Eigen  3.4.90 (git rev a4098ac676528a83cfb73d4d26ce1b42ec05f47c)
EigenSolver.h
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#include "./InternalHeaderCheck.h"
17
18namespace Eigen {
19
66template<typename MatrixType_> class EigenSolver
67{
68 public:
69
71 typedef MatrixType_ MatrixType;
72
73 enum {
74 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
75 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
76 Options = MatrixType::Options,
77 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
78 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
79 };
80
82 typedef typename MatrixType::Scalar Scalar;
83 typedef typename NumTraits<Scalar>::Real RealScalar;
85
92 typedef std::complex<RealScalar> ComplexScalar;
93
100
107
115 EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_eigenvectorsOk(false), m_realSchur(), m_matT(), m_tmp() {}
116
123 explicit EigenSolver(Index size)
124 : m_eivec(size, size),
125 m_eivalues(size),
126 m_isInitialized(false),
127 m_eigenvectorsOk(false),
128 m_realSchur(size),
129 m_matT(size, size),
130 m_tmp(size)
131 {}
132
148 template<typename InputType>
149 explicit EigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true)
150 : m_eivec(matrix.rows(), matrix.cols()),
151 m_eivalues(matrix.cols()),
152 m_isInitialized(false),
153 m_eigenvectorsOk(false),
154 m_realSchur(matrix.cols()),
155 m_matT(matrix.rows(), matrix.cols()),
156 m_tmp(matrix.cols())
157 {
158 compute(matrix.derived(), computeEigenvectors);
159 }
160
182
202 {
203 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
204 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
205 return m_eivec;
206 }
207
227
247 {
248 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
249 return m_eivalues;
250 }
251
279 template<typename InputType>
280 EigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true);
281
284 {
285 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
286 return m_info;
287 }
288
291 {
292 m_realSchur.setMaxIterations(maxIters);
293 return *this;
294 }
295
298 {
299 return m_realSchur.getMaxIterations();
300 }
301
302 private:
303 void doComputeEigenvectors();
304
305 protected:
306
307 static void check_template_parameters()
308 {
309 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
310 EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
311 }
312
313 MatrixType m_eivec;
314 EigenvalueType m_eivalues;
315 bool m_isInitialized;
316 bool m_eigenvectorsOk;
317 ComputationInfo m_info;
318 RealSchur<MatrixType> m_realSchur;
319 MatrixType m_matT;
320
321 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
322 ColumnVectorType m_tmp;
323};
324
325template<typename MatrixType>
327{
328 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
329 const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon();
330 Index n = m_eivalues.rows();
331 MatrixType matD = MatrixType::Zero(n,n);
332 for (Index i=0; i<n; ++i)
333 {
334 if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i)), precision))
335 matD.coeffRef(i,i) = numext::real(m_eivalues.coeff(i));
336 else
337 {
338 matD.template block<2,2>(i,i) << numext::real(m_eivalues.coeff(i)), numext::imag(m_eivalues.coeff(i)),
339 -numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i));
340 ++i;
341 }
342 }
343 return matD;
344}
345
346template<typename MatrixType>
348{
349 eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
350 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
351 const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon();
352 Index n = m_eivec.cols();
353 EigenvectorsType matV(n,n);
354 for (Index j=0; j<n; ++j)
355 {
356 if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j)), precision) || j+1==n)
357 {
358 // we have a real eigen value
359 matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
360 matV.col(j).normalize();
361 }
362 else
363 {
364 // we have a pair of complex eigen values
365 for (Index i=0; i<n; ++i)
366 {
367 matV.coeffRef(i,j) = ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1));
368 matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
369 }
370 matV.col(j).normalize();
371 matV.col(j+1).normalize();
372 ++j;
373 }
374 }
375 return matV;
376}
377
378template<typename MatrixType>
379template<typename InputType>
381EigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeEigenvectors)
382{
383 check_template_parameters();
384
385 using std::sqrt;
386 using std::abs;
387 using numext::isfinite;
388 eigen_assert(matrix.cols() == matrix.rows());
389
390 // Reduce to real Schur form.
391 m_realSchur.compute(matrix.derived(), computeEigenvectors);
392
393 m_info = m_realSchur.info();
394
395 if (m_info == Success)
396 {
397 m_matT = m_realSchur.matrixT();
398 if (computeEigenvectors)
399 m_eivec = m_realSchur.matrixU();
400
401 // Compute eigenvalues from matT
402 m_eivalues.resize(matrix.cols());
403 Index i = 0;
404 while (i < matrix.cols())
405 {
406 if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0))
407 {
408 m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
409 if(!(isfinite)(m_eivalues.coeffRef(i)))
410 {
411 m_isInitialized = true;
412 m_eigenvectorsOk = false;
413 m_info = NumericalIssue;
414 return *this;
415 }
416 ++i;
417 }
418 else
419 {
420 Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
421 Scalar z;
422 // Compute z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
423 // without overflow
424 {
425 Scalar t0 = m_matT.coeff(i+1, i);
426 Scalar t1 = m_matT.coeff(i, i+1);
427 Scalar maxval = numext::maxi<Scalar>(abs(p),numext::maxi<Scalar>(abs(t0),abs(t1)));
428 t0 /= maxval;
429 t1 /= maxval;
430 Scalar p0 = p/maxval;
431 z = maxval * sqrt(abs(p0 * p0 + t0 * t1));
432 }
433
434 m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
435 m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
436 if(!((isfinite)(m_eivalues.coeffRef(i)) && (isfinite)(m_eivalues.coeffRef(i+1))))
437 {
438 m_isInitialized = true;
439 m_eigenvectorsOk = false;
440 m_info = NumericalIssue;
441 return *this;
442 }
443 i += 2;
444 }
445 }
446
447 // Compute eigenvectors.
448 if (computeEigenvectors)
449 doComputeEigenvectors();
450 }
451
452 m_isInitialized = true;
453 m_eigenvectorsOk = computeEigenvectors;
454
455 return *this;
456}
457
458
459template<typename MatrixType>
460void EigenSolver<MatrixType>::doComputeEigenvectors()
461{
462 using std::abs;
463 const Index size = m_eivec.cols();
464 const Scalar eps = NumTraits<Scalar>::epsilon();
465
466 // inefficient! this is already computed in RealSchur
467 Scalar norm(0);
468 for (Index j = 0; j < size; ++j)
469 {
470 norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
471 }
472
473 // Backsubstitute to find vectors of upper triangular form
474 if (norm == Scalar(0))
475 {
476 return;
477 }
478
479 for (Index n = size-1; n >= 0; n--)
480 {
481 Scalar p = m_eivalues.coeff(n).real();
482 Scalar q = m_eivalues.coeff(n).imag();
483
484 // Scalar vector
485 if (q == Scalar(0))
486 {
487 Scalar lastr(0), lastw(0);
488 Index l = n;
489
490 m_matT.coeffRef(n,n) = Scalar(1);
491 for (Index i = n-1; i >= 0; i--)
492 {
493 Scalar w = m_matT.coeff(i,i) - p;
494 Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
495
496 if (m_eivalues.coeff(i).imag() < Scalar(0))
497 {
498 lastw = w;
499 lastr = r;
500 }
501 else
502 {
503 l = i;
504 if (m_eivalues.coeff(i).imag() == Scalar(0))
505 {
506 if (w != Scalar(0))
507 m_matT.coeffRef(i,n) = -r / w;
508 else
509 m_matT.coeffRef(i,n) = -r / (eps * norm);
510 }
511 else // Solve real equations
512 {
513 Scalar x = m_matT.coeff(i,i+1);
514 Scalar y = m_matT.coeff(i+1,i);
515 Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
516 Scalar t = (x * lastr - lastw * r) / denom;
517 m_matT.coeffRef(i,n) = t;
518 if (abs(x) > abs(lastw))
519 m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
520 else
521 m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
522 }
523
524 // Overflow control
525 Scalar t = abs(m_matT.coeff(i,n));
526 if ((eps * t) * t > Scalar(1))
527 m_matT.col(n).tail(size-i) /= t;
528 }
529 }
530 }
531 else if (q < Scalar(0) && n > 0) // Complex vector
532 {
533 Scalar lastra(0), lastsa(0), lastw(0);
534 Index l = n-1;
535
536 // Last vector component imaginary so matrix is triangular
537 if (abs(m_matT.coeff(n,n-1)) > abs(m_matT.coeff(n-1,n)))
538 {
539 m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
540 m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
541 }
542 else
543 {
544 ComplexScalar cc = ComplexScalar(Scalar(0),-m_matT.coeff(n-1,n)) / ComplexScalar(m_matT.coeff(n-1,n-1)-p,q);
545 m_matT.coeffRef(n-1,n-1) = numext::real(cc);
546 m_matT.coeffRef(n-1,n) = numext::imag(cc);
547 }
548 m_matT.coeffRef(n,n-1) = Scalar(0);
549 m_matT.coeffRef(n,n) = Scalar(1);
550 for (Index i = n-2; i >= 0; i--)
551 {
552 Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
553 Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
554 Scalar w = m_matT.coeff(i,i) - p;
555
556 if (m_eivalues.coeff(i).imag() < Scalar(0))
557 {
558 lastw = w;
559 lastra = ra;
560 lastsa = sa;
561 }
562 else
563 {
564 l = i;
565 if (m_eivalues.coeff(i).imag() == RealScalar(0))
566 {
567 ComplexScalar cc = ComplexScalar(-ra,-sa) / ComplexScalar(w,q);
568 m_matT.coeffRef(i,n-1) = numext::real(cc);
569 m_matT.coeffRef(i,n) = numext::imag(cc);
570 }
571 else
572 {
573 // Solve complex equations
574 Scalar x = m_matT.coeff(i,i+1);
575 Scalar y = m_matT.coeff(i+1,i);
576 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;
577 Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
578 if ((vr == Scalar(0)) && (vi == Scalar(0)))
579 vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw));
580
581 ComplexScalar cc = ComplexScalar(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra) / ComplexScalar(vr,vi);
582 m_matT.coeffRef(i,n-1) = numext::real(cc);
583 m_matT.coeffRef(i,n) = numext::imag(cc);
584 if (abs(x) > (abs(lastw) + abs(q)))
585 {
586 m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
587 m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
588 }
589 else
590 {
591 cc = ComplexScalar(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n)) / ComplexScalar(lastw,q);
592 m_matT.coeffRef(i+1,n-1) = numext::real(cc);
593 m_matT.coeffRef(i+1,n) = numext::imag(cc);
594 }
595 }
596
597 // Overflow control
598 Scalar t = numext::maxi<Scalar>(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n)));
599 if ((eps * t) * t > Scalar(1))
600 m_matT.block(i, n-1, size-i, 2) /= t;
601
602 }
603 }
604
605 // We handled a pair of complex conjugate eigenvalues, so need to skip them both
606 n--;
607 }
608 else
609 {
610 eigen_assert(0 && "Internal bug in EigenSolver (INF or NaN has not been detected)"); // this should not happen
611 }
612 }
613
614 // Back transformation to get eigenvectors of original matrix
615 for (Index j = size-1; j >= 0; j--)
616 {
617 m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
618 m_eivec.col(j) = m_tmp;
619 }
620}
621
622} // end namespace Eigen
623
624#endif // EIGEN_EIGENSOLVER_H
Computes eigenvalues and eigenvectors of general matrices.
Definition: EigenSolver.h:67
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: EigenSolver.h:82
EigenSolver(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Constructor; computes eigendecomposition of given matrix.
Definition: EigenSolver.h:149
EigenSolver(Index size)
Default constructor with memory preallocation.
Definition: EigenSolver.h:123
MatrixType pseudoEigenvalueMatrix() const
Returns the block-diagonal matrix in the pseudo-eigendecomposition.
Definition: EigenSolver.h:326
ComputationInfo info() const
Definition: EigenSolver.h:283
EigenvectorsType eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: EigenSolver.h:347
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Definition: EigenSolver.h:92
EigenSolver()
Default constructor.
Definition: EigenSolver.h:115
EigenSolver & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: EigenSolver.h:290
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > EigenvalueType
Type for vector of eigenvalues as returned by eigenvalues().
Definition: EigenSolver.h:99
Eigen::Index Index
Definition: EigenSolver.h:84
const MatrixType & pseudoEigenvectors() const
Returns the pseudo-eigenvectors of given matrix.
Definition: EigenSolver.h:201
EigenSolver & compute(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Computes eigendecomposition of given matrix.
const EigenvalueType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: EigenSolver.h:246
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: EigenSolver.h:297
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorsType
Type for matrix of eigenvectors as returned by eigenvectors().
Definition: EigenSolver.h:106
MatrixType_ MatrixType
Synonym for the template parameter MatrixType_.
Definition: EigenSolver.h:71
Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:187
ComputationInfo
Definition: Constants.h:442
@ NumericalIssue
Definition: Constants.h:446
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: B01_Experimental.dox:1
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_isfinite_op< typename Derived::Scalar >, const Derived > isfinite(const Eigen::ArrayBase< Derived > &x)
Definition: EigenBase.h:32
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:65
Derived & derived()
Definition: EigenBase.h:48
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:62
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:235