C---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
      PROGRAM LINSYS

      !Program:  Solves for vector X in the linear equation AX=B where
      !A is an N by N non-singular matrix.  You must first input the
      !elements for matrix A followed by the elements of column vector
      !B.  The method used is Gauss-Jordan.  The parameter "NMAX=80"
      !below is arbitrary and represents the largest matrix A linear
      !system that can be solved.
      !
      !Example (three independent linear variables, X1, X2, and X3):
      !
      !             7 X1    -2 X2 +    X3 = -5
      !            -2 X1 +  10 X2 -  5 X3 = 20
      !             4 X1 - 400 X2 + 17 X3 = -200
      !
      !The solution is X1=-0.151515, X2=0.361815, X3=-3.21576 to six
      !significant digits.
      !
      ! Author:  Jerry R. Ziemke
      
      PARAMETER(NMAX=80)
      DOUBLE PRECISION B(NMAX),A(NMAX,NMAX+1),W(NMAX,NMAX+1),X(NMAX)

      WRITE(*,*)'This program solves for vector X in the linear'
      WRITE(*,*)'equation AX=B where A is an N by N non-singular'
      WRITE(*,*)'matrix.  You must first input the elements for'
      WRITE(*,*)'matrix A followed by the elements of column vector B.'
      WRITE(*,*)'Method used: Basic Gauss-Jordan.'
      WRITE(*,*)'Enter N for N X N matrix A:'
      READ(*,*) N
      DO ICOL=1,N
	WRITE(*,*)'Column:',ICOL
	DO IROW=1,N
	  WRITE(*,'(1X,A8,I2,A1,I2,A2)') 'Enter A(',IROW,',',ICOL,'):'
	  READ(*,*) A(IROW,ICOL)
	ENDDO
      ENDDO
      DO IROW=1,N
	WRITE(*,'(1X,A8,I2,A2)') 'Enter B(',IROW,'):'
	READ(*,*) B(IROW)
      ENDDO
      CALL LNSYS(N,A,W,B,1,1,X)

      STOP
      END

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
      SUBROUTINE LNSYS(N,A,W,B,IPRINT,ICHECK,X)

C Subroutine: Computes the solution vector column "X" for the linear
C	equation AX=B where A is a (real) N by N matrix.  Gauss-Jordan
C       method is used.
C Arrays, etc. (NOTE: You must declare in your calling program that A,
C       W, B and X arrays are DOUBLE PRECISION):
C   N = Dimension of linear system to be solved.  Note that the arrays
C	A and W are dimensioned N rows by N+1 columns (input).
C   A = Matrix containing data coefficients (input).
C   W = Just a "work" array entirely internal to LNSYS.  This array
C	must be declared in an external call and with the same
C       dimensions as array A.
C   B = Column vector of constants (input).
C   IPRINT = 1 to print out inverse elements of A, and any other
C       integer to skip doing this (input).
C   ICHECK = 1 for an error check; any other integer skips doing this.
C	This error analysis computes the product AX (X is the computed
C       solution), yielding a column of numbers which SHOULD be equal
C       to the original column vector B (input).
C   X = Solution column vector (output).

      PARAMETER(NMAX=80)
      DOUBLE PRECISION A(NMAX,NMAX+1),W(NMAX,NMAX+1),X(NMAX),B(NMAX)

      DO 20 I=1,N
        DO 10 J=1,N
          W(I,J)=A(I,J)
 10     CONTINUE
 20   CONTINUE
      DO 30 J=1,N
        A(J,N+1)=B(J)
 30   CONTINUE
      DO 60 K=2,N
        IF (A(K-1,K-1).EQ.0.0) CALL ADDEM(K,N,A)
        DO 50 I=K,N
          U=-A(I,K-1)/A(K-1,K-1)
          DO 40 J=K,N+1
            A(I,J)=A(I,J)+U*A(K-1,J)
 40       CONTINUE
 50     CONTINUE
 60   CONTINUE
      X(N)=A(N,N+1)/A(N,N)
      DO 80 JJ=1,N-1
	J=N-JJ
	X(J)=0.0
        DO 70 I=J+1,N
          X(J)=X(J)+A(J,I)*X(I)
 70     CONTINUE
        X(J)=(A(J,N+1)-X(J))/A(J,J)
 80   CONTINUE
      IF (IPRINT.EQ.1) THEN
        WRITE(*,*) 'SOLUTION X(I)...'
        WRITE(*,*) 'I, X(I):'
        DO 90 I=1,N
          WRITE(*,*) I,X(I)
 90     CONTINUE
      ENDIF
      IF (ICHECK.EQ.1) THEN
	WRITE(*,*) 'Solution (X) error check...'
	WRITE(*,*) '(Derived column vector B should equal actual B)'
	WRITE(*,*) 'Row index I, Actual B(I), Derived B(I) (from AX=B):'
        DO 110 I=1,N
          SS=0.0
          DO 100 J=1,N
            SS=SS+W(I,J)*X(J)
 100      CONTINUE
          WRITE(*,*) I,B(I),SS
 110    CONTINUE
      ENDIF

      RETURN
      END

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7--
      SUBROUTINE ADDEM(K,N,A)
      
      DOUBLE PRECISION A(N,N+1)

      L=1
 10   IF ((L.LT.N).AND.(A(L,K-1).EQ.0.0)) L=L+1
      IF ((L.LT.N).AND.(A(L,K-1).EQ.0.0)) GOTO 10
      DO 20 M=K-1,N+1
        A(K-1,M)=A(K-1,M)+A(L,M)
 20   CONTINUE

      RETURN
      END