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