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