C
C      ________________________________________________________
C     |                                                        |
C     | COMPUTE ALL EIGENVALUES OF A COMPLEX HESSENBERG MATRIX |
C     |                                                        |
C     |    INPUT:                                              |
C     |                                                        |
C     |         A     --COEFFICIENTS OF HESSENBERG MATRIX      |
C     |                 PACKED AT START OF COMPLEX ARRAY       |
C     |                                                        |
C     |         N     --DIMENSION OF MATRIX STORED IN A        |
C     |                                                        |
C     |         W     --COMPLEX WORK ARRAY WITH AT LEAST       |
C     |                 N ELEMENTS                             |
C     |                                                        |
C     |         D,P   --REAL WORK ARRAYS WITH AT LEAST N       |
C     |                 ELEMENTS IN EACH ARRAY                 |
C     |                                                        |
C     |    OUTPUT:                                             |
C     |                                                        |
C     |         E     --COMPLEX ARRAY OF EIGENVALUES           |
C     |                                                        |
C     |    BUILTIN FUNCTIONS: AMAX1,CONJG,SQRT                 |
C     |    PACKAGE FUNCTIONS: MAG,SQR                          |
C     |    PACKAGE SUBROUTINES: CEIG                           |
C     |________________________________________________________|
C
      SUBROUTINE VLS(E,A,N,D,P,W)
      COMPLEX A(1),E(1),W(1),X,Y,Z,Z1,Z2,Z3,Z4
      REAL D(1),P(1),B,Q,R,S,T,T0,T1,T2,U,MAG,SQR
      INTEGER F,F0,G,G0,H,I,J,K,KL,KM,KS,K0,L,LL,L0,L1,M,N,NS,O
      DO 10 I = 1,N
10         D(I) = 1.
      IF ( N .EQ. 1 ) GOTO 330
      B = 65536.**(-3)
C     ---------------------------------
C     |*** COMPUTE MACHINE EPSILON ***|
C     ---------------------------------
      T = 1.
20    T = .5*T
      S = 1. + T
      IF ( S .GT. 1. ) GOTO 20
      T0 = T + T
      T2 = T0*T0
      NS = 50*N
      LL = 0
      KL = (N*(N+3))/2 - 1
      M = N
      KM = KL - M
30    I = KM
      J = M
C     --------------------------------------------
C     |*** CHECK FOR ZERO SUBDIAGONAL ELEMENT ***|
C     --------------------------------------------
40    IF ( MAG(A(I)) .EQ. 0. ) GOTO 50
      I = I - J
      J = J - 1
      IF ( J .GT. 1 ) GOTO 40
50    K0 = J
      L = J + 1
      G = I + L
      F = G - 1
      F0 = F
      G0 = G
      L0 = L
      IF ( J .EQ. M ) GOTO 320
      DO 60 I = J,M
60         P(I) = SQRT(D(I))
      S = 0.
      K = J
C     --------------------------------
C     |*** COMPUTE NORM OF MATRIX ***|
C     --------------------------------
70    H = F - J
      T = P(K)
      DO 80 I = F,G
80         S = S + T*MAG(A(I))*P(I-H)
      K = L
      L = L + 1
      F = F + K
      G = G + L
      IF ( K .LT. M ) GOTO 70
      IF ( K .GT. M ) GOTO 90
      G = G - 1
      GOTO 70
90    IF ( S .EQ. 0. ) S = 1.
      T1 = 1./(S*T0)
      L1 = 0
      GOTO 300
100   LL = LL + 1
      IF ( LL .GT. NS ) GOTO 340
      T = 1.
      DO 110 I = K0,M
110        IF ( D(I) .LT. T ) T = D(I)
      IF ( T .GT. B ) GOTO 150
C     ----------------------------
C     |*** RESCALE THE MATRIX ***|
C     ----------------------------
      DO 120 I = K0,M
           P(I) = SQRT(D(I))
120        D(I) = 1.
      F = F0
      G = G0
      L = L0
      K = K0
130   H = F - K0
      T = P(K)
      DO 140 I = F,G
140        A(I) = A(I)*T*P(I-H)
      K = L
      L = L + 1
      F = F + K
      G = G + L
      IF ( K .LT. M ) GOTO 130
      IF ( K .GT. M ) GOTO 150
      G = G - 1
      GOTO 130
C     -----------------------
C     |*** COMPUTE SHIFT ***|
C     -----------------------
150   S = D(M-1)
      T = D(M)
      Z1 = A(KM-1)*S
      Z2 = A(KM)*S
      Z3 = A(KL-1)*T
      Z4 = A(KL)*T
      CALL CEIG(Z,Z,Z1,Z2,Z3,Z4)
160   K = K0
      L = L0
      F = F0
      G = G0
      KS = K
      Z1 = D(K)*A(F) - Z
      Z2 = D(K)*A(G)
C     ---------------------------------
C     |*** COMPUTE GIVENS ROTATION ***|
C     ---------------------------------
170   T = MAG(Z1) + MAG(Z2)
      IF ( T .EQ. 0. ) GOTO 190
      T = 1./T
      Q = SQR(Z1,T)
      R = SQR(Z2,T)
      Q = D(K)*Q
      R = D(L)*R
      IF ( Q .GT. R ) GOTO 180
      Z4 = Z1/Z2
      Z3 = (D(K)/D(L))*CONJG(Z4)
      E(K) = Z3
      W(K) = Z4
      P(K) = 1.
      S = R/(Q+R)
      R = D(K)
      D(K) = D(L)*S
      Q = D(L)
      D(L) = R*S
      GOTO 200
180   Z4 = Z2/Z1
      Z3 = (D(L)/D(K))*CONJG(Z4)
      E(K) = Z3
      W(K) = Z4
      P(K) = 0.
      S = Q/(Q+R)
      D(K) = D(K)*S
      Q = D(L)
      D(L) = D(L)*S
      GOTO 200
190   P(K) = 0.
      Z3 = (0.,0.)
      Z4 = (0.,0.)
      E(K) = Z3
      W(K) = Z4
      Q = D(L)
200   IF ( K .GT. KS ) GOTO 220
C     ------------------------------------------
C     |*** PROCESS INITIAL COLUMN FROM LEFT ***|
C     ------------------------------------------
      IF ( P(K) .EQ. 0. ) GOTO 210
      Y = A(G) + Z3*A(F)
      A(G) = Z4*A(G) - A(F)
      A(F) = Y
      GOTO 240
210   Y = A(F) + Z3*A(G)
      A(G) = A(G) - Z4*A(F)
      A(F) = Y
      GOTO 240
C     -----------------------------------
C     |*** PROCESS COLUMNS FROM LEFT ***|
C     -----------------------------------
220   I = G - 1
      IF ( P(K) .EQ. 0. ) GOTO 230
      X = A(H)*Z3 + X
      A(H) = X
      Y = A(G) + Z3*A(I)
      A(G) = Z4*A(G) - A(I)
      A(I) = Y
      IF ( MAG(X) .NE. 0. ) GOTO 240
      K0 = K
      L0 = L
      F0 = I
      G0 = G
      GOTO 240
230   X = X*Z3 + A(H)
      A(H) = X
      Y = A(I) + Z3*A(G)
      A(G) = A(G) - Z4*A(I)
      A(I) = Y
      IF ( MAG(X) .NE. 0. ) GOTO 240
      K0 = K
      L0 = L
      F0 = I
      G0 = G
240   J = K
      K = L
      L = L + 1
      DO 260 I = KS,J
           O = G + I
           H = O + 1
           IF ( P(I) .EQ. 0. ) GOTO 250
           Y = A(H) + E(I)*A(O)
           A(H) = W(I)*A(H) - A(O)
           A(O) = Y
           GOTO 260
250        Y = A(O) + E(I)*A(H)
           A(H) = A(H) - W(I)*A(O)
           A(O) = Y
260   CONTINUE
C     ------------------------------------
C     |*** PROCESS COLUMNS FROM RIGHT ***|
C     ------------------------------------
      Z3 = CONJG(Z3)
      Z4 = CONJG(Z4)
      IF ( P(J) .EQ. 0. ) GOTO 280
      Z1 = Q*A(G+K) - Z*W(J)
      DO 270 I = F,G
           H = I + K
           Y = A(H) + Z3*A(I)
           A(H) = Z4*A(H) - A(I)
270        A(I) = Y
      IF ( L .GT. M ) GOTO 300
      F = F + K
      H = G
      G = G + L
      Z2 = Q*A(G)
      X = A(G)
      A(G) = Z4*A(G)
      GOTO 170
280   Z1 = Q*A(G+K) - Z
      DO 290 I = F,G
           H = I + K
           Y = A(I) + Z3*A(H)
           A(H) = A(H) - Z4*A(I)
290        A(I) = Y
      IF ( L .GT. M ) GOTO 300
      F = F + K
      H = G
      G = G + L
      Z2 = Q*A(G)
      X = Z3*A(G)
      GOTO 170
300   T = D(M)
      S = D(M-1)
      Q = MAG(A(KM))
C     ------------------------------
C     |*** TEST FOR CONVERGENCE ***|
C     ------------------------------
      IF ( ((T*Q)*T1)*((S*Q)*T1) .GT. 1. ) GOTO 100
      L1 = L1 + 1
      IF ( L1 .GT. 30 ) GOTO 310
      R = MAG(A(KL))
      U = AMAX1(Q,R)
      IF ( U .EQ. 0. ) GOTO 310
      IF ( (S*Q)*(Q/U) .GT. T2*(T*R)*(R/U) ) GOTO 100
310   L1 = 0
      E(M) = A(KL)*D(M)
      KL = KM - 1
      KM = KM - M
      M = M - 1
      IF ( M .GT. K0 ) GOTO 300
320   E(M) = A(KL)*D(M)
      KL = KM - 1
      KM = KM - M
      M = M - 1
      IF ( M .GT. 1 ) GOTO 30
330   E(1) = A(1)*D(1)
      P(1) = N
      RETURN
340   M = N - M + 1
      WRITE(6,*) 'SINCE THE STOPPING CRITERION NOT SATISFIED'
      WRITE(6,*) 'AFTER',NS,'ITERATIONS, WE STOP WHILE COMPUTING'
      WRITE(6,*) 'EIGENVALUE NUMBER',M
      P(1) = M
      RETURN
      END