
C     TITLE: TEST SUBROUTINE TRIDGL
C
C     THIS TESTS THE SUBROUTINE 'TRIDGL' GIVEN IN SECTION 8.4.
C     
      DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MAX_N=200)
      DIMENSION A(MAX_N), B(MAX_N), C(MAX_N), F(MAX_N)
C
C     INPUT THE SYSTEM TO BE SOLVED.
C
      PRINT *, ' WHAT IS ORDER N OF SYSTEM?'
      READ *, N
      PRINT *, ' GIVE THE EQUATION COEFFICIENTS.'
      PRINT *, ' GIVE THEM ONE EQUATION AT A TIME.'
      PRINT *, ' GIVE ONLY THE ELEMENTS OF THE '
      PRINT *, ' SUBDIAGONAL, DIAGONAL, AND '
      PRINT *, ' SUPERDIAGONAL, AS APPROPRIATE.'
      PRINT *, ' ALSO INCLUDE THE EQUATIONS RIGHT SIDE.'
      PRINT *, ' '
      PRINT *, ' GIVE B(1), C(1), RHS(1) FOR EQUATION', 1
      READ *, B(1), C(1), F(1)
      DO 10 I=2,N-1
        PRINT *, ' GIVE A(I), B(I), C(I), RHS(I) FOR EQUATION',I
10      READ *, A(I), B(I), C(I), F(I)
      PRINT *, ' GIVE A(N), B(N), RHS(N) FOR EQUATION',N
      READ *, A(N), B(N), F(N)
C
C     CALL SUBROUTINE TRIDGL TO SOLVE THE TRIDIAGONAL SYSTEM.
      IFLAG = 0
      CALL TRIDGL(A,B,C,F,N,IFLAG,IER)
C
C     OUTPUT THE SOLUTION OF THE SYSTEM.
      PRINT *, ' SOLUTIONS'
      PRINT *, ' I      X(I)'
      DO 20 J=1,N
20      PRINT *, J,F(J)
      STOP
      END

      SUBROUTINE TRIDGL(A,B,C,F,N,IFLAG,IER)
C     -----------------
C
C     THIS SOLVES A TRIDIAGONAL SYSTEM OF LINEAR EQUATIONS  M*X=F.
C
C     INPUT:
C     THE ORDER OF THE LINEAR SYSTEM IS GIVEN BY  N.
C     THE SUBDIAGONAL,DIAGONAL, AND SUPERDIAGONAL OF  M  ARE GIVEN
C     BY THE ARRAYS  A, B, AND C, RESPECTIVELY. MORE PRECISELY,
C        M(I,I-1) = A(I),  I=2,...,N
C        M(I,I)   = B(I),  I=1,...,N
C        M(I,I+1) = C(I),  I=1,...,N-1
C     IFLAG=0 MEANS THAT THE ORIGINAL MATRIX  M  IS GIVEN AS
C     SPECIFIED ABOVE.
C     IFLAG=1 MEANS THAT THE LU FACTORIZATION OF  M  IS ALREADY KNOWN
C     AND IS STORED IN  A,B, AND C. THIS WILL HAVE BEEN ACCOMPLISHED
C     BY A PREVIOUS CALL TO THIS SUBROUTINE.
C
C     OUTPUT:
C     UPON EXIT, THE LU FACTORIZATION OF  M  WILL BE STORED
C     IN  A,B, AND C.
C     THE SOLUTION VECTOR  X  WILL BE RETURNED IN  F.
C     IER=0 MEANS THE PROGRAM WAS COMPLETED SATISFACTORILY.
C     IER=1 MEANS THAT A ZERO PIVOT ELEMENT WAS ENCOUNTERED, AND
C     NO SOLUTION WAS ATTEMPTED.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (ZERO=0.0D0)
      DIMENSION A(*), B(*), C(*), F(*)
C
      IF(IFLAG .EQ. 0) THEN
C         COMPUTE LU FACTORIZATION OF MATRIX M.
          DO 10 J=2,N
            IF(B(J-1) .EQ. ZERO) THEN
                IER = 1
                RETURN
            END IF
            A(J) = A(J)/B(J-1)
10          B(J) = B(J) - A(J)*C(J-1)
          IF(B(N) .EQ. ZERO) THEN
              IER = 1
              RETURN
          END IF
      END IF
C
C     COMPUTE SOLUTION X TO M*X=F.
      DO 20 J=2,N
20      F(J) = F(J) - A(J)*F(J-1)
      F(N) = F(N)/B(N)
      DO 30 J=N-1,1,-1
30      F(J) = (F(J) - C(J)*F(J+1))/B(J)
      IER = 0
      RETURN
      END
