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 Eigen::Index Index; ///< \deprecated since Eigen 3.3
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    explicit 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    template<typename InputType>
126    explicit ComplexEigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true)
127            : m_eivec(matrix.rows(),matrix.cols()),
128              m_eivalues(matrix.cols()),
129              m_schur(matrix.rows()),
130              m_isInitialized(false),
131              m_eigenvectorsOk(false),
132              m_matX(matrix.rows(),matrix.cols())
133    {
134      compute(matrix.derived(), computeEigenvectors);
135    }
136
137    /** \brief Returns the eigenvectors of given matrix.
138      *
139      * \returns  A const reference to the matrix whose columns are the eigenvectors.
140      *
141      * \pre Either the constructor
142      * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
143      * function compute(const MatrixType& matrix, bool) has been called before
144      * to compute the eigendecomposition of a matrix, and
145      * \p computeEigenvectors was set to true (the default).
146      *
147      * This function returns a matrix whose columns are the eigenvectors. Column
148      * \f$ k \f$ is an eigenvector corresponding to eigenvalue number \f$ k
149      * \f$ as returned by eigenvalues().  The eigenvectors are normalized to
150      * have (Euclidean) norm equal to one. The matrix returned by this
151      * function is the matrix \f$ V \f$ in the eigendecomposition \f$ A = V D
152      * V^{-1} \f$, if it exists.
153      *
154      * Example: \include ComplexEigenSolver_eigenvectors.cpp
155      * Output: \verbinclude ComplexEigenSolver_eigenvectors.out
156      */
157    const EigenvectorType& eigenvectors() const
158    {
159      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
160      eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
161      return m_eivec;
162    }
163
164    /** \brief Returns the eigenvalues of given matrix.
165      *
166      * \returns A const reference to the column vector containing the eigenvalues.
167      *
168      * \pre Either the constructor
169      * ComplexEigenSolver(const MatrixType& matrix, bool) or the member
170      * function compute(const MatrixType& matrix, bool) has been called before
171      * to compute the eigendecomposition of a matrix.
172      *
173      * This function returns a column vector containing the
174      * eigenvalues. Eigenvalues are repeated according to their
175      * algebraic multiplicity, so there are as many eigenvalues as
176      * rows in the matrix. The eigenvalues are not sorted in any particular
177      * order.
178      *
179      * Example: \include ComplexEigenSolver_eigenvalues.cpp
180      * Output: \verbinclude ComplexEigenSolver_eigenvalues.out
181      */
182    const EigenvalueType& eigenvalues() const
183    {
184      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
185      return m_eivalues;
186    }
187
188    /** \brief Computes eigendecomposition of given matrix.
189      *
190      * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
191      * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
192      *    eigenvalues are computed; if false, only the eigenvalues are
193      *    computed.
194      * \returns    Reference to \c *this
195      *
196      * This function computes the eigenvalues of the complex matrix \p matrix.
197      * The eigenvalues() function can be used to retrieve them.  If
198      * \p computeEigenvectors is true, then the eigenvectors are also computed
199      * and can be retrieved by calling eigenvectors().
200      *
201      * The matrix is first reduced to Schur form using the
202      * ComplexSchur class. The Schur decomposition is then used to
203      * compute the eigenvalues and eigenvectors.
204      *
205      * The cost of the computation is dominated by the cost of the
206      * Schur decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$
207      * is the size of the matrix.
208      *
209      * Example: \include ComplexEigenSolver_compute.cpp
210      * Output: \verbinclude ComplexEigenSolver_compute.out
211      */
212    template<typename InputType>
213    ComplexEigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true);
214
215    /** \brief Reports whether previous computation was successful.
216      *
217      * \returns \c Success if computation was succesful, \c NoConvergence otherwise.
218      */
219    ComputationInfo info() const
220    {
221      eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
222      return m_schur.info();
223    }
224
225    /** \brief Sets the maximum number of iterations allowed. */
226    ComplexEigenSolver& setMaxIterations(Index maxIters)
227    {
228      m_schur.setMaxIterations(maxIters);
229      return *this;
230    }
231
232    /** \brief Returns the maximum number of iterations. */
233    Index getMaxIterations()
234    {
235      return m_schur.getMaxIterations();
236    }
237
238  protected:
239
240    static void check_template_parameters()
241    {
242      EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
243    }
244
245    EigenvectorType m_eivec;
246    EigenvalueType m_eivalues;
247    ComplexSchur<MatrixType> m_schur;
248    bool m_isInitialized;
249    bool m_eigenvectorsOk;
250    EigenvectorType m_matX;
251
252  private:
253    void doComputeEigenvectors(RealScalar matrixnorm);
254    void sortEigenvalues(bool computeEigenvectors);
255};
256
257
258template<typename MatrixType>
259template<typename InputType>
260ComplexEigenSolver<MatrixType>&
261ComplexEigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeEigenvectors)
262{
263  check_template_parameters();
264
265  // this code is inspired from Jampack
266  eigen_assert(matrix.cols() == matrix.rows());
267
268  // Do a complex Schur decomposition, A = U T U^*
269  // The eigenvalues are on the diagonal of T.
270  m_schur.compute(matrix.derived(), computeEigenvectors);
271
272  if(m_schur.info() == Success)
273  {
274    m_eivalues = m_schur.matrixT().diagonal();
275    if(computeEigenvectors)
276      doComputeEigenvectors(m_schur.matrixT().norm());
277    sortEigenvalues(computeEigenvectors);
278  }
279
280  m_isInitialized = true;
281  m_eigenvectorsOk = computeEigenvectors;
282  return *this;
283}
284
285
286template<typename MatrixType>
287void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm)
288{
289  const Index n = m_eivalues.size();
290
291  matrixnorm = numext::maxi(matrixnorm,(std::numeric_limits<RealScalar>::min)());
292
293  // Compute X such that T = X D X^(-1), where D is the diagonal of T.
294  // The matrix X is unit triangular.
295  m_matX = EigenvectorType::Zero(n, n);
296  for(Index k=n-1 ; k>=0 ; k--)
297  {
298    m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
299    // Compute X(i,k) using the (i,k) entry of the equation X T = D X
300    for(Index i=k-1 ; i>=0 ; i--)
301    {
302      m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
303      if(k-i-1>0)
304        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();
305      ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
306      if(z==ComplexScalar(0))
307      {
308        // If the i-th and k-th eigenvalue are equal, then z equals 0.
309        // Use a small value instead, to prevent division by zero.
310        numext::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
311      }
312      m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
313    }
314  }
315
316  // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
317  m_eivec.noalias() = m_schur.matrixU() * m_matX;
318  // .. and normalize the eigenvectors
319  for(Index k=0 ; k<n ; k++)
320  {
321    m_eivec.col(k).normalize();
322  }
323}
324
325
326template<typename MatrixType>
327void ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors)
328{
329  const Index n =  m_eivalues.size();
330  for (Index i=0; i<n; i++)
331  {
332    Index k;
333    m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
334    if (k != 0)
335    {
336      k += i;
337      std::swap(m_eivalues[k],m_eivalues[i]);
338      if(computeEigenvectors)
339	m_eivec.col(i).swap(m_eivec.col(k));
340    }
341  }
342}
343
344} // end namespace Eigen
345
346#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H
347