Eigen::EigenSolver< MatrixType_ > Class Template Reference

Computes eigenvalues and eigenvectors of general matrices. More...

#include <EigenSolver.h>

Public Types

enum  {
  RowsAtCompileTime = MatrixType::RowsAtCompileTime , ColsAtCompileTime = MatrixType::ColsAtCompileTime , Options = internal::traits<MatrixType>::Options , MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime ,
  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
}
 
typedef MatrixType_ MatrixType
 Synonym for the template parameter MatrixType_. More...
 
typedef MatrixType::Scalar Scalar
 Scalar type for matrices of type MatrixType. More...
 
typedef NumTraits< Scalar >::Real RealScalar
 
typedef Eigen::Index Index
 
typedef std::complex< RealScalarComplexScalar
 Complex scalar type for MatrixType. More...
 
typedef Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > EigenvalueType
 Type for vector of eigenvalues as returned by eigenvalues(). More...
 
typedef Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTimeEigenvectorsType
 Type for matrix of eigenvectors as returned by eigenvectors(). More...
 

Public Member Functions

 EigenSolver ()
 Default constructor. More...
 
 EigenSolver (Index size)
 Default constructor with memory preallocation. More...
 
template<typename InputType >
 EigenSolver (const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
 Constructor; computes eigendecomposition of given matrix. More...
 
EigenvectorsType eigenvectors () const
 Returns the eigenvectors of given matrix. More...
 
const MatrixTypepseudoEigenvectors () const
 Returns the pseudo-eigenvectors of given matrix. More...
 
MatrixType pseudoEigenvalueMatrix () const
 Returns the block-diagonal matrix in the pseudo-eigendecomposition. More...
 
const EigenvalueTypeeigenvalues () const
 Returns the eigenvalues of given matrix. More...
 
template<typename InputType >
EigenSolvercompute (const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
 Computes eigendecomposition of given matrix. More...
 
ComputationInfo info () const
 
EigenSolversetMaxIterations (Index maxIters)
 Sets the maximum number of iterations allowed. More...
 
Index getMaxIterations ()
 Returns the maximum number of iterations. More...
 
template<typename InputType >
EigenSolver< MatrixType > & compute (const EigenBase< InputType > &matrix, bool computeEigenvectors)
 

Protected Types

typedef Matrix< Scalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > ColumnVectorType
 

Static Protected Member Functions

static void check_template_parameters ()
 

Protected Attributes

MatrixType m_eivec
 
EigenvalueType m_eivalues
 
bool m_isInitialized
 
bool m_eigenvectorsOk
 
ComputationInfo m_info
 
RealSchur< MatrixTypem_realSchur
 
MatrixType m_matT
 
ColumnVectorType m_tmp
 

Private Member Functions

void doComputeEigenvectors ()
 

Detailed Description

template<typename MatrixType_>
class Eigen::EigenSolver< MatrixType_ >

Computes eigenvalues and eigenvectors of general matrices.

\eigenvalues_module

Template Parameters
MatrixType_the type of the matrix of which we are computing the eigendecomposition; this is expected to be an instantiation of the Matrix class template. Currently, only real matrices are supported.

The eigenvalues and eigenvectors of a matrix \( A \) are scalars \( \lambda \) and vectors \( v \) such that \( Av = \lambda v \). If \( D \) is a diagonal matrix with the eigenvalues on the diagonal, and \( V \) is a matrix with the eigenvectors as its columns, then \( A V = V D \). The matrix \( V \) is almost always invertible, in which case we have \( A = V D V^{-1} \). This is called the eigendecomposition.

The eigenvalues and eigenvectors of a matrix may be complex, even when the matrix is real. However, we can choose real matrices \( V \) and \( D \) satisfying \( A V = V D \), just like the eigendecomposition, if the matrix \( D \) is not required to be diagonal, but if it is allowed to have blocks of the form

\[ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \]

(where \( u \) and \( v \) are real numbers) on the diagonal. These blocks correspond to complex eigenvalue pairs \( u \pm iv \). We call this variant of the eigendecomposition the pseudo-eigendecomposition.

Call the function compute() to compute the eigenvalues and eigenvectors of a given matrix. Alternatively, you can use the EigenSolver(const MatrixType&, bool) constructor which computes the eigenvalues and eigenvectors at construction time. Once the eigenvalue and eigenvectors are computed, they can be retrieved with the eigenvalues() and eigenvectors() functions. The pseudoEigenvalueMatrix() and pseudoEigenvectors() methods allow the construction of the pseudo-eigendecomposition.

The documentation for EigenSolver(const MatrixType&, bool) contains an example of the typical use of this class.

Note
The implementation is adapted from JAMA (public domain). Their code is based on EISPACK.
See also
MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver

Member Typedef Documentation

◆ ColumnVectorType

template<typename MatrixType_ >
typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> Eigen::EigenSolver< MatrixType_ >::ColumnVectorType
protected

◆ ComplexScalar

template<typename MatrixType_ >
typedef std::complex<RealScalar> Eigen::EigenSolver< MatrixType_ >::ComplexScalar

Complex scalar type for MatrixType.

This is std::complex<Scalar> if Scalar is real (e.g., float or double) and just Scalar if Scalar is complex.

◆ EigenvalueType

template<typename MatrixType_ >
typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> Eigen::EigenSolver< MatrixType_ >::EigenvalueType

Type for vector of eigenvalues as returned by eigenvalues().

This is a column vector with entries of type ComplexScalar. The length of the vector is the size of MatrixType.

◆ EigenvectorsType

Type for matrix of eigenvectors as returned by eigenvectors().

This is a square matrix with entries of type ComplexScalar. The size is the same as the size of MatrixType.

◆ Index

template<typename MatrixType_ >
typedef Eigen::Index Eigen::EigenSolver< MatrixType_ >::Index
Deprecated:
since Eigen 3.3

◆ MatrixType

template<typename MatrixType_ >
typedef MatrixType_ Eigen::EigenSolver< MatrixType_ >::MatrixType

Synonym for the template parameter MatrixType_.

◆ RealScalar

template<typename MatrixType_ >
typedef NumTraits<Scalar>::Real Eigen::EigenSolver< MatrixType_ >::RealScalar

◆ Scalar

template<typename MatrixType_ >
typedef MatrixType::Scalar Eigen::EigenSolver< MatrixType_ >::Scalar

Scalar type for matrices of type MatrixType.

Member Enumeration Documentation

◆ anonymous enum

template<typename MatrixType_ >
anonymous enum
Enumerator
RowsAtCompileTime 
ColsAtCompileTime 
Options 
MaxRowsAtCompileTime 
MaxColsAtCompileTime 
73  {
74  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
75  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
76  Options = internal::traits<MatrixType>::Options,
77  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
78  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
79  };
@ MaxRowsAtCompileTime
Definition: EigenSolver.h:77
@ MaxColsAtCompileTime
Definition: EigenSolver.h:78
@ ColsAtCompileTime
Definition: EigenSolver.h:75
@ Options
Definition: EigenSolver.h:76
@ RowsAtCompileTime
Definition: EigenSolver.h:74

Constructor & Destructor Documentation

◆ EigenSolver() [1/3]

template<typename MatrixType_ >
Eigen::EigenSolver< MatrixType_ >::EigenSolver ( )
inline

Default constructor.

The default constructor is useful in cases in which the user intends to perform decompositions via EigenSolver::compute(const MatrixType&, bool).

See also
compute() for an example.
118  : m_eivec(), m_eivalues(), m_isInitialized(false), m_eigenvectorsOk(false), m_realSchur(), m_matT(), m_tmp() {}
MatrixType m_eivec
Definition: EigenSolver.h:306
RealSchur< MatrixType > m_realSchur
Definition: EigenSolver.h:311
bool m_eigenvectorsOk
Definition: EigenSolver.h:309
bool m_isInitialized
Definition: EigenSolver.h:308
MatrixType m_matT
Definition: EigenSolver.h:312
EigenvalueType m_eivalues
Definition: EigenSolver.h:307
ColumnVectorType m_tmp
Definition: EigenSolver.h:315

◆ EigenSolver() [2/3]

template<typename MatrixType_ >
Eigen::EigenSolver< MatrixType_ >::EigenSolver ( Index  size)
inlineexplicit

Default constructor with memory preallocation.

Like the default constructor but with preallocation of the internal data according to the specified problem size.

See also
EigenSolver()
127  : m_eivec(size, size),
128  m_eivalues(size),
129  m_isInitialized(false),
130  m_eigenvectorsOk(false),
131  m_realSchur(size),
132  m_matT(size, size),
133  m_tmp(size) {}
Scalar Scalar int size
Definition: benchVecAdd.cpp:17

◆ EigenSolver() [3/3]

template<typename MatrixType_ >
template<typename InputType >
Eigen::EigenSolver< MatrixType_ >::EigenSolver ( const EigenBase< InputType > &  matrix,
bool  computeEigenvectors = true 
)
inlineexplicit

Constructor; computes eigendecomposition of given matrix.

Parameters
[in]matrixSquare matrix whose eigendecomposition is to be computed.
[in]computeEigenvectorsIf true, both the eigenvectors and the eigenvalues are computed; if false, only the eigenvalues are computed.

This constructor calls compute() to compute the eigenvalues and eigenvectors.

Example:

MatrixXd A = MatrixXd::Random(6, 6);
cout << "Here is a random 6x6 matrix, A:" << endl << A << endl << endl;
EigenSolver<MatrixXd> es(A);
cout << "The eigenvalues of A are:" << endl << es.eigenvalues() << endl;
cout << "The matrix of eigenvectors, V, is:" << endl << es.eigenvectors() << endl << endl;
complex<double> lambda = es.eigenvalues()[0];
cout << "Consider the first eigenvalue, lambda = " << lambda << endl;
VectorXcd v = es.eigenvectors().col(0);
cout << "If v is the corresponding eigenvector, then lambda * v = " << endl << lambda * v << endl;
cout << "... and A * v = " << endl << A.cast<complex<double> >() * v << endl << endl;
MatrixXcd D = es.eigenvalues().asDiagonal();
MatrixXcd V = es.eigenvectors();
cout << "Finally, V * D * V^(-1) = " << endl << V * D * V.inverse() << endl;
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
cout<< "The eigenvalues of A are:"<< endl<< ces.eigenvalues()<< endl;cout<< "The matrix of eigenvectors, V, is:"<< endl<< ces.eigenvectors()<< endl<< endl;complex< float > lambda
Definition: ComplexEigenSolver_compute.cpp:9
dominoes D
Definition: Domino.cpp:55
MatrixXcd V
Definition: EigenSolver_EigenSolver_MatrixType.cpp:15
EigenSolver< MatrixXf > es
Definition: EigenSolver_compute.cpp:1
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:47
Definition: datatypes.h:12

Output:

See also
compute()
152  : m_eivec(matrix.rows(), matrix.cols()),
153  m_eivalues(matrix.cols()),
154  m_isInitialized(false),
155  m_eigenvectorsOk(false),
156  m_realSchur(matrix.cols()),
157  m_matT(matrix.rows(), matrix.cols()),
158  m_tmp(matrix.cols()) {
159  compute(matrix.derived(), computeEigenvectors);
160  }
EigenSolver & compute(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Computes eigendecomposition of given matrix.
Eigen::Map< Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor >, 0, Eigen::OuterStride<> > matrix(T *data, int rows, int cols, int stride)
Definition: common.h:85

References Eigen::EigenSolver< MatrixType_ >::compute(), and matrix().

Member Function Documentation

◆ check_template_parameters()

template<typename MatrixType_ >
static void Eigen::EigenSolver< MatrixType_ >::check_template_parameters ( )
inlinestaticprotected
301  {
303  EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
304  }
#define EIGEN_STATIC_ASSERT(X, MSG)
Definition: StaticAssert.h:26
#define EIGEN_STATIC_ASSERT_NON_INTEGER(TYPE)
Definition: StaticAssert.h:74
SCALAR Scalar
Definition: bench_gemm.cpp:45
@ IsComplex
Definition: NumTraits.h:176

References EIGEN_STATIC_ASSERT, and EIGEN_STATIC_ASSERT_NON_INTEGER.

◆ compute() [1/2]

template<typename MatrixType_ >
template<typename InputType >
EigenSolver<MatrixType>& Eigen::EigenSolver< MatrixType_ >::compute ( const EigenBase< InputType > &  matrix,
bool  computeEigenvectors 
)
373  {
375 
376  using numext::isfinite;
377  using std::abs;
378  using std::sqrt;
379  eigen_assert(matrix.cols() == matrix.rows());
380 
381  // Reduce to real Schur form.
382  m_realSchur.compute(matrix.derived(), computeEigenvectors);
383 
384  m_info = m_realSchur.info();
385 
386  if (m_info == Success) {
388  if (computeEigenvectors) m_eivec = m_realSchur.matrixU();
389 
390  // Compute eigenvalues from matT
391  m_eivalues.resize(matrix.cols());
392  Index i = 0;
393  while (i < matrix.cols()) {
394  if (i == matrix.cols() - 1 || m_matT.coeff(i + 1, i) == Scalar(0)) {
395  m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
396  if (!(isfinite)(m_eivalues.coeffRef(i))) {
397  m_isInitialized = true;
398  m_eigenvectorsOk = false;
400  return *this;
401  }
402  ++i;
403  } else {
404  Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i + 1, i + 1));
405  Scalar z;
406  // Compute z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
407  // without overflow
408  {
409  Scalar t0 = m_matT.coeff(i + 1, i);
410  Scalar t1 = m_matT.coeff(i, i + 1);
411  Scalar maxval = numext::maxi<Scalar>(abs(p), numext::maxi<Scalar>(abs(t0), abs(t1)));
412  t0 /= maxval;
413  t1 /= maxval;
414  Scalar p0 = p / maxval;
415  z = maxval * sqrt(abs(p0 * p0 + t0 * t1));
416  }
417 
418  m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i + 1, i + 1) + p, z);
419  m_eivalues.coeffRef(i + 1) = ComplexScalar(m_matT.coeff(i + 1, i + 1) + p, -z);
420  if (!((isfinite)(m_eivalues.coeffRef(i)) && (isfinite)(m_eivalues.coeffRef(i + 1)))) {
421  m_isInitialized = true;
422  m_eigenvectorsOk = false;
424  return *this;
425  }
426  i += 2;
427  }
428  }
429 
430  // Compute eigenvectors.
431  if (computeEigenvectors) doComputeEigenvectors();
432  }
433 
434  m_isInitialized = true;
435  m_eigenvectorsOk = computeEigenvectors;
436 
437  return *this;
438 }
AnnoyingScalar abs(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:135
AnnoyingScalar sqrt(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:134
int i
Definition: BiCGSTAB_step_by_step.cpp:9
#define eigen_assert(x)
Definition: Macros.h:910
Vector3f p0
Definition: MatrixBase_all.cpp:2
float * p
Definition: Tutorial_Map_using.cpp:9
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: EigenSolver.h:82
void doComputeEigenvectors()
Definition: EigenSolver.h:441
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Definition: EigenSolver.h:92
static void check_template_parameters()
Definition: EigenSolver.h:301
Eigen::Index Index
Definition: EigenSolver.h:84
ComputationInfo m_info
Definition: EigenSolver.h:310
EIGEN_DEVICE_FUNC constexpr EIGEN_STRONG_INLINE Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:217
EIGEN_DEVICE_FUNC constexpr EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:294
const MatrixType & matrixU() const
Returns the orthogonal matrix in the Schur decomposition.
Definition: RealSchur.h:128
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealSchur.h:194
const MatrixType & matrixT() const
Returns the quasi-triangular matrix in the Schur decomposition.
Definition: RealSchur.h:144
RealSchur & compute(const EigenBase< InputType > &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
@ NumericalIssue
Definition: Constants.h:442
@ Success
Definition: Constants.h:440
#define isfinite(X)
Definition: main.h:111
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool() isfinite(const Eigen::bfloat16 &h)
Definition: BFloat16.h:752

References abs(), eigen_assert, i, Eigen::numext::isfinite(), isfinite, matrix(), Eigen::NumericalIssue, p, p0, sqrt(), and Eigen::Success.

◆ compute() [2/2]

template<typename MatrixType_ >
template<typename InputType >
EigenSolver& Eigen::EigenSolver< MatrixType_ >::compute ( const EigenBase< InputType > &  matrix,
bool  computeEigenvectors = true 
)

Computes eigendecomposition of given matrix.

Parameters
[in]matrixSquare matrix whose eigendecomposition is to be computed.
[in]computeEigenvectorsIf true, both the eigenvectors and the eigenvalues are computed; if false, only the eigenvalues are computed.
Returns
Reference to *this

This function computes the eigenvalues of the real matrix matrix. The eigenvalues() function can be used to retrieve them. If computeEigenvectors is true, then the eigenvectors are also computed and can be retrieved by calling eigenvectors().

The matrix is first reduced to real Schur form using the RealSchur class. The Schur decomposition is then used to compute the eigenvalues and eigenvectors.

The cost of the computation is dominated by the cost of the Schur decomposition, which is very approximately \( 25n^3 \) (where \( n \) is the size of the matrix) if computeEigenvectors is true, and \( 10n^3 \) if computeEigenvectors is false.

This method reuses of the allocated data in the EigenSolver object.

Example:

EigenSolver<MatrixXf> es;
MatrixXf A = MatrixXf::Random(4, 4);
es.compute(A, /* computeEigenvectors = */ false);
cout << "The eigenvalues of A are: " << es.eigenvalues().transpose() << endl;
es.compute(A + MatrixXf::Identity(4, 4), false); // re-use es to compute eigenvalues of A+I
cout << "The eigenvalues of A+I are: " << es.eigenvalues().transpose() << endl;

Output:

 

Referenced by __attribute__(), ctms_decompositions(), Eigen::EigenSolver< MatrixType_ >::EigenSolver(), eigensolver(), and eigensolver_verify_assert().

◆ doComputeEigenvectors()

template<typename MatrixType >
void Eigen::EigenSolver< MatrixType >::doComputeEigenvectors
private
441  {
442  using std::abs;
443  const Index size = m_eivec.cols();
445 
446  // inefficient! this is already computed in RealSchur
447  Scalar norm(0);
448  for (Index j = 0; j < size; ++j) {
449  norm += m_matT.row(j).segment((std::max)(j - 1, Index(0)), size - (std::max)(j - 1, Index(0))).cwiseAbs().sum();
450  }
451 
452  // Backsubstitute to find vectors of upper triangular form
453  if (norm == Scalar(0)) {
454  return;
455  }
456 
457  for (Index n = size - 1; n >= 0; n--) {
458  Scalar p = m_eivalues.coeff(n).real();
459  Scalar q = m_eivalues.coeff(n).imag();
460 
461  // Scalar vector
462  if (q == Scalar(0)) {
463  Scalar lastr(0), lastw(0);
464  Index l = n;
465 
466  m_matT.coeffRef(n, n) = Scalar(1);
467  for (Index i = n - 1; i >= 0; i--) {
468  Scalar w = m_matT.coeff(i, i) - p;
469  Scalar r = m_matT.row(i).segment(l, n - l + 1).dot(m_matT.col(n).segment(l, n - l + 1));
470 
471  if (m_eivalues.coeff(i).imag() < Scalar(0)) {
472  lastw = w;
473  lastr = r;
474  } else {
475  l = i;
476  if (m_eivalues.coeff(i).imag() == Scalar(0)) {
477  if (w != Scalar(0))
478  m_matT.coeffRef(i, n) = -r / w;
479  else
480  m_matT.coeffRef(i, n) = -r / (eps * norm);
481  } else // Solve real equations
482  {
483  Scalar x = m_matT.coeff(i, i + 1);
484  Scalar y = m_matT.coeff(i + 1, i);
485  Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) +
486  m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
487  Scalar t = (x * lastr - lastw * r) / denom;
488  m_matT.coeffRef(i, n) = t;
489  if (abs(x) > abs(lastw))
490  m_matT.coeffRef(i + 1, n) = (-r - w * t) / x;
491  else
492  m_matT.coeffRef(i + 1, n) = (-lastr - y * t) / lastw;
493  }
494 
495  // Overflow control
496  Scalar t = abs(m_matT.coeff(i, n));
497  if ((eps * t) * t > Scalar(1)) m_matT.col(n).tail(size - i) /= t;
498  }
499  }
500  } else if (q < Scalar(0) && n > 0) // Complex vector
501  {
502  Scalar lastra(0), lastsa(0), lastw(0);
503  Index l = n - 1;
504 
505  // Last vector component imaginary so matrix is triangular
506  if (abs(m_matT.coeff(n, n - 1)) > abs(m_matT.coeff(n - 1, n))) {
507  m_matT.coeffRef(n - 1, n - 1) = q / m_matT.coeff(n, n - 1);
508  m_matT.coeffRef(n - 1, n) = -(m_matT.coeff(n, n) - p) / m_matT.coeff(n, n - 1);
509  } else {
510  ComplexScalar cc =
511  ComplexScalar(Scalar(0), -m_matT.coeff(n - 1, n)) / ComplexScalar(m_matT.coeff(n - 1, n - 1) - p, q);
512  m_matT.coeffRef(n - 1, n - 1) = numext::real(cc);
513  m_matT.coeffRef(n - 1, n) = numext::imag(cc);
514  }
515  m_matT.coeffRef(n, n - 1) = Scalar(0);
516  m_matT.coeffRef(n, n) = Scalar(1);
517  for (Index i = n - 2; i >= 0; i--) {
518  Scalar ra = m_matT.row(i).segment(l, n - l + 1).dot(m_matT.col(n - 1).segment(l, n - l + 1));
519  Scalar sa = m_matT.row(i).segment(l, n - l + 1).dot(m_matT.col(n).segment(l, n - l + 1));
520  Scalar w = m_matT.coeff(i, i) - p;
521 
522  if (m_eivalues.coeff(i).imag() < Scalar(0)) {
523  lastw = w;
524  lastra = ra;
525  lastsa = sa;
526  } else {
527  l = i;
528  if (m_eivalues.coeff(i).imag() == RealScalar(0)) {
529  ComplexScalar cc = ComplexScalar(-ra, -sa) / ComplexScalar(w, q);
530  m_matT.coeffRef(i, n - 1) = numext::real(cc);
531  m_matT.coeffRef(i, n) = numext::imag(cc);
532  } else {
533  // Solve complex equations
534  Scalar x = m_matT.coeff(i, i + 1);
535  Scalar y = m_matT.coeff(i + 1, i);
536  Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) +
537  m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
538  Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
539  if ((vr == Scalar(0)) && (vi == Scalar(0)))
540  vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw));
541 
542  ComplexScalar cc = ComplexScalar(x * lastra - lastw * ra + q * sa, x * lastsa - lastw * sa - q * ra) /
543  ComplexScalar(vr, vi);
544  m_matT.coeffRef(i, n - 1) = numext::real(cc);
545  m_matT.coeffRef(i, n) = numext::imag(cc);
546  if (abs(x) > (abs(lastw) + abs(q))) {
547  m_matT.coeffRef(i + 1, n - 1) = (-ra - w * m_matT.coeff(i, n - 1) + q * m_matT.coeff(i, n)) / x;
548  m_matT.coeffRef(i + 1, n) = (-sa - w * m_matT.coeff(i, n) - q * m_matT.coeff(i, n - 1)) / x;
549  } else {
550  cc = ComplexScalar(-lastra - y * m_matT.coeff(i, n - 1), -lastsa - y * m_matT.coeff(i, n)) /
551  ComplexScalar(lastw, q);
552  m_matT.coeffRef(i + 1, n - 1) = numext::real(cc);
553  m_matT.coeffRef(i + 1, n) = numext::imag(cc);
554  }
555  }
556 
557  // Overflow control
558  Scalar t = numext::maxi<Scalar>(abs(m_matT.coeff(i, n - 1)), abs(m_matT.coeff(i, n)));
559  if ((eps * t) * t > Scalar(1)) m_matT.block(i, n - 1, size - i, 2) /= t;
560  }
561  }
562 
563  // We handled a pair of complex conjugate eigenvalues, so need to skip them both
564  n--;
565  } else {
566  eigen_assert(0 && "Internal bug in EigenSolver (INF or NaN has not been detected)"); // this should not happen
567  }
568  }
569 
570  // Back transformation to get eigenvectors of original matrix
571  for (Index j = size - 1; j >= 0; j--) {
572  m_tmp.noalias() = m_eivec.leftCols(j + 1) * m_matT.col(j).segment(0, j + 1);
573  m_eivec.col(j) = m_tmp;
574  }
575 }
AnnoyingScalar imag(const AnnoyingScalar &)
Definition: AnnoyingScalar.h:132
const unsigned n
Definition: CG3DPackingUnitTest.cpp:11
RowVector3d w
Definition: Matrix_resize_int.cpp:3
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:46
EIGEN_DEVICE_FUNC constexpr EIGEN_STRONG_INLINE const Scalar & coeff(Index rowId, Index colId) const
Definition: PlainObjectBase.h:198
float real
Definition: datatypes.h:10
#define max(a, b)
Definition: datatypes.h:23
Scalar * y
Definition: level1_cplx_impl.h:128
double eps
Definition: crbond_bessel.cc:24
EIGEN_DEVICE_FUNC const Scalar & q
Definition: SpecialFunctionsImpl.h:2019
r
Definition: UniformPSDSelfTest.py:20
double epsilon
Definition: osc_ring_sarah_asymptotics.h:43
list x
Definition: plotDoE.py:28
t
Definition: plotPSD.py:36
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2

References abs(), eigen_assert, CRBond_Bessel::eps, oomph::SarahBL::epsilon, i, imag(), j, max, n, p, Eigen::numext::q, UniformPSDSelfTest::r, size, plotPSD::t, w, plotDoE::x, and y.

◆ eigenvalues()

template<typename MatrixType_ >
const EigenvalueType& Eigen::EigenSolver< MatrixType_ >::eigenvalues ( ) const
inline

Returns the eigenvalues of given matrix.

Returns
A const reference to the column vector containing the eigenvalues.
Precondition
Either the constructor EigenSolver(const MatrixType&,bool) or the member function compute(const MatrixType&, bool) has been called before.

The eigenvalues are repeated according to their algebraic multiplicity, so there are as many eigenvalues as rows in the matrix. The eigenvalues are not sorted in any particular order.

Example:

MatrixXd ones = MatrixXd::Ones(3, 3);
EigenSolver<MatrixXd> es(ones, false);
cout << "The eigenvalues of the 3x3 matrix of ones are:" << endl << es.eigenvalues() << endl;
MatrixXcf ones
Definition: ComplexEigenSolver_eigenvalues.cpp:1

Output:

See also
eigenvectors(), pseudoEigenvalueMatrix(), MatrixBase::eigenvalues()
246  {
247  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
248  return m_eivalues;
249  }

References eigen_assert, Eigen::EigenSolver< MatrixType_ >::m_eivalues, and Eigen::EigenSolver< MatrixType_ >::m_isInitialized.

Referenced by eigensolver(), eigensolver_generic_extra(), eigensolver_verify_assert(), and Eigen::internal::eigenvalues_selector< Derived, false >::run().

◆ eigenvectors()

template<typename MatrixType >
EigenSolver< MatrixType >::EigenvectorsType Eigen::EigenSolver< MatrixType >::eigenvectors

Returns the eigenvectors of given matrix.

Returns
Matrix whose columns are the (possibly complex) eigenvectors.
Precondition
Either the constructor EigenSolver(const MatrixType&,bool) or the member function compute(const MatrixType&, bool) has been called before, and computeEigenvectors was set to true (the default).

Column \( k \) of the returned matrix is an eigenvector corresponding to eigenvalue number \( k \) as returned by eigenvalues(). The eigenvectors are normalized to have (Euclidean) norm equal to one. The matrix returned by this function is the matrix \( V \) in the eigendecomposition \( A = V D V^{-1} \), if it exists.

Example:

MatrixXd ones = MatrixXd::Ones(3, 3);
EigenSolver<MatrixXd> es(ones);
cout << "The first eigenvector of the 3x3 matrix of ones is:" << endl << es.eigenvectors().col(0) << endl;

Output:

See also
eigenvalues(), pseudoEigenvectors()
344  {
345  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
346  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
347  const RealScalar precision = RealScalar(2) * NumTraits<RealScalar>::epsilon();
348  Index n = m_eivec.cols();
349  EigenvectorsType matV(n, n);
350  for (Index j = 0; j < n; ++j) {
352  j + 1 == n) {
353  // we have a real eigen value
354  matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
355  matV.col(j).normalize();
356  } else {
357  // we have a pair of complex eigen values
358  for (Index i = 0; i < n; ++i) {
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 }
NumTraits< Scalar >::Real RealScalar
Definition: EigenSolver.h:83
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorsType
Type for matrix of eigenvectors as returned by eigenvectors().
Definition: EigenSolver.h:108
EIGEN_DEVICE_FUNC bool isMuchSmallerThan(const Scalar &x, const OtherScalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
Definition: MathFunctions.h:1916

References Eigen::Matrix< Scalar_, Rows_, Cols_, Options_, MaxRows_, MaxCols_ >::coeffRef(), eigen_assert, oomph::SarahBL::epsilon, i, imag(), Eigen::internal::isMuchSmallerThan(), j, and n.

Referenced by __attribute__(), eigensolver(), eigensolver_generic_extra(), and eigensolver_verify_assert().

◆ getMaxIterations()

template<typename MatrixType_ >
Index Eigen::EigenSolver< MatrixType_ >::getMaxIterations ( )
inline

Returns the maximum number of iterations.

295 { return m_realSchur.getMaxIterations(); }
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: RealSchur.h:210

References Eigen::RealSchur< MatrixType_ >::getMaxIterations(), and Eigen::EigenSolver< MatrixType_ >::m_realSchur.

Referenced by eigensolver().

◆ info()

template<typename MatrixType_ >
ComputationInfo Eigen::EigenSolver< MatrixType_ >::info ( ) const
inline
Returns
NumericalIssue if the input contains INF or NaN values or overflow occurred. Returns Success otherwise.
283  {
284  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
285  return m_info;
286  }

References eigen_assert, Eigen::EigenSolver< MatrixType_ >::m_info, and Eigen::EigenSolver< MatrixType_ >::m_isInitialized.

Referenced by eigensolver().

◆ pseudoEigenvalueMatrix()

template<typename MatrixType >
MatrixType Eigen::EigenSolver< MatrixType >::pseudoEigenvalueMatrix

Returns the block-diagonal matrix in the pseudo-eigendecomposition.

Returns
A block-diagonal matrix.
Precondition
Either the constructor EigenSolver(const MatrixType&,bool) or the member function compute(const MatrixType&, bool) has been called before.

The matrix \( D \) returned by this function is real and block-diagonal. The blocks on the diagonal are either 1-by-1 or 2-by-2 blocks of the form \( \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \). These blocks are not sorted in any particular order. The matrix \( D \) and the matrix \( V \) returned by pseudoEigenvectors() satisfy \( AV = VD \).

See also
pseudoEigenvectors() for an example, eigenvalues()
319  {
320  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
321  const RealScalar precision = RealScalar(2) * NumTraits<RealScalar>::epsilon();
322  const Index n = m_eivalues.rows();
323  MatrixType matD = MatrixType::Zero(n, n);
324  Index i = 0;
325  for (; i < n - 1; ++i) {
328  matD.coeffRef(i, i) = real;
329  if (!internal::isMuchSmallerThan(imag, real, precision)) {
330  matD.coeffRef(i, i + 1) = imag;
331  matD.coeffRef(i + 1, i) = -imag;
332  matD.coeffRef(i + 1, i + 1) = real;
333  ++i;
334  }
335  }
336  if (i == n - 1) {
337  matD.coeffRef(i, i) = numext::real(m_eivalues.coeff(i));
338  }
339 
340  return matD;
341 }
MatrixXf MatrixType
Definition: benchmark-blocking-sizes.cpp:52
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: PlainObjectBase.h:191
const AutoDiffScalar< DerType > & real(const AutoDiffScalar< DerType > &x)
Definition: AutoDiffScalar.h:486
DerType::Scalar imag(const AutoDiffScalar< DerType > &)
Definition: AutoDiffScalar.h:490
double Zero
Definition: pseudosolid_node_update_elements.cc:35
Definition: main.h:116
Definition: main.h:115

References eigen_assert, oomph::SarahBL::epsilon, i, imag(), Eigen::imag(), Eigen::internal::isMuchSmallerThan(), n, Eigen::real(), and oomph::PseudoSolidHelper::Zero.

Referenced by eigensolver(), eigensolver_generic_extra(), and eigensolver_verify_assert().

◆ pseudoEigenvectors()

template<typename MatrixType_ >
const MatrixType& Eigen::EigenSolver< MatrixType_ >::pseudoEigenvectors ( ) const
inline

Returns the pseudo-eigenvectors of given matrix.

Returns
Const reference to matrix whose columns are the pseudo-eigenvectors.
Precondition
Either the constructor EigenSolver(const MatrixType&,bool) or the member function compute(const MatrixType&, bool) has been called before, and computeEigenvectors was set to true (the default).

The real matrix \( V \) returned by this function and the block-diagonal matrix \( D \) returned by pseudoEigenvalueMatrix() satisfy \( AV = VD \).

Example:

MatrixXd A = MatrixXd::Random(6, 6);
cout << "Here is a random 6x6 matrix, A:" << endl << A << endl << endl;
EigenSolver<MatrixXd> es(A);
MatrixXd D = es.pseudoEigenvalueMatrix();
MatrixXd V = es.pseudoEigenvectors();
cout << "The pseudo-eigenvalue matrix D is:" << endl << D << endl;
cout << "The pseudo-eigenvector matrix V is:" << endl << V << endl;
cout << "Finally, V * D * V^(-1) = " << endl << V * D * V.inverse() << endl;

Output:

See also
pseudoEigenvalueMatrix(), eigenvectors()
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  }

References eigen_assert, Eigen::EigenSolver< MatrixType_ >::m_eigenvectorsOk, Eigen::EigenSolver< MatrixType_ >::m_eivec, and Eigen::EigenSolver< MatrixType_ >::m_isInitialized.

Referenced by eigensolver(), eigensolver_generic_extra(), and eigensolver_verify_assert().

◆ setMaxIterations()

template<typename MatrixType_ >
EigenSolver& Eigen::EigenSolver< MatrixType_ >::setMaxIterations ( Index  maxIters)
inline

Sets the maximum number of iterations allowed.

289  {
290  m_realSchur.setMaxIterations(maxIters);
291  return *this;
292  }
RealSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: RealSchur.h:204

References Eigen::EigenSolver< MatrixType_ >::m_realSchur, and Eigen::RealSchur< MatrixType_ >::setMaxIterations().

Referenced by eigensolver().

Member Data Documentation

◆ m_eigenvectorsOk

template<typename MatrixType_ >
bool Eigen::EigenSolver< MatrixType_ >::m_eigenvectorsOk
protected

◆ m_eivalues

template<typename MatrixType_ >
EigenvalueType Eigen::EigenSolver< MatrixType_ >::m_eivalues
protected

◆ m_eivec

template<typename MatrixType_ >
MatrixType Eigen::EigenSolver< MatrixType_ >::m_eivec
protected

◆ m_info

template<typename MatrixType_ >
ComputationInfo Eigen::EigenSolver< MatrixType_ >::m_info
protected

◆ m_isInitialized

◆ m_matT

template<typename MatrixType_ >
MatrixType Eigen::EigenSolver< MatrixType_ >::m_matT
protected

◆ m_realSchur

template<typename MatrixType_ >
RealSchur<MatrixType> Eigen::EigenSolver< MatrixType_ >::m_realSchur
protected

◆ m_tmp

template<typename MatrixType_ >
ColumnVectorType Eigen::EigenSolver< MatrixType_ >::m_tmp
protected

The documentation for this class was generated from the following file: