program simq11
c
c -- solve simultaneous equations with complex coefficients
c -- by gauss-jordan elimination
c -- subroutines gaussj and swap required
c -- figure 4.4
c     
      logical error
      integer n, m, maxr, maxc, out, i, j, index(8,3), nvec
      real a(8,8), y(8), coef(8), b(8,8)
      common /inout/ out, maxr, maxc, error
      data nvec/1/
c
      out = 6
      maxr = 4
      maxc = 4
      write(out, 101)
10    call input(a, y, n, m)
      if (n .lt. 2) goto 100
      do 30 i = 1, n
       do 20 j = 1, n
        b(i,j) = a(i,j)
20     continue
       coef(i) = y(i)
30    continue
      call gaussj(b, coef, index, n, maxr*2, nvec, error, out)
      if (.not. error) call output(a, y, coef, n)
      goto 10
100   stop
101   format('1 simultaneous solution',
     * ' with complex coefficients')
      end
      subroutine input(a, y, n, m)
c --  get values for n and arrays a, y
c
      integer n, out, i, j, m, maxr, k, l
      real a(8,8), y(8), dreal(4,4), dimag(4,4), v(4,2)
      common /inout/ out, maxr, maxc, error
c
5     write(out, 105)
      read(*, 106) n
      m = n
      if (n .gt. maxr) goto 5
      if (n .lt. 2) return
      do 20 i = 1, n
       write(out, 101) i
       k = 0
       l = 2 * i - 1
       do 10 j = 1, n
        k = k + 1
        write(out, 102) j
        read(*, 103) dreal(i,j)
        a(l,k) = dreal(i,j)
        a(l+1, k+1) = dreal(i,j)
        k = k + 1
        write(out, 107) j
        read(*, 103) dimag(i,j)
        a(l,k) = -dimag(i,j)
        a(l+1, k-1) = dimag(i,j)
10     continue
       write(out, 104)
       read(*, 103) v(i,1)
       y(l) = v(i,1)
       write(out, 108)
       read(*, 103) v(i,2)
       y(l+1) = v(i,2)
20    continue
c print original matrix
      do 30 i = 1, n
       write(out, 109) (dreal(i,j), dimag(i,j), j = 1, n),
     * v(i,1), v(i,2)
30    continue
      write(out, 110)
      n = 2 * n
      m = n
      return
101   format(' equation ', i3/)
102   format('+real ', i3, ': ' )
103   format(f10.0)
104   format('+real constant ' )
105   format(' how many equations? ' )
106   format(i2)
107   format('+imag ', i3, ': ' )
108   format('+imag constant ' )
109   format(1p6e12.4)
110   format(/)
      end
      subroutine output(a, y, coef, n)
c
c -- print the answers
c
      logical error
      integer n, out, i, j, maxr, maxc, n2
      real a(8,8), y(8), coef(8), re, im, mag, ang, p180
      common /inout/ out, maxr, maxc, error
      data p180/57.29578/
c
      do 10 i = 1, n
       write(out, 101) (a(i,j), j = 1, n), y(i)
10    continue
      if (error) return
      n2 = n / 2
      write(out, 102)
      do 20 i = 1, n2
       j = 2 * i - 1
       re = coef(j)
       im = coef(j + 1)
       mag = sqrt(re*re + im*im)
       ang = atan2(im, re) * p180
       write(out, 103) re, im, mag, ang
20    continue
      return
101   format(1p6e12.4)
102   format(/ '    real       imaginary',
     * '   magnitude    angle')
103   format(1p4e12.4)
      end