474 Chapter 11.Eigensystems and use the scaled variables for the transformation.(A Householder transformation depends only on the ratios of the elements. Note that when dealing with a matrix whose elements vary over many orders of magnitude,it is important that the matrix be permuted,insofar as possible,so that the smaller elements are in the top left-hand corner.This is because the reduction is performed starting from the bottom right-hand corner,and a mixture of small and large elements there can lead to considerable rounding errors. The routine tred2 is designed for use with the routine tqli of the next section. tgli finds the eigenvalues and eigenvectors of a symmetric,tridiagonal matrix. The combination of tred2 and tqli is the most efficient known technique for 8 finding all the eigenvalues and eigenvectors (or just all the eigenvalues)of a real. symmetric matrix. nted for In the listing below,the statements indicated by comments are required only for subsequent computation of eigenvectors.If only eigenvalues are required,omission of the commented statements speeds up the execution time of tred2 by a factor of 2 for large n.In the limit of large n,the operation count of the Householder reduction R以 is 2n3/3 for eigenvalues only,and 4n3/3 for both eigenvalues and eigenvectors. 令 #include <math.h> void tred2(float *a,int n,float d[],float e[]) Americ computer, Householder reduction of a real,symmetric matrix a[1..n][1..n].On output,a is replaced by the orthogonal matrix Q effecting the transformation.d[1..n]returns the diagonal ele- ments of the tridiagonal matrix,and e[1..n]the off-diagonal elements,with e[1]=0.Several statements,as noted in comments,can be omitted if only eigenvalues are to be found,in which case a contains no useful information on output.Otherwise they are to be included. SCIENTIFIC int 1,k,j,i; float scale,hh,h,g,f; to dir for(1=n;1>=2;i--){ 1=1-1； h=scale=0.0; 1f(1>1) 1788-1982 COMPUTING(ISBN for(k=1;k<=1;k++) scale +fabs(a[i][k]); if (scale ==0.0) Skip transformation Further reproduction, Numerical Recipes 10-621 e[i]=a[i][1]; else for(k=1;k<=1;k++){ 43108 a[i][k]/scale; Use scaled a's for transformation. h+=a[i]k]*a[i][k]; Form o in h. (outside f=a[i][1]: g=(f >=0.0 ?-sqrt(h)sgrt(h)); Software. e[i]=scale*g; h -f*g; Now h is equation (11.2.4). a[1][1]=f-g Store u in the ith row of a. f=0.0; for(j=1;j<=1;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]ta[i][k] for(k=j+1;k<=1;k++) g+=a[k][j]*a[i]k]: e[j]=g/h; Form element of p in temporarily unused element of e.474 Chapter 11. Eigensystems Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). and use the scaled variables for the transformation. (A Householder transformation depends only on the ratios of the elements.) Note that when dealing with a matrix whose elements vary over many orders of magnitude, it is important that the matrix be permuted, insofar as possible, so that the smaller elements are in the top left-hand corner. This is because the reduction is performed starting from the bottom right-hand corner, and a mixture of small and large 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 for finding 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 for subsequent computation of eigenvectors. If only eigenvalues are required, omission of the commented statements speeds up the execution time of tred2 by a factor of 2 for large n. In the limit of large n, the operation count of the Householder reduction is 2n3/3 for eigenvalues only, and 4n3/3 for both eigenvalues and eigenvectors. #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 replaced by the orthogonal matrix Q effecting the transformation. d[1..n] returns the diagonal ele￾ments of the tridiagonal matrix, and e[1..n] the off-diagonal elements, with e[1]=0. Several statements, as noted in comments, can be omitted if only eigenvalues are to be found, in which case 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 unused element of e