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)
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.");
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
356 typename TridiagonalizationType::SubDiagonalType m_subdiag;
358 bool m_isInitialized;
359 bool m_eigenvectorsOk;
379 template<
int StorageOrder,
typename RealScalar,
typename Scalar,
typename Index>
380 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
383 template<
typename MatrixType>
387 eigen_assert(matrix.cols() == matrix.rows());
388 eigen_assert((options&~(EigVecMask|GenEigMask))==0
389 && (options&EigVecMask)!=EigVecMask
390 &&
"invalid option parameter");
392 Index n = matrix.cols();
393 m_eivalues.resize(n,1);
397 m_eivalues.coeffRef(0,0) = internal::real(matrix.coeff(0,0));
398 if(computeEigenvectors)
399 m_eivec.setOnes(n,n);
401 m_isInitialized =
true;
402 m_eigenvectorsOk = computeEigenvectors;
408 MatrixType& mat = m_eivec;
411 RealScalar scale = matrix.cwiseAbs().maxCoeff();
413 mat = matrix / scale;
414 m_subdiag.resize(n-1);
415 internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
423 for (Index i = start; i<end; ++i)
424 if (internal::isMuchSmallerThan(internal::abs(m_subdiag[i]),(internal::abs(diag[i])+internal::abs(diag[i+1]))))
428 while (end>0 && m_subdiag[end-1]==0)
437 if(iter > m_maxIterations * n)
break;
440 while (start>0 && m_subdiag[start-1]!=0)
443 internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (
Scalar*)0, n);
446 if (iter <= m_maxIterations * n)
456 for (Index i = 0; i < n-1; ++i)
459 m_eivalues.segment(i,n-i).minCoeff(&k);
462 std::swap(m_eivalues[i], m_eivalues[k+i]);
463 if(computeEigenvectors)
464 m_eivec.col(i).swap(m_eivec.col(k+i));
472 m_isInitialized =
true;
473 m_eigenvectorsOk = computeEigenvectors;
480 template<
typename SolverType,
int Size,
bool IsComplex>
struct direct_selfadjoint_eigenvalues
482 static inline void run(SolverType& eig,
const typename SolverType::MatrixType& A,
int options)
483 { eig.compute(A,options); }
486 template<
typename SolverType>
struct direct_selfadjoint_eigenvalues<SolverType,3,false>
488 typedef typename SolverType::MatrixType MatrixType;
489 typedef typename SolverType::RealVectorType VectorType;
490 typedef typename SolverType::Scalar Scalar;
492 static inline void computeRoots(
const MatrixType& m, VectorType& roots)
498 const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0);
499 const Scalar s_sqrt3 = sqrt(Scalar(3.0));
504 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);
505 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);
506 Scalar c2 = m(0,0) + m(1,1) + m(2,2);
510 Scalar c2_over_3 = c2*s_inv3;
511 Scalar a_over_3 = (c1 - c2*c2_over_3)*s_inv3;
512 if (a_over_3 > Scalar(0))
513 a_over_3 = Scalar(0);
515 Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1));
517 Scalar q = half_b*half_b + a_over_3*a_over_3*a_over_3;
522 Scalar rho = sqrt(-a_over_3);
523 Scalar theta = atan2(sqrt(-q),half_b)*s_inv3;
524 Scalar cos_theta = cos(theta);
525 Scalar sin_theta = sin(theta);
526 roots(0) = c2_over_3 + Scalar(2)*rho*cos_theta;
527 roots(1) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta);
528 roots(2) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta);
531 if (roots(0) >= roots(1))
532 std::swap(roots(0),roots(1));
533 if (roots(1) >= roots(2))
535 std::swap(roots(1),roots(2));
536 if (roots(0) >= roots(1))
537 std::swap(roots(0),roots(1));
541 static inline void run(SolverType& solver,
const MatrixType& mat,
int options)
544 eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows());
545 eigen_assert((options&~(EigVecMask|GenEigMask))==0
546 && (options&EigVecMask)!=EigVecMask
547 &&
"invalid option parameter");
550 MatrixType& eivecs = solver.m_eivec;
551 VectorType& eivals = solver.m_eivalues;
554 Scalar scale = mat.cwiseAbs().maxCoeff();
555 MatrixType scaledMat = mat / scale;
558 computeRoots(scaledMat,eivals);
561 if(computeEigenvectors)
564 safeNorm2 *= safeNorm2;
567 eivecs.setIdentity();
571 scaledMat = scaledMat.template selfadjointView<Lower>();
575 Scalar d0 = eivals(2) - eivals(1);
576 Scalar d1 = eivals(1) - eivals(0);
577 int k = d0 > d1 ? 2 : 0;
578 d0 = d0 > d1 ? d1 : d0;
580 tmp.diagonal().array () -= eivals(k);
583 n = (cross = tmp.row(0).cross(tmp.row(1))).squaredNorm();
586 eivecs.col(k) = cross / sqrt(n);
589 n = (cross = tmp.row(0).cross(tmp.row(2))).squaredNorm();
592 eivecs.col(k) = cross / sqrt(n);
595 n = (cross = tmp.row(1).cross(tmp.row(2))).squaredNorm();
598 eivecs.col(k) = cross / sqrt(n);
607 solver.m_isInitialized =
true;
608 solver.m_eigenvectorsOk = computeEigenvectors;
615 tmp.diagonal().array() -= eivals(1);
618 eivecs.col(1) = eivecs.col(k).unitOrthogonal();
621 n = (cross = eivecs.col(k).cross(tmp.row(0).normalized())).squaredNorm();
623 eivecs.col(1) = cross / sqrt(n);
626 n = (cross = eivecs.col(k).cross(tmp.row(1))).squaredNorm();
628 eivecs.col(1) = cross / sqrt(n);
631 n = (cross = eivecs.col(k).cross(tmp.row(2))).squaredNorm();
633 eivecs.col(1) = cross / sqrt(n);
638 eivecs.col(1) = eivecs.col(k).unitOrthogonal();
644 Scalar d = eivecs.col(1).dot(eivecs.col(k));
645 eivecs.col(1) = (eivecs.col(1) - d * eivecs.col(k)).normalized();
648 eivecs.col(k==2 ? 0 : 2) = eivecs.col(k).cross(eivecs.col(1)).normalized();
655 solver.m_isInitialized =
true;
656 solver.m_eigenvectorsOk = computeEigenvectors;
661 template<
typename SolverType>
struct direct_selfadjoint_eigenvalues<SolverType,2,false>
663 typedef typename SolverType::MatrixType MatrixType;
664 typedef typename SolverType::RealVectorType VectorType;
665 typedef typename SolverType::Scalar Scalar;
667 static inline void computeRoots(
const MatrixType& m, VectorType& roots)
670 const Scalar t0 = Scalar(0.5) * sqrt( abs2(m(0,0)-m(1,1)) + Scalar(4)*m(1,0)*m(1,0));
671 const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1));
676 static inline void run(SolverType& solver,
const MatrixType& mat,
int options)
678 eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
679 eigen_assert((options&~(EigVecMask|GenEigMask))==0
680 && (options&EigVecMask)!=EigVecMask
681 &&
"invalid option parameter");
684 MatrixType& eivecs = solver.m_eivec;
685 VectorType& eivals = solver.m_eivalues;
688 Scalar scale = mat.cwiseAbs().maxCoeff();
689 scale = (std::max)(scale,Scalar(1));
690 MatrixType scaledMat = mat / scale;
693 computeRoots(scaledMat,eivals);
696 if(computeEigenvectors)
698 scaledMat.diagonal().array () -= eivals(1);
699 Scalar a2 = abs2(scaledMat(0,0));
700 Scalar c2 = abs2(scaledMat(1,1));
701 Scalar b2 = abs2(scaledMat(1,0));
704 eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0);
705 eivecs.col(1) /= sqrt(a2+b2);
709 eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0);
710 eivecs.col(1) /= sqrt(c2+b2);
713 eivecs.col(0) << eivecs.col(1).unitOrthogonal();
720 solver.m_isInitialized =
true;
721 solver.m_eigenvectorsOk = computeEigenvectors;
727 template<
typename MatrixType>
731 internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>::run(*
this,matrix,options);
736 template<
int StorageOrder,
typename RealScalar,
typename Scalar,
typename Index>
737 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
739 RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
740 RealScalar e = subdiag[end-1];
746 RealScalar mu = diag[end];
751 RealScalar e2 = abs2(subdiag[end-1]);
752 RealScalar h = hypot(td,e);
753 if(e2==0) mu -= (e / (td + (td>0 ? 1 : -1))) * (e / h);
754 else mu -= e2 / (td + (td>0 ? h : -h));
757 RealScalar x = diag[start] - mu;
758 RealScalar z = subdiag[start];
759 for (Index k = start; k < end; ++k)
761 JacobiRotation<RealScalar> rot;
762 rot.makeGivens(x, z);
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];
788 Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n);
789 q.applyOnTheRight(k,k+1,rot);
798 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H