26 #ifndef EIGEN_SELFADJOINTEIGENSOLVER_H
27 #define EIGEN_SELFADJOINTEIGENSOLVER_H
33 template<
typename _MatrixType>
34 class GeneralizedSelfAdjointEigenSolver;
37 template<
typename SolverType,
int Size,
bool IsComplex>
struct direct_selfadjoint_eigenvalues;
89 Size = MatrixType::RowsAtCompileTime,
96 typedef typename MatrixType::Scalar
Scalar;
97 typedef typename MatrixType::Index
Index;
114 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type
RealVectorType;
169 :
m_eivec(matrix.rows(), matrix.cols()),
171 m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
338 #ifdef EIGEN2_SUPPORT
340 :
m_eivec(matrix.rows(), matrix.cols()),
342 m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
345 compute(matrix, computeEigenvectors);
349 :
m_eivec(matA.cols(), matA.cols()),
351 m_subdiag(matA.cols() > 1 ? matA.cols() - 1 : 1),
366 #endif // EIGEN2_SUPPORT
394 template<
int StorageOrder,
typename RealScalar,
typename Scalar,
typename Index>
395 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
398 template<
typename MatrixType>
405 &&
"invalid option parameter");
407 Index n = matrix.cols();
408 m_eivalues.resize(n,1);
413 if(computeEigenvectors)
414 m_eivec.setOnes(n,n);
416 m_isInitialized =
true;
417 m_eigenvectorsOk = computeEigenvectors;
426 RealScalar scale = matrix.cwiseAbs().maxCoeff();
428 mat = matrix / scale;
429 m_subdiag.resize(n-1);
438 for (
Index i = start; i<end; ++i)
443 while (end>0 && m_subdiag[end-1]==0)
452 if(iter > m_maxIterations * n)
break;
455 while (start>0 && m_subdiag[start-1]!=0)
458 internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (
Scalar*)0, n);
461 if (iter <= m_maxIterations * n)
471 for (
Index i = 0; i < n-1; ++i)
474 m_eivalues.segment(i,n-i).minCoeff(&k);
477 std::swap(m_eivalues[i], m_eivalues[k+i]);
478 if(computeEigenvectors)
479 m_eivec.col(i).swap(m_eivec.col(k+i));
487 m_isInitialized =
true;
488 m_eigenvectorsOk = computeEigenvectors;
495 template<
typename SolverType,
int Size,
bool IsComplex>
struct direct_selfadjoint_eigenvalues
497 static inline void run(SolverType& eig,
const typename SolverType::MatrixType& A,
int options)
498 { eig.compute(A,options); }
501 template<
typename SolverType>
struct direct_selfadjoint_eigenvalues<SolverType,3,false>
503 typedef typename SolverType::MatrixType MatrixType;
504 typedef typename SolverType::RealVectorType VectorType;
505 typedef typename SolverType::Scalar Scalar;
507 static inline void computeRoots(
const MatrixType& m, VectorType& roots)
513 const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0);
514 const Scalar s_sqrt3 =
sqrt(Scalar(3.0));
519 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);
520 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);
521 Scalar c2 = m(0,0) + m(1,1) + m(2,2);
525 Scalar c2_over_3 = c2*s_inv3;
526 Scalar a_over_3 = (c1 - c2*c2_over_3)*s_inv3;
527 if (a_over_3 > Scalar(0))
528 a_over_3 = Scalar(0);
530 Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1));
532 Scalar
q = half_b*half_b + a_over_3*a_over_3*a_over_3;
537 Scalar rho =
sqrt(-a_over_3);
538 Scalar theta = atan2(
sqrt(-q),half_b)*s_inv3;
539 Scalar cos_theta =
cos(theta);
540 Scalar sin_theta =
sin(theta);
541 roots(0) = c2_over_3 + Scalar(2)*rho*cos_theta;
542 roots(1) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta);
543 roots(2) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta);
546 if (roots(0) >= roots(1))
547 std::swap(roots(0),roots(1));
548 if (roots(1) >= roots(2))
550 std::swap(roots(1),roots(2));
551 if (roots(0) >= roots(1))
552 std::swap(roots(0),roots(1));
556 static inline void run(SolverType& solver,
const MatrixType& mat,
int options)
559 eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows());
562 &&
"invalid option parameter");
565 MatrixType& eivecs = solver.m_eivec;
566 VectorType& eivals = solver.m_eivalues;
569 Scalar scale = mat.cwiseAbs().maxCoeff();
570 MatrixType scaledMat = mat / scale;
573 computeRoots(scaledMat,eivals);
576 if(computeEigenvectors)
579 safeNorm2 *= safeNorm2;
582 eivecs.setIdentity();
586 scaledMat = scaledMat.template selfadjointView<Lower>();
590 Scalar d0 = eivals(2) - eivals(1);
591 Scalar d1 = eivals(1) - eivals(0);
592 int k = d0 > d1 ? 2 : 0;
593 d0 = d0 > d1 ? d1 : d0;
595 tmp.diagonal().array () -= eivals(k);
598 n = (cross = tmp.row(0).cross(tmp.row(1))).squaredNorm();
601 eivecs.col(k) = cross /
sqrt(n);
604 n = (cross = tmp.row(0).cross(tmp.row(2))).squaredNorm();
607 eivecs.col(k) = cross /
sqrt(n);
610 n = (cross = tmp.row(1).cross(tmp.row(2))).squaredNorm();
613 eivecs.col(k) = cross /
sqrt(n);
622 solver.m_isInitialized =
true;
623 solver.m_eigenvectorsOk = computeEigenvectors;
630 tmp.diagonal().array() -= eivals(1);
633 eivecs.col(1) = eivecs.col(k).unitOrthogonal();
636 n = (cross = eivecs.col(k).cross(tmp.row(0).normalized())).squaredNorm();
638 eivecs.col(1) = cross /
sqrt(n);
641 n = (cross = eivecs.col(k).cross(tmp.row(1))).squaredNorm();
643 eivecs.col(1) = cross /
sqrt(n);
646 n = (cross = eivecs.col(k).cross(tmp.row(2))).squaredNorm();
648 eivecs.col(1) = cross /
sqrt(n);
653 eivecs.col(1) = eivecs.col(k).unitOrthogonal();
659 Scalar d = eivecs.col(1).dot(eivecs.col(k));
660 eivecs.col(1) = (eivecs.col(1) - d * eivecs.col(k)).normalized();
663 eivecs.col(k==2 ? 0 : 2) = eivecs.col(k).cross(eivecs.col(1)).normalized();
670 solver.m_isInitialized =
true;
671 solver.m_eigenvectorsOk = computeEigenvectors;
676 template<
typename SolverType>
struct direct_selfadjoint_eigenvalues<SolverType,2,false>
678 typedef typename SolverType::MatrixType MatrixType;
679 typedef typename SolverType::RealVectorType VectorType;
680 typedef typename SolverType::Scalar Scalar;
682 static inline void computeRoots(
const MatrixType& m, VectorType& roots)
685 const Scalar t0 = Scalar(0.5) *
sqrt(
abs2(m(0,0)-m(1,1)) + Scalar(4)*m(1,0)*m(1,0));
686 const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1));
691 static inline void run(SolverType& solver,
const MatrixType& mat,
int options)
693 eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
696 &&
"invalid option parameter");
699 MatrixType& eivecs = solver.m_eivec;
700 VectorType& eivals = solver.m_eivalues;
703 Scalar scale = mat.cwiseAbs().maxCoeff();
704 scale = (std::max)(scale,Scalar(1));
705 MatrixType scaledMat = mat / scale;
708 computeRoots(scaledMat,eivals);
711 if(computeEigenvectors)
713 scaledMat.diagonal().array () -= eivals(1);
714 Scalar a2 =
abs2(scaledMat(0,0));
715 Scalar c2 =
abs2(scaledMat(1,1));
716 Scalar b2 =
abs2(scaledMat(1,0));
719 eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0);
720 eivecs.col(1) /=
sqrt(a2+b2);
724 eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0);
725 eivecs.col(1) /=
sqrt(c2+b2);
728 eivecs.col(0) << eivecs.col(1).unitOrthogonal();
735 solver.m_isInitialized =
true;
736 solver.m_eigenvectorsOk = computeEigenvectors;
742 template<
typename MatrixType>
751 template<
int StorageOrder,
typename RealScalar,
typename Scalar,
typename Index>
752 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
760 RealScalar
td = (diag[end-1] - diag[end])*RealScalar(0.5);
761 RealScalar e2 =
abs2(subdiag[end-1]);
762 RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) *
sqrt(td*td + e2));
763 RealScalar x = diag[start] - mu;
764 RealScalar z = subdiag[start];
765 for (Index k = start; k < end; ++k)
771 RealScalar sdk = rot.
s() * diag[k] + rot.
c() * subdiag[k];
772 RealScalar dkp1 = rot.
s() * subdiag[k] + rot.
c() * diag[k+1];
774 diag[k] = rot.
c() * (rot.
c() * diag[k] - rot.
s() * subdiag[k]) - rot.
s() * (rot.
c() * subdiag[k] - rot.
s() * diag[k+1]);
775 diag[k+1] = rot.
s() * sdk + rot.
c() * dkp1;
776 subdiag[k] = rot.
c() * sdk - rot.
s() * dkp1;
780 subdiag[k - 1] = rot.
c() * subdiag[k-1] - rot.
s() * z;
786 z = -rot.
s() * subdiag[k+1];
787 subdiag[k + 1] = rot.
c() * subdiag[k+1];
794 Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> >
q(matrixQ,n,n);
795 q.applyOnTheRight(k,k+1,rot);
804 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H