c11-6 (Numerical Recipes in C), страница 2

There are twopossible ways of terminating the iteration for an eigenvalue. First, if an,n−1 becomes“negligible,” then ann is an eigenvalue. We can then delete the nth row and columnSample 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).where H is upper Hessenberg. Thus490Chapter 11.Eigensystems|am,m−1 |(|q| + |r|) |p|(|am+1,m+1 | + |amm | + |am−1,m−1 |)(11.6.26)Very rarely, the procedure described so far will fail to converge.

On suchmatrices, experience shows that if one double step is performed with any shiftsthat are of order the norm of the matrix, convergence is subsequently very rapid.Accordingly, if ten iterations occur without determining an eigenvalue, the usualshifts are replaced for the next iteration by shifts defined byks + ks+1 = 1.5 × (|an,n−1| + |an−1,n−2|)ks ks+1 = (|an,n−1| + |an−1,n−2|)2(11.6.27)The factor 1.5 was arbitrarily chosen to lessen the likelihood of an “unfortunate”choice of shifts. This strategy is repeated after 20 unsuccessful iterations. After 30unsuccessful iterations, the routine reports failure.The operation count for the QR algorithm described here is ∼ 5k 2 per iteration,where k is the current size of the matrix.

The typical average number of iterations pereigenvalue is ∼ 1.8, so the total operation count for all the eigenvalues is ∼ 3n3 . Thisestimate neglects any possible efficiency due to splitting or sparseness of the matrix.The following routine hqr is based algorithmically on the above description,in turn following the implementations in [1,2].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).of the matrix and look for the next eigenvalue. Alternatively, an−1,n−2 may becomenegligible. In this case the eigenvalues of the 2 × 2 matrix in the lower right-handcorner may be taken to be eigenvalues.

We delete the nth and (n − 1)st rows andcolumns of the matrix and continue.The test for convergence to an eigenvalue is combined with a test for negligiblesubdiagonal elements that allows splitting of the matrix into submatrices. We findthe largest i such that ai,i−1 is negligible. If i = n, we have found a singleeigenvalue. If i = n − 1, we have found two eigenvalues. Otherwise we continuethe iteration on the submatrix in rows i to n (i being set to unity if there is nosmall subdiagonal element).After determining i, the submatrix in rows i to n is examined to see if theproduct of any two consecutive subdiagonal elements is small enough that wecan work with an even smaller submatrix, starting say in row m. We start withm = n − 2 and decrement it down to i + 1, computing p, q, and r according toequations (11.6.23) with 1 replaced by m and 2 by m + 1.

If these were indeed theelements of the special “first” Householder matrix in a double QR step, then applyingthe Householder matrix would lead to nonzero elements in positions (m + 1, m − 1),(m + 2, m − 1), and (m + 2, m). We require that the first two of these elements besmall compared with the local diagonal elements am−1,m−1 , amm and am+1,m+1 .A satisfactory approximate criterion is11.6 The QR Algorithm for Real Hessenberg Matrices491#include <math.h>#include "nrutil.h"anorm=0.0;Compute matrix norm for possible use in lofor (i=1;i<=n;i++)cating single small subdiagonal element.for (j=IMAX(i-1,1);j<=n;j++)anorm += fabs(a[i][j]);nn=n;t=0.0;Gets changed only by an exceptional shift.while (nn >= 1) {Begin search for next eigenvalue.its=0;do {for (l=nn;l>=2;l--) {Begin iteration: look for single small subdis=fabs(a[l-1][l-1])+fabs(a[l][l]);agonal element.if (s == 0.0) s=anorm;if ((float)(fabs(a[l][l-1]) + s) == s) break;}x=a[nn][nn];if (l == nn) {One root found.wr[nn]=x+t;wi[nn--]=0.0;} else {y=a[nn-1][nn-1];w=a[nn][nn-1]*a[nn-1][nn];if (l == (nn-1)) {Two roots found...p=0.5*(y-x);q=p*p+w;z=sqrt(fabs(q));x += t;if (q >= 0.0) {...a real pair.z=p+SIGN(z,p);wr[nn-1]=wr[nn]=x+z;if (z) wr[nn]=x-w/z;wi[nn-1]=wi[nn]=0.0;} else {...a complex pair.wr[nn-1]=wr[nn]=x+p;wi[nn-1]= -(wi[nn]=z);}nn -= 2;} else {No roots found.

Continue iteration.if (its == 30) nrerror("Too many iterations in hqr");if (its == 10 || its == 20) {Form exceptional shift.t += x;for (i=1;i<=nn;i++) a[i][i] -= x;s=fabs(a[nn][nn-1])+fabs(a[nn-1][nn-2]);y=x=0.75*s;w = -0.4375*s*s;}++its;for (m=(nn-2);m>=l;m--) {Form shift and then look forz=a[m][m];2 consecutive small subr=x-z;diagonal elements.s=y-z;p=(r*s-w)/a[m+1][m]+a[m][m+1];Equation (11.6.23).q=a[m+1][m+1]-z-r-s;r=a[m+2][m+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).void hqr(float **a, int n, float wr[], float wi[])Finds all eigenvalues of an upper Hessenberg matrix a[1..n][1..n]. On input a can beexactly as output from elmhes §11.5; on output it is destroyed.

The real and imaginary partsof the eigenvalues are returned in wr[1..n] and wi[1..n], respectively.{int nn,m,l,k,j,its,i,mmin;float z,y,x,w,v,u,t,s,r,q,p,anorm;492Chapter 11.Eigensystems}for (i=m+2;i<=nn;i++) {a[i][i-2]=0.0;if (i != (m+2)) a[i][i-3]=0.0;}for (k=m;k<=nn-1;k++) {Double QR step on rows l to nn and columns m to nn.if (k != m) {p=a[k][k-1];Begin setup of Householderq=a[k+1][k-1];vector.r=0.0;if (k != (nn-1)) r=a[k+2][k-1];if ((x=fabs(p)+fabs(q)+fabs(r)) != 0.0) {p /= x;Scale to prevent overflow orq /= x;underflow.r /= x;}}if ((s=SIGN(sqrt(p*p+q*q+r*r),p)) != 0.0) {if (k == m) {if (l != m)a[k][k-1] = -a[k][k-1];} elsea[k][k-1] = -s*x;p += s;Equations (11.6.24).x=p/s;y=q/s;z=r/s;q /= p;r /= p;for (j=k;j<=nn;j++) {Row modification.p=a[k][j]+q*a[k+1][j];if (k != (nn-1)) {p += r*a[k+2][j];a[k+2][j] -= p*z;}a[k+1][j] -= p*y;a[k][j] -= p*x;}mmin = nn<k+3 ? nn : k+3;for (i=l;i<=mmin;i++) {Column modification.p=x*a[i][k]+y*a[i][k+1];if (k != (nn-1)) {p += z*a[i][k+2];a[i][k+2] -= p*r;}a[i][k+1] -= p*q;a[i][k] -= p;}}}}}} while (l < nn-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).s=fabs(p)+fabs(q)+fabs(r);Scale to prevent overflow orp /= s;underflow.q /= s;r /= s;if (m == l) break;u=fabs(a[m][m-1])*(fabs(q)+fabs(r));v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1]));if ((float)(u+v) == v) break;Equation (11.6.26).49311.7 Eigenvalues or Eigenvectors by Inverse IterationCITED 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]Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns HopkinsUniversity Press), §7.5.11.7 Improving Eigenvalues and/or FindingEigenvectors by Inverse IterationThe basic idea behind inverse iteration is quite simple.

Let y be the solutionof the linear system(A − τ 1) · y = b(11.7.1)where b is a random vector and τ is close to some eigenvalue λ of A. Then thesolution y will be close to the eigenvector corresponding to λ. The procedure canbe iterated: Replace b by y and solve for a new y, which will be even closer tothe true eigenvector.We can see why this works by expanding both y and b as linear combinationsof the eigenvectors xj of A:y=Xα j xjb=Xjβj xj(11.7.2)jThen (11.7.1) givesXαj (λj − τ )xj =Xjβj xj(11.7.3)jso thatαj =andy=βjλj − τX βj xjλj − τ(11.7.4)(11.7.5)jIf τ is close to λn , say, then provided βn is not accidentally too small, y will beapproximately xn , up to a normalization.

Moreover, each iteration of this proceduregives another power of λj − τ in the denominator of (11.7.5). Thus the convergenceis rapid for well-separated eigenvalues.Suppose at the kth stage of iteration we are solving the equation(A − τk 1) · y = bk(11.7.6)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).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].