11.6 The QR Algorithm for Real Hessenberg Matrices 491 int nn,m,1,k,j,its,i,mmin; float z,y,x,w,v,u,t,s,r,q,p,anormi anorm=0.0; Compute matrix norm for possible use in lo- for(i=1:i<=n:1++) cating single small subdiagonal element. for (j=IMAX(i-1,1);j<=n;j++) anorm +fabs(a[i][j]); nn=n; t=0.0; Gets changed only by an exceptional shift. while (nn >=1){ Begin search for next eigenvalue. 1ts=0; do for(1=nn;1>=2;1-)[ Begin iteration:look for single small subdi- s=fabs(a[1-1][1-1])+fabs(a[1][1]); agonal element. if (s ==0.0)s=anorm; if((f1oat)(fabs(a[1][1-1])+s)=s){ 18881992 a[1][1-1]=0.0; break; 100 (including this one) 2 872 x=a[nn][nn]; if (1 ==nn){ One root found. 7423 from NUMERICAL RECIPES IN wr[nn]=x+t; wi[nn--]=0.0; else y=a[nn-1][nn-1]; w=a[nn][nn-1]*a[nn-1][nn]; (North America to to to make one paper THE if(1==(nn-1))[ Two roots found... /Cambridge University Press.Programs ART p=0.5*(y-x); q-p*p+v; z=sqrt(fabs(q)); x+t; send 1f(q>=0.0){ ...a real pair. z=p+SIGN(z,p); wr [nn-1]=wr [nn]=x+z; Copyright (C) if (z)wr[nn]=x-w/z; wi[nn-1]=wi[nn]=0.0; else ...a complex pair. wr [nn-1]=wr [nn]=x+p; to directcustser 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING(ISBN wi[nn-1]=-(wi[nn]=z); nn-=2; else No roots found.Continue iteration v@cambr 10-621 if (its =30)nrerror("Too many iterations in hgr"); 1f(1ts==10I11ts==20){ Form exceptional shift. t+=X; -43108 for (i=1;i<=mn;i++)a[i][i]-=x; s=fabs(a[nn][nn-1])+fabs(a[nn-1][nn-2]); y=x=0.75*s; (outside W=-0.4375**8; North Software. ++1ts; for(m=(nn-2);m>=1;m-){ Form shift and then look for z=a[m][m]; 2 consecutive small sub- 工=X=Z: diagonal elements. 8=y-z; p=(r*s-w)/a[m+1][m]+a[m][m+1]: Equation (11.6.23) qa[m+1][m+1]-z-r-s; r=a[m+2][m+1]; s=fabs(p)+fabs(q)+fabs(r); Scale to prevent overflow or p/=s; underflow. 9/as; r /s; if (m ==1)break;11.6 The QR Algorithm for Real Hessenberg Matrices 491 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). int nn,m,l,k,j,its,i,mmin; float z,y,x,w,v,u,t,s,r,q,p,anorm; anorm=0.0; Compute matrix norm for possible use in lo￾for (i=1;i<=n;i++) cating single small subdiagonal element. for (j=IMAX(i-1,1);j<=n;j++) anorm += fabs(a[i][j]); nn=n; t=0.0; Gets changed only by an exceptional shift. while (nn >= 1) { Begin search for next eigenvalue. its=0; do { for (l=nn;l>=2;l--) { Begin iteration: look for single small subdi￾s=fabs(a[l-1][l-1])+fabs(a[l][l]); agonal element. if (s == 0.0) s=anorm; if ((float)(fabs(a[l][l-1]) + s) == s) { a[l][l-1]=0.0; break; } } x=a[nn][nn]; if (l == nn) { One root found. wr[nn]=x+t; wi[nn--]=0.0; } else { y=a[nn-1][nn-1]; w=a[nn][nn-1]*a[nn-1][nn]; if (l == (nn-1)) { Two roots found... p=0.5*(y-x); q=p*p+w; z=sqrt(fabs(q)); x += t; if (q >= 0.0) { ...a real pair. z=p+SIGN(z,p); wr[nn-1]=wr[nn]=x+z; if (z) wr[nn]=x-w/z; wi[nn-1]=wi[nn]=0.0; } else { ...a complex pair. wr[nn-1]=wr[nn]=x+p; wi[nn-1]= -(wi[nn]=z); } nn -= 2; } else { No roots found. Continue iteration. if (its == 30) nrerror("Too many iterations in hqr"); if (its == 10 || its == 20) { Form exceptional shift. t += x; for (i=1;i<=nn;i++) a[i][i] -= x; s=fabs(a[nn][nn-1])+fabs(a[nn-1][nn-2]); y=x=0.75*s; w = -0.4375*s*s; } ++its; for (m=(nn-2);m>=l;m--) { Form shift and then look for 2 consecutive small sub￾diagonal elements. z=a[m][m]; r=x-z; s=y-z; p=(r*s-w)/a[m+1][m]+a[m][m+1]; Equation (11.6.23). q=a[m+1][m+1]-z-r-s; r=a[m+2][m+1]; s=fabs(p)+fabs(q)+fabs(r); Scale to prevent overflow or p /= s; underflow. q /= s; r /= s; if (m == l) break;