C
C      ________________________________________________________
C     |                                                        |
C     |     COMPUTE EIGENVECTOR CORRESPONDING TO GIVEN REAL    |
C     |         EIGENVALUE FOR A REAL HESSENBERG MATRIX        |
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         E     --ESTIMATE OF EIGENVALUE                 |
C     |                                                        |
C     |         A     --HESSENBERG MATRIX PACKED AT START OF   |
C     |                 ARRAY WITH LENGTH AT LEAST N(N+7)/2 - 2|
C     |                 (ALGORITHM DESTROYS ELEMENTS)          |
C     |                                                        |
C     |         N     --DIMENSION OF MATRIX STORED IN A        |
C     |                                                        |
C     |    OUTPUT:                                             |
C     |                                                        |
C     |         E     --IMPROVED ESTIMATE FOR EIGENVALUE       |
C     |                                                        |
C     |         X     --ARRAY CONTAINING EIGENVECTOR           |
C     |                                                        |
C     |   BUILTIN FUNCTIONS: ABS                               |
C     |________________________________________________________|
C
      SUBROUTINE EVECT(E,X,A,N)
      REAL A(1),X(1),E,R,S,T,U,V
      INTEGER F,G,H,I,J,K,L,M,N,O,P,Q
      IF ( N .GT. 1 ) GOTO 10
      E = A(1)
      X(1) = 1.
      RETURN
10    G = (N*(N+1))/2
      F = G
      H = G + N - 2
      M = 1
      L = 2
      J = 0
      K = 1
C     ------------------------------------------------------
C     |*** MAKE ROOM FOR PIVOTS AND SUBTRACT EIGENVALUE ***|
C     ------------------------------------------------------
20    IF ( K .GT. N ) GOTO 40
      P = L - 1
      A(P) = A(P) - E
      IF ( K .EQ. N ) L = P
      DO 30 I = M,L
30         A(I-J) = A(I)
      J = K
      K = K + 1
      A(H+K) = A(L)
      M = L + 1
      L = M + K
      GOTO 20
40    I = G + 1
      J = H + 2
      K = H + N
50    IF ( I .GT. K ) GOTO 60
      A(I) = A(J)
      I = I + 2
      J = J + 1
      GOTO 50
60    K = N
      M = G
      L = G - N + 1
      G = H + N - 1
      R = ABS(A(M)) + ABS(A(G))
70    IF ( K .EQ. 1 ) GOTO 140
      J = K
      K = K - 1
      S = A(M)
      T = A(G)
      IF ( ABS(S) .GE. ABS(T) ) GOTO 100
C     --------------------------------------------
C     |*** PIVOT AND PERFORM ELIMINATION STEP ***|
C     --------------------------------------------
      A(M) = T
      IF ( R .LT. ABS(T) ) GOTO 80
      R = ABS(T)
      P = J
      O = M
80    T = S/T
      A(G) = T
      A(G+1) = 1.
      G = G - 2
      L = L - K
      M = M - J
      DO 90 I = L,M
           Q = I + K
           S = A(I)
           A(I) = A(Q) - T*S
90         A(Q) = S
      GOTO 70
C     --------------------------
C     |*** ELIMINATION STEP ***|
C     --------------------------
100   IF ( R .LT. ABS(T) ) GOTO 110
      R = ABS(T)
      P = J
      O = M
110   IF ( S .EQ. 0. ) GOTO 130
      T = T/S
      L = L - K
      M = M - J
      DO 120 I = L,M
120        A(I) = A(I) - T*A(I+K)
130   A(G) = T
      A(G+1) = 0.
      G = G - 2
      GOTO 70
140   T = ABS(A(1))
      IF ( T .GT. R ) GOTO 150
      R = T
      P = K
C     ------------------------------------------------------
C     |*** COMPUTE INITIAL APPROXIMATION TO EIGENVECTOR ***|
C     ------------------------------------------------------
150   DO 160 I = 1,N
160        X(I) = 0.
      T = 1.
      L = F
      J = O - P
      K = P
      GOTO 180
170   T = X(K)/A(J+K)
180   X(K) = T
      IF ( K .EQ. 1 ) GOTO 200
      K = K - 1
      DO 190 I = 1,K
190        X(I) = X(I) - T*A(I+J)
      J = J - K
      GOTO 170
200   IF ( K .EQ. N ) GOTO 210
      J = K
      K = K + 1
      L = L + 2
      X(K) = X(K) - A(L-1)*X(J)
      IF ( A(L) .EQ. 0. ) GOTO 200
      T = X(K)
      X(K) = X(J)
      X(J) = T
      GOTO 200
210   IF ( R .EQ. 0. ) GOTO 310
C     -----------------------------------------------------
C     |*** APPLY ONE ITERATION OF INVERSE POWER METHOD ***|
C     -----------------------------------------------------
      S = 0.
      DO 220 I = 1,N
           T = ABS(X(I))
220        IF ( T .GT. S ) S = T
      P = H + N
      R = 0.
      S = 1./S
      DO 230 I = 1,N
           T = S*X(I)
           R = R + T*T
           X(I) = T
230        A(I+P) = T
      V = 0.
      L = F
      J = L - N
      K = N
240   T = X(K)/A(J+K)
      X(K) = T
      IF ( K .EQ. 1 ) GOTO 260
      K = K - 1
      DO 250 I = 1,K
250        X(I) = X(I) - T*A(I+J)
      J = J - K
      GOTO 240
260   U = ABS(X(1))
270   IF ( K .EQ. N ) GOTO 290
      J = K
      K = K + 1
      L = L + 2
      T = X(J)
      S = X(K) - A(L-1)*T
      X(K) = S
      IF ( ABS(S) .GT. U ) U = ABS(S)
      IF ( A(L) .EQ. 1. ) GOTO 280
      V = V + A(P+J)*T
      GOTO 270
280   X(J) = S
      X(K) = T
      V = V + A(P+J)*X(J)
      GOTO 270
290   V = V + A(P+N)*X(N)
      IF ( V .NE. 0. ) V = R/V
      S = 0.
      T = 1./U
      DO 300 I = 1,N
           S = S + (V*X(I))**2
300        X(I) = T*X(I)
C     --------------------------------------------------------------
C     |*** USE RAYLEIGH QUOTIENT TO IMPROVE EIGENVALUE ESTIMATE ***|
C     --------------------------------------------------------------
      IF ( R+R .GE. S ) E = E + V
      RETURN
310   U = 0.
      DO 320 I = 1,N
           T = ABS(X(I))
320        IF ( T .GT. U ) U = T
      T = 1./U
      DO 330 I = 1,N
330        X(I) = T*X(I)
      RETURN
      END