11 #ifndef EIGEN_EIGENSOLVER_H 12 #define EIGEN_EIGENSOLVER_H 14 #include "./RealSchur.h" 72 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
73 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
74 Options = MatrixType::Options,
75 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
76 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
80 typedef typename MatrixType::Scalar
Scalar;
113 EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
122 : m_eivec(size, size),
124 m_isInitialized(false),
125 m_eigenvectorsOk(false),
146 template<
typename InputType>
148 : m_eivec(matrix.rows(), matrix.cols()),
149 m_eivalues(matrix.cols()),
150 m_isInitialized(false),
151 m_eigenvectorsOk(false),
152 m_realSchur(matrix.cols()),
153 m_matT(matrix.rows(), matrix.cols()),
201 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
202 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
246 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
277 template<
typename InputType>
283 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
301 void doComputeEigenvectors();
305 static void check_template_parameters()
307 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
312 EigenvalueType m_eivalues;
313 bool m_isInitialized;
314 bool m_eigenvectorsOk;
320 ColumnVectorType m_tmp;
323 template<
typename MatrixType>
326 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
327 Index n = m_eivalues.rows();
329 for (
Index i=0; i<n; ++i)
331 if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i))))
332 matD.coeffRef(i,i) = numext::real(m_eivalues.coeff(i));
335 matD.template block<2,2>(i,i) << numext::real(m_eivalues.coeff(i)), numext::imag(m_eivalues.coeff(i)),
336 -numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i));
343 template<
typename MatrixType>
346 eigen_assert(m_isInitialized &&
"EigenSolver is not initialized.");
347 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
348 Index n = m_eivec.cols();
350 for (
Index j=0; j<n; ++j)
352 if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j))) || j+1==n)
355 matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
356 matV.col(j).normalize();
361 for (
Index i=0; i<n; ++i)
363 matV.coeffRef(i,j) =
ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1));
364 matV.coeffRef(i,j+1) =
ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
366 matV.col(j).normalize();
367 matV.col(j+1).normalize();
374 template<
typename MatrixType>
375 template<
typename InputType>
379 check_template_parameters();
383 using numext::isfinite;
384 eigen_assert(matrix.
cols() == matrix.
rows());
389 m_info = m_realSchur.
info();
393 m_matT = m_realSchur.
matrixT();
394 if (computeEigenvectors)
395 m_eivec = m_realSchur.
matrixU();
400 while (i < matrix.
cols())
402 if (i == matrix.
cols() - 1 || m_matT.coeff(i+1, i) ==
Scalar(0))
404 m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
405 if(!(isfinite)(m_eivalues.coeffRef(i)))
407 m_isInitialized =
true;
408 m_eigenvectorsOk =
false;
416 Scalar p =
Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
421 Scalar t0 = m_matT.coeff(i+1, i);
422 Scalar t1 = m_matT.coeff(i, i+1);
423 Scalar maxval = numext::maxi<Scalar>(abs(p),numext::maxi<Scalar>(abs(t0),abs(t1)));
427 z = maxval * sqrt(abs(p0 * p0 + t0 * t1));
430 m_eivalues.coeffRef(i) =
ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
431 m_eivalues.coeffRef(i+1) =
ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
432 if(!((isfinite)(m_eivalues.coeffRef(i)) && (isfinite)(m_eivalues.coeffRef(i+1))))
434 m_isInitialized =
true;
435 m_eigenvectorsOk =
false;
444 if (computeEigenvectors)
445 doComputeEigenvectors();
448 m_isInitialized =
true;
449 m_eigenvectorsOk = computeEigenvectors;
455 template<
typename Scalar>
460 if (abs(yr) > abs(yi))
464 return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
470 return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
475 template<
typename MatrixType>
479 const Index size = m_eivec.cols();
484 for (
Index j = 0; j < size; ++j)
486 norm += m_matT.row(j).segment((std::max)(j-1,
Index(0)), size-(std::max)(j-1,
Index(0))).cwiseAbs().sum();
495 for (
Index n = size-1; n >= 0; n--)
497 Scalar p = m_eivalues.coeff(n).real();
498 Scalar q = m_eivalues.coeff(n).imag();
503 Scalar lastr(0), lastw(0);
506 m_matT.coeffRef(n,n) = 1.0;
507 for (
Index i = n-1; i >= 0; i--)
509 Scalar w = m_matT.coeff(i,i) - p;
510 Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
512 if (m_eivalues.coeff(i).imag() <
Scalar(0))
520 if (m_eivalues.coeff(i).imag() ==
Scalar(0))
523 m_matT.coeffRef(i,n) = -r / w;
525 m_matT.coeffRef(i,n) = -r / (eps * norm);
529 Scalar x = m_matT.coeff(i,i+1);
530 Scalar y = m_matT.coeff(i+1,i);
531 Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
532 Scalar t = (x * lastr - lastw * r) / denom;
533 m_matT.coeffRef(i,n) = t;
534 if (abs(x) > abs(lastw))
535 m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
537 m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
541 Scalar t = abs(m_matT.coeff(i,n));
542 if ((eps * t) * t >
Scalar(1))
543 m_matT.col(n).tail(size-i) /= t;
547 else if (q <
Scalar(0) && n > 0)
549 Scalar lastra(0), lastsa(0), lastw(0);
553 if (abs(m_matT.coeff(n,n-1)) > abs(m_matT.coeff(n-1,n)))
555 m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
556 m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
560 std::complex<Scalar> cc = cdiv<Scalar>(
Scalar(0),-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
561 m_matT.coeffRef(n-1,n-1) = numext::real(cc);
562 m_matT.coeffRef(n-1,n) = numext::imag(cc);
564 m_matT.coeffRef(n,n-1) =
Scalar(0);
565 m_matT.coeffRef(n,n) =
Scalar(1);
566 for (
Index i = n-2; i >= 0; i--)
568 Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
569 Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
570 Scalar w = m_matT.coeff(i,i) - p;
572 if (m_eivalues.coeff(i).imag() <
Scalar(0))
581 if (m_eivalues.coeff(i).imag() == RealScalar(0))
583 std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
584 m_matT.coeffRef(i,n-1) = numext::real(cc);
585 m_matT.coeffRef(i,n) = numext::imag(cc);
590 Scalar x = m_matT.coeff(i,i+1);
591 Scalar y = m_matT.coeff(i+1,i);
592 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;
593 Scalar vi = (m_eivalues.coeff(i).real() - p) *
Scalar(2) * q;
595 vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw));
597 std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
598 m_matT.coeffRef(i,n-1) = numext::real(cc);
599 m_matT.coeffRef(i,n) = numext::imag(cc);
600 if (abs(x) > (abs(lastw) + abs(q)))
602 m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
603 m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
607 cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q);
608 m_matT.coeffRef(i+1,n-1) = numext::real(cc);
609 m_matT.coeffRef(i+1,n) = numext::imag(cc);
614 Scalar t = numext::maxi<Scalar>(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n)));
615 if ((eps * t) * t >
Scalar(1))
616 m_matT.block(i, n-1, size-i, 2) /= t;
626 eigen_assert(0 &&
"Internal bug in EigenSolver (INF or NaN has not been detected)");
631 for (
Index j = size-1; j >= 0; j--)
633 m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
634 m_eivec.col(j) = m_tmp;
640 #endif // EIGEN_EIGENSOLVER_H Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorsType
Type for matrix of eigenvectors as returned by eigenvectors().
Definition: EigenSolver.h:104
EigenSolver & compute(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Computes eigendecomposition of given matrix.
EigenvectorsType eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: EigenSolver.h:344
EigenSolver(Index size)
Default constructor with memory preallocation.
Definition: EigenSolver.h:121
const MatrixType & matrixU() const
Returns the orthogonal matrix in the Schur decomposition.
Definition: RealSchur.h:127
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: EigenSolver.h:295
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealSchur.h:195
Eigen::Index Index
Definition: EigenSolver.h:82
Derived & derived()
Definition: EigenBase.h:44
ComputationInfo info() const
Definition: EigenSolver.h:281
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:252
Index rows() const
Definition: EigenBase.h:58
const MatrixType & matrixT() const
Returns the quasi-triangular matrix in the Schur decomposition.
Definition: RealSchur.h:144
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: RealSchur.h:213
Definition: EigenBase.h:28
EigenSolver(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Constructor; computes eigendecomposition of given matrix.
Definition: EigenSolver.h:147
MatrixType pseudoEigenvalueMatrix() const
Returns the block-diagonal matrix in the pseudo-eigendecomposition.
Definition: EigenSolver.h:324
const MatrixType & pseudoEigenvectors() const
Returns the pseudo-eigenvectors of given matrix.
Definition: EigenSolver.h:199
Definition: Constants.h:434
RealSchur & compute(const EigenBase< InputType > &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
EigenSolver & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: EigenSolver.h:288
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Definition: EigenSolver.h:90
Definition: Constants.h:432
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
Definition: EigenSolver.h:69
EigenSolver()
Default constructor.
Definition: EigenSolver.h:113
Index cols() const
Definition: EigenBase.h:61
RealSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: RealSchur.h:206
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: EigenSolver.h:80
Computes eigenvalues and eigenvectors of general matrices.
Definition: EigenSolver.h:64
ComputationInfo
Definition: Constants.h:430
const EigenvalueType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: EigenSolver.h:244
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > EigenvalueType
Type for vector of eigenvalues as returned by eigenvalues().
Definition: EigenSolver.h:97