c  Applied Numerical Methods, Example 5.1
c  Gauss-Jordan Reduction
c
c  This program finds M solution vectors corresponding to a set of N
c  simultaneous linear equations using the Gauss-Jordan reduction
c  algorithm with the diagonal pivot strategy.  The N by N matrix of 
c  coefficients appears in the first N columns of the matrix A.  M 
c  right-hand side vectors are appended as the (N+1)th to (N+M)th 
c  columns of A.  On the K(th) pass of the elimination scheme, K(th) 
c  row elements are normalized by dividing by the pivot element A(K,K).
c  DETER, the determinant of the coefficient matrix, is updated
c  prior to elimination of all nonzero elements (except for the pivot
c  element) in the K(th) column.  Should a pivot element be encountered
c  which is smaller in magniture than EPS, computation is disco,ntinued
c  and an appropriate comment is printed.  The solution vectors 
c  corresponding to the original right-hand side vectors are stored
c  in columns N+I to N+M of A.
c
      Implicit Real *8(a-h,o-z)
      dimension A(50,51)
c
1     read(5,100)N,M, EPS
      NPLUSM=N+M
      write(6,200)N,M, EPS
      do 2 I=1,N
      read(5,101)(A(I,J),J=1,NPLUSM)
      write(6,201)(A(I,J),J=1,NPLUSM)
2     continue
c
c  ....Begin Elimination Procedure
c
      DETER=1.
      do 9 K=1,N
c
c  ....Update the determinant value
      DETER=DETER*A(K,K)
c
c  ....Check for pivot element too small
      If (DABS(A(K,K)).GT.EPS) GO TO 5
      write(6,202)
      GO TO 1
c
c  ....Normalize the pivot row
5     KP1=K+1
      do 6 J=KP1,NPLUSM
6     A(K,J)=A(K,J)/A(K,K)
      A(K,K)=1.
c
c  ....Eliminate K(th) column elements except for pivot
      do 9 I=1,N
      IF(I.EQ.K .OR. A(I,K).EQ.0.) GO TO 9
          do 8 J=KP1,NPLUSM
8         A(I,J)=A(I,J)-A(I,K)*A(K,J)
          A(I,K)=0.
9     continue
c
      write(6,203)DETER,N,NPLUSM
      do 10 I=1,N
10    write(6,201)(A(I,J),J=1,NPLUSM)
      go to 1
c
c  ....Formats for input and output statements
100   Format( 2(10x,I2,8x),10x,E7.1)
101   Format( 20x,5F10.5)
200   Format( ' N = ',I8 / ' M = ', I8 / ' EPS = ',1PE14.1,
     & '  A(1,1)...A(..', / )
201   Format(1H ,7F13.7)
202   Format(' Small pivot - Matrix may be singular')
203   Format(' DETER = ', E115.6/ '       A(1,1)....A(',I2, 1H,,I2,1H))
c
      end
