SelfAdjointEigenSolver.h
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_SELFADJOINTEIGENSOLVER_H
00012 #define EIGEN_SELFADJOINTEIGENSOLVER_H
00013 
00014 #include "./Tridiagonalization.h"
00015 
00016 namespace Eigen { 
00017 
00018 template<typename _MatrixType>
00019 class GeneralizedSelfAdjointEigenSolver;
00020 
00021 namespace internal {
00022 template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues;
00023 }
00024 
00068 template<typename _MatrixType> class SelfAdjointEigenSolver
00069 {
00070   public:
00071 
00072     typedef _MatrixType MatrixType;
00073     enum {
00074       Size = MatrixType::RowsAtCompileTime,
00075       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00076       Options = MatrixType::Options,
00077       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00078     };
00079     
00081     typedef typename MatrixType::Scalar Scalar;
00082     typedef typename MatrixType::Index Index;
00083 
00090     typedef typename NumTraits<Scalar>::Real RealScalar;
00091     
00092     friend struct internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>;
00093 
00099     typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
00100     typedef Tridiagonalization<MatrixType> TridiagonalizationType;
00101 
00112     SelfAdjointEigenSolver()
00113         : m_eivec(),
00114           m_eivalues(),
00115           m_subdiag(),
00116           m_isInitialized(false)
00117     { }
00118 
00131     SelfAdjointEigenSolver(Index size)
00132         : m_eivec(size, size),
00133           m_eivalues(size),
00134           m_subdiag(size > 1 ? size - 1 : 1),
00135           m_isInitialized(false)
00136     {}
00137 
00153     SelfAdjointEigenSolver(const MatrixType& matrix, int options = ComputeEigenvectors)
00154       : m_eivec(matrix.rows(), matrix.cols()),
00155         m_eivalues(matrix.cols()),
00156         m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
00157         m_isInitialized(false)
00158     {
00159       compute(matrix, options);
00160     }
00161 
00192     SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors);
00193     
00208     SelfAdjointEigenSolver& computeDirect(const MatrixType& matrix, int options = ComputeEigenvectors);
00209 
00228     const MatrixType& eigenvectors() const
00229     {
00230       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00231       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00232       return m_eivec;
00233     }
00234 
00250     const RealVectorType& eigenvalues() const
00251     {
00252       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00253       return m_eivalues;
00254     }
00255 
00274     MatrixType operatorSqrt() const
00275     {
00276       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00277       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00278       return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
00279     }
00280 
00299     MatrixType operatorInverseSqrt() const
00300     {
00301       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00302       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00303       return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
00304     }
00305 
00310     ComputationInfo info() const
00311     {
00312       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
00313       return m_info;
00314     }
00315 
00321     static const int m_maxIterations = 30;
00322 
00323     #ifdef EIGEN2_SUPPORT
00324     SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors)
00325       : m_eivec(matrix.rows(), matrix.cols()),
00326         m_eivalues(matrix.cols()),
00327         m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
00328         m_isInitialized(false)
00329     {
00330       compute(matrix, computeEigenvectors);
00331     }
00332     
00333     SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
00334         : m_eivec(matA.cols(), matA.cols()),
00335           m_eivalues(matA.cols()),
00336           m_subdiag(matA.cols() > 1 ? matA.cols() - 1 : 1),
00337           m_isInitialized(false)
00338     {
00339       static_cast<GeneralizedSelfAdjointEigenSolver<MatrixType>*>(this)->compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
00340     }
00341     
00342     void compute(const MatrixType& matrix, bool computeEigenvectors)
00343     {
00344       compute(matrix, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
00345     }
00346 
00347     void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
00348     {
00349       compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
00350     }
00351     #endif // EIGEN2_SUPPORT
00352 
00353   protected:
00354     MatrixType m_eivec;
00355     RealVectorType m_eivalues;
00356     typename TridiagonalizationType::SubDiagonalType m_subdiag;
00357     ComputationInfo m_info;
00358     bool m_isInitialized;
00359     bool m_eigenvectorsOk;
00360 };
00361 
00378 namespace internal {
00379 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
00380 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
00381 }
00382 
00383 template<typename MatrixType>
00384 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
00385 ::compute(const MatrixType& matrix, int options)
00386 {
00387   eigen_assert(matrix.cols() == matrix.rows());
00388   eigen_assert((options&~(EigVecMask|GenEigMask))==0
00389           && (options&EigVecMask)!=EigVecMask
00390           && "invalid option parameter");
00391   bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
00392   Index n = matrix.cols();
00393   m_eivalues.resize(n,1);
00394 
00395   if(n==1)
00396   {
00397     m_eivalues.coeffRef(0,0) = internal::real(matrix.coeff(0,0));
00398     if(computeEigenvectors)
00399       m_eivec.setOnes(n,n);
00400     m_info = Success;
00401     m_isInitialized = true;
00402     m_eigenvectorsOk = computeEigenvectors;
00403     return *this;
00404   }
00405 
00406   // declare some aliases
00407   RealVectorType& diag = m_eivalues;
00408   MatrixType& mat = m_eivec;
00409 
00410   // map the matrix coefficients to [-1:1] to avoid over- and underflow.
00411   RealScalar scale = matrix.cwiseAbs().maxCoeff();
00412   if(scale==RealScalar(0)) scale = RealScalar(1);
00413   mat = matrix / scale;
00414   m_subdiag.resize(n-1);
00415   internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
00416   
00417   Index end = n-1;
00418   Index start = 0;
00419   Index iter = 0; // total number of iterations
00420 
00421   while (end>0)
00422   {
00423     for (Index i = start; i<end; ++i)
00424       if (internal::isMuchSmallerThan(internal::abs(m_subdiag[i]),(internal::abs(diag[i])+internal::abs(diag[i+1]))))
00425         m_subdiag[i] = 0;
00426 
00427     // find the largest unreduced block
00428     while (end>0 && m_subdiag[end-1]==0)
00429     {
00430       end--;
00431     }
00432     if (end<=0)
00433       break;
00434 
00435     // if we spent too many iterations, we give up
00436     iter++;
00437     if(iter > m_maxIterations * n) break;
00438 
00439     start = end - 1;
00440     while (start>0 && m_subdiag[start-1]!=0)
00441       start--;
00442 
00443     internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n);
00444   }
00445 
00446   if (iter <= m_maxIterations * n)
00447     m_info = Success;
00448   else
00449     m_info = NoConvergence;
00450 
00451   // Sort eigenvalues and corresponding vectors.
00452   // TODO make the sort optional ?
00453   // TODO use a better sort algorithm !!
00454   if (m_info == Success)
00455   {
00456     for (Index i = 0; i < n-1; ++i)
00457     {
00458       Index k;
00459       m_eivalues.segment(i,n-i).minCoeff(&k);
00460       if (k > 0)
00461       {
00462         std::swap(m_eivalues[i], m_eivalues[k+i]);
00463         if(computeEigenvectors)
00464           m_eivec.col(i).swap(m_eivec.col(k+i));
00465       }
00466     }
00467   }
00468   
00469   // scale back the eigen values
00470   m_eivalues *= scale;
00471 
00472   m_isInitialized = true;
00473   m_eigenvectorsOk = computeEigenvectors;
00474   return *this;
00475 }
00476 
00477 
00478 namespace internal {
00479   
00480 template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues
00481 {
00482   static inline void run(SolverType& eig, const typename SolverType::MatrixType& A, int options)
00483   { eig.compute(A,options); }
00484 };
00485 
00486 template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3,false>
00487 {
00488   typedef typename SolverType::MatrixType MatrixType;
00489   typedef typename SolverType::RealVectorType VectorType;
00490   typedef typename SolverType::Scalar Scalar;
00491   
00492   static inline void computeRoots(const MatrixType& m, VectorType& roots)
00493   {
00494     using std::sqrt;
00495     using std::atan2;
00496     using std::cos;
00497     using std::sin;
00498     const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0);
00499     const Scalar s_sqrt3 = sqrt(Scalar(3.0));
00500 
00501     // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0.  The
00502     // eigenvalues are the roots to this equation, all guaranteed to be
00503     // real-valued, because the matrix is symmetric.
00504     Scalar c0 = m(0,0)*m(1,1)*m(2,2) + Scalar(2)*m(1,0)*m(2,0)*m(2,1) - m(0,0)*m(2,1)*m(2,1) - m(1,1)*m(2,0)*m(2,0) - m(2,2)*m(1,0)*m(1,0);
00505     Scalar c1 = m(0,0)*m(1,1) - m(1,0)*m(1,0) + m(0,0)*m(2,2) - m(2,0)*m(2,0) + m(1,1)*m(2,2) - m(2,1)*m(2,1);
00506     Scalar c2 = m(0,0) + m(1,1) + m(2,2);
00507 
00508     // Construct the parameters used in classifying the roots of the equation
00509     // and in solving the equation for the roots in closed form.
00510     Scalar c2_over_3 = c2*s_inv3;
00511     Scalar a_over_3 = (c1 - c2*c2_over_3)*s_inv3;
00512     if (a_over_3 > Scalar(0))
00513       a_over_3 = Scalar(0);
00514 
00515     Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1));
00516 
00517     Scalar q = half_b*half_b + a_over_3*a_over_3*a_over_3;
00518     if (q > Scalar(0))
00519       q = Scalar(0);
00520 
00521     // Compute the eigenvalues by solving for the roots of the polynomial.
00522     Scalar rho = sqrt(-a_over_3);
00523     Scalar theta = atan2(sqrt(-q),half_b)*s_inv3;
00524     Scalar cos_theta = cos(theta);
00525     Scalar sin_theta = sin(theta);
00526     roots(0) = c2_over_3 + Scalar(2)*rho*cos_theta;
00527     roots(1) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta);
00528     roots(2) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta);
00529 
00530     // Sort in increasing order.
00531     if (roots(0) >= roots(1))
00532       std::swap(roots(0),roots(1));
00533     if (roots(1) >= roots(2))
00534     {
00535       std::swap(roots(1),roots(2));
00536       if (roots(0) >= roots(1))
00537         std::swap(roots(0),roots(1));
00538     }
00539   }
00540   
00541   static inline void run(SolverType& solver, const MatrixType& mat, int options)
00542   {
00543     using std::sqrt;
00544     eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows());
00545     eigen_assert((options&~(EigVecMask|GenEigMask))==0
00546             && (options&EigVecMask)!=EigVecMask
00547             && "invalid option parameter");
00548     bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
00549     
00550     MatrixType& eivecs = solver.m_eivec;
00551     VectorType& eivals = solver.m_eivalues;
00552   
00553     // map the matrix coefficients to [-1:1] to avoid over- and underflow.
00554     Scalar scale = mat.cwiseAbs().maxCoeff();
00555     MatrixType scaledMat = mat / scale;
00556 
00557     // compute the eigenvalues
00558     computeRoots(scaledMat,eivals);
00559 
00560     // compute the eigen vectors
00561     if(computeEigenvectors)
00562     {
00563       Scalar safeNorm2 = Eigen::NumTraits<Scalar>::epsilon();
00564       safeNorm2 *= safeNorm2;
00565       if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::epsilon())
00566       {
00567         eivecs.setIdentity();
00568       }
00569       else
00570       {
00571         scaledMat = scaledMat.template selfadjointView<Lower>();
00572         MatrixType tmp;
00573         tmp = scaledMat;
00574 
00575         Scalar d0 = eivals(2) - eivals(1);
00576         Scalar d1 = eivals(1) - eivals(0);
00577         int k =  d0 > d1 ? 2 : 0;
00578         d0 = d0 > d1 ? d1 : d0;
00579 
00580         tmp.diagonal().array () -= eivals(k);
00581         VectorType cross;
00582         Scalar n;
00583         n = (cross = tmp.row(0).cross(tmp.row(1))).squaredNorm();
00584 
00585         if(n>safeNorm2)
00586           eivecs.col(k) = cross / sqrt(n);
00587         else
00588         {
00589           n = (cross = tmp.row(0).cross(tmp.row(2))).squaredNorm();
00590 
00591           if(n>safeNorm2)
00592             eivecs.col(k) = cross / sqrt(n);
00593           else
00594           {
00595             n = (cross = tmp.row(1).cross(tmp.row(2))).squaredNorm();
00596 
00597             if(n>safeNorm2)
00598               eivecs.col(k) = cross / sqrt(n);
00599             else
00600             {
00601               // the input matrix and/or the eigenvaues probably contains some inf/NaN,
00602               // => exit
00603               // scale back to the original size.
00604               eivals *= scale;
00605 
00606               solver.m_info = NumericalIssue;
00607               solver.m_isInitialized = true;
00608               solver.m_eigenvectorsOk = computeEigenvectors;
00609               return;
00610             }
00611           }
00612         }
00613 
00614         tmp = scaledMat;
00615         tmp.diagonal().array() -= eivals(1);
00616 
00617         if(d0<=Eigen::NumTraits<Scalar>::epsilon())
00618           eivecs.col(1) = eivecs.col(k).unitOrthogonal();
00619         else
00620         {
00621           n = (cross = eivecs.col(k).cross(tmp.row(0).normalized())).squaredNorm();
00622           if(n>safeNorm2)
00623             eivecs.col(1) = cross / sqrt(n);
00624           else
00625           {
00626             n = (cross = eivecs.col(k).cross(tmp.row(1))).squaredNorm();
00627             if(n>safeNorm2)
00628               eivecs.col(1) = cross / sqrt(n);
00629             else
00630             {
00631               n = (cross = eivecs.col(k).cross(tmp.row(2))).squaredNorm();
00632               if(n>safeNorm2)
00633                 eivecs.col(1) = cross / sqrt(n);
00634               else
00635               {
00636                 // we should never reach this point,
00637                 // if so the last two eigenvalues are likely to ve very closed to each other
00638                 eivecs.col(1) = eivecs.col(k).unitOrthogonal();
00639               }
00640             }
00641           }
00642 
00643           // make sure that eivecs[1] is orthogonal to eivecs[2]
00644           Scalar d = eivecs.col(1).dot(eivecs.col(k));
00645           eivecs.col(1) = (eivecs.col(1) - d * eivecs.col(k)).normalized();
00646         }
00647 
00648         eivecs.col(k==2 ? 0 : 2) = eivecs.col(k).cross(eivecs.col(1)).normalized();
00649       }
00650     }
00651     // Rescale back to the original size.
00652     eivals *= scale;
00653     
00654     solver.m_info = Success;
00655     solver.m_isInitialized = true;
00656     solver.m_eigenvectorsOk = computeEigenvectors;
00657   }
00658 };
00659 
00660 // 2x2 direct eigenvalues decomposition, code from Hauke Heibel
00661 template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2,false>
00662 {
00663   typedef typename SolverType::MatrixType MatrixType;
00664   typedef typename SolverType::RealVectorType VectorType;
00665   typedef typename SolverType::Scalar Scalar;
00666   
00667   static inline void computeRoots(const MatrixType& m, VectorType& roots)
00668   {
00669     using std::sqrt;
00670     const Scalar t0 = Scalar(0.5) * sqrt( abs2(m(0,0)-m(1,1)) + Scalar(4)*m(1,0)*m(1,0));
00671     const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1));
00672     roots(0) = t1 - t0;
00673     roots(1) = t1 + t0;
00674   }
00675   
00676   static inline void run(SolverType& solver, const MatrixType& mat, int options)
00677   {
00678     eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
00679     eigen_assert((options&~(EigVecMask|GenEigMask))==0
00680             && (options&EigVecMask)!=EigVecMask
00681             && "invalid option parameter");
00682     bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
00683     
00684     MatrixType& eivecs = solver.m_eivec;
00685     VectorType& eivals = solver.m_eivalues;
00686   
00687     // map the matrix coefficients to [-1:1] to avoid over- and underflow.
00688     Scalar scale = mat.cwiseAbs().maxCoeff();
00689     scale = (std::max)(scale,Scalar(1));
00690     MatrixType scaledMat = mat / scale;
00691     
00692     // Compute the eigenvalues
00693     computeRoots(scaledMat,eivals);
00694     
00695     // compute the eigen vectors
00696     if(computeEigenvectors)
00697     {
00698       scaledMat.diagonal().array () -= eivals(1);
00699       Scalar a2 = abs2(scaledMat(0,0));
00700       Scalar c2 = abs2(scaledMat(1,1));
00701       Scalar b2 = abs2(scaledMat(1,0));
00702       if(a2>c2)
00703       {
00704         eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0);
00705         eivecs.col(1) /= sqrt(a2+b2);
00706       }
00707       else
00708       {
00709         eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0);
00710         eivecs.col(1) /= sqrt(c2+b2);
00711       }
00712 
00713       eivecs.col(0) << eivecs.col(1).unitOrthogonal();
00714     }
00715     
00716     // Rescale back to the original size.
00717     eivals *= scale;
00718     
00719     solver.m_info = Success;
00720     solver.m_isInitialized = true;
00721     solver.m_eigenvectorsOk = computeEigenvectors;
00722   }
00723 };
00724 
00725 }
00726 
00727 template<typename MatrixType>
00728 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
00729 ::computeDirect(const MatrixType& matrix, int options)
00730 {
00731   internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>::run(*this,matrix,options);
00732   return *this;
00733 }
00734 
00735 namespace internal {
00736 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
00737 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
00738 {
00739   RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
00740   RealScalar e = subdiag[end-1];
00741   // Note that thanks to scaling, e^2 or td^2 cannot overflow, however they can still
00742   // underflow thus leading to inf/NaN values when using the following commented code:
00743 //   RealScalar e2 = abs2(subdiag[end-1]);
00744 //   RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
00745   // This explain the following, somewhat more complicated, version:
00746   RealScalar mu = diag[end];
00747   if(td==0)
00748     mu -= abs(e);
00749   else
00750   {
00751     RealScalar e2 = abs2(subdiag[end-1]);
00752     RealScalar h = hypot(td,e);
00753     if(e2==0)  mu -= (e / (td + (td>0 ? 1 : -1))) * (e / h);
00754     else       mu -= e2 / (td + (td>0 ? h : -h));
00755   }
00756   
00757   RealScalar x = diag[start] - mu;
00758   RealScalar z = subdiag[start];
00759   for (Index k = start; k < end; ++k)
00760   {
00761     JacobiRotation<RealScalar> rot;
00762     rot.makeGivens(x, z);
00763 
00764     // do T = G' T G
00765     RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
00766     RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
00767 
00768     diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
00769     diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
00770     subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
00771     
00772 
00773     if (k > start)
00774       subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
00775 
00776     x = subdiag[k];
00777 
00778     if (k < end - 1)
00779     {
00780       z = -rot.s() * subdiag[k+1];
00781       subdiag[k + 1] = rot.c() * subdiag[k+1];
00782     }
00783     
00784     // apply the givens rotation to the unit matrix Q = Q * G
00785     if (matrixQ)
00786     {
00787       // FIXME if StorageOrder == RowMajor this operation is not very efficient
00788       Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n);
00789       q.applyOnTheRight(k,k+1,rot);
00790     }
00791   }
00792 }
00793 
00794 } // end namespace internal
00795 
00796 } // end namespace Eigen
00797 
00798 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:25:52