C
C      ________________________________________________________
C     |                                                        |
C     |     COMPUTE EIGENVECTOR CORRESPONDING TO GIVEN REAL    |
C     |        EIGENVALUE FOR A REAL TRIDIAGONAL MATRIX        |
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         E     --EIGENVALUE                             |
C     |                                                        |
C     |         L     --SUBDIAGONAL                            |
C     |                                                        |
C     |         D     --DIAGONAL                               |
C     |                                                        |
C     |         U     --SUPERDIAGONAL                          |
C     |                                                        |
C     |         N     --MATRIX DIMENSION                       |
C     |                                                        |
C     |         W     --WORK ARRAY (LENGTH AT LEAST 4N)        |
C     |                                                        |
C     |    OUTPUT:                                             |
C     |                                                        |
C     |         E     --IMPROVED ESTIMATE FOR EIGENVALUE       |
C     |                                                        |
C     |         X     --EIGENVECTOR                            |
C     |                                                        |
C     |    BUILTIN FUNCTIONS: ABS,SQRT                         |
C     |________________________________________________________|
C
      SUBROUTINE TVECT(E,X,L,D,U,N,W)
      REAL D(1),L(1),U(1),W(1),X(1),E,O,P,Q,R,S,T,V,Y,Z
      INTEGER F,G,H,I,J,K,M,N
      IF ( N .GT. 1 ) GOTO 10
      E = D(1)
      X(1) = 1.
      RETURN
10    M = N - 1
      J = 2
C     ---------------------------------------------------------
C     |*** STORE MATRIX IN W ARRAY AND SUBTRACT EIGENVALUE ***|
C     ---------------------------------------------------------
      DO 20 I = 1,M
           X(I) = 0.
           W(J) = U(I)
           W(J+1) = D(I) - E
           W(J+2) = L(I)
20         J = J + 4
      W(J) = 0.
      W(J+1) = D(N) - E
      O = 65536.**(-4)
      T = .5*O
      S = T
30    T = .5*T
      P = S
      S = S + T
      IF ( S .GE. O ) GOTO 40
      IF ( S+T .GT. S ) GOTO 30
40    R = W(3)
      V = ABS(R) + ABS(W(4))
      F = 4
      M = 4*N - 3
      G = -2
C     ---------------------------
C     |*** FACTOR THE MATRIX ***|
C     ---------------------------
50    G = G + 4
      IF ( G .GT. M ) GOTO 110
      H = G - 1
      I = G + 2
      J = G + 5
C     --------------------
C     |*** FIND PIVOT ***|
C     --------------------
      Q = W(I)
      Y = ABS(Q)
      Z = ABS(R)
      IF ( Z .GE. Y ) GOTO 80
C     -------------------
C     |*** SWAP ROWS ***|
C     -------------------
      IF ( V .GE. Y ) GOTO 60
      V = Y
      F = I
60    T = W(G)
      W(G) = W(J)
      W(J) = T
      T = R/Q
      K = G + 1
      W(K) = Q
      K = J - 1
      S = W(K)
      IF ( S .EQ. 0. ) GOTO 70
      IF ( S .EQ. O ) S = P
      W(K) = -S*T
      W(H) = S
      GOTO 100
70    W(K) = S
      W(H) = O
      GOTO 100
80    W(H) = 0.
      IF ( V .GE. Z ) GOTO 90
      V = Z
      F = I
90    IF ( R .EQ. 0. ) GOTO 120
      T = Q/R
C     -------------------
C     |*** ELIMINATE ***|
C     -------------------
100   R = W(J) - T*W(G)
      W(J) = R
      W(I) = T
      GOTO 50
110   IF ( ABS(R) .GE. V ) GOTO 120
      V = R
      F = J + 1
C     ---------------------------------------------------
C     |*** COMPUTE INITIAL EIGENVECTOR APPROXIMATION ***|
C     ---------------------------------------------------
120   J = F/4
      X(J) = 1.
      IF ( J .EQ. 1 ) GOTO 140
      K = F - 5
      J = J - 1
      X(J) = (X(J)-W(K-1)*X(J+1))/W(K)
130   IF ( J .EQ. 1 ) GOTO 140
      J = J - 1
      K = K - 4
      T = W(K-2)
      IF ( T .EQ. O ) T = 0.
      X(J) = (X(J)-W(K-1)*X(J+1)-T*X(J+2))/W(K)
      GOTO 130
140   IF ( V .EQ. 0. ) GOTO 230
      S = 0.
      DO 150 I = 1,N
           T = ABS(X(I))
150        IF ( T .GT. S ) S = T
      R = 0.
      S = 1./S
      DO 160 I = 1,N
           T = S*X(I)
           R = R + T*T
160        X(I) = T
      K = 0
      J = 1
      Y = X(1)
C     -----------------------------------------------------
C     |*** APPLY ONE ITERATION OF INVERSE POWER METHOD ***|
C     -----------------------------------------------------
170   K = K + 4
      I = J
      J = J + 1
      S = W(K)
      W(K) = Y
      Y = X(J)
      IF ( W(K-3) .EQ. 0. ) GOTO 180
      T = X(J)
      X(J) = X(I)
      X(I) = T
180   X(J) = X(J) - S*X(I)
      IF ( J .LT. N ) GOTO 170
C     ---------------------------
C     |*** BACK SUBSTITUTION ***|
C     ---------------------------
      S = X(J)/W(K+3)
      X(J) = S
      T = ABS(S)
      V = S*Y
      J = J - 1
      K = K - 1
      S = (X(J)-W(K-1)*S)/W(K)
      X(J) = S
      IF ( ABS(S) .GT. T ) T = ABS(S)
      V = V + S*W(K+1)
190   IF ( J .EQ. 1 ) GOTO 200
      J = J - 1
      K = K - 4
      Z = W(K-2)
      IF ( Z .EQ. O ) Z = 0.
      S = (X(J)-W(K-1)*S-Z*X(J+2))/W(K)
      X(J) = S
      IF ( ABS(S) .GT. T ) T = ABS(S)
      V = V + S*W(K+1)
      GOTO 190
200   IF ( V .NE. 0. ) V = R/V
      S = 0.
      T = 1./T
      Z = 0.
      DO 210 I = 1,N
           S = S + (X(I)*V)**2
210        Z = Z + (T*X(I))**2
      T = T/SQRT(Z)
      DO 220 I = 1,N
220        X(I) = T*X(I)
C     --------------------------------------------------------------
C     |*** USE RAYLEIGH QUOTIENT TO IMPROVE EIGENVALUE ESTIMATE ***|
C     --------------------------------------------------------------
      IF ( R+R .GE. S ) E = E + V
      RETURN
230   T = 0.
      DO 240 I = 1,N
           S = ABS(X(I))
240        IF ( S .GT. T ) T = S
      T = 1./T
      Z = 0.
      DO 250 I = 1,N
250        Z = Z + (T*X(I))**2
      T = T/SQRT(Z)
      DO 260 I = 1,N
260        X(I) = T*X(I)
      RETURN
      END