#include #define NRANSI #include "nrutil.h" void hqr(float **a, int n, float wr[], float wi[]) { int nn,m,l,k,j,its,i,mmin; float z,y,x,w,v,u,t,s,r,q,p,anorm; anorm=fabs(a[1][1]); for (i=2;i<=n;i++) for (j=(i-1);j<=n;j++) anorm += fabs(a[i][j]); nn=n; t=0.0; while (nn >= 1) { its=0; do { for (l=nn;l>=2;l--) { s=fabs(a[l-1][l-1])+fabs(a[l][l]); if (s == 0.0) s=anorm; if ((float)(fabs(a[l][l-1]) + s) == s) break; } x=a[nn][nn]; if (l == nn) { 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)) { p=0.5*(y-x); q=p*p+w; z=sqrt(fabs(q)); x += t; if (q >= 0.0) { 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 { wr[nn-1]=wr[nn]=x+p; wi[nn-1]= -(wi[nn]=z); } nn -= 2; } else { if (its == 30) nrerror("Too many iterations in hqr"); if (its == 10 || its == 20) { 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--) { z=a[m][m]; r=x-z; s=y-z; p=(r*s-w)/a[m+1][m]+a[m][m+1]; q=a[m+1][m+1]-z-r-s; r=a[m+2][m+1]; s=fabs(p)+fabs(q)+fabs(r); p /= s; q /= s; r /= s; if (m == l) break; 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; } 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++) { if (k != m) { p=a[k][k-1]; q=a[k+1][k-1]; 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; q /= x; 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; x=p/s; y=q/s; z=r/s; q /= p; r /= p; for (j=k;j<=nn;j++) { 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