' ******************************************************************
' Eigenvalues of a real upper Hessenberg matrix by the QR method
' ******************************************************************

' ------------------------------------------------------------------
' 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 function
' ------------------------------------------------------------------

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

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

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

SUB Hqr(H()      AS DOUBLE,_ 
        Lo       AS INTEGER,_ 
        Hi       AS INTEGER,_ 
        Lambda() AS Complex)
' ------------------------------------------------------------------
' This subroutine is a translation of the EISPACK subroutine hqr.
'
' This subroutine finds the eigenvalues 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)
'
' 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.
'
'   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, Z, 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 ErrCode = 0: EXIT SUB

  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
    Z = H(M, M)
    R = X - Z
    S = Y - Z
    P = (R * S - W) / H(M + 1, M) + H(M, M + 1)
    Q = H(M + 1, M + 1) - Z - 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(Z) + _
           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
    Z = R / S
    Q = Q / P
    R = R / P

    IF NotLas THEN GOTO L225

    ' Row modification
    FOR J = K TO En
       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 = L 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

    GOTO L260

L225:

    ' Row modification
    FOR J = K TO En
      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 * Z
    NEXT J

    J = IMin(En, K + 3)

    ' Column modification
    FOR I = L TO J
      P = X * H(I, K) + Y * H(I, K + 1) + Z * 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

L260:

  NEXT K

  GOTO L70

L270:

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

L280:

  ' Two roots found
  P = .5 * (Y - X)
  Q = P * P + W
  Z = SQR(ABS(Q))
  X = X + T
  IF Q < 0 THEN GOTO L320

  ' Real pair
  Z = P + DSgn(Z, P)
  Lambda(Na).X = X + Z
  Lambda(En).X = Lambda(Na).X
  IF Z <> 0 THEN Lambda(En).X = X - W / Z
  Lambda(Na).Y = 0
  Lambda(En).Y = 0
  GOTO L330

L320:

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

L330:

  En = Enm2
  GOTO L60

END SUB
