program simq12
c
c -- solve simultaneous equations by gauss-seidel
c -- figure 4.15
c     
      logical error
      integer length, maxr, maxc, out
      real a(8,8), y(8), coef(8)
      common /inout/ out, maxr, maxc
c
      out = 6
      maxr = 8
      maxc = 8
      write(out, 101)
10    call input(a, y, length)
      if (length .lt. 2) stop
      call seid(a, y, coef, length, error)
      if (.not. error) call output(a, y, coef, length)
      goto 10
101   format('1 simultaneous solution by gauss-seidel')
      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, maxc
c
5     write(out, 105)
      read(*, 106) n
      m = n
      if (n .gt. maxr) goto 5
      if (n .eq. 0) return
      if (n .lt. 0) goto 30
      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
c -- use matrix from previous problem
30    n = -n
      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 seid(a, y, coef, ncol, error)
c
c -- solve simultaneous equations by gauss-seidel
c -- apr 28, 81
c
      logical error, done
      integer ncol, i, j, out, n, n1, k, l, max
      real a(8,8), y(8), coef(8), tol, lambda
      real b(8,8), w(8), big, ab, sum, nextc
      common /inout/ out, maxr, maxc
      data tol/ 1.0e-4/, max/100/
c
      error = .false.
5     write(out, 101)
      read(*,102) lambda
      if ((lambda .lt. 0.0).or.(lambda .gt. 2.0)) goto 5
      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 99
       if (l .eq. i) goto 70
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))
70    continue
      if (b(n,n) .eq. 0.0) goto 99
c -- initial values
      do 75 i = 1, n
      coef(i) = 0.0
75    continue
      i = 0
80    i = i + 1
      done = .true.
      do 90 j = 1, n
       sum = y(j)
       do 85 k = 1, n
        if (j .ne. k) sum = sum - b(j,k) * coef(k)
85     continue
       nextc = sum / b(j,j)
       if (abs(nextc - coef(j)) .lt. tol) goto 95
       done = .false.
       if (nextc * coef(j) .lt. 0.0)
     *    nextc = (coef(j) + nextc) / 2
       coef(j) = lambda * nextc + (1 - lambda)*coef(j)
       write(out, 103) i, j, coef(j)
90    continue
      if ((i .le. 100).and.(.not. done)) goto 80
      write(out,*) ' no solution'
      error = .true.
95    return
99    write(out,*) ' error--matrix singular'
      error = .true.
      return
101   format(' relaxation factor? ' )
102   format(e10.0)
103   format(i4, '; coef(', i3, ') =', 1pe12.4)
      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, maxc
c
      do 10 i = 1, n
       write(out, 101) (a(i,j), j = 1, n), y(i)
10    continue
      write(out,102)
      write(out, 101) (coef(i), i = 1, n)
      return
101   format(1p6e12.4)
102   format(/ '   solution')
      end