EigenSolver.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 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_EIGENSOLVER_H
00012 #define EIGEN_EIGENSOLVER_H
00013 
00014 #include "./RealSchur.h"
00015 
00016 namespace Eigen { 
00017 
00064 template<typename _MatrixType> class EigenSolver
00065 {
00066   public:
00067 
00069     typedef _MatrixType MatrixType;
00070 
00071     enum {
00072       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00073       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00074       Options = MatrixType::Options,
00075       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00076       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00077     };
00078 
00080     typedef typename MatrixType::Scalar Scalar;
00081     typedef typename NumTraits<Scalar>::Real RealScalar;
00082     typedef typename MatrixType::Index Index;
00083 
00090     typedef std::complex<RealScalar> ComplexScalar;
00091 
00097     typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
00098 
00104     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;
00105 
00113  EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
00114 
00121     EigenSolver(Index size)
00122       : m_eivec(size, size),
00123         m_eivalues(size),
00124         m_isInitialized(false),
00125         m_eigenvectorsOk(false),
00126         m_realSchur(size),
00127         m_matT(size, size), 
00128         m_tmp(size)
00129     {}
00130 
00146     EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
00147       : m_eivec(matrix.rows(), matrix.cols()),
00148         m_eivalues(matrix.cols()),
00149         m_isInitialized(false),
00150         m_eigenvectorsOk(false),
00151         m_realSchur(matrix.cols()),
00152         m_matT(matrix.rows(), matrix.cols()), 
00153         m_tmp(matrix.cols())
00154     {
00155       compute(matrix, computeEigenvectors);
00156     }
00157 
00178     EigenvectorsType eigenvectors() const;
00179 
00198     const MatrixType& pseudoEigenvectors() const
00199     {
00200       eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00201       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00202       return m_eivec;
00203     }
00204 
00223     MatrixType pseudoEigenvalueMatrix() const;
00224 
00243     const EigenvalueType& eigenvalues() const
00244     {
00245       eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00246       return m_eivalues;
00247     }
00248 
00276     EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
00277 
00278     ComputationInfo info() const
00279     {
00280       eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00281       return m_realSchur.info();
00282     }
00283 
00284   private:
00285     void doComputeEigenvectors();
00286 
00287   protected:
00288     MatrixType m_eivec;
00289     EigenvalueType m_eivalues;
00290     bool m_isInitialized;
00291     bool m_eigenvectorsOk;
00292     RealSchur<MatrixType> m_realSchur;
00293     MatrixType m_matT;
00294 
00295     typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
00296     ColumnVectorType m_tmp;
00297 };
00298 
00299 template<typename MatrixType>
00300 MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
00301 {
00302   eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00303   Index n = m_eivalues.rows();
00304   MatrixType matD = MatrixType::Zero(n,n);
00305   for (Index i=0; i<n; ++i)
00306   {
00307     if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i))))
00308       matD.coeffRef(i,i) = internal::real(m_eivalues.coeff(i));
00309     else
00310     {
00311       matD.template block<2,2>(i,i) <<  internal::real(m_eivalues.coeff(i)), internal::imag(m_eivalues.coeff(i)),
00312                                        -internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i));
00313       ++i;
00314     }
00315   }
00316   return matD;
00317 }
00318 
00319 template<typename MatrixType>
00320 typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const
00321 {
00322   eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00323   eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00324   Index n = m_eivec.cols();
00325   EigenvectorsType matV(n,n);
00326   for (Index j=0; j<n; ++j)
00327   {
00328     if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(j)), internal::real(m_eivalues.coeff(j))) || j+1==n)
00329     {
00330       // we have a real eigen value
00331       matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
00332       matV.col(j).normalize();
00333     }
00334     else
00335     {
00336       // we have a pair of complex eigen values
00337       for (Index i=0; i<n; ++i)
00338       {
00339         matV.coeffRef(i,j)   = ComplexScalar(m_eivec.coeff(i,j),  m_eivec.coeff(i,j+1));
00340         matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
00341       }
00342       matV.col(j).normalize();
00343       matV.col(j+1).normalize();
00344       ++j;
00345     }
00346   }
00347   return matV;
00348 }
00349 
00350 template<typename MatrixType>
00351 EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
00352 {
00353   assert(matrix.cols() == matrix.rows());
00354 
00355   // Reduce to real Schur form.
00356   m_realSchur.compute(matrix, computeEigenvectors);
00357   if (m_realSchur.info() == Success)
00358   {
00359     m_matT = m_realSchur.matrixT();
00360     if (computeEigenvectors)
00361       m_eivec = m_realSchur.matrixU();
00362   
00363     // Compute eigenvalues from matT
00364     m_eivalues.resize(matrix.cols());
00365     Index i = 0;
00366     while (i < matrix.cols()) 
00367     {
00368       if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) 
00369       {
00370         m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
00371         ++i;
00372       }
00373       else
00374       {
00375         Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
00376         Scalar z = internal::sqrt(internal::abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
00377         m_eivalues.coeffRef(i)   = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
00378         m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
00379         i += 2;
00380       }
00381     }
00382     
00383     // Compute eigenvectors.
00384     if (computeEigenvectors)
00385       doComputeEigenvectors();
00386   }
00387 
00388   m_isInitialized = true;
00389   m_eigenvectorsOk = computeEigenvectors;
00390 
00391   return *this;
00392 }
00393 
00394 // Complex scalar division.
00395 template<typename Scalar>
00396 std::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi)
00397 {
00398   Scalar r,d;
00399   if (internal::abs(yr) > internal::abs(yi))
00400   {
00401       r = yi/yr;
00402       d = yr + r*yi;
00403       return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
00404   }
00405   else
00406   {
00407       r = yr/yi;
00408       d = yi + r*yr;
00409       return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
00410   }
00411 }
00412 
00413 
00414 template<typename MatrixType>
00415 void EigenSolver<MatrixType>::doComputeEigenvectors()
00416 {
00417   const Index size = m_eivec.cols();
00418   const Scalar eps = NumTraits<Scalar>::epsilon();
00419 
00420   // inefficient! this is already computed in RealSchur
00421   Scalar norm(0);
00422   for (Index j = 0; j < size; ++j)
00423   {
00424     norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
00425   }
00426   
00427   // Backsubstitute to find vectors of upper triangular form
00428   if (norm == 0.0)
00429   {
00430     return;
00431   }
00432 
00433   for (Index n = size-1; n >= 0; n--)
00434   {
00435     Scalar p = m_eivalues.coeff(n).real();
00436     Scalar q = m_eivalues.coeff(n).imag();
00437 
00438     // Scalar vector
00439     if (q == Scalar(0))
00440     {
00441       Scalar lastr(0), lastw(0);
00442       Index l = n;
00443 
00444       m_matT.coeffRef(n,n) = 1.0;
00445       for (Index i = n-1; i >= 0; i--)
00446       {
00447         Scalar w = m_matT.coeff(i,i) - p;
00448         Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
00449 
00450         if (m_eivalues.coeff(i).imag() < 0.0)
00451         {
00452           lastw = w;
00453           lastr = r;
00454         }
00455         else
00456         {
00457           l = i;
00458           if (m_eivalues.coeff(i).imag() == 0.0)
00459           {
00460             if (w != 0.0)
00461               m_matT.coeffRef(i,n) = -r / w;
00462             else
00463               m_matT.coeffRef(i,n) = -r / (eps * norm);
00464           }
00465           else // Solve real equations
00466           {
00467             Scalar x = m_matT.coeff(i,i+1);
00468             Scalar y = m_matT.coeff(i+1,i);
00469             Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
00470             Scalar t = (x * lastr - lastw * r) / denom;
00471             m_matT.coeffRef(i,n) = t;
00472             if (internal::abs(x) > internal::abs(lastw))
00473               m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
00474             else
00475               m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
00476           }
00477 
00478           // Overflow control
00479           Scalar t = internal::abs(m_matT.coeff(i,n));
00480           if ((eps * t) * t > Scalar(1))
00481             m_matT.col(n).tail(size-i) /= t;
00482         }
00483       }
00484     }
00485     else if (q < Scalar(0) && n > 0) // Complex vector
00486     {
00487       Scalar lastra(0), lastsa(0), lastw(0);
00488       Index l = n-1;
00489 
00490       // Last vector component imaginary so matrix is triangular
00491       if (internal::abs(m_matT.coeff(n,n-1)) > internal::abs(m_matT.coeff(n-1,n)))
00492       {
00493         m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
00494         m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
00495       }
00496       else
00497       {
00498         std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
00499         m_matT.coeffRef(n-1,n-1) = internal::real(cc);
00500         m_matT.coeffRef(n-1,n) = internal::imag(cc);
00501       }
00502       m_matT.coeffRef(n,n-1) = 0.0;
00503       m_matT.coeffRef(n,n) = 1.0;
00504       for (Index i = n-2; i >= 0; i--)
00505       {
00506         Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
00507         Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
00508         Scalar w = m_matT.coeff(i,i) - p;
00509 
00510         if (m_eivalues.coeff(i).imag() < 0.0)
00511         {
00512           lastw = w;
00513           lastra = ra;
00514           lastsa = sa;
00515         }
00516         else
00517         {
00518           l = i;
00519           if (m_eivalues.coeff(i).imag() == RealScalar(0))
00520           {
00521             std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
00522             m_matT.coeffRef(i,n-1) = internal::real(cc);
00523             m_matT.coeffRef(i,n) = internal::imag(cc);
00524           }
00525           else
00526           {
00527             // Solve complex equations
00528             Scalar x = m_matT.coeff(i,i+1);
00529             Scalar y = m_matT.coeff(i+1,i);
00530             Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
00531             Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
00532             if ((vr == 0.0) && (vi == 0.0))
00533               vr = eps * norm * (internal::abs(w) + internal::abs(q) + internal::abs(x) + internal::abs(y) + internal::abs(lastw));
00534 
00535             std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
00536             m_matT.coeffRef(i,n-1) = internal::real(cc);
00537             m_matT.coeffRef(i,n) = internal::imag(cc);
00538             if (internal::abs(x) > (internal::abs(lastw) + internal::abs(q)))
00539             {
00540               m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
00541               m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
00542             }
00543             else
00544             {
00545               cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q);
00546               m_matT.coeffRef(i+1,n-1) = internal::real(cc);
00547               m_matT.coeffRef(i+1,n) = internal::imag(cc);
00548             }
00549           }
00550 
00551           // Overflow control
00552           using std::max;
00553           Scalar t = (max)(internal::abs(m_matT.coeff(i,n-1)),internal::abs(m_matT.coeff(i,n)));
00554           if ((eps * t) * t > Scalar(1))
00555             m_matT.block(i, n-1, size-i, 2) /= t;
00556 
00557         }
00558       }
00559       
00560       // We handled a pair of complex conjugate eigenvalues, so need to skip them both
00561       n--;
00562     }
00563     else
00564     {
00565       eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen
00566     }
00567   }
00568 
00569   // Back transformation to get eigenvectors of original matrix
00570   for (Index j = size-1; j >= 0; j--)
00571   {
00572     m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
00573     m_eivec.col(j) = m_tmp;
00574   }
00575 }
00576 
00577 } // end namespace Eigen
00578 
00579 #endif // EIGEN_EIGENSOLVER_H


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