《数字信号处理》教学参考资料（Numerical Recipes in C，The Art of Scientific Computing Second Edition）Chapter 11.1

11.1 Jacobi Transformations of a Symmetric Matrix 463 CITED REFERENCES AND FURTHER READING: Stoer,J.,and Bulirsch,R.1980,Introduction to Numerical Analysis(New York:Springer-Verlag). Chapter 6.[1] Wilkinson,J.H.,and Reinsch,C.1971,Linear Algebra,vol.ll of Handbook for Automatic Com- putation (New York:Springer-Verlag).[2] 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).[3] IMSL Math/Library Users Manual (IMSL Inc.,2500 CityWest Boulevard,Houston TX 77042).[4] NAG Fortran Library (Numerical Algorithms Group,256 Banbury Road,Oxford OX27DE,U.K.). Chapter F02.[5] Golub,G.H.,and Van Loan,C.F.1989,Matrix Computations,2nd ed.(Baltimore:Johns Hopkins University Press),87.7.[6] Wilkinson,J.H.1965,The Algebraic Eigenvalue Problem(New York:Oxford University Press).[7] Acton,F.S.1970,Numerica/Methods That Work;1990,corrected edition (Washington:Mathe- matical Association of America),Chapter 13. 骨家 Horn,R.A.,and Johnson,C.R.1985,Matrix Analysis(Cambridge:Cambridge University Press). RECIPES 11.1 Jacobi Transformations of a Symmetric Matrix 县% 令 The Jacobi method consists of a sequence of orthogonal similarity transforma- tions of the form of equation (11.0.14).Each transformation (a Jacobi rotation)is just a plane rotation designed to annihilate one of the off-diagonal matrix elements. 家 SCIENTIFIC Successive transformations undo previously set zeros,but the off-diagonal elements nevertheless get smaller and smaller,until the matrix is diagonal to machine preci- 可 sion.Accumulating the product of the transformations as you go gives the matrix of eigenvectors,equation(11.0.15),while the elements of the final diagonal matrix are the eigenvalues. The Jacobi method is absolutely foolproof for all real symmetric matrices.For matrices of order greater than about 10,say,the algorithm is slower,by a significant 10.621 constant factor,than the OR method we shall give in $11.3.However,the Jacobi algorithm is much simpler than the more efficient methods.We thus recommend it E喜 Numerical Recipes 43106 for matrices of moderate order,where expense is not a major consideration. The basic Jacobi rotation Pp is a matrix of the form (outside North Software. c (11.1.1) Here all the diagonal elements are unity except for the two elements c in rows (and columns)p and g.All off-diagonal elements are zero except the two elements s and -s.The numbers c and s are the cosine and sine of a rotation angle o,so c2+s2 =1

11.1 Jacobi Transformations of a Symmetric Matrix 463 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). CITED REFERENCES AND FURTHER READING: Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag), Chapter 6. [1] Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Com￾putation (New York: Springer-Verlag). [2] 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). [3] IMSL Math/Library Users Manual (IMSL Inc., 2500 CityWest Boulevard, Houston TX 77042). [4] NAG Fortran Library (Numerical Algorithms Group, 256 Banbury Road, Oxford OX27DE, U.K.), Chapter F02. [5] Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns Hopkins University Press), §7.7. [6] Wilkinson, J.H. 1965, The Algebraic Eigenvalue Problem (New York: Oxford University Press). [7] Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathe￾matical Association of America), Chapter 13. Horn, R.A., and Johnson, C.R. 1985, Matrix Analysis (Cambridge: Cambridge University Press). 11.1 Jacobi Transformations of a Symmetric Matrix The Jacobi method consists of a sequence of orthogonal similarity transforma￾tions of the form of equation (11.0.14). Each transformation (a Jacobi rotation) is just a plane rotation designed to annihilate one of the off-diagonal matrix elements. Successive transformations undo previously set zeros, but the off-diagonal elements nevertheless get smaller and smaller, until the matrix is diagonal to machine preci￾sion. Accumulating the product of the transformations as you go gives the matrix of eigenvectors, equation (11.0.15), while the elements of the final diagonal matrix are the eigenvalues. The Jacobi method is absolutely foolproof for all real symmetric matrices. For matrices of order greater than about 10, say, the algorithm is slower, by a significant constant factor, than the QR method we shall give in §11.3. However, the Jacobi algorithm is much simpler than the more efficient methods. We thus recommend it for matrices of moderate order, where expense is not a major consideration. The basic Jacobi rotation Ppq is a matrix of the form Ppq =           1 ··· c ··· s . . . 1 . . . −s ··· c ··· 1           (11.1.1) Here all the diagonal elements are unity except for the two elements c in rows (and columns) p and q. All off-diagonal elements are zero except the two elements s and −s. The numbers c and s are the cosine and sine of a rotation angle φ, so c 2 +s2 = 1

464 Chapter 11.Eigensystems A plane rotation such as (11.1.1)is used to transform the matrix A according to A'=PZ·APg (11.1.2) Now,PT.A changes only rows p and q ofA,while A.Pp changes only columns p and g.Notice that the subscripts p and g do not denote components of Pp but rather label which kind of rotation the matrix is,i.e.,which rows and columns it affects.Thus the changed elements of A in(11.1.2)are only in the p and g rows and columns indicated below: …a1p…a1g… aplappapg A'= (11.1.3) Cam ag…dgp d RECIPES ge。 o? 9 Multiplying out equation(11.1.2)and using the symmetry of A,we get the explicit formulas arp=carp-sara 9 r卡p,r卡q (11.1.4) e62A arg=cara +sarp amp=ca+sa-2scam (11.1.5) 6 ang =s2app+c2ag +2scam (11.1.6) apa=(c2-s2)am+sc(app-ana) (11.1.7) The idea of the Jacobi method is to try to zero the off-diagonal elements by a series of plane rotations.Accordingly,to set,equation(11.1.7)gives the following expression for the rotation angle Numerical 10621 c2-s2 agg-app 0=cot 26=2sc 431 (11.1.8) 2apq Recipes If we let t=s/c,the definition of can be rewritten t2+2t0-1=0 (11.1.9) The smaller root of this equation corresponds to a rotation angle less than /4 in magnitude;this choice at each stage gives the most stable reduction.Using the form of the quadratic formula with the discriminant in the denominator,we can write this smaller root as sgn(0) t=1l+√0+T (11.1.10)

464 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). A plane rotation such as (11.1.1) is used to transform the matrix A according to A
= PT pq · A · Ppq (11.1.2) Now, PT pq · A changes only rows p and q of A, while A · Ppq changes only columns p and q. Notice that the subscripts p and q do not denote components of Ppq, but rather label which kind of rotation the matrix is, i.e., which rows and columns it affects. Thus the changed elements of A in (11.1.2) are only in the p and q rows and columns indicated below: A
=              ··· a
1p ··· a
1q ··· . . . . . . . . . . . . a
p1 ··· a
pp ··· a
pq ··· a
pn . . . . . . . . . . . . a
q1 ··· a
qp ··· a
qq ··· a
qn . . . . . . . . . . . . ··· a
np ··· a
nq ···              (11.1.3) Multiplying out equation (11.1.2) and using the symmetry of A, we get the explicit formulas a
rp = carp − sarq a
rq = carq + sarp r = p, r = q (11.1.4) a
pp = c2app + s2aqq − 2scapq (11.1.5) a
qq = s2app + c2aqq + 2scapq (11.1.6) a
pq = (c2 − s2)apq + sc(app − aqq) (11.1.7) The idea of the Jacobi method is to try to zero the off-diagonal elements by a series of plane rotations. Accordingly, to set a
pq = 0, equation (11.1.7) gives the following expression for the rotation angle φ θ ≡ cot 2φ ≡ c2 − s2 2sc = aqq − app 2apq (11.1.8) If we let t ≡ s/c, the definition of θ can be rewritten t 2 + 2tθ − 1=0 (11.1.9) The smaller root of this equation corresponds to a rotation angle less than π/4 in magnitude; this choice at each stage gives the most stable reduction. Using the form of the quadratic formula with the discriminant in the denominator, we can write this smaller root as t = sgn(θ) |θ| + √ θ2 + 1 (11.1.10)

11.1 Jacobi Transformations of a Symmetric Matrix 465 If is so large that 02 would overflow on the computer,we set t =1/(20).It now follows that 1 c=1 (11.1.11) V2+1 s=tc (11.1.12) When we actually use equations(11.1.4)(11.1.7)numerically,we rewrite them to minimize roundoff error.Equation(11.1.7)is replaced by apg =0 (11.1.13) 菲 The idea in the remaining equations is to set the new quantity equal to the old 、三麦&州 quantity plus a small correction.Thus we can use(11.1.7)and(11.1.13)to eliminate aag from (11.1.5),giving RECIPES dpp=app -tapa (11.1.14) Similarly, dna=ana tapa (11.1.15) arp arp-s(arq+Tarp) (11.1.16) arg=arg+s(arp-Targ） (11.1.17) where T (tan o/2)is defined by T三1+C (11.1.18) One can see the convergence of the Jacobi method by considering the sum of the squares of the off-diagonal elements S=>lar2 (11.1.19) 及 r≠8 Numerical 10521 43126 Equations (11.1.4)-(11.1.7)imply that S'=S-2lapa2 (11.1.20) (Since the transformation is orthogonal,the sum of the squares of the diagonal elements increases correspondingly by 22.)The sequence of S's thus decreases monotonically.Since the sequence is bounded below by zero,and since we can choose apa to be whatever element we want,the sequence can be made to converge to zero. Eventually one obtains a matrix D that is diagonal to machine precision.The diagonal elements give the eigenvalues of the original matrix A,since D-VT.A.V (11.1.21)

11.1 Jacobi Transformations of a Symmetric Matrix 465 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). If θ is so large that θ2 would overflow on the computer, we set t = 1/(2θ). It now follows that c = 1 √ t2 + 1 (11.1.11) s = tc (11.1.12) When we actually use equations (11.1.4)–(11.1.7) numerically, we rewrite them to minimize roundoff error. Equation (11.1.7) is replaced by a
pq =0 (11.1.13) The idea in the remaining equations is to set the new quantity equal to the old quantity plus a small correction. Thus we can use (11.1.7) and (11.1.13) to eliminate aqq from (11.1.5), giving a
pp = app − tapq (11.1.14) Similarly, a
qq = aqq + tapq (11.1.15) a
rp = arp − s(arq + τarp) (11.1.16) a
rq = arq + s(arp − τarq) (11.1.17) where τ (= tan φ/2) is defined by τ ≡ s 1 + c (11.1.18) One can see the convergence of the Jacobi method by considering the sum of the squares of the off-diagonal elements S =  r=s |ars| 2 (11.1.19) Equations (11.1.4)–(11.1.7) imply that S
= S − 2|apq| 2 (11.1.20) (Since the transformation is orthogonal, the sum of the squares of the diagonal elements increases correspondingly by 2|apq| 2.) The sequence of S’s thus decreases monotonically. Since the sequence is bounded below by zero, and since we can choose apq to be whatever element we want, the sequence can be made to converge to zero. Eventually one obtains a matrix D that is diagonal to machine precision. The diagonal elements give the eigenvalues of the original matrix A, since D = VT · A · V (11.1.21)

466 Chapter 11. Eigensystems where V=P1·P2·P3… (11.1.22) the Pi's being the successive Jacobi rotation matrices.The columns of V are the eigenvectors(since A.V=V.D).They can be computed by applying V=V.Pi (11.1.23) at each stage of calculation,where initially V is the identity matrix.In detail, equation (11.1.23)is =rs(s≠p,s≠q vhp =curp-sUrg (11.1.24) 3会 vrg=surp +curg RECIPES We rewrite these equations in terms of r as in equations(11.1.16)and (11.1.17) to minimize roundoff. The only remaining question is the strategy one should adopt for the order in which the elements are to be annihilated.Jacobi's original algorithm of 1846 searched 36 the whole upper triangle at each stage and set the largest off-diagonal element to zero. This is a reasonable strategy for hand calculation,but it is prohibitive on a computer since the search alone makes each Jacobi rotation a process oforder N2 instead of N. S￥%/A 9 A better strategy for our purposes is the cyclic Jacobi method,where one SCIENTIFIC( annihilates elements in strict order.For example,one can simply proceed down the rows:P12,P13,...,Pin;then P23,P24,etc.One can show that convergence is generally quadratic for both the original or the cyclic Jacobi methods,for nondegenerate eigenvalues.One such set of n(n-1)/2 Jacobi rotations is called a sweep. The program below,based on the implementations in [1.21,uses two further Numerica 10.621 refinements: In the first three sweeps,we carry out the pq rotation only if lapal >e 431 for some threshold value 1S0 (outside Recipes E= (11.1.25) 5n2 where So is the sum of the off-diagonal moduli. So=∑larsl (11.1.26) r<s .After four sweeps,if lapal<lappl and lapal<lagal,we set lapal =0 and skip the rotation.The criterion used in the comparison isp< 10-(where D is the number of significant decimal digits on the machine,and similarly for lagl

466 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). where V = P1 · P2 · P3 ··· (11.1.22) the Pi’s being the successive Jacobi rotation matrices. The columns of V are the eigenvectors (since A · V = V · D). They can be computed by applying V
= V · Pi (11.1.23) at each stage of calculation, where initially V is the identity matrix. In detail, equation (11.1.23) is v
rs = vrs (s = p, s = q) v
rp = cvrp − svrq v
rq = svrp + cvrq (11.1.24) We rewrite these equations in terms of τ as in equations (11.1.16) and (11.1.17) to minimize roundoff. The only remaining question is the strategy one should adopt for the order in which the elements are to be annihilated. Jacobi’s original algorithm of 1846 searched the whole upper triangle at each stage and set the largest off-diagonal element to zero. This is a reasonable strategy for hand calculation, but it is prohibitive on a computer since the search alone makes each Jacobi rotation a process of order N 2 instead of N. A better strategy for our purposes is the cyclic Jacobi method, where one annihilates elements in strict order. For example, one can simply proceed down the rows: P12, P13, ..., P1n; then P23, P24, etc. One can show that convergence is generally quadratic for both the original or the cyclic Jacobi methods, for nondegenerate eigenvalues. One such set of n(n − 1)/2 Jacobi rotations is called a sweep. The program below, based on the implementations in [1,2], uses two further refinements: • In the first three sweeps, we carry out the pq rotation only if |a pq| >  for some threshold value  = 1 5 S0 n2 (11.1.25) where S0 is the sum of the off-diagonal moduli, S0 =  r<s |ars| (11.1.26) • After four sweeps, if |apq| |app| and |apq| |aqq|, we set |apq| = 0 and skip the rotation. The criterion used in the comparison is |a pq| < 10−(D+2)|app|, where D is the number of significant decimal digits on the machine, and similarly for |aqq|

11.1 Jacobi Transformations of a Symmetric Matrix 467 In the following routine the nxn symmetric matrix a is stored as a [1..n] [1..n].On output,the superdiagonal elements of a are destroyed,but the diagonal and subdiagonal are unchanged and give full information on the original symmetric matrix a.The vector d[1..n]returns the eigenvalues of a.During the computation, it contains the current diagonal of a.The matrix v[1..n][1..n]outputs the normalized eigenvector belonging to d[k]in its kth column.The parameter nrot is the number of Jacobi rotations that were needed to achieve convergence. Typical matrices require 6 to 10 sweeps to achieve convergence,or 3n2 to 5n2 Jacobi rotations.Each rotation requires of order 4n operations,each consisting of a multiply and an add,so the total labor is of order 12n3 to 20n3operations. 8 Calculation of the eigenvectors as well as the eigenvalues changes the operation count from 4n to 6n per rotation,which is only a 50 percent overhead. #include #include "nrutil.h" #define ROTATE(a,i,j,k,1)g=a[i][j];h=a[k][1];a[i][j]=g-s*(h+g*tau);\ a[k][1]=h+s*(g-h*tau); 沙 void jacobi(float *a,int n,float d[],float **v,int *nrot) Computes all eigenvalues and eigenvectors of a real symmetric matrix a[1..n][1..n].On output,elements of a above the diagonal are destroyed.d[1..n]returns the eigenvalues of a. v[1..n][1..n]is a matrix whose columns contain,on output,the normalized eigenvectors of a.nrot returns the number of Jacobi rotations that were required. merica Press. ART 1ntj,1q,1p,1; float tresh,theta,tau,t,sm,s,h,g,c,*b,*z; Programs b=vector(1,n); z=vector(1,n); SCIENTIFIC for (ip=1;ip<=n;ip++){ Initialize to the identity matrix. for(1q=1;1q<=n;iq++)v[ip][1q]=0.0; v[ip][ip]=1.0; for (ip=1;ip<=n;ip++) Initialize b and d to the diagonal b[ip]=d[ip]=a[ip][ip]; of a. z[ip]=0.0; This vector will accumulate terms y 1920 COMPUTING(ISBN of the form tape as in equa- *nrot=0; tion(11.1.14). for(1=1;1<=50;1++)[ sm=0.0; Numerica 10621 for (ip=1;ip<=n-1;ip++){ Sum off-diagonal elements. for(1q=1p+1;1q<=n:1q++) sm +fabs(a[ip][ig]); (outside Recipes 43106 1f(sm==0.0)( The normal return,which relies free_vector(z,1,n); on quadratic convergence to Software. free_vector(b,1,n); machine underflow. North return; 2 1f(1<4) tresh=0.2*sm/(n*n); ..on the first three sweeps else tresh=0.0; ...thereafter for (ip=1;ip<=n-1;ip++) for (iq=ip+1;iq<=n;ig++){ g=100.0*fabs(a[ip][ig]); After four sweeps,skip the rotation if the off-diagonal element is small. if (i 4&&(float)(fabs(d[ip])+g)==(float)fabs(d[ip]) &&(float)(fabs(d[ig])+g)==(float)fabs(d[ig])) a[ip][iq]=0.0;

11.1 Jacobi Transformations of a Symmetric Matrix 467 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). In the following routine the n×n symmetric matrix a is stored as a[1..n] [1..n]. On output, the superdiagonal elements of a are destroyed, but the diagonal and subdiagonal are unchanged and give full information on the original symmetric matrix a. The vector d[1..n] returns the eigenvalues of a. During the computation, it contains the current diagonal of a. The matrix v[1..n][1..n] outputs the normalized eigenvector belonging to d[k] in its kth column. The parameter nrot is the number of Jacobi rotations that were needed to achieve convergence. Typical matrices require 6 to 10 sweeps to achieve convergence, or 3n 2 to 5n2 Jacobi rotations. Each rotation requires of order 4n operations, each consisting of a multiply and an add, so the total labor is of order 12n3 to 20n3 operations. Calculation of the eigenvectors as well as the eigenvalues changes the operation count from 4n to 6n per rotation, which is only a 50 percent overhead. #include #include "nrutil.h" #define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\ a[k][l]=h+s*(g-h*tau); void jacobi(float **a, int n, float d[], float **v, int *nrot) Computes all eigenvalues and eigenvectors of a real symmetric matrix a[1..n][1..n]. On output, elements of a above the diagonal are destroyed. d[1..n] returns the eigenvalues of a. v[1..n][1..n] is a matrix whose columns contain, on output, the normalized eigenvectors of a. nrot returns the number of Jacobi rotations that were required. { int j,iq,ip,i; float tresh,theta,tau,t,sm,s,h,g,c,*b,*z; b=vector(1,n); z=vector(1,n); for (ip=1;ip 4 && (float)(fabs(d[ip])+g) == (float)fabs(d[ip]) && (float)(fabs(d[iq])+g) == (float)fabs(d[iq])) a[ip][iq]=0.0;

468 Chapter 11.Eigensystems else if (fabs(a[ip][ig])>tresh){ h=d[ig]-d[ip]; if ((float)(fabs(h)+g)==(float)fabs(h)) t=(a[ip][ig])/h; t=1/(28) else theta=0.5*h/(a[ip][iq]); Equation (11.1.10). t=1.0/(fabs(theta)+sqrt(1.0+theta*theta)); if (theta 0.0)t =-t; y c=1.0/sqrt(1+t*t); s=t*c; tau=s/(1.0+c); http://www.nr. h=t*a[ip][ig]; z[ip]-h; z[1q]+=h 83g d[ip]-=h; granted for (including 19881992 d[ig]+=h; a[ip][1q]=0.0; 11.800 for (j-1;j<=ip-1;j++){ Case of rotations 1<j<p. ROTATE(a,j,ip,j,iq) for (j=ip+1;j<=iq-1;j++){ Case of rotations p<j<q. ROTATE(a,ip,j,j,iq) (Nort server for (j=iq+1;j<=n;j++) Case of rotations g<j<n. ROTATE(a,ip,j,iq,j) America computer University Press. 令 THE for(j=1;j<=n;j++){ one paper ART ROTATE(v,j,ip,j,iq) ++(*nrot); Programs for (ip=1;ip<=n;ip++){ b[ip]+z[ip]; Update d with the sum of tapa, to dir d[ip]=b[ip]; z[ip]=0.0; and reinitialize z. 188 OF SCIENTIFIC COMPUTING(ISBN nrerror("Too many iterations in routine jacobi"); 1920 v@cam Note that the above routine assumes that underflows are set to zero.On machines where this is not true,the program must be modified. Numerical Recipes 021 43108 The eigenvalues are not ordered on output.If sorting is desired,the following routine can be invoked to reorder the output of jacobi or of later routines in this (outside chapter.(The method,straight insertion,is N2 rather than Nlog N;but since you Software. have just done an N3 procedure to get the eigenvalues,you can afford yourself North this little indulgence. void eigsrt(f1oatd☐，float**v,intn) Given the eigenvalues d[1..n]and eigenvectors v[1..n][1..n]as output from jacobi ($11.1)or tqli(11.3).this routine sorts the eigenvalues into descending order,and rearranges the columns of v correspondingly.The method is straight insertion. int k,j,i; float pi for(1i=1;i<n;1++)[ p=d[k=i];

468 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). else if (fabs(a[ip][iq]) > tresh) { h=d[iq]-d[ip]; if ((float)(fabs(h)+g) == (float)fabs(h)) t=(a[ip][iq])/h; t = 1/(2θ) else { theta=0.5*h/(a[ip][iq]); Equation (11.1.10). t=1.0/(fabs(theta)+sqrt(1.0+theta*theta)); if (theta < 0.0) t = -t; } c=1.0/sqrt(1+t*t); s=t*c; tau=s/(1.0+c); h=t*a[ip][iq]; z[ip] -= h; z[iq] += h; d[ip] -= h; d[iq] += h; a[ip][iq]=0.0; for (j=1;j<=ip-1;j++) { Case of rotations 1 ≤ j<p. ROTATE(a,j,ip,j,iq) } for (j=ip+1;j<=iq-1;j++) { Case of rotations p<j<q. ROTATE(a,ip,j,j,iq) } for (j=iq+1;j<=n;j++) { Case of rotations q<j ≤ n. ROTATE(a,ip,j,iq,j) } for (j=1;j<=n;j++) { ROTATE(v,j,ip,j,iq) } ++(*nrot); } } } for (ip=1;ip<=n;ip++) { b[ip] += z[ip]; d[ip]=b[ip]; Update d with the sum of tapq, z[ip]=0.0; and reinitialize z. } } nrerror("Too many iterations in routine jacobi"); } Note that the above routine assumes that underflows are set to zero. On machines where this is not true, the program must be modified. The eigenvalues are not ordered on output. If sorting is desired, the following routine can be invoked to reorder the output of jacobi or of later routines in this chapter. (The method, straight insertion, is N 2 rather than N log N; but since you have just done an N 3 procedure to get the eigenvalues, you can afford yourself this little indulgence.) void eigsrt(float d[], float **v, int n) Given the eigenvalues d[1..n] and eigenvectors v[1..n][1..n] as output from jacobi (§11.1) or tqli (§11.3), this routine sorts the eigenvalues into descending order, and rearranges the columns of v correspondingly. The method is straight insertion. { int k,j,i; float p; for (i=1;i<n;i++) { p=d[k=i];

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 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). 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 Com￾putation (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