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;
24 template<
typename MatrixType,
typename DiagType,
typename SubDiagType>
26 ComputationInfo 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)
181 compute(matrix.
derived(), options);
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.");
372 static const int m_maxIterations = 30;
375 static EIGEN_DEVICE_FUNC
376 void check_template_parameters()
378 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
381 EigenvectorsType m_eivec;
386 bool m_isInitialized;
387 bool m_eigenvectorsOk;
411 template<
int StorageOrder,
typename RealScalar,
typename Scalar,
typename Index>
416 template<
typename MatrixType>
417 template<
typename InputType>
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;
452 mat = matrix.template triangularView<Lower>();
455 mat.template triangularView<Lower>() /= scale;
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;
470 template<
typename MatrixType>
479 if (computeEigenvectors)
483 m_info = internal::computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
485 m_isInitialized =
true;
486 m_eigenvectorsOk = computeEigenvectors;
502 template<
typename MatrixType,
typename DiagType,
typename SubDiagType>
504 ComputationInfo 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));
575 template<
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); }
582 template<
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;
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;
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;
740 template<
typename SolverType>
741 struct direct_selfadjoint_eigenvalues<SolverType,2,false>
743 typedef typename SolverType::MatrixType MatrixType;
744 typedef typename SolverType::RealVectorType VectorType;
745 typedef typename SolverType::Scalar
Scalar;
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;
825 template<
typename MatrixType>
830 internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>::run(*
this,matrix,options);
837 template<
int StorageOrder,
typename RealScalar,
typename Scalar,
typename Index>
851 mu -= numext::abs(e);
856 mu -= e / ((td + (td>
RealScalar(0) ? h : -h)) / e);
858 mu -= e2 / (td + (td>
RealScalar(0) ? h : -h));
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];
895 q.applyOnTheRight(k,k+1,rot);
904 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H MatrixType::Scalar Scalar
Scalar type for matrices of type _MatrixType.
Definition: SelfAdjointEigenSolver.h:89
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
const RealVectorType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: SelfAdjointEigenSolver.h:300
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:76
Namespace containing all symbols from the Eigen library.
Definition: Core:141
Rotation given by a cosine-sine pair.
Definition: Jacobi.h:34
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
Tridiagonal decomposition of a selfadjoint matrix.
Definition: Tridiagonalization.h:64
Derived & setIdentity()
Definition: CwiseNullaryOp.h:873
Derived & derived()
Definition: EigenBase.h:46
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
Definition: EigenBase.h:29
Eigen::Index Index
Definition: SelfAdjointEigenSolver.h:90
const EigenvectorsType & eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: SelfAdjointEigenSolver.h:277
Definition: Constants.h:405
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:74
SelfAdjointEigenSolver & computeFromTridiagonal(const RealVectorType &diag, const SubDiagonalType &subdiag, int options=ComputeEigenvectors)
Computes the eigen decomposition from a tridiagonal symmetric matrix.
Definition: SelfAdjointEigenSolver.h:472
internal::plain_col_type< MatrixType, RealScalar >::type RealVectorType
Type for vector of eigenvalues as returned by eigenvalues().
Definition: SelfAdjointEigenSolver.h:109
NumTraits< Scalar >::Real RealScalar
Real scalar type for _MatrixType.
Definition: SelfAdjointEigenSolver.h:100
Definition: Constants.h:442
SelfAdjointEigenSolver(const EigenBase< InputType > &matrix, int options=ComputeEigenvectors)
Constructor; computes eigendecomposition of given matrix.
Definition: SelfAdjointEigenSolver.h:173
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:281
SelfAdjointEigenSolver(Index size)
Constructor, pre-allocates memory for dynamic-size matrices.
Definition: SelfAdjointEigenSolver.h:147
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:324
SelfAdjointEigenSolver & computeDirect(const MatrixType &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix using a closed-form algorithm.
Definition: SelfAdjointEigenSolver.h:828
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.
ComputationInfo
Definition: Constants.h:440
Definition: Constants.h:446
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sin_op< typename Derived::Scalar >, const Derived > sin(const Eigen::ArrayBase< Derived > &x)
void makeGivens(const Scalar &p, const Scalar &q, Scalar *r=0)
Definition: Jacobi.h:162
MatrixType operatorInverseSqrt() const
Computes the inverse square root of the matrix.
Definition: SelfAdjointEigenSolver.h:349