11#ifndef EIGEN_SELFADJOINTEIGENSOLVER_H
12#define EIGEN_SELFADJOINTEIGENSOLVER_H
14#include "./Tridiagonalization.h"
18template<
typename _MatrixType>
19class GeneralizedSelfAdjointEigenSolver;
22template<
typename SolverType,
int Size,
bool IsComplex>
struct direct_selfadjoint_eigenvalues;
24template<
typename MatrixType,
typename DiagType,
typename SubDiagType>
26ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag,
const Index maxIterations,
bool computeEigenvectors, MatrixType& eivec);
80 typedef _MatrixType MatrixType;
82 Size = MatrixType::RowsAtCompileTime,
83 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
84 Options = MatrixType::Options,
85 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
89 typedef typename MatrixType::Scalar
Scalar;
109 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type
RealVectorType;
129 m_info(InvalidInput),
130 m_isInitialized(false),
131 m_eigenvectorsOk(false)
148 : m_eivec(size, size),
150 m_subdiag(size > 1 ? size - 1 : 1),
151 m_hcoeffs(size > 1 ? size - 1 : 1),
152 m_isInitialized(false),
153 m_eigenvectorsOk(false)
171 template<
typename InputType>
174 : m_eivec(matrix.rows(), matrix.cols()),
175 m_eivalues(matrix.cols()),
176 m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
177 m_hcoeffs(matrix.cols() > 1 ? matrix.cols() - 1 : 1),
178 m_isInitialized(false),
179 m_eigenvectorsOk(false)
214 template<
typename InputType>
279 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
280 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
302 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
326 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
327 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
328 return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
351 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
352 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
353 return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
363 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
375 static EIGEN_DEVICE_FUNC
376 void check_template_parameters()
378 EIGEN_STATIC_ASSERT_NON_INTEGER(
Scalar);
381 EigenvectorsType m_eivec;
383 typename TridiagonalizationType::SubDiagonalType m_subdiag;
384 typename TridiagonalizationType::CoeffVectorType m_hcoeffs;
386 bool m_isInitialized;
387 bool m_eigenvectorsOk;
411template<
int StorageOrder,
typename RealScalar,
typename Scalar,
typename Index>
413static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag,
Index start,
Index end, Scalar* matrixQ,
Index n);
416template<
typename MatrixType>
417template<
typename InputType>
419SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
420::compute(
const EigenBase<InputType>& a_matrix,
int options)
422 check_template_parameters();
424 const InputType &matrix(a_matrix.derived());
426 EIGEN_USING_STD(abs);
427 eigen_assert(matrix.cols() == matrix.rows());
428 eigen_assert((options&~(EigVecMask|GenEigMask))==0
429 && (options&EigVecMask)!=EigVecMask
430 &&
"invalid option parameter");
432 Index n = matrix.cols();
433 m_eivalues.resize(n,1);
438 m_eivalues.coeffRef(0,0) = numext::real(m_eivec.coeff(0,0));
439 if(computeEigenvectors)
440 m_eivec.setOnes(n,n);
442 m_isInitialized =
true;
443 m_eigenvectorsOk = computeEigenvectors;
448 RealVectorType& diag = m_eivalues;
449 EigenvectorsType& mat = m_eivec;
452 mat = matrix.template triangularView<Lower>();
453 RealScalar scale = mat.cwiseAbs().maxCoeff();
454 if(scale==RealScalar(0)) scale = RealScalar(1);
455 mat.template triangularView<Lower>() /= scale;
456 m_subdiag.resize(n-1);
457 m_hcoeffs.resize(n-1);
458 internal::tridiagonalization_inplace(mat, diag, m_subdiag, m_hcoeffs, computeEigenvectors);
460 m_info = internal::computeFromTridiagonal_impl(diag, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
465 m_isInitialized =
true;
466 m_eigenvectorsOk = computeEigenvectors;
470template<
typename MatrixType>
471SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
479 if (computeEigenvectors)
481 m_eivec.setIdentity(diag.size(), diag.size());
483 m_info = internal::computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
485 m_isInitialized =
true;
486 m_eigenvectorsOk = computeEigenvectors;
502template<
typename MatrixType,
typename DiagType,
typename SubDiagType>
504ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag,
const Index maxIterations,
bool computeEigenvectors, MatrixType& eivec)
507 typedef typename MatrixType::Scalar Scalar;
509 Index n = diag.size();
514 typedef typename DiagType::RealScalar RealScalar;
515 const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
519 for (
Index i = start; i<end; ++i) {
520 if (numext::abs(subdiag[i]) < considerAsZero) {
521 subdiag[i] = RealScalar(0);
525 const RealScalar scaled_subdiag = precision_inv * subdiag[i];
526 if (scaled_subdiag * scaled_subdiag <= (numext::abs(diag[i])+numext::abs(diag[i+1]))) {
527 subdiag[i] = RealScalar(0);
533 while (end>0 && subdiag[end-1]==RealScalar(0))
542 if(iter > maxIterations * n)
break;
545 while (start>0 && subdiag[start-1]!=0)
548 internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), subdiag.data(), start, end, computeEigenvectors ? eivec.data() : (Scalar*)0, n);
550 if (iter <= maxIterations * n)
560 for (
Index i = 0; i < n-1; ++i)
563 diag.segment(i,n-i).minCoeff(&k);
566 numext::swap(diag[i], diag[k+i]);
567 if(computeEigenvectors)
568 eivec.col(i).swap(eivec.col(k+i));
575template<
typename SolverType,
int Size,
bool IsComplex>
struct direct_selfadjoint_eigenvalues
578 static inline void run(SolverType& eig,
const typename SolverType::MatrixType& A,
int options)
579 { eig.compute(A,options); }
582template<
typename SolverType>
struct direct_selfadjoint_eigenvalues<SolverType,3,false>
584 typedef typename SolverType::MatrixType MatrixType;
585 typedef typename SolverType::RealVectorType VectorType;
586 typedef typename SolverType::Scalar Scalar;
587 typedef typename SolverType::EigenvectorsType EigenvectorsType;
595 static inline void computeRoots(
const MatrixType& m, VectorType& roots)
597 EIGEN_USING_STD(sqrt)
598 EIGEN_USING_STD(atan2)
601 const Scalar s_inv3 = Scalar(1)/Scalar(3);
602 const Scalar s_sqrt3 = sqrt(Scalar(3));
607 Scalar c0 = m(0,0)*m(1,1)*m(2,2) + Scalar(2)*m(1,0)*m(2,0)*m(2,1) - m(0,0)*m(2,1)*m(2,1) - m(1,1)*m(2,0)*m(2,0) - m(2,2)*m(1,0)*m(1,0);
608 Scalar c1 = m(0,0)*m(1,1) - m(1,0)*m(1,0) + m(0,0)*m(2,2) - m(2,0)*m(2,0) + m(1,1)*m(2,2) - m(2,1)*m(2,1);
609 Scalar c2 = m(0,0) + m(1,1) + m(2,2);
613 Scalar c2_over_3 = c2*s_inv3;
614 Scalar a_over_3 = (c2*c2_over_3 - c1)*s_inv3;
615 a_over_3 = numext::maxi(a_over_3, Scalar(0));
617 Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1));
619 Scalar q = a_over_3*a_over_3*a_over_3 - half_b*half_b;
620 q = numext::maxi(q, Scalar(0));
623 Scalar rho = sqrt(a_over_3);
624 Scalar theta = atan2(sqrt(q),half_b)*s_inv3;
625 Scalar cos_theta = cos(theta);
626 Scalar sin_theta = sin(theta);
628 roots(0) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta);
629 roots(1) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta);
630 roots(2) = c2_over_3 + Scalar(2)*rho*cos_theta;
634 static inline bool extract_kernel(MatrixType& mat, Ref<VectorType> res, Ref<VectorType> representative)
636 EIGEN_USING_STD(abs);
637 EIGEN_USING_STD(sqrt);
640 mat.diagonal().cwiseAbs().maxCoeff(&i0);
643 representative = mat.col(i0);
646 n0 = (c0 = representative.cross(mat.col((i0+1)%3))).squaredNorm();
647 n1 = (c1 = representative.cross(mat.col((i0+2)%3))).squaredNorm();
648 if(n0>n1) res = c0/sqrt(n0);
649 else res = c1/sqrt(n1);
655 static inline void run(SolverType& solver,
const MatrixType& mat,
int options)
657 eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows());
658 eigen_assert((options&~(EigVecMask|GenEigMask))==0
659 && (options&EigVecMask)!=EigVecMask
660 &&
"invalid option parameter");
663 EigenvectorsType& eivecs = solver.m_eivec;
664 VectorType& eivals = solver.m_eivalues;
667 Scalar shift = mat.trace() / Scalar(3);
669 MatrixType scaledMat = mat.template selfadjointView<Lower>();
670 scaledMat.diagonal().array() -= shift;
671 Scalar scale = scaledMat.cwiseAbs().maxCoeff();
672 if(scale > 0) scaledMat /= scale;
675 computeRoots(scaledMat,eivals);
678 if(computeEigenvectors)
683 eivecs.setIdentity();
691 Scalar d0 = eivals(2) - eivals(1);
692 Scalar d1 = eivals(1) - eivals(0);
702 tmp.diagonal().array () -= eivals(k);
704 extract_kernel(tmp, eivecs.col(k), eivecs.col(l));
712 eivecs.col(l) -= eivecs.col(k).dot(eivecs.col(l))*eivecs.col(l);
713 eivecs.col(l).normalize();
718 tmp.diagonal().array () -= eivals(l);
721 extract_kernel(tmp, eivecs.col(l), dummy);
725 eivecs.col(1) = eivecs.col(2).cross(eivecs.col(0)).normalized();
731 eivals.array() += shift;
734 solver.m_isInitialized =
true;
735 solver.m_eigenvectorsOk = computeEigenvectors;
740template<
typename SolverType>
741struct direct_selfadjoint_eigenvalues<SolverType,2,false>
743 typedef typename SolverType::MatrixType MatrixType;
744 typedef typename SolverType::RealVectorType VectorType;
745 typedef typename SolverType::Scalar Scalar;
746 typedef typename SolverType::EigenvectorsType EigenvectorsType;
749 static inline void computeRoots(
const MatrixType& m, VectorType& roots)
751 EIGEN_USING_STD(sqrt);
752 const Scalar t0 = Scalar(0.5) * sqrt( numext::abs2(m(0,0)-m(1,1)) + Scalar(4)*numext::abs2(m(1,0)));
753 const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1));
759 static inline void run(SolverType& solver,
const MatrixType& mat,
int options)
761 EIGEN_USING_STD(sqrt);
762 EIGEN_USING_STD(abs);
764 eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
765 eigen_assert((options&~(EigVecMask|GenEigMask))==0
766 && (options&EigVecMask)!=EigVecMask
767 &&
"invalid option parameter");
770 EigenvectorsType& eivecs = solver.m_eivec;
771 VectorType& eivals = solver.m_eivalues;
774 Scalar shift = mat.trace() / Scalar(2);
775 MatrixType scaledMat = mat;
776 scaledMat.coeffRef(0,1) = mat.coeff(1,0);
777 scaledMat.diagonal().array() -= shift;
778 Scalar scale = scaledMat.cwiseAbs().maxCoeff();
779 if(scale > Scalar(0))
783 computeRoots(scaledMat,eivals);
786 if(computeEigenvectors)
790 eivecs.setIdentity();
794 scaledMat.diagonal().array () -= eivals(1);
795 Scalar a2 = numext::abs2(scaledMat(0,0));
796 Scalar c2 = numext::abs2(scaledMat(1,1));
797 Scalar b2 = numext::abs2(scaledMat(1,0));
800 eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0);
801 eivecs.col(1) /= sqrt(a2+b2);
805 eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0);
806 eivecs.col(1) /= sqrt(c2+b2);
809 eivecs.col(0) << eivecs.col(1).unitOrthogonal();
815 eivals.array() += shift;
818 solver.m_isInitialized =
true;
819 solver.m_eigenvectorsOk = computeEigenvectors;
825template<
typename MatrixType>
827SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
828::computeDirect(
const MatrixType& matrix,
int options)
830 internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>::run(*
this,matrix,options);
837template<
int StorageOrder,
typename RealScalar,
typename Scalar,
typename Index>
839static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag,
Index start,
Index end, Scalar* matrixQ,
Index n)
842 RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
843 RealScalar e = subdiag[end-1];
849 RealScalar mu = diag[end];
850 if(td==RealScalar(0)) {
851 mu -= numext::abs(e);
852 }
else if (e != RealScalar(0)) {
853 const RealScalar e2 = numext::abs2(e);
854 const RealScalar h = numext::hypot(td,e);
855 if(e2 == RealScalar(0)) {
856 mu -= e / ((td + (td>RealScalar(0) ? h : -h)) / e);
858 mu -= e2 / (td + (td>RealScalar(0) ? h : -h));
862 RealScalar x = diag[start] - mu;
863 RealScalar z = subdiag[start];
866 for (
Index k = start; k < end && z != RealScalar(0); ++k)
868 JacobiRotation<RealScalar> rot;
869 rot.makeGivens(x, z);
872 RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
873 RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
875 diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
876 diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
877 subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
880 subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
886 z = -rot.s() * subdiag[k+1];
887 subdiag[k + 1] = rot.c() * subdiag[k+1];
894 Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n);
895 q.applyOnTheRight(k,k+1,rot);
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Computes eigenvalues and eigenvectors of selfadjoint matrices.
Definition: SelfAdjointEigenSolver.h:77
MatrixType::Scalar Scalar
Scalar type for matrices of type _MatrixType.
Definition: SelfAdjointEigenSolver.h:89
SelfAdjointEigenSolver & computeFromTridiagonal(const RealVectorType &diag, const SubDiagonalType &subdiag, int options=ComputeEigenvectors)
Computes the eigen decomposition from a tridiagonal symmetric matrix.
Definition: SelfAdjointEigenSolver.h:472
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SelfAdjointEigenSolver.h:361
SelfAdjointEigenSolver & compute(const EigenBase< InputType > &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix.
const RealVectorType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: SelfAdjointEigenSolver.h:300
NumTraits< Scalar >::Real RealScalar
Real scalar type for _MatrixType.
Definition: SelfAdjointEigenSolver.h:100
const EigenvectorsType & eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: SelfAdjointEigenSolver.h:277
SelfAdjointEigenSolver(const EigenBase< InputType > &matrix, int options=ComputeEigenvectors)
Constructor; computes eigendecomposition of given matrix.
Definition: SelfAdjointEigenSolver.h:173
internal::plain_col_type< MatrixType, RealScalar >::type RealVectorType
Type for vector of eigenvalues as returned by eigenvalues().
Definition: SelfAdjointEigenSolver.h:109
SelfAdjointEigenSolver(Index size)
Constructor, pre-allocates memory for dynamic-size matrices.
Definition: SelfAdjointEigenSolver.h:147
Eigen::Index Index
Definition: SelfAdjointEigenSolver.h:90
static const int m_maxIterations
Maximum number of iterations.
Definition: SelfAdjointEigenSolver.h:372
MatrixType operatorInverseSqrt() const
Computes the inverse square root of the matrix.
Definition: SelfAdjointEigenSolver.h:349
MatrixType operatorSqrt() const
Computes the positive-definite square root of the matrix.
Definition: SelfAdjointEigenSolver.h:324
SelfAdjointEigenSolver & computeDirect(const MatrixType &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix using a closed-form algorithm.
Definition: SelfAdjointEigenSolver.h:828
Tridiagonal decomposition of a selfadjoint matrix.
Definition: Tridiagonalization.h:65
ComputationInfo
Definition: Constants.h:440
@ Success
Definition: Constants.h:442
@ NoConvergence
Definition: Constants.h:446
@ ComputeEigenvectors
Definition: Constants.h:405
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Definition: EigenBase.h:30
Derived & derived()
Definition: EigenBase.h:46
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:233