c11-5 (Numerical Recipes in C)

482Chapter 11.Eigensystemsis equivalent to the 2n × 2n real problem A −Buu·=λB Avv(11.4.2)is also an eigenvector, as you can verify by writing out the two matrix equations implied by (11.4.2). Thus if λ1 , λ2 , . . . , λn are the eigenvalues of C,then the 2n eigenvalues of the augmented problem (11.4.2) are λ1 , λ1 , λ2 , λ2 , . . .

,λn , λn ; each, in other words, is repeated twice. The eigenvectors are pairs of the formu + iv and i(u + iv); that is, they are the same up to an inessential phase. Thus wesolve the augmented problem (11.4.2), and choose one eigenvalue and eigenvectorfrom each pair. These give the eigenvalues and eigenvectors of the original matrix C.Working with the augmented matrix requires a factor of 2 more storage than theoriginal complex matrix.

In principle, a complex algorithm is also a factor of 2 moreefficient in computer time than is the solution of the augmented problem.CITED REFERENCES AND FURTHER READING:Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag). [1]Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide, 2nd ed., vol. 6 ofLecture Notes in Computer Science (New York: Springer-Verlag). [2]11.5 Reduction of a General Matrix toHessenberg FormThe algorithms for symmetric matrices, given in the preceding sections, arehighly satisfactory in practice. By contrast, it is impossible to design equallysatisfactory algorithms for the nonsymmetric case.

There are two reasons for this.First, the eigenvalues of a nonsymmetric matrix can be very sensitive to small changesin the matrix elements. Second, the matrix itself can be defective, so that there isno complete set of eigenvectors. We emphasize that these difficulties are intrinsicproperties of certain nonsymmetric matrices, and no numerical procedure can “cure”them. The best we can hope for are procedures that don’t exacerbate such problems.The presence of rounding error can only make the situation worse. With finiteprecision arithmetic, one cannot even design a foolproof algorithm to determinewhether a given matrix is defective or not. Thus current algorithms generally try tofind a complete set of eigenvectors, and rely on the user to inspect the results. If anyeigenvectors are almost parallel, the matrix is probably defective.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.

Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Note that the 2n × 2n matrix in (11.4.2) is symmetric: AT = A and BT = −Bif C is Hermitian.Corresponding to a given eigenvalue λ, the vector−v(11.4.3)u11.5 Reduction of a General Matrix to Hessenberg Form483Apart from referring you to the literature, and to the collected routines in [1,2], weare going to sidestep the problem of eigenvectors, giving algorithms for eigenvaluesonly.

If you require just a few eigenvectors, you can read §11.7 and consider findingthem by inverse iteration. We consider the problem of finding all eigenvectors of anonsymmetric matrix as lying beyond the scope of this book.The sensitivity of eigenvalues to rounding errors during the execution ofsome algorithms can be reduced by the procedure of balancing. The errors inthe eigensystem found by a numerical procedure are generally proportional to theEuclidean norm of the matrix, that is, to the square root of the sum of the squaresof the elements.

The idea of balancing is to use similarity transformations tomake corresponding rows and columns of the matrix have comparable norms, thusreducing the overall norm of the matrix while leaving the eigenvalues unchanged.A symmetric matrix is already balanced.Balancing is a procedure with of order N 2 operations.

Thus, the time takenby the procedure balanc, given below, should never be more than a few percentof the total time required to find the eigenvalues. It is therefore recommended thatyou always balance nonsymmetric matrices. It never hurts, and it can substantiallyimprove the accuracy of the eigenvalues computed for a badly balanced matrix.The actual algorithm used is due to Osborne, as discussed in [1].

It consists of asequence of similarity transformations by diagonal matrices D. To avoid introducingrounding errors during the balancing process, the elements of D are restricted to beexact powers of the radix base employed for floating-point arithmetic (i.e., 2 formost machines, but 16 for IBM mainframe architectures).

The output is a matrixthat is balanced in the norm given by summing the absolute magnitudes of thematrix elements. This is more efficient than using the Euclidean norm, and equallyeffective: A large reduction in one norm implies a large reduction in the other.Note that if the off-diagonal elements of any row or column of a matrix areall zero, then the diagonal element is an eigenvalue.

If the eigenvalue happens tobe ill-conditioned (sensitive to small changes in the matrix elements), it will haverelatively large errors when determined by the routine hqr (§11.6). Had we merelyinspected the matrix beforehand, we could have determined the isolated eigenvalueexactly and then deleted the corresponding row and column from the matrix.

Youshould consider whether such a pre-inspection might be useful in your application.(For symmetric matrices, the routines we gave will determine isolated eigenvaluesaccurately in all cases.)The routine balanc does not keep track of the accumulated similarity transformation of the original matrix, since we will only be concerned with findingeigenvalues of nonsymmetric matrices, not eigenvectors. Consult [1-3] if you wantto keep track of the transformation.#include <math.h>#define RADIX 2.0void balanc(float **a, int n)Given a matrix a[1..n][1..n], this routine replaces it by a balanced matrix with identicaleigenvalues. A symmetric matrix is already balanced and is unaffected by this procedure.

Theparameter RADIX should be the machine’s floating-point radix.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.

To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Balancing484Chapter 11.Eigensystems{int last,j,i;float s,r,g,f,c,sqrdx;}Reduction to Hessenberg FormThe strategy for finding the eigensystem of a general matrix parallels that of thesymmetric case.

First we reduce the matrix to a simpler form, and then we performan iterative procedure on the simplified matrix. The simpler structure we use here iscalled Hessenberg form. An upper Hessenberg matrix has zeros everywhere belowthe diagonal except for the first subdiagonal row. For example, in the 6 × 6 case,the nonzero elements are:× × × × × ×× × × × × ×× × × × ×× × × ×× × ×× ×By now you should be able to tell at a glance that such a structure can beachieved by a sequence of Householder transformations, each one zeroing theSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.

Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).sqrdx=RADIX*RADIX;last=0;while (last == 0) {last=1;for (i=1;i<=n;i++) {Calculate row and column norms.r=c=0.0;for (j=1;j<=n;j++)if (j != i) {c += fabs(a[j][i]);r += fabs(a[i][j]);}if (c && r) {If both are nonzero,g=r/RADIX;f=1.0;s=c+r;while (c<g) {find the integer power of the machine radix thatf *= RADIX;comes closest to balancing the matrix.c *= sqrdx;}g=r*RADIX;while (c>g) {f /= RADIX;c /= sqrdx;}if ((c+r)/f < 0.95*s) {last=0;g=1.0/f;for (j=1;j<=n;j++) a[i][j] *= g;Apply similarity transformafor (j=1;j<=n;j++) a[j][i] *= f;tion.}}}}11.5 Reduction of a General Matrix to Hessenberg Form485ni,r+1 ≡airar+1,rSubtract ni,r+1 times row r + 1 from row i.

To make the elimination asimilarity transformation, also add ni,r+1 times column i to column r + 1.A total of N − 2 such stages are required.When the magnitudes of the matrix elements vary over many orders, you shouldtry to rearrange the matrix so that the largest elements are in the top left-hand corner.This reduces the roundoff error, since the reduction proceeds from left to right.Since we are concerned only with eigenvalues, the routine elmhes does notkeep track of the accumulated similarity transformation. The operation count isabout 5N 3 /6 for large N .#include <math.h>#define SWAP(g,h) {y=(g);(g)=(h);(h)=y;}void elmhes(float **a, int n)Reduction to Hessenberg form by the elimination method. The real, nonsymmetric matrixa[1..n][1..n] is replaced by an upper Hessenberg matrix with identical eigenvalues.