1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2009 Claire Maurice
5// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
6// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
7//
8// This Source Code Form is subject to the terms of the Mozilla
9// Public License v. 2.0. If a copy of the MPL was not distributed
10// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
11
12#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
13#define EIGEN_COMPLEX_EIGEN_SOLVER_H
14
15#include "./ComplexSchur.h"
16
17namespace Eigen {
18
19/** \eigenvalues_module \ingroup Eigenvalues_Module
20  *
21  *
22  * \class ComplexEigenSolver
23  *
24  * \brief Computes eigenvalues and eigenvectors of general complex matrices
25  *
26  * \tparam _MatrixType the type of the matrix of which we are
27  * computing the eigendecomposition; this is expected to be an
28  * instantiation of the Matrix class template.
29  *
30  * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars
31  * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v
32  * \f$.  If \f$ D \f$ is a diagonal matrix with the eigenvalues on
33  * the diagonal, and \f$ V \f$ is a matrix with the eigenvectors as
34  * its columns, then \f$ A V = V D \f$. The matrix \f$ V \f$ is
35  * almost always invertible, in which case we have \f$ A = V D V^{-1}
36  * \f$. This is called the eigendecomposition.
37  *
38  * The main function in this class is compute(), which computes the
39  * eigenvalues and eigenvectors of a given function. The
40  * documentation for that function contains an example showing the
41  * main features of the class.
42  *
43  * \sa class EigenSolver, class SelfAdjointEigenSolver
44  */
45template<typename _MatrixType> class ComplexEigenSolver
46{
47  public:
48
49    /** \brief Synonym for the template parameter \p _MatrixType. */
50    typedef _MatrixType MatrixType;
51
52    enum {
53      RowsAtCompileTime = MatrixType::RowsAtCompileTime,
54      ColsAtCompileTime = MatrixType::ColsAtCompileTime,
55      Options = MatrixType::Options,
56      MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
57      MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
58    };
59
60    /** \brief Scalar type for matrices of type #MatrixType. */
61    typedef typename MatrixType::Scalar Scalar;
62    typedef typename NumTraits<Scalar>::Real RealScalar;
63    typedef typename MatrixType::Index Index;
64
65    /** \brief Complex scalar type for #MatrixType.
66      *
67      * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
68      * \c float or \c double) and just \c Scalar if #Scalar is
69      * complex.
70      */
71    typedef std::complex<RealScalar> ComplexScalar;
72
73    /** \brief Type for vector of eigenvalues as returned by eigenvalues().
74      *
75      * This is a column vector with entries of type #ComplexScalar.
76      * The length of the vector is the size of #MatrixType.
77      */
78    typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options&(~RowMajor), MaxColsAtCompileTime, 1> EigenvalueType;
79
80    /** \brief Type for matrix of eigenvectors as returned by eigenvectors().
81      *
82      * This is a square matrix with entries of type #ComplexScalar.
83      * The size is the same as the size of #MatrixType.
84      */
85    typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorType;
86
87    /** \brief Default constructor.
88      *
89      * The default constructor is useful in cases in which the user intends to
90      * perform decompositions via compute().
91      */
92    ComplexEigenSolver()
93            : m_eivec(),
94              m_eivalues(),
95              m_schur(),
96              m_isInitialized(false),
97              m_eigenvectorsOk(false),
98              m_matX()
99    {}
100
101    /** \brief Default Constructor with memory preallocation
102      *
103      * Like the default constructor but with preallocation of the internal data
104      * according to the specified problem \a size.
105      * \sa ComplexEigenSolver()
106      */
107    ComplexEigenSolver(Index size)
108            : m_eivec(size, size),
109              m_eivalues(size),
110              m_schur(size),
111              m_isInitialized(false),
112              m_eigenvectorsOk(false),
113              m_matX(size, size)
114    {}
115
116    /** \brief Constructor; computes eigendecomposition of given matrix.
117      *
118      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
119      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
120      *    eigenvalues are computed; if false, only the eigenvalues are
121      *    computed.
122      *
123      * This constructor calls compute() to compute the eigendecomposition.
124      */
125      ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
126            : m_eivec(matrix.rows(),matrix.cols()),
127              m_eivalues(matrix.cols()),
128              m_schur(matrix.rows()),
129              m_isInitialized(false),
130              m_eigenvectorsOk(false),
131              m_matX(matrix.rows(),matrix.cols())
132    {
133      compute(matrix, computeEigenvectors);
134    }
135
136    /** \brief Returns the eigenvectors of given matrix.
137      *
138      * \returns  A const reference to the matrix whose columns are the eigenvectors.
139      *
140      * \pre Either the constructor
141      * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
142      * function compute(const MatrixType& matrix, bool) has been called before
143      * to compute the eigendecomposition of a matrix, and
144      * \p computeEigenvectors was set to true (the default).
145      *
146      * This function returns a matrix whose columns are the eigenvectors. Column
147      * \f$ k \f$ is an eigenvector corresponding to eigenvalue number \f$ k
148      * \f$ as returned by eigenvalues().  The eigenvectors are normalized to
149      * have (Euclidean) norm equal to one. The matrix returned by this
150      * function is the matrix \f$ V \f$ in the eigendecomposition \f$ A = V D
151      * V^{-1} \f$, if it exists.
152      *
153      * Example: \include ComplexEigenSolver_eigenvectors.cpp
154      * Output: \verbinclude ComplexEigenSolver_eigenvectors.out
155      */
156    const EigenvectorType& eigenvectors() const
157    {
158      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
159      eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
160      return m_eivec;
161    }
162
163    /** \brief Returns the eigenvalues of given matrix.
164      *
165      * \returns A const reference to the column vector containing the eigenvalues.
166      *
167      * \pre Either the constructor
168      * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
169      * function compute(const MatrixType& matrix, bool) has been called before
170      * to compute the eigendecomposition of a matrix.
171      *
172      * This function returns a column vector containing the
173      * eigenvalues. Eigenvalues are repeated according to their
174      * algebraic multiplicity, so there are as many eigenvalues as
175      * rows in the matrix. The eigenvalues are not sorted in any particular
176      * order.
177      *
178      * Example: \include ComplexEigenSolver_eigenvalues.cpp
179      * Output: \verbinclude ComplexEigenSolver_eigenvalues.out
180      */
181    const EigenvalueType& eigenvalues() const
182    {
183      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
184      return m_eivalues;
185    }
186
187    /** \brief Computes eigendecomposition of given matrix.
188      *
189      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
190      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
191      *    eigenvalues are computed; if false, only the eigenvalues are
192      *    computed.
193      * \returns    Reference to \c *this
194      *
195      * This function computes the eigenvalues of the complex matrix \p matrix.
196      * The eigenvalues() function can be used to retrieve them.  If
197      * \p computeEigenvectors is true, then the eigenvectors are also computed
198      * and can be retrieved by calling eigenvectors().
199      *
200      * The matrix is first reduced to Schur form using the
201      * ComplexSchur class. The Schur decomposition is then used to
202      * compute the eigenvalues and eigenvectors.
203      *
204      * The cost of the computation is dominated by the cost of the
205      * Schur decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$
206      * is the size of the matrix.
207      *
208      * Example: \include ComplexEigenSolver_compute.cpp
209      * Output: \verbinclude ComplexEigenSolver_compute.out
210      */
211    ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
212
213    /** \brief Reports whether previous computation was successful.
214      *
215      * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
216      */
217    ComputationInfo info() const
218    {
219      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
220      return m_schur.info();
221    }
222
223    /** \brief Sets the maximum number of iterations allowed. */
224    ComplexEigenSolver& setMaxIterations(Index maxIters)
225    {
226      m_schur.setMaxIterations(maxIters);
227      return *this;
228    }
229
230    /** \brief Returns the maximum number of iterations. */
231    Index getMaxIterations()
232    {
233      return m_schur.getMaxIterations();
234    }
235
236  protected:
237    EigenvectorType m_eivec;
238    EigenvalueType m_eivalues;
239    ComplexSchur<MatrixType> m_schur;
240    bool m_isInitialized;
241    bool m_eigenvectorsOk;
242    EigenvectorType m_matX;
243
244  private:
245    void doComputeEigenvectors(const RealScalar& matrixnorm);
246    void sortEigenvalues(bool computeEigenvectors);
247};
248
249
250template<typename MatrixType>
251ComplexEigenSolver<MatrixType>&
252ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
253{
254  // this code is inspired from Jampack
255  eigen_assert(matrix.cols() == matrix.rows());
256
257  // Do a complex Schur decomposition, A = U T U^*
258  // The eigenvalues are on the diagonal of T.
259  m_schur.compute(matrix, computeEigenvectors);
260
261  if(m_schur.info() == Success)
262  {
263    m_eivalues = m_schur.matrixT().diagonal();
264    if(computeEigenvectors)
265      doComputeEigenvectors(matrix.norm());
266    sortEigenvalues(computeEigenvectors);
267  }
268
269  m_isInitialized = true;
270  m_eigenvectorsOk = computeEigenvectors;
271  return *this;
272}
273
274
275template<typename MatrixType>
276void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(const RealScalar& matrixnorm)
277{
278  const Index n = m_eivalues.size();
279
280  // Compute X such that T = X D X^(-1), where D is the diagonal of T.
281  // The matrix X is unit triangular.
282  m_matX = EigenvectorType::Zero(n, n);
283  for(Index k=n-1 ; k>=0 ; k--)
284  {
285    m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
286    // Compute X(i,k) using the (i,k) entry of the equation X T = D X
287    for(Index i=k-1 ; i>=0 ; i--)
288    {
289      m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
290      if(k-i-1>0)
291        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();
292      ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
293      if(z==ComplexScalar(0))
294      {
295        // If the i-th and k-th eigenvalue are equal, then z equals 0.
296        // Use a small value instead, to prevent division by zero.
297        numext::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
298      }
299      m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
300    }
301  }
302
303  // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
304  m_eivec.noalias() = m_schur.matrixU() * m_matX;
305  // .. and normalize the eigenvectors
306  for(Index k=0 ; k<n ; k++)
307  {
308    m_eivec.col(k).normalize();
309  }
310}
311
312
313template<typename MatrixType>
314void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors)
315{
316  const Index n =  m_eivalues.size();
317  for (Index i=0; i<n; i++)
318  {
319    Index k;
320    m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
321    if (k != 0)
322    {
323      k += i;
324      std::swap(m_eivalues[k],m_eivalues[i]);
325      if(computeEigenvectors)
326	m_eivec.col(i).swap(m_eivec.col(k));
327    }
328  }
329}
330
331} // end namespace Eigen
332
333#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H
334