Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
SelfAdjointEigenSolver.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
5// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
6//
7// This Source Code Form is subject to the terms of the Mozilla
8// Public License v. 2.0. If a copy of the MPL was not distributed
9// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11#ifndef EIGEN_SELFADJOINTEIGENSOLVER_H
12#define EIGEN_SELFADJOINTEIGENSOLVER_H
13
14#include "./Tridiagonalization.h"
15
16namespace Eigen {
17
18template<typename _MatrixType>
19class GeneralizedSelfAdjointEigenSolver;
20
21namespace internal {
22template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues;
23
24template<typename MatrixType, typename DiagType, typename SubDiagType>
25EIGEN_DEVICE_FUNC
26ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec);
27}
28
76template<typename _MatrixType> class SelfAdjointEigenSolver
77{
78 public:
79
80 typedef _MatrixType MatrixType;
81 enum {
82 Size = MatrixType::RowsAtCompileTime,
83 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
84 Options = MatrixType::Options,
85 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
86 };
87
89 typedef typename MatrixType::Scalar Scalar;
91
93
100 typedef typename NumTraits<Scalar>::Real RealScalar;
101
102 friend struct internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>;
103
109 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
110 typedef Tridiagonalization<MatrixType> TridiagonalizationType;
112
123 EIGEN_DEVICE_FUNC
125 : m_eivec(),
126 m_eivalues(),
127 m_subdiag(),
128 m_hcoeffs(),
129 m_info(InvalidInput),
130 m_isInitialized(false),
131 m_eigenvectorsOk(false)
132 { }
133
146 EIGEN_DEVICE_FUNC
148 : m_eivec(size, size),
149 m_eivalues(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)
154 {}
155
171 template<typename InputType>
172 EIGEN_DEVICE_FUNC
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)
180 {
181 compute(matrix.derived(), options);
182 }
183
214 template<typename InputType>
215 EIGEN_DEVICE_FUNC
217
236 EIGEN_DEVICE_FUNC
237 SelfAdjointEigenSolver& computeDirect(const MatrixType& matrix, int options = ComputeEigenvectors);
238
252
276 EIGEN_DEVICE_FUNC
278 {
279 eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
280 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
281 return m_eivec;
282 }
283
299 EIGEN_DEVICE_FUNC
301 {
302 eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
303 return m_eivalues;
304 }
305
323 EIGEN_DEVICE_FUNC
324 MatrixType operatorSqrt() const
325 {
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();
329 }
330
348 EIGEN_DEVICE_FUNC
349 MatrixType operatorInverseSqrt() const
350 {
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();
354 }
355
360 EIGEN_DEVICE_FUNC
362 {
363 eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
364 return m_info;
365 }
366
372 static const int m_maxIterations = 30;
373
374 protected:
375 static EIGEN_DEVICE_FUNC
376 void check_template_parameters()
377 {
378 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
379 }
380
381 EigenvectorsType m_eivec;
382 RealVectorType m_eivalues;
383 typename TridiagonalizationType::SubDiagonalType m_subdiag;
384 typename TridiagonalizationType::CoeffVectorType m_hcoeffs;
385 ComputationInfo m_info;
386 bool m_isInitialized;
387 bool m_eigenvectorsOk;
388};
389
390namespace internal {
411template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
412EIGEN_DEVICE_FUNC
413static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
414}
415
416template<typename MatrixType>
417template<typename InputType>
418EIGEN_DEVICE_FUNC
419SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
420::compute(const EigenBase<InputType>& a_matrix, int options)
421{
422 check_template_parameters();
423
424 const InputType &matrix(a_matrix.derived());
425
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");
431 bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
432 Index n = matrix.cols();
433 m_eivalues.resize(n,1);
434
435 if(n==1)
436 {
437 m_eivec = matrix;
438 m_eivalues.coeffRef(0,0) = numext::real(m_eivec.coeff(0,0));
439 if(computeEigenvectors)
440 m_eivec.setOnes(n,n);
441 m_info = Success;
442 m_isInitialized = true;
443 m_eigenvectorsOk = computeEigenvectors;
444 return *this;
445 }
446
447 // declare some aliases
448 RealVectorType& diag = m_eivalues;
449 EigenvectorsType& mat = m_eivec;
450
451 // map the matrix coefficients to [-1:1] to avoid over- and underflow.
452 mat = matrix.template triangularView<Lower>();
453 RealScalar scale = mat.cwiseAbs().maxCoeff();
454 if(scale==RealScalar(0)) scale = RealScalar(1);
455 mat.template triangularView<Lower>() /= scale;
456 m_subdiag.resize(n-1);
457 m_hcoeffs.resize(n-1);
458 internal::tridiagonalization_inplace(mat, diag, m_subdiag, m_hcoeffs, computeEigenvectors);
459
460 m_info = internal::computeFromTridiagonal_impl(diag, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
461
462 // scale back the eigen values
463 m_eivalues *= scale;
464
465 m_isInitialized = true;
466 m_eigenvectorsOk = computeEigenvectors;
467 return *this;
468}
469
470template<typename MatrixType>
471SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
472::computeFromTridiagonal(const RealVectorType& diag, const SubDiagonalType& subdiag , int options)
473{
474 //TODO : Add an option to scale the values beforehand
475 bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
476
477 m_eivalues = diag;
478 m_subdiag = subdiag;
479 if (computeEigenvectors)
480 {
481 m_eivec.setIdentity(diag.size(), diag.size());
482 }
483 m_info = internal::computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
484
485 m_isInitialized = true;
486 m_eigenvectorsOk = computeEigenvectors;
487 return *this;
488}
489
490namespace internal {
502template<typename MatrixType, typename DiagType, typename SubDiagType>
503EIGEN_DEVICE_FUNC
504ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec)
505{
506 ComputationInfo info;
507 typedef typename MatrixType::Scalar Scalar;
508
509 Index n = diag.size();
510 Index end = n-1;
511 Index start = 0;
512 Index iter = 0; // total number of iterations
513
514 typedef typename DiagType::RealScalar RealScalar;
515 const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
516 const RealScalar precision_inv = RealScalar(1)/NumTraits<RealScalar>::epsilon();
517 while (end>0)
518 {
519 for (Index i = start; i<end; ++i) {
520 if (numext::abs(subdiag[i]) < considerAsZero) {
521 subdiag[i] = RealScalar(0);
522 } else {
523 // abs(subdiag[i]) <= epsilon * sqrt(abs(diag[i]) + abs(diag[i+1]))
524 // Scaled to prevent underflows.
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);
528 }
529 }
530 }
531
532 // find the largest unreduced block at the end of the matrix.
533 while (end>0 && subdiag[end-1]==RealScalar(0))
534 {
535 end--;
536 }
537 if (end<=0)
538 break;
539
540 // if we spent too many iterations, we give up
541 iter++;
542 if(iter > maxIterations * n) break;
543
544 start = end - 1;
545 while (start>0 && subdiag[start-1]!=0)
546 start--;
547
548 internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), subdiag.data(), start, end, computeEigenvectors ? eivec.data() : (Scalar*)0, n);
549 }
550 if (iter <= maxIterations * n)
551 info = Success;
552 else
553 info = NoConvergence;
554
555 // Sort eigenvalues and corresponding vectors.
556 // TODO make the sort optional ?
557 // TODO use a better sort algorithm !!
558 if (info == Success)
559 {
560 for (Index i = 0; i < n-1; ++i)
561 {
562 Index k;
563 diag.segment(i,n-i).minCoeff(&k);
564 if (k > 0)
565 {
566 numext::swap(diag[i], diag[k+i]);
567 if(computeEigenvectors)
568 eivec.col(i).swap(eivec.col(k+i));
569 }
570 }
571 }
572 return info;
573}
574
575template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues
576{
577 EIGEN_DEVICE_FUNC
578 static inline void run(SolverType& eig, const typename SolverType::MatrixType& A, int options)
579 { eig.compute(A,options); }
580};
581
582template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3,false>
583{
584 typedef typename SolverType::MatrixType MatrixType;
585 typedef typename SolverType::RealVectorType VectorType;
586 typedef typename SolverType::Scalar Scalar;
587 typedef typename SolverType::EigenvectorsType EigenvectorsType;
588
589
594 EIGEN_DEVICE_FUNC
595 static inline void computeRoots(const MatrixType& m, VectorType& roots)
596 {
597 EIGEN_USING_STD(sqrt)
598 EIGEN_USING_STD(atan2)
599 EIGEN_USING_STD(cos)
600 EIGEN_USING_STD(sin)
601 const Scalar s_inv3 = Scalar(1)/Scalar(3);
602 const Scalar s_sqrt3 = sqrt(Scalar(3));
603
604 // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0. The
605 // eigenvalues are the roots to this equation, all guaranteed to be
606 // real-valued, because the matrix is symmetric.
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);
610
611 // Construct the parameters used in classifying the roots of the equation
612 // and in solving the equation for the roots in closed form.
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));
616
617 Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1));
618
619 Scalar q = a_over_3*a_over_3*a_over_3 - half_b*half_b;
620 q = numext::maxi(q, Scalar(0));
621
622 // Compute the eigenvalues by solving for the roots of the polynomial.
623 Scalar rho = sqrt(a_over_3);
624 Scalar theta = atan2(sqrt(q),half_b)*s_inv3; // since sqrt(q) > 0, atan2 is in [0, pi] and theta is in [0, pi/3]
625 Scalar cos_theta = cos(theta);
626 Scalar sin_theta = sin(theta);
627 // roots are already sorted, since cos is monotonically decreasing on [0, pi]
628 roots(0) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta); // == 2*rho*cos(theta+2pi/3)
629 roots(1) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta); // == 2*rho*cos(theta+ pi/3)
630 roots(2) = c2_over_3 + Scalar(2)*rho*cos_theta;
631 }
632
633 EIGEN_DEVICE_FUNC
634 static inline bool extract_kernel(MatrixType& mat, Ref<VectorType> res, Ref<VectorType> representative)
635 {
636 EIGEN_USING_STD(abs);
637 EIGEN_USING_STD(sqrt);
638 Index i0;
639 // Find non-zero column i0 (by construction, there must exist a non zero coefficient on the diagonal):
640 mat.diagonal().cwiseAbs().maxCoeff(&i0);
641 // mat.col(i0) is a good candidate for an orthogonal vector to the current eigenvector,
642 // so let's save it:
643 representative = mat.col(i0);
644 Scalar n0, n1;
645 VectorType c0, c1;
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);
650
651 return true;
652 }
653
654 EIGEN_DEVICE_FUNC
655 static inline void run(SolverType& solver, const MatrixType& mat, int options)
656 {
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");
661 bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
662
663 EigenvectorsType& eivecs = solver.m_eivec;
664 VectorType& eivals = solver.m_eivalues;
665
666 // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow.
667 Scalar shift = mat.trace() / Scalar(3);
668 // TODO Avoid this copy. Currently it is necessary to suppress bogus values when determining maxCoeff and for computing the eigenvectors later
669 MatrixType scaledMat = mat.template selfadjointView<Lower>();
670 scaledMat.diagonal().array() -= shift;
671 Scalar scale = scaledMat.cwiseAbs().maxCoeff();
672 if(scale > 0) scaledMat /= scale; // TODO for scale==0 we could save the remaining operations
673
674 // compute the eigenvalues
675 computeRoots(scaledMat,eivals);
676
677 // compute the eigenvectors
678 if(computeEigenvectors)
679 {
680 if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::epsilon())
681 {
682 // All three eigenvalues are numerically the same
683 eivecs.setIdentity();
684 }
685 else
686 {
687 MatrixType tmp;
688 tmp = scaledMat;
689
690 // Compute the eigenvector of the most distinct eigenvalue
691 Scalar d0 = eivals(2) - eivals(1);
692 Scalar d1 = eivals(1) - eivals(0);
693 Index k(0), l(2);
694 if(d0 > d1)
695 {
696 numext::swap(k,l);
697 d0 = d1;
698 }
699
700 // Compute the eigenvector of index k
701 {
702 tmp.diagonal().array () -= eivals(k);
703 // By construction, 'tmp' is of rank 2, and its kernel corresponds to the respective eigenvector.
704 extract_kernel(tmp, eivecs.col(k), eivecs.col(l));
705 }
706
707 // Compute eigenvector of index l
709 {
710 // If d0 is too small, then the two other eigenvalues are numerically the same,
711 // and thus we only have to ortho-normalize the near orthogonal vector we saved above.
712 eivecs.col(l) -= eivecs.col(k).dot(eivecs.col(l))*eivecs.col(l);
713 eivecs.col(l).normalize();
714 }
715 else
716 {
717 tmp = scaledMat;
718 tmp.diagonal().array () -= eivals(l);
719
720 VectorType dummy;
721 extract_kernel(tmp, eivecs.col(l), dummy);
722 }
723
724 // Compute last eigenvector from the other two
725 eivecs.col(1) = eivecs.col(2).cross(eivecs.col(0)).normalized();
726 }
727 }
728
729 // Rescale back to the original size.
730 eivals *= scale;
731 eivals.array() += shift;
732
733 solver.m_info = Success;
734 solver.m_isInitialized = true;
735 solver.m_eigenvectorsOk = computeEigenvectors;
736 }
737};
738
739// 2x2 direct eigenvalues decomposition, code from Hauke Heibel
740template<typename SolverType>
741struct direct_selfadjoint_eigenvalues<SolverType,2,false>
742{
743 typedef typename SolverType::MatrixType MatrixType;
744 typedef typename SolverType::RealVectorType VectorType;
745 typedef typename SolverType::Scalar Scalar;
746 typedef typename SolverType::EigenvectorsType EigenvectorsType;
747
748 EIGEN_DEVICE_FUNC
749 static inline void computeRoots(const MatrixType& m, VectorType& roots)
750 {
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));
754 roots(0) = t1 - t0;
755 roots(1) = t1 + t0;
756 }
757
758 EIGEN_DEVICE_FUNC
759 static inline void run(SolverType& solver, const MatrixType& mat, int options)
760 {
761 EIGEN_USING_STD(sqrt);
762 EIGEN_USING_STD(abs);
763
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");
768 bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
769
770 EigenvectorsType& eivecs = solver.m_eivec;
771 VectorType& eivals = solver.m_eivalues;
772
773 // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow.
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))
780 scaledMat /= scale;
781
782 // Compute the eigenvalues
783 computeRoots(scaledMat,eivals);
784
785 // compute the eigen vectors
786 if(computeEigenvectors)
787 {
788 if((eivals(1)-eivals(0))<=abs(eivals(1))*Eigen::NumTraits<Scalar>::epsilon())
789 {
790 eivecs.setIdentity();
791 }
792 else
793 {
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));
798 if(a2>c2)
799 {
800 eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0);
801 eivecs.col(1) /= sqrt(a2+b2);
802 }
803 else
804 {
805 eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0);
806 eivecs.col(1) /= sqrt(c2+b2);
807 }
808
809 eivecs.col(0) << eivecs.col(1).unitOrthogonal();
810 }
811 }
812
813 // Rescale back to the original size.
814 eivals *= scale;
815 eivals.array() += shift;
816
817 solver.m_info = Success;
818 solver.m_isInitialized = true;
819 solver.m_eigenvectorsOk = computeEigenvectors;
820 }
821};
822
823}
824
825template<typename MatrixType>
826EIGEN_DEVICE_FUNC
827SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
828::computeDirect(const MatrixType& matrix, int options)
829{
830 internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>::run(*this,matrix,options);
831 return *this;
832}
833
834namespace internal {
835
836// Francis implicit QR step.
837template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
838EIGEN_DEVICE_FUNC
839static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
840{
841 // Wilkinson Shift.
842 RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
843 RealScalar e = subdiag[end-1];
844 // Note that thanks to scaling, e^2 or td^2 cannot overflow, however they can still
845 // underflow thus leading to inf/NaN values when using the following commented code:
846 // RealScalar e2 = numext::abs2(subdiag[end-1]);
847 // RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
848 // This explain the following, somewhat more complicated, version:
849 RealScalar mu = diag[end];
850 if(td==RealScalar(0)) {
851 mu -= numext::abs(e);
852 } else if (e != RealScalar(0)) {
853 const RealScalar e2 = numext::abs2(e);
854 const RealScalar h = numext::hypot(td,e);
855 if(e2 == RealScalar(0)) {
856 mu -= e / ((td + (td>RealScalar(0) ? h : -h)) / e);
857 } else {
858 mu -= e2 / (td + (td>RealScalar(0) ? h : -h));
859 }
860 }
861
862 RealScalar x = diag[start] - mu;
863 RealScalar z = subdiag[start];
864 // If z ever becomes zero, the Givens rotation will be the identity and
865 // z will stay zero for all future iterations.
866 for (Index k = start; k < end && z != RealScalar(0); ++k)
867 {
868 JacobiRotation<RealScalar> rot;
869 rot.makeGivens(x, z);
870
871 // do T = G' T G
872 RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
873 RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
874
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;
878
879 if (k > start)
880 subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
881
882 // "Chasing the bulge" to return to triangular form.
883 x = subdiag[k];
884 if (k < end - 1)
885 {
886 z = -rot.s() * subdiag[k+1];
887 subdiag[k + 1] = rot.c() * subdiag[k+1];
888 }
889
890 // apply the givens rotation to the unit matrix Q = Q * G
891 if (matrixQ)
892 {
893 // FIXME if StorageOrder == RowMajor this operation is not very efficient
894 Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n);
895 q.applyOnTheRight(k,k+1,rot);
896 }
897 }
898}
899
900} // end namespace internal
901
902} // end namespace Eigen
903
904#endif // EIGEN_SELFADJOINTEIGENSOLVER_H
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Computes eigenvalues and eigenvectors of selfadjoint matrices.
Definition: SelfAdjointEigenSolver.h:77
MatrixType::Scalar Scalar
Scalar type for matrices of type _MatrixType.
Definition: SelfAdjointEigenSolver.h:89
SelfAdjointEigenSolver & computeFromTridiagonal(const RealVectorType &diag, const SubDiagonalType &subdiag, int options=ComputeEigenvectors)
Computes the eigen decomposition from a tridiagonal symmetric matrix.
Definition: SelfAdjointEigenSolver.h:472
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.
const RealVectorType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: SelfAdjointEigenSolver.h:300
NumTraits< Scalar >::Real RealScalar
Real scalar type for _MatrixType.
Definition: SelfAdjointEigenSolver.h:100
const EigenvectorsType & eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: SelfAdjointEigenSolver.h:277
SelfAdjointEigenSolver(const EigenBase< InputType > &matrix, int options=ComputeEigenvectors)
Constructor; computes eigendecomposition of given matrix.
Definition: SelfAdjointEigenSolver.h:173
internal::plain_col_type< MatrixType, RealScalar >::type RealVectorType
Type for vector of eigenvalues as returned by eigenvalues().
Definition: SelfAdjointEigenSolver.h:109
SelfAdjointEigenSolver(Index size)
Constructor, pre-allocates memory for dynamic-size matrices.
Definition: SelfAdjointEigenSolver.h:147
Eigen::Index Index
Definition: SelfAdjointEigenSolver.h:90
static const int m_maxIterations
Maximum number of iterations.
Definition: SelfAdjointEigenSolver.h:372
MatrixType operatorInverseSqrt() const
Computes the inverse square root of the matrix.
Definition: SelfAdjointEigenSolver.h:349
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
Tridiagonal decomposition of a selfadjoint matrix.
Definition: Tridiagonalization.h:65
ComputationInfo
Definition: Constants.h:440
@ Success
Definition: Constants.h:442
@ NoConvergence
Definition: Constants.h:446
@ ComputeEigenvectors
Definition: Constants.h:405
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Definition: EigenBase.h:30
Derived & derived()
Definition: EigenBase.h:46
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:233