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,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 typename MatrixType::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     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       ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
00126             : m_eivec(matrix.rows(),matrix.cols()),
00127               m_eivalues(matrix.cols()),
00128               m_schur(matrix.rows()),
00129               m_isInitialized(false),
00130               m_eigenvectorsOk(false),
00131               m_matX(matrix.rows(),matrix.cols())
00132     {
00133       compute(matrix, computeEigenvectors);
00134     }
00135 
00156     const EigenvectorType& eigenvectors() const
00157     {
00158       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00159       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00160       return m_eivec;
00161     }
00162 
00181     const EigenvalueType& eigenvalues() const
00182     {
00183       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00184       return m_eivalues;
00185     }
00186 
00211     ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
00212 
00217     ComputationInfo info() const
00218     {
00219       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00220       return m_schur.info();
00221     }
00222 
00224     ComplexEigenSolver& setMaxIterations(Index maxIters)
00225     {
00226       m_schur.setMaxIterations(maxIters);
00227       return *this;
00228     }
00229 
00231     Index getMaxIterations()
00232     {
00233       return m_schur.getMaxIterations();
00234     }
00235 
00236   protected:
00237     EigenvectorType m_eivec;
00238     EigenvalueType m_eivalues;
00239     ComplexSchur<MatrixType> m_schur;
00240     bool m_isInitialized;
00241     bool m_eigenvectorsOk;
00242     EigenvectorType m_matX;
00243 
00244   private:
00245     void doComputeEigenvectors(const RealScalar& matrixnorm);
00246     void sortEigenvalues(bool computeEigenvectors);
00247 };
00248 
00249 
00250 template<typename MatrixType>
00251 ComplexEigenSolver<MatrixType>& 
00252 ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
00253 {
00254   // this code is inspired from Jampack
00255   eigen_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(const 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         numext::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 } // end namespace Eigen
00332 
00333 #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:36:54