Eigen  3.3.3
ComplexEigenSolver.h
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,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
00007 //
00008 // This Source Code Form is subject to the terms of the Mozilla
00009 // Public License v. 2.0. If a copy of the MPL was not distributed
00010 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00011 
00012 #ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
00013 #define EIGEN_COMPLEX_EIGEN_SOLVER_H
00014 
00015 #include "./ComplexSchur.h"
00016 
00017 namespace Eigen { 
00018 
00045 template<typename _MatrixType> class ComplexEigenSolver
00046 {
00047   public:
00048 
00050     typedef _MatrixType MatrixType;
00051 
00052     enum {
00053       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00054       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00055       Options = MatrixType::Options,
00056       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00057       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00058     };
00059 
00061     typedef typename MatrixType::Scalar Scalar;
00062     typedef typename NumTraits<Scalar>::Real RealScalar;
00063     typedef Eigen::Index Index; 
00064 
00071     typedef std::complex<RealScalar> ComplexScalar;
00072 
00078     typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options&(~RowMajor), MaxColsAtCompileTime, 1> EigenvalueType;
00079 
00085     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorType;
00086 
00092     ComplexEigenSolver()
00093             : m_eivec(),
00094               m_eivalues(),
00095               m_schur(),
00096               m_isInitialized(false),
00097               m_eigenvectorsOk(false),
00098               m_matX()
00099     {}
00100 
00107     explicit ComplexEigenSolver(Index size)
00108             : m_eivec(size, size),
00109               m_eivalues(size),
00110               m_schur(size),
00111               m_isInitialized(false),
00112               m_eigenvectorsOk(false),
00113               m_matX(size, size)
00114     {}
00115 
00125     template<typename InputType>
00126     explicit ComplexEigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true)
00127             : m_eivec(matrix.rows(),matrix.cols()),
00128               m_eivalues(matrix.cols()),
00129               m_schur(matrix.rows()),
00130               m_isInitialized(false),
00131               m_eigenvectorsOk(false),
00132               m_matX(matrix.rows(),matrix.cols())
00133     {
00134       compute(matrix.derived(), computeEigenvectors);
00135     }
00136 
00157     const EigenvectorType& eigenvectors() const
00158     {
00159       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00160       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00161       return m_eivec;
00162     }
00163 
00182     const EigenvalueType& eigenvalues() const
00183     {
00184       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00185       return m_eivalues;
00186     }
00187 
00212     template<typename InputType>
00213     ComplexEigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true);
00214 
00219     ComputationInfo info() const
00220     {
00221       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00222       return m_schur.info();
00223     }
00224 
00226     ComplexEigenSolver& setMaxIterations(Index maxIters)
00227     {
00228       m_schur.setMaxIterations(maxIters);
00229       return *this;
00230     }
00231 
00233     Index getMaxIterations()
00234     {
00235       return m_schur.getMaxIterations();
00236     }
00237 
00238   protected:
00239     
00240     static void check_template_parameters()
00241     {
00242       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00243     }
00244     
00245     EigenvectorType m_eivec;
00246     EigenvalueType m_eivalues;
00247     ComplexSchur<MatrixType> m_schur;
00248     bool m_isInitialized;
00249     bool m_eigenvectorsOk;
00250     EigenvectorType m_matX;
00251 
00252   private:
00253     void doComputeEigenvectors(RealScalar matrixnorm);
00254     void sortEigenvalues(bool computeEigenvectors);
00255 };
00256 
00257 
00258 template<typename MatrixType>
00259 template<typename InputType>
00260 ComplexEigenSolver<MatrixType>& 
00261 ComplexEigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeEigenvectors)
00262 {
00263   check_template_parameters();
00264   
00265   // this code is inspired from Jampack
00266   eigen_assert(matrix.cols() == matrix.rows());
00267 
00268   // Do a complex Schur decomposition, A = U T U^*
00269   // The eigenvalues are on the diagonal of T.
00270   m_schur.compute(matrix.derived(), computeEigenvectors);
00271 
00272   if(m_schur.info() == Success)
00273   {
00274     m_eivalues = m_schur.matrixT().diagonal();
00275     if(computeEigenvectors)
00276       doComputeEigenvectors(m_schur.matrixT().norm());
00277     sortEigenvalues(computeEigenvectors);
00278   }
00279 
00280   m_isInitialized = true;
00281   m_eigenvectorsOk = computeEigenvectors;
00282   return *this;
00283 }
00284 
00285 
00286 template<typename MatrixType>
00287 void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm)
00288 {
00289   const Index n = m_eivalues.size();
00290 
00291   matrixnorm = numext::maxi(matrixnorm,(std::numeric_limits<RealScalar>::min)());
00292 
00293   // Compute X such that T = X D X^(-1), where D is the diagonal of T.
00294   // The matrix X is unit triangular.
00295   m_matX = EigenvectorType::Zero(n, n);
00296   for(Index k=n-1 ; k>=0 ; k--)
00297   {
00298     m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
00299     // Compute X(i,k) using the (i,k) entry of the equation X T = D X
00300     for(Index i=k-1 ; i>=0 ; i--)
00301     {
00302       m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
00303       if(k-i-1>0)
00304         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();
00305       ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
00306       if(z==ComplexScalar(0))
00307       {
00308         // If the i-th and k-th eigenvalue are equal, then z equals 0.
00309         // Use a small value instead, to prevent division by zero.
00310         numext::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
00311       }
00312       m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
00313     }
00314   }
00315 
00316   // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
00317   m_eivec.noalias() = m_schur.matrixU() * m_matX;
00318   // .. and normalize the eigenvectors
00319   for(Index k=0 ; k<n ; k++)
00320   {
00321     m_eivec.col(k).normalize();
00322   }
00323 }
00324 
00325 
00326 template<typename MatrixType>
00327 void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors)
00328 {
00329   const Index n =  m_eivalues.size();
00330   for (Index i=0; i<n; i++)
00331   {
00332     Index k;
00333     m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
00334     if (k != 0)
00335     {
00336       k += i;
00337       std::swap(m_eivalues[k],m_eivalues[i]);
00338       if(computeEigenvectors)
00339         m_eivec.col(i).swap(m_eivec.col(k));
00340     }
00341   }
00342 }
00343 
00344 } // end namespace Eigen
00345 
00346 #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends