program simq1
c
c -- fortran program to solve three
c -- simultaneous equations by cramer's rule
c -- figure 4.3
c     
      logical error
      integer length, out
      character*1 yesno
      real a(3,3), y(3), coef(3)
      common /inout/ out
c
      out = 6
10    write(out,*)' simultaneous solution by cramers rule'
      call input(a, y, length)
      call cramer(a, y, coef, error)
      if (.not. error) call output(a, y, coef, length)
      write(out, 102)
      read(*, 103) yesno
      if (yesno .eq. 'y' .or. yesno .eq. 'Y') goto 10
      stop
102   format(/' more? ' )
103   format(a1)
      end
      subroutine input(a, y, n)
c -- get values for n and arrays a, y
c
      integer n, out, i, j
      real a(3,3), y(3)
      common /inout/ out
c
      n = 3
      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: ' )
      end
      subroutine cramer(a, y, coef, error)
c
c -- solve 3 simultaneous equations by cramer-s rule
c
      logical error
      integer n, i, j, out
      real a(3,3), y(3), coef(3), b(3,3), det
      common /inout/ out
c
      error = .false.
      n = 3
      do 20 i = 1, n
       do 10 j = 1, n
        b(i,j) = a(i,j)
10     continue
20    continue
      det = determ(b)
      if (det .eq. 0.0) goto 90
      call setup(det, a, b, y, coef,1)
      call setup(det, a, b, y, coef,2)
      call setup(det, a, b, y, coef,3)
      return
c
90    error = .true.
      write(out,*) ' error--matrix singular'
      return
      end
      function determ(x)
c
c -- find the determinant of the 3-by-3 matrix x
c
      real x(3,3), sum
c
      sum = x(1,1) * (x(2,2)*x(3,3) - x(3,2)*x(2,3))
      sum = sum - x(1,2)* (x(2,1)*x(3,3) - x(3,1)*x(2,3))
      sum = sum + x(1,3)* (x(2,1)*x(3,2) - x(3,1)*x(2,2))
      determ = sum
      return
      end
      subroutine setup(det, a, b, y, coef, j)
c
c -- set up the numerator matrices
c
      integer i, j
      real a(3,3), b(3,3), y(3), coef(3), det
c
      do 10 i = 1, 3
       b(i,j) = y(i)
       if (j .gt. 1) b(i, j-1) = a(i, j-1)
10    continue
      coef(j) = determ(b) / det
      return
      end
      subroutine output(a, y, coef, n)
c
c -- print out the answers
c
      integer n, out, i, j
      real a(n,n), y(n), coef(n)
      common /inout/ out
c
      write(out, 101) ((a(i,j), j = 1, n), y(i), i = 1, n)
      write(out,*) '   solution'
      write(out, 101) (coef(i), i = 1, n)
      return
101   format(1p4e12.4)
      end