ComplexEigenSolver.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) 2009 Claire Maurice
00005 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
00007 //
00008 // Eigen is free software; you can redistribute it and/or
00009 // modify it under the terms of the GNU Lesser General Public
00010 // License as published by the Free Software Foundation; either
00011 // version 3 of the License, or (at your option) any later version.
00012 //
00013 // Alternatively, you can redistribute it and/or
00014 // modify it under the terms of the GNU General Public License as
00015 // published by the Free Software Foundation; either version 2 of
00016 // the License, or (at your option) any later version.
00017 //
00018 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00019 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00020 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00021 // GNU General Public License for more details.
00022 //
00023 // You should have received a copy of the GNU Lesser General Public
00024 // License and a copy of the GNU General Public License along with
00025 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 #ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
00028 #define EIGEN_COMPLEX_EIGEN_SOLVER_H
00029 
00030 #include "./ComplexSchur.h"
00031 
00032 namespace Eigen { 
00033 
00060 template<typename _MatrixType> class ComplexEigenSolver
00061 {
00062   public:
00063 
00065     typedef _MatrixType MatrixType;
00066 
00067     enum {
00068       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00069       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00070       Options = MatrixType::Options,
00071       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00072       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00073     };
00074 
00076     typedef typename MatrixType::Scalar Scalar;
00077     typedef typename NumTraits<Scalar>::Real RealScalar;
00078     typedef typename MatrixType::Index Index;
00079 
00086     typedef std::complex<RealScalar> ComplexScalar;
00087 
00093     typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options&(~RowMajor), MaxColsAtCompileTime, 1> EigenvalueType;
00094 
00100     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorType;
00101 
00107     ComplexEigenSolver()
00108             : m_eivec(),
00109               m_eivalues(),
00110               m_schur(),
00111               m_isInitialized(false),
00112               m_eigenvectorsOk(false),
00113               m_matX()
00114     {}
00115 
00122     ComplexEigenSolver(Index size)
00123             : m_eivec(size, size),
00124               m_eivalues(size),
00125               m_schur(size),
00126               m_isInitialized(false),
00127               m_eigenvectorsOk(false),
00128               m_matX(size, size)
00129     {}
00130 
00140       ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
00141             : m_eivec(matrix.rows(),matrix.cols()),
00142               m_eivalues(matrix.cols()),
00143               m_schur(matrix.rows()),
00144               m_isInitialized(false),
00145               m_eigenvectorsOk(false),
00146               m_matX(matrix.rows(),matrix.cols())
00147     {
00148       compute(matrix, computeEigenvectors);
00149     }
00150 
00171     const EigenvectorType& eigenvectors() const
00172     {
00173       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00174       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00175       return m_eivec;
00176     }
00177 
00196     const EigenvalueType& eigenvalues() const
00197     {
00198       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00199       return m_eivalues;
00200     }
00201 
00226     ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
00227 
00232     ComputationInfo info() const
00233     {
00234       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00235       return m_schur.info();
00236     }
00237 
00238   protected:
00239     EigenvectorType m_eivec;
00240     EigenvalueType m_eivalues;
00241     ComplexSchur<MatrixType> m_schur;
00242     bool m_isInitialized;
00243     bool m_eigenvectorsOk;
00244     EigenvectorType m_matX;
00245 
00246   private:
00247     void doComputeEigenvectors(RealScalar matrixnorm);
00248     void sortEigenvalues(bool computeEigenvectors);
00249 };
00250 
00251 
00252 template<typename MatrixType>
00253 ComplexEigenSolver<MatrixType>& ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
00254 {
00255   // this code is inspired from Jampack
00256   assert(matrix.cols() == matrix.rows());
00257 
00258   // Do a complex Schur decomposition, A = U T U^*
00259   // The eigenvalues are on the diagonal of T.
00260   m_schur.compute(matrix, computeEigenvectors);
00261 
00262   if(m_schur.info() == Success)
00263   {
00264     m_eivalues = m_schur.matrixT().diagonal();
00265     if(computeEigenvectors)
00266       doComputeEigenvectors(matrix.norm());
00267     sortEigenvalues(computeEigenvectors);
00268   }
00269 
00270   m_isInitialized = true;
00271   m_eigenvectorsOk = computeEigenvectors;
00272   return *this;
00273 }
00274 
00275 
00276 template<typename MatrixType>
00277 void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm)
00278 {
00279   const Index n = m_eivalues.size();
00280 
00281   // Compute X such that T = X D X^(-1), where D is the diagonal of T.
00282   // The matrix X is unit triangular.
00283   m_matX = EigenvectorType::Zero(n, n);
00284   for(Index k=n-1 ; k>=0 ; k--)
00285   {
00286     m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
00287     // Compute X(i,k) using the (i,k) entry of the equation X T = D X
00288     for(Index i=k-1 ; i>=0 ; i--)
00289     {
00290       m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
00291       if(k-i-1>0)
00292         m_matX.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * m_matX.col(k).segment(i+1,k-i-1)).value();
00293       ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
00294       if(z==ComplexScalar(0))
00295       {
00296         // If the i-th and k-th eigenvalue are equal, then z equals 0.
00297         // Use a small value instead, to prevent division by zero.
00298         internal::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
00299       }
00300       m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
00301     }
00302   }
00303 
00304   // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
00305   m_eivec.noalias() = m_schur.matrixU() * m_matX;
00306   // .. and normalize the eigenvectors
00307   for(Index k=0 ; k<n ; k++)
00308   {
00309     m_eivec.col(k).normalize();
00310   }
00311 }
00312 
00313 
00314 template<typename MatrixType>
00315 void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors)
00316 {
00317   const Index n =  m_eivalues.size();
00318   for (Index i=0; i<n; i++)
00319   {
00320     Index k;
00321     m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
00322     if (k != 0)
00323     {
00324       k += i;
00325       std::swap(m_eivalues[k],m_eivalues[i]);
00326       if(computeEigenvectors)
00327         m_eivec.col(i).swap(m_eivec.col(k));
00328     }
00329   }
00330 }
00331 
00332 } // end namespace Eigen
00333 
00334 #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H