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 // 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 
00223   protected:
00224     EigenvectorType m_eivec;
00225     EigenvalueType m_eivalues;
00226     ComplexSchur<MatrixType> m_schur;
00227     bool m_isInitialized;
00228     bool m_eigenvectorsOk;
00229     EigenvectorType m_matX;
00230 
00231   private:
00232     void doComputeEigenvectors(RealScalar matrixnorm);
00233     void sortEigenvalues(bool computeEigenvectors);
00234 };
00235 
00236 
00237 template<typename MatrixType>
00238 ComplexEigenSolver<MatrixType>& ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
00239 {
00240   // this code is inspired from Jampack
00241   assert(matrix.cols() == matrix.rows());
00242 
00243   // Do a complex Schur decomposition, A = U T U^*
00244   // The eigenvalues are on the diagonal of T.
00245   m_schur.compute(matrix, computeEigenvectors);
00246 
00247   if(m_schur.info() == Success)
00248   {
00249     m_eivalues = m_schur.matrixT().diagonal();
00250     if(computeEigenvectors)
00251       doComputeEigenvectors(matrix.norm());
00252     sortEigenvalues(computeEigenvectors);
00253   }
00254 
00255   m_isInitialized = true;
00256   m_eigenvectorsOk = computeEigenvectors;
00257   return *this;
00258 }
00259 
00260 
00261 template<typename MatrixType>
00262 void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm)
00263 {
00264   const Index n = m_eivalues.size();
00265 
00266   // Compute X such that T = X D X^(-1), where D is the diagonal of T.
00267   // The matrix X is unit triangular.
00268   m_matX = EigenvectorType::Zero(n, n);
00269   for(Index k=n-1 ; k>=0 ; k--)
00270   {
00271     m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
00272     // Compute X(i,k) using the (i,k) entry of the equation X T = D X
00273     for(Index i=k-1 ; i>=0 ; i--)
00274     {
00275       m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
00276       if(k-i-1>0)
00277         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();
00278       ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
00279       if(z==ComplexScalar(0))
00280       {
00281         // If the i-th and k-th eigenvalue are equal, then z equals 0.
00282         // Use a small value instead, to prevent division by zero.
00283         internal::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
00284       }
00285       m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
00286     }
00287   }
00288 
00289   // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
00290   m_eivec.noalias() = m_schur.matrixU() * m_matX;
00291   // .. and normalize the eigenvectors
00292   for(Index k=0 ; k<n ; k++)
00293   {
00294     m_eivec.col(k).normalize();
00295   }
00296 }
00297 
00298 
00299 template<typename MatrixType>
00300 void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors)
00301 {
00302   const Index n =  m_eivalues.size();
00303   for (Index i=0; i<n; i++)
00304   {
00305     Index k;
00306     m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
00307     if (k != 0)
00308     {
00309       k += i;
00310       std::swap(m_eivalues[k],m_eivalues[i]);
00311       if(computeEigenvectors)
00312         m_eivec.col(i).swap(m_eivec.col(k));
00313     }
00314   }
00315 }
00316 
00317 } // end namespace Eigen
00318 
00319 #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H


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