1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
2 | // for linear algebra.
|
---|
3 | //
|
---|
4 | // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
|
---|
5 | // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
|
---|
6 | //
|
---|
7 | // This Source Code Form is subject to the terms of the Mozilla
|
---|
8 | // Public License v. 2.0. If a copy of the MPL was not distributed
|
---|
9 | // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
---|
10 |
|
---|
11 | #ifndef EIGEN_EIGENSOLVER_H
|
---|
12 | #define EIGEN_EIGENSOLVER_H
|
---|
13 |
|
---|
14 | #include "./RealSchur.h"
|
---|
15 |
|
---|
16 | namespace Eigen {
|
---|
17 |
|
---|
18 | /** \eigenvalues_module \ingroup Eigenvalues_Module
|
---|
19 | *
|
---|
20 | *
|
---|
21 | * \class EigenSolver
|
---|
22 | *
|
---|
23 | * \brief Computes eigenvalues and eigenvectors of general matrices
|
---|
24 | *
|
---|
25 | * \tparam _MatrixType the type of the matrix of which we are computing the
|
---|
26 | * eigendecomposition; this is expected to be an instantiation of the Matrix
|
---|
27 | * class template. Currently, only real matrices are supported.
|
---|
28 | *
|
---|
29 | * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars
|
---|
30 | * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v \f$. If
|
---|
31 | * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and
|
---|
32 | * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V =
|
---|
33 | * V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we
|
---|
34 | * have \f$ A = V D V^{-1} \f$. This is called the eigendecomposition.
|
---|
35 | *
|
---|
36 | * The eigenvalues and eigenvectors of a matrix may be complex, even when the
|
---|
37 | * matrix is real. However, we can choose real matrices \f$ V \f$ and \f$ D
|
---|
38 | * \f$ satisfying \f$ A V = V D \f$, just like the eigendecomposition, if the
|
---|
39 | * matrix \f$ D \f$ is not required to be diagonal, but if it is allowed to
|
---|
40 | * have blocks of the form
|
---|
41 | * \f[ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f]
|
---|
42 | * (where \f$ u \f$ and \f$ v \f$ are real numbers) on the diagonal. These
|
---|
43 | * blocks correspond to complex eigenvalue pairs \f$ u \pm iv \f$. We call
|
---|
44 | * this variant of the eigendecomposition the pseudo-eigendecomposition.
|
---|
45 | *
|
---|
46 | * Call the function compute() to compute the eigenvalues and eigenvectors of
|
---|
47 | * a given matrix. Alternatively, you can use the
|
---|
48 | * EigenSolver(const MatrixType&, bool) constructor which computes the
|
---|
49 | * eigenvalues and eigenvectors at construction time. Once the eigenvalue and
|
---|
50 | * eigenvectors are computed, they can be retrieved with the eigenvalues() and
|
---|
51 | * eigenvectors() functions. The pseudoEigenvalueMatrix() and
|
---|
52 | * pseudoEigenvectors() methods allow the construction of the
|
---|
53 | * pseudo-eigendecomposition.
|
---|
54 | *
|
---|
55 | * The documentation for EigenSolver(const MatrixType&, bool) contains an
|
---|
56 | * example of the typical use of this class.
|
---|
57 | *
|
---|
58 | * \note The implementation is adapted from
|
---|
59 | * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
|
---|
60 | * Their code is based on EISPACK.
|
---|
61 | *
|
---|
62 | * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver
|
---|
63 | */
|
---|
64 | template<typename _MatrixType> class EigenSolver
|
---|
65 | {
|
---|
66 | public:
|
---|
67 |
|
---|
68 | /** \brief Synonym for the template parameter \p _MatrixType. */
|
---|
69 | typedef _MatrixType MatrixType;
|
---|
70 |
|
---|
71 | enum {
|
---|
72 | RowsAtCompileTime = MatrixType::RowsAtCompileTime,
|
---|
73 | ColsAtCompileTime = MatrixType::ColsAtCompileTime,
|
---|
74 | Options = MatrixType::Options,
|
---|
75 | MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
|
---|
76 | MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
|
---|
77 | };
|
---|
78 |
|
---|
79 | /** \brief Scalar type for matrices of type #MatrixType. */
|
---|
80 | typedef typename MatrixType::Scalar Scalar;
|
---|
81 | typedef typename NumTraits<Scalar>::Real RealScalar;
|
---|
82 | typedef typename MatrixType::Index Index;
|
---|
83 |
|
---|
84 | /** \brief Complex scalar type for #MatrixType.
|
---|
85 | *
|
---|
86 | * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
|
---|
87 | * \c float or \c double) and just \c Scalar if #Scalar is
|
---|
88 | * complex.
|
---|
89 | */
|
---|
90 | typedef std::complex<RealScalar> ComplexScalar;
|
---|
91 |
|
---|
92 | /** \brief Type for vector of eigenvalues as returned by eigenvalues().
|
---|
93 | *
|
---|
94 | * This is a column vector with entries of type #ComplexScalar.
|
---|
95 | * The length of the vector is the size of #MatrixType.
|
---|
96 | */
|
---|
97 | typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
|
---|
98 |
|
---|
99 | /** \brief Type for matrix of eigenvectors as returned by eigenvectors().
|
---|
100 | *
|
---|
101 | * This is a square matrix with entries of type #ComplexScalar.
|
---|
102 | * The size is the same as the size of #MatrixType.
|
---|
103 | */
|
---|
104 | typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType;
|
---|
105 |
|
---|
106 | /** \brief Default constructor.
|
---|
107 | *
|
---|
108 | * The default constructor is useful in cases in which the user intends to
|
---|
109 | * perform decompositions via EigenSolver::compute(const MatrixType&, bool).
|
---|
110 | *
|
---|
111 | * \sa compute() for an example.
|
---|
112 | */
|
---|
113 | EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
|
---|
114 |
|
---|
115 | /** \brief Default constructor with memory preallocation
|
---|
116 | *
|
---|
117 | * Like the default constructor but with preallocation of the internal data
|
---|
118 | * according to the specified problem \a size.
|
---|
119 | * \sa EigenSolver()
|
---|
120 | */
|
---|
121 | EigenSolver(Index size)
|
---|
122 | : m_eivec(size, size),
|
---|
123 | m_eivalues(size),
|
---|
124 | m_isInitialized(false),
|
---|
125 | m_eigenvectorsOk(false),
|
---|
126 | m_realSchur(size),
|
---|
127 | m_matT(size, size),
|
---|
128 | m_tmp(size)
|
---|
129 | {}
|
---|
130 |
|
---|
131 | /** \brief Constructor; computes eigendecomposition of given matrix.
|
---|
132 | *
|
---|
133 | * \param[in] matrix Square matrix whose eigendecomposition is to be computed.
|
---|
134 | * \param[in] computeEigenvectors If true, both the eigenvectors and the
|
---|
135 | * eigenvalues are computed; if false, only the eigenvalues are
|
---|
136 | * computed.
|
---|
137 | *
|
---|
138 | * This constructor calls compute() to compute the eigenvalues
|
---|
139 | * and eigenvectors.
|
---|
140 | *
|
---|
141 | * Example: \include EigenSolver_EigenSolver_MatrixType.cpp
|
---|
142 | * Output: \verbinclude EigenSolver_EigenSolver_MatrixType.out
|
---|
143 | *
|
---|
144 | * \sa compute()
|
---|
145 | */
|
---|
146 | EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
|
---|
147 | : m_eivec(matrix.rows(), matrix.cols()),
|
---|
148 | m_eivalues(matrix.cols()),
|
---|
149 | m_isInitialized(false),
|
---|
150 | m_eigenvectorsOk(false),
|
---|
151 | m_realSchur(matrix.cols()),
|
---|
152 | m_matT(matrix.rows(), matrix.cols()),
|
---|
153 | m_tmp(matrix.cols())
|
---|
154 | {
|
---|
155 | compute(matrix, computeEigenvectors);
|
---|
156 | }
|
---|
157 |
|
---|
158 | /** \brief Returns the eigenvectors of given matrix.
|
---|
159 | *
|
---|
160 | * \returns %Matrix whose columns are the (possibly complex) eigenvectors.
|
---|
161 | *
|
---|
162 | * \pre Either the constructor
|
---|
163 | * EigenSolver(const MatrixType&,bool) or the member function
|
---|
164 | * compute(const MatrixType&, bool) has been called before, and
|
---|
165 | * \p computeEigenvectors was set to true (the default).
|
---|
166 | *
|
---|
167 | * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
|
---|
168 | * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The
|
---|
169 | * eigenvectors are normalized to have (Euclidean) norm equal to one. The
|
---|
170 | * matrix returned by this function is the matrix \f$ V \f$ in the
|
---|
171 | * eigendecomposition \f$ A = V D V^{-1} \f$, if it exists.
|
---|
172 | *
|
---|
173 | * Example: \include EigenSolver_eigenvectors.cpp
|
---|
174 | * Output: \verbinclude EigenSolver_eigenvectors.out
|
---|
175 | *
|
---|
176 | * \sa eigenvalues(), pseudoEigenvectors()
|
---|
177 | */
|
---|
178 | EigenvectorsType eigenvectors() const;
|
---|
179 |
|
---|
180 | /** \brief Returns the pseudo-eigenvectors of given matrix.
|
---|
181 | *
|
---|
182 | * \returns Const reference to matrix whose columns are the pseudo-eigenvectors.
|
---|
183 | *
|
---|
184 | * \pre Either the constructor
|
---|
185 | * EigenSolver(const MatrixType&,bool) or the member function
|
---|
186 | * compute(const MatrixType&, bool) has been called before, and
|
---|
187 | * \p computeEigenvectors was set to true (the default).
|
---|
188 | *
|
---|
189 | * The real matrix \f$ V \f$ returned by this function and the
|
---|
190 | * block-diagonal matrix \f$ D \f$ returned by pseudoEigenvalueMatrix()
|
---|
191 | * satisfy \f$ AV = VD \f$.
|
---|
192 | *
|
---|
193 | * Example: \include EigenSolver_pseudoEigenvectors.cpp
|
---|
194 | * Output: \verbinclude EigenSolver_pseudoEigenvectors.out
|
---|
195 | *
|
---|
196 | * \sa pseudoEigenvalueMatrix(), eigenvectors()
|
---|
197 | */
|
---|
198 | const MatrixType& pseudoEigenvectors() const
|
---|
199 | {
|
---|
200 | eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
|
---|
201 | eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
|
---|
202 | return m_eivec;
|
---|
203 | }
|
---|
204 |
|
---|
205 | /** \brief Returns the block-diagonal matrix in the pseudo-eigendecomposition.
|
---|
206 | *
|
---|
207 | * \returns A block-diagonal matrix.
|
---|
208 | *
|
---|
209 | * \pre Either the constructor
|
---|
210 | * EigenSolver(const MatrixType&,bool) or the member function
|
---|
211 | * compute(const MatrixType&, bool) has been called before.
|
---|
212 | *
|
---|
213 | * The matrix \f$ D \f$ returned by this function is real and
|
---|
214 | * block-diagonal. The blocks on the diagonal are either 1-by-1 or 2-by-2
|
---|
215 | * blocks of the form
|
---|
216 | * \f$ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f$.
|
---|
217 | * These blocks are not sorted in any particular order.
|
---|
218 | * The matrix \f$ D \f$ and the matrix \f$ V \f$ returned by
|
---|
219 | * pseudoEigenvectors() satisfy \f$ AV = VD \f$.
|
---|
220 | *
|
---|
221 | * \sa pseudoEigenvectors() for an example, eigenvalues()
|
---|
222 | */
|
---|
223 | MatrixType pseudoEigenvalueMatrix() const;
|
---|
224 |
|
---|
225 | /** \brief Returns the eigenvalues of given matrix.
|
---|
226 | *
|
---|
227 | * \returns A const reference to the column vector containing the eigenvalues.
|
---|
228 | *
|
---|
229 | * \pre Either the constructor
|
---|
230 | * EigenSolver(const MatrixType&,bool) or the member function
|
---|
231 | * compute(const MatrixType&, bool) has been called before.
|
---|
232 | *
|
---|
233 | * The eigenvalues are repeated according to their algebraic multiplicity,
|
---|
234 | * so there are as many eigenvalues as rows in the matrix. The eigenvalues
|
---|
235 | * are not sorted in any particular order.
|
---|
236 | *
|
---|
237 | * Example: \include EigenSolver_eigenvalues.cpp
|
---|
238 | * Output: \verbinclude EigenSolver_eigenvalues.out
|
---|
239 | *
|
---|
240 | * \sa eigenvectors(), pseudoEigenvalueMatrix(),
|
---|
241 | * MatrixBase::eigenvalues()
|
---|
242 | */
|
---|
243 | const EigenvalueType& eigenvalues() const
|
---|
244 | {
|
---|
245 | eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
|
---|
246 | return m_eivalues;
|
---|
247 | }
|
---|
248 |
|
---|
249 | /** \brief Computes eigendecomposition of given matrix.
|
---|
250 | *
|
---|
251 | * \param[in] matrix Square matrix whose eigendecomposition is to be computed.
|
---|
252 | * \param[in] computeEigenvectors If true, both the eigenvectors and the
|
---|
253 | * eigenvalues are computed; if false, only the eigenvalues are
|
---|
254 | * computed.
|
---|
255 | * \returns Reference to \c *this
|
---|
256 | *
|
---|
257 | * This function computes the eigenvalues of the real matrix \p matrix.
|
---|
258 | * The eigenvalues() function can be used to retrieve them. If
|
---|
259 | * \p computeEigenvectors is true, then the eigenvectors are also computed
|
---|
260 | * and can be retrieved by calling eigenvectors().
|
---|
261 | *
|
---|
262 | * The matrix is first reduced to real Schur form using the RealSchur
|
---|
263 | * class. The Schur decomposition is then used to compute the eigenvalues
|
---|
264 | * and eigenvectors.
|
---|
265 | *
|
---|
266 | * The cost of the computation is dominated by the cost of the
|
---|
267 | * Schur decomposition, which is very approximately \f$ 25n^3 \f$
|
---|
268 | * (where \f$ n \f$ is the size of the matrix) if \p computeEigenvectors
|
---|
269 | * is true, and \f$ 10n^3 \f$ if \p computeEigenvectors is false.
|
---|
270 | *
|
---|
271 | * This method reuses of the allocated data in the EigenSolver object.
|
---|
272 | *
|
---|
273 | * Example: \include EigenSolver_compute.cpp
|
---|
274 | * Output: \verbinclude EigenSolver_compute.out
|
---|
275 | */
|
---|
276 | EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
|
---|
277 |
|
---|
278 | ComputationInfo info() const
|
---|
279 | {
|
---|
280 | eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
|
---|
281 | return m_realSchur.info();
|
---|
282 | }
|
---|
283 |
|
---|
284 | /** \brief Sets the maximum number of iterations allowed. */
|
---|
285 | EigenSolver& setMaxIterations(Index maxIters)
|
---|
286 | {
|
---|
287 | m_realSchur.setMaxIterations(maxIters);
|
---|
288 | return *this;
|
---|
289 | }
|
---|
290 |
|
---|
291 | /** \brief Returns the maximum number of iterations. */
|
---|
292 | Index getMaxIterations()
|
---|
293 | {
|
---|
294 | return m_realSchur.getMaxIterations();
|
---|
295 | }
|
---|
296 |
|
---|
297 | private:
|
---|
298 | void doComputeEigenvectors();
|
---|
299 |
|
---|
300 | protected:
|
---|
301 |
|
---|
302 | static void check_template_parameters()
|
---|
303 | {
|
---|
304 | EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
|
---|
305 | EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
|
---|
306 | }
|
---|
307 |
|
---|
308 | MatrixType m_eivec;
|
---|
309 | EigenvalueType m_eivalues;
|
---|
310 | bool m_isInitialized;
|
---|
311 | bool m_eigenvectorsOk;
|
---|
312 | RealSchur<MatrixType> m_realSchur;
|
---|
313 | MatrixType m_matT;
|
---|
314 |
|
---|
315 | typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
|
---|
316 | ColumnVectorType m_tmp;
|
---|
317 | };
|
---|
318 |
|
---|
319 | template<typename MatrixType>
|
---|
320 | MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
|
---|
321 | {
|
---|
322 | eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
|
---|
323 | Index n = m_eivalues.rows();
|
---|
324 | MatrixType matD = MatrixType::Zero(n,n);
|
---|
325 | for (Index i=0; i<n; ++i)
|
---|
326 | {
|
---|
327 | if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i))))
|
---|
328 | matD.coeffRef(i,i) = numext::real(m_eivalues.coeff(i));
|
---|
329 | else
|
---|
330 | {
|
---|
331 | matD.template block<2,2>(i,i) << numext::real(m_eivalues.coeff(i)), numext::imag(m_eivalues.coeff(i)),
|
---|
332 | -numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i));
|
---|
333 | ++i;
|
---|
334 | }
|
---|
335 | }
|
---|
336 | return matD;
|
---|
337 | }
|
---|
338 |
|
---|
339 | template<typename MatrixType>
|
---|
340 | typename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const
|
---|
341 | {
|
---|
342 | eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
|
---|
343 | eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
|
---|
344 | Index n = m_eivec.cols();
|
---|
345 | EigenvectorsType matV(n,n);
|
---|
346 | for (Index j=0; j<n; ++j)
|
---|
347 | {
|
---|
348 | if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j))) || j+1==n)
|
---|
349 | {
|
---|
350 | // we have a real eigen value
|
---|
351 | matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
|
---|
352 | matV.col(j).normalize();
|
---|
353 | }
|
---|
354 | else
|
---|
355 | {
|
---|
356 | // we have a pair of complex eigen values
|
---|
357 | for (Index i=0; i<n; ++i)
|
---|
358 | {
|
---|
359 | matV.coeffRef(i,j) = ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1));
|
---|
360 | matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
|
---|
361 | }
|
---|
362 | matV.col(j).normalize();
|
---|
363 | matV.col(j+1).normalize();
|
---|
364 | ++j;
|
---|
365 | }
|
---|
366 | }
|
---|
367 | return matV;
|
---|
368 | }
|
---|
369 |
|
---|
370 | template<typename MatrixType>
|
---|
371 | EigenSolver<MatrixType>&
|
---|
372 | EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
|
---|
373 | {
|
---|
374 | check_template_parameters();
|
---|
375 |
|
---|
376 | using std::sqrt;
|
---|
377 | using std::abs;
|
---|
378 | eigen_assert(matrix.cols() == matrix.rows());
|
---|
379 |
|
---|
380 | // Reduce to real Schur form.
|
---|
381 | m_realSchur.compute(matrix, computeEigenvectors);
|
---|
382 |
|
---|
383 | if (m_realSchur.info() == Success)
|
---|
384 | {
|
---|
385 | m_matT = m_realSchur.matrixT();
|
---|
386 | if (computeEigenvectors)
|
---|
387 | m_eivec = m_realSchur.matrixU();
|
---|
388 |
|
---|
389 | // Compute eigenvalues from matT
|
---|
390 | m_eivalues.resize(matrix.cols());
|
---|
391 | Index i = 0;
|
---|
392 | while (i < matrix.cols())
|
---|
393 | {
|
---|
394 | if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0))
|
---|
395 | {
|
---|
396 | m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
|
---|
397 | ++i;
|
---|
398 | }
|
---|
399 | else
|
---|
400 | {
|
---|
401 | Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
|
---|
402 | Scalar z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
|
---|
403 | m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
|
---|
404 | m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
|
---|
405 | i += 2;
|
---|
406 | }
|
---|
407 | }
|
---|
408 |
|
---|
409 | // Compute eigenvectors.
|
---|
410 | if (computeEigenvectors)
|
---|
411 | doComputeEigenvectors();
|
---|
412 | }
|
---|
413 |
|
---|
414 | m_isInitialized = true;
|
---|
415 | m_eigenvectorsOk = computeEigenvectors;
|
---|
416 |
|
---|
417 | return *this;
|
---|
418 | }
|
---|
419 |
|
---|
420 | // Complex scalar division.
|
---|
421 | template<typename Scalar>
|
---|
422 | std::complex<Scalar> cdiv(const Scalar& xr, const Scalar& xi, const Scalar& yr, const Scalar& yi)
|
---|
423 | {
|
---|
424 | using std::abs;
|
---|
425 | Scalar r,d;
|
---|
426 | if (abs(yr) > abs(yi))
|
---|
427 | {
|
---|
428 | r = yi/yr;
|
---|
429 | d = yr + r*yi;
|
---|
430 | return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
|
---|
431 | }
|
---|
432 | else
|
---|
433 | {
|
---|
434 | r = yr/yi;
|
---|
435 | d = yi + r*yr;
|
---|
436 | return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
|
---|
437 | }
|
---|
438 | }
|
---|
439 |
|
---|
440 |
|
---|
441 | template<typename MatrixType>
|
---|
442 | void EigenSolver<MatrixType>::doComputeEigenvectors()
|
---|
443 | {
|
---|
444 | using std::abs;
|
---|
445 | const Index size = m_eivec.cols();
|
---|
446 | const Scalar eps = NumTraits<Scalar>::epsilon();
|
---|
447 |
|
---|
448 | // inefficient! this is already computed in RealSchur
|
---|
449 | Scalar norm(0);
|
---|
450 | for (Index j = 0; j < size; ++j)
|
---|
451 | {
|
---|
452 | norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
|
---|
453 | }
|
---|
454 |
|
---|
455 | // Backsubstitute to find vectors of upper triangular form
|
---|
456 | if (norm == 0.0)
|
---|
457 | {
|
---|
458 | return;
|
---|
459 | }
|
---|
460 |
|
---|
461 | for (Index n = size-1; n >= 0; n--)
|
---|
462 | {
|
---|
463 | Scalar p = m_eivalues.coeff(n).real();
|
---|
464 | Scalar q = m_eivalues.coeff(n).imag();
|
---|
465 |
|
---|
466 | // Scalar vector
|
---|
467 | if (q == Scalar(0))
|
---|
468 | {
|
---|
469 | Scalar lastr(0), lastw(0);
|
---|
470 | Index l = n;
|
---|
471 |
|
---|
472 | m_matT.coeffRef(n,n) = 1.0;
|
---|
473 | for (Index i = n-1; i >= 0; i--)
|
---|
474 | {
|
---|
475 | Scalar w = m_matT.coeff(i,i) - p;
|
---|
476 | Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
|
---|
477 |
|
---|
478 | if (m_eivalues.coeff(i).imag() < 0.0)
|
---|
479 | {
|
---|
480 | lastw = w;
|
---|
481 | lastr = r;
|
---|
482 | }
|
---|
483 | else
|
---|
484 | {
|
---|
485 | l = i;
|
---|
486 | if (m_eivalues.coeff(i).imag() == 0.0)
|
---|
487 | {
|
---|
488 | if (w != 0.0)
|
---|
489 | m_matT.coeffRef(i,n) = -r / w;
|
---|
490 | else
|
---|
491 | m_matT.coeffRef(i,n) = -r / (eps * norm);
|
---|
492 | }
|
---|
493 | else // Solve real equations
|
---|
494 | {
|
---|
495 | Scalar x = m_matT.coeff(i,i+1);
|
---|
496 | Scalar y = m_matT.coeff(i+1,i);
|
---|
497 | Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
|
---|
498 | Scalar t = (x * lastr - lastw * r) / denom;
|
---|
499 | m_matT.coeffRef(i,n) = t;
|
---|
500 | if (abs(x) > abs(lastw))
|
---|
501 | m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
|
---|
502 | else
|
---|
503 | m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
|
---|
504 | }
|
---|
505 |
|
---|
506 | // Overflow control
|
---|
507 | Scalar t = abs(m_matT.coeff(i,n));
|
---|
508 | if ((eps * t) * t > Scalar(1))
|
---|
509 | m_matT.col(n).tail(size-i) /= t;
|
---|
510 | }
|
---|
511 | }
|
---|
512 | }
|
---|
513 | else if (q < Scalar(0) && n > 0) // Complex vector
|
---|
514 | {
|
---|
515 | Scalar lastra(0), lastsa(0), lastw(0);
|
---|
516 | Index l = n-1;
|
---|
517 |
|
---|
518 | // Last vector component imaginary so matrix is triangular
|
---|
519 | if (abs(m_matT.coeff(n,n-1)) > abs(m_matT.coeff(n-1,n)))
|
---|
520 | {
|
---|
521 | m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
|
---|
522 | m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
|
---|
523 | }
|
---|
524 | else
|
---|
525 | {
|
---|
526 | std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
|
---|
527 | m_matT.coeffRef(n-1,n-1) = numext::real(cc);
|
---|
528 | m_matT.coeffRef(n-1,n) = numext::imag(cc);
|
---|
529 | }
|
---|
530 | m_matT.coeffRef(n,n-1) = 0.0;
|
---|
531 | m_matT.coeffRef(n,n) = 1.0;
|
---|
532 | for (Index i = n-2; i >= 0; i--)
|
---|
533 | {
|
---|
534 | Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
|
---|
535 | Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
|
---|
536 | Scalar w = m_matT.coeff(i,i) - p;
|
---|
537 |
|
---|
538 | if (m_eivalues.coeff(i).imag() < 0.0)
|
---|
539 | {
|
---|
540 | lastw = w;
|
---|
541 | lastra = ra;
|
---|
542 | lastsa = sa;
|
---|
543 | }
|
---|
544 | else
|
---|
545 | {
|
---|
546 | l = i;
|
---|
547 | if (m_eivalues.coeff(i).imag() == RealScalar(0))
|
---|
548 | {
|
---|
549 | std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
|
---|
550 | m_matT.coeffRef(i,n-1) = numext::real(cc);
|
---|
551 | m_matT.coeffRef(i,n) = numext::imag(cc);
|
---|
552 | }
|
---|
553 | else
|
---|
554 | {
|
---|
555 | // Solve complex equations
|
---|
556 | Scalar x = m_matT.coeff(i,i+1);
|
---|
557 | Scalar y = m_matT.coeff(i+1,i);
|
---|
558 | Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
|
---|
559 | Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
|
---|
560 | if ((vr == 0.0) && (vi == 0.0))
|
---|
561 | vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw));
|
---|
562 |
|
---|
563 | std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
|
---|
564 | m_matT.coeffRef(i,n-1) = numext::real(cc);
|
---|
565 | m_matT.coeffRef(i,n) = numext::imag(cc);
|
---|
566 | if (abs(x) > (abs(lastw) + abs(q)))
|
---|
567 | {
|
---|
568 | m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
|
---|
569 | m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
|
---|
570 | }
|
---|
571 | else
|
---|
572 | {
|
---|
573 | cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q);
|
---|
574 | m_matT.coeffRef(i+1,n-1) = numext::real(cc);
|
---|
575 | m_matT.coeffRef(i+1,n) = numext::imag(cc);
|
---|
576 | }
|
---|
577 | }
|
---|
578 |
|
---|
579 | // Overflow control
|
---|
580 | using std::max;
|
---|
581 | Scalar t = (max)(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n)));
|
---|
582 | if ((eps * t) * t > Scalar(1))
|
---|
583 | m_matT.block(i, n-1, size-i, 2) /= t;
|
---|
584 |
|
---|
585 | }
|
---|
586 | }
|
---|
587 |
|
---|
588 | // We handled a pair of complex conjugate eigenvalues, so need to skip them both
|
---|
589 | n--;
|
---|
590 | }
|
---|
591 | else
|
---|
592 | {
|
---|
593 | eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen
|
---|
594 | }
|
---|
595 | }
|
---|
596 |
|
---|
597 | // Back transformation to get eigenvectors of original matrix
|
---|
598 | for (Index j = size-1; j >= 0; j--)
|
---|
599 | {
|
---|
600 | m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
|
---|
601 | m_eivec.col(j) = m_tmp;
|
---|
602 | }
|
---|
603 | }
|
---|
604 |
|
---|
605 | } // end namespace Eigen
|
---|
606 |
|
---|
607 | #endif // EIGEN_EIGENSOLVER_H
|
---|