program simq2
c
c -- fortran program to solve
c -- simultaneous equations by gauss elimination
c -- figure 4.4
c     
      logical error
      integer length, maxr, out
      real a(8,8), y(8), coef(8)
      common /inout/ out, maxr
c
      out = 6
      maxr = 8
      write(out, 101)
10    call input(a, y, length)
      if (length .lt. 2) stop
      call gauss(a, y, coef, length, error)
      if (.not. error) call output(a, y, coef, length)
      goto 10
101   format('1 simultaneous solution by gauss elimination')
      end
      subroutine input(a, y, n)
c -- get values for n and arrays a, y
c
      integer n, out, i, j, m
      real a(8,8), y(8)
      common /inout/ out, maxr
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
       do 10 j = 1, n
        write(out, 102) j
        read(*, 103) a(i,j)
10     continue
       write(out, 104)
       read(*, 103) y(i)
20    continue
      return
101   format(' equation ', i3/)
102   format('+',i4, ': ' )
103   format(f10.0)
104   format('+ c: ' )
105   format(/'+how many equations? ' )
106   format(i2)
      end
      subroutine gauss(a, y, coef, ncol, error)
c
c -- simultaneous solution by gauss elimination
c
      logical error
      integer ncol, i, j, out, n, n1, i1, k, l
      real a(8,8), y(8), coef(8)
      real b(8,8), w(8), big, ab, sum, t
      common /inout/ out, maxr
c
      error = .false.
      n = ncol
      do 20 i = 1, n
       do 10 j = 1, n
        b(i,j) = a(i,j)
10     continue
       w(i) = y(i)
20    continue
      n1 = n - 1
      do 70 i = 1, n1
       big = abs(b(i,i))
       l = i
       i1 = i + 1
       do 30 j = i1, n
        ab = abs(b(j,i))
        if (ab .le. big) goto 30
         big = ab
         l = j
30     continue
       if (big .eq. 0.0) goto 70
       if (l .eq. i) goto 50
c -- interchange rows to put largest element on diagonal
       do 40 j = 1, n
        call swap(b(l,j), b(i,j))
40     continue
       call swap(w(i), w(l))
50     continue
       do 60 j = i1, n
        t = b(j,i) / b(i,i)
        do 55 k = i1, n
         b(j,k) = b(j,k) - t * b(i,k)
55      continue
        w(j) = w(j) - t * w(i)
60     continue
70    continue
      if (b(n,n) .eq. 0.0) goto 99
      coef(n) = w(n) / b(n,n)
      i = n - 1
c -- back substitution
80    sum = 0
       i1 = i + 1
       do 75 j = i1, n
        sum = sum + b(i,j) * coef(j)
75     continue
       coef(i) = (w(i) - sum) / b(i,i)
       i = i - 1
      if (i .gt. 0) goto 80
      return
99    write(out,*)' error--matrix singular'
      return
      end
      subroutine swap(a,b)
c -- interchange two values
      real a, b, hold
c
      hold = a
      a = b
      b = hold
      return
      end
      subroutine output(a, y, coef, n)
c
c -- print out the answers
c
      integer n, out, i, j
      real a(8,8), y(8), coef(8)
      common /inout/ out, maxr
c
      do 10 i = 1, n
       write(out, 101) (a(i,j), j = 1, n), y(i)
10    continue
      write(out,*) '   solution'
      write(out, 101) (coef(i), i = 1, n)
      return
101   format(1p6e12.4)
      end