'*===================================================================== '* 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