[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/eigensystem.hxx

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*                  Copyright 2004 by Ullrich Koethe                    */
00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
00005 /*                                                                      */
00006 /*    This file is part of the VIGRA computer vision library.           */
00007 /*    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)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
VIGRA 1.6.0 (5 Nov 2009)