      SUBROUTINE FACTOR( N, D, E, IERR, RV1, RV2, RV3, RV4 )
C
C     MICHAEL B. PORTER 7/1/85
C
C     THIS SUBROUTINE IS BASED ON TINVIT IN EISPACK
C     USES GAUSSIAN ELIMINATION WITH PIVOTING TO FACTOR A
C        SYMMETRIC TRIDIAGONAL LINEAR SYSTEM
C     
C     ON INPUT
C
C        N IS THE ORDER OF THE MATRIX.
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
C
C
C     ON OUTPUT
C
C        RV1, RV2, RV3 AND RV4 CONTAIN THE FACTORED MATRIX
C
      COMPLEX U, V, XU, D(N), E(N), RV1(N), RV2(N), RV3(N), RV4(N)

C     LU DECOMPOSITION WITHOUT INTERCHANGES

      U = D(1)
      V = E(2)

      DO I = 2, N-1
         XU = E( I ) / U
         RV4( I   ) = XU
         RV1( I-1 ) = 1.0 / U
         RV2( I-1 ) = V
         U = D(I) - XU * V
         V = E(I+1)
      END DO

      XU       = E( N ) / U
      RV4(N)   = XU
      RV1(N-1) = 1.0 / U
      RV2(N-1) = V
      U        = D(N) - XU * V

      IF( U .EQ. 0.0 ) WRITE( *, * ) 'SINGULAR MATRIX'
      RV1(N) = 1.0 / U
      RV2(N) = 0.0

      RETURN
      END

      SUBROUTINE BACKSB(N,D,E,IERR,RV1,RV2,RV3,RV4,B)
C
C     MICHAEL B. PORTER 7/1/85
C
C     THIS SUBROUTINE IS BASED ON TINVIT IN EISPACK
C     PERFORMS BACK-SUBSTITUTION FOR A
C        SYMMETRIC TRIDIAGONAL LINEAR SYSTEM
C     
C     ON INPUT
C
C        N IS THE ORDER OF THE MATRIX.
C        D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX.
C        E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
C          IN ITS LAST N-1 POSITIONS.  E(1) IS ARBITRARY.
C        B CONTAINS THE RIGHT HAND SIDE (AX=B)
C
C     ON OUTPUT
C
C        ALL INPUT ARRAYS ARE UNALTERED EXCEPT
C        B WHICH CONTAINS THE SOLUTION
C
C     RV1, RV2, RV3 AND RV4 ARE TEMPORARY STORAGE ARRAYS.

      COMPLEX D(N),E(N),RV1(N),RV2(N),RV3(N),RV4(N),B(N)
      DATA SMALL / 1.0E-37 /   ! smallest pos. machine constant

C     *** FORWARD ELIMINATION ***

      DO I = 2, N

         B( I ) = B( I ) - RV4( I ) * B( I - 1 )

C        --- prevent underflows from propagating:
         IF ( ABS(  REAL( B( I ) ) ) .LE. SMALL .AND.
     &        ABS( AIMAG( B( I ) ) ) .LE. SMALL ) B( I ) = 0.0
      END DO

C     *** BACK SUBSTITUTION (RESULT IN B) ***

      B( N ) = B( N ) * RV1( N )

      IF ( N .GE. 2 ) THEN
      DO I = N - 1, 1, -1
         B( I ) = ( B( I ) - B( I + 1 ) * RV2( I ) ) * RV1( I )
      END DO
      END IF

      RETURN
      END
