      SUBROUTINE FACTOR( N, A, IERR, RV, RVU )
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        A CONTAINS THE INPUT MATRIX.
C
C
C     ON OUTPUT
C
C        RV1, RV2, RV3 AND RV4 CONTAIN THE FACTORED MATRIX
C
      COMPLEX U, XU, A(2, N), RV(2, N), RVU(2, N)

C     LU DECOMPOSITION WITHOUT INTERCHANGES

      XU = 1.0
      U = A(1,1)

      IF ( N .GT. 1 ) THEN
      DO 580 I = 2, N-1
         XU = A( 2, I ) / U
         RV( 2, I   ) = XU
         RVU( 1, I-1 ) = 1.0 / U
         RVU( 2, I-1 ) = A( 2, I )
         U = A( 1, I   ) - XU * A( 2, I )
  580 CONTINUE
      ENDIF

      XU = A( 2, N ) / U
      RV( 2, N   ) = XU
      RVU( 1, N-1 ) = 1.0 / U
      RVU( 2, N-1 ) = A( 2, N )
      U = A( 1, N ) - XU * A( 2, N )

      IF( U .EQ. 0.0 ) WRITE( *, * ) 'SINGULAR MATRIX'
      RVU(1,N) = 1.0 / U
      RVU(2,N) = 0.0

      RETURN
      END

      SUBROUTINE BACKSB(N,A,IERR,RV,RVU,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     RV IS A TEMPORARY STORAGE ARRAY.

      COMPLEX U,V,A(2,N),RV(2,N), RVU(2,N), B(N)

C     *** FORWARD ELIMINATION ***

      DO 820 I = 2, N
         B( I ) = B( I ) - RV( 2, I ) * B( I-1 )

  820 CONTINUE

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

      B( N ) = B( N ) * RV( 1, N )

      IF ( N .GT. 1 ) THEN
      V = RVU( 2, N-1 ) * B( N )
      DO 620 I = N-1, 1, -1


C         B( I ) = ( B( I ) - B( I + 1 ) * RVU( 2, I ) ) * RVU( 1, I )
         B( I ) = ( B( I ) - B( I + 1 ) * RVU( 2, I ) )
C         B( I ) = ( B( I ) - V ) * RVU( 1, I )
C         B( I ) = ( RVU( 2, I ) - V ) * RVU( 1, I )
C         V = RVU( 2, I - 1 ) * B( I )

  620 CONTINUE
      ENDIF

      RETURN
      END
