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 
00019 /** \eigenvalues_module \ingroup Eigenvalues_Module
00020   *
00021   *
00022   * \class ComplexEigenSolver
00023   *
00024   * \brief Computes eigenvalues and eigenvectors of general complex matrices
00025   *
00026   * \tparam _MatrixType the type of the matrix of which we are
00027   * computing the eigendecomposition; this is expected to be an
00028   * instantiation of the Matrix class template.
00029   *
00030   * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars
00031   * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v
00032   * \f$.  If \f$ D \f$ is a diagonal matrix with the eigenvalues on
00033   * the diagonal, and \f$ V \f$ is a matrix with the eigenvectors as
00034   * its columns, then \f$ A V = V D \f$. The matrix \f$ V \f$ is
00035   * almost always invertible, in which case we have \f$ A = V D V^{-1}
00036   * \f$. This is called the eigendecomposition.
00037   *
00038   * The main function in this class is compute(), which computes the
00039   * eigenvalues and eigenvectors of a given function. The
00040   * documentation for that function contains an example showing the
00041   * main features of the class.
00042   *
00043   * \sa class EigenSolver, class SelfAdjointEigenSolver
00044   */
00045 template<typename _MatrixType> class ComplexEigenSolver 
00046 {
00047   public:
00048 
00049     /** \brief Synonym for the template parameter \p _MatrixType. */
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 
00060     /** \brief Scalar type for matrices of type #MatrixType. */
00061     typedef typename MatrixType::Scalar Scalar;
00062     typedef typename NumTraits<Scalar>::Real RealScalar;
00063     typedef typename MatrixType::Index Index;
00064 
00065     /** \brief Complex scalar type for #MatrixType.
00066       *
00067       * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
00068       * \c float or \c double) and just \c Scalar if #Scalar is
00069       * complex.
00070       */
00071     typedef std::complex<RealScalar> ComplexScalar;
00072 
00073     /** \brief Type for vector of eigenvalues as returned by eigenvalues().
00074       *
00075       * This is a column vector with entries of type #ComplexScalar.
00076       * The length of the vector is the size of #MatrixType.
00077       */
00078     typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options&(~RowMajor), MaxColsAtCompileTime, 1>  EigenvalueType;
00079 
00080     /** \brief Type for matrix of eigenvectors as returned by eigenvectors().
00081       *
00082       * This is a square matrix with entries of type #ComplexScalar.
00083       * The size is the same as the size of #MatrixType.
00084       */
00085     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime>  EigenvectorType;
00086 
00087     /** \brief Default constructor.
00088       *
00089       * The default constructor is useful in cases in which the user intends to
00090       * perform decompositions via compute().
00091       */
00092     ComplexEigenSolver()
00093             : m_eivec(),
00094               m_eivalues(),
00095               m_schur(),
00096               m_isInitialized(false),
00097               m_eigenvectorsOk(false),
00098               m_matX()
00099     {}
00100 
00101     /** \brief Default Constructor with memory preallocation
00102       *
00103       * Like the default constructor but with preallocation of the internal data
00104       * according to the specified problem \a size.
00105       * \sa ComplexEigenSolver()
00106       */
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 
00116     /** \brief Constructor; computes eigendecomposition of given matrix.
00117       *
00118       * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
00119       * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
00120       *    eigenvalues are computed; if false, only the eigenvalues are
00121       *    computed.
00122       *
00123       * This constructor calls compute() to compute the eigendecomposition.
00124       */
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 
00136     /** \brief Returns the eigenvectors of given matrix.
00137       *
00138       * \returns  A const reference to the matrix whose columns are the eigenvectors.
00139       *
00140       * \pre Either the constructor
00141       * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
00142       * function compute(const MatrixType& matrix, bool) has been called before
00143       * to compute the eigendecomposition of a matrix, and
00144       * \p computeEigenvectors was set to true (the default).
00145       *
00146       * This function returns a matrix whose columns are the eigenvectors. Column
00147       * \f$ k \f$ is an eigenvector corresponding to eigenvalue number \f$ k
00148       * \f$ as returned by eigenvalues().  The eigenvectors are normalized to
00149       * have (Euclidean) norm equal to one. The matrix returned by this
00150       * function is the matrix \f$ V \f$ in the eigendecomposition \f$ A = V D
00151       * V^{-1} \f$, if it exists.
00152       *
00153       * Example: \include ComplexEigenSolver_eigenvectors.cpp
00154       * Output: \verbinclude ComplexEigenSolver_eigenvectors.out
00155       */
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 
00163     /** \brief Returns the eigenvalues of given matrix.
00164       *
00165       * \returns A const reference to the column vector containing the eigenvalues.
00166       *
00167       * \pre Either the constructor
00168       * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
00169       * function compute(const MatrixType& matrix, bool) has been called before
00170       * to compute the eigendecomposition of a matrix.
00171       *
00172       * This function returns a column vector containing the
00173       * eigenvalues. Eigenvalues are repeated according to their
00174       * algebraic multiplicity, so there are as many eigenvalues as
00175       * rows in the matrix. The eigenvalues are not sorted in any particular
00176       * order.
00177       *
00178       * Example: \include ComplexEigenSolver_eigenvalues.cpp
00179       * Output: \verbinclude ComplexEigenSolver_eigenvalues.out
00180       */
00181     const EigenvalueType & eigenvalues() const
00182     {
00183       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00184       return m_eivalues;
00185     }
00186 
00187     /** \brief Computes eigendecomposition of given matrix.
00188       *
00189       * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
00190       * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
00191       *    eigenvalues are computed; if false, only the eigenvalues are
00192       *    computed.
00193       * \returns    Reference to \c *this
00194       *
00195       * This function computes the eigenvalues of the complex matrix \p matrix.
00196       * The eigenvalues() function can be used to retrieve them.  If
00197       * \p computeEigenvectors is true, then the eigenvectors are also computed
00198       * and can be retrieved by calling eigenvectors().
00199       *
00200       * The matrix is first reduced to Schur form using the
00201       * ComplexSchur class. The Schur decomposition is then used to
00202       * compute the eigenvalues and eigenvectors.
00203       *
00204       * The cost of the computation is dominated by the cost of the
00205       * Schur decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$
00206       * is the size of the matrix.
00207       *
00208       * Example: \include ComplexEigenSolver_compute.cpp
00209       * Output: \verbinclude ComplexEigenSolver_compute.out
00210       */
00211     ComplexEigenSolver & compute(const MatrixType& matrix, bool computeEigenvectors = true);
00212 
00213     /** \brief Reports whether previous computation was successful.
00214       *
00215       * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
00216       */
00217     ComputationInfo info() const
00218     {
00219       eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
00220       return m_schur.info();
00221     }
00222 
00223     /** \brief Sets the maximum number of iterations allowed. */
00224     ComplexEigenSolver & setMaxIterations(Index maxIters)
00225     {
00226       m_schur.setMaxIterations(maxIters);
00227       return *this;
00228     }
00229 
00230     /** \brief Returns the maximum number of iterations. */
00231     Index getMaxIterations()
00232     {
00233       return m_schur.getMaxIterations();
00234     }
00235 
00236   protected:
00237     
00238     static void check_template_parameters()
00239     {
00240       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00241     }
00242     
00243     EigenvectorType m_eivec;
00244     EigenvalueType m_eivalues;
00245     ComplexSchur<MatrixType> m_schur;
00246     bool m_isInitialized;
00247     bool m_eigenvectorsOk;
00248     EigenvectorType m_matX;
00249 
00250   private:
00251     void doComputeEigenvectors(const RealScalar& matrixnorm);
00252     void sortEigenvalues(bool computeEigenvectors);
00253 };
00254 
00255 
00256 template<typename MatrixType>
00257 ComplexEigenSolver<MatrixType>& 
00258 ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
00259 {
00260   check_template_parameters();
00261   
00262   // this code is inspired from Jampack
00263   eigen_assert(matrix.cols() == matrix.rows());
00264 
00265   // Do a complex Schur decomposition, A = U T U^*
00266   // The eigenvalues are on the diagonal of T.
00267   m_schur.compute(matrix, computeEigenvectors);
00268 
00269   if(m_schur.info() == Success)
00270   {
00271     m_eivalues = m_schur.matrixT().diagonal();
00272     if(computeEigenvectors)
00273       doComputeEigenvectors(matrix.norm());
00274     sortEigenvalues(computeEigenvectors);
00275   }
00276 
00277   m_isInitialized = true;
00278   m_eigenvectorsOk = computeEigenvectors;
00279   return *this;
00280 }
00281 
00282 
00283 template<typename MatrixType>
00284 void ComplexEigenSolver<MatrixType>::doComputeEigenvectors (const RealScalar& matrixnorm)
00285 {
00286   const Index n = m_eivalues.size();
00287 
00288   // Compute X such that T = X D X^(-1), where D is the diagonal of T.
00289   // The matrix X is unit triangular.
00290   m_matX = EigenvectorType::Zero(n, n);
00291   for(Index k=n-1 ; k>=0 ; k--)
00292   {
00293     m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
00294     // Compute X(i,k) using the (i,k) entry of the equation X T = D X
00295     for(Index i=k-1 ; i>=0 ; i--)
00296     {
00297       m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
00298       if(k-i-1>0)
00299         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();
00300       ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
00301       if(z==ComplexScalar(0))
00302       {
00303         // If the i-th and k-th eigenvalue are equal, then z equals 0.
00304         // Use a small value instead, to prevent division by zero.
00305         numext::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
00306       }
00307       m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
00308     }
00309   }
00310 
00311   // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
00312   m_eivec.noalias() = m_schur.matrixU() * m_matX;
00313   // .. and normalize the eigenvectors
00314   for(Index k=0 ; k<n ; k++)
00315   {
00316     m_eivec.col(k).normalize();
00317   }
00318 }
00319 
00320 
00321 template<typename MatrixType>
00322 void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors)
00323 {
00324   const Index n =  m_eivalues.size();
00325   for (Index i=0; i<n; i++)
00326   {
00327     Index k;
00328     m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
00329     if (k != 0)
00330     {
00331       k += i;
00332       std::swap(m_eivalues[k],m_eivalues[i]);
00333       if(computeEigenvectors)
00334     m_eivec.col(i).swap(m_eivec.col(k));
00335     }
00336   }
00337 }
00338 
00339 } // end namespace Eigen
00340 
00341 #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H