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