[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 2004 by Ullrich Koethe */ 00004 /* Cognitive Systems Group, University of Hamburg, Germany */ 00005 /* */ 00006 /* This file is part of the VIGRA computer vision library. */ 00007 /* ( Version 1.6.0, Aug 13 2008 ) */ 00008 /* The VIGRA Website is */ 00009 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00010 /* Please direct questions, bug reports, and contributions to */ 00011 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00012 /* vigra@informatik.uni-hamburg.de */ 00013 /* */ 00014 /* Permission is hereby granted, free of charge, to any person */ 00015 /* obtaining a copy of this software and associated documentation */ 00016 /* files (the "Software"), to deal in the Software without */ 00017 /* restriction, including without limitation the rights to use, */ 00018 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00019 /* sell copies of the Software, and to permit persons to whom the */ 00020 /* Software is furnished to do so, subject to the following */ 00021 /* conditions: */ 00022 /* */ 00023 /* The above copyright notice and this permission notice shall be */ 00024 /* included in all copies or substantial portions of the */ 00025 /* Software. */ 00026 /* */ 00027 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00028 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00029 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00030 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00031 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00032 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00033 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00034 /* OTHER DEALINGS IN THE SOFTWARE. */ 00035 /* */ 00036 /************************************************************************/ 00037 00038 00039 #ifndef VIGRA_EIGENSYSTEM_HXX 00040 #define VIGRA_EIGENSYSTEM_HXX 00041 00042 #include <algorithm> 00043 #include <complex> 00044 #include "matrix.hxx" 00045 #include "array_vector.hxx" 00046 #include "polynomial.hxx" 00047 00048 namespace vigra 00049 { 00050 00051 namespace linalg 00052 { 00053 00054 namespace detail 00055 { 00056 00057 // code adapted from JAMA 00058 // a and b will be overwritten 00059 template <class T, class C1, class C2> 00060 void 00061 housholderTridiagonalization(MultiArrayView<2, T, C1> &a, MultiArrayView<2, T, C2> &b) 00062 { 00063 const unsigned int n = rowCount(a); 00064 vigra_precondition(n == columnCount(a), 00065 "housholderTridiagonalization(): matrix must be square."); 00066 vigra_precondition(n == rowCount(b) && 2 <= columnCount(b), 00067 "housholderTridiagonalization(): matrix size mismatch."); 00068 00069 MultiArrayView<1, T, C2> d = b.bindOuter(0); 00070 MultiArrayView<1, T, C2> e = b.bindOuter(1); 00071 00072 for(unsigned int j = 0; j < n; ++j) 00073 { 00074 d(j) = a(n-1, j); 00075 } 00076 00077 // Householder reduction to tridiagonalMatrix form. 00078 00079 for(int i = n-1; i > 0; --i) 00080 { 00081 // Scale to avoid under/overflow. 00082 00083 T scale = 0.0; 00084 T h = 0.0; 00085 for(int k = 0; k < i; ++k) 00086 { 00087 scale = scale + abs(d(k)); 00088 } 00089 if(scale == 0.0) 00090 { 00091 e(i) = d(i-1); 00092 for(int j = 0; j < i; ++j) 00093 { 00094 d(j) = a(i-1, j); 00095 a(i, j) = 0.0; 00096 a(j, i) = 0.0; 00097 } 00098 } 00099 else 00100 { 00101 // Generate Householder vector. 00102 00103 for(int k = 0; k < i; ++k) 00104 { 00105 d(k) /= scale; 00106 h += sq(d(k)); 00107 } 00108 T f = d(i-1); 00109 T g = VIGRA_CSTD::sqrt(h); 00110 if(f > 0) { 00111 g = -g; 00112 } 00113 e(i) = scale * g; 00114 h -= f * g; 00115 d(i-1) = f - g; 00116 for(int j = 0; j < i; ++j) 00117 { 00118 e(j) = 0.0; 00119 } 00120 00121 // Apply similarity transformation to remaining columns. 00122 00123 for(int j = 0; j < i; ++j) 00124 { 00125 f = d(j); 00126 a(j, i) = f; 00127 g = e(j) + a(j, j) * f; 00128 for(int k = j+1; k <= i-1; ++k) 00129 { 00130 g += a(k, j) * d(k); 00131 e(k) += a(k, j) * f; 00132 } 00133 e(j) = g; 00134 } 00135 f = 0.0; 00136 for(int j = 0; j < i; ++j) 00137 { 00138 e(j) /= h; 00139 f += e(j) * d(j); 00140 } 00141 T hh = f / (h + h); 00142 for(int j = 0; j < i; ++j) 00143 { 00144 e(j) -= hh * d(j); 00145 } 00146 for(int j = 0; j < i; ++j) 00147 { 00148 f = d(j); 00149 g = e(j); 00150 for(int k = j; k <= i-1; ++k) 00151 { 00152 a(k, j) -= (f * e(k) + g * d(k)); 00153 } 00154 d(j) = a(i-1, j); 00155 a(i, j) = 0.0; 00156 } 00157 } 00158 d(i) = h; 00159 } 00160 00161 // Accumulate transformations. 00162 00163 for(unsigned int i = 0; i < n-1; ++i) 00164 { 00165 a(n-1, i) = a(i, i); 00166 a(i, i) = 1.0; 00167 T h = d(i+1); 00168 if(h != 0.0) 00169 { 00170 for(unsigned int k = 0; k <= i; ++k) 00171 { 00172 d(k) = a(k, i+1) / h; 00173 } 00174 for(unsigned int j = 0; j <= i; ++j) 00175 { 00176 T g = 0.0; 00177 for(unsigned int k = 0; k <= i; ++k) 00178 { 00179 g += a(k, i+1) * a(k, j); 00180 } 00181 for(unsigned int k = 0; k <= i; ++k) 00182 { 00183 a(k, j) -= g * d(k); 00184 } 00185 } 00186 } 00187 for(unsigned int k = 0; k <= i; ++k) 00188 { 00189 a(k, i+1) = 0.0; 00190 } 00191 } 00192 for(unsigned int j = 0; j < n; ++j) 00193 { 00194 d(j) = a(n-1, j); 00195 a(n-1, j) = 0.0; 00196 } 00197 a(n-1, n-1) = 1.0; 00198 e(0) = 0.0; 00199 } 00200 00201 // code adapted from JAMA 00202 // de and z will be overwritten 00203 template <class T, class C1, class C2> 00204 bool 00205 tridiagonalMatrixEigensystem(MultiArrayView<2, T, C1> &de, MultiArrayView<2, T, C2> &z) 00206 { 00207 const unsigned int n = rowCount(z); 00208 vigra_precondition(n == columnCount(z), 00209 "tridiagonalMatrixEigensystem(): matrix must be square."); 00210 vigra_precondition(n == rowCount(de) && 2 <= columnCount(de), 00211 "tridiagonalMatrixEigensystem(): matrix size mismatch."); 00212 00213 MultiArrayView<1, T, C2> d = de.bindOuter(0); 00214 MultiArrayView<1, T, C2> e = de.bindOuter(1); 00215 00216 for(unsigned int i = 1; i < n; i++) { 00217 e(i-1) = e(i); 00218 } 00219 e(n-1) = 0.0; 00220 00221 T f = 0.0; 00222 T tst1 = 0.0; 00223 T eps = VIGRA_CSTD::pow(2.0,-52.0); 00224 for(unsigned int l = 0; l < n; ++l) 00225 { 00226 // Find small subdiagonalMatrix element 00227 00228 tst1 = std::max(tst1, abs(d(l)) + abs(e(l))); 00229 unsigned int m = l; 00230 00231 // Original while-loop from Java code 00232 while(m < n) 00233 { 00234 if(abs(e(m)) <= eps*tst1) 00235 break; 00236 ++m; 00237 } 00238 00239 // If m == l, d(l) is an eigenvalue, 00240 // otherwise, iterate. 00241 00242 if(m > l) 00243 { 00244 int iter = 0; 00245 do 00246 { 00247 if(++iter > 50) 00248 return false; // too many iterations 00249 00250 // Compute implicit shift 00251 00252 T g = d(l); 00253 T p = (d(l+1) - g) / (2.0 * e(l)); 00254 T r = hypot(p,1.0); 00255 if(p < 0) 00256 { 00257 r = -r; 00258 } 00259 d(l) = e(l) / (p + r); 00260 d(l+1) = e(l) * (p + r); 00261 T dl1 = d(l+1); 00262 T h = g - d(l); 00263 for(unsigned int i = l+2; i < n; ++i) 00264 { 00265 d(i) -= h; 00266 } 00267 f = f + h; 00268 00269 // Implicit QL transformation. 00270 00271 p = d(m); 00272 T c = 1.0; 00273 T c2 = c; 00274 T c3 = c; 00275 T el1 = e(l+1); 00276 T s = 0.0; 00277 T s2 = 0.0; 00278 for(int i = m-1; i >= (int)l; --i) 00279 { 00280 c3 = c2; 00281 c2 = c; 00282 s2 = s; 00283 g = c * e(i); 00284 h = c * p; 00285 r = hypot(p,e(i)); 00286 e(i+1) = s * r; 00287 s = e(i) / r; 00288 c = p / r; 00289 p = c * d(i) - s * g; 00290 d(i+1) = h + s * (c * g + s * d(i)); 00291 00292 // Accumulate transformation. 00293 00294 for(unsigned int k = 0; k < n; ++k) 00295 { 00296 h = z(k, i+1); 00297 z(k, i+1) = s * z(k, i) + c * h; 00298 z(k, i) = c * z(k, i) - s * h; 00299 } 00300 } 00301 p = -s * s2 * c3 * el1 * e(l) / dl1; 00302 e(l) = s * p; 00303 d(l) = c * p; 00304 00305 // Check for convergence. 00306 00307 } while(abs(e(l)) > eps*tst1); 00308 } 00309 d(l) = d(l) + f; 00310 e(l) = 0.0; 00311 } 00312 00313 // Sort eigenvalues and corresponding vectors. 00314 00315 for(unsigned int i = 0; i < n-1; ++i) 00316 { 00317 int k = i; 00318 T p = abs(d(i)); 00319 for(unsigned int j = i+1; j < n; ++j) 00320 { 00321 T p1 = abs(d(j)); 00322 if(p < p1) 00323 { 00324 k = j; 00325 p = p1; 00326 } 00327 } 00328 if(k != i) 00329 { 00330 std::swap(d(k), d(i)); 00331 for(unsigned int j = 0; j < n; ++j) 00332 { 00333 std::swap(z(j, i), z(j, k)); 00334 } 00335 } 00336 } 00337 return true; 00338 } 00339 00340 // Nonsymmetric reduction to Hessenberg form. 00341 00342 template <class T, class C1, class C2> 00343 void nonsymmetricHessenbergReduction(MultiArrayView<2, T, C1> & H, MultiArrayView<2, T, C2> & V) 00344 { 00345 // This is derived from the Algol procedures orthes and ortran, 00346 // by Martin and Wilkinson, Handbook for Auto. Comp., 00347 // Vol.ii-Linear Algebra, and the corresponding 00348 // Fortran subroutines in EISPACK. 00349 00350 int n = rowCount(H); 00351 int low = 0; 00352 int high = n-1; 00353 ArrayVector<T> ort(n); 00354 00355 for(int m = low+1; m <= high-1; ++m) 00356 { 00357 // Scale column. 00358 00359 T scale = 0.0; 00360 for(int i = m; i <= high; ++i) 00361 { 00362 scale = scale + abs(H(i, m-1)); 00363 } 00364 if(scale != 0.0) 00365 { 00366 00367 // Compute Householder transformation. 00368 00369 T h = 0.0; 00370 for(int i = high; i >= m; --i) 00371 { 00372 ort[i] = H(i, m-1)/scale; 00373 h += sq(ort[i]); 00374 } 00375 T g = VIGRA_CSTD::sqrt(h); 00376 if(ort[m] > 0) 00377 { 00378 g = -g; 00379 } 00380 h = h - ort[m] * g; 00381 ort[m] = ort[m] - g; 00382 00383 // Apply Householder similarity transformation 00384 // H = (I-u*u'/h)*H*(I-u*u')/h) 00385 00386 for(int j = m; j < n; ++j) 00387 { 00388 T f = 0.0; 00389 for(int i = high; i >= m; --i) 00390 { 00391 f += ort[i]*H(i, j); 00392 } 00393 f = f/h; 00394 for(int i = m; i <= high; ++i) 00395 { 00396 H(i, j) -= f*ort[i]; 00397 } 00398 } 00399 00400 for(int i = 0; i <= high; ++i) 00401 { 00402 T f = 0.0; 00403 for(int j = high; j >= m; --j) 00404 { 00405 f += ort[j]*H(i, j); 00406 } 00407 f = f/h; 00408 for(int j = m; j <= high; ++j) 00409 { 00410 H(i, j) -= f*ort[j]; 00411 } 00412 } 00413 ort[m] = scale*ort[m]; 00414 H(m, m-1) = scale*g; 00415 } 00416 } 00417 00418 // Accumulate transformations (Algol's ortran). 00419 00420 for(int i = 0; i < n; ++i) 00421 { 00422 for(int j = 0; j < n; ++j) 00423 { 00424 V(i, j) = (i == j ? 1.0 : 0.0); 00425 } 00426 } 00427 00428 for(int m = high-1; m >= low+1; --m) 00429 { 00430 if(H(m, m-1) != 0.0) 00431 { 00432 for(int i = m+1; i <= high; ++i) 00433 { 00434 ort[i] = H(i, m-1); 00435 } 00436 for(int j = m; j <= high; ++j) 00437 { 00438 T g = 0.0; 00439 for(int i = m; i <= high; ++i) 00440 { 00441 g += ort[i] * V(i, j); 00442 } 00443 // Double division avoids possible underflow 00444 g = (g / ort[m]) / H(m, m-1); 00445 for(int i = m; i <= high; ++i) 00446 { 00447 V(i, j) += g * ort[i]; 00448 } 00449 } 00450 } 00451 } 00452 } 00453 00454 00455 // Complex scalar division. 00456 00457 template <class T> 00458 void cdiv(T xr, T xi, T yr, T yi, T & cdivr, T & cdivi) 00459 { 00460 T r,d; 00461 if(abs(yr) > abs(yi)) 00462 { 00463 r = yi/yr; 00464 d = yr + r*yi; 00465 cdivr = (xr + r*xi)/d; 00466 cdivi = (xi - r*xr)/d; 00467 } 00468 else 00469 { 00470 r = yr/yi; 00471 d = yi + r*yr; 00472 cdivr = (r*xr + xi)/d; 00473 cdivi = (r*xi - xr)/d; 00474 } 00475 } 00476 00477 template <class T, class C> 00478 int hessenbergQrDecompositionHelper(MultiArrayView<2, T, C> & H, int n, int l, double eps, 00479 T & p, T & q, T & r, T & s, T & w, T & x, T & y, T & z) 00480 { 00481 int m = n-2; 00482 while(m >= l) 00483 { 00484 z = H(m, m); 00485 r = x - z; 00486 s = y - z; 00487 p = (r * s - w) / H(m+1, m) + H(m, m+1); 00488 q = H(m+1, m+1) - z - r - s; 00489 r = H(m+2, m+1); 00490 s = abs(p) + abs(q) + abs(r); 00491 p = p / s; 00492 q = q / s; 00493 r = r / s; 00494 if(m == l) 00495 { 00496 break; 00497 } 00498 if(abs(H(m, m-1)) * (abs(q) + abs(r)) < 00499 eps * (abs(p) * (abs(H(m-1, m-1)) + abs(z) + 00500 abs(H(m+1, m+1))))) 00501 { 00502 break; 00503 } 00504 --m; 00505 } 00506 return m; 00507 } 00508 00509 00510 00511 // Nonsymmetric reduction from Hessenberg to real Schur form. 00512 00513 template <class T, class C1, class C2, class C3> 00514 bool hessenbergQrDecomposition(MultiArrayView<2, T, C1> & H, MultiArrayView<2, T, C2> & V, 00515 MultiArrayView<2, T, C3> & de) 00516 { 00517 00518 // This is derived from the Algol procedure hqr2, 00519 // by Martin and Wilkinson, Handbook for Auto. Comp., 00520 // Vol.ii-Linear Algebra, and the corresponding 00521 // Fortran subroutine in EISPACK. 00522 00523 // Initialize 00524 MultiArrayView<1, T, C3> d = de.bindOuter(0); 00525 MultiArrayView<1, T, C3> e = de.bindOuter(1); 00526 00527 int nn = rowCount(H); 00528 int n = nn-1; 00529 int low = 0; 00530 int high = nn-1; 00531 T eps = VIGRA_CSTD::pow(2.0, sizeof(T) == sizeof(float) 00532 ? -23.0 00533 : -52.0); 00534 T exshift = 0.0; 00535 T p=0,q=0,r=0,s=0,z=0,t,w,x,y; 00536 T norm = vigra::norm(H); 00537 00538 // Outer loop over eigenvalue index 00539 int iter = 0; 00540 while(n >= low) 00541 { 00542 00543 // Look for single small sub-diagonal element 00544 int l = n; 00545 while (l > low) 00546 { 00547 s = abs(H(l-1, l-1)) + abs(H(l, l)); 00548 if(s == 0.0) 00549 { 00550 s = norm; 00551 } 00552 if(abs(H(l, l-1)) < eps * s) 00553 { 00554 break; 00555 } 00556 --l; 00557 } 00558 00559 // Check for convergence 00560 // One root found 00561 if(l == n) 00562 { 00563 H(n, n) = H(n, n) + exshift; 00564 d(n) = H(n, n); 00565 e(n) = 0.0; 00566 --n; 00567 iter = 0; 00568 00569 // Two roots found 00570 00571 } 00572 else if(l == n-1) 00573 { 00574 w = H(n, n-1) * H(n-1, n); 00575 p = (H(n-1, n-1) - H(n, n)) / 2.0; 00576 q = p * p + w; 00577 z = VIGRA_CSTD::sqrt(abs(q)); 00578 H(n, n) = H(n, n) + exshift; 00579 H(n-1, n-1) = H(n-1, n-1) + exshift; 00580 x = H(n, n); 00581 00582 // Real pair 00583 00584 if(q >= 0) 00585 { 00586 if(p >= 0) 00587 { 00588 z = p + z; 00589 } 00590 else 00591 { 00592 z = p - z; 00593 } 00594 d(n-1) = x + z; 00595 d(n) = d(n-1); 00596 if(z != 0.0) 00597 { 00598 d(n) = x - w / z; 00599 } 00600 e(n-1) = 0.0; 00601 e(n) = 0.0; 00602 x = H(n, n-1); 00603 s = abs(x) + abs(z); 00604 p = x / s; 00605 q = z / s; 00606 r = VIGRA_CSTD::sqrt(p * p+q * q); 00607 p = p / r; 00608 q = q / r; 00609 00610 // Row modification 00611 00612 for(int j = n-1; j < nn; ++j) 00613 { 00614 z = H(n-1, j); 00615 H(n-1, j) = q * z + p * H(n, j); 00616 H(n, j) = q * H(n, j) - p * z; 00617 } 00618 00619 // Column modification 00620 00621 for(int i = 0; i <= n; ++i) 00622 { 00623 z = H(i, n-1); 00624 H(i, n-1) = q * z + p * H(i, n); 00625 H(i, n) = q * H(i, n) - p * z; 00626 } 00627 00628 // Accumulate transformations 00629 00630 for(int i = low; i <= high; ++i) 00631 { 00632 z = V(i, n-1); 00633 V(i, n-1) = q * z + p * V(i, n); 00634 V(i, n) = q * V(i, n) - p * z; 00635 } 00636 00637 // Complex pair 00638 00639 } 00640 else 00641 { 00642 d(n-1) = x + p; 00643 d(n) = x + p; 00644 e(n-1) = z; 00645 e(n) = -z; 00646 } 00647 n = n - 2; 00648 iter = 0; 00649 00650 // No convergence yet 00651 00652 } 00653 else 00654 { 00655 00656 // Form shift 00657 00658 x = H(n, n); 00659 y = 0.0; 00660 w = 0.0; 00661 if(l < n) 00662 { 00663 y = H(n-1, n-1); 00664 w = H(n, n-1) * H(n-1, n); 00665 } 00666 00667 // Wilkinson's original ad hoc shift 00668 00669 if(iter == 10) 00670 { 00671 exshift += x; 00672 for(int i = low; i <= n; ++i) 00673 { 00674 H(i, i) -= x; 00675 } 00676 s = abs(H(n, n-1)) + abs(H(n-1, n-2)); 00677 x = y = 0.75 * s; 00678 w = -0.4375 * s * s; 00679 } 00680 00681 // MATLAB's new ad hoc shift 00682 00683 if(iter == 30) 00684 { 00685 s = (y - x) / 2.0; 00686 s = s * s + w; 00687 if(s > 0) 00688 { 00689 s = VIGRA_CSTD::sqrt(s); 00690 if(y < x) 00691 { 00692 s = -s; 00693 } 00694 s = x - w / ((y - x) / 2.0 + s); 00695 for(int i = low; i <= n; ++i) 00696 { 00697 H(i, i) -= s; 00698 } 00699 exshift += s; 00700 x = y = w = 0.964; 00701 } 00702 } 00703 00704 iter = iter + 1; 00705 if(iter > 60) 00706 return false; 00707 00708 // Look for two consecutive small sub-diagonal elements 00709 int m = hessenbergQrDecompositionHelper(H, n, l, eps, p, q, r, s, w, x, y, z); 00710 for(int i = m+2; i <= n; ++i) 00711 { 00712 H(i, i-2) = 0.0; 00713 if(i > m+2) 00714 { 00715 H(i, i-3) = 0.0; 00716 } 00717 } 00718 00719 // Double QR step involving rows l:n and columns m:n 00720 00721 for(int k = m; k <= n-1; ++k) 00722 { 00723 int notlast = (k != n-1); 00724 if(k != m) { 00725 p = H(k, k-1); 00726 q = H(k+1, k-1); 00727 r = (notlast ? H(k+2, k-1) : 0.0); 00728 x = abs(p) + abs(q) + abs(r); 00729 if(x != 0.0) 00730 { 00731 p = p / x; 00732 q = q / x; 00733 r = r / x; 00734 } 00735 } 00736 if(x == 0.0) 00737 { 00738 break; 00739 } 00740 s = VIGRA_CSTD::sqrt(p * p + q * q + r * r); 00741 if(p < 0) 00742 { 00743 s = -s; 00744 } 00745 if(s != 0) 00746 { 00747 if(k != m) 00748 { 00749 H(k, k-1) = -s * x; 00750 } 00751 else if(l != m) 00752 { 00753 H(k, k-1) = -H(k, k-1); 00754 } 00755 p = p + s; 00756 x = p / s; 00757 y = q / s; 00758 z = r / s; 00759 q = q / p; 00760 r = r / p; 00761 00762 // Row modification 00763 00764 for(int j = k; j < nn; ++j) 00765 { 00766 p = H(k, j) + q * H(k+1, j); 00767 if(notlast) 00768 { 00769 p = p + r * H(k+2, j); 00770 H(k+2, j) = H(k+2, j) - p * z; 00771 } 00772 H(k, j) = H(k, j) - p * x; 00773 H(k+1, j) = H(k+1, j) - p * y; 00774 } 00775 00776 // Column modification 00777 00778 for(int i = 0; i <= std::min(n,k+3); ++i) 00779 { 00780 p = x * H(i, k) + y * H(i, k+1); 00781 if(notlast) 00782 { 00783 p = p + z * H(i, k+2); 00784 H(i, k+2) = H(i, k+2) - p * r; 00785 } 00786 H(i, k) = H(i, k) - p; 00787 H(i, k+1) = H(i, k+1) - p * q; 00788 } 00789 00790 // Accumulate transformations 00791 00792 for(int i = low; i <= high; ++i) 00793 { 00794 p = x * V(i, k) + y * V(i, k+1); 00795 if(notlast) 00796 { 00797 p = p + z * V(i, k+2); 00798 V(i, k+2) = V(i, k+2) - p * r; 00799 } 00800 V(i, k) = V(i, k) - p; 00801 V(i, k+1) = V(i, k+1) - p * q; 00802 } 00803 } // (s != 0) 00804 } // k loop 00805 } // check convergence 00806 } // while (n >= low) 00807 00808 // Backsubstitute to find vectors of upper triangular form 00809 00810 if(norm == 0.0) 00811 { 00812 return false; 00813 } 00814 00815 for(n = nn-1; n >= 0; --n) 00816 { 00817 p = d(n); 00818 q = e(n); 00819 00820 // Real vector 00821 00822 if(q == 0) 00823 { 00824 int l = n; 00825 H(n, n) = 1.0; 00826 for(int i = n-1; i >= 0; --i) 00827 { 00828 w = H(i, i) - p; 00829 r = 0.0; 00830 for(int j = l; j <= n; ++j) 00831 { 00832 r = r + H(i, j) * H(j, n); 00833 } 00834 if(e(i) < 0.0) 00835 { 00836 z = w; 00837 s = r; 00838 } 00839 else 00840 { 00841 l = i; 00842 if(e(i) == 0.0) 00843 { 00844 if(w != 0.0) 00845 { 00846 H(i, n) = -r / w; 00847 } 00848 else 00849 { 00850 H(i, n) = -r / (eps * norm); 00851 } 00852 00853 // Solve real equations 00854 00855 } 00856 else 00857 { 00858 x = H(i, i+1); 00859 y = H(i+1, i); 00860 q = (d(i) - p) * (d(i) - p) + e(i) * e(i); 00861 t = (x * s - z * r) / q; 00862 H(i, n) = t; 00863 if(abs(x) > abs(z)) 00864 { 00865 H(i+1, n) = (-r - w * t) / x; 00866 } 00867 else 00868 { 00869 H(i+1, n) = (-s - y * t) / z; 00870 } 00871 } 00872 00873 // Overflow control 00874 00875 t = abs(H(i, n)); 00876 if((eps * t) * t > 1) 00877 { 00878 for(int j = i; j <= n; ++j) 00879 { 00880 H(j, n) = H(j, n) / t; 00881 } 00882 } 00883 } 00884 } 00885 00886 // Complex vector 00887 00888 } 00889 else if(q < 0) 00890 { 00891 int l = n-1; 00892 00893 // Last vector component imaginary so matrix is triangular 00894 00895 if(abs(H(n, n-1)) > abs(H(n-1, n))) 00896 { 00897 H(n-1, n-1) = q / H(n, n-1); 00898 H(n-1, n) = -(H(n, n) - p) / H(n, n-1); 00899 } 00900 else 00901 { 00902 cdiv(0.0,-H(n-1, n),H(n-1, n-1)-p,q, H(n-1, n-1), H(n-1, n)); 00903 } 00904 H(n, n-1) = 0.0; 00905 H(n, n) = 1.0; 00906 for(int i = n-2; i >= 0; --i) 00907 { 00908 T ra,sa,vr,vi; 00909 ra = 0.0; 00910 sa = 0.0; 00911 for(int j = l; j <= n; ++j) 00912 { 00913 ra = ra + H(j, n-1) * H(i, j); 00914 sa = sa + H(j, n) * H(i, j); 00915 } 00916 w = H(i, i) - p; 00917 00918 if(e(i) < 0.0) 00919 { 00920 z = w; 00921 r = ra; 00922 s = sa; 00923 } 00924 else 00925 { 00926 l = i; 00927 if(e(i) == 0) 00928 { 00929 cdiv(-ra,-sa,w,q, H(i, n-1), H(i, n)); 00930 } 00931 else 00932 { 00933 // Solve complex equations 00934 00935 x = H(i, i+1); 00936 y = H(i+1, i); 00937 vr = (d(i) - p) * (d(i) - p) + e(i) * e(i) - q * q; 00938 vi = (d(i) - p) * 2.0 * q; 00939 if((vr == 0.0) && (vi == 0.0)) 00940 { 00941 vr = eps * norm * (abs(w) + abs(q) + 00942 abs(x) + abs(y) + abs(z)); 00943 } 00944 cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi, H(i, n-1), H(i, n)); 00945 if(abs(x) > (abs(z) + abs(q))) 00946 { 00947 H(i+1, n-1) = (-ra - w * H(i, n-1) + q * H(i, n)) / x; 00948 H(i+1, n) = (-sa - w * H(i, n) - q * H(i, n-1)) / x; 00949 } 00950 else 00951 { 00952 cdiv(-r-y*H(i, n-1),-s-y*H(i, n),z,q, H(i+1, n-1), H(i+1, n)); 00953 } 00954 } 00955 00956 // Overflow control 00957 00958 t = std::max(abs(H(i, n-1)),abs(H(i, n))); 00959 if((eps * t) * t > 1) 00960 { 00961 for(int j = i; j <= n; ++j) 00962 { 00963 H(j, n-1) = H(j, n-1) / t; 00964 H(j, n) = H(j, n) / t; 00965 } 00966 } 00967 } 00968 } 00969 } 00970 } 00971 00972 // Back transformation to get eigenvectors of original matrix 00973 00974 for(int j = nn-1; j >= low; --j) 00975 { 00976 for(int i = low; i <= high; ++i) 00977 { 00978 z = 0.0; 00979 for(int k = low; k <= std::min(j,high); ++k) 00980 { 00981 z = z + V(i, k) * H(k, j); 00982 } 00983 V(i, j) = z; 00984 } 00985 } 00986 return true; 00987 } 00988 00989 } // namespace detail 00990 00991 /** \addtogroup MatrixAlgebra 00992 */ 00993 //@{ 00994 /** Compute the eigensystem of a symmetric matrix. 00995 00996 \a a is a real symmetric matrix, \a ew is a single-column matrix 00997 holding the eigenvalues, and \a ev is a matrix of the same size as 00998 \a a whose columns are the corresponding eigenvectors. Eigenvalues 00999 will be sorted from largest to smallest magnitude. 01000 The algorithm returns <tt>false</tt> when it doesn't 01001 converge. It can be applied in-place, i.e. <tt>&a == &ev</tt> is allowed. 01002 The code of this function was adapted from JAMA. 01003 01004 <b>\#include</b> <<a href="eigensystem_8hxx-source.html">vigra/eigensystem.hxx</a>> or<br> 01005 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01006 Namespaces: vigra and vigra::linalg 01007 */ 01008 template <class T, class C1, class C2, class C3> 01009 bool 01010 symmetricEigensystem(MultiArrayView<2, T, C1> const & a, 01011 MultiArrayView<2, T, C2> & ew, MultiArrayView<2, T, C3> & ev) 01012 { 01013 vigra_precondition(isSymmetric(a), 01014 "symmetricEigensystem(): symmetric input matrix required."); 01015 unsigned int acols = columnCount(a); 01016 vigra_precondition(1 == columnCount(ew) && acols == rowCount(ew) && 01017 acols == columnCount(ev) && acols == rowCount(ev), 01018 "symmetricEigensystem(): matrix shape mismatch."); 01019 01020 ev.copy(a); // does nothing if &ev == &a 01021 Matrix<T> de(acols, 2); 01022 detail::housholderTridiagonalization(ev, de); 01023 if(!detail::tridiagonalMatrixEigensystem(de, ev)) 01024 return false; 01025 01026 ew.copy(columnVector(de, 0)); 01027 return true; 01028 } 01029 01030 /** Compute the eigensystem of a square, but 01031 not necessarily symmetric matrix. 01032 01033 \a a is a real square matrix, \a ew is a single-column matrix 01034 holding the possibly complex eigenvalues, and \a ev is a matrix of 01035 the same size as \a a whose columns are the corresponding eigenvectors. 01036 Eigenvalues will be sorted from largest to smallest magnitude. 01037 The algorithm returns <tt>false</tt> when it doesn't 01038 converge. It can be applied in-place, i.e. <tt>&a == &ev</tt> is allowed. 01039 The code of this function was adapted from JAMA. 01040 01041 <b>\#include</b> <<a href="eigensystem_8hxx-source.html">vigra/eigensystem.hxx</a>> or<br> 01042 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01043 Namespaces: vigra and vigra::linalg 01044 */ 01045 template <class T, class C1, class C2, class C3> 01046 bool 01047 nonsymmetricEigensystem(MultiArrayView<2, T, C1> const & a, 01048 MultiArrayView<2, std::complex<T>, C2> & ew, MultiArrayView<2, T, C3> & ev) 01049 { 01050 unsigned int acols = columnCount(a); 01051 vigra_precondition(acols == rowCount(a), 01052 "nonsymmetricEigensystem(): square input matrix required."); 01053 vigra_precondition(1 == columnCount(ew) && acols == rowCount(ew) && 01054 acols == columnCount(ev) && acols == rowCount(ev), 01055 "nonsymmetricEigensystem(): matrix shape mismatch."); 01056 01057 Matrix<T> H(a); 01058 Matrix<T> de(acols, 2); 01059 detail::nonsymmetricHessenbergReduction(H, ev); 01060 if(!detail::hessenbergQrDecomposition(H, ev, de)) 01061 return false; 01062 01063 for(unsigned int i=0; i < acols; ++i) 01064 { 01065 ew(i,0) = std::complex<T>(de(i, 0), de(i, 1)); 01066 } 01067 return true; 01068 } 01069 01070 /** Compute the roots of a polynomial using the eigenvalue method. 01071 01072 \a poly is a real polynomial (compatible to \ref vigra::PolynomialView), 01073 and \a roots a complex valued vector (compatible to <tt>std::vector</tt> 01074 with a <tt>value_type</tt> compatible to the type <tt>POLYNOMIAL::Complex</tt>) to which 01075 the roots are appended. The function calls \ref nonsymmetricEigensystem() with the standard 01076 companion matrix yielding the roots as eigenvalues. It returns <tt>false</tt> if 01077 it fails to converge. 01078 01079 <b>\#include</b> <<a href="eigensystem_8hxx-source.html">vigra/eigensystem.hxx</a>> or<br> 01080 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01081 Namespaces: vigra and vigra::linalg 01082 01083 \see polynomialRoots(), vigra::Polynomial 01084 */ 01085 template <class POLYNOMIAL, class VECTOR> 01086 bool polynomialRootsEigenvalueMethod(POLYNOMIAL const & poly, VECTOR & roots, bool polishRoots) 01087 { 01088 typedef typename POLYNOMIAL::value_type T; 01089 typedef typename POLYNOMIAL::Real Real; 01090 typedef typename POLYNOMIAL::Complex Complex; 01091 typedef Matrix<T> TMatrix; 01092 typedef Matrix<Complex> ComplexMatrix; 01093 01094 int const degree = poly.order(); 01095 double const eps = poly.epsilon(); 01096 01097 TMatrix inMatrix(degree, degree); 01098 for(int i = 0; i < degree; ++i) 01099 inMatrix(0, i) = -poly[degree - i - 1] / poly[degree]; 01100 for(int i = 0; i < degree - 1; ++i) 01101 inMatrix(i + 1, i) = NumericTraits<T>::one(); 01102 ComplexMatrix ew(degree, 1); 01103 TMatrix ev(degree, degree); 01104 bool success = nonsymmetricEigensystem(inMatrix, ew, ev); 01105 if(!success) 01106 return false; 01107 for(int i = 0; i < degree; ++i) 01108 { 01109 if(polishRoots) 01110 vigra::detail::laguerre1Root(poly, ew(i,0), 1); 01111 roots.push_back(vigra::detail::deleteBelowEpsilon(ew(i,0), eps)); 01112 } 01113 std::sort(roots.begin(), roots.end(), vigra::detail::PolynomialRootCompare<Real>(eps)); 01114 return true; 01115 } 01116 01117 template <class POLYNOMIAL, class VECTOR> 01118 bool polynomialRootsEigenvalueMethod(POLYNOMIAL const & poly, VECTOR & roots) 01119 { 01120 return polynomialRootsEigenvalueMethod(poly, roots, true); 01121 } 01122 01123 /** Compute the real roots of a real polynomial using the eigenvalue method. 01124 01125 \a poly is a real polynomial (compatible to \ref vigra::PolynomialView), 01126 and \a roots a real valued vector (compatible to <tt>std::vector</tt> 01127 with a <tt>value_type</tt> compatible to the type <tt>POLYNOMIAL::Real</tt>) to which 01128 the roots are appended. The function calls \ref polynomialRootsEigenvalueMethod() and 01129 throws away all complex roots. It returns <tt>false</tt> if it fails to converge. 01130 01131 <b>\#include</b> <<a href="eigensystem_8hxx-source.html">vigra/eigensystem.hxx</a>> or<br> 01132 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 01133 Namespaces: vigra and vigra::linalg 01134 01135 \see polynomialRealRoots(), vigra::Polynomial 01136 */ 01137 template <class POLYNOMIAL, class VECTOR> 01138 bool polynomialRealRootsEigenvalueMethod(POLYNOMIAL const & p, VECTOR & roots, bool polishRoots) 01139 { 01140 typedef typename NumericTraits<typename VECTOR::value_type>::ComplexPromote Complex; 01141 ArrayVector<Complex> croots; 01142 if(!polynomialRootsEigenvalueMethod(p, croots)) 01143 return false; 01144 for(unsigned int i = 0; i < croots.size(); ++i) 01145 if(croots[i].imag() == 0.0) 01146 roots.push_back(croots[i].real()); 01147 return true; 01148 } 01149 01150 template <class POLYNOMIAL, class VECTOR> 01151 bool polynomialRealRootsEigenvalueMethod(POLYNOMIAL const & p, VECTOR & roots) 01152 { 01153 return polynomialRealRootsEigenvalueMethod(p, roots, true); 01154 } 01155 01156 01157 //@} 01158 01159 } // namespace linalg 01160 01161 using linalg::symmetricEigensystem; 01162 using linalg::nonsymmetricEigensystem; 01163 using linalg::polynomialRootsEigenvalueMethod; 01164 using linalg::polynomialRealRootsEigenvalueMethod; 01165 01166 } // namespace vigra 01167 01168 #endif // VIGRA_EIGENSYSTEM_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|