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