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)
160 template<
typename InputType>
163 : m_eivec(matrix.rows(), matrix.cols()),
164 m_eivalues(matrix.cols()),
165 m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
166 m_isInitialized(false)
168 compute(matrix.
derived(), options);
201 template<
typename InputType>
261 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
262 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
284 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
308 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
309 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
310 return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
333 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
334 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
335 return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
345 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
354 static const int m_maxIterations = 30;
357 static void check_template_parameters()
359 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
362 EigenvectorsType m_eivec;
366 bool m_isInitialized;
367 bool m_eigenvectorsOk;
391 template<
int StorageOrder,
typename RealScalar,
typename Scalar,
typename Index>
396 template<
typename MatrixType>
397 template<
typename InputType>
402 check_template_parameters();
404 const InputType &matrix(a_matrix.
derived());
407 eigen_assert(matrix.cols() == matrix.rows());
408 eigen_assert((options&~(EigVecMask|GenEigMask))==0
409 && (options&EigVecMask)!=EigVecMask
410 &&
"invalid option parameter");
412 Index n = matrix.cols();
413 m_eivalues.resize(n,1);
418 m_eivalues.coeffRef(0,0) = numext::real(m_eivec.coeff(0,0));
419 if(computeEigenvectors)
420 m_eivec.setOnes(n,n);
422 m_isInitialized =
true;
423 m_eigenvectorsOk = computeEigenvectors;
432 mat = matrix.template triangularView<Lower>();
435 mat.template triangularView<Lower>() /= scale;
437 internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
439 m_info = internal::computeFromTridiagonal_impl(diag, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
444 m_isInitialized =
true;
445 m_eigenvectorsOk = computeEigenvectors;
449 template<
typename MatrixType>
458 if (computeEigenvectors)
462 m_info = internal::computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
464 m_isInitialized =
true;
465 m_eigenvectorsOk = computeEigenvectors;
481 template<
typename MatrixType,
typename DiagType,
typename SubDiagType>
482 ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag,
const Index maxIterations,
bool computeEigenvectors, MatrixType& eivec)
487 typedef typename MatrixType::Scalar
Scalar;
489 Index n = diag.size();
494 typedef typename DiagType::RealScalar
RealScalar;
495 const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
500 for (
Index i = start; i<end; ++i)
501 if (internal::isMuchSmallerThan(
abs(subdiag[i]),(
abs(diag[i])+
abs(diag[i+1])),precision) ||
abs(subdiag[i]) <= considerAsZero)
505 while (end>0 && subdiag[end-1]==RealScalar(0))
514 if(iter > maxIterations * n)
break;
517 while (start>0 && subdiag[start-1]!=0)
520 internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), subdiag.data(), start, end, computeEigenvectors ? eivec.data() : (Scalar*)0, n);
522 if (iter <= maxIterations * n)
532 for (
Index i = 0; i < n-1; ++i)
535 diag.segment(i,n-i).minCoeff(&k);
538 std::swap(diag[i], diag[k+i]);
539 if(computeEigenvectors)
540 eivec.col(i).swap(eivec.col(k+i));
547 template<
typename SolverType,
int Size,
bool IsComplex>
struct direct_selfadjoint_eigenvalues
550 static inline void run(SolverType& eig,
const typename SolverType::MatrixType& A,
int options)
551 { eig.compute(A,options); }
554 template<
typename SolverType>
struct direct_selfadjoint_eigenvalues<SolverType,3,false>
556 typedef typename SolverType::MatrixType MatrixType;
557 typedef typename SolverType::RealVectorType VectorType;
558 typedef typename SolverType::Scalar
Scalar;
567 static inline void computeRoots(
const MatrixType& m, VectorType& roots)
569 EIGEN_USING_STD_MATH(
sqrt)
570 EIGEN_USING_STD_MATH(atan2)
571 EIGEN_USING_STD_MATH(
cos)
572 EIGEN_USING_STD_MATH(
sin)
573 const Scalar s_inv3 = Scalar(1)/Scalar(3);
574 const Scalar s_sqrt3 =
sqrt(Scalar(3));
579 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);
580 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);
581 Scalar c2 = m(0,0) + m(1,1) + m(2,2);
585 Scalar c2_over_3 = c2*s_inv3;
586 Scalar a_over_3 = (c2*c2_over_3 - c1)*s_inv3;
587 a_over_3 = numext::maxi(a_over_3, Scalar(0));
589 Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1));
591 Scalar q = a_over_3*a_over_3*a_over_3 - half_b*half_b;
592 q = numext::maxi(q, Scalar(0));
595 Scalar rho =
sqrt(a_over_3);
596 Scalar theta = atan2(
sqrt(q),half_b)*s_inv3;
597 Scalar cos_theta =
cos(theta);
598 Scalar sin_theta =
sin(theta);
600 roots(0) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta);
601 roots(1) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta);
602 roots(2) = c2_over_3 + Scalar(2)*rho*cos_theta;
611 mat.diagonal().cwiseAbs().maxCoeff(&i0);
614 representative = mat.col(i0);
617 n0 = (c0 = representative.cross(mat.col((i0+1)%3))).squaredNorm();
618 n1 = (c1 = representative.cross(mat.col((i0+2)%3))).squaredNorm();
619 if(n0>n1) res = c0/std::sqrt(n0);
620 else res = c1/std::sqrt(n1);
626 static inline void run(SolverType& solver,
const MatrixType& mat,
int options)
628 eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows());
629 eigen_assert((options&~(EigVecMask|GenEigMask))==0
630 && (options&EigVecMask)!=EigVecMask
631 &&
"invalid option parameter");
634 EigenvectorsType& eivecs = solver.m_eivec;
635 VectorType& eivals = solver.m_eivalues;
638 Scalar shift = mat.trace() / Scalar(3);
640 MatrixType scaledMat = mat.template selfadjointView<Lower>();
641 scaledMat.diagonal().array() -= shift;
642 Scalar scale = scaledMat.cwiseAbs().maxCoeff();
643 if(scale > 0) scaledMat /= scale;
646 computeRoots(scaledMat,eivals);
649 if(computeEigenvectors)
654 eivecs.setIdentity();
662 Scalar d0 = eivals(2) - eivals(1);
663 Scalar d1 = eivals(1) - eivals(0);
673 tmp.diagonal().array () -= eivals(k);
675 extract_kernel(tmp, eivecs.col(k), eivecs.col(l));
683 eivecs.col(l) -= eivecs.col(k).dot(eivecs.col(l))*eivecs.col(l);
684 eivecs.col(l).normalize();
689 tmp.diagonal().array () -= eivals(l);
692 extract_kernel(tmp, eivecs.col(l), dummy);
696 eivecs.col(1) = eivecs.col(2).cross(eivecs.col(0)).normalized();
702 eivals.array() += shift;
705 solver.m_isInitialized =
true;
706 solver.m_eigenvectorsOk = computeEigenvectors;
711 template<
typename SolverType>
712 struct direct_selfadjoint_eigenvalues<SolverType,2,false>
714 typedef typename SolverType::MatrixType MatrixType;
715 typedef typename SolverType::RealVectorType VectorType;
716 typedef typename SolverType::Scalar
Scalar;
720 static inline void computeRoots(
const MatrixType& m, VectorType& roots)
723 const Scalar t0 = Scalar(0.5) *
sqrt( numext::abs2(m(0,0)-m(1,1)) + Scalar(4)*numext::abs2(m(1,0)));
724 const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1));
730 static inline void run(SolverType& solver,
const MatrixType& mat,
int options)
732 EIGEN_USING_STD_MATH(
sqrt);
733 EIGEN_USING_STD_MATH(
abs);
735 eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
736 eigen_assert((options&~(EigVecMask|GenEigMask))==0
737 && (options&EigVecMask)!=EigVecMask
738 &&
"invalid option parameter");
741 EigenvectorsType& eivecs = solver.m_eivec;
742 VectorType& eivals = solver.m_eivalues;
745 Scalar shift = mat.trace() / Scalar(2);
746 MatrixType scaledMat = mat;
747 scaledMat.coeffRef(0,1) = mat.coeff(1,0);
748 scaledMat.diagonal().array() -= shift;
749 Scalar scale = scaledMat.cwiseAbs().maxCoeff();
750 if(scale > Scalar(0))
754 computeRoots(scaledMat,eivals);
757 if(computeEigenvectors)
761 eivecs.setIdentity();
765 scaledMat.diagonal().array () -= eivals(1);
766 Scalar a2 = numext::abs2(scaledMat(0,0));
767 Scalar c2 = numext::abs2(scaledMat(1,1));
768 Scalar b2 = numext::abs2(scaledMat(1,0));
771 eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0);
772 eivecs.col(1) /=
sqrt(a2+b2);
776 eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0);
777 eivecs.col(1) /=
sqrt(c2+b2);
780 eivecs.col(0) << eivecs.col(1).unitOrthogonal();
786 eivals.array() += shift;
789 solver.m_isInitialized =
true;
790 solver.m_eigenvectorsOk = computeEigenvectors;
796 template<
typename MatrixType>
801 internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>::run(*
this,matrix,options);
806 template<
int StorageOrder,
typename RealScalar,
typename Scalar,
typename Index>
826 else mu -= e2 / (td + (td>RealScalar(0) ? h : -h));
831 for (
Index k = start; k < end; ++k)
837 RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
838 RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
840 diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
841 diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
842 subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
846 subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
852 z = -rot.s() * subdiag[k+1];
853 subdiag[k + 1] = rot.c() * subdiag[k+1];
861 q.applyOnTheRight(k,k+1,rot);
870 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H MatrixType::Scalar Scalar
Scalar type for matrices of type _MatrixType.
Definition: SelfAdjointEigenSolver.h:83
void makeGivens(const Scalar &p, const Scalar &q, Scalar *z=0)
Definition: Jacobi.h:149
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:88
const RealVectorType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: SelfAdjointEigenSolver.h:282
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Computes eigenvalues and eigenvectors of selfadjoint matrices.
Definition: SelfAdjointEigenSolver.h:70
Namespace containing all symbols from the Eigen library.
Definition: Core:287
Rotation given by a cosine-sine pair.
Definition: ForwardDeclarations.h:263
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
Tridiagonal decomposition of a selfadjoint matrix.
Definition: Tridiagonalization.h:63
Derived & setIdentity()
Definition: CwiseNullaryOp.h:774
Derived & derived()
Definition: EigenBase.h:45
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:279
Definition: EigenBase.h:29
Eigen::Index Index
Definition: SelfAdjointEigenSolver.h:84
const EigenvectorsType & eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: SelfAdjointEigenSolver.h:259
Definition: Constants.h:395
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_cos_op< typename Derived::Scalar >, const Derived > cos(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
SelfAdjointEigenSolver & computeFromTridiagonal(const RealVectorType &diag, const SubDiagonalType &subdiag, int options=ComputeEigenvectors)
Computes the eigen decomposition from a tridiagonal symmetric matrix.
Definition: SelfAdjointEigenSolver.h:451
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:432
SelfAdjointEigenSolver(const EigenBase< InputType > &matrix, int options=ComputeEigenvectors)
Constructor; computes eigendecomposition of given matrix.
Definition: SelfAdjointEigenSolver.h:162
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:190
SelfAdjointEigenSolver(Index size)
Constructor, pre-allocates memory for dynamic-size matrices.
Definition: SelfAdjointEigenSolver.h:138
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: Eigen_Colamd.h:50
MatrixType operatorSqrt() const
Computes the positive-definite square root of the matrix.
Definition: SelfAdjointEigenSolver.h:306
SelfAdjointEigenSolver & computeDirect(const MatrixType &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix using a closed-form algorithm.
Definition: SelfAdjointEigenSolver.h:799
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SelfAdjointEigenSolver.h:343
SelfAdjointEigenSolver & compute(const EigenBase< InputType > &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix.
ComputationInfo
Definition: Constants.h:430
Definition: Constants.h:436
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sin_op< typename Derived::Scalar >, const Derived > sin(const Eigen::ArrayBase< Derived > &x)
MatrixType operatorInverseSqrt() const
Computes the inverse square root of the matrix.
Definition: SelfAdjointEigenSolver.h:331