11#ifndef EIGEN_EIGENSOLVER_H
12#define EIGEN_EIGENSOLVER_H
80 typedef typename MatrixType::Scalar
Scalar;
146 template<
typename InputType>
148 :
m_eivec(matrix.rows(), matrix.cols()),
153 m_matT(matrix.rows(), matrix.cols()),
277 template<
typename InputType>
301 void doComputeEigenvectors();
323template<
typename MatrixType>
326 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
328 Index n = m_eivalues.rows();
330 for (
Index i=0; i<n; ++i)
344template<
typename MatrixType>
347 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
348 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
350 Index n = m_eivec.cols();
352 for (
Index j=0; j<n; ++j)
357 matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
358 matV.col(j).normalize();
363 for (
Index i=0; i<n; ++i)
368 matV.col(j).normalize();
369 matV.col(j+1).normalize();
376template<
typename MatrixType>
377template<
typename InputType>
381 check_template_parameters();
389 m_realSchur.compute(matrix.
derived(), computeEigenvectors);
391 m_info = m_realSchur.info();
395 m_matT = m_realSchur.matrixT();
396 if (computeEigenvectors)
397 m_eivec = m_realSchur.matrixU();
400 m_eivalues.resize(matrix.
cols());
402 while (i < matrix.
cols())
404 if (i == matrix.
cols() - 1 || m_matT.coeff(i+1, i) ==
Scalar(0))
406 m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
407 if(!(
isfinite)(m_eivalues.coeffRef(i)))
409 m_isInitialized =
true;
410 m_eigenvectorsOk =
false;
418 Scalar p =
Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
423 Scalar t0 = m_matT.coeff(i+1, i);
424 Scalar t1 = m_matT.coeff(i, i+1);
425 Scalar maxval = numext::maxi<Scalar>(
abs(p),numext::maxi<Scalar>(
abs(t0),
abs(t1)));
429 z = maxval *
sqrt(
abs(p0 * p0 + t0 * t1));
432 m_eivalues.coeffRef(i) =
ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
433 m_eivalues.coeffRef(i+1) =
ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
434 if(!((
isfinite)(m_eivalues.coeffRef(i)) && (
isfinite)(m_eivalues.coeffRef(i+1))))
436 m_isInitialized =
true;
437 m_eigenvectorsOk =
false;
446 if (computeEigenvectors)
447 doComputeEigenvectors();
450 m_isInitialized =
true;
451 m_eigenvectorsOk = computeEigenvectors;
457template<
typename MatrixType>
472 if (norm == Scalar(0))
479 Scalar p = m_eivalues.coeff(n).real();
480 Scalar q = m_eivalues.coeff(n).imag();
485 Scalar lastr(0), lastw(0);
488 m_matT.coeffRef(n,n) = Scalar(1);
489 for (
Index i = n-1; i >= 0; i--)
491 Scalar w = m_matT.coeff(i,i) - p;
492 Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
494 if (m_eivalues.coeff(i).imag() < Scalar(0))
502 if (m_eivalues.coeff(i).imag() == Scalar(0))
505 m_matT.coeffRef(i,n) = -r / w;
507 m_matT.coeffRef(i,n) = -r / (eps * norm);
511 Scalar x = m_matT.coeff(i,i+1);
512 Scalar
y = m_matT.coeff(i+1,i);
513 Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
514 Scalar t = (x * lastr - lastw * r) / denom;
515 m_matT.coeffRef(i,n) = t;
517 m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
519 m_matT.coeffRef(i+1,n) = (-lastr -
y * t) / lastw;
523 Scalar t =
abs(m_matT.coeff(i,n));
524 if ((eps * t) * t > Scalar(1))
525 m_matT.col(n).tail(
size-i) /= t;
529 else if (q < Scalar(0) && n > 0)
531 Scalar lastra(0), lastsa(0), lastw(0);
535 if (
abs(m_matT.coeff(n,n-1)) >
abs(m_matT.coeff(n-1,n)))
537 m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
538 m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
542 ComplexScalar cc = ComplexScalar(Scalar(0),-m_matT.coeff(n-1,n)) / ComplexScalar(m_matT.coeff(n-1,n-1)-p,q);
546 m_matT.coeffRef(n,n-1) = Scalar(0);
547 m_matT.coeffRef(n,n) = Scalar(1);
548 for (
Index i = n-2; i >= 0; i--)
550 Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
551 Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
552 Scalar w = m_matT.coeff(i,i) - p;
554 if (m_eivalues.coeff(i).imag() < Scalar(0))
563 if (m_eivalues.coeff(i).imag() == RealScalar(0))
565 ComplexScalar cc = ComplexScalar(-ra,-sa) / ComplexScalar(w,q);
572 Scalar x = m_matT.coeff(i,i+1);
573 Scalar
y = m_matT.coeff(i+1,i);
574 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;
575 Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
576 if ((vr == Scalar(0)) && (vi == Scalar(0)))
579 ComplexScalar cc = ComplexScalar(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra) / ComplexScalar(vr,vi);
584 m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
585 m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
589 cc = ComplexScalar(-lastra-
y*m_matT.coeff(i,n-1),-lastsa-
y*m_matT.coeff(i,n)) / ComplexScalar(lastw,q);
596 Scalar t = numext::maxi<Scalar>(
abs(m_matT.coeff(i,n-1)),
abs(m_matT.coeff(i,n)));
597 if ((eps * t) * t > Scalar(1))
598 m_matT.block(i, n-1,
size-i, 2) /= t;
608 eigen_assert(0 &&
"Internal bug in EigenSolver (INF or NaN has not been detected)");
615 m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
616 m_eivec.col(j) = m_tmp;
EIGEN_DEVICE_FUNC RealReturnType real() const
Definition: CommonCwiseUnaryOps.h:100
EIGEN_DEVICE_FUNC const ImagReturnType imag() const
Definition: CommonCwiseUnaryOps.h:109
#define eigen_assert(x)
Definition: Macros.h:1047
#define EIGEN_STATIC_ASSERT_NON_INTEGER(TYPE)
Definition: StaticAssert.h:187
#define EIGEN_STATIC_ASSERT(CONDITION, MSG)
Definition: StaticAssert.h:127
\eigenvalues_module
Definition: EigenSolver.h:65
@ Options
Definition: EigenSolver.h:74
@ RowsAtCompileTime
Definition: EigenSolver.h:72
@ ColsAtCompileTime
Definition: EigenSolver.h:73
@ MaxRowsAtCompileTime
Definition: EigenSolver.h:75
@ MaxColsAtCompileTime
Definition: EigenSolver.h:76
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: EigenSolver.h:80
ColumnVectorType m_tmp
Definition: EigenSolver.h:320
bool m_eigenvectorsOk
Definition: EigenSolver.h:314
bool m_isInitialized
Definition: EigenSolver.h:313
EigenSolver()
Default constructor.
Definition: EigenSolver.h:113
EigenSolver & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: EigenSolver.h:288
RealSchur< MatrixType > m_realSchur
Definition: EigenSolver.h:316
MatrixType pseudoEigenvalueMatrix() const
Returns the block-diagonal matrix in the pseudo-eigendecomposition.
Definition: EigenSolver.h:324
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Definition: EigenSolver.h:90
ComputationInfo m_info
Definition: EigenSolver.h:315
Eigen::Index Index
Definition: EigenSolver.h:82
EigenvalueType m_eivalues
Definition: EigenSolver.h:312
EigenvectorsType eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: EigenSolver.h:345
EigenSolver(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Constructor; computes eigendecomposition of given matrix.
Definition: EigenSolver.h:147
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
Definition: EigenSolver.h:69
const MatrixType & pseudoEigenvectors() const
Returns the pseudo-eigenvectors of given matrix.
Definition: EigenSolver.h:199
static void check_template_parameters()
Definition: EigenSolver.h:305
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorsType
Type for matrix of eigenvectors as returned by eigenvectors().
Definition: EigenSolver.h:104
EigenSolver(Index size)
Default constructor with memory preallocation.
Definition: EigenSolver.h:121
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: EigenSolver.h:295
MatrixType m_matT
Definition: EigenSolver.h:317
Matrix< Scalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > ColumnVectorType
Definition: EigenSolver.h:319
NumTraits< Scalar >::Real RealScalar
Definition: EigenSolver.h:81
MatrixType m_eivec
Definition: EigenSolver.h:311
ComputationInfo info() const
Definition: EigenSolver.h:281
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > EigenvalueType
Type for vector of eigenvalues as returned by eigenvalues().
Definition: EigenSolver.h:97
const EigenvalueType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: EigenSolver.h:244
EigenSolver & compute(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Computes eigendecomposition of given matrix.
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar & coeffRef(Index rowId, Index colId)
This is an overloaded version of DenseCoeffsBase<Derived,WriteAccessors>::coeffRef(Index,...
Definition: PlainObjectBase.h:175
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: RealSchur.h:213
RealSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: RealSchur.h:206
UnitType abs(const UnitType x) noexcept
Compute absolute value.
Definition: math.h:721
auto sqrt(const UnitType &value) noexcept -> unit_t< square_root< typename units::traits::unit_t_traits< UnitType >::unit_type >, typename units::traits::unit_t_traits< UnitType >::underlying_type, linear_scale >
computes the square root of value
Definition: math.h:483
ComputationInfo
Enum for reporting the status of a computation.
Definition: Constants.h:440
@ NumericalIssue
The provided data did not satisfy the prerequisites.
Definition: Constants.h:444
@ Success
Computation was successful.
Definition: Constants.h:442
constexpr common_t< T1, T2 > max(const T1 x, const T2 y) noexcept
Compile-time pairwise maximum function.
Definition: max.hpp:35
const Scalar & y
Definition: MathFunctions.h:821
EIGEN_DEVICE_FUNC bool isMuchSmallerThan(const Scalar &x, const OtherScalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
Definition: MathFunctions.h:1940
EIGEN_CONSTEXPR Index size(const T &x)
Definition: Meta.h:479
EIGEN_DEVICE_FUNC bool() isfinite(const T &x)
Definition: MathFunctions.h:1372
Namespace containing all symbols from the Eigen library.
Definition: MatrixExponential.h:16
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Common base class for all classes T such that MatrixBase has an operator=(T) and a constructor Matrix...
Definition: EigenBase.h:30
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:63
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:46
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
Holds information about the various numeric (i.e.
Definition: NumTraits.h:233