program simq5
c
c -- solve simultaneous equations by
c -- gauss-jordan elimination with multiple constant vectors
c -- subroutines gaussj and swap required
c -- figure 4.9
c     
      logical error
      integer size, maxr, maxc, out, nvec, index(8,3)
      real a(8,8), y(8,3), coef(8,3), b(8,8)
      common /inout/ out, maxr, maxc, error
c
      out = 6
      maxr = 8
      maxc = 8
      write(out, 101)
10    call input(a, y, size, nvec)
      if (size .lt. 2) goto 100
      do 30 i = 1, size
       do 20 j = 1, size
        b(i,j) = a(i,j)
20     continue
       do 30 j = 1, nvec
        coef(i,j) = y(i,j)
30    continue
      call gaussj(b, coef, index, size, maxr, nvec, error, out)
      if (.not. error)
     *  call output(a, y, coef, b, size, nvec)
      goto 10
100   stop
101   format('1 simultaneous solution by',
     * ' gauss-jordan elimination')
      end
      subroutine input(a, y, n, nvec)
c --  get values for n and arrays a, y
c
      integer n, out, i, j, m, maxr, nvec
      real a(8,8), y(8,3)
      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
 7    write(out, 107)
      read(*, 106) nvec
      if ((nvec .gt. 3) .or. (nvec.lt.0)) goto 7
      do 30 i = 1, n
       write(out, 101) i
       do 10 j = 1, n
        write(out, 102) j
        read(*, 103) a(i,j)
10     continue
       if (nvec .eq. 0) goto 30
       do 20 j = 1, nvec
        write(out, 104)
        read(*, 103) y(i,j)
20     continue
30    continue
      return
101   format(' equation ', i3/)
102   format('+',i4, ': ' )
103   format(f10.0)
104   format('+ c: ' )
105   format(/' how many equations? ' )
106   format(i2)
107   format(' how many constant vectors? ' )
      end
      subroutine output(a, y, coef, inver, n, nvec)
c
c -- print the answers
c
      logical error
      integer n, out, i, j, maxr, maxc, nvec
      real a(8,8), y(8,3), coef(8,3), inver(8,8)
      common /inout/ out, maxr, maxc, error
c
      if (nvec .eq. 0) goto 30
c -- print matrix and solution vector
      do 10 i = 1, n
       write(out, 101)
     *   (a(i,j), j = 1, n), (y(i,j), j = 1, nvec)
10    continue
      if (error) return
      write(out, 102)
      do 20 i = 1, n
       write(out, 101) (coef(i,j), j = 1, nvec)
20    continue
      return
c -- print inverse
30    write(out, 103)
      do 40 i = 1, n
       write(out, 101)
     *   (inver(i,j), j = 1, n)
40    continue
      return
101   format(1p6e12.4)
102   format(/ '   solution')
103   format(/ '   matrix inverse')
      end