' ******************************************************************
' Eigenvalues and eigenvectors of a real upper Hessenberg matrix
' ******************************************************************

' ------------------------------------------------------------------
' Type definition
' ------------------------------------------------------------------

TYPE Complex
  X AS DOUBLE
  Y AS DOUBLE
END TYPE

' ------------------------------------------------------------------
' Global variable
' ------------------------------------------------------------------

COMMON SHARED ErrCode AS INTEGER
' Error code from the last function evaluation

' ------------------------------------------------------------------
' External functions
' ------------------------------------------------------------------

DECLARE FUNCTION Max (X AS DOUBLE, Y AS DOUBLE) AS DOUBLE
' Maximum of 2 reals

DECLARE FUNCTION IMin (X AS INTEGER, Y AS INTEGER) AS INTEGER
' Minimum of 2 integers

' ******************************************************************

SUB Cdiv(      Ar AS DOUBLE,       Ai AS DOUBLE, _ 
               Br AS DOUBLE,       Bi AS DOUBLE, _
         BYREF Cr AS DOUBLE, BYREF Ci AS DOUBLE)
' Complex division, (Cr,Ci) = (Ar,Ai)/(Br,Bi)

  DIM AS DOUBLE S, Ars, Ais, Brs, Bis
  
  S = ABS(Br) + ABS(Bi)
  Ars = Ar / S
  Ais = Ai / S
  Brs = Br / S
  Bis = Bi / S
  S = Brs * Brs + Bis * Bis
  Cr = (Ars * Brs + Ais * Bis) / S
  Ci = (Ais * Brs - Ars * Bis) / S
END SUB

FUNCTION DSgn(A AS DOUBLE, B AS DOUBLE) AS DOUBLE
  IF B < 0 THEN DSgn = -ABS(A) ELSE DSgn = ABS(A)
END FUNCTION

SUB Hqr2(H()      AS DOUBLE,_ 
         Lo       AS INTEGER,_ 
         Hi       AS INTEGER,_ 
         Lambda() AS Complex,_ 
         Z()      AS DOUBLE)
' ------------------------------------------------------------------
' This function is a translation of the EISPACK subroutine hqr2
'
' This procedure finds the eigenvalues and eigenvectors of a real
' upper Hessenberg matrix by the QR method.
'
' On input:
'
'   H contains the upper Hessenberg matrix.
'
'   Lo and Hi are integers determined by the balancing subroutine
'   Balance. If Balance has not been used, set Lo=Lbound(H),
'   Hi=Ubound(H)
'
'   Z contains the transformation matrix produced by Eltran after
'   the reduction by Elmhes, or by Ortran after the reduction by
'   Orthes, if performed. If the eigenvectors of the Hessenberg
'   matrix are desired, Z must contain the identity matrix.
'
' On output:
'
'   H has been destroyed.
'
'   Lambda contains the eigenvalues. They are unordered except that
'   complex conjugate pairs of values appear consecutively with the
'   eigenvalue having the positive imaginary part first.
'
'   Z contains the real and imaginary parts of the eigenvectors. If
'   the i-th eigenvalue is real, the i-th column of Z contains its
'   eigenvector. If the i-th eigenvalue is complex with positive
'   imaginary part, the i-th and (i+1)-th columns of Z contain the
'   real and imaginary parts of its eigenvector. The eigenvectors
'   are unnormalized. If an error exit is made, none of the
'   eigenvectors has been found.
'
'   The global variable ErrCode contains an error code:
'      zero  for normal return,
'      -j    if the limit of 30*N iterations is exhausted
'            while the j-th eigenvalue is being sought
'            (N being the size of the matrix). The eigenvalues
'            should be correct for indices j+1,...,Ubound(H).
' ------------------------------------------------------------------
' Note: This is a crude translation. Many of the original goto's 
' have been kept !
' ------------------------------------------------------------------

  DIM AS INTEGER Lb, Ub, I, J, K, L, M, N, En, Na
  DIM AS INTEGER Itn, Its, Mp2, Enm2, NotLas
  DIM AS DOUBLE  P, Q, R, S, T, W, X, Y, Zz 
  DIM AS DOUBLE  Ra, Sa, Vi, Vr, Norm, Tst1, Tst2

  Lb = LBOUND(H, 1)
  Ub = UBOUND(H, 1)

  ' Store roots isolated by Balance and compute matrix norm
  K = Lb
  Norm = 0
  FOR I = Lb TO Ub
    FOR J = K TO Ub
      Norm = Norm + ABS(H(I, J))
    NEXT J
    K = I
    IF I < Lo OR I > Hi THEN
      Lambda(I).X = H(I, I)
      Lambda(I).Y = 0       
    END IF
  NEXT I

  N = Ub - Lb + 1
  Itn = 30 * N
  En = Hi
  T = 0

L60: ' Search for next eigenvalues

  IF En < Lo THEN GOTO L340
  
  Its = 0
  Na = En - 1
  Enm2 = Na - 1

L70: ' Look for single small sub-diagonal element

  FOR L = En TO Lo STEP -1
    IF L = Lo THEN GOTO L100
    S = ABS(H(L - 1, L - 1)) + ABS(H(L, L))
    IF S = 0 THEN S = Norm
    Tst1 = S
    Tst2 = Tst1 + ABS(H(L, L - 1))
    IF Tst2 = Tst1 THEN GOTO L100
  NEXT L

L100: ' Form shift

  X = H(En, En)
  IF L = En THEN GOTO L270
  Y = H(Na, Na)
  W = H(En, Na) * H(Na, En)
  IF L = Na THEN GOTO L280

  ' Set error -- all eigenvalues have not
  ' converged after 30*N iterations
  IF Itn = 0 THEN ErrCode = -En: EXIT SUB

  IF Its <> 10 AND Its <> 20 THEN GOTO L130

  ' Form exceptional shift
  T = T + X

  FOR I = Lo TO En
    H(I, I) = H(I, I) - X
  NEXT I

  S = ABS(H(En, Na)) + ABS(H(Na, Enm2))
  X = 0.75 * S
  Y = X
  W = -0.4375 * S * S

L130:

  Its = Its + 1
  Itn = Itn - 1

  ' Look for two consecutive small sub-diagonal elements
  FOR M = Enm2 TO L STEP -1
    Zz = H(M, M)
    R = X - Zz
    S = Y - Zz
    P = (R * S - W) / H(M + 1, M) + H(M, M + 1)
    Q = H(M + 1, M + 1) - Zz - R - S
    R = H(M + 2, M + 1)
    S = ABS(P) + ABS(Q) + ABS(R)
    P = P / S
    Q = Q / S
    R = R / S
    IF M = L THEN GOTO L150
    Tst1 = ABS(P) * (ABS(H(M - 1, M - 1)) + ABS(Zz) + _
           ABS(H(M + 1, M + 1)))
    Tst2 = Tst1 + ABS(H(M, M - 1)) * (ABS(Q) + ABS(R))
    IF Tst2 = Tst1 THEN GOTO L150
  NEXT M

L150:

  Mp2 = M + 2

  FOR I = Mp2 TO En
    H(I, I - 2) = 0
    IF I <> Mp2 THEN H(I, I - 3) = 0
  NEXT I

  ' Double QR step involving rows L to En and columns M to En
  FOR K = M TO Na
    NotLas = (K <> Na)
    IF K = M THEN GOTO L170
    
    P = H(K, K - 1)
    Q = H(K + 1, K - 1)
    R = 0
    
    IF NotLas THEN R = H(K + 2, K - 1)
    
    X = ABS(P) + ABS(Q) + ABS(R)
    IF X = 0 THEN GOTO L260
    
    P = P / X
    Q = Q / X
    R = R / X

L170:

    S = DSgn(SQR(P * P + Q * Q + R * R), P)
    
    IF K <> M THEN
      H(K, K - 1) = -S * X
    ELSEIF L <> M THEN
      H(K, K - 1) = -H(K, K - 1)
    END IF
    
    P = P + S
    X = P / S
    Y = Q / S
    Zz = R / S
    Q = Q / P
    R = R / P
    
    IF NotLas THEN GOTO L225

    ' Row modification
    FOR J = K TO Ub
      P = H(K, J) + Q * H(K + 1, J)
      H(K, J) = H(K, J) - P * X
      H(K + 1, J) = H(K + 1, J) - P * Y
    NEXT J

    J = IMin(En, K + 3)

    ' Column modification
    FOR I = Lb TO J
      P = X * H(I, K) + Y * H(I, K + 1)
      H(I, K) = H(I, K) - P
      H(I, K + 1) = H(I, K + 1) - P * Q
    NEXT I

    ' Accumulate transformations
    FOR I = Lo TO Hi
      P = X * Z(I, K) + Y * Z(I, K + 1)
      Z(I, K) = Z(I, K) - P
      Z(I, K + 1) = Z(I, K + 1) - P * Q
    NEXT I
    
    GOTO L260

L225:

    ' Row modification
    FOR J = K TO Ub
      P = H(K, J) + Q * H(K + 1, J) + R * H(K + 2, J)
      H(K, J) = H(K, J) - P * X
      H(K + 1, J) = H(K + 1, J) - P * Y
      H(K + 2, J) = H(K + 2, J) - P * Zz
    NEXT J

    J = IMin(En, K + 3)

    ' Column modification
    FOR I = Lb TO J
      P = X * H(I, K) + Y * H(I, K + 1) + Zz * H(I, K + 2)
      H(I, K) = H(I, K) - P
      H(I, K + 1) = H(I, K + 1) - P * Q
      H(I, K + 2) = H(I, K + 2) - P * R
    NEXT I

    ' Accumulate transformations
    FOR I = Lo TO Hi
      P = X * Z(I, K) + Y * Z(I, K + 1) + Zz * Z(I, K + 2)
      Z(I, K) = Z(I, K) - P
      Z(I, K + 1) = Z(I, K + 1) - P * Q
      Z(I, K + 2) = Z(I, K + 2) - P * R
    NEXT I

L260:

  NEXT K

  GOTO L70

L270:

  ' One root found
  H(En, En) = X + T
  Lambda(En).X = H(En, En)
  Lambda(En).Y = 0
  En = Na
  GOTO L60

L280:

  ' Two roots found
  P = 0.5 * (Y - X)
  Q = P * P + W
  Zz = SQR(ABS(Q))
  H(En, En) = X + T
  X = H(En, En)
  H(Na, Na) = Y + T
  IF Q < 0 THEN GOTO L320

  ' Real pair
  Zz = P + DSgn(Zz, P)
  Lambda(Na).X = X + Zz
  Lambda(En).X = Lambda(Na).X
  IF Zz <> 0 THEN Lambda(En).X = X - W / Zz
  Lambda(Na).Y = 0
  Lambda(En).Y = 0
  X = H(En, Na)
  S = ABS(X) + ABS(Zz)
  P = X / S
  Q = Zz / S
  R = SQR(P * P + Q * Q)
  P = P / R
  Q = Q / R

  ' Row modification
  FOR J = Na TO Ub
    Zz = H(Na, J)
    H(Na, J) = Q * Zz + P * H(En, J)
    H(En, J) = Q * H(En, J) - P * Zz
  NEXT J

  ' Column modification
  FOR I = Lb TO En
    Zz = H(I, Na)
    H(I, Na) = Q * Zz + P * H(I, En)
    H(I, En) = Q * H(I, En) - P * Zz
  NEXT I

  ' Accumulate transformations
  FOR I = Lo TO Hi
    Zz = Z(I, Na)
    Z(I, Na) = Q * Zz + P * Z(I, En)
    Z(I, En) = Q * Z(I, En) - P * Zz
  NEXT I

  GOTO L330

L320:

  ' Complex pair
  Lambda(Na).X = X + P
  Lambda(En).X = X + P
  Lambda(Na).Y = Zz
  Lambda(En).Y = -Zz

L330:

  En = Enm2
  GOTO L60

L340:

  IF Norm = 0 THEN ErrCode = 0 : EXIT SUB

  ' All roots found. Backsubstitute to find
  ' vectors of upper triangular form
  FOR En = Ub TO Lb STEP -1
    P = Lambda(En).X
    Q = Lambda(En).Y
    Na = En - 1
    IF Q < 0 THEN
      GOTO L710
    ELSEIF Q = 0 THEN
      GOTO L600
    ELSE
      GOTO L800
    END IF

L600:

    ' Real vector
    M = En
    H(En, En) = 1
    IF Na < Lb THEN GOTO L800

    FOR I = Na TO Lb STEP -1
      W = H(I, I) - P
      R = 0

      FOR J = M TO En
        R = R + H(I, J) * H(J, En)
      NEXT J

      IF Lambda(I).Y >= 0 THEN GOTO L630
      Zz = W
      S = R
      GOTO L700
L630:
      M = I
      IF Lambda(I).Y <> 0 THEN GOTO L640
      T = W
      IF T <> 0 THEN GOTO L635
      Tst1 = Norm
      T = Tst1

      DO
        T = 0.01 * T
        Tst2 = Norm + T
      LOOP UNTIL Tst2 <= Tst1
L635:
      H(I, En) = -R / T
      GOTO L680
L640:
      ' Solve real equations
      X = H(I, I + 1)
      Y = H(I + 1, I)
      Q = (Lambda(I).X - P) * (Lambda(I).X - P) + _
           Lambda(I).Y * Lambda(I).Y
      T = (X * S - Zz * R) / Q
      H(I, En) = T
      IF ABS(X) > ABS(Zz) THEN
        H(I + 1, En) = (-R - W * T) / X
      ELSE
        H(I + 1, En) = (-S - Y * T) / Zz
      END IF
L680:
      ' Overflow control
      T = ABS(H(I, En))
      IF T = 0 THEN GOTO L700
      Tst1 = T
      Tst2 = Tst1 + 1 / Tst1
      IF Tst2 > Tst1 THEN GOTO L700
      FOR J = I TO En
        H(J, En) = H(J, En) / T
      NEXT J
L700:
    NEXT I

    ' End real vector
    GOTO L800
    
L710: ' Complex vector
    M = Na

    ' Last vector component chosen imaginary so that
    ' eigenvector matrix is triangular
    IF ABS(H(En, Na)) > ABS(H(Na, En)) THEN
      H(Na, Na) = Q / H(En, Na)
      H(Na, En) = -(H(En, En) - P) / H(En, Na)
    ELSE
      Cdiv 0, -H(Na, En), H(Na, Na) - P, Q, H(Na, Na), H(Na, En)
    END IF

    H(En, Na) = 0
    H(En, En) = 1
    Enm2 = Na - 1
    IF Enm2 < Lb THEN GOTO L800

    FOR I = Enm2 TO Lb STEP -1
      W = H(I, I) - P
      Ra = 0
      Sa = 0

      FOR J = M TO En
        Ra = Ra + H(I, J) * H(J, Na)
        Sa = Sa + H(I, J) * H(J, En)
      NEXT J

      IF Lambda(I).Y >= 0 THEN GOTO L770
      Zz = W
      R = Ra
      S = Sa
      GOTO L795
L770:
      M = I
      IF Lambda(I).Y <> 0 THEN GOTO L780
      Cdiv -Ra, -Sa, W, Q, H(I, Na), H(I, En)
      GOTO L790
L780:
      ' Solve complex equations
      X = H(I, I + 1)
      Y = H(I + 1, I)
      Vr = (Lambda(I).X - P) * (Lambda(I).X - P) + _
            Lambda(I).Y * Lambda(I).Y - Q * Q 
      Vi = (Lambda(I).X - P) * 2 * Q
      IF Vr = 0 AND Vi = 0 THEN
        Tst1 = Norm * (ABS(W) + ABS(Q) + ABS(X) + ABS(Y) + ABS(Zz))
        Vr = Tst1
        DO
          Vr = 0.01 * Vr
          Tst2 = Tst1 + Vr
        LOOP UNTIL Tst2 <= Tst1
      END IF

      Cdiv X * R - Zz * Ra + Q * Sa, X * S - Zz * Sa - Q * Ra, _
           Vr, Vi, H(I, Na), H(I, En)

      IF ABS(X) > ABS(Zz) + ABS(Q) THEN
        H(I + 1, Na) = (-Ra - W * H(I, Na) + Q * H(I, En)) / X
        H(I + 1, En) = (-Sa - W * H(I, En) - Q * H(I, Na)) / X
      ELSE
        Cdiv -R - Y * H(I, Na), -S - Y * H(I, En), Zz, Q, _
             H(I + 1, Na), H(I + 1, En)
      END IF
L790:
      ' Overflow control
      T = Max(ABS(H(I, Na)), ABS(H(I, En)))
      IF T = 0 THEN GOTO L795

      Tst1 = T
      Tst2 = Tst1 + 1 / Tst1
      IF Tst2 > Tst1 THEN GOTO L795

      FOR J = I TO En
        H(J, Na) = H(J, Na) / T
        H(J, En) = H(J, En) / T
      NEXT J
L795:
    NEXT I

    ' End complex vector
L800:
  NEXT En

  ' End back substitution.
  ' Vectors of isolated roots
  FOR I = Lb TO Ub
    IF I < Lo OR I > Hi THEN
      FOR J = I TO Ub
        Z(I, J) = H(I, J)
      NEXT J
    END IF
  NEXT I

  ' Multiply by transformation matrix to give
  ' vectors of original full matrix.
  FOR J = Ub TO Lo STEP -1
    M = IMin(J, Hi)
    FOR I = Lo TO Hi
      Zz = 0
      FOR K = Lo TO M
        Zz = Zz + Z(I, K) * H(K, J)
      NEXT K
      Z(I, J) = Zz
    NEXT I
  NEXT J

  ErrCode = 0
END SUB
