[ 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(unsigned int 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(unsigned int k = 0; k <= i; ++k)
00169             {
00170                 d(k) = a(k, i+1) / h;
00171             }
00172             for(unsigned int j = 0; j <= i; ++j)
00173             {
00174                 T g = 0.0;
00175                 for(unsigned int k = 0; k <= i; ++k)
00176                 {
00177                     g += a(k, i+1) * a(k, j);
00178                 }
00179                 for(unsigned int k = 0; k <= i; ++k)
00180                 {
00181                     a(k, j) -= g * d(k);
00182                 }
00183             }
00184         }
00185         for(unsigned int k = 0; k <= i; ++k)
00186         {
00187             a(k, i+1) = 0.0;
00188         }
00189     }
00190     for(unsigned int 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         unsigned int 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(unsigned int 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(unsigned int 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(unsigned int i = 0; i < n-1; ++i)
00314     {
00315         unsigned int k = i;
00316         T p = d(i);
00317         for(unsigned int 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(unsigned int 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> <<a href="eigensystem_8hxx-source.html">vigra/eigensystem.hxx</a>> or<br>
01003     <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><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     /** Compute the eigensystem of a square, but
01029         not necessarily symmetric matrix.
01030 
01031         \a a is a real square matrix, \a ew is a single-column matrix
01032         holding the possibly complex eigenvalues, and \a ev is a matrix of
01033         the same size as \a a whose columns are the corresponding eigenvectors.
01034         Eigenvalues will be sorted from largest to smallest magnitude.
01035         The algorithm returns <tt>false</tt> when it doesn't
01036         converge. It can be applied in-place, i.e. <tt>&a == &ev</tt> is allowed.
01037         The code of this function was adapted from JAMA.
01038 
01039     <b>\#include</b> <<a href="eigensystem_8hxx-source.html">vigra/eigensystem.hxx</a>> or<br>
01040     <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
01041         Namespaces: vigra and vigra::linalg
01042      */
01043 template <class T, class C1, class C2, class C3>
01044 bool
01045 nonsymmetricEigensystem(MultiArrayView<2, T, C1> const & a,
01046          MultiArrayView<2, std::complex<T>, C2> & ew, MultiArrayView<2, T, C3> & ev)
01047 {
01048     const MultiArrayIndex acols = columnCount(a);
01049     vigra_precondition(acols == rowCount(a),
01050         "nonsymmetricEigensystem(): square input matrix required.");
01051     vigra_precondition(1 == columnCount(ew) && acols == rowCount(ew) &&
01052                        acols == columnCount(ev) && acols == rowCount(ev),
01053         "nonsymmetricEigensystem(): matrix shape mismatch.");
01054 
01055     Matrix<T> H(a);
01056     Matrix<T> de(acols, 2);
01057     detail::nonsymmetricHessenbergReduction(H, ev);
01058     if(!detail::hessenbergQrDecomposition(H, ev, de))
01059         return false;
01060 
01061     for(MultiArrayIndex i = 0; i < acols; ++i)
01062     {
01063         ew(i,0) = std::complex<T>(de(i, 0), de(i, 1));
01064     }
01065     return true;
01066 }
01067 
01068     /** Compute the roots of a polynomial using the eigenvalue method.
01069 
01070         \a poly is a real polynomial (compatible to \ref vigra::PolynomialView),
01071         and \a roots a complex valued vector (compatible to <tt>std::vector</tt>
01072         with a <tt>value_type</tt> compatible to the type <tt>POLYNOMIAL::Complex</tt>) to which
01073         the roots are appended. The function calls \ref nonsymmetricEigensystem() with the standard
01074         companion matrix yielding the roots as eigenvalues. It returns <tt>false</tt> if
01075         it fails to converge.
01076 
01077         <b>\#include</b> <<a href="eigensystem_8hxx-source.html">vigra/eigensystem.hxx</a>> or<br>
01078         <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
01079         Namespaces: vigra and vigra::linalg
01080 
01081         \see polynomialRoots(), vigra::Polynomial
01082      */
01083 template <class POLYNOMIAL, class VECTOR>
01084 bool polynomialRootsEigenvalueMethod(POLYNOMIAL const & poly, VECTOR & roots, bool polishRoots)
01085 {
01086     typedef typename POLYNOMIAL::value_type T;
01087     typedef typename POLYNOMIAL::Real    Real;
01088     typedef typename POLYNOMIAL::Complex Complex;
01089     typedef Matrix<T> TMatrix;
01090     typedef Matrix<Complex> ComplexMatrix;
01091 
01092     int const degree = poly.order();
01093     double const eps = poly.epsilon();
01094 
01095     TMatrix inMatrix(degree, degree);
01096     for(int i = 0; i < degree; ++i)
01097         inMatrix(0, i) = -poly[degree - i - 1] / poly[degree];
01098     for(int i = 0; i < degree - 1; ++i)
01099         inMatrix(i + 1, i) = NumericTraits<T>::one();
01100     ComplexMatrix ew(degree, 1);
01101     TMatrix ev(degree, degree);
01102     bool success = nonsymmetricEigensystem(inMatrix, ew, ev);
01103     if(!success)
01104         return false;
01105     for(int i = 0; i < degree; ++i)
01106     {
01107         if(polishRoots)
01108             vigra::detail::laguerre1Root(poly, ew(i,0), 1);
01109         roots.push_back(vigra::detail::deleteBelowEpsilon(ew(i,0), eps));
01110     }
01111     std::sort(roots.begin(), roots.end(), vigra::detail::PolynomialRootCompare<Real>(eps));
01112     return true;
01113 }
01114 
01115 template <class POLYNOMIAL, class VECTOR>
01116 bool polynomialRootsEigenvalueMethod(POLYNOMIAL const & poly, VECTOR & roots)
01117 {
01118     return polynomialRootsEigenvalueMethod(poly, roots, true);
01119 }
01120 
01121     /** Compute the real roots of a real polynomial using the eigenvalue method.
01122 
01123         \a poly is a real polynomial (compatible to \ref vigra::PolynomialView),
01124         and \a roots a real valued vector (compatible to <tt>std::vector</tt>
01125         with a <tt>value_type</tt> compatible to the type <tt>POLYNOMIAL::Real</tt>) to which
01126         the roots are appended. The function calls \ref polynomialRootsEigenvalueMethod() and
01127         throws away all complex roots. It returns <tt>false</tt> if it fails to converge.
01128         The parameter <tt>polishRoots</tt> is ignored (it is only here for syntax compatibility 
01129         with polynomialRealRoots()).
01130         
01131 
01132         <b>\#include</b> <<a href="eigensystem_8hxx-source.html">vigra/eigensystem.hxx</a>> or<br>
01133         <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br>
01134         Namespaces: vigra and vigra::linalg
01135 
01136         \see polynomialRealRoots(), vigra::Polynomial
01137      */
01138 template <class POLYNOMIAL, class VECTOR>
01139 bool polynomialRealRootsEigenvalueMethod(POLYNOMIAL const & p, VECTOR & roots, bool /* polishRoots */)
01140 {
01141     typedef typename NumericTraits<typename VECTOR::value_type>::ComplexPromote Complex;
01142     ArrayVector<Complex> croots;
01143     if(!polynomialRootsEigenvalueMethod(p, croots))
01144         return false;
01145     for(unsigned int i = 0; i < croots.size(); ++i)
01146         if(croots[i].imag() == 0.0)
01147             roots.push_back(croots[i].real());
01148     return true;
01149 }
01150 
01151 template <class POLYNOMIAL, class VECTOR>
01152 bool polynomialRealRootsEigenvalueMethod(POLYNOMIAL const & p, VECTOR & roots)
01153 {
01154     return polynomialRealRootsEigenvalueMethod(p, roots, true);
01155 }
01156 
01157 
01158 //@}
01159 
01160 } // namespace linalg
01161 
01162 using linalg::symmetricEigensystem;
01163 using linalg::nonsymmetricEigensystem;
01164 using linalg::polynomialRootsEigenvalueMethod;
01165 using linalg::polynomialRealRootsEigenvalueMethod;
01166 
01167 } // namespace vigra
01168 
01169 #endif // VIGRA_EIGENSYSTEM_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.7.0 (Thu Aug 25 2011)