'*=====================================================================
'*            Calculation of eigenvalues and eigenvectors             *
'*               of a real non symmetric square matrix                *
'*                         by QR algorithm                            *
'* ------------------------------------------------------------------ *
'* SAMPLE RUN:                                                        *
'*                                                                    *
'* Input file Hqr5.dat contains:                                      *
'*                                                                    *
'* 5                                                                  *
'*  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 Hqr5.lst contains:                                     *
'*                                                                    *
'* --------------------------------------------------------------     *
'*  Eigenvalues and Eigenvectors by QR algorithm                      *
'* --------------------------------------------------------------     *
'*  Dimension of the input matrix =  5                                * 
'*                                                                    *
'*  Input matrix:                                                     *
'*   1.000000  2.000000  3.000000 -7.000000 12.000000                 *
'*   2.000000  4.000000  7.000000  3.000000 -1.000000                 *
'*   3.000000  7.000000 10.000000  8.000000  4.000000                 *
'*  -7.000000  3.000000  8.000000 -0.750000 -9.000000                 *
'*  12.000000 -1.000000  4.000000 -9.000000 10.000000                 *
'*                                                                    *
'*  Normalized Eigenvectors (in columns):                             *
'*   0.467172  1.000000  0.705820  0.093053  0.208812                 *
'*  -0.007947 -0.326100 -0.012677  0.594911  1.000000                 *
'*  -0.507827  0.209620  0.131486  1.000000 -0.478027                 *
'*   1.000000 -0.147689 -0.527499  0.455666 -0.275000                 *
'*   0.264432 -0.815422  1.000000  0.050740 -0.216915                 *
'*                                                                    *
'*  Eigenvalues:               Iterations:                            *
'*  -10.486545 +   0.0000000 i      0                                 *
'*   -7.774580 +   0.0000000 i      4                                 *
'*   23.755955 +   0.0000000 i     -4                                 *
'*   18.291821 +   0.0000000 i      5                                 *
'*    0.463350 +   0.0000000 i     -5                                 *
'*                                                                    *
'*  Check sum =  2.783592651625817D-13                                *
'*  (must be approximately 0).                                        *
'* --------------------------------------------------------------     *
'* Reference: "Numerical Algorithms with C By G. Engeln-Mueller and   *
'*             F. Uhlig, Springer-Verlag, 1996" [BIBLI 11].           *
'*                                                                    *
'*                               Basic version by J-P Moreau, Paris.  *
'*                                        (www.jpmoreau.fr)           *
'*====================================================================*
'PROGRAM Test_HQR
DEFDBL A-H, O-Z
DEFINT I-N

OPTION BASE 0                  'start index from zero (translation from C)

'global constants
NSIZE = 4: MAXIT = 50
ONE = 1#: TWO = 2#: ZERO = 0#: XMACHEPS = 2.22E-16

DIM xmat(NSIZE, NSIZE)         'input real square matrix
DIM A(NSIZE, NSIZE)            'copy of the input matrix (for check sum)
DIM ev(NSIZE, NSIZE)           'Eigenvectors (if ivec <> 0)
DIM wr(NSIZE)                  'Eigenvalues (Real part)
DIM wi(NSIZE)                  'Eigenvalues (Imaginary parts)
DIM scal(NSIZE)                'see subroutine 500
DIM icnt(NSIZE)                'Iteration counter
DIM iperm(NSIZE)               'Permutation vector

'integer irc                    Return Code
'integer ivec                   flag for eigenvectors (=0 -> none)
'integer ievnorm                flag for normalization of Eigenvectors
			       ' 0 = no normalization
'double  xnorm                  norm of input matrix (must be <> 0)

  CLS
  ivec = 1: ievnorm = 1

' open input/output files
  OPEN "hqr5.dat" FOR INPUT AS #1       'read mode
  OPEN "hqr5.lst" FOR OUTPUT AS #2      'print mode

  PRINT #2, "--------------------------------------------------------------"
  PRINT #2, " Eigenvalues and Eigenvectors by QR algorithm"
  PRINT #2, "--------------------------------------------------------------"

' read size of matrix from input file
  INPUT #1, n

  IF n < 1 THEN
    PRINT #2, " Dimension must be > 0."
    GOTO 10   'go to ending section
  END IF

' read matrix from input file
  FOR i = 0 TO n - 1
    FOR j = 0 TO n - 1
      INPUT #1, xmat(i, j)
      A(i, j) = xmat(i, j)
    NEXT j
  NEXT i
  CLOSE #1

' print data to output file
  PRINT #2, " Dimension of the input matrix = "; n
  PRINT #2,
  PRINT #2, " Input matrix:"
  GOSUB 1400  'print matrix xmat
  PRINT #2,

  'see subroutine 2000
  GOSUB 2000 'call eigen(ivec, ievnorm, n, xmat, ev, wr, wi, icnt, irc)

  IF irc <> 0 THEN
    PRINT #2, " Error in procedure Eigen, rc="; irc
    GOTO 10   'go to ending section
  END IF

' If ivec <> 0, print eigenvectors
  IF ivec <> 0 THEN
    IF ievnorm <> 0 THEN
      PRINT #2, " Normalized Eigenvectors (in columns):"
    ELSE
      PRINT #2, " Not normalized Eigenvectors (in columns):"
    END IF
    GOSUB 1410  'print matrix ev
  END IF

  PRINT #2,
  PRINT #2, " Eigenvalues:               Iterations:"
  FOR i = 0 TO n - 1
    PRINT #2, USING " ###.###### + ###.####### i     ##"; wr(i); wi(i); icnt(i)
  NEXT i

' Check result: sum of L1 norms of Matrix*Eigenvector - Eigenvalue*Eigenvector
' (this must be nearly 0).
  IF ivec <> 0 THEN
    xnorm = ZERO: k = 0
    WHILE (k <= n - 1)
      IF wi(k) = ZERO THEN
	FOR i = 0 TO n - 1
	  w = ZERO
	  FOR j = 0 TO n - 1
	    w = w + A(i, j) * ev(j, k)
	  NEXT j
	  w = w - wr(k) * ev(i, k)
	  xnorm = xnorm + ABS(w)
	NEXT i
      ELSE
	FOR i = 0 TO n - 1
	  w = ZERO
	  FOR j = 0 TO n - 1
	    w = w + A(i, j) * ev(j, k)
	  NEXT j
	  w = w - wr(k) * ev(i, k) - wi(k) * ev(i, k + 1)
	  V = ZERO
	  FOR j = 0 TO n - 1
	    V = V + A(i, j) * ev(j, k + 1)
	  NEXT j
	  V = V - wr(k) * ev(i, k + 1) + wi(k) * ev(i, k)
	  xnorm = xnorm + 2# * SQR(V * V + w * w)
	NEXT i
	k = k + 1
      END IF
      k = k + 1
    WEND 'while k<=n-1
    PRINT #2, " "
    PRINT #2, " Check sum = "; xnorm
    PRINT #2, " (must be approximately 0)."
  END IF

10 PRINT #2, "--------------------------------------------------------------"
  CLOSE #2
  PRINT
  PRINT " Results in hqr5.lst."
  PRINT

END 'of main program

'* ------------------------ QR Subroutines -------------------------- *
'* Reference: "Numerical Algorithms with C by G. Engeln-Mueller and   *
'*             F. Uhlig, Springer-Verlag, 1996"                       *
'*                                                                    *
'*                               Basic version by J-P Moreau, Paris.  *
'* ------------------------------------------------------------------ *

'Subroutine balance        "balance a matrix
'                (n,       "size of matrix
'                 xmat,    "input matrix
'                 scal,    "Scaling data
'                 low,     "first relevant row index
'                 ihigh    "last relevant row index
'                )
'integer ihigh
'*====================================================================*
'*                                                                    *
'* balance balances the matrix xmat so that the rows with zero entries*
'* off the diagonal are isolated and the remaining columns and rows   *
'* are resized to have one norm close to 1.                           *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Input parameters:                                                *
'*   ================                                                 *
'*      n        integer;  ( n > 0 )                                  *
'*               Dimension of xmat                                    *
'*      xmat     n x n input matrix                                   *
'*                                                                    *
'*   Output parameters:                                               *
'*   ==================                                               *
'*      xmat     n x n scaled matrix                                  *
'*      low      integer;                                             *
'*      ihigh    integer;                                             *
'*               the rows 0 to low-1 and those from ihigh to n-1      *
'*               contain isolated eigenvalues (only nonzero entry on  *
'*               the diagonal)                                        *
'*      scal     vector of size n                                     *
'*               the vector scal contains the isolated eigenvalues in *
'*               the positions 0 to low-1 and ihigh to n-1, its other *
'*               components contain the scaling factors for           *
'*               transforming mat.                                    *
'*                                                                    *
'*====================================================================*
500 ibasis = 2       'Floating point basis of used CPU
'double b2, r, c, f, g, s

  b2 = ibasis * ibasis
  m = 0
  k = n - 1
  iter = 1
  WHILE (iter = 1)
    iter = 0
    FOR j = k TO 0 STEP -1
      r = ZERO
      FOR i = 0 TO k
	IF i <> j THEN r = r + ABS(xmat(j, i))
      NEXT i
      IF r = ZERO THEN
	scal(k) = j
	IF j <> k THEN
	  FOR i = 0 TO k
	    'CALL SWAP(xmat(i, j), xmat(i, k))
	    temp = xmat(i, k): xmat(i, k) = xmat(i, j): xmat(i, j) = temp
	  NEXT i
	  FOR i = m TO n - 1
	    'CALL SWAP(xmat(j, i), xmat(k, i))
	    temp = xmat(k, i): xmat(k, i) = xmat(j, i): xmat(j, i) = temp
	  NEXT i
	END IF
	k = k - 1
	iter = 1
      END IF
    NEXT j
  WEND 'while iter=1

  iter = 1
  WHILE (iter = 1)
    iter = 0
    FOR j = m TO k
      c = ZERO
      FOR i = m TO k
	IF i <> j THEN c = c + ABS(xmat(i, j))
      NEXT i
      IF c = ZERO THEN
	scal(m) = j
	IF j <> m THEN
	  FOR i = 0 TO k
	    'CALL SWAP(xmat(i, j), xmat(i, m))
	    temp = xmat(i, m): xmat(i, m) = xmat(i, j): xmat(i, j) = temp
	  NEXT i
	  FOR i = m TO n - 1
	    'CALL SWAP(xmat(j, i), xmat(m, i))
	    temp = xmat(m, i): xmat(m, i) = xmat(j, i): xmat(j, i) = temp
	  NEXT i
	END IF
	m = m + 1
	iter = 1
      END IF
    NEXT j
  WEND 'while iter=1

  low = m
  ihigh = k
  FOR i = m TO k
    scal(i) = ONE
  NEXT i

  iter = 1
  WHILE (iter = 1)
    iter = 0
    FOR i = m TO k
      c = ZERO: r = ZERO
      FOR j = m TO k
	IF j <> i THEN
	  c = c + ABS(xmat(j, i))
	  r = r + ABS(xmat(i, j))
	END IF
      NEXT j
      g = r / ibasis
      f = ONE
      s = c + r

      WHILE (c < g)
	f = f * ibasis
	c = c * b2
      WEND

      g = r * ibasis
      WHILE (c >= g)
	f = f / ibasis
	c = c / b2
      WEND

      IF ((c + r) / f < .95 * s) THEN
	g = ONE / f
	scal(i) = scal(i) * f
	iter = 1
	FOR j = m TO n - 1
	  xmat(i, j) = xmat(i, j) * g
	NEXT j
	FOR j = 0 TO k
	  xmat(j, i) = xmat(j, i) * f
	NEXT j
      END IF
    NEXT i
  WEND 'while iter=1
RETURN

'Subroutine balback           "reverse balancing ...........
'                  (n,        "Dimension of matrix .........
'                   low,      "first nonzero row ...........
'                   ihigh,    "last nonzero row ............
'                   scal,     "Scaling data ................
'                   ev        "Eigenvectors ................
'                  )
'integer ihigh
'*====================================================================*
'*                                                                    *
'*  balback reverses the balancing of balance for the eigenvactors.   *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Input parameters:                                                *
'*   ----------------                                                 *
'*      n        integer:  ( n > 0 )                                  *
'*               Dimension of xmat                                    *
'*      low      integer:                                             *
'*      ihigh    integer:   see balance                               *
'*      ev       n x n matrix of eigenvectors, as computed in  qr2    *
'*      scal     vector of size n:                                    *
'*               Scaling data from  balance                           *
'*                                                                    *
'*   Output parameter:                                                *
'*   ----------------                                                 *
'*      ev       n x n matrix:                                        *
'*               Non-normalized eigenvectors of the original matrix   *
'*                                                                    *
'*====================================================================*
600 'double s

  FOR i = low TO ihigh
    s = scal(i)
    FOR j = 0 TO n - 1
      ev(i, j) = ev(i, j) * s
    NEXT j
  NEXT i

  FOR i = low - 1 TO 0 STEP -1
    k = INT(scal(i))
    IF k <> i THEN
      FOR j = 0 TO n - 1
	'CALL SWAP(ev(i, j), ev(k, j))
	temp = ev(k, j): ev(k, j) = ev(i, j): ev(i, j) = temp
      NEXT j
    END IF
  NEXT i

  FOR i = ihigh + 1 TO n - 1
    k = INT(scal(i))
    IF k <> i THEN
      FOR j = 0 TO n - 1
	'CALL SWAP(ev(i, j), ev(k, j))
	temp = ev(k, j): ev(k, j) = ev(i, j): ev(i, j) = temp
      NEXT j
    END IF
  NEXT i
RETURN

'Subroutine elmhes         "reduce matrix to upper Hessenberg form
'                (n,       "Dimension of matrix .........
'                 low,     "first nonzero row ...........
'                 ihigh,   "last nonzero row ............
'                 xmat,    "input/output matrix .........
'                 iperm    "Permutation vector ..........
'                )
'integer ihigh
'*====================================================================*
'*                                                                    *
'*  elmhes transforms the matrix mat to upper Hessenberg form.        *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Input parameters:                                                *
'*   ================                                                 *
'*      n        integer:  ( n > 0 )                                  *
'*               Dimension of mat                                     *
'*      low      integer:                                             *
'*      ihigh    integer: see  balance                                *
'*      xmat     n x n matrix                                         *
'*                                                                    *
'*   Output parameter:                                                *
'*   =================                                                *
'*      xmat     n x n matrix:                                        *
'*               upper Hessenberg matrix: additional information on   *
'*               the transformation is stored in the lower triangle   *
'*      iperm    integer vector of size n:                            *
'*               Permutation vector for elmtrans                      *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Subroutine used:   None                                          *
'*   ===============                                                  *
'*                                                                    *
'*====================================================================*
700 'double x, y
  FOR m = low + 1 TO ihigh - 1
    i = m
    x = ZERO
    FOR j = m TO ihigh
      IF ABS(xmat(j, m - 1)) > ABS(x) THEN
	x = xmat(j, m - 1)
	i = j
      END IF
    NEXT j

    iperm(m) = i
    IF i <> m THEN
      FOR j = m - 1 TO n - 1
	'CALL SWAP(xmat(i, j), xmat(m, j))
	temp = xmat(m, j): xmat(m, j) = xmat(i, j): xmat(i, j) = temp
      NEXT j
      FOR j = 0 TO ihigh
	'CALL SWAP(xmat(j, i), xmat(j, m))
	temp = xmat(j, m): xmat(j, m) = xmat(j, i): xmat(j, i) = temp
      NEXT j
    END IF

    IF x <> ZERO THEN
      FOR i = m + 1 TO ihigh
	y = xmat(i, m - 1)
	IF y <> ZERO THEN
	  y = y / x
	  xmat(i, m - 1) = y
	  FOR j = m TO n - 1
	    xmat(i, j) = xmat(i, j) - y * xmat(m, j)
	  NEXT j
	  FOR j = 0 TO ihigh
	    xmat(j, m) = xmat(j, m) + y * xmat(j, i)
	  NEXT j
	END IF
      NEXT i
    END IF 'x <> ZERO
  NEXT m
RETURN

'Subroutine elmtrans        "initialize eigenvectors......
'                  (n,      "Dimension of matrix .........
'                   low,    "first nonzero row ...........
'                   high,   "last nonzero row ............
'                   xmat,   "input matrix ................
'                   iperm,  "row permutations ............
'                   ev      "Eigenvectors matrix .........
'                   )
'integer ihigh, iperm(0:n)
'*====================================================================*
'*                                                                    *
'*  elmtrans initializes the eigenvectors matrix.                     *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Input parameters:                                                *
'*   ================                                                 *
'*      n        integer:  ( n > 0 )                                  *
'*               Dimension of xmat and ev                             *
'*      low      integer:                                             *
'*      ihigh    integer: see  balance                                *
'*      xmat     n x n input matrix                                   *
'*      iperm    Integer vector of size n:                            *
'*               Permutation data from  elmhes                        *
'*                                                                    *
'*   Output parameter:                                                *
'*   ================                                                 *
'*      ev       n x n matrix:                                        *
'*               Eigenvectors matrix                                  *
'*                                                                    *
'*====================================================================*
800 'begin elmtrans
  FOR i = 0 TO n - 1
    FOR k = 0 TO n - 1
      ev(i, k) = ZERO
    NEXT k
    ev(i, i) = ONE
  NEXT i

  FOR i = ihigh - 1 TO low + 1 STEP -1
    j = iperm(i)
    FOR k = i + 1 TO ihigh
      ev(k, i) = xmat(k, i - 1)
    NEXT k
    IF i <> j THEN
      FOR k = i TO ihigh
	ev(i, k) = ev(j, k)
	ev(j, k) = ZERO
      NEXT k
      ev(j, i) = ONE
    END IF
  NEXT i
RETURN

'Subroutine Comdiv (        "Complex division ................
'            ar,            "Real part of numerator ..........
'            ai,            "Imaginary part of numerator .....
'            br,            "Real part of denominator ........
'            bi,            "Imaginary part of denominator ...
'            cr,            "Real part of quotient ...........
'            ci,            "Imaginary part of quotient ......
'            irc            "return code .....................
'           )
'double ar,ai,br,bi,cr,ci
'integer irc
'*====================================================================*
'*                                                                    *
'*  Complex division  c = a / b                                       *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Input parameters:                                                *
'*   ================                                                 *
'*      ar,ai    double;                                              *
'*               real, imaginary parts of numerator                   *
'*      br,bi    double:                                              *
'*               real, imaginary parts of denominator                 *
'*                                                                    *
'*   Output parameters:                                               *
'*   ==================                                               *
'*      cr,ci    double:                                              *
'*               real , imaginary parts of the quotient               *
'*                                                                    *
'*   Return code value:                                               *
'*   =================                                                *
'*      = 0      ok                                                   *
'*      = 1      division by 0                                        *
'*                                                                    *
'*====================================================================*
900 'double tmp
  IF br = ZERO AND bi = ZERO THEN
    irc = 1
    RETURN
  END IF
  IF ABS(br) > ABS(bi) THEN
    tmp = bi / br
    br = tmp * bi + br
    cr = (ar + tmp * ai) / br
    ci = (ai - tmp * ar) / br
  ELSE
    tmp = br / bi
    bi = tmp * br + bi
    cr = (tmp * ar + ai) / bi
    ci = (tmp * ai - ar) / bi
  END IF
 irc = 0
RETURN

'double FUNCTION Comabs(         "Complex absolute value ..........
'                       ar,      "Real part .......................
'                       ai       "Imaginary part ..................
'                      )
'double ar,ai
'*====================================================================*
'*                                                                    *
'*  Complex absolute value of a complex number                        *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Input parameters:                                                *
'*   ================                                                 *
'*      ar,ai    double;                                              *
'*               Real, imaginary parts of  a                          *
'*                                                                    *
'*   Return value :                                                   *
'*   =============                                                    *
'*      Absolute value of a (double)                                  *
'*                                                                    *
'*====================================================================*
1000 'begin Comabs
  IF ar = ZERO AND ai = ZERO THEN
    Comabs = ZERO
    RETURN
  END IF

  ar = ABS(ar)
  ai = ABS(ai)

  IF ai > ar THEN                                  'Switch  ai and ar
    'CALL SWAP(ai, ar)
    temp = ar: ar = ai: ai = temp
  END IF

  IF ai = ZERO THEN
    Comabs = ar
  ELSE
    Comabs = ar * SQR(ONE + ai / ar * ai / ar)
  END IF
RETURN

'Subroutine  hqrvec          "compute eigenvectors ......
'                  (n,       "Dimension of matrix .......
'                   low,     "first nonzero row .........
'                   ihigh,   "last nonzero row ..........
'                   xmat,    "upper Hessenberg matrix ...
'                   wr,      "Real parts of evalues .....
'                   wi,      "Imaginary parts of evalues
'                   ev,      "Eigenvectors ..............
'                   irc      "return code ...............
'                  )
'integer ihigh, irc
'*====================================================================*
'*                                                                    *
'* hqrvec computes the eigenvectors for the eigenvalues found in hqr2 *
'* (ev must have been initialized before by subroutine 800 Elmtrans ) *  
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Input parameters:                                                *
'*   ================                                                 *
'*      n        integer;  ( n > 0 )                                  *
'*               Dimension of xmat and ev, number of eigenvalues.     *
'*      low      int low;                                             *
'*      ihigh    integer; see  balance                                *
'*      xmat     n x n upper Hessenberg matrix                        *
'*      wr       vector of size n;                                    *
'*               Real parts of the n eigenvalues.                     *
'*      wi       vector of size n;                                    *
'*               Imaginary parts of the n eigenvalues.                *
'*                                                                    *
'*   Output parameter:                                                *
'*   ================                                                 *
'*      ev       n x n matrix, whose columns are the eigenvectors     *
'*                                                                    *
'*   Return value :                                                   *
'*   =============                                                    *
'*      =  0     all ok                                               *
'*      =  1     h is the zero matrix.                                *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Subroutine in use:                                               *
'*   =================                                                *
'*                                                                    *
'*      comdiv(): complex division                                    *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Constant used:    XMACHEPS                                       *
'*   =============                                                    *
'*                                                                    *
'*====================================================================*
1200 'integer ien

  r = ZERO: s = ZERO: z = ZERO: xnorm = ZERO

  FOR i = 0 TO n - 1                          'find norm of xmat
    FOR j = i TO n - 1
      xnorm = xnorm + ABS(xmat(i, j))
    NEXT j
  NEXT i

  IF xnorm = ZERO THEN
    irc = 1                                   'zero matrix
    RETURN
  END IF

  FOR ien = n - 1 TO 0 STEP -1                'transform back
    p = wr(ien)
    q = wi(ien)
    na = ien - 1
    IF q = ZERO THEN
      m = ien
      xmat(ien, ien) = ONE
      FOR i = na TO 0 STEP -1
	w = xmat(i, i) - p
	r = xmat(i, ien)
	FOR j = m TO na
	  r = r + xmat(i, j) * xmat(j, ien)
	NEXT j
	IF wi(i) < ZERO THEN
	  z = w
	  s = r
	ELSE
	  m = i
	  IF wi(i) = ZERO THEN
	    IF w <> ZERO THEN
	      temp = w
	    ELSE
	      temp = XMACHEPS * xnorm
	    END IF
	    xmat(i, ien) = -r / temp
	  ELSE
	    'Solve the linear system:
	    '| w   x |  | xmat(i,ien)   |   | -r |
	    '|       |  |               | = |    |
	    '| y   z |  | xmat(i+1,ien) |   | -s |
	    x = xmat(i, i + 1)
	    y = xmat(i + 1, i)
	    q = (wr(i) - p) ^ 2 + wi(i) ^ 2
	    xmat(i, ien) = (x * s - z * r) / q
	    t = xmat(i, ien)
	    IF ABS(x) > ABS(z) THEN
	      temp = (-r - w * t) / x
	    ELSE
	      temp = (-s - y * t) / z
	    END IF
	    xmat(i + 1, ien) = temp
	  END IF
	END IF 'wi(i) < 0
      NEXT i
    ELSEIF q < ZERO THEN
      m = na
      IF ABS(xmat(ien, na)) > ABS(xmat(na, ien)) THEN
	xmat(na, na) = -(xmat(ien, ien) - p) / xmat(ien, na)
	xmat(na, ien) = -q / xmat(ien, na)
      ELSE
 'CALL Comdiv(-xmat(na,ien),0,xmat(na,na)-p,q,xmat(na,na),xmat(na,ien),code)
 'Note: return code not used here (must be 0).
	ar = -xmat(na, ien): ai = 0#: br = xmat(na, na) - p: bi = q
	GOSUB 900: xmat(na, na) = cr: xmat(na, ien) = ci
      END IF
      xmat(ien, na) = ONE
      xmat(ien, ien) = ZERO
      FOR i = na - 1 TO 0 STEP -1
	w = xmat(i, i) - p
	ra = xmat(i, ien)
	sa = ZERO
	FOR j = m TO na
	  ra = ra + xmat(i, j) * xmat(j, na)
	  sa = sa + xmat(i, j) * xmat(j, ien)
	NEXT j

	IF wi(i) < ZERO THEN
	  z = w
	  r = ra
	  s = sa
	ELSE
	  m = i
	  IF wi(i) = ZERO THEN
	    'CALL Comdiv(-ra,-sa,w,q,xmat(i,na),xmat(i,ien),code)
	    ar = -ra: ai = -sa: br = w: bi = q
	    GOSUB 900: xmat(i, na) = cr: xmat(i, ien) = ci
	  ELSE
	  '  Solve complex linear system:
	  '| w+i*q     x | | xmat(i,na) + i*xmat(i,ien)  |   | -ra+i*sa |
	  '|             | |                             | = |          |
	  '|   y    z+i*q| | xmat(i+1,na)+i*xmat(i+1,en) |   | -r+i*s   |
	    x = xmat(i, i + 1)
	    y = xmat(i + 1, i)
	    vr = (wr(i) - p) ^ 2 + wi(i) ^ 2 - q * q
	    vi = TWO * q * (wr(i) - p)
	    IF vr = ZERO AND vi = ZERO THEN
	      vr = XMACHEPS
	    END IF
'CALL Comdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi,xmat(i,na),xmat(i,ien),code)
	    ar = x * r - z * ra + q * sa: ai = x * s - z * sa - q * ra
	    br = vr: bi = vi: GOSUB 900: xmat(i, na) = cr: xmat(i, ien) = ci
	    IF ABS(x) > ABS(z) + ABS(q) THEN
	    xmat(i + 1, na) = (-ra - w * xmat(i, na) + q * xmat(i, ien)) / x
	    xmat(i + 1, ien) = (-sa - w * xmat(i, ien) - q * xmat(i, na)) / x
	    ELSE
'CALL Comdiv(-r-y*h(i,na),-s-y*xmat(i,ien),z,q,xmat(i+1,na),xmat(i+1,ien),code)
	      ar = -r - y * xmat(i, na): ai = -s - y * xmat(i, ien)
	      br = z: bi = q: GOSUB 900
	      xmat(i + 1, na) = cr: xmat(i + 1, ien) = ci
	    END IF
	  END IF 'wi(i) = 0
	END IF 'wi(i) < 0
      NEXT i
    END IF 'else if q < 0
  NEXT ien

  FOR i = 0 TO n - 1                   'Eigenvectors for the evalues for
    IF i < low OR i > ihigh THEN       'rows < low  and rows  > high
      FOR k = i + 1 TO n - 1
	ev(i, k) = xmat(i, k)
      NEXT k
    END IF
  NEXT i

  j = n - 1
  WHILE (j >= low)
    IF j <= ihigh THEN
      m = j
    ELSE
      j = ihigh
    END IF
    IF wi(j) < ZERO THEN
      l = j - 1
      FOR i = low TO ihigh
	y = ZERO: z = ZERO
	FOR k = low TO m
	  y = y + ev(i, k) * xmat(k, l)
	  z = z + ev(i, k) * xmat(k, j)
	NEXT k
	ev(i, l) = y
	ev(i, j) = z
      NEXT i
    ELSE
      IF wi(j) = ZERO THEN
	FOR i = low TO ihigh
	  z = ZERO
	  FOR k = low TO m
	    z = z + ev(i, k) * xmat(k, j)
	  NEXT k
	  ev(i, j) = z
	NEXT i
      END IF
    END IF
    j = j - 1
  WEND 'while j>=low

  irc = 0
RETURN 'hqrvec

'utility subroutines to print matrices
1400 FOR i = 0 TO n - 1
       FOR j = 0 TO n - 1
	 PRINT #2, USING "###.######"; xmat(i, j);
       NEXT j
       PRINT #2,
     NEXT i
RETURN

1410 FOR i = 0 TO n - 1
       FOR j = 0 TO n - 1
	 PRINT #2, USING "###.######"; ev(i, j);
       NEXT j
       PRINT #2,
     NEXT i
RETURN

'Subroutine hqr2           "compute eigenvalues .........
'                (ivec     "switch for computing evectors
'                 n,       "Dimension of matrix .........
'                 low,     "first nonzero row ...........
'                 ihigh,   "last nonzero row ............
'                 xmat,    "Hessenberg matrix ...........
'                 wr,      "Real parts of eigenvalues ...
'                 wi,      "Imaginary parts of evalues ..
'                 ev       "Matrix of eigenvectors ......
'                 icnt,    "Iteration counter ...........
'                 irc      "return code .................
'                )
'integer ivec, ihigh, icnt(NSIZE), irc
'**********************************************************************
'*                                                                    *
'* hqr2 computes the eigenvalues and (if ivec <> 0) the eigenvectors  *
'* of an  n * n upper Hessenberg matrix.                              *
'*                                                                    *
'* ------------------------------------------------------------------ *
'*                                                                    *
'*   Control parameter:                                               *
'*   -----------------                                                *
'*      ivec     integer;                                             *
'*       = 0     compute eigenvalues only                             *
'*       = 1     compute all eigenvalues and eigenvectors             *
'*                                                                    *
'*   Input parameters:                                                *
'*   ----------------                                                 *
'*      n        integer;  ( n > 0 )                                  *
'*               Dimension of xmat and ev,                            *
'*               length of the real parts vector  wr and of the       *
'*               imaginary parts vector  wi of the eigenvalues.       *
'*      low      integer;                                             *
'*      ihigh    integer;  see balance                                *
'*      xmat     n x n matrix:                                        *
'*               upper Hessenberg matrix as output of Elmhes          *
'*               (destroyed in the process).                          *
'*                                                                    *
'*   Output parameters:                                               *
'*   -----------------                                                *
'*      ev       n x n matrix;  (only if ivec <> 0)                   *
'*               Matrix, which for ivec <> 0 contains the             *
'*               eigenvectors as follows:                             *
'*               For real eigebvalues the corresponding column        *
'*               contains the corresponding eigenvector, while for    *
'*               complex eigenvalues the corresponding column contains*
'*               the real part of the eigenvector with its imaginary  *
'*               part is stored in the subsequent column of ev.       *
'*               The eigenvector for the complex conjugate eigenvector*
'*               is given by the complex conjugate eigenvector.       *
'*      wr       vector of size n;                                    *
'*               Real part of the n eigenvalues.                      *
'*      wi       vector of size n;                                    *
'*               Imaginary parts of the eigenvalues                   *
'*      icnt     Integer vector of size n;                            *
'*               vector of iterations used for each eigenvalue.       *
'*               For a complex conjugate eigenvalue pair the second   *
'*               entry is negative.                                   *
'*                                                                    *
'*   Return code:                                                     *
'*   ------------                                                     *
'*      =   0    all ok                                               *
'*      = 4xx    Iteration maximum exceeded when computing evalue xx  *
'*      =  99    zero  matrix                                         *
'*                                                                    *
'* ------------------------------------------------------------------ *
'*                                                                    *
'*   Subroutine in use:                                               *
'*   -----------------                                                *
'*                                                                    *
'*   1200 hqrvec:    reverse transform for eigenvectors               *
'*                                                                    *
'* ------------------------------------------------------------------ *
'*                                                                    *
'*   Constants used:     XMACHEPS, MAXIT                              *
'*   --------------                                                   *
'*                                                                    *
'**********************************************************************
'Labels: 1510, 1512, 1515, 1530
'integer idebug, ien
'double p, q, r, s, t, w, x, y, z
  
1500 idebug = 0
  IF idebug <> 0 THEN
    PRINT #2, " Matrix ev at begin hqr2:"
    GOSUB 1410
  END IF

  p = ZERO: q = ZERO: r = ZERO
  FOR i = 0 TO n - 1
    IF i < low OR i > ihigh THEN
      wr(i) = xmat(i, i)
      wi(i) = ZERO
      icnt(i) = 0
    END IF
  NEXT i

  ien = ihigh
  t = ZERO

  WHILE (ien >= low)
    iter = 0
    na = ien - 1

    WHILE (1 < 2)
      ll = 999
      FOR l = ien TO low + 1 STEP -1            'search for small
						'subdiagonal element
	IF ABS(xmat(l, l - 1)) <= XMACHEPS * (ABS(xmat(l - 1, l - 1) + ABS(xmat(l, l)))) THEN
	  ll = l                                'save current index
	  GOTO 1510                             'exit l loop
	END IF
      NEXT l
1510  IF ll <> 999 THEN
	l = ll
      ELSE
	l = 0                                   'restore l
      END IF

      x = xmat(ien, ien)
      IF l = ien THEN                           'found one evalue
	wr(ien) = x + t
	xmat(ien, ien) = x + t
	wi(ien) = ZERO
	icnt(ien) = iter
	ien = ien - 1
	GOTO 1515                               'exit from loop while(1<2)
      END IF

      y = xmat(na, na)
      w = xmat(ien, na) * xmat(na, ien)

      IF l = na THEN                           'found two evalues
	p = (y - x) * .5#
	q = p * p + w
	z = SQR(ABS(q))
	x = x + t
	xmat(ien, ien) = x + t
	xmat(na, na) = y + t
	icnt(ien) = -iter
	icnt(na) = iter
	IF q >= ZERO THEN                       'real eigenvalues
	  IF p < ZERO THEN
	    z = p - z
	  ELSE
	    z = p + z
	  END IF
	  wr(na) = x + z
	  wr(ien) = x - w / z
	  s = w - w / z
	  wi(na) = ZERO
	  wi(ien) = ZERO
	  x = xmat(ien, na)
	  r = SQR(x * x + z * z)

	  IF ivec <> 0 THEN
	    p = x / r
	    q = z / r
	    FOR j = na TO n - 1
	      z = xmat(na, j)
	      xmat(na, j) = q * z + p * xmat(ien, j)
	      xmat(ien, j) = q * xmat(ien, j) - p * z
	    NEXT j

	    FOR i = 0 TO ien
	      z = xmat(i, na)
	      xmat(i, na) = q * z + p * xmat(i, ien)
	      xmat(i, ien) = q * xmat(i, ien) - p * z
	    NEXT i

	    FOR i = low TO ihigh
	      z = ev(i, na)
	      ev(i, na) = q * z + p * ev(i, ien)
	      ev(i, ien) = q * ev(i, ien) - p * z
	    NEXT i
	  END IF 'if ivec <> ZERO
	ELSE                                  'pair of complex
	  wr(na) = x + p
	  wr(ien) = x + p
	  wi(na) = z
	  wi(ien) = -z
	END IF 'if q>=ZERO

	ien = ien - 2
	GOTO 1515                             'exit while(1<2)
      END IF 'if l = na

      IF iter >= MAXIT THEN
	icnt(ien) = MAXIT + 1
	irc = ien
	PRINT " stop at iter >= MAXIT."
	RETURN                               'MAXIT Iterations: return
      END IF

      IF iter <> 0 AND (iter MOD 10) = 0 THEN
	t = t + x
	FOR i = low TO ien
	  xmat(i, i) = xmat(i, i) - x
	NEXT i
	s = ABS(xmat(ien, na)) + ABS(xmat(na, ien - 2))
	x = .75# * s: y = x
	w = -.4375# * s * s
      END IF

      iter = iter + 1

      FOR m = ien - 2 TO l STEP -1
	z = xmat(m, m)
	r = x - z
	s = y - z
	p = (r * s - w) / xmat(m + 1, m) + xmat(m, m + 1)
	q = xmat(m + 1, m + 1) - z - r - s
	r = xmat(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 1512
	IF ABS(xmat(m, m - 1)) * (ABS(q) + ABS(r)) <= XMACHEPS * ABS(p) * (ABS(xmat(m - 1, m - 1)) + ABS(z) + ABS(h(m + 1, m + 1))) THEN
	  GOTO 1512                'exit m loop
	END IF
      NEXT m

1512  FOR i = m + 2 TO ien
	xmat(i, i - 2) = ZERO
      NEXT i
      FOR i = m + 3 TO ien
	xmat(i, i - 3) = ZERO
      NEXT i

      FOR k = m TO na
	IF k <> m THEN           'double  QR step, for rows l to ien
				 'and columns m to ien
	  p = xmat(k, k - 1)
	  q = xmat(k + 1, k - 1)
	  IF k <> na THEN
	    r = xmat(k + 2, k - 1)
	  ELSE
	    r = ZERO
	  END IF
	  x = ABS(p) + ABS(q) + ABS(r)
	  IF x = ZERO THEN GOTO 1530               'go to next k
	  p = p / x
	  q = q / x
	  r = r / x
	END IF
	s = SQR(p * p + q * q + r * r)
	IF p < ZERO THEN s = -s

	IF k <> m THEN
	  xmat(k, k - 1) = -s * x
	ELSEIF l <> m THEN
	  xmat(k, k - 1) = -xmat(k, k - 1)
	END IF
	p = p + s
	x = p / s
	y = q / s
	z = r / s
	q = q / p
	r = r / p

	FOR j = k TO n - 1                        'modify rows
	  p = xmat(k, j) + q * xmat(k + 1, j)
	  IF k <> na THEN
	    p = p + r * xmat(k + 2, j)
	    xmat(k + 2, j) = xmat(k + 2, j) - p * z
	  END IF
	  xmat(k + 1, j) = xmat(k + 1, j) - p * y
	  xmat(k, j) = xmat(k, j) - p * x
	NEXT j

	IF k + 3 < ien THEN
	  j = k + 3
	ELSE
	  j = ien
	END IF
	FOR i = 0 TO j                            'modify columns
	  p = x * xmat(i, k) + y * xmat(i, k + 1)
	  IF k <> na THEN
	    p = p + z * xmat(i, k + 2)
	    xmat(i, k + 2) = xmat(i, k + 2) - p * r
	  END IF
	  xmat(i, k + 1) = xmat(i, k + 1) - p * q
	  xmat(i, k) = xmat(i, k) - p
	NEXT i

	IF ivec <> 0 THEN        'if eigenvectors are needed
	  FOR i = low TO ihigh
	    p = x * ev(i, k) + y * ev(i, k + 1)
	    IF k <> na THEN
	      p = p + z * ev(i, k + 2)
	      ev(i, k + 2) = ev(i, k + 2) - p * r
	    END IF
	    ev(i, k + 1) = ev(i, k + 1) - p * q
	    ev(i, k) = ev(i, k) - p
	  NEXT i
	END IF
1530  NEXT k

     WEND 'while(1<2)

1515 WEND 'while ien >= low                     'All evalues found

  IF ivec <> 0 THEN                             'transform evectors back
    IF idebug <> 0 THEN
      PRINT #2, " Matrix ev before hqrvec:"
      GOSUB 1410
    END IF
    'CALL hqrvec(n, low, ihigh, xmat, wr, wi, ev, irc)
    GOSUB 1200
    IF idebug <> 0 THEN
      PRINT #2, " irc="; irc
      PRINT #2, " Matrix ev after hqrvec:"
      GOSUB 1410
    END IF
  ELSE
    irc = 0
  END IF

RETURN

'Subroutine norm_1 (n,      "Dimension of matrix ...........
'                   ev,     "Matrix with eigenvectors ......
'                   wi,     "Imaginary parts of evalues ....
'                   irc     "return code ...................
'                  )
'integer irc
'**********************************************************************
'*                                                                    *
'*  norm_1 normalizes the one norm of the column vectors in ev.       *
'*  (special attention to complex vectors in ev is given)             *
'*                                                                    *
'* ------------------------------------------------------------------ *
'*                                                                    *
'*   Input parameters:                                                *
'*   ----------------                                                 *
'*      n        integer; ( n > 0 )                                   *
'*               Dimension of matrix ev                               *
'*      ev       n x n matrix;                                        *
'*               Matrix of (not normalized) eigenvectors              *
'*      wi       vector of size n:                                    *
'*               Imaginary parts of the eigenvalues                   *
'*                                                                    *
'*   Output parameter:                                                *
'*   ----------------                                                 *
'*      ev       n x n matrix;                                        *
'*               Matrix with normalized eigenvectors                  *
'*               (in columns).                                        *
'*   Return code:                                                     *
'*   -----------                                                      *
'*      = 0      all ok                                               *
'*      = 1      n < 1                                                *
'*                                                                    *
'* ------------------------------------------------------------------ *
'*                                                                    *
'*   function or subroutine used:                                     *
'*   ---------------------------                                      *
'*      double   comabs():  complex absolute value                    *
'*               comdiv():  complex division                          *
'*                                                                    *
'**********************************************************************
1600 'double xmaxi, tr, ti

  IF n < 1 THEN
    irc = 1
    RETURN
  END IF

  j = 0
  WHILE (j < n)
    IF wi(j) = ZERO THEN
      xmaxi = ev(0, j)
      FOR i = 1 TO n - 1
	IF ABS(ev(i, j)) > ABS(xmaxi) THEN xmaxi = ev(i, j)
      NEXT i

      IF xmaxi <> ZERO THEN
	xmaxi = ONE / xmaxi
	FOR i = 0 TO n - 1
	  ev(i, j) = ev(i, j) * xmaxi
	NEXT i
      END IF
    ELSE
      tr = ev(0, j)
      ti = ev(0, j + 1)
      FOR i = 1 TO n - 1
	ar = ev(i, j): ai = ev(i, j + 1)
	GOSUB 1000: temp = Comabs     'temp=Comabs(ar,ai)
	ar = tr: ai = ti: GOSUB 1000
	IF temp > Comabs THEN
	  tr = ev(i, j)
	  ti = ev(i, j + 1)
	END IF
      NEXT i

      IF tr <> ZERO OR ti <> ZERO THEN
	FOR i = 0 TO n - 1
'CALL Comdiv(ev(i, j), ev(i, j + 1), tr, ti, ev(i, j), ev(i, j + 1), irc)
	ar = ev(i, j): ai = ev(i, j + 1): br = tr: bi = ti
	GOSUB 900: ev(i, j) = cr: ev(i, j + 1) = ci
	NEXT i
      END IF
      j = j + 1                                        'raise j by two
    END IF
    j = j + 1
  WEND 'j loop
  irc = 0
RETURN

'Subroutine eigen (        "Compute all evalues/evectors of a matrix
'           ivec,          "switch for computing evectors ...
'           ievnorm        "<> 0 for normalized eigenvectors
'           n,             "size of matrix ..................
'           xmat,          "input matrix ....................
'           ev,            "Eigenvectors ....................
'           wr,            "real parts of eigenvalues .......
'           wi,            "imaginary parts of eigenvalues ..
'           icnt,          "Iteration counter ...............
'           irc            "return code .....................
'          )
'**********************************************************************
'*                                                                    *
'* The subroutine eigen  determines all eigenvalues and (if desired)  *
'* all eigenvectors of a real square  n * n  matrix via the QR method *
'* in the version of Martin, Parlett, Peters, Reinsch and Wilkinson.  *
'*                                                                    *
'* ------------------------------------------------------------------ *
'*                                                                    *
'*   Litterature:                                                     *
'*   -----------                                                      *
'*      1) Peters, Wilkinson: Eigenvectors of real and complex        *
'*         matrices by LR and QR triangularisations,                  *
'*         Num. Math. 16, p.184-204, (1970): [PETE70]: contribution   *
'*         II/15, p. 372 - 395 in [WILK71].                           *
'*      2) Martin, Wilkinson: Similarity reductions of a general      *
'*         matrix to Hessenberg form, Num. Math. 12, p. 349-368,(1968)*
'*         [MART 68]: contribution II,13, p. 339 - 358 in [WILK71].   *
'*      3) Parlett, Reinsch: Balancing a matrix for calculations of   *
'*         eigenvalues and eigenvectors, Num. Math. 13, p. 293-304,   *
'*         (1969): [PARL69]: contribution II/11, p.315 - 326 in       *
'*         [WILK71].                                                  *
'*                                                                    *
'* ------------------------------------------------------------------ *
'*                                                                    *
'*   Control parameters:                                              *
'*   ------------------                                               *
'*      ivec     integer;                                             *
'*               call for eigen :                                     *
'*      = 0      compute eigenvalues only                             *
'*      = 1      compute all eigenvalues and eigenvectors             *
'*      ortho    boolean flag that shows if transformation of mat to  *
'*               Hessenberg form shall be done orthogonally by        *
'*               "orthes" (flag=1) or elementarily by "elmhes"        *
'*               (flag=0). The Householder matrices used in           *
'*               orthogonal transformation have the advantage of      *
'*               preserving the symmetry of input matrices.           *
'*               (only elmhes is implemented here).                   *
'*      ievnorm  integer flag that shows if Eigenvectors shall be     *
'*               normalized (flag=1) or not (flag=0).                 *
'*                                                                    *
'*   Input parameters:                                                *
'*   ----------------                                                 *
'*      n        integer: ( n > 0 )                                   *
'*               size of matrix, number of eigenvalues                *
'*      xmat     n x n matrix:                                        *
'*               input matrix                                         *
'*                                                                    *
'*   Output parameters:                                               *
'*   -----------------                                                *
'*      ev       n x n matrix;     (only if ivec <> 0)                *
'*               matrix, if ivec<>0 that holds the eigenvectors thus: *
'*               If the jth eigenvalue of the matrix is real then the *
'*               jth column is the corresponding real eigenvector:    *
'*               if the jth eigenvalue is complex then the jth column *
'*               of ev contains the real part of the eigenvector while*
'*               its imaginary part is in column j+1.                 *
'*               (the j+1st eigenvector is the complex conjugate      *
'*               vector.)                                             *
'*      wr       vector of size n;                                    *
'*               Real parts of the eigenvalues.                       *
'*      wi       vector of size n;                                    *
'*               Imaginary parts of the eigenvalues                   *
'*      icnt     Integer vector of size n:                            *
'*               vector containing the number of iterations for each  *
'*               eigenvalue. (for a complex conjugate pair the second *
'*               entry is negative).                                  *
'*                                                                    *
'*   Return code:                                                     *
'*   ------------                                                     *
'*      =   0    all ok                                               *
'*      =   1    n < 1 or other invalid input parameter               *
'*      =   2    insufficient memory                                  *
'*      = 10x    error x from balance()                               *
'*      = 20x    error x from elmh()                                  *
'*      = 30x    error x from elmtrans()   (for vec = Trie only)      *
'*      = 4xx    error xx from hqr2()                                 *
'*      = 50x    error x from balback()    (for vec = True only)      *
'*      = 60x    error x from norm_1()     (for vec = True only)      *
'*                                                                    *
'* ------------------------------------------------------------------ *
'*                                                                    *
'*  Subroutines used:                                                 *
'*  ----------------                                                  *
'*                                                                    *
'*   500 balance() : Balancing of an  n x n  matrix                   *
'*   600 balback() : Reverse balancing to obtain eigenvectors         *
'*   700 elmhes()  : Transformation to upper Hessenberg form          *
'*   800 elmtrans(): Intialize eigenvectors                           *
'*   900 comdiv    : Divide two complex numbers                       *
'*  1000 comabs    : This function returns the absolute value         *
'*                   of a complex number                              *
'*  1200 hqrvec    : Compute eigenvectors (called by 1500)            *
'*  1500 hqr2()    : Compute eigenvalues/eigenvectors                 *
'*  1600 norm_1()  : Normalize eigenvectors                           *
'*                                                                    *
'**********************************************************************
'integer ihigh
 
2000 idebug = 0
  
  IF n < 1 THEN                            'case n < 1 .............
    irc = 1
    RETURN
  END IF

  FOR i = 0 TO n - 1
    icnt(i) = 0
  NEXT i

  IF n = 1 THEN                               'case n = 1 ...........
    ev(0, 0) = ONE
    wr(0) = xmat(0, 0)
    wi(0) = ZERO
    irc = 0
    RETURN
  END IF
					      'balance xmat for nearly
  'CALL Balance(n, xmat, scale, low, ihigh)    equal row and column
  GOSUB 500                                   'one norms
  
  IF idebug <> 0 THEN
    PRINT #2, " Matrix xmat after Balance:"
    GOSUB 1400
  END IF

  'CALL Elmhes(n, low, ihigh, xmat, iperm)     reduce xmat to upper
  GOSUB 700                                   'Hessenberg form

  IF idebug <> 0 THEN
    PRINT #2, " Matrix xmat after Elmhes:"
    GOSUB 1400
  END IF

  IF ivec <> 0 THEN                           'initialize eigenvectors
   'CALL elmtrans(n, low, ihigh, xmat, iperm, ev)
    GOSUB 800
  END IF

  IF idebug <> 0 THEN
    PRINT #2, " Matrix xmat after Elmtrans:"
    GOSUB 1400
    PRINT #2, " Vector iperm after Elmtrans:"
    FOR i = 0 TO n - 1
      PRINT #2, USING "####"; iperm(i);
    NEXT i
    PRINT #2,
    PRINT #2, " Matrix ev after Elmtrans:"
    GOSUB 1410
  END IF

  'CALL hqr2 (ivec, n, low, ihigh, xmat,    'execute Francis QR algorithm
  '           wr, wi, ev, icnt, irc)         to obtain the eigenvalues and
  GOSUB 1500                                'the eigenvectors of input matrix

  IF irc <> 0 THEN
    irc = 2
    RETURN
  END IF

  IF ivec <> 0 THEN

    'CALL balback(n, low, ihigh, scal, ev)      'reverse balancing if
    GOSUB 600                                   'eigenvaectors are to
						'be determined
    IF ievnorm <> 0 THEN
      'CALL norm_1                              'normalize eigenvectors
      GOSUB 1600
    END IF

    IF irc <> 0 THEN
      irc = 3
      RETURN
    END IF
  END IF
  irc = 0
RETURN

'end of file thqr.bas