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