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