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 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 */
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 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 */
eigenvectors()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 */
eigenvalues()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 */
info()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. */
setMaxIterations(Index maxIters)224 ComplexEigenSolver& setMaxIterations(Index maxIters)
225 {
226 m_schur.setMaxIterations(maxIters);
227 return *this;
228 }
229
230 /** \brief Returns the maximum number of iterations. */
getMaxIterations()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
250 template<typename MatrixType>
251 ComplexEigenSolver<MatrixType>&
compute(const MatrixType & matrix,bool computeEigenvectors)252 ComplexEigenSolver<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
275 template<typename MatrixType>
doComputeEigenvectors(const RealScalar & matrixnorm)276 void 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
313 template<typename MatrixType>
sortEigenvalues(bool computeEigenvectors)314 void 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