• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
17 namespace 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   */
45 template<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       */
ComplexEigenSolver()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       */
ComplexEigenSolver(Index size)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       */
eigenvectors()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       */
eigenvalues()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       */
info()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. */
setMaxIterations(Index maxIters)226     ComplexEigenSolver& setMaxIterations(Index maxIters)
227     {
228       m_schur.setMaxIterations(maxIters);
229       return *this;
230     }
231 
232     /** \brief Returns the maximum number of iterations. */
getMaxIterations()233     Index getMaxIterations()
234     {
235       return m_schur.getMaxIterations();
236     }
237 
238   protected:
239 
check_template_parameters()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 
258 template<typename MatrixType>
259 template<typename InputType>
260 ComplexEigenSolver<MatrixType>&
compute(const EigenBase<InputType> & matrix,bool computeEigenvectors)261 ComplexEigenSolver<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 
286 template<typename MatrixType>
doComputeEigenvectors(RealScalar matrixnorm)287 void 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 
326 template<typename MatrixType>
sortEigenvalues(bool computeEigenvectors)327 void 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