' *****************************************************************
' Balances a matrix and tries to isolate eigenvalues
' *****************************************************************

#DEFINE True -1
#DEFINE False 0

DIM SHARED Lb AS INTEGER, Ub AS INTEGER

SUB Exchange (A()     AS DOUBLE,_
              Lo      AS INTEGER,_
              Hi      AS INTEGER,_
              J       AS INTEGER,_
              M       AS INTEGER,_
              Scale() AS DOUBLE)
' Row and column exchange

  DIM I AS INTEGER

  Scale(M) = J : IF J = M THEN EXIT SUB

  FOR I = Lb TO Hi 
    SWAP A(I, J), A(I, M)
  NEXT I

  FOR I = Lo TO Ub 
    SWAP A(J, I), A(M, I)
  NEXT I
END SUB

SUB Balance(      A()     AS DOUBLE,  _ 
            BYREF Lo      AS INTEGER, _ 
            BYREF Hi      AS INTEGER, _ 
                  Scale() AS DOUBLE)
' ------------------------------------------------------------------
' This procedure is a translation of the EISPACK procedure Balanc.
'
' This procedure balances a real matrix and isolates eigenvalues
' whenever possible.
'
' On input:
'
'   A contains the input matrix to be balanced.
'
' On output:
'
'   A contains the balanced matrix.
'
'   Lo and Hi are two integers such that A(i,j)
'   is equal to zero if
'     (1) i is greater than j and
'     (2) j=Lbound(A),...,Lo-1 or i=Hi+1,...,Ubound(A)
'
'   Scale contains information determining the permutations
'   and scaling factors used.
'
'   Suppose that the principal submatrix in rows Lo through Hi
'   has been balanced, that P(j) denotes the index interchanged
'   with j during the permutation step, and that the elements
'   of the diagonal matrix used are denoted by D(i,j).  then
'       Scale(j) = P(j),    for j = Lbound(A),...,Lo-1
'                = D(j,j),      j = Lo,...,Hi
'                = P(j)         j = Hi+1,...,Ubound(A).
'   the order in which the interchanges are made is
'   Ubound(A) to Hi+1, then Lbound(A) to Lo-1.
'
'   Note that Lbound(A) is returned for Hi if Hi is < Lbound(A)
'   formally
' ------------------------------------------------------------------

  CONST RADIX  = 2              ' Base for floating point numbers
  CONST RADIX2 = RADIX * RADIX
  
  DIM AS INTEGER I, J, M, Flag , Found, Conv
  DIM AS DOUBLE  C, F, G, R, S
  
  Lb = LBOUND(A, 1)
  Ub = UBOUND(A, 1)

  Lo = Lb
  Hi = Ub

  ' Search for rows isolating an eigenvalue and push them down
  DO
    J = Hi
    DO
      I = Lb
      DO
        Flag = (I <> J) AND (A(J, I) <> 0)
        I = I + 1
      LOOP UNTIL Flag OR (I > Hi)
      Found = NOT Flag
      IF Found THEN
        M = Hi
        Exchange A(), Lo, Hi, J, M, Scale()
        Hi = Hi - 1
      END IF
      J = J - 1
    LOOP UNTIL Found OR (J < Lb)
  LOOP UNTIL (NOT Found) OR (Hi < Lb)

  IF Hi < Lb THEN Hi = Lb
  IF Hi = Lb THEN EXIT SUB

  ' Search for columns isolating an eigenvalue and push them left
  DO
    J = Lo
    DO
      I = Lo
      DO
        Flag = (I <> J) AND (A(I, J) <> 0)
        I = I + 1
      LOOP UNTIL Flag OR (I > Hi)
      Found = NOT Flag
      IF Found THEN
        M = Lo
        Exchange A(), Lo, Hi, J, M, Scale()
        Lo = Lo + 1
      END IF
      J = J + 1
    LOOP UNTIL Found OR (J > Hi)
  LOOP UNTIL (NOT Found)

  ' Now balance the submatrix in rows Lo to Hi
  FOR I = Lo TO Hi
    Scale(I) = 1
  NEXT I

  ' Iterative loop for norm reduction
  DO
    Conv = True

    FOR I = Lo TO Hi
      C = 0
      R = 0
      FOR J = Lo TO Hi
        IF J <> I THEN 
          C = C + ABS(A(J, I))
          R = R + ABS(A(I, J))
        END IF  
      NEXT J

      ' Guard against zero C or R due to underflow
      IF (C <> 0) AND (R <> 0) THEN
        G = R / RADIX
        F = 1
        S = C + R

        DO WHILE C < G
          F = F * RADIX
          C = C * RADIX2
        LOOP

        G = R * RADIX

        DO WHILE C >= G
          F = F / RADIX
          C = C / RADIX2
        LOOP

        ' Now balance
        IF (C + R) / F < .95 * S THEN
          G = 1 / F
          Scale(I) = Scale(I) * F
          Conv = False
          FOR J = Lo TO Ub
            A(I, J) = A(I, J) * G
          NEXT J
          FOR J = Lb TO Hi
            A(J, I) = A(J, I) * F
          NEXT J
        END IF
      END IF
    NEXT I
  LOOP UNTIL Conv
END SUB
