' ******************************************************************
' Reduction of a square matrix to upper Hessenberg form
' ******************************************************************

SUB ElmHes(A()     AS DOUBLE,  _ 
           Lo      AS INTEGER, _ 
           Hi      AS INTEGER, _ 
           Inter() AS INTEGER)
' ------------------------------------------------------------------
' This subroutine is a translation of the EISPACK subroutine Elmhes
'
' Given a real general matrix, this subroutine reduces a submatrix
' situated in rows and columns Lo through Hi to upper Hessenberg
' form by stabilized elementary similarity transformations.
'
' On input:
'
'   A() contains the input matrix.
'
'   Lo and Hi are integers determined by the balancing procedure
'   Balance. If Balance has not been used, set Lo=Lbound(A),
'   Hi=Ubound(A).
'
' On output:
'
'   A() contains the Hessenberg matrix. The multipliers which were
'   used in the reduction are stored in the remaining triangle under
'   the Hessenberg matrix.
'
'   Inter() contains information on the rows and columns interchanged
'   in the reduction. Only elements Lo through Hi are used.
' ------------------------------------------------------------------

  DIM AS INTEGER Lb, Ub, I, J, K, L, M, Mm1
  DIM AS DOUBLE  X, Y

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

  L = Hi - 1
  K = Lo + 1
  IF L < K THEN EXIT SUB

  FOR M = K TO L
    Mm1 = M - 1
    X = 0
    I = M

    FOR J = M TO Hi
      IF ABS(A(J, Mm1)) > ABS(X) THEN
        X = A(J, Mm1)
        I = J
      END IF
    NEXT J

    Inter(M) = I

    ' Interchange rows and columns of A
    IF I <> M THEN
      FOR J = Mm1 TO Ub
        SWAP A(I, J), A(M, J)
      NEXT J
      FOR J = Lb TO Hi
        SWAP A(J, I), A(J, M)
      NEXT J
    END IF

    IF X <> 0 THEN
      FOR I = M + 1 TO Hi
        Y = A(I, Mm1)
        IF Y <> 0 THEN
          Y = Y / X
          A(I, Mm1) = Y
          FOR J = M TO Ub
            A(I, J) = A(I, J) - Y * A(M, J)
          NEXT J
          FOR J = Lb TO Hi
            A(J, M) = A(J, M) + Y * A(J, I)
          NEXT J
        END IF
      NEXT I
    END IF
  NEXT M
END SUB
