template<class MatrixT, class VectorT>
| void Impala::Core::Matrix::Tqli | ( | VectorT * | d, | |
| VectorT * | e, | |||
| MatrixT * | z | |||
| ) |
Definition at line 122 of file MatKLM.h.
References MatE(), SIGN, VecE(), and VecNrElem().
Referenced by MatKLM().
00123 { 00124 00125 int m,l,iter,i,k; 00126 double s,r,p,g,f,dd,c,b; 00127 00128 int n = VecNrElem(e); 00129 00130 for (i=1; i<n; i++) 00131 *VecE(e, i-1) = *VecE(e, i); //renumber 00132 *VecE(e, n-1) = 0.0; 00133 00134 00135 for (l=0; l<n; l++) 00136 { 00137 iter=0; 00138 do { 00139 for (m=l; m<n-1; m++) 00140 { 00141 dd = fabs(*VecE(d, m)) + fabs(*VecE(d, m+1)); 00142 if ( (fabs(*VecE(e, m)) + dd) == dd) 00143 break; 00144 } 00145 if (m != l) 00146 { 00147 if (iter++ == 30) 00148 std::cout << "Tqli: Too many iterations" << std::endl; 00149 g = (*VecE(d, l+1) - *VecE(d, l)) / (2.0 * *VecE(e, l)); 00150 r = sqrt(g*g + 1.0);//=pythag(g,1.0); 00151 g = *VecE(d, m) - *VecE(d, l) + *VecE(e, l)/(g+SIGN(r,g)); 00152 s = c = 1.0; 00153 p = 0.0; 00154 for (i=m-1; i>=l; i--) 00155 { 00156 f = s * *VecE(e, i); 00157 b = c * *VecE(e, i); 00158 *VecE(e, i+1) = (r=sqrt(f*f+g*g));//=(r=pythag(f,g)); 00159 if (r == 0.0) 00160 { 00161 *VecE(d, i+1) -= p; 00162 *VecE(e, m) = 0.0; 00163 break; 00164 } 00165 s = f/r; 00166 c = g/r; 00167 g = *VecE(d, i+1)-p; 00168 r = (*VecE(d, i)-g)*s + 2.0*c*b; 00169 *VecE(d, i+1) = g+(p=s*r); 00170 g = c*r-b; 00171 00172 for (k=0; k<n; k++) 00173 { 00174 f = *MatE(z, k, i+1); 00175 *MatE(z, k, i+1) = s * *MatE(z, k, i) + c*f; 00176 *MatE(z, k, i) = c * *MatE(z, k, i) - s*f; 00177 } 00178 } 00179 if (r==0.0 && i>=l) 00180 continue; 00181 *VecE(d, l) -= p; 00182 *VecE(e, l) = g; 00183 *VecE(e, m) = 0.0; 00184 } 00185 } while (m != l); 00186 } 00187 00188 //eigenvalues d[0..n-1], eigenvectors z[0..n-1][0..n-1], sort them now 00189 int j; 00190 for (i=0; i<n-1; i++) 00191 { 00192 p = *VecE(d, k=i); 00193 for (j=i+1; j<n; j++) 00194 if (*VecE(d, j) >= p) 00195 p = *VecE(d, k=j); 00196 if (k!=i) 00197 { //insert 00198 *VecE(d, k) = *VecE(d, i); 00199 *VecE(d, i) = p; 00200 //exchange corresponding eigenvectors 00201 for (j=0; j<n; j++) 00202 { 00203 p = *MatE(z, j, i); 00204 *MatE(z, j, i) = *MatE(z, j, k); 00205 *MatE(z, j, k) = p; 00206 } 00207 } 00208 } 00209 }
Here is the call graph for this function:
