ComplexEigenSolver.h 11.9 KB
Newer Older
LM's avatar
LM committed
1 2 3 4 5
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Claire Maurice
// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
Don Gagne's avatar
Don Gagne committed
6
// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
LM's avatar
LM committed
7
//
Don Gagne's avatar
Don Gagne committed
8 9 10
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
LM's avatar
LM committed
11 12 13 14 15 16

#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
#define EIGEN_COMPLEX_EIGEN_SOLVER_H

#include "./ComplexSchur.h"

Don Gagne's avatar
Don Gagne committed
17 18
namespace Eigen { 

LM's avatar
LM committed
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222
/** \eigenvalues_module \ingroup Eigenvalues_Module
  *
  *
  * \class ComplexEigenSolver
  *
  * \brief Computes eigenvalues and eigenvectors of general complex matrices
  *
  * \tparam _MatrixType the type of the matrix of which we are
  * computing the eigendecomposition; this is expected to be an
  * instantiation of the Matrix class template.
  *
  * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars
  * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v
  * \f$.  If \f$ D \f$ is a diagonal matrix with the eigenvalues on
  * the diagonal, and \f$ V \f$ is a matrix with the eigenvectors as
  * its columns, then \f$ A V = V D \f$. The matrix \f$ V \f$ is
  * almost always invertible, in which case we have \f$ A = V D V^{-1}
  * \f$. This is called the eigendecomposition.
  *
  * The main function in this class is compute(), which computes the
  * eigenvalues and eigenvectors of a given function. The
  * documentation for that function contains an example showing the
  * main features of the class.
  *
  * \sa class EigenSolver, class SelfAdjointEigenSolver
  */
template<typename _MatrixType> class ComplexEigenSolver
{
  public:

    /** \brief Synonym for the template parameter \p _MatrixType. */
    typedef _MatrixType MatrixType;

    enum {
      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
      Options = MatrixType::Options,
      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
    };

    /** \brief Scalar type for matrices of type #MatrixType. */
    typedef typename MatrixType::Scalar Scalar;
    typedef typename NumTraits<Scalar>::Real RealScalar;
    typedef typename MatrixType::Index Index;

    /** \brief Complex scalar type for #MatrixType.
      *
      * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
      * \c float or \c double) and just \c Scalar if #Scalar is
      * complex.
      */
    typedef std::complex<RealScalar> ComplexScalar;

    /** \brief Type for vector of eigenvalues as returned by eigenvalues().
      *
      * This is a column vector with entries of type #ComplexScalar.
      * The length of the vector is the size of #MatrixType.
      */
    typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options&(~RowMajor), MaxColsAtCompileTime, 1> EigenvalueType;

    /** \brief Type for matrix of eigenvectors as returned by eigenvectors().
      *
      * This is a square matrix with entries of type #ComplexScalar.
      * The size is the same as the size of #MatrixType.
      */
    typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorType;

    /** \brief Default constructor.
      *
      * The default constructor is useful in cases in which the user intends to
      * perform decompositions via compute().
      */
    ComplexEigenSolver()
            : m_eivec(),
              m_eivalues(),
              m_schur(),
              m_isInitialized(false),
              m_eigenvectorsOk(false),
              m_matX()
    {}

    /** \brief Default Constructor with memory preallocation
      *
      * Like the default constructor but with preallocation of the internal data
      * according to the specified problem \a size.
      * \sa ComplexEigenSolver()
      */
    ComplexEigenSolver(Index size)
            : m_eivec(size, size),
              m_eivalues(size),
              m_schur(size),
              m_isInitialized(false),
              m_eigenvectorsOk(false),
              m_matX(size, size)
    {}

    /** \brief Constructor; computes eigendecomposition of given matrix.
      *
      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
      *    eigenvalues are computed; if false, only the eigenvalues are
      *    computed.
      *
      * This constructor calls compute() to compute the eigendecomposition.
      */
      ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
            : m_eivec(matrix.rows(),matrix.cols()),
              m_eivalues(matrix.cols()),
              m_schur(matrix.rows()),
              m_isInitialized(false),
              m_eigenvectorsOk(false),
              m_matX(matrix.rows(),matrix.cols())
    {
      compute(matrix, computeEigenvectors);
    }

    /** \brief Returns the eigenvectors of given matrix.
      *
      * \returns  A const reference to the matrix whose columns are the eigenvectors.
      *
      * \pre Either the constructor
      * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
      * function compute(const MatrixType& matrix, bool) has been called before
      * to compute the eigendecomposition of a matrix, and
      * \p computeEigenvectors was set to true (the default).
      *
      * This function returns a matrix whose columns are the eigenvectors. Column
      * \f$ k \f$ is an eigenvector corresponding to eigenvalue number \f$ k
      * \f$ as returned by eigenvalues().  The eigenvectors are normalized to
      * have (Euclidean) norm equal to one. The matrix returned by this
      * function is the matrix \f$ V \f$ in the eigendecomposition \f$ A = V D
      * V^{-1} \f$, if it exists.
      *
      * Example: \include ComplexEigenSolver_eigenvectors.cpp
      * Output: \verbinclude ComplexEigenSolver_eigenvectors.out
      */
    const EigenvectorType& eigenvectors() const
    {
      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
      eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
      return m_eivec;
    }

    /** \brief Returns the eigenvalues of given matrix.
      *
      * \returns A const reference to the column vector containing the eigenvalues.
      *
      * \pre Either the constructor
      * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
      * function compute(const MatrixType& matrix, bool) has been called before
      * to compute the eigendecomposition of a matrix.
      *
      * This function returns a column vector containing the
      * eigenvalues. Eigenvalues are repeated according to their
      * algebraic multiplicity, so there are as many eigenvalues as
      * rows in the matrix. The eigenvalues are not sorted in any particular
      * order.
      *
      * Example: \include ComplexEigenSolver_eigenvalues.cpp
      * Output: \verbinclude ComplexEigenSolver_eigenvalues.out
      */
    const EigenvalueType& eigenvalues() const
    {
      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
      return m_eivalues;
    }

    /** \brief Computes eigendecomposition of given matrix.
      *
      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
      *    eigenvalues are computed; if false, only the eigenvalues are
      *    computed.
      * \returns    Reference to \c *this
      *
      * This function computes the eigenvalues of the complex matrix \p matrix.
      * The eigenvalues() function can be used to retrieve them.  If
      * \p computeEigenvectors is true, then the eigenvectors are also computed
      * and can be retrieved by calling eigenvectors().
      *
      * The matrix is first reduced to Schur form using the
      * ComplexSchur class. The Schur decomposition is then used to
      * compute the eigenvalues and eigenvectors.
      *
      * The cost of the computation is dominated by the cost of the
      * Schur decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$
      * is the size of the matrix.
      *
      * Example: \include ComplexEigenSolver_compute.cpp
      * Output: \verbinclude ComplexEigenSolver_compute.out
      */
    ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);

    /** \brief Reports whether previous computation was successful.
      *
      * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
      */
    ComputationInfo info() const
    {
      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
      return m_schur.info();
    }

Don Gagne's avatar
Don Gagne committed
223 224 225 226 227 228 229 230 231 232 233 234 235
    /** \brief Sets the maximum number of iterations allowed. */
    ComplexEigenSolver& setMaxIterations(Index maxIters)
    {
      m_schur.setMaxIterations(maxIters);
      return *this;
    }

    /** \brief Returns the maximum number of iterations. */
    Index getMaxIterations()
    {
      return m_schur.getMaxIterations();
    }

LM's avatar
LM committed
236 237 238 239 240 241 242 243 244
  protected:
    EigenvectorType m_eivec;
    EigenvalueType m_eivalues;
    ComplexSchur<MatrixType> m_schur;
    bool m_isInitialized;
    bool m_eigenvectorsOk;
    EigenvectorType m_matX;

  private:
Don Gagne's avatar
Don Gagne committed
245
    void doComputeEigenvectors(const RealScalar& matrixnorm);
LM's avatar
LM committed
246 247 248 249 250
    void sortEigenvalues(bool computeEigenvectors);
};


template<typename MatrixType>
Don Gagne's avatar
Don Gagne committed
251 252
ComplexEigenSolver<MatrixType>& 
ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
LM's avatar
LM committed
253 254
{
  // this code is inspired from Jampack
Don Gagne's avatar
Don Gagne committed
255
  eigen_assert(matrix.cols() == matrix.rows());
LM's avatar
LM committed
256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275

  // Do a complex Schur decomposition, A = U T U^*
  // The eigenvalues are on the diagonal of T.
  m_schur.compute(matrix, computeEigenvectors);

  if(m_schur.info() == Success)
  {
    m_eivalues = m_schur.matrixT().diagonal();
    if(computeEigenvectors)
      doComputeEigenvectors(matrix.norm());
    sortEigenvalues(computeEigenvectors);
  }

  m_isInitialized = true;
  m_eigenvectorsOk = computeEigenvectors;
  return *this;
}


template<typename MatrixType>
Don Gagne's avatar
Don Gagne committed
276
void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(const RealScalar& matrixnorm)
LM's avatar
LM committed
277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296
{
  const Index n = m_eivalues.size();

  // Compute X such that T = X D X^(-1), where D is the diagonal of T.
  // The matrix X is unit triangular.
  m_matX = EigenvectorType::Zero(n, n);
  for(Index k=n-1 ; k>=0 ; k--)
  {
    m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
    // Compute X(i,k) using the (i,k) entry of the equation X T = D X
    for(Index i=k-1 ; i>=0 ; i--)
    {
      m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
      if(k-i-1>0)
        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();
      ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
      if(z==ComplexScalar(0))
      {
        // If the i-th and k-th eigenvalue are equal, then z equals 0.
        // Use a small value instead, to prevent division by zero.
Don Gagne's avatar
Don Gagne committed
297
        numext::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
LM's avatar
LM committed
298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330
      }
      m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
    }
  }

  // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
  m_eivec.noalias() = m_schur.matrixU() * m_matX;
  // .. and normalize the eigenvectors
  for(Index k=0 ; k<n ; k++)
  {
    m_eivec.col(k).normalize();
  }
}


template<typename MatrixType>
void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors)
{
  const Index n =  m_eivalues.size();
  for (Index i=0; i<n; i++)
  {
    Index k;
    m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
    if (k != 0)
    {
      k += i;
      std::swap(m_eivalues[k],m_eivalues[i]);
      if(computeEigenvectors)
	m_eivec.col(i).swap(m_eivec.col(k));
    }
  }
}

Don Gagne's avatar
Don Gagne committed
331
} // end namespace Eigen
LM's avatar
LM committed
332 333

#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H