492 Chapter 11.Eigensystems 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). for(i=m+2;1<=nn;i++)[ a[i][i-2]=0.0; 1f(11=(m+2))a[i][i-3]=0.0; for (k=m;k<=nn-1;k++){ Double QR step on rows 1 to nn and columns m to nn. if (k !m){ p=a[k][k-1]: Begin setup of Householder http://www.nr. q=a0k+1][k-1]; vector. r=0.0; 1f(kI=(nn-1))r=a[k+2][k-1]; 83 if ((x=fabs(p)+fabs(g)+fabs(r))!=0.0){ 鱼 198891992 p/=x; Scale to prevent overflow or q/=x; underflow. 1.800 (including this one) r/=x; if ((s=SIGN(sqrt(p*p+q*q+r*r),p))!-0.0){ 7423 to to tusers to make one paper by Cambridge University Press.Programs from NUMERICAL RECIPES IN 1f(k=m){ 1f(11=m) ak][k-1]=-ak]k-1] else a0k][k-1]=-g*x; (North America THE p+=s; Equations (11.6.24). x=p/s; y=q/s: only). z=r/s; q/=p: send r /=pi for (i=k;i<=nn;++) Row modification. strictly prohibited. p=a[k][j]+q*ak+1][j]; if (k !=(nn-1)){ email to directcustsen Copyright(C) p+=rta[k+2][j]; ak+2][j]-=p*z; ak+1][]-=p*y ART OF SCIENTIFIC COMPUTING(ISBN 0-521 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[1][k+1]; if (k !(nn-1)){ @cambridge.org(outside North America). To order Numerical Recipes books or 1988-1992 by Numerical Recipes p+=z*a[i][k+2] 1-431085 a[1][k+2]-=p*r; a[1][k+1]-=p*q a[i][k]-=p; Software. ying of machine while (1<nn-1) 22492 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). 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). } 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 Householder q=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 or q /= 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]; } else a[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); } }