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

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).Now choose the vector x for the second Householder matrix to be the bottomn − 2 elements of the second column, and from it construct1 00···00 10···0(11.2.9)P2 ≡  0 0. . .. ..(n−2)P20 011.2 Reduction of a Symmetric Matrix to Tridiagonal Form473Variables are thus computed in the following order: σ, u, H, p, K, q, A0 .

At anystage m, A is tridiagonal in its last m − 1 rows and columns.If the eigenvectors of the final tridiagonal matrix are found (for example, by theroutine in the next section), then the eigenvectors of A can be obtained by applyingthe accumulated transformation(11.2.17)to those eigenvectors. We therefore form Q by recursion after all the P’s havebeen determined:Qn−2 = Pn−2Qj = Pj · Qj+1 ,j = n − 3, . . . , 1(11.2.18)Q = Q1Input for the routine below is the real, symmetric matrix a[1..n][1..n]. Onoutput, a contains the elements of the orthogonal matrix q. The vector d[1..n] isset to the diagonal elements of the tridiagonal matrix A0 , while the vector e[1..n]is set to the off-diagonal elements in its components 2 through n, with e[1]=0.Note that since a is overwritten, you should copy it before calling the routine, if itis required for subsequent computations.No extra storage arrays are needed for the intermediate results.

At stage m, thevectors p and q are nonzero only in elements 1, . . . , i (recall that i = n − m + 1),while u is nonzero only in elements 1, . . . , i − 1. The elements of the vector e arebeing determined in the order n, n − 1, . . . , so we can store p in the elements of enot already determined. The vector q can overwrite p once p is no longer needed.We store u in the ith row of a and u/H in the ith column of a. Once the reductionis complete, we compute the matrices Qj using the quantities u and u/H that havebeen stored in a. Since Qj is an identity matrix in the last n − j + 1 rows andcolumns, we only need compute its elements up to row and column n − j.

Thesecan overwrite the u’s and u/H’s in the corresponding rows and columns of a, whichare no longer required for subsequent Q’s.The routine tred2, given below, includes one further refinement. If the quantityσ is zero or “small” at any stage, one can skip the corresponding transformation.A simple criterion, such asσ<smallest positive number representable on machinemachine precisionwould be fine most of the time. A more careful criterion is actually used. Definethe quantity=i−1X|aik |(11.2.19)k=1If = 0 to machine precision, we skip the transformation.

Otherwise we redefineaikbecomesaik /(11.2.20)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).Q = P1 · P2 · · · Pn−2474Chapter 11.Eigensystems#include <math.h>void tred2(float **a, int n, float d[], float e[])Householder reduction of a real, symmetric matrix a[1..n][1..n]. On output, a is replacedby the orthogonal matrix Q effecting the transformation. d[1..n] returns the diagonal elements of the tridiagonal matrix, and e[1..n] the off-diagonal elements, with e[1]=0. Severalstatements, as noted in comments, can be omitted if only eigenvalues are to be found, in whichcase a contains no useful information on output.

Otherwise they are to be included.{int l,k,j,i;float scale,hh,h,g,f;for (i=n;i>=2;i--) {l=i-1;h=scale=0.0;if (l > 1) {for (k=1;k<=l;k++)scale += fabs(a[i][k]);if (scale == 0.0)Skip transformation.e[i]=a[i][l];else {for (k=1;k<=l;k++) {a[i][k] /= scale;Use scaled a’s for transformation.h += a[i][k]*a[i][k];Form σ in h.}f=a[i][l];g=(f >= 0.0 ? -sqrt(h) : sqrt(h));e[i]=scale*g;h -= f*g;Now h is equation (11.2.4).a[i][l]=f-g;Store u in the ith row of a.f=0.0;for (j=1;j<=l;j++) {/* Next statement can be omitted if eigenvectors not wanted */a[j][i]=a[i][j]/h;Store u/H in ith column of a.g=0.0;Form an element of A · u in g.for (k=1;k<=j;k++)g += a[j][k]*a[i][k];for (k=j+1;k<=l;k++)g += a[k][j]*a[i][k];e[j]=g/h;Form element of p in temporarily unusedelement of e.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).and use the scaled variables for the transformation. (A Householder transformationdepends only on the ratios of the elements.)Note that when dealing with a matrix whose elements vary over many ordersof magnitude, it is important that the matrix be permuted, insofar as possible, so thatthe smaller elements are in the top left-hand corner.

This is because the reductionis performed starting from the bottom right-hand corner, and a mixture of small andlarge elements there can lead to considerable rounding errors.The routine tred2 is designed for use with the routine tqli of the next section.tqli finds the eigenvalues and eigenvectors of a symmetric, tridiagonal matrix.The combination of tred2 and tqli is the most efficient known technique forfinding all the eigenvalues and eigenvectors (or just all the eigenvalues) of a real,symmetric matrix.In the listing below, the statements indicated by comments are required only forsubsequent computation of eigenvectors.

If only eigenvalues are required, omissionof the commented statements speeds up the execution time of tred2 by a factor of 2for large n. In the limit of large n, the operation count of the Householder reductionis 2n3 /3 for eigenvalues only, and 4n3 /3 for both eigenvalues and eigenvectors.11.3 Eigenvalues and Eigenvectors of a Tridiagonal Matrix475}} elsee[i]=a[i][l];d[i]=h;}/* Next statement can be omitted if eigenvectors not wanted */d[1]=0.0;e[1]=0.0;/* Contents of this loop can be omitted if eigenvectors notwanted except for statement d[i]=a[i][i]; */for (i=1;i<=n;i++) {Begin accumulation of transformation mal=i-1;trices.if (d[i]) {This block skipped when i=1.for (j=1;j<=l;j++) {g=0.0;for (k=1;k<=l;k++)Use u and u/H stored in a to form P·Q.g += a[i][k]*a[k][j];for (k=1;k<=l;k++)a[k][j] -= g*a[k][i];}}d[i]=a[i][i];This statement remains.a[i][i]=1.0;Reset row and column of a to identityfor (j=1;j<=l;j++) a[j][i]=a[i][j]=0.0;matrix for next iteration.}}CITED REFERENCES AND FURTHER READING:Golub, G.H., and Van Loan, C.F.

1989, Matrix Computations, 2nd ed. (Baltimore: Johns HopkinsUniversity Press), §5.1. [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).Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag). [2]11.3 Eigenvalues and Eigenvectors of aTridiagonal MatrixEvaluation of the Characteristic PolynomialOnce our original, real, symmetric matrix has been reduced to tridiagonal form,one possible way to determine its eigenvalues is to find the roots of the characteristicpolynomial pn (λ) directly. The characteristic polynomial of a tridiagonal matrix canbe evaluated for any trial value of λ by an efficient recursion relation (see [1], forexample).

The polynomials of lower degree produced during the recurrence form aSample 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).f += e[j]*a[i][j];}hh=f/(h+h);Form K, equation (11.2.11).for (j=1;j<=l;j++) {Form q and store in e overwriting p.f=a[i][j];e[j]=g=e[j]-hh*f;for (k=1;k<=j;k++)Reduce a, equation (11.2.13).a[j][k] -= (f*e[k]+g*a[i][k]);}.