C
C      ________________________________________________________
C     |                                                        |
C     |    COMPUTE EIGENVECTOR CORRESPONDING TO GIVEN REAL     |
C     |      EIGENVALUE FOR A REAL SYMMETRIC BAND MATRIX       |
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         E     --EIGENVALUE                             |
C     |                                                        |
C     |         A     --ARRAY CONTAINING BANDS ON DIAGONAL AND |
C     |                 BELOW DIAGONAL                         |
C     |                                                        |
C     |         LA    --LEADING (ROW) DIMENSION OF ARRAY A     |
C     |                                                        |
C     |         N     --MATRIX DIMENSION                       |
C     |                                                        |
C     |         H     --HALF BANDWIDTH                         |
C     |                                                        |
C     |         W     --WORK ARRAY (LENGTH AT LEAST (3H+2)N)   |
C     |                                                        |
C     |    OUTPUT:                                             |
C     |                                                        |
C     |         E     --IMPROVED ESTIMATE FOR EIGENVALUE       |
C     |                                                        |
C     |         X     --EIGENVECTOR                            |
C     |                                                        |
C     |    BUILTIN FUNCTIONS: ABS,MAX0,MIN0,SQRT               |
C     |________________________________________________________|
C
      SUBROUTINE HVECT(E,X,A,LA,N,H,W)
      INTEGER B,C,D,F,G,H,I,J,K,L,LA,M,N,O,P,Q,Y,Z
      REAL A(LA,1),W(1),X(1),E,R,S,T,U,V
      IF ( H .GT. 0 ) GOTO 20
      S = ABS(A(1,1)-E)
      DO 10 J = 1,N
           X(J) = 0.
           T = ABS(A(1,J)-E)
           IF ( T .GT. S ) GOTO 10
           S = T
           K = J
10    CONTINUE
      X(K) = 1.
      E = A(1,K)
      RETURN
C     ---------------------------------------------------------
C     |*** STORE MATRIX IN W ARRAY AND SUBTRACT EIGENVALUE ***|
C     ---------------------------------------------------------
20    C = 1 + 3*H
      D = H + H
      P = H + 1
      M = 1 - H
      DO 70 J = 1,N
           X(J) = 0.
           L = M + D
           M = M + P
           DO 30 I = M,L
30              W(I) = 0.
           M = L + P
           W(M) = A(1,J) - E
           F = M
           L = N - J
           IF ( L .GE. H ) GOTO 40
           IF ( L .EQ. 0 ) GOTO 70
           L = L + M
           GOTO 50
40         L = M + H
50         M = M + 1
           K = 2 - M
           DO 60 I = M,L
                T = A(K+I,J)
                W(I) = T
                F = F + C
60              W(F) = T
70    CONTINUE
      R = 0.
      Q = P + P
      L = C + 1
      DO 80 I = Q,L
           T = ABS(W(I))
           IF ( T .GT. R ) R = T
80    CONTINUE
      I = -H
      Z = P + H
C     ---------------------------
C     |*** FACTOR THE MATRIX ***|
C     ---------------------------
      K = 0
90    K = K + 1
      I = I + L
      IF ( K .EQ. N ) GOTO 160
      M = I + 1
      Q = I
      O = MIN0(H,N-K)
      P = I + O
C     ---------------------------------------
C     |*** FIND PIVOT AND START ROW SWAP ***|
C     ---------------------------------------
      S = ABS(W(I))
      DO 100 J = M,P
           T = ABS(W(J))
           IF ( S .GT. T ) GOTO 100
           Q = J
           S = T
100   CONTINUE
      IF ( R .LT. S ) GOTO 110
      R = S
      Y = K
      IF ( R .EQ. 0 ) GOTO 170
110   J = I - Z
      B = Q - I
      W(J) = K + B
      T = W(Q)
      W(Q) = W(I)
      W(I) = T
C     -----------------------------
C     |*** COMPUTE MULTIPLIERS ***|
C     -----------------------------
      DO 120 J = M,P
120        W(J) = W(J)/T
      F = I + C*MIN0(D,N-K)
      G = C - O
130   M = P + G
      P = M + B
      T = W(P)
      W(P) = W(M)
      W(M) = T
      P = M + O
      IF ( T .EQ. 0. ) GOTO 150
      Q = I - M
      M = M + 1
C     ------------------------------
C     |*** ELIMINATE BY COLUMNS ***|
C     ------------------------------
      DO 140 J = M,P
140        W(J) = W(J) - T*W(J+Q)
150   IF ( P .LT. F ) GOTO 130
      GOTO 90
160   IF ( ABS(W(I)) .GE. R ) GOTO 170
      R = W(I)
      Y = N
C     ---------------------------------------------------
C     |*** COMPUTE INITIAL EIGENVECTOR APPROXIMATION ***|
C     ---------------------------------------------------
170   K = Y
      J = C*K - H
      T = 1.
      GOTO 190
180   T = X(K)/W(J+K)
190   X(K) = T
      IF ( K .EQ. 1 ) GOTO 210
      Q = MAX0(1,K-D)
      K = K - 1
      DO 200 I = Q,K
200        X(I) = X(I) - T*W(I+J)
      J = J - C
      GOTO 180
210   IF ( R .EQ. 0. ) GOTO 360
      S = 0.
      DO 220 I = 1,N
           T = ABS(X(I))
220        IF ( T .GT. S ) S = T
      R = 0.
      S = 1./S
      DO 230 I = 1,N
           T = S*X(I)
           R = R + T*T
230        X(I) = T
      J = -H
      Q = 1
      P = H + 1
      DO 240 I = 1,P
240        W(I+1) = X(I)
C     -----------------------------------------------------
C     |*** APPLY ONE ITERATION OF INVERSE POWER METHOD ***|
C     -----------------------------------------------------
250   J = J + C
      K = Q
      IF ( K .EQ. N ) GOTO 270
      Q = K + 1
      I = W(J+K-Z)
      T = X(I)
      X(I) = X(K)
      X(K) = T
      O = MIN0(K+H,N)
      DO 260 I = Q,O
260        X(I) = X(I) - T*W(I+J)
      IF ( O .EQ. N ) GOTO 250
      W(Q+J) = X(O+1)
      GOTO 250
C     ---------------------------
C     |*** BACK SUBSTITUTION ***|
C     ---------------------------
270   U = 0.
280   T = X(K)/W(J+K)
      IF ( ABS(T) .GT. U ) U = ABS(T)
      X(K) = T
      IF ( K .EQ. 1 ) GOTO 300
      Q = MAX0(1,K-D)
      K = K - 1
      DO 290 I = Q,K
290        X(I) = X(I) - T*W(I+J)
      J = J - C
      GOTO 280
300   V = 0.
      DO 310 I = 1,P
310        V = V + X(I)*W(I+1)
      IF ( P .EQ. N ) GOTO 330
      J = Z + 1
      K = P + 1
      DO 320 I = K,N
           V = V + X(I)*W(J)
320        J = J + L
330   IF ( V .NE. 0. ) V = R/V
      S = 0.
      T = 1./U
      U = 0.
      DO 340 I = 1,N
           S = S + (X(I)*V)**2
340        U = U + (T*X(I))**2
      T = T/SQRT(U)
      DO 350 I = 1,N
350        X(I) = T*X(I)
C     --------------------------------------------------------------
C     |*** USE RAYLEIGH QUOTIENT TO IMPROVE EIGENVALUE ESTIMATE ***|
C     --------------------------------------------------------------
      IF ( R+R .GE. S ) E = E + V
      RETURN
360   U = 0.
      DO 370 I = 1,N
           T = ABS(X(I))
370        IF ( T .GT. U ) U = T
      T = 1./U
      U = 0.
      DO 380 I = 1,N
380        U = U + (T*X(I))**2
      T = T/SQRT(U)
      DO 390 I = 1,N
390        X(I) = T*X(I)
      RETURN
      END