'************************************************************************
'*          Eigenvalues of a real nonsymmetric square matrix            *
'*                      using the HQR method                            *
'*                                                                      *
'*                              Basic version by J-P Moreau, Paris      *
'*                                      (www.jpmoreau.fr)               *
'* -------------------------------------------------------------------- *
'* SAMPLE RUN:                                                          *
'*                                                                      *
'* Input matrix (symmetric for comparison with Jacobi):                 *
'*                                                                      *
'*        1,  2,  3, -7,   12                                           *
'*        2,  4,  7,  3,   -1                                           *
'*        3,  7, 10,  8,    4                                           *
'*       -7,  3,  8, -0.75,-9                                           *
'*       12, -1,  4, -9,   10                                           *
'*                                                                      *
'*                                                                      *
'* Output file (test_hqr.lst):                                          *
'*                                                                      *
'* =============================================================        *
'*   TESTING THE QR ALGORITHM TO FIND THE REAL OR COMPLEX               *
'*     EIGENVALUES OF A REAL SQUARE NON SYMMETRIC MATRIX                *
'*            IN DOUBLE PRECISION (16 digits)                           *
'* =============================================================        *
'*                                                                      *
'* GIVEN MATRIX A:                                                      *
'*                                                                      *
'*    1.00    2.00    3.00   -7.00   12.00                              *
'*    2.00    4.00    7.00    3.00   -1.00                              *
'*    3.00    7.00   10.00    8.00    4.00                              *
'*   -7.00    3.00    8.00   -0.75   -9.00                              *
'*   12.00   -1.00    4.00   -9.00   10.00                              *
'*                                                                      *
'* MATRIX A AFTER BALANC:                                               *
'*                                                                      *
'*    1.0000    2.0000    3.0000   -7.0000   12.0000                    *
'*    2.0000    4.0000    7.0000    3.0000   -1.0000                    *
'*    3.0000    7.0000   10.0000    8.0000    4.0000                    *
'*   -7.0000    3.0000    8.0000   -0.7500   -9.0000                    *
'*   12.0000   -1.0000    4.0000   -9.0000   10.0000                    *
'*                                                                      *
'* ElmHes done.                                                         *
'*                                                                      *
'* MATRIX A IN HESSENBERG FORM:                                         *
'*                                                                      *
'*    1.000000   17.166667   -9.738494    2.674367    3.000000          *
'*   12.000000   16.083333   -9.322176   -0.100844    4.000000          *
'*    0.250000    3.319444  -11.372036    4.739487   10.333333          *
'*   -0.583333   -0.307531  -11.556061    9.893547   15.715481          *
'*    0.166667   -0.907950    0.224789    8.506681    8.645156          *
'*                                                                      *
'*                                                                      *
'* HQR done.                                                            *
'*                                                                      *
'* EIGEN VALUES, REAL AND IMAGINARY PARTS:                              *
'*                                                                      *
'*  -10.486545    0.000000                                              *
'*   -7.774580    0.000000                                              *
'*   23.755955    0.000000                                              *
'*   18.291821    0.000000                                              *
'*    0.463350    0.000000                                              *
'*                                                                      *
'* End of file Test_hqr.lst.                                            *
'*                                                                      *
'* -------------------------------------------------------------------- *
'* Subroutines used:                                                    *
'*   1000  Balance:  Balances the elements of the input matrix          *
'*   2000  ElmHes :  Reduction of input matrix to Hessenberg form       *
'*   3000  HQR    :  Finds all the eigenvalues of an upper triangle     *
'                    Hessenberg matrix.                                 *
'************************************************************************
DEFINT I-N
DEFDBL A-H, O-Z

ZERO = 0#
n = 5        'size of square matrix

DIM A(n, n)

A(1, 1) = 1#: A(1, 2) = 2#: A(1, 3) = 3#: A(1, 4) = -7#: A(1, 5) = 12#
A(2, 1) = 2#: A(2, 2) = 4#: A(2, 3) = 7#: A(2, 4) = 3#: A(2, 5) = -1#
A(3, 1) = 3#: A(3, 2) = 7#: A(3, 3) = 10#: A(3, 4) = 8#: A(3, 5) = 4#
A(4, 1) = -7#: A(4, 2) = 3#: A(4, 3) = 8#: A(4, 4) = -.75#: A(4, 5) = -9#
A(5, 1) = 12#: A(5, 2) = -1#: A(5, 3) = 4#: A(5, 4) = -9#: A(5, 5) = 10#

DIM d(n), wr(n), wi(n)

OPEN "test_hqr.lst" FOR OUTPUT AS #1

CLS
PRINT #1,
PRINT #1, " ==========================================================="
PRINT #1, "   TESTING THE QR ALGORITHM TO FIND THE REAL OR COMPLEX"
PRINT #1, "     EIGENVALUES OF A REAL SQUARE NON SYMMETRIC MATRIX"
PRINT #1, "            IN DOUBLE PRECISION (16 digits)"
PRINT #1, " ==========================================================="
PRINT #1,

PRINT #1, " GIVEN MATRIX A:"
PRINT #1,

FOR i = 1 TO n
  FOR j = 1 TO n
    PRINT #1, USING " ####.##"; A(i, j);
  NEXT j
  PRINT #1,
NEXT i
PRINT #1,

low = 1
ihi = n

'Balanc(A,n, d, low, hi)
GOSUB 1000

PRINT
PRINT " Balanc done."
PRINT #1,
PRINT #1, " MATRIX A AFTER BALANC:"
PRINT #1,
FOR i = 1 TO n
  FOR j = 1 TO n
    PRINT #1, USING " ####.####"; A(i, j);
  NEXT j
  PRINT #1,
NEXT i
PRINT #1,

'ElmHes(A,n)
GOSUB 2000

PRINT
PRINT " ElmHes done."

PRINT #1,
PRINT #1, " MATRIX A IN HESSENBERG FORM:"
PRINT #1,
FOR i = 1 TO n
    FOR j = 1 TO n
      PRINT #1, USING " ####.######"; A(i, j);
    NEXT j
    PRINT #1,
NEXT i
PRINT #1,

'HQR_MR(A,n,wr,wi)
GOSUB 3000

PRINT
PRINT " HQR done."

PRINT #1,
PRINT #1, " EIGEN VALUES, REAL AND IMAGINARY PARTS:"
PRINT #1,

FOR i = 1 TO n
  PRINT #1, USING "  ###.######  ###.######"; wr(i); wi(i)
NEXT i
PRINT #1,
PRINT #1, " End of file Test_hqr.lst."
CLOSE #1
PRINT
PRINT " Results in file test_hqr.lst."
END 'of main program

'----------------------------------------------------------------------
900 'Subroutine Exc(m : INTEGER)  internal subroutine of Balanc
  d(m) = 1# * j
  IF j <> m THEN
    FOR ii = 1 TO k
      f = A(ii, j)
      A(ii, j) = A(ii, m)
      A(ii, m) = f
    NEXT ii
    FOR ii = l TO n
      f = A(j, ii)
      A(j, ii) = A(m, ii)
      A(m, ii) = f
    NEXT ii
  END IF
RETURN
'--------------------------------------------------------------

'***********************************************************************
'* Balanc:   From Pascal version translated from ALGOL by J-P Dumont   *
'*                                  Basic version by J-P Moreau        *
'* ------------------------------------------------------------------- *
'* Ref.: "Handbook FOR automatic computation Volume II Linear Algebra  *
'*        J.H. Wilkinson et C.Reinsch, Springer Verlag, 1971"          *
'* ------------------------------------------------------------------- *
'*   This procedure balances the elements of a non symmetric square    *
'* matrix. Given a matrix n x n stored in a table of size n x n, that  *
'* matrix is replaced by a balanced matrix having the same eigenvalues *
'* A symmetric matrix is not affected by this procedure.               *
'* ------------------------------------------------------------------- *
'* Inputs:                                                             *
'*       n       Order of non symmetric square matrix A (n x n)        *
'*       A       Table n*n storing the  elements of A                  *
'* ------------------------------------------------------------------- *
'* Outputs:                                                            *
'*       A       The original A matrix is replaced by the balanced     *
'*               matrix. However, it is possible to recover exactly    *
'*               the original matrix.                                  *
'* low,ihi       Two integers such as A(i,j) = 0 if                    *
'*               (1) i > j  AND                                        *
'*               (2) j =1,...,low-1 OR i = ihi+1,...,n                 *
'*       D       Vector(1..n) containing enough information to trace   *
'*               the done permutations and the used scale factors.     *
'***********************************************************************
1000 'Subroutine Balanc
'LABELS: Iteration(1100),L1(1200),L2(1300),L3(1400),L4(1500)

b = 2    'Floating point basis of used CPU

'i,j,k,l        : INTEGER
'b2,c,f,g,r,s   : DOUBLE
'noconv         : INTEGER (0 or 1)

b2 = b * b
l = 1
k = n
'Search lines isolating an eigenvalue and shift downwards
1200 FOR j = k TO 1 STEP -1
  r = ZERO
  FOR i = 1 TO j - 1
    r = r + ABS(A(j, i))
  NEXT i
  FOR i = j + 1 TO k
    r = r + ABS(A(j, i))
  NEXT i
  IF (r = 0) THEN
    m = k: GOSUB 900 'Exc(k)
    k = k - 1
    GOTO 1200
  END IF
NEXT j
'Search columns isolating an eigenvalue and shift to the left
1300 FOR j = l TO k
  c = ZERO
  FOR i = l TO j - 1
    c = c + ABS(A(i, j))
  NEXT i
  FOR i = j + 1 TO k
    c = c + ABS(A(i, j))
  NEXT i
  IF (c = 0) THEN
    m = l: GOSUB 900 'Exc(l)
    l = l + 1
    GOTO 1300
  END IF
NEXT j
'Now balance submatrix from ligne l to line k
  low = l
  ihi = k
  FOR i = 1 TO k
    d(i) = 1#
  NEXT i
1100  noconv = 0
  FOR i = l TO k
    c = ZERO
    r = c
    FOR j = l TO i - 1
      c = c + ABS(A(j, i))
      r = r + ABS(A(i, j))
    NEXT j
    FOR j = i + 1 TO k
      c = c + ABS(A(j, i))
      r = r + ABS(A(i, j))
    NEXT j
    g = r / b
    f = 1#
    s = c + r
1400 IF c < g THEN
      f = f * b
      c = c * b2
      GOTO 1400
    END IF
    g = r * b
1500 IF c >= g THEN
      f = f / b
      c = c / b2
      GOTO 1500
    END IF

'Balancing the elements of submatrix
    IF (c + r) / f < .95 * s THEN
      g = 1# / f
      d(i) = d 'i) * f
      noconv = 1
      FOR j = l TO n
        A(i, j) = A(i, j) * g
      NEXT j
      FOR j = 1 TO k
        A(j, i) = A(j, i) * f
      NEXT j
    END IF
  NEXT i
  IF noconv = 1 THEN GOTO 1100

RETURN 'from Balanc


'{-------------------------Documentation-------------------------------}
'{               VAR A: Square_Matrix; n: INTEGER);                    }
'{ Reduction of a non symmetric real matrix to Hessenberg form by the  }
'{ elimination method. Matrix A (n x n), stored in a table of size     }
'{ n x n, is replaced by a superior Hessenberg matrix having the same  }
'{ eigenvalues. It is recommanded to call previously the Balanc        }
'{ subroutine. In output, the Hessenberg matrix has elements A(i,j)    }
'{ with i<=j+1. The elements for i>j+1, that in theory equal zero, are }
'{ actually filled with random (not used) values.                      }
'{---------------------------------------------------------------------}
'{ From Pascal version by J.P.DUMONT - Extracted from reference:       }
'{       "William H.PRESS, Brian P.FLANNERY, Saul A.TEUKOLSKY AND      }
'{        William T.VETTERLING                                         }
'{              N U M E R I C A L  R E C I P E S                       }
'{              The Art OF Scientific Computing                        }
'{              CAMBRIDGE UNIVERSITY PRESS 1987"                       }
'{                                       Basic version by J-P Moreau   }
'{/////////////////////////////////////////////////////////////////////}
2000 'Subroutine ElmHes
  IF n > 2 THEN
    FOR m = 2 TO n - 1
      x = ZERO
      i = m
      FOR j = m TO n
        IF ABS(A(j, m - 1)) > ABS(x) THEN
          x = A(j, m - 1)
          i = j
        END IF 'IF Abs
      NEXT j 'FOR j= m TO n
      IF i <> m THEN
        FOR j = m - 1 TO n
          y = A(i, j)
          A(i, j) = A(m, j)
          A(m, j) = y
        NEXT j
        FOR j = 1 TO n
          y = A(j, i)
          A(j, i) = A(j, m)
          A(j, m) = y
        NEXT j
      END IF 'IF i <> m
      IF x <> ZERO THEN
        FOR i = m + 1 TO n
          y = A(i, m - 1)
          IF y <> ZERO THEN
            y = y / x
            A(i, m - 1) = y
            FOR j = m TO n
              A(i, j) = A(i, j) - y * A(m, j)
            NEXT j
            FOR j = 1 TO n
              A(j, m) = A(j, m) + y * A(j, i)
            NEXT j
          END IF  'IF y
        NEXT i  'FOR i
      END IF  'IF x
    NEXT m  'FOR m
  END IF 'if n>2
RETURN 'from ElmHes

'---------------------------------------------------------------------
2900 'FUNCTION Sign(aa,bb) used by HQR
  IF bb < 0 THEN
    Sign = -ABS(aa)
  ELSE
    Sign = ABS(aa)
  END IF
RETURN
'---------------------------------------------------------------------

'{-------------------------Documentation--------------------------------
'                   VAR A : Square_Matrix; n: INTEGER;
'                   VAR wr,wi : Real_Vector);}
'{/////////////////////////////////////////////////////////////////////}
'{               QR algorithm for real Hessenberg matrices             }
'{                                                                     }
'{   This procedure finds all the eigenvalues of an upper triangle     }
'{ Hessenberg matrix of size n x n, stored in a table of size n x n.   }
'{ The input matrix A has been previously produced by the Helmes       }
'{ subroutine. In output this A matrix is destroyed. Tables wr and wi  }
'{ respectively return real parts and imaginary parts of the eigen-    }
'{ values.                                                             }
'{/////////////////////////////////////////////////////////////////////}
'{ Reference (Pascal version):                                         }
'{        William H.PRESS, Brian P.FLANNERY, Saul A.TEUKOLSKY AND      }
'{        William T.VETTERLING                                         }
'{              N U M E R I C A L  R E C I P E S                       }
'{              The Art OF Scientific Computing                        }
'{              CAMBRIDGE UNIVERSITY PRESS 1987                        }
'{                                       Basic version by J-P Moreau   }
'{/////////////////////////////////////////////////////////////////////}
3000 'Subroutine HQR
'LABELS: 3100,3200,3300,3400

  itsmx = 30

  anorm = ABS(A(1, 1))
  FOR i = 2 TO n
    FOR j = i - 1 TO n
      anorm = anorm + ABS(A(i, j))
    NEXT j
  NEXT i
  nn = n
  t = ZERO

3100 its = 0
3200 FOR l = nn TO 2 STEP -1
    s = ABS(A(l - 1, l - 1)) + ABS(A(l, l))
    IF (s = 0!) THEN s = anorm
    IF ((ABS(A(l, l - 1)) + s) = s) THEN GOTO 3300
  NEXT l
  l = 1
3300 x = A(nn, nn)
  IF (l = nn) THEN
    wr(nn) = x + t
    wi(nn) = ZERO
    nn = nn - 1
  ELSE
    y = A(nn - 1, nn - 1)
    w = A(nn, nn - 1) * A(nn - 1, nn)
    IF (l = nn - 1) THEN
      p = .5# * (y - x)
      q = p * p + w
      z = SQR(ABS(q))
      x = x + t
      IF q >= ZERO THEN
        aa = z: bb = p: GOSUB 2900
        z = p + Sign
        wr(nn) = x + z
        wr(nn - 1) = wr(nn)
        IF z <> ZERO THEN wr(nn) = x - w / z
        wi(nn) = ZERO
        wi(nn - 1) = ZERO
      ELSE
        wr(nn) = x + p
        wr(nn - 1) = wr(nn)
        wi(nn) = z
        wi(nn - 1) = -z
      END IF
      nn = nn - 2
    ELSE
      IF its = itsmx THEN
        PRINT " Pause in HQR subroutine."
        PRINT " Too many iterations!"
        INPUT r$
      END IF
      IF (its = 10) OR (its = 20) THEN
        t = t + x
        FOR i = 1 TO nn
          A(i, i) = A(i, i) - x
        NEXT i
        s = ABS(A(nn, nn - 1)) + ABS(A(nn - 1, nn - 2))
        x = .75# * s
        y = x
        w = -.4375# * s * s
      END IF
      its = its + 1
      FOR m = nn - 2 TO 1 STEP -1
        z = A(m, m)
        r = x - z
        s = y - z
        p = (r * s - w) / A(m + 1, m) + A(m, m + 1)
        q = A(m + 1, m + 1) - z - r - s
        r = A(m + 2, m + 1)
        s = ABS(p) + ABS(q) + ABS(r)
        p = p / s
        q = q / s
        r = r / s
        IF m = 1 THEN GOTO 3400
        u = ABS(A(m, m - 1)) * (ABS(q) + ABS(r))
        v = ABS(p) * (ABS(A(m - 1, m - 1)) + ABS(z) + ABS(A(m + 1, m + 1)))
        IF u + v = v THEN GOTO 3400
      NEXT m
3400  FOR i = m + 2 TO nn
        A(i, i - 2) = ZERO
        IF i <> (m + 2) THEN A(i, i - 3) = ZERO
      NEXT i
      FOR k = m TO nn - 1
        IF k <> m THEN
          p = A(k, k - 1)
          q = A(k + 1, k - 1)
          r = ZERO
          IF k <> (nn - 1) THEN r = A(k + 2, k - 1)
          x = ABS(p) + ABS(q) + ABS(r)
          IF x <> ZERO THEN
            p = p / x
            q = q / x
            r = r / x
          END IF
        END IF
        aa = SQR(p * p + q * q + r * r): bb = p: GOSUB 2900
        s = Sign
        IF s <> ZERO THEN
          IF k = m THEN
            IF l <> m THEN A(k, k - 1) = -A(k, k - 1)
          ELSE
            A(k, k - 1) = -s * x
          END IF
          p = p + s
          x = p / s
          y = q / s
          z = r / s
          q = q / p
          r = r / p
          FOR j = k TO nn
            p = A(k, j) + q * A(k + 1, j)
            IF k <> (nn - 1) THEN
              p = p + r * A(k + 2, j)
              A(k + 2, j) = A(k + 2, j) - p * z
            END IF
            A(k + 1, j) = A(k + 1, j) - p * y
            A(k, j) = A(k, j) - p * x
          NEXT j
          IF nn < k + 3 THEN
            mmin = nn
          ELSE
            mmin = k + 3
          END IF
          FOR i = 1 TO mmin
            p = x * A(i, k) + y * A(i, k + 1)
            IF k <> (nn - 1) THEN
              p = p + z * A(i, k + 2)
              A(i, k + 2) = A(i, k + 2) - p * r
            END IF
            A(i, k + 1) = A(i, k + 1) - p * q
            A(i, k) = A(i, k) - p
          NEXT i
        END IF
      NEXT k
      GOTO 3200
    END IF
  END IF
  IF nn >= 1 THEN GOTO 3100
RETURN 'from Subroutine HQR

'End of file Test_hqr.bas