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 "./EigenvaluesCommon.h"
00031 #include "./ComplexSchur.h"
00032 
00059 template<typename _MatrixType> class ComplexEigenSolver
00060 {
00061   public:
00062 
00064     typedef _MatrixType MatrixType;
00065 
00066     enum {
00067       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00068       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00069       Options = MatrixType::Options,
00070       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00071       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00072     };
00073 
00075     typedef typename MatrixType::Scalar Scalar;
00076     typedef typename NumTraits<Scalar>::Real RealScalar;
00077     typedef typename MatrixType::Index Index;
00078 
00085     typedef std::complex<RealScalar> ComplexScalar;
00086 
00092     typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options&(~RowMajor), MaxColsAtCompileTime, 1> EigenvalueType;
00093 
00099     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorType;
00100 
00106     ComplexEigenSolver()
00107             : m_eivec(),
00108               m_eivalues(),
00109               m_schur(),
00110               m_isInitialized(false),
00111               m_eigenvectorsOk(false),
00112               m_matX()
00113     {}
00114 
00121     ComplexEigenSolver(Index size)
00122             : m_eivec(size, size),
00123               m_eivalues(size),
00124               m_schur(size),
00125               m_isInitialized(false),
00126               m_eigenvectorsOk(false),
00127               m_matX(size, size)
00128     {}
00129 
00139       ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
00140             : m_eivec(matrix.rows(),matrix.cols()),
00141               m_eivalues(matrix.cols()),
00142               m_schur(matrix.rows()),
00143               m_isInitialized(false),
00144               m_eigenvectorsOk(false),
00145               m_matX(matrix.rows(),matrix.cols())
00146     {
00147       compute(matrix, computeEigenvectors);
00148     }
00149 
00170     const EigenvectorType& eigenvectors() const
00171     {
00172       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00173       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00174       return m_eivec;
00175     }
00176 
00195     const EigenvalueType& eigenvalues() const
00196     {
00197       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00198       return m_eivalues;
00199     }
00200 
00225     ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
00226 
00231     ComputationInfo info() const
00232     {
00233       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00234       return m_schur.info();
00235     }
00236 
00237   protected:
00238     EigenvectorType m_eivec;
00239     EigenvalueType m_eivalues;
00240     ComplexSchur<MatrixType> m_schur;
00241     bool m_isInitialized;
00242     bool m_eigenvectorsOk;
00243     EigenvectorType m_matX;
00244 
00245   private:
00246     void doComputeEigenvectors(RealScalar matrixnorm);
00247     void sortEigenvalues(bool computeEigenvectors);
00248 };
00249 
00250 
00251 template<typename MatrixType>
00252 ComplexEigenSolver<MatrixType>& ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
00253 {
00254   // this code is inspired from Jampack
00255   assert(matrix.cols() == matrix.rows());
00256 
00257   // Do a complex Schur decomposition, A = U T U^*
00258   // The eigenvalues are on the diagonal of T.
00259   m_schur.compute(matrix, computeEigenvectors);
00260 
00261   if(m_schur.info() == Success)
00262   {
00263     m_eivalues = m_schur.matrixT().diagonal();
00264     if(computeEigenvectors)
00265       doComputeEigenvectors(matrix.norm());
00266     sortEigenvalues(computeEigenvectors);
00267   }
00268 
00269   m_isInitialized = true;
00270   m_eigenvectorsOk = computeEigenvectors;
00271   return *this;
00272 }
00273 
00274 
00275 template<typename MatrixType>
00276 void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm)
00277 {
00278   const Index n = m_eivalues.size();
00279 
00280   // Compute X such that T = X D X^(-1), where D is the diagonal of T.
00281   // The matrix X is unit triangular.
00282   m_matX = EigenvectorType::Zero(n, n);
00283   for(Index k=n-1 ; k>=0 ; k--)
00284   {
00285     m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
00286     // Compute X(i,k) using the (i,k) entry of the equation X T = D X
00287     for(Index i=k-1 ; i>=0 ; i--)
00288     {
00289       m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
00290       if(k-i-1>0)
00291         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();
00292       ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
00293       if(z==ComplexScalar(0))
00294       {
00295         // If the i-th and k-th eigenvalue are equal, then z equals 0.
00296         // Use a small value instead, to prevent division by zero.
00297         internal::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
00298       }
00299       m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
00300     }
00301   }
00302 
00303   // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
00304   m_eivec.noalias() = m_schur.matrixU() * m_matX;
00305   // .. and normalize the eigenvectors
00306   for(Index k=0 ; k<n ; k++)
00307   {
00308     m_eivec.col(k).normalize();
00309   }
00310 }
00311 
00312 
00313 template<typename MatrixType>
00314 void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors)
00315 {
00316   const Index n =  m_eivalues.size();
00317   for (Index i=0; i<n; i++)
00318   {
00319     Index k;
00320     m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
00321     if (k != 0)
00322     {
00323       k += i;
00324       std::swap(m_eivalues[k],m_eivalues[i]);
00325       if(computeEigenvectors)
00326         m_eivec.col(i).swap(m_eivec.col(k));
00327     }
00328   }
00329 }
00330 
00331 
00332 #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H


re_vision
Author(s): Dorian Galvez-Lopez
autogenerated on Sun Jan 5 2014 11:30:58