11 #ifndef EIGEN_SELFADJOINTEIGENSOLVER_H
12 #define EIGEN_SELFADJOINTEIGENSOLVER_H
14 #include "./Tridiagonalization.h"
18 template<
typename _MatrixType>
19 class GeneralizedSelfAdjointEigenSolver;
22 template<
typename SolverType,
int Size,
bool IsComplex>
struct direct_selfadjoint_eigenvalues;
23 template<
typename MatrixType,
typename DiagType,
typename SubDiagType>
24 ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag,
const Index maxIterations,
bool computeEigenvectors, MatrixType& eivec);
74 typedef _MatrixType MatrixType;
76 Size = MatrixType::RowsAtCompileTime,
77 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
78 Options = MatrixType::Options,
79 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
83 typedef typename MatrixType::Scalar
Scalar;
103 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type
RealVectorType;
122 m_isInitialized(false)
139 : m_eivec(size, size),
141 m_subdiag(size > 1 ? size - 1 : 1),
142 m_isInitialized(false)
162 : m_eivec(matrix.rows(), matrix.cols()),
163 m_eivalues(matrix.cols()),
164 m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
165 m_isInitialized(false)
258 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
259 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
281 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
306 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
307 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
308 return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
332 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
333 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
334 return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
344 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
356 static void check_template_parameters()
358 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
361 EigenvectorsType m_eivec;
363 typename TridiagonalizationType::SubDiagonalType m_subdiag;
365 bool m_isInitialized;
366 bool m_eigenvectorsOk;
386 template<
int StorageOrder,
typename RealScalar,
typename Scalar,
typename Index>
388 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
391 template<
typename MatrixType>
396 check_template_parameters();
399 eigen_assert(matrix.cols() == matrix.rows());
400 eigen_assert((options&~(EigVecMask|GenEigMask))==0
401 && (options&EigVecMask)!=EigVecMask
402 &&
"invalid option parameter");
404 Index n = matrix.cols();
405 m_eivalues.resize(n,1);
409 m_eivalues.coeffRef(0,0) = numext::real(matrix.coeff(0,0));
410 if(computeEigenvectors)
411 m_eivec.setOnes(n,n);
413 m_isInitialized =
true;
414 m_eigenvectorsOk = computeEigenvectors;
423 mat = matrix.template triangularView<Lower>();
426 mat.template triangularView<Lower>() /= scale;
428 internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
430 m_info = internal::computeFromTridiagonal_impl(diag, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
435 m_isInitialized =
true;
436 m_eigenvectorsOk = computeEigenvectors;
440 template<
typename MatrixType>
449 if (computeEigenvectors)
453 m_info = computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
455 m_isInitialized =
true;
456 m_eigenvectorsOk = computeEigenvectors;
471 template<
typename MatrixType,
typename DiagType,
typename SubDiagType>
472 ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag,
const Index maxIterations,
bool computeEigenvectors, MatrixType& eivec)
477 typedef typename MatrixType::Scalar Scalar;
479 Index n = diag.size();
484 typedef typename DiagType::RealScalar RealScalar;
485 const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
489 for (Index i = start; i<end; ++i)
490 if (internal::isMuchSmallerThan(abs(subdiag[i]),(abs(diag[i])+abs(diag[i+1]))) || abs(subdiag[i]) <= considerAsZero)
494 while (end>0 && subdiag[end-1]==0)
503 if(iter > maxIterations * n)
break;
506 while (start>0 && subdiag[start-1]!=0)
509 internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), subdiag.data(), start, end, computeEigenvectors ? eivec.data() : (Scalar*)0, n);
511 if (iter <= maxIterations * n)
521 for (Index i = 0; i < n-1; ++i)
524 diag.segment(i,n-i).minCoeff(&k);
527 std::swap(diag[i], diag[k+i]);
528 if(computeEigenvectors)
529 eivec.col(i).swap(eivec.col(k+i));
536 template<
typename SolverType,
int Size,
bool IsComplex>
struct direct_selfadjoint_eigenvalues
539 static inline void run(SolverType& eig,
const typename SolverType::MatrixType& A,
int options)
540 { eig.compute(A,options); }
543 template<
typename SolverType>
struct direct_selfadjoint_eigenvalues<SolverType,3,false>
545 typedef typename SolverType::MatrixType MatrixType;
546 typedef typename SolverType::RealVectorType VectorType;
547 typedef typename SolverType::Scalar Scalar;
548 typedef typename SolverType::EigenvectorsType EigenvectorsType;
556 static inline void computeRoots(
const MatrixType& m, VectorType& roots)
558 EIGEN_USING_STD_MATH(sqrt)
559 EIGEN_USING_STD_MATH(atan2)
560 EIGEN_USING_STD_MATH(cos)
561 EIGEN_USING_STD_MATH(sin)
562 const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0);
563 const Scalar s_sqrt3 = sqrt(Scalar(3.0));
568 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);
569 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);
570 Scalar c2 = m(0,0) + m(1,1) + m(2,2);
574 Scalar c2_over_3 = c2*s_inv3;
575 Scalar a_over_3 = (c2*c2_over_3 - c1)*s_inv3;
576 a_over_3 = numext::maxi(a_over_3, Scalar(0));
578 Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1));
580 Scalar q = a_over_3*a_over_3*a_over_3 - half_b*half_b;
581 q = numext::maxi(q, Scalar(0));
584 Scalar rho = sqrt(a_over_3);
585 Scalar theta = atan2(sqrt(q),half_b)*s_inv3;
586 Scalar cos_theta = cos(theta);
587 Scalar sin_theta = sin(theta);
589 roots(0) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta);
590 roots(1) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta);
591 roots(2) = c2_over_3 + Scalar(2)*rho*cos_theta;
595 static inline
bool extract_kernel(MatrixType& mat, Ref<VectorType> res, Ref<VectorType> representative)
600 mat.diagonal().cwiseAbs().maxCoeff(&i0);
603 representative = mat.col(i0);
606 n0 = (c0 = representative.cross(mat.col((i0+1)%3))).squaredNorm();
607 n1 = (c1 = representative.cross(mat.col((i0+2)%3))).squaredNorm();
608 if(n0>n1) res = c0/std::sqrt(n0);
609 else res = c1/std::sqrt(n1);
615 static inline void run(SolverType& solver,
const MatrixType& mat,
int options)
617 eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows());
618 eigen_assert((options&~(EigVecMask|GenEigMask))==0
619 && (options&EigVecMask)!=EigVecMask
620 &&
"invalid option parameter");
623 EigenvectorsType& eivecs = solver.m_eivec;
624 VectorType& eivals = solver.m_eivalues;
627 Scalar shift = mat.trace() / Scalar(3);
629 MatrixType scaledMat = mat.template selfadjointView<Lower>();
630 scaledMat.diagonal().array() -= shift;
631 Scalar scale = scaledMat.cwiseAbs().maxCoeff();
632 if(scale > 0) scaledMat /= scale;
635 computeRoots(scaledMat,eivals);
638 if(computeEigenvectors)
643 eivecs.setIdentity();
651 Scalar d0 = eivals(2) - eivals(1);
652 Scalar d1 = eivals(1) - eivals(0);
662 tmp.diagonal().array () -= eivals(k);
664 extract_kernel(tmp, eivecs.col(k), eivecs.col(l));
672 eivecs.col(l) -= eivecs.col(k).dot(eivecs.col(l))*eivecs.col(l);
673 eivecs.col(l).normalize();
678 tmp.diagonal().array () -= eivals(l);
681 extract_kernel(tmp, eivecs.col(l), dummy);
685 eivecs.col(1) = eivecs.col(2).cross(eivecs.col(0)).normalized();
691 eivals.array() += shift;
694 solver.m_isInitialized =
true;
695 solver.m_eigenvectorsOk = computeEigenvectors;
700 template<
typename SolverType>
701 struct direct_selfadjoint_eigenvalues<SolverType,2,false>
703 typedef typename SolverType::MatrixType MatrixType;
704 typedef typename SolverType::RealVectorType VectorType;
705 typedef typename SolverType::Scalar Scalar;
706 typedef typename SolverType::EigenvectorsType EigenvectorsType;
709 static inline void computeRoots(
const MatrixType& m, VectorType& roots)
712 const Scalar t0 = Scalar(0.5) * sqrt( numext::abs2(m(0,0)-m(1,1)) + Scalar(4)*numext::abs2(m(1,0)));
713 const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1));
719 static inline void run(SolverType& solver,
const MatrixType& mat,
int options)
721 EIGEN_USING_STD_MATH(sqrt);
722 EIGEN_USING_STD_MATH(abs);
724 eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
725 eigen_assert((options&~(EigVecMask|GenEigMask))==0
726 && (options&EigVecMask)!=EigVecMask
727 &&
"invalid option parameter");
730 EigenvectorsType& eivecs = solver.m_eivec;
731 VectorType& eivals = solver.m_eivalues;
734 Scalar scale = mat.cwiseAbs().maxCoeff();
735 scale = numext::maxi(scale,Scalar(1));
736 MatrixType scaledMat = mat / scale;
739 computeRoots(scaledMat,eivals);
742 if(computeEigenvectors)
746 eivecs.setIdentity();
750 scaledMat.diagonal().array () -= eivals(1);
751 Scalar a2 = numext::abs2(scaledMat(0,0));
752 Scalar c2 = numext::abs2(scaledMat(1,1));
753 Scalar b2 = numext::abs2(scaledMat(1,0));
756 eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0);
757 eivecs.col(1) /= sqrt(a2+b2);
761 eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0);
762 eivecs.col(1) /= sqrt(c2+b2);
765 eivecs.col(0) << eivecs.col(1).unitOrthogonal();
773 solver.m_isInitialized =
true;
774 solver.m_eigenvectorsOk = computeEigenvectors;
780 template<
typename MatrixType>
785 internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>::run(*
this,matrix,options);
790 template<
int StorageOrder,
typename RealScalar,
typename Scalar,
typename Index>
792 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
795 RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
796 RealScalar e = subdiag[end-1];
802 RealScalar mu = diag[end];
807 RealScalar e2 = numext::abs2(subdiag[end-1]);
808 RealScalar h = numext::hypot(td,e);
809 if(e2==0) mu -= (e / (td + (td>0 ? 1 : -1))) * (e / h);
810 else mu -= e2 / (td + (td>0 ? h : -h));
813 RealScalar x = diag[start] - mu;
814 RealScalar z = subdiag[start];
815 for (Index k = start; k < end; ++k)
817 JacobiRotation<RealScalar> rot;
818 rot.makeGivens(x, z);
821 RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
822 RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
824 diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
825 diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
826 subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
830 subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
836 z = -rot.s() * subdiag[k+1];
837 subdiag[k + 1] = rot.c() * subdiag[k+1];
844 Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n);
845 q.applyOnTheRight(k,k+1,rot);
854 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H
MatrixType operatorSqrt() const
Computes the positive-definite square root of the matrix.
Definition: SelfAdjointEigenSolver.h:304
static const int m_maxIterations
Maximum number of iterations.
Definition: SelfAdjointEigenSolver.h:353
MatrixType::Scalar Scalar
Scalar type for matrices of type _MatrixType.
Definition: SelfAdjointEigenSolver.h:83
Computes eigenvalues and eigenvectors of selfadjoint matrices.
Definition: SelfAdjointEigenSolver.h:70
SelfAdjointEigenSolver(const MatrixType &matrix, int options=ComputeEigenvectors)
Constructor; computes eigendecomposition of given matrix.
Definition: SelfAdjointEigenSolver.h:161
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
Tridiagonal decomposition of a selfadjoint matrix.
Definition: Tridiagonalization.h:63
Derived & setIdentity()
Definition: CwiseNullaryOp.h:779
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:252
Eigen::Index Index
Definition: SelfAdjointEigenSolver.h:84
Definition: Constants.h:387
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SelfAdjointEigenSolver.h:342
SelfAdjointEigenSolver & computeFromTridiagonal(const RealVectorType &diag, const SubDiagonalType &subdiag, int options=ComputeEigenvectors)
Computes the eigen decomposition from a tridiagonal symmetric matrix.
Definition: SelfAdjointEigenSolver.h:442
MatrixType operatorInverseSqrt() const
Computes the inverse square root of the matrix.
Definition: SelfAdjointEigenSolver.h:330
internal::plain_col_type< MatrixType, RealScalar >::type RealVectorType
Type for vector of eigenvalues as returned by eigenvalues().
Definition: SelfAdjointEigenSolver.h:103
NumTraits< Scalar >::Real RealScalar
Real scalar type for _MatrixType.
Definition: SelfAdjointEigenSolver.h:94
Definition: Constants.h:424
SelfAdjointEigenSolver(Index size)
Constructor, pre-allocates memory for dynamic-size matrices.
Definition: SelfAdjointEigenSolver.h:138
Definition: Eigen_Colamd.h:54
SelfAdjointEigenSolver & compute(const MatrixType &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix.
Definition: SelfAdjointEigenSolver.h:394
SelfAdjointEigenSolver & computeDirect(const MatrixType &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix using a closed-form algorithm.
Definition: SelfAdjointEigenSolver.h:783
ComputationInfo
Definition: Constants.h:422
Definition: Constants.h:428
const RealVectorType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: SelfAdjointEigenSolver.h:279
const EigenvectorsType & eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: SelfAdjointEigenSolver.h:256