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

details vigra/eigensystem.hxx VIGRA

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

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.3.3 (18 Aug 2005)