11.2 Reduction of a Symmetric Matrix to Tridiagonal Form 469 for (j=i+1;jp)p=d[k=j]; 1f(kI=1)[ d[k]=d[i]; d[i]=p; for(j=1;j=n;j+){ p=v[j][i]; v[j][i]v[j][k]: v[j][]=p; 222 83 19881992 CITED REFERENCES AND FURTHER READING: 111800 Golub,G.H.,and Van Loan,C.F.1989,Matrix Computations,2nd ed.(Baltimore:Johns Hopkins 州N University Press),$8.4. Smith,B.T.,et al.1976,Matrix Eigensystem Routines-EISPACK Guide,2nd ed.,vol.6 of NUMERICAL RECIPES I Lecture Notes in Computer Science (New York:Springer-Verlag).[1] Wilkinson,J.H.,and Reinsch,C.1971,Linear Algebra,vol.Il of Handbook for Automatic Com- putation (New York:Springer-Verlag).[2] 11.2 Reduction of a Symmetric Matrix 9 Programs to Tridiagonal Form:Givens and IENTIFIC Householder Reductions to dir As already mentioned,the optimum strategy for finding eigenvalues and eigenvectors is,first,to reduce the matrix to a simple form,only then beginning an COMPUTING (ISBN 198918920 iterative procedure.For symmetric matrices,the preferred simple form is tridiagonal. The Givens reduction is a modification of the Jacobi method.Instead of trying to Recipes Numerica 10621 reduce the matrix all the way to diagonal form,we are content to stop when the matrix is tridiagonal.This allows the procedure to be carried out in a finite number Recipes 43106 of steps,unlike the Jacobi method,which requires iteration to convergence. Givens Method (outside North For the Givens method,we choose the rotation angle in equation(11.1.1)so as to zero an element that is not at one of the four "corners,"i.e.,not app,apa, or a in equation (11.1.3).Specifically,we first choose P23 to annihilate a31 (and,by symmetry,a13).Then we choose P24 to annihilate a41.In general,we choose the sequence P23,P24,P2n;P34,,P3n;.;Pn-1,n where Pik annihilates ak.j-1.The method works because elements such as ap and a withrprg,are linear combinations of the old quantities arp and ar by
11.2 Reduction of a Symmetric Matrix to Tridiagonal Form 469 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 machinereadable 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). for (j=i+1;j= p) p=d[k=j]; if (k != i) { d[k]=d[i]; d[i]=p; for (j=1;j<=n;j++) { p=v[j][i]; v[j][i]=v[j][k]; v[j][k]=p; } } } } CITED REFERENCES AND FURTHER READING: Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns Hopkins University Press), §8.4. Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide, 2nd ed., vol. 6 of Lecture Notes in Computer Science (New York: Springer-Verlag). [1] Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag). [2] 11.2 Reduction of a Symmetric Matrix to Tridiagonal Form: Givens and Householder Reductions As already mentioned, the optimum strategy for finding eigenvalues and eigenvectors is, first, to reduce the matrix to a simple form, only then beginning an iterative procedure. For symmetric matrices, the preferred simple form is tridiagonal. The Givens reduction is a modification of the Jacobi method. Instead of trying to reduce the matrix all the way to diagonal form, we are content to stop when the matrix is tridiagonal. This allows the procedure to be carried out in a finite number of steps, unlike the Jacobi method, which requires iteration to convergence. Givens Method For the Givens method, we choose the rotation angle in equation (11.1.1) so as to zero an element that is not at one of the four “corners,” i.e., not a pp, apq, or aqq in equation (11.1.3). Specifically, we first choose P23 to annihilate a31 (and, by symmetry, a13). Then we choose P24 to annihilate a41. In general, we choose the sequence P23, P24,..., P2n; P34,..., P3n; ... ; Pn−1,n where Pjk annihilates ak,j−1. The method works because elements such as a rp and a rq, with r = p r = q, are linear combinations of the old quantities a rp and arq, by
470 Chapter 11.Eigensystems equation(11.1.4).Thus,if arp and ar have already been set to zero,they remain zero as the reduction proceeds.Evidently,of order n2/2 rotations are required, and the number of multiplications in a straightforward implementation is of order 4n3/3,not counting those for keeping track of the product of the transformation matrices,required for the eigenvectors. The Householder method,to be discussed next,is just as stable as the Givens reduction and it is a factor of 2 more efficient,so the Givens method is not generally used.Recent work(see [11)has shown that the Givens reduction can be reformulated to reduce the number of operations by a factor of 2,and also avoid the necessity of taking square roots.This appears to make the algorithm competitive with the Householder reduction.However,this"fast Givens"reduction has to be monitored to avoid overflows,and the variables have to be periodically rescaled.There does not seem to be any compelling reason to prefer the Givens reduction over the Householder method. Householder Method RECIPESI 2 The Householder algorithm reduces an n x n symmetric matrix A to tridiagonal form by n-2 orthogonal transformations.Each transformation annihilates the required part of a whole column and whole corresponding row.The basic ingredient Press. is a Householder matrix P,which has the form P=1-2w.w (11.2.1) OF SCIENTIFIC( where w is a real vector with w2=1.(In the present notation,the outer or matrix product of two vectors,a and b is written a.bT,while the inner or scalar product of 6 the vectors is written as aT.b.)The matrix P is orthogonal,because P2=(1-2w.w).(1-2w.w) =1-4w.w+4w.(wT.w).wT (11.2.2) 10-621 =1 秀N时 Numerical Recipes 43108 Therefore P=P-1.But PT=P,and so PT=P-1,proving orthogonality. Rewrite P as (outside North P=1、u·uT H (11.2.3) where the scalar H is H (11.2.4) and u can now be any vector.Suppose x is the vector composed of the first column of A.Choose u=x干lxe1 (11.2.5)
470 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 machinereadable 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). equation (11.1.4). Thus, if arp and arq have already been set to zero, they remain zero as the reduction proceeds. Evidently, of order n2/2 rotations are required, and the number of multiplications in a straightforward implementation is of order 4n3/3, not counting those for keeping track of the product of the transformation matrices, required for the eigenvectors. The Householder method, to be discussed next, is just as stable as the Givens reduction and it is a factor of 2 more efficient, so the Givens method is not generally used. Recent work (see [1]) has shown that the Givens reduction can be reformulated to reduce the number of operations by a factor of 2, and also avoid the necessity of taking square roots. This appears to make the algorithm competitive with the Householder reduction. However, this “fast Givens” reduction has to be monitored to avoid overflows, and the variables have to be periodically rescaled. There does not seem to be any compelling reason to prefer the Givens reduction over the Householder method. Householder Method The Householder algorithm reduces an n×n symmetric matrix A to tridiagonal form by n − 2 orthogonal transformations. Each transformation annihilates the required part of a whole column and whole corresponding row. The basic ingredient is a Householder matrix P, which has the form P = 1 − 2w · wT (11.2.1) where w is a real vector with |w| 2 = 1. (In the present notation, the outer or matrix product of two vectors, a and b is written a · bT , while the inner or scalar product of the vectors is written as aT · b.) The matrix P is orthogonal, because P2 = (1 − 2w · wT ) · (1 − 2w · wT ) = 1 − 4w · wT + 4w · (wT · w) · wT = 1 (11.2.2) Therefore P = P−1. But PT = P, and so PT = P−1, proving orthogonality. Rewrite P as P =1 − u · uT H (11.2.3) where the scalar H is H ≡ 1 2 |u| 2 (11.2.4) and u can now be any vector. Suppose x is the vector composed of the first column of A. Choose u = x ∓ |x|e1 (11.2.5)
11.2 Reduction of a Symmetric Matrix to Tridiagonal Form 471 where er is the unit vector [1,0,...,0]T,and the choice of signs will be made later.Then p=x-音年kePx 2u·(x2干lxz1) =X-22干21 (11.2.6) =X一u =±xei 8= This shows that the Householder matrix P acts on a given vector x to zero all its elements except the first one. To reduce a symmetric matrix A to tridiagonal form,we choose the vector x for the first Householder matrix to be the lower n-1 elements of the first column. Then the lower n-2 elements will be zeroed: (Nort server [1 00 0 a11 a12a13· make one paper 9 computer, Americ University Press. THE 0 ART a21 是 P1·A= 0 a31 Programs :. (m-1)P1 irrelevant 10 a11 a12a13 0 192 OF SCIENTIFIC COMPUTING(ISBN (11.2.7) irrelevant 0 Fuurggoglrion Numerical Recipes 10-521 Here we have written the matrices in partitioned form,with (n-1)P denoting a 4310-5 Householder matrix with dimensions(n-1)x(n-1).The quantity k is simply plus or minus the magnitude of the vector (outside The complete orthogonal transformation is now North Software. a11 k0 0 显 A'=PA·P= 0 (11.2.8) irrelevant We have used the fact that PT =P
11.2 Reduction of a Symmetric Matrix to Tridiagonal Form 471 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 machinereadable 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). where e1 is the unit vector [1, 0,..., 0]T , and the choice of signs will be made later. Then P · x = x − u H · (x ∓ |x|e1) T · x = x − 2u · (|x| 2 ∓ |x|x1) 2|x| 2 ∓ 2|x|x1 = x − u = ±|x|e1 (11.2.6) This shows that the Householder matrix P acts on a given vector x to zero all its elements except the first one. To reduce a symmetric matrix A to tridiagonal form, we choose the vector x for the first Householder matrix to be the lower n − 1 elements of the first column. Then the lower n − 2 elements will be zeroed: P1 · A = 1 0 0 ··· 0 0 0 . . . (n−1)P1 0 · a11 a12 a13 ··· a1n a21 a31 . . . irrelevant an1 = a11 a12 a13 ··· a1n k 0 . . . irrelevant 0 (11.2.7) Here we have written the matrices in partitioned form, with (n−1)P denoting a Householder matrix with dimensions (n − 1) × (n − 1). The quantity k is simply plus or minus the magnitude of the vector [a 21,...,an1] T . The complete orthogonal transformation is now A = P · A · P = a11 k 0 ··· 0 k 0 . . . irrelevant 0 (11.2.8) We have used the fact that PT = P.
472 Chapter 11. Eigensystems Now choose the vector x for the second Householder matrix to be the bottom n-2 elements of the second column,and from it construct 10 0 0 01 0 0 P2三 00 (11.2.9) (m-2)P2 0 0 The identity block in the upper left corner insures that the tridiagonalization achieved .com or in the first step will not be spoiled by this one,while the (n-2)-dimensional Householder matrix(-2)P2 creates one additional row and column ofthe tridiagonal output.Clearly,a sequence of n-2 such transformations will reduce the matrix A to tridiagonal form. Instead of actually carrying out the matrix multiplications in P.A.P,we 彩 compute a vector ⊙ 2 p (11.2.10) Then 三9 Press. A.P=A(1-")=A-pu7 Program 9 A'=P.A.P=A-p.uT-u.pT+2Ku.uT where the scalar K is defined by CIENTIFIC 6 ur·p K= 2H (11.2.11) If we write q≡p-Ku (11.2.12) COMPUTING (ISBN 1992 then we have A'=A-q.uT-u.qT (11.2.13) 、君0m Numerical 10621 This is the computationally useful formula. 43106 Following [2],the routine for Householder reduction given below actually starts in the nth column of A,not the first as in the explanation above.In detail,the equations are as follows:At stage m (m =1,2,...,n-2)the vector u has the form uT=[a1,02,,ai,i-2,ai,i-1±VG,0,,0] (11.2.14) Here i≡n-m+1=n,n-1,…,3 (11.2.15) and the quantity o(2 in our earlier notation)is 0=(a)2+…+(a.-1)2 (11.2.16) We choose the sign of a in (11.2.14)to be the same as the sign of ai.-1 to lessen roundoff error
472 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 machinereadable 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). Now choose the vector x for the second Householder matrix to be the bottom n − 2 elements of the second column, and from it construct P2 ≡ 1 0 0 ··· 0 0 1 0 ··· 0 0 0 . . . . . . (n−2)P2 0 0 (11.2.9) The identity block in the upper left corner insures that the tridiagonalization achieved in the first step will not be spoiled by this one, while the (n − 2)-dimensional Householder matrix (n−2)P2 creates one additional row and column of the tridiagonal output. Clearly, a sequence of n − 2 such transformations will reduce the matrix A to tridiagonal form. Instead of actually carrying out the matrix multiplications in P · A · P, we compute a vector p ≡ A · u H (11.2.10) Then A · P = A · (1 − u · uT H ) = A − p · uT A = P · A · P = A − p · uT − u · pT + 2Ku · uT where the scalar K is defined by K = uT · p 2H (11.2.11) If we write q ≡ p − Ku (11.2.12) then we have A = A − q · uT − u · qT (11.2.13) This is the computationally useful formula. Following [2], the routine for Householder reduction given below actually starts in the nth column of A, not the first as in the explanation above. In detail, the equations are as follows: At stage m (m = 1, 2,...,n−2) the vector u has the form uT = [ai1, ai2,...,ai,i−2, ai,i−1 ± √σ, 0,..., 0] (11.2.14) Here i ≡ n − m +1= n, n − 1,..., 3 (11.2.15) and the quantity σ (|x| 2 in our earlier notation) is σ = (ai1) 2 + ··· + (ai,i−1) 2 (11.2.16) We choose the sign of σ in (11.2.14) to be the same as the sign of a i,i−1 to lessen roundoff error.
11.2 Reduction of a Symmetric Matrix to Tridiagonal Form 473 Variables are thus computed in the following order:o,u,H,p,K,q,A'.At any stage 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 the routine in the next section),then the eigenvectors of A can be obtained by applying the accumulated transformation Q=P1·P2…Pn-2 (11.2.17) to those eigenvectors.We therefore form Q by recursion after all the P's have been determined: Qn-2=Pn-2 菲 Q1=P)Qj+1, j=n-3,1 (11.2.18) Q=Q1 Input for the routine below is the real,symmetric matrix a[1..n][1..n].On output,a contains the elements of the orthogonal matrix q.The vector d[1..n]is set to the diagonal elements of the tridiagonal matrix A'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 it 里4 is required for subsequent computations. No extra storage arrays are needed for the intermediate results.At stage m,the vectors p and q are nonzero only in elements 1,...,i(recall that i=n-m +1), 总3的 while u is nonzero only in elements 1,...,i-1.The elements of the vector e are being determined in the order n,n-1,...,so we can store p in the elements of e SCIENTIFIC( not already determined.The vector q can overwrite p once p is no longer needed. 61 We store u in the ith row of a and u/H in the ith column of a.Once the reduction is complete,we compute the matrices Q;using the quantities u and u/H that have been stored in a.Since Q;is an identity matrix in the last n-j+1 rows and columns,we only need compute its elements up to row and column n-j.These 9 can overwrite the u's and u/H's in the corresponding rows and columns of a,which are no longer required for subsequent Q's. The routine tred2,given below,includes one further refinement.If the quantity Recipes Numerica 10621 o is zero or "small"at any stage,one can skip the corresponding transformation. 431 A simple criterion,such as (outside Recipes smallest positive number representable on machine machine precision North would be fine most of the time.A more careful criterion is actually used.Define the quantity (11.2.19) k=1 If e=0 to machine precision,we skip the transformation.Otherwise we redefine aik becomes aik/e (11.2.20)
11.2 Reduction of a Symmetric Matrix to Tridiagonal Form 473 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 machinereadable 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). Variables are thus computed in the following order: σ, u, H, p, K, q, A . At any stage 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 the routine in the next section), then the eigenvectors of A can be obtained by applying the accumulated transformation Q = P1 · P2 ··· Pn−2 (11.2.17) to those eigenvectors. We therefore form Q by recursion after all the P’s have been determined: Qn−2 = Pn−2 Qj = Pj · Qj+1, j = n − 3,..., 1 Q = Q1 (11.2.18) Input for the routine below is the real, symmetric matrix a[1..n][1..n]. On output, a contains the elements of the orthogonal matrix q. The vector d[1..n] is set to the diagonal elements of the tridiagonal matrix A , 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 it is required for subsequent computations. No extra storage arrays are needed for the intermediate results. At stage m, the vectors 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 are being determined in the order n, n − 1,... , so we can store p in the elements of e not 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 reduction is complete, we compute the matrices Qj using the quantities u and u/H that have been stored in a. Since Qj is an identity matrix in the last n − j + 1 rows and columns, we only need compute its elements up to row and column n − j. These can overwrite the u’s and u/H’s in the corresponding rows and columns of a, which are 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 machine machine precision would be fine most of the time. A more careful criterion is actually used. Define the quantity = i−1 k=1 |aik| (11.2.19) If = 0 to machine precision, we skip the transformation. Otherwise we redefine aik becomes aik/ (11.2.20)
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 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=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 machinereadable 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 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 elements 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= 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
11.3 Eigenvalues and Eigenvectors of a Tridiagonal Matrix 475 f+=e[j]*a[i][j]; hh=f/(h+h); Form K,equation (11.2.11). for(j=1;j<1;j+)[ Form q and store in e overwriting p. f=a[i][j]; e[j]=gme[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]); else e[i]=a[i][1]; d[i]=h; /Next statement can be omitted if eigenvectors not wanted*/ 83g d[1]=0.0: granted for 19881992 e[1]=0.0: /Contents of this loop can be omitted if eigenvectors not 1.200 wanted except for statement d[i]=a[i][i];/ for (i=1;i<=n;i++){ Begin accumulation of transformation ma- 1=1-1; trices. if (a[i]){ This block skipped when i=1 for(j=1;j<=1;j+){ g=0.0; for(k=1;k<=1;k++) Use u and u/H stored in a to form P.Q. (North g+=a[i]k]*a[k][]: for(k=1;k<=1;k++) America computer, THE ak][j]-=g*a[k][1]; ART d[i]=a[i][i]; This statement remains. 9 Programs a[i][i]=1.0; Reset row and column of a to identity for(j=1;j<=1;j++)a[j][1]=a[i][j]=0.0; matrix for next iteration. IENTIFIC to dir CITED REFERENCES AND FURTHER READING: Golub,G.H.,and Van Loan,C.F.1989,Matrix Computations,2nd ed.(Baltimore:Johns Hopkins University Press),85.1.[1] COMPUTING(ISBN 198918920 Smith,B.T.,et al.1976,Matrix Eigensystem Routines-EISPACK Guide,2nd ed.,vol.6 of Lecture Notes in Computer Science(New York:Springer-Verlag). 10-621 Wilkinson,J.H.,and Reinsch,C.1971,Linear Algebra,vol.Il of Handbook for Automatic Com- putation (New York:Springer-Verlag).[2] Numerical Recipes 43108 11.3 Eigenvalues and Eigenvectors of a (outside Tridiagonal Matrix North Software. Amer ying of visit website Evaluation of the Characteristic Polynomial machine Once 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 characteristic polynomial pn(A)directly.The characteristic polynomial of a tridiagonal matrix can be evaluated for any trial value of A by an efficient recursion relation(see [1],for example).The polynomials of lower degree produced during the recurrence form a
11.3 Eigenvalues and Eigenvectors of a Tridiagonal Matrix 475 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 machinereadable 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). 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]); } } } else e[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 not wanted 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 identity for (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 Hopkins University Press), §5.1. [1] Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide, 2nd ed., vol. 6 of Lecture 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 a Tridiagonal Matrix Evaluation of the Characteristic Polynomial Once 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 characteristic polynomial pn(λ) directly. The characteristic polynomial of a tridiagonal matrix can be evaluated for any trial value of λ by an efficient recursion relation (see [1], for example). The polynomials of lower degree produced during the recurrence form a