c11-4 (Numerical Recipes in C)

11.4 Hermitian Matrices481}} while (m != l);}}CITED REFERENCES AND FURTHER READING:Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), pp. 331–335. [1]Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag). [2]Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide, 2nd ed., vol. 6 ofLecture Notes in Computer Science (New York: Springer-Verlag). [3]Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),§6.6.6. [4]11.4 Hermitian MatricesThe complex analog of a real, symmetric matrix is a Hermitian matrix,satisfying equation (11.0.4).

Jacobi transformations can be used to find eigenvaluesand eigenvectors, as also can Householder reduction to tridiagonal form followed byQL iteration. Complex versions of the previous routines jacobi, tred2, and tqliare quite analogous to their real counterparts. For working routines, consult [1,2] .An alternative, using the routines in this book, is to convert the Hermitianproblem to a real, symmetric one: If C = A + iB is a Hermitian matrix, then then × n complex eigenvalue problem(A + iB) · (u + iv) = λ(u + iv)(11.4.1)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).for (i=m-1;i>=l;i--) {A plane rotation as in the origif=s*e[i];nal QL, followed by Givensb=c*e[i];rotations to restore tridiage[i+1]=(r=pythag(f,g));onal form.if (r == 0.0) {Recover from underflow.d[i+1] -= p;e[m]=0.0;break;}s=f/r;c=g/r;g=d[i+1]-p;r=(d[i]-g)*s+2.0*c*b;d[i+1]=g+(p=s*r);g=c*r-b;/* Next loop can be omitted if eigenvectors not wanted*/for (k=1;k<=n;k++) {Form eigenvectors.f=z[k][i+1];z[k][i+1]=s*z[k][i]+c*f;z[k][i]=c*z[k][i]-s*f;}}if (r == 0.0 && i >= l) continue;d[l] -= p;e[l]=g;e[m]=0.0;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)u.