11#ifndef EIGEN_EIGENSOLVER_H
12#define EIGEN_EIGENSOLVER_H
14#include "./RealSchur.h"
16#include "./InternalHeaderCheck.h"
74 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
75 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
76 Options = MatrixType::Options,
77 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
78 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
82 typedef typename MatrixType::Scalar
Scalar;
115 EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_eigenvectorsOk(false), m_realSchur(), m_matT(), m_tmp() {}
124 : m_eivec(size, size),
126 m_isInitialized(false),
127 m_eigenvectorsOk(false),
148 template<
typename InputType>
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()),
203 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
204 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
248 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
279 template<
typename InputType>
285 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
292 m_realSchur.setMaxIterations(maxIters);
299 return m_realSchur.getMaxIterations();
303 void doComputeEigenvectors();
307 static void check_template_parameters()
309 EIGEN_STATIC_ASSERT_NON_INTEGER(
Scalar);
315 bool m_isInitialized;
316 bool m_eigenvectorsOk;
318 RealSchur<MatrixType> m_realSchur;
321 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
322 ColumnVectorType m_tmp;
325template<
typename MatrixType>
328 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
330 Index n = m_eivalues.rows();
332 for (
Index i=0; i<n; ++i)
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));
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));
346template<
typename MatrixType>
349 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
350 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
352 Index n = m_eivec.cols();
354 for (
Index j=0; j<n; ++j)
356 if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j)), precision) || j+1==n)
359 matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
360 matV.col(j).normalize();
365 for (
Index i=0; i<n; ++i)
370 matV.col(j).normalize();
371 matV.col(j+1).normalize();
378template<
typename MatrixType>
379template<
typename InputType>
383 check_template_parameters();
387 using numext::isfinite;
388 eigen_assert(matrix.
cols() == matrix.
rows());
391 m_realSchur.compute(matrix.
derived(), computeEigenvectors);
393 m_info = m_realSchur.info();
397 m_matT = m_realSchur.matrixT();
398 if (computeEigenvectors)
399 m_eivec = m_realSchur.matrixU();
402 m_eivalues.resize(matrix.
cols());
404 while (i < matrix.
cols())
406 if (i == matrix.
cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0))
408 m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
409 if(!(
isfinite)(m_eivalues.coeffRef(i)))
411 m_isInitialized =
true;
412 m_eigenvectorsOk =
false;
420 Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
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)));
430 Scalar p0 = p/maxval;
431 z = maxval *
sqrt(
abs(p0 * p0 + t0 * t1));
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))))
438 m_isInitialized =
true;
439 m_eigenvectorsOk =
false;
448 if (computeEigenvectors)
449 doComputeEigenvectors();
452 m_isInitialized =
true;
453 m_eigenvectorsOk = computeEigenvectors;
459template<
typename MatrixType>
460void EigenSolver<MatrixType>::doComputeEigenvectors()
463 const Index size = m_eivec.cols();
464 const Scalar eps = NumTraits<Scalar>::epsilon();
468 for (
Index j = 0; j < size; ++j)
470 norm += m_matT.row(j).segment((std::max)(j-1,
Index(0)), size-(std::max)(j-1,
Index(0))).cwiseAbs().sum();
474 if (norm == Scalar(0))
479 for (
Index n = size-1; n >= 0; n--)
481 Scalar p = m_eivalues.coeff(n).real();
482 Scalar q = m_eivalues.coeff(n).imag();
487 Scalar lastr(0), lastw(0);
490 m_matT.coeffRef(n,n) = Scalar(1);
491 for (
Index i = n-1; i >= 0; i--)
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));
496 if (m_eivalues.coeff(i).imag() < Scalar(0))
504 if (m_eivalues.coeff(i).imag() == Scalar(0))
507 m_matT.coeffRef(i,n) = -r / w;
509 m_matT.coeffRef(i,n) = -r / (eps * norm);
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;
519 m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
521 m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
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;
531 else if (q < Scalar(0) && n > 0)
533 Scalar lastra(0), lastsa(0), lastw(0);
537 if (
abs(m_matT.coeff(n,n-1)) >
abs(m_matT.coeff(n-1,n)))
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);
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);
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--)
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;
556 if (m_eivalues.coeff(i).imag() < Scalar(0))
565 if (m_eivalues.coeff(i).imag() == RealScalar(0))
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);
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)))
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);
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;
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);
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;
610 eigen_assert(0 &&
"Internal bug in EigenSolver (INF or NaN has not been detected)");
615 for (
Index j = size-1; j >= 0; j--)
617 m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
618 m_eivec.col(j) = m_tmp;
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