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;
72 typedef _MatrixType MatrixType;
74 Size = MatrixType::RowsAtCompileTime,
75 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
76 Options = MatrixType::Options,
77 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
81 typedef typename MatrixType::Scalar
Scalar;
82 typedef typename MatrixType::Index Index;
99 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type
RealVectorType;
116 m_isInitialized(false)
132 : m_eivec(size, size),
134 m_subdiag(size > 1 ? size - 1 : 1),
135 m_isInitialized(false)
154 : m_eivec(matrix.rows(), matrix.cols()),
155 m_eivalues(matrix.cols()),
156 m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
157 m_isInitialized(false)
159 compute(matrix, options);
230 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
231 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
252 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
276 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
277 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
278 return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
301 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
302 eigen_assert(m_eigenvectorsOk &&
"The eigenvectors have not been computed together with the eigenvalues.");
303 return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
312 eigen_assert(m_isInitialized &&
"SelfAdjointEigenSolver is not initialized.");
321 static const int m_maxIterations = 30;
323 #ifdef EIGEN2_SUPPORT 325 : m_eivec(matrix.rows(), matrix.cols()),
326 m_eivalues(matrix.cols()),
327 m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
328 m_isInitialized(false)
330 compute(matrix, computeEigenvectors);
334 : m_eivec(matA.cols(), matA.cols()),
335 m_eivalues(matA.cols()),
336 m_subdiag(matA.cols() > 1 ? matA.cols() - 1 : 1),
337 m_isInitialized(
false)
342 void compute(
const MatrixType& matrix,
bool computeEigenvectors)
347 void compute(
const MatrixType& matA,
const MatrixType& matB,
bool computeEigenvectors =
true)
351 #endif // EIGEN2_SUPPORT 354 static void check_template_parameters()
356 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
363 bool m_isInitialized;
364 bool m_eigenvectorsOk;
384 template<
int StorageOrder,
typename RealScalar,
typename Scalar,
typename Index>
388 template<
typename MatrixType>
392 check_template_parameters();
395 eigen_assert(matrix.cols() == matrix.rows());
396 eigen_assert((options&~(EigVecMask|GenEigMask))==0
397 && (options&EigVecMask)!=EigVecMask
398 &&
"invalid option parameter");
400 Index n = matrix.cols();
401 m_eivalues.resize(n,1);
405 m_eivalues.coeffRef(0,0) = numext::real(matrix.coeff(0,0));
406 if(computeEigenvectors)
407 m_eivec.setOnes(n,n);
409 m_isInitialized =
true;
410 m_eigenvectorsOk = computeEigenvectors;
416 MatrixType& mat = m_eivec;
419 mat = matrix.template triangularView<Lower>();
422 mat.template triangularView<Lower>() /= scale;
423 m_subdiag.resize(n-1);
424 internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
432 for (Index i = start; i<end; ++i)
433 if (internal::isMuchSmallerThan(abs(m_subdiag[i]),(abs(diag[i])+abs(diag[i+1]))))
437 while (end>0 && m_subdiag[end-1]==0)
446 if(iter > m_maxIterations * n)
break;
449 while (start>0 && m_subdiag[start-1]!=0)
452 internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (
Scalar*)0, n);
455 if (iter <= m_maxIterations * n)
465 for (Index i = 0; i < n-1; ++i)
468 m_eivalues.segment(i,n-i).minCoeff(&k);
471 std::swap(m_eivalues[i], m_eivalues[k+i]);
472 if(computeEigenvectors)
473 m_eivec.col(i).swap(m_eivec.col(k+i));
481 m_isInitialized =
true;
482 m_eigenvectorsOk = computeEigenvectors;
489 template<
typename SolverType,
int Size,
bool IsComplex>
struct direct_selfadjoint_eigenvalues
491 static inline void run(SolverType& eig,
const typename SolverType::MatrixType& A,
int options)
492 { eig.compute(A,options); }
495 template<
typename SolverType>
struct direct_selfadjoint_eigenvalues<SolverType,3,false>
497 typedef typename SolverType::MatrixType MatrixType;
498 typedef typename SolverType::RealVectorType VectorType;
499 typedef typename SolverType::Scalar
Scalar;
500 typedef typename MatrixType::Index Index;
506 static inline void computeRoots(
const MatrixType& m, VectorType& roots)
512 const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0);
513 const Scalar s_sqrt3 = sqrt(Scalar(3.0));
518 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);
519 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);
520 Scalar c2 = m(0,0) + m(1,1) + m(2,2);
524 Scalar c2_over_3 = c2*s_inv3;
525 Scalar a_over_3 = (c2*c2_over_3 - c1)*s_inv3;
526 if(a_over_3<Scalar(0))
527 a_over_3 = Scalar(0);
529 Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1));
531 Scalar q = a_over_3*a_over_3*a_over_3 - half_b*half_b;
536 Scalar rho = sqrt(a_over_3);
537 Scalar theta = atan2(sqrt(q),half_b)*s_inv3;
538 Scalar cos_theta = cos(theta);
539 Scalar sin_theta = sin(theta);
541 roots(0) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta);
542 roots(1) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta);
543 roots(2) = c2_over_3 + Scalar(2)*rho*cos_theta;
551 mat.diagonal().cwiseAbs().maxCoeff(&i0);
554 representative = mat.col(i0);
557 n0 = (c0 = representative.cross(mat.col((i0+1)%3))).squaredNorm();
558 n1 = (c1 = representative.cross(mat.col((i0+2)%3))).squaredNorm();
559 if(n0>n1) res = c0/std::sqrt(n0);
560 else res = c1/std::sqrt(n1);
565 static inline void run(SolverType& solver,
const MatrixType& mat,
int options)
567 eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows());
568 eigen_assert((options&~(EigVecMask|GenEigMask))==0
569 && (options&EigVecMask)!=EigVecMask
570 &&
"invalid option parameter");
573 MatrixType& eivecs = solver.m_eivec;
574 VectorType& eivals = solver.m_eivalues;
577 Scalar shift = mat.trace() / Scalar(3);
579 MatrixType scaledMat = mat.template selfadjointView<Lower>();
580 scaledMat.diagonal().array() -= shift;
581 Scalar scale = scaledMat.cwiseAbs().maxCoeff();
582 if(scale > 0) scaledMat /= scale;
585 computeRoots(scaledMat,eivals);
588 if(computeEigenvectors)
593 eivecs.setIdentity();
601 Scalar d0 = eivals(2) - eivals(1);
602 Scalar d1 = eivals(1) - eivals(0);
612 tmp.diagonal().array () -= eivals(k);
614 extract_kernel(tmp, eivecs.col(k), eivecs.col(l));
622 eivecs.col(l) -= eivecs.col(k).dot(eivecs.col(l))*eivecs.col(l);
623 eivecs.col(l).normalize();
628 tmp.diagonal().array () -= eivals(l);
631 extract_kernel(tmp, eivecs.col(l), dummy);
635 eivecs.col(1) = eivecs.col(2).cross(eivecs.col(0)).normalized();
641 eivals.array() += shift;
644 solver.m_isInitialized =
true;
645 solver.m_eigenvectorsOk = computeEigenvectors;
650 template<
typename SolverType>
struct direct_selfadjoint_eigenvalues<SolverType,2,false>
652 typedef typename SolverType::MatrixType MatrixType;
653 typedef typename SolverType::RealVectorType VectorType;
654 typedef typename SolverType::Scalar
Scalar;
656 static inline void computeRoots(
const MatrixType& m, VectorType& roots)
659 const Scalar t0 = Scalar(0.5) * sqrt( numext::abs2(m(0,0)-m(1,1)) + Scalar(4)*numext::abs2(m(1,0)));
660 const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1));
665 static inline void run(SolverType& solver,
const MatrixType& mat,
int options)
670 eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
671 eigen_assert((options&~(EigVecMask|GenEigMask))==0
672 && (options&EigVecMask)!=EigVecMask
673 &&
"invalid option parameter");
676 MatrixType& eivecs = solver.m_eivec;
677 VectorType& eivals = solver.m_eivalues;
680 Scalar scale = mat.cwiseAbs().maxCoeff();
681 scale = (std::max)(scale,Scalar(1));
682 MatrixType scaledMat = mat / scale;
685 computeRoots(scaledMat,eivals);
688 if(computeEigenvectors)
692 eivecs.setIdentity();
696 scaledMat.diagonal().array () -= eivals(1);
697 Scalar a2 = numext::abs2(scaledMat(0,0));
698 Scalar c2 = numext::abs2(scaledMat(1,1));
699 Scalar b2 = numext::abs2(scaledMat(1,0));
702 eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0);
703 eivecs.col(1) /= sqrt(a2+b2);
707 eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0);
708 eivecs.col(1) /= sqrt(c2+b2);
711 eivecs.col(0) << eivecs.col(1).unitOrthogonal();
719 solver.m_isInitialized =
true;
720 solver.m_eigenvectorsOk = computeEigenvectors;
726 template<
typename MatrixType>
730 internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>::run(*
this,matrix,options);
735 template<
int StorageOrder,
typename RealScalar,
typename Scalar,
typename Index>
753 if(e2==0) mu -= (e / (td + (td>0 ? 1 : -1))) * (e / h);
754 else mu -= e2 / (td + (td>0 ? h : -h));
759 for (Index k = start; k < end; ++k)
765 RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
766 RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
768 diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
769 diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
770 subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
774 subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
780 z = -rot.s() * subdiag[k+1];
781 subdiag[k + 1] = rot.c() * subdiag[k+1];
789 q.applyOnTheRight(k,k+1,rot);
798 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SelfAdjointEigenSolver.h:310
MatrixType::Scalar Scalar
Scalar type for matrices of type _MatrixType.
Definition: SelfAdjointEigenSolver.h:81
MatrixType operatorInverseSqrt() const
Computes the inverse square root of the matrix.
Definition: SelfAdjointEigenSolver.h:299
Definition: Constants.h:339
const RealVectorType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: SelfAdjointEigenSolver.h:250
SelfAdjointEigenSolver(const MatrixType &matrix, int options=ComputeEigenvectors)
Constructor; computes eigendecomposition of given matrix.
Definition: SelfAdjointEigenSolver.h:153
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:104
SelfAdjointEigenSolver(Index size)
Constructor, pre-allocates memory for dynamic-size matrices.
Definition: SelfAdjointEigenSolver.h:131
Computes eigenvalues and eigenvectors of selfadjoint matrices.
Definition: SelfAdjointEigenSolver.h:68
Rotation given by a cosine-sine pair.
Definition: ForwardDeclarations.h:228
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
Tridiagonal decomposition of a selfadjoint matrix.
Definition: Tridiagonalization.h:61
SelfAdjointEigenSolver & computeDirect(const MatrixType &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix using a direct algorithm.
Definition: SelfAdjointEigenSolver.h:728
void makeGivens(const Scalar &p, const Scalar &q, Scalar *z=0)
Definition: Jacobi.h:148
Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem.
Definition: GeneralizedSelfAdjointEigenSolver.h:48
MatrixType operatorSqrt() const
Computes the positive-definite square root of the matrix.
Definition: SelfAdjointEigenSolver.h:274
NumTraits< Scalar >::Real RealScalar
Real scalar type for _MatrixType.
Definition: SelfAdjointEigenSolver.h:90
internal::plain_col_type< MatrixType, RealScalar >::type RealVectorType
Type for vector of eigenvalues as returned by eigenvalues().
Definition: SelfAdjointEigenSolver.h:99
A matrix or vector expression mapping an existing expressions.
Definition: Ref.h:17
Definition: Eigen_Colamd.h:54
const MatrixType & eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: SelfAdjointEigenSolver.h:228
Definition: Constants.h:380
Definition: Constants.h:376
SelfAdjointEigenSolver & compute(const MatrixType &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix.
Definition: SelfAdjointEigenSolver.h:390
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:127
ComputationInfo
Definition: Constants.h:374
Definition: Constants.h:336