[ 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) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.9.0 (Tue Nov 6 2012)