source: pacpussensors/trunk/Vislab/lib3dv/eigen/Eigen/src/Eigenvalues/EigenSolver.h@ 136

Last change on this file since 136 was 136, checked in by ldecherf, 7 years ago

Doc
File size: 21.3 KB
Line 
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
16namespace 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 */
64template<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
319template<typename MatrixType>
320MatrixType 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
339template<typename MatrixType>
340typename 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
370template<typename MatrixType>
371EigenSolver<MatrixType>&
372EigenSolver<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.
421template<typename Scalar>
422std::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
441template<typename MatrixType>
442void 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
Note: See TracBrowser for help on using the repository browser.