      PROGRAM BACSUBS
C     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
C     To accompany the text:
C     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
C     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
C     This free software is complements of the author.
C
C     Algorithm 3.1 (Back Substitution).
C     Section 3.3, Upper-Triangular Linear Systems, Page 145
C
      PARAMETER(MaxR=10)
      INTEGER InRC,N
      REAL A,B,DET,X
      CHARACTER Ans*1,Ach*1,Bch*1,Xch*1
      DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),X(1:MaxR)
10    CALL INPUTS(A,B,N,InRC,DET)
      CALL BACKSUB(A,B,N,X,DET)
      CALL RESULTS(A,B,N,X,DET)
      WRITE(9,*)' '
      WRITE(9,*)'WANT TO SOLVE ANOTHER UPPER-TRIANGULAR SYSTEM ? <Y/N> '
      READ(9,'(A)') Ans
      IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') GOTO 10
      STOP
      END

      SUBROUTINE BACKSUB(A,B,N,X,DET)
      PARAMETER(MaxR=10)
      INTEGER J,K,N,R
      REAL A,B,DET,SUM,X
      DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),X(1:MaxR)
      DO K=1,N
      IF (A(K,K).EQ.0) THEN
         DET=0
         RETURN
      ENDIF
      ENDDO
      DET=A(N,N)
      X(N)=B(N)/A(N,N)
      DO R=(N-1),1,-1
        DET=DET*A(R,R)
        SUM= 0
        DO J=(R+1),N
          SUM=SUM+A(R,J)*X(J)
        ENDDO
      X(R)= (B(R)-SUM)/A(R,R)
      ENDDO
      RETURN
      END

      SUBROUTINE XBACKSUB(A,B,N,X,DET)
C     This subroutine uses simulated WHILE loop(s).
      PARAMETER(MaxR=10)
      INTEGER J,K,N,R
      REAL A,B,DET,SUM,X
      DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),X(1:MaxR)
      DO 10 K=1,N
      IF (A(K,K).EQ.0) THEN
         DET=0
         RETURN
      ENDIF
10    CONTINUE
      DET=A(N,N)
      X(N)=B(N)/A(N,N)
      DO 30 R=(N-1),1,-1
        DET=DET*A(R,R)
        SUM= 0
        DO 20 J=(R+1),N
          SUM=SUM+A(R,J)*X(J)
20      CONTINUE
      X(R)= (B(R)-SUM)/A(R,R)
30    CONTINUE
      RETURN
      END

      SUBROUTINE INPUTS(A,B,N,InRC,DET)
      PARAMETER(MaxR=10)
      INTEGER C,I,InRC,N,R
      REAL A,B,DET
      CHARACTER Ans*1,Ach*1
      DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR)
      DO 10 I=1,15
        WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)'SOLUTION OF AN UPPER-TRIANGULAR LINEAR SYSTEM:'
      WRITE(9,*)' '
      WRITE(9,*)'      A*X = B'
      WRITE(9,*)' '
      WRITE(9,*)'A  IS AN  N BY N  UPPER-TRIANGULAR MATRIX.'
      WRITE(9,*)' '
      WRITE(9,*)'B  IS AN  N  DIMENSIONAL VECTOR OF CONSTANTS.'
      WRITE(9,*)' '
      WRITE(9,*)'X  IS THE N DIMENSIONAL SOLUTION VECTOR.'
      WRITE(9,*)' '
      WRITE(9,*)'ENTER THE NUMBER OF EQUATIONS N = '
      READ(9,*) N
      InRC=0
      WRITE(9,*)'WANT TO INPUT COLUMNS ? <Y/N> '
      READ(9,'(A)') Ans
      IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') THEN
        InRC=1
      ENDIF
      WRITE(9,*)' '
      WRITE(9,*)'ENTER THE UPPER-TRIANGULAR MATRIX A(I,J):'
      CALL MATIN('A',A,N,InRC)
      WRITE(9,*)' '
      WRITE(9,*)'ENTER THE COLUMN VECTOR B(J):'
      CALL VECIN('B',B,N)
      RETURN
      END

      SUBROUTINE MATIN(Ach,A,N,InRC)
      PARAMETER(MaxR=10)
      INTEGER Col,InRC,N,Row
      REAL A
      CHARACTER Ach*1
      DIMENSION A(1:MaxR,1:MaXR)
      DO 20 Row=1,N
      DO 10 Col=1,N
        A(Row,Col)=0
10    CONTINUE
20    CONTINUE
      WRITE(9,*)'GIVE THE ELEMENTS OF THE MATRIX ',Ach
      WRITE(9,*)' '
      IF (InRC.EQ.0) THEN
        DO 40 Row=1,N
          WRITE(9,*)'GIVE THE ELEMENTS OF ROW ',Row
          WRITE(9,*)' '
        DO 30 Col=Row,N
          WRITE(9,*) Ach,'(',Row,',',Col,') = '
          READ(9,*) A(Row,Col)
30      CONTINUE
40      CONTINUE
      ELSE
        DO 60 Col=1,N
          WRITE(9,*)'GIVE THE ELEMENTS OF COLUMN ',Col
          WRITE(9,*)' '
        DO 50 Row=1,Col
          WRITE(9,*) Ach,'(',Row,',',Col,') = '
          READ(9,*) A(Row,Col)
50      CONTINUE
60      CONTINUE
      ENDIF
      RETURN
      END

      SUBROUTINE VECIN(Bch,B,N)
      PARAMETER(MaxR=10)
      INTEGER N,Row
      REAL B
      CHARACTER Bch*1
      DIMENSION B(1:MaxR)
      WRITE(9,*)'GIVE THE ELEMENTS OF THE COLUMN VECTOR ',Bch
      WRITE(9,*)' '
      DO 10 Row=1,N
        WRITE(9,*) Bch,'(',Row,') = '
        READ(9,*) B(Row)
10    CONTINUE
      RETURN
      END

      SUBROUTINE MATPRINT(A,N)
      PARAMETER(MaxR=10)
      INTEGER Col,N,Row
      REAL A
      DIMENSION A(1:MaxR,1:MaxR)
      DO 10 Row=1,N
        WRITE(9,999) ( A(Row,Col), Col=1,N )
10    CONTINUE
999   FORMAT(1X,5F15.7,4X)
      RETURN
      END

      SUBROUTINE VECPRINT(Bch,B,Xch,X,N)
      PARAMETER(MaxR=10)
      INTEGER N,Row
      REAL B,X
      CHARACTER Bch*1,Xch*1
      DIMENSION B(1:MaxR),X(1:MaxR)
      WRITE(9,*) Bch,' COEFFICIENT VECTOR                ',Xch,
     +' SOLUTION VECTOR'
      WRITE(9,*)
      DO 10 Row=1,N
        WRITE(9,*) Bch,'(',Row,') =',B(Row),'                     ',
     +Xch,'(',Row,') =',X(Row)
10    CONTINUE
      RETURN
      END

      SUBROUTINE RESULTS(A,B,N,X,DET)
      PARAMETER(MaxR=10)
      INTEGER Col,I,N,Row
      REAL A,B,DET,X
      CHARACTER Bch*1,Xch*1
      DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),X(1:MaxR)
      DO 10 I=1,15
        WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)'THE UPPER-TRIANGULAR MATRIX A(I,J) IS:'
      WRITE(9,*)' '
      CALL MATPRINT(A,N)
      IF (DET.EQ.0) THEN
        WRITE(9,*)'THE MATRIX IS SINGULAR.'
        WRITE(9,*)
        WRITE(9,*)'THE BACK-SUBSTITUTION ALGORITHM DOES NOT APPLY.'
      ELSE
        WRITE(9,*)
        CALL VECPRINT('B',B,'X',X,N)
      ENDIF
      WRITE(9,*)' '
      WRITE(9,*)'THE DETERMINANT`S VALUE IS  DET A = ',DET
      RETURN
      END

      PROGRAM GAUSJORD
C     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
C     To accompany the text:
C     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
C     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
C     This free software is complements of the author.
C
C     Algorithm 3.2 (Upper-Triangularization Followed by Back Substitution).
C     Section 3.4, Gaussian Elimination and Pivoting, Page 156
C
      PARAMETER(MaxR=10,MaxC=11)
      INTEGER InRC,N
      REAL A,A1,DET,X
      CHARACTER Ans*1,Ach*1,Bch*1,Xch*1
      DIMENSION A(1:MaxR,1:MaxC),A1(1:MaxR,1:MaxC),X(1:MaxR)
10    CALL INPUTS(A,A1,N,InRC)
      CALL GAUSSJO(A,X,N,DET)
      CALL RESULTS(A1,X,N,DET)
      WRITE(9,*)' '
      WRITE(9,*)'WANT TO SOLVE ANOTHER LINEAR SYSTEM ? <Y/N> '
      READ(9,'(A)') Ans
      IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') GOTO 10
      STOP
      END

      SUBROUTINE GAUSSJO(A,X,N,DET)
      PARAMETER(MaxR=10,MaxC=11)
      INTEGER C,J,K,N,P,R,Row,T
      REAL A,DET,M,SUM,X
      DIMENSION A(1:MaxR,1:MaxC),Row(1:MaxR),X(1:MaxR)
      DO J=1,N
        Row(J)=J
      ENDDO
      DO P=1,(N-1)
        DO K=(P+1),N
          IF (ABS(A(Row(K),P)).GT.ABS(A(Row(P),P))) THEN
            T = Row(P)
            Row(P)=Row(K)
            Row(K)=T
          ENDIF
        ENDDO
        IF (A(Row(P),P).EQ.0) THEN
            DET=0
            RETURN
        ENDIF
        DO K=(P+1),N
          M = A(Row(K),P)/A(Row(P),P)
          DO C=(P+1),(N+1)
            A(Row(K),C)=A(Row(K),C)-M*A(Row(P),C)
          ENDDO
        ENDDO
      ENDDO
      DET=A(Row(N),N)
      IF (DET.NE.0) THEN
        X(N)=A(Row(N),(N+1))/A(Row(N),N)
        DO K=(N-1),1,-1
          DET=DET*A(Row(K),K)
          SUM=0
          DO C=(K+1),N
            SUM=SUM+A(Row(K),C)*X(C)
          ENDDO
          X(K)=(A(Row(K),(N+1))-SUM)/A(Row(K),K)
        ENDDO
      ENDIF
      RETURN
      END

      SUBROUTINE XGAUSSJO(A,X,N,DET)
C     This subroutine uses simulated WHILE loop(s).
      PARAMETER(MaxR=10,MaxC=11)
      INTEGER C,J,K,N,P,R,Row,T
      REAL A,DET,M,SUM,X
      DIMENSION A(1:MaxR,1:MaxC),Row(1:MaxR),X(1:MaxR)
      DO 10 J=1,N
        Row(J)=J
10    CONTINUE
      DO 50 P=1,(N-1)
      DO 20 K=(P+1),N
        IF (ABS(A(Row(K),P)).GT.ABS(A(Row(P),P))) THEN
          T = Row(P)
          Row(P)=Row(K)
          Row(K)=T
        ENDIF
20    CONTINUE
        IF (A(Row(P),P).EQ.0) THEN
            DET=0
            RETURN
        ENDIF
        DO 40 K=(P+1),N
          M = A(Row(K),P)/A(Row(P),P)
          DO 30 C=(P+1),(N+1)
          A(Row(K),C)=A(Row(K),C)-M*A(Row(P),C)
30        CONTINUE
40      CONTINUE
50    CONTINUE
      DET=A(Row(N),N)
      IF (DET.NE.0) THEN
        X(N)=A(Row(N),(N+1))/A(Row(N),N)
        DO 70 K=(N-1),1,-1
          DET=DET*A(Row(K),K)
          SUM=0
          DO 60 C=(K+1),N
            SUM=SUM+A(Row(K),C)*X(C)
60        CONTINUE
          X(K)=(A(Row(K),(N+1))-SUM)/A(Row(K),K)
70      CONTINUE
      ENDIF
      RETURN
      END

      SUBROUTINE INPUTS(A,A1,N,InRC)
      PARAMETER(MaxR=10,MaxC=11)
      INTEGER Col,I,InRC,N,Row
      REAL A
      CHARACTER Ans*1,Ach*1
      DIMENSION A(1:MaxR,1:MaxC),A1(1:MaxR,1:MaxC)
      DO 10 I=1,15
        WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)'SOLUTION OF A LINEAR SYSTEM:'
      WRITE(9,*)' '
      WRITE(9,*)'      A*X = B'
      WRITE(9,*)' '
      WRITE(9,*)'GAUSSIAN ELIMINATION IS USED TO UPPER-TRIANGULARIZE,'
      WRITE(9,*)' '
      WRITE(9,*)'AND THEN BACK-SUBSTITUTION IS USED TO SOLVE FOR  X.'
      WRITE(9,*)' '
      WRITE(9,*)'A  IS AN  N BY N  NONSINGULAR MATRIX.'
      WRITE(9,*)' '
      WRITE(9,*)'B  IS AN  N  DIMENSIONAL VECTOR OF CONSTANTS.'
      WRITE(9,*)' '
      WRITE(9,*)'X  IS THE N DIMENSIONAL SOLUTION VECTOR.'
      WRITE(9,*)' '
      WRITE(9,*)'ENTER THE NUMBER OF EQUATIONS N = '
      READ(9,*) N
      InRC=0
      WRITE(9,*)'WANT TO INPUT COLUMNS ? <Y/N> '
      READ(9,'(A)') Ans
      IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') THEN
        InRC=1
      ENDIF
      WRITE(9,*)' '
      WRITE(9,*)'ENTER THE COEFFICIENT MATRIX  A:'
      CALL MATIN('A',A,N,InRC)
      WRITE(9,*)' '
      WRITE(9,*)'ENTER THE COLUMN VECTOR  B:'
      CALL VECIN('B',A,N,InRC)
      DO 30 Row=1,N
      DO 20 Col=1,N+1
        A1(Row,Col)=A(Row,Col)
20    CONTINUE
30    CONTINUE
      RETURN
      END

      SUBROUTINE MATIN(Ach,A,N,InRC)
      PARAMETER(MaxR=10,MaxC=11)
      INTEGER Col,InRC,N,Row
      REAL A
      CHARACTER Ach*1
      DIMENSION A(1:MaxR,1:MaxC)
      DO 20 Row=1,N
        DO 10 Col=1,(N+1)
          A(Row,Col)=0
10      CONTINUE
20    CONTINUE
      WRITE(9,*)'GIVE THE ELEMENTS OF THE MATRIX ',Ach
      WRITE(9,*)' '
      IF (InRC.EQ.0) THEN
        DO 40 Row=1,N
          WRITE(9,*)'GIVE THE ELEMENTS OF ROW ',Row
          WRITE(9,*)' '
          DO 30 Col=1,N
            WRITE(9,*) Ach,'(',Row,',',Col,') = '
            READ(9,*) A(Row,Col)
30        CONTINUE
40      CONTINUE
      ELSE
        DO 60 Col=1,N
          WRITE(9,*)'GIVE THE ELEMENTS OF COLUMN ',Col
          WRITE(9,*)' '
          DO 50 Row=1,N
            WRITE(9,*) Ach,'(',Row,',',Col,') = '
            READ(9,*) A(Row,Col)
50        CONTINUE
60      CONTINUE
      ENDIF
      RETURN
      END

      SUBROUTINE VECIN(Ach,A,N,InRC)
      PARAMETER(MaxR=10,MaxC=11)
      INTEGER Col,InRC,N,Row
      REAL A
      CHARACTER Ach*1
      DIMENSION A(1:MaxR,1:MaxC)
      WRITE(9,*)'GIVE THE ELEMENTS OF THE COLUMN VECTOR ',Ach
      WRITE(9,*)' '
      Col=N+1
      DO 10 Row=1,N
        WRITE(9,*) Ach,'(',Row,') = '
        READ(9,*) A(Row,Col)
10    CONTINUE
      RETURN
      END

      SUBROUTINE MATPRINT(A,N)
      PARAMETER(MaxR=10,MaxC=11)
      INTEGER Col,N,Row
      REAL A
      DIMENSION A(1:MaxR,1:MaxC)
      DO 10 Row=1,N
        WRITE(9,999) ( A(Row,Col), Col=1,N )
10    CONTINUE
999   FORMAT(1X,5F15.7,4X)
      RETURN
      END

      SUBROUTINE VECPRINT(Bch,Xch,A,X,N)
      PARAMETER(MaxR=10,MaxC=11)
      INTEGER Col,N,Row
      REAL A,X
      CHARACTER Bch*1,Xch*1
      DIMENSION A(1:MaxR,MaxC),X(1:MaxR)
      Col=N+1
      WRITE(9,*) Bch,' COEFFICIENT VECTOR                ',Xch,
     +' SOLUTION VECTOR'
      WRITE(9,*)
      DO 10 Row=1,N
        WRITE(9,*) Bch,'(',Row,') =',A(Row,Col),'                     ',
     +Xch,'(',Row,') =',X(Row)
10    CONTINUE
      RETURN
      END

      SUBROUTINE RESULTS(A,X,N,DET)
      PARAMETER(MaxR=10,MaxC=11)
      INTEGER I,N
      REAL A,DET,X
      CHARACTER Bch*1,Xch*1
      DIMENSION A(1:MaxR,1:MaxC),X(1:MaxR)
      DO 10 I=1,15
        WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)'COMPUTATION OF THE SOLUTION FOR THE LINEAR SYSTEM A*X =
     + B.'
      WRITE(9,*)' '
      WRITE(9,*)'GAUSSIAN ELIMINATION WAS USED TO UPPER-TRIANGULARIZE,'
      WRITE(9,*)' '
      WRITE(9,*)'AND THEN BACK-SUBSTITUTION IS USED TO SOLVE FOR  X.'
      WRITE(9,*)' '
      WRITE(9,*)'THE COEFFICIENT MATRIX  A  IS:'
      WRITE(9,*)' '
      CALL MATPRINT(A,N)
      IF (DET.EQ.0) THEN
        WRITE(9,*)
        WRITE(9,*)'THE MATRIX IS SINGULAR.'
        WRITE(9,*)
        WRITE(9,*)'A ZERO PIVOT ELEMENT WAS ENCOUNTERED.'
        WRITE(9,*)' '
        WRITE(9,*)'THE METHOD DOES NOT APPLY.'
      ELSE
        WRITE(9,*)
        CALL VECPRINT('B','X',A,X,N)
      ENDIF
      WRITE(9,*)' '
      WRITE(9,*)'THE DETERMINANT`S VALUE IS  DET A = ',DET
      RETURN
      END

      PROGRAM LUFACTOR
C     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
C     To accompany the text:
C     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
C     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
C     This free software is complements of the author.
C
C     Algorithm 3.3 (PA = LU Factorization with Pivoting).
C     Section 3.6, Triangular Factorizaton, Page 175
C
      PARAMETER(MaxR=10)
      INTEGER InRC,N,Row
      REAL A,A1,B,DET,X
      CHARACTER Ans*1,Ach*1,Bch*1,Xch*1
      DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR)
      DIMENSION B(1:MaxR),Row(1:MaxR),X(1:MaxR)
10    CALL INPUTS(A,A1,B,N,InRC)
      CALL FACTOR(A,Row,N,DET)
20    WRITE(9,*)' '
      WRITE(9,*)'ENTER THE COLUMN VECTOR  B:'
      CALL VECIN('B',B,N)
      IF (DET.NE.0) THEN
        CALL SOLVE(A,Row,B,X,N)
      ENDIF
        CALL RESULTS(A1,B,X,N,DET)
      IF (DET.NE.0) THEN
        WRITE(9,*)'WANT TO SOLVE A*X=B WITH A NEW VECTOR B ? <Y/N> '
        READ(9,'(A)') Ans
        IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') GOTO 20
      ENDIF
      WRITE(9,*)'WANT TO SOLVE ANOTHER LINEAR SYSTEM ? <Y/N> '
      READ(9,'(A)') Ans
      IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') GOTO 10
      END

      SUBROUTINE FACTOR(A,Row,N,DET)
      PARAMETER(MaxR=10)
      INTEGER C,J,InRC,K,N,P,Row,RowK,RowP,T
      REAL A,A1,B,DET,X
      CHARACTER Ans*1,Ach*1,Bch*1,Xch*1
      DIMENSION A(1:MaxR,1:MaxR),Row(1:MaxR)
      DET=1
      DO J=1,N
        Row(J)=J
      ENDDO
      DO P=1,(N-1)
        DO K=(P+1),N
          IF (ABS(A(Row(K),P)).GT.ABS(A(Row(P),P))) THEN
            T=Row(P)
            Row(P)=Row(K)
            Row(K)=T
            DET=-DET
          ENDIF
        ENDDO
        DET=DET*A(Row(P),P)
        IF (DET.EQ.0) RETURN
        DO K=(P+1),N
          RowK=Row(K)
          RowP=Row(P)
          A(RowK,P)=A(RowK,P)/A(RowP,P)
          DO C=(P+1),N
            A(RowK,C)=A(RowK,C)-A(RowK,P)*A(RowP,C)
          ENDDO
        ENDDO
      ENDDO
      DET=DET*A(Row(N),N)
      RETURN
      END

      SUBROUTINE SOLVE(A,Row,B,X,N)
      PARAMETER(MaxR=10)
      INTEGER C,K,N,Row,RowK
      REAL A,B,Sum,X
      DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),Row(1:MaxR),X(1:MaxR)
      DO K=1,N
        IF (A(Row(K),K).EQ.0) RETURN
      ENDDO
      X(1)=B(Row(1))
      DO K=2,N
        Sum=0
        RowK=Row(K)
        DO C=1,(K-1)
          Sum=Sum+A(RowK,C)*X(C)
        ENDDO
        X(K)=B(RowK)-Sum
      ENDDO
      X(N)=X(N)/A(Row(N),N)
      DO K=(N-1),1,-1
        Sum=0
        RowK=Row(K)
        DO C=(K+1),N
          Sum=Sum+A(RowK,C)*X(C)
        ENDDO
        X(K)=(X(K)-Sum)/A(RowK,K)
      ENDDO
      RETURN
      END

      SUBROUTINE XFACTOR(A,Row,N,DET)
C     This subroutine uses simulated WHILE loop(s).
      PARAMETER(MaxR=10)
      INTEGER C,J,InRC,K,N,P,Row,RowK,RowP,T
      REAL A,A1,B,DET,X
      CHARACTER Ans*1,Ach*1,Bch*1,Xch*1
      DIMENSION A(1:MaxR,1:MaxR),Row(1:MaxR)
      DET=1
      DO 10 J=1,N
        Row(J)=J
10    CONTINUE
      DO 50 P=1,(N-1)
        DO 20 K=(P+1),N
        IF (ABS(A(Row(K),P)).GT.ABS(A(Row(P),P))) THEN
          T=Row(P)
          Row(P)=Row(K)
          Row(K)=T
          DET=-DET
        ENDIF
20      CONTINUE
        DET=DET*A(Row(P),P)
        IF (DET.EQ.0) RETURN
        DO 40 K=(P+1),N
          RowK=Row(K)
          RowP=Row(P)
          A(RowK,P)=A(RowK,P)/A(RowP,P)
          DO 30 C=(P+1),N
            A(RowK,C)=A(RowK,C)-A(RowK,P)*A(RowP,C)
30        CONTINUE
40      CONTINUE
50    CONTINUE
      DET=DET*A(Row(N),N)
      RETURN
      END

      SUBROUTINE XSOLVE(A,Row,B,X,N)
C     This subroutine uses simulated WHILE loop(s).
      PARAMETER(MaxR=10)
      INTEGER C,K,N,Row,RowK
      REAL A,B,Sum,X
      DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),Row(1:MaxR),X(1:MaxR)
      DO 10 K=1,N
        IF (A(Row(K),K).EQ.0) RETURN
10    CONTINUE
      X(1)=B(Row(1))
      DO 30 K=2,N
        Sum=0
        RowK=Row(K)
        DO 20 C=1,(K-1)
          Sum=Sum+A(RowK,C)*X(C)
20      CONTINUE
        X(K)=B(RowK)-Sum
30    CONTINUE
      X(N)=X(N)/A(Row(N),N)
      DO 50 K=(N-1),1,-1
        Sum=0
        RowK=Row(K)
        DO 40 C=(K+1),N
          Sum=Sum+A(RowK,C)*X(C)
40      CONTINUE
        X(K)=(X(K)-Sum)/A(RowK,K)
50    CONTINUE
      RETURN
      END

      SUBROUTINE INPUTS(A,A1,B,N,InRC)
      PARAMETER(MaxR=10)
      INTEGER Col,I,InRC,N,Row
      REAL A,B
      CHARACTER Ans*1,Ach*1
      DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR),B(1:MaxR)
      DO 10 I=1,12
        WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)'SOLUTION OF A LINEAR SYSTEM:'
      WRITE(9,*)' '
      WRITE(9,*)'      A*X = B'
      WRITE(9,*)' '
      WRITE(9,*)'THE TRIANGULAR FACTORIZATION L*U = P*A IS CONSTRUCTED'
      WRITE(9,*)' '
      WRITE(9,*)'FIRST,  THE SOLUTION  Y  TO  L*Y = P*B  IS FOUND,'
      WRITE(9,*)' '
      WRITE(9,*)'SECOND, THE SOLUTION  X  TO  U*X = Y    IS FOUND.'
      WRITE(9,*)' '
      WRITE(9,*)'A  IS AN  N BY N  NONSINGULAR MATRIX.'
      WRITE(9,*)' '
      WRITE(9,*)'B  IS AN  N  DIMENSIONAL VECTOR OF CONSTANTS.'
      WRITE(9,*)' '
      WRITE(9,*)'X  IS THE N DIMENSIONAL SOLUTION VECTOR.'
      WRITE(9,*)' '
      WRITE(9,*)'ENTER THE NUMBER OF EQUATIONS N = '
      READ(9,*) N
      InRC=0
      WRITE(9,*)'WANT TO INPUT COLUMNS ? <Y/N> '
      READ(9,'(A)') Ans
      IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') THEN
        InRC=1
      ENDIF
      WRITE(9,*)' '
      WRITE(9,*)'ENTER THE COEFFICIENT MATRIX  A:'
      CALL MATIN('A',A,N,InRC)
      DO 30 Row=1,N
      DO 20 Col=1,N
        A1(Row,Col)=A(Row,Col)
20    CONTINUE
30    CONTINUE
      RETURN
      END

      SUBROUTINE MATIN(Ach,A,N,InRC)
      PARAMETER(MaxR=10)
      INTEGER Col,InRC,N,Row
      REAL A
      CHARACTER Ach*1
      DIMENSION A(1:MaxR,1:MaxR)
      DO 20 Row=1,N
      DO 10 Col=1,N
        A(Row,Col)=0
10    CONTINUE
20    CONTINUE
      WRITE(9,*)'GIVE THE ELEMENTS OF THE MATRIX ',Ach
      WRITE(9,*)' '
      IF (InRC.EQ.0) THEN
        DO 40 Row=1,N
          WRITE(9,*)'GIVE THE ELEMENTS OF ROW ',Row
          WRITE(9,*)' '
        DO 30 Col=1,N
          WRITE(9,*) Ach,'(',Row,',',Col,') = '
          READ(9,*) A(Row,Col)
30      CONTINUE
40      CONTINUE
      ELSE
        DO 60 Col=1,N
          WRITE(9,*)'GIVE THE ELEMENTS OF COLUMN ',Col
          WRITE(9,*)' '
        DO 50 Row=1,N
          WRITE(9,*) Ach,'(',Row,',',Col,') = '
          READ(9,*) A(Row,Col)
50      CONTINUE
60      CONTINUE
      ENDIF
      RETURN
      END

      SUBROUTINE VECIN(Bch,B,N)
      PARAMETER(MaxR=10)
      INTEGER N,Row
      REAL A
      CHARACTER Bch*1
      DIMENSION B(1:MaxR)
      WRITE(9,*)'GIVE THE ELEMENTS OF THE COLUMN VECTOR ',Bch
      WRITE(9,*)' '
      DO 10 Row=1,N
        WRITE(9,*) Bch,'(',Row,') = '
        READ(9,*) B(Row)
10    CONTINUE
      RETURN
      END

      SUBROUTINE MATPRINT(A,N)
      PARAMETER(MaxR=10)
      INTEGER Col,N,Row
      REAL A
      DIMENSION A(1:MaxR,1:MaxR)
      DO 10 Row=1,N
        WRITE(9,999) ( A(Row,Col), Col=1,N )
10    CONTINUE
999   FORMAT(1X,5F15.7,4X)
      RETURN
      END

      SUBROUTINE VECPRINT(Bch,Xch,B,X,N)
      PARAMETER(MaxR=10)
      INTEGER N,Row
      REAL B,X
      CHARACTER Bch*1,Xch*1
      DIMENSION B(1:MaxR),X(1:MaxR)
      WRITE(9,*) Bch,' COEFFICIENT VECTOR                ',Xch,
     +' SOLUTION VECTOR'
      WRITE(9,*)
      DO 10 Row=1,N
        WRITE(9,*) Bch,'(',Row,') =',B(Row),'                     ',
     +Xch,'(',Row,') =',X(Row)
10    CONTINUE
      RETURN
      END

      SUBROUTINE RESULTS(A,B,X,N,DET)
      PARAMETER(MaxR=10)
      INTEGER I,N
      REAL A,B,DET,X
      CHARACTER Bch*1,Xch*1
      DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),X(1:MaxR)
      DO 10 I=1,15
        WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)'COMPUTATION OF THE SOLUTION FOR THE LINEAR SYSTEM A*X =
     + B.'
      IF (DET.NE.0) THEN
      WRITE(9,*)'THE TRIANGULAR FACTORIZATION L*U = P*A WAS CONSTRUCTED'
      WRITE(9,*)'FIRST,  THE SOLUTION  Y  TO  L*Y = P*B WAS FOUND,'
      WRITE(9,*)'SECOND, THE SOLUTION  X  TO  U*X = Y   WAS FOUND.'
      ENDIF
      WRITE(9,*)
      CALL MATPRINT(A,N)
      IF (DET.EQ.0) THEN
        WRITE(9,*)
        WRITE(9,*)'THE MATRIX IS SINGULAR.'
        WRITE(9,*)
        WRITE(9,*)'A ZERO PIVOT ELEMENT WAS ENCOUNTERED.'
        WRITE(9,*)' '
        WRITE(9,*)'THE MATRIX DOES NOT HAVE TRIANGULAR FACTORIZATION.'
        WRITE(9,*)' '
        WRITE(9,*)'THE METHOD DOES NOT APPLY.'
      ELSE
        WRITE(9,*)
        CALL VECPRINT('B','X',A,X,N)
      ENDIF
      WRITE(9,*)' '
      WRITE(9,*)'THE DETERMINANT`S VALUE IS  DET A = ',DET
      RETURN
      END

      PROGRAM JACOBI
C     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
C     To accompany the text:
C     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
C     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
C     This free software is complements of the author.
C
C     Algorithm 3.4 (Jacobi Iteration).
C     Section 3.7, Iterative Methods for Linear Systems, Page 186
C
      PARAMETER(MaxR=25,Max=99,Tol=1E-6)
      INTEGER Cond,Count,InRC,N,Row
      REAL A,A1,B,DET,Sep,X
      CHARACTER Ans*1,Ach*1,Bch*1,Xch*1
      DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR)
      DIMENSION B(1:MaxR),Row(1:MaxR),X(1:MaxR)
10    CALL INPUTS(A,A1,B,N,InRC)
      WRITE(9,*)'ENTER the starting vector.'
      Xch='X'
      CALL VECIN(Xch,X,N)
      CALL ITERATE(A,B,X,N,Max,Tol,Sep,Cond,Count)
      CALL RESULTS(A,B,X,N,Tol,Sep,Cond,Count)
      WRITE(9,*)' '
      WRITE(9,*)'Want to solve another linear system ? <Y/N> '
      READ(9,'(A)') Ans
      IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') GOTO 10
      END

      SUBROUTINE ITERATE(A,B,X,N,Max,Tol,Sep,Cond,Count)
      PARAMETER(MaxR=25)
      INTEGER C,Cond,Count,J,K,N,Max,R
      REAL A,A1,B,DET,X,Xold,Row,Sep,Sum,T,Tol
      DIMENSION A(1:MaxR,1:MaxR)
      DIMENSION B(1:MaxR),X(1:MaxR),Xold(1:MaxR)
      Sep=1
      K=0
      Cond=1
      DO J=1,N
        Xold(J)=X(J)
      ENDDO
      DO R=1,N
        Row=0
        DO C=1,N
          Row=Row+ABS(A(R,C))
        ENDDO
        IF (Row.GE.2*ABS(A(R,R))) Cond=0
      ENDDO
      IF (Cond.EQ.0) RETURN
      WRITE(9,1000) K, ( Xold(J), J=1,N )
      WHILE ((K.LT.Max).AND.(Sep.GT.Tol))
        DO R=1,N
          Sum=B(R)
          DO C=1,N
            IF (C.NE.R) Sum=Sum-A(R,C)*Xold(C)
          ENDDO
          X(R)=Sum/A(R,R)
        ENDDO
        Sep=0
        DO J=1,N
          Sep=Sep+ABS(X(J)-Xold(J))
          Xold(J)=X(J)
        ENDDO
        K=K+1
        WRITE(9,1000) K, ( Xold(J), J=1,N )
      REPEAT
      Count=K
1000  FORMAT(I2,4X,5F15.7)
      PAUSE
      RETURN
      END

      SUBROUTINE XITERATE(A,B,X,N,Max,Tol,Sep,Cond,Count)
C     This subroutine uses simulated WHILE loop(s).
      PARAMETER(MaxR=25)
      INTEGER C,Cond,Count,J,K,N,Max,R
      REAL A,A1,B,DET,X,Xold,Row,Sep,Sum,T,Tol
      CHARACTER Ans*1,Ach*1,Bch*1,Xch*1
      DIMENSION A(1:MaxR,1:MaxR)
      DIMENSION B(1:MaxR),X(1:MaxR),Xold(1:MaxR)
      Sep=1
      K=0
      Cond=1
      WRITE(9,*)'ENTER the starting vector.'
      Xch='X'
      CALL VECIN(Xch,X,N)
      DO 10 J=1,N
        Xold(J)=X(J)
10    CONTINUE
      DO 30 R=1,N
        Row=0
        DO 20 C=1,N
          Row=Row+ABS(A(R,C))
20      CONTINUE
        IF (Row.GE.2*ABS(A(R,R))) Cond=0
30    CONTINUE
      IF (Cond.EQ.0) RETURN
      WRITE(9,1000) K, ( Xold(J), J=1,N )
50    IF ((K.LT.Max).AND.(Sep.GT.Tol)) THEN
        DO 70 R=1,N
          Sum=B(R)
          DO 60 C=1,N
            IF (C.NE.R) Sum=Sum-A(R,C)*Xold(C)
60        CONTINUE
          X(R)=Sum/A(R,R)
70      CONTINUE
        Sep=0
        DO 80 J=1,N
          Sep=Sep+ABS(X(J)-Xold(J))
80      CONTINUE
        DO 90 J=1,N
          Xold(J)=X(J)
90      CONTINUE
        K=K+1
        WRITE(9,1000) K, ( Xold(J), J=1,N )
        GOTO 50
      ENDIF
      Count=K
1000  FORMAT(I2,4X,5F15.7)
      PAUSE
      RETURN
      END

      SUBROUTINE INPUTS(A,A1,B,N,InRC)
      PARAMETER(MaxR=25)
      INTEGER Col,I,InRC,N,Row
      REAL A,B
      CHARACTER Ans*1,Ach*1,Bch*1
      DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR),B(1:MaxR)
      DO 10 I=1,12
        WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)' '
      WRITE(9,*)'       Jacobi iteration is is performed to find t
     +he solution to the linear'
      WRITE(9,*)' '
      WRITE(9,*)'     system  A*X = B.   The matrix must be diagonally d
     +ominant.   Starting with'
      WRITE(9,*)' '
      WRITE(9,*)'                                             (k)  (k)
     +      (k)    '
      WRITE(9,*)'     X  = (0,0,...,0),  a sequence  { X  = (x   ,x   ,.
     +..,x   ) }  is generated'
      WRITE(9,*)'      0                                k     1    2
     +        N      '
      WRITE(9,*)' '
      WRITE(9,*)'     which converges to the solution.'
      WRITE(9,*)' '
      WRITE(9,*)'                           (k+1)'
      WRITE(9,*)'     the j-th coordinate  x       of  X     is computed
     + using the formula:'
      WRITE(9,*)'                           j           k+1'
      WRITE(9,*)' '
      WRITE(9,*)'      (k+1)                (k+1)              (k+1)'
      WRITE(9,*)'     x       =  [b  - a   x      -...- a     x     '
      WRITE(9,*)'      j           j    j,1 1            j,j-1 j-1  '
      WRITE(9,*)' '
      WRITE(9,*)'                             (k)            (k)      '
      WRITE(9,*)'                    - a     x    -...- a   x   ]/a
     +       for  j=1,2,...,N. '
      WRITE(9,*)'                       j,j+1 j+1        j,N N     j,j'
      WRITE(9,*)' '
      WRITE(9,*)' '
      WRITE(9,*)'     give the dimension  N = '
      READ(9,*) N
      InRC=0
      WRITE(9,*)'Want to input columns ? <Y/N> '
      READ(9,'(A)') Ans
      IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') THEN
        InRC=1
      ENDIF
      WRITE(9,*)' '
      WRITE(9,*)'ENTER the coefficient matrix  A:'
      Ach='A'
      CALL MATIN(Ach,A,N,InRC)
      WRITE(9,*)'The matrix  A  is:'
      CALL MATPRINT(A,N)
      Bch='B'
      CALL VECIN(Bch,B,N)
      RETURN
      END

      SUBROUTINE MATIN(Ach,A,N,InRC)
      PARAMETER(MaxR=25)
      INTEGER Col,InRC,N,Row
      REAL A
      CHARACTER Ach*1
      DIMENSION A(1:MaxR,1:MaxR)
      DO 20 Row=1,N
      DO 10 Col=1,N
        A(Row,Col)=0
10    CONTINUE
20    CONTINUE
      WRITE(9,*)'Give the elements of the matrix ',Ach
      WRITE(9,*)' '
      IF (InRC.EQ.0) THEN
        DO 40 Row=1,N
          WRITE(9,*)'ENTER the elements of row ',Row
          WRITE(9,*)' '
        DO 30 Col=1,N
          WRITE(9,*) Ach,'(',Row,',',Col,') = '
          READ(9,*) A(Row,Col)
30      CONTINUE
40      CONTINUE
      ELSE
        DO 60 Col=1,N
          WRITE(9,*)'ENTER the elements of column ',Col
          WRITE(9,*)' '
        DO 50 Row=1,N
          WRITE(9,*) Ach,'(',Row,',',Col,') = '
          READ(9,*) A(Row,Col)
50      CONTINUE
60      CONTINUE
      ENDIF
      RETURN
      END

      SUBROUTINE VECIN(Bch,B,N)
      PARAMETER(MaxR=25)
      INTEGER N,Row
      REAL B
      CHARACTER Bch*1
      DIMENSION B(1:MaxR)
      WRITE(9,*)'ENTER the elements of the column vector ',Bch
      WRITE(9,*)' '
      DO 10 Row=1,N
        WRITE(9,*) Bch,'(',Row,') = '
        READ(9,*) B(Row)
10    CONTINUE
      RETURN
      END

      SUBROUTINE MATPRINT(A,N)
      PARAMETER(MaxR=25)
      INTEGER Col,N,Row
      REAL A
      DIMENSION A(1:MaxR,1:MaxR)
      DO 10 Row=1,N
        WRITE(9,999) ( A(Row,Col), Col=1,N )
10    CONTINUE
999   FORMAT(1X,5F15.7,4X)
      RETURN
      END

      SUBROUTINE VECPRINT(Bch,Xch,B,X,N)
      PARAMETER(MaxR=25)
      INTEGER N,Row
      REAL B,X
      CHARACTER Bch*1,Xch*1
      DIMENSION B(1:MaxR),X(1:MaxR)
      WRITE(9,*) Bch,' COEFFICIENT VECTOR                ',Xch,
     +' SOLUTION VECTOR'
      WRITE(9,*)
      DO 10 Row=1,N
        WRITE(9,*) Bch,'(',Row,') =',B(Row),'    ',Xch,'(',Row,') ='
     +,X(Row)
10    CONTINUE
      RETURN
      END

      SUBROUTINE RESULTS(A,B,X,N,Tol,Sep,Cond,Count)
      PARAMETER(MaxR=25)
      INTEGER C,Cond,Count,K,N,R
      REAL A,B,X,Tol,Sep
      CHARACTER Bch*1,Xch*1
      DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),X(1:MaxR)
      DO 10 I=1,15
        WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)' '
      WRITE(9,*)'The matrix A(I,J) is:'
      IF (N.LT.6) THEN
        WRITE(9,*)' '
      ENDIF
      CALL MATPRINT(A,N)
      WRITE(9,*)' '
      WRITE(9,*)'The vector B is:'
      IF (N.LT.6) THEN
        WRITE(9,*)' '
      ENDIF
      DO 20 K=1,N
        WRITE(9,*)'      B(',K,') = ',B(K)
20    CONTINUE
      WRITE(9,*)' '
      IF (Cond.EQ.1) THEN
        IF (Sep.LT.Tol) THEN
          WRITE(9,*)'   Jacobi iteration solved  A*X = B.'
          WRITE(9,*)'It took',Count,' iterations to find the solution:'
        ELSE
          WRITE(9,*)'    Jacobi iteration did not converge.'
          WRITE(9,*)'The status after ',Count,' iterations is'
        ENDIF
        IF (N.LT.6) THEN
          WRITE(9,*)' '
        ENDIF
        WRITE(9,*)' '
        DO 30 K=1,N
            WRITE(9,*)'      X(',K,') = ',X(K)
30      CONTINUE
        WRITE(9,*)' '
        WRITE(9,*)'           +-',Sep,'  is the estimated accuracy.'
      ELSE
        WRITE(9,*)' '
        WRITE(9,*)'The matrix is not diagonally dominant.'
        WRITE(9,*)' '
        WRITE(9,*)'Jacobi iteration cannot be used.'
      ENDIF
      RETURN
      END

      PROGRAM GAUSSIT
C     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
C     To accompany the text:
C     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
C     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
C     This free software is complements of the author.
C
C     Algorithm 3.5 (Gauss-Seidel Iteration).
C     Section 3.7, Iterative Methods for Linear Systems, Page 187
C
      PARAMETER(MaxR=25,Max=99,Tol=1E-7)
      INTEGER Cond,Count,InRC,N,Row
      REAL A,A1,B,DET,Sep,X
      CHARACTER Ans*1,Ach*1,Bch*1,Xch*1
      DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR)
      DIMENSION B(1:MaxR),Row(1:MaxR),X(1:MaxR)
10    CALL INPUTS(A,A1,B,N,InRC)
      WRITE(9,*)'ENTER the starting vector.'
      Xch='X'
      CALL VECIN(Xch,X,N)
      CALL ITERATE(A,B,X,N,Max,Tol,Sep,Cond,Count)
      CALL RESULTS(A,B,X,N,Tol,Sep,Cond,Count)
      WRITE(9,*)' '
      WRITE(9,*)'Want to solve another linear system ? <Y/N> '
      READ(9,'(A)') Ans
      IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') GOTO 10
      END

      SUBROUTINE ITERATE(A,B,X,N,Max,Tol,Sep,Cond,Count)
      PARAMETER(MaxR=25)
      INTEGER C,Cond,Count,J,K,N,Max,R
      REAL A,A1,B,DET,X,Xold,Row,Sep,Sum,T,Tol
      CHARACTER Ans*1,Ach*1,Bch*1,Xch*1
      DIMENSION A(1:MaxR,1:MaxR)
      DIMENSION B(1:MaxR),X(1:MaxR),Xold(1:MaxR)
      Sep=1
      K=0
      Cond=1
      DO J=1,N
        Xold(J)=X(J)
      ENDDO
      DO R=1,N
        Row=0
        DO C=1,N
          Row=Row+ABS(A(R,C))
        ENDDO
        IF (Row.GE.2*ABS(A(R,R))) Cond=0
      ENDDO
      IF (Cond.EQ.0) RETURN
      WRITE(9,1000) K, ( Xold(J), J=1,N )
      WHILE ((K.LT.Max).AND.(Sep.GT.Tol))
        DO R=1,N
          Sum=B(R)
          DO C=1,N
            IF (C.NE.R) Sum=Sum-A(R,C)*X(C)
          ENDDO
          X(R)=Sum/A(R,R)
        ENDDO
        Sep=0
        DO J=1,N
          Sep=Sep+ABS(X(J)-Xold(J))
          Xold(J)=X(J)
        ENDDO
        K=K+1
        WRITE(9,1000) K, ( Xold(J), J=1,N )
      REPEAT
      Count=K
1000  FORMAT(I2,4X,5F15.7)
      PAUSE
      RETURN
      END

      SUBROUTINE XITERATE(A,B,X,N,Max,Tol,Sep,Cond,Count)
C     This subroutine uses simulated WHILE loop(s).
      PARAMETER(MaxR=25)
      INTEGER C,Cond,Count,J,K,N,Max,R
      REAL A,A1,B,DET,X,Xold,Row,Sep,Sum,T,Tol
      CHARACTER Ans*1,Ach*1,Bch*1,Xch*1
      DIMENSION A(1:MaxR,1:MaxR)
      DIMENSION B(1:MaxR),X(1:MaxR),Xold(1:MaxR)
      Sep=1
      K=0
      Cond=1
      DO 10 J=1,N
        Xold(J)=X(J)
10    CONTINUE
      DO 30 R=1,N
        Row=0
        DO 20 C=1,N
          Row=Row+ABS(A(R,C))
20      CONTINUE
        IF (Row.GE.2*ABS(A(R,R))) Cond=0
30    CONTINUE
      IF (Cond.EQ.0) RETURN
      DO 40 J=1,10
      WRITE(9,*) ' '
40    CONTINUE
      WRITE(9,1000) K, ( Xold(J), J=1,N )
50    IF ((K.LT.Max).AND.(Sep.GT.Tol)) THEN
        DO 70 R=1,N
          Sum=B(R)
          DO 60 C=1,N
            IF (C.NE.R) Sum=Sum-A(R,C)*X(C)
60        CONTINUE
          X(R)=Sum/A(R,R)
70      CONTINUE
        Sep=0
        DO 80 J=1,N
          Sep=Sep+ABS(X(J)-Xold(J))
          Xold(J)=X(J)
80      CONTINUE
        K=K+1
        WRITE(9,1000) K, ( Xold(J), J=1,N )
        GOTO 50
      ENDIF
      Count=K
1000  FORMAT(I2,4X,5F15.7)
      PAUSE
      RETURN
      END

      SUBROUTINE INPUTS(A,A1,B,N,InRC)
      PARAMETER(MaxR=25)
      INTEGER Col,I,InRC,N,Row
      REAL A,B
      CHARACTER Ans*1,Ach*1,Bch*1
      DIMENSION A(1:MaxR,1:MaxR),A1(1:MaxR,1:MaxR),B(1:MaxR)
      DO 10 I=1,12
        WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)' '
      WRITE(9,*)'       Gauss-Seidel iteration is is performed to find t
     +he solution to the linear'
      WRITE(9,*)' '
      WRITE(9,*)'     system  A*X = B.   The matrix must be diagonally d
     +ominant.   Starting with'
      WRITE(9,*)' '
      WRITE(9,*)'                                             (k)  (k)
     +      (k)    '
      WRITE(9,*)'     X  = (0,0,...,0),  a sequence  { X  = (x   ,x   ,.
     +..,x   ) }  is generated'
      WRITE(9,*)'      0                                k     1    2
     +        N      '
      WRITE(9,*)' '
      WRITE(9,*)'     which converges to the solution.'
      WRITE(9,*)' '
      WRITE(9,*)'                           (k+1)'
      WRITE(9,*)'     the j-th coordinate  x       of  X     is computed
     + using the formula:'
      WRITE(9,*)'                           j           k+1'
      WRITE(9,*)' '
      WRITE(9,*)'      (k+1)                (k+1)              (k+1)'
      WRITE(9,*)'     x       =  [b  - a   x      -...- a     x     '
      WRITE(9,*)'      j           j    j,1 1            j,j-1 j-1  '
      WRITE(9,*)' '
      WRITE(9,*)'                             (k)            (k)      '
      WRITE(9,*)'                    - a     x    -...- a   x   ]/a
     +       for  j=1,2,...,N. '
      WRITE(9,*)'                       j,j+1 j+1        j,N N     j,j'
      WRITE(9,*)' '
      WRITE(9,*)' '
      WRITE(9,*)'     give the dimension  N = '
      READ(9,*) N
      InRC=0
      WRITE(9,*)'Want to input columns ? <Y/N> '
      READ(9,'(A)') Ans
      IF (Ans.EQ.'Y' .OR. Ans.EQ.'y') THEN
        InRC=1
      ENDIF
      WRITE(9,*)' '
      WRITE(9,*)'ENTER the coefficient matrix  A:'
      Ach='A'
      CALL MATIN(Ach,A,N,InRC)
      WRITE(9,*)'The matrix  A  is:'
      CALL MATPRINT(A,N)
      Bch='B'
      CALL VECIN(Bch,B,N)
      RETURN
      END

      SUBROUTINE MATIN(Ach,A,N,InRC)
      PARAMETER(MaxR=25)
      INTEGER Col,InRC,N,Row
      REAL A
      CHARACTER Ach*1
      DIMENSION A(1:MaxR,1:MaxR)
      DO 20 Row=1,N
      DO 10 Col=1,N
        A(Row,Col)=0
10    CONTINUE
20    CONTINUE
      WRITE(9,*)'Give the elements of the matrix ',Ach
      WRITE(9,*)' '
      IF (InRC.EQ.0) THEN
        DO 40 Row=1,N
          WRITE(9,*)'ENTER the elements of row ',Row
          WRITE(9,*)' '
        DO 30 Col=1,N
          WRITE(9,*) Ach,'(',Row,',',Col,') = '
          READ(9,*) A(Row,Col)
30      CONTINUE
40      CONTINUE
      ELSE
        DO 60 Col=1,N
          WRITE(9,*)'ENTER the elements of column ',Col
          WRITE(9,*)' '
        DO 50 Row=1,N
          WRITE(9,*) Ach,'(',Row,',',Col,') = '
          READ(9,*) A(Row,Col)
50      CONTINUE
60      CONTINUE
      ENDIF
      RETURN
      END

      SUBROUTINE VECIN(Bch,B,N)
      PARAMETER(MaxR=25)
      INTEGER N,Row
      REAL B
      CHARACTER Bch*1
      DIMENSION B(1:MaxR)
      WRITE(9,*)'ENTER the elements of the column vector ',Bch
      WRITE(9,*)' '
      DO 10 Row=1,N
        WRITE(9,*) Bch,'(',Row,') = '
        READ(9,*) B(Row)
10    CONTINUE
      RETURN
      END

      SUBROUTINE MATPRINT(A,N)
      PARAMETER(MaxR=25)
      INTEGER Col,N,Row
      REAL A
      DIMENSION A(1:MaxR,1:MaxR)
      DO 10 Row=1,N
        WRITE(9,999) ( A(Row,Col), Col=1,N )
10    CONTINUE
999   FORMAT(1X,5F15.7,4X)
      RETURN
      END

      SUBROUTINE VECPRINT(Bch,Xch,B,X,N)
      PARAMETER(MaxR=25)
      INTEGER N,Row
      REAL B,X
      CHARACTER Bch*1,Xch*1
      DIMENSION B(1:MaxR),X(1:MaxR)
      WRITE(9,*) Bch,' COEFFICIENT VECTOR                ',Xch,
     +' SOLUTION VECTOR'
      WRITE(9,*)
      DO 10 Row=1,N
        WRITE(9,*) Bch,'(',Row,') =',B(Row),'    ',Xch,'(',Row,') ='
     +,X(Row)
10    CONTINUE
      RETURN
      END

      SUBROUTINE RESULTS(A,B,X,N,Tol,Sep,Cond,Count)
      PARAMETER(MaxR=25)
      INTEGER C,Cond,Count,K,N,R
      REAL A,B,X,Tol,Sep
      CHARACTER Bch*1,Xch*1
      DIMENSION A(1:MaxR,1:MaxR),B(1:MaxR),X(1:MaxR)
      DO 10 I=1,15
        WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)' '
      WRITE(9,*)'The matrix A(I,J) is:'
      IF (N.LT.6) THEN
        WRITE(9,*)' '
      ENDIF
      CALL MATPRINT(A,N)
      WRITE(9,*)' '
      WRITE(9,*)'The vector B is:'
      IF (N.LT.6) THEN
        WRITE(9,*)' '
      ENDIF
      DO 20 K=1,N
        WRITE(9,*)'      B(',K,') = ',B(K)
20    CONTINUE
      WRITE(9,*)' '
      IF (Cond.EQ.1) THEN
        IF (Sep.LT.Tol) THEN
          WRITE(9,*)'Gauss-Seidel iteration solved  A*X = B.'
          WRITE(9,*)'It took',Count,' iterations to find the solution:'
        ELSE
          WRITE(9,*)'Gauss-Seidel iteration did not converge.'
          WRITE(9,*)'The status after ',Count,' iterations is'
        ENDIF
        IF (N.LT.6) THEN
          WRITE(9,*)' '
        ENDIF
        WRITE(9,*)' '
        DO 30 K=1,N
            WRITE(9,*)'      X(',K,') = ',X(K)
30      CONTINUE
        WRITE(9,*)' '
        WRITE(9,*)'           +-',Sep,'  is the estimated accuracy.'
      ELSE
        WRITE(9,*)' '
        WRITE(9,*)'The matrix is not diagonally dominant.'
        WRITE(9,*)' '
        WRITE(9,*)'Gauss-Seidel iteration cannot be used.'
      ENDIF
      RETURN
      END