|
template<class C>
Definition at line 715 of file Matrix.h. References SIGN. Referenced by CxMatrixTem< C >::klm(). 00716 { 00717 00718 int m,l,iter,i,k; 00719 double s,r,p,g,f,dd,c,b; 00720 00721 int n = e.Size() ; 00722 00723 for (i=1; i<n; i++) e[i-1]=e[i]; //renumber 00724 e[n-1] = 0.0; 00725 00726 00727 for (l=0; l<n; l++) 00728 { 00729 iter=0; 00730 do 00731 { 00732 for (m=l; m<n-1; m++) 00733 { 00734 dd = fabs(d[m]) + fabs(d[m+1]); 00735 if ( (fabs(e[m]) + dd) == dd) break; 00736 } 00737 if (m != l) 00738 { 00739 if (iter++ == 30) printf("Too many iterations in tqli\n"); 00740 g = (d[l+1]-d[l])/(2.0*e[l]); 00741 r = sqrt(g*g + 1.0);//=pythag(g,1.0); 00742 g = d[m] - d[l] + e[l]/(g+SIGN(r,g)); 00743 s = c = 1.0; 00744 p = 0.0; 00745 for (i=m-1; i>=l; i--) 00746 { 00747 f = s*e[i]; 00748 b = c*e[i]; 00749 e[i+1] = (r=sqrt(f*f+g*g));//=(r=pythag(f,g)); 00750 if (r == 0.0) 00751 { 00752 d[i+1] -= p; 00753 e[m] = 0.0; 00754 break; 00755 } 00756 s = f/r; 00757 c = g/r; 00758 g = d[i+1]-p; 00759 r = (d[i]-g)*s + 2.0*c*b; 00760 d[i+1] = g+(p=s*r); 00761 g = c*r-b; 00762 00763 for (k=0; k<n; k++) 00764 { 00765 f = z[k][i+1]; 00766 z[k][i+1] = s*z[k][i] + c*f; 00767 z[k][i] = c*z[k][i] - s*f; 00768 } 00769 } 00770 if (r==0.0 && i>=l) continue; 00771 d[l] -= p; 00772 e[l] = g; 00773 e[m] = 0.0; 00774 } 00775 } while (m != l); 00776 } 00777 00778 //eigenvalues d[0..n-1], eigenvectors z[0..n-1][0..n-1], sort them now 00779 int j; 00780 for (i=0; i<n-1; i++) 00781 { 00782 p = d[k=i]; 00783 for (j=i+1; j<n; j++) if (d[j]>=p) p = d[k=j]; 00784 if (k!=i) { //insert 00785 d[k] = d[i]; 00786 d[i] = p; 00787 //exchange corresponding eigenvectors 00788 for (j=0; j<n; j++) 00789 { 00790 p = z[j][i]; 00791 z[j][i] = z[j][k]; 00792 z[j][k] = p; 00793 } 00794 } 00795 } 00796 }
|
1.5.1