Eigen  3.4.0 (git rev e3e74001f7c4bf95f0dde572e8a08c5b2918a3ab)
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 
16 namespace Eigen {
17 
18 template<typename _MatrixType>
19 class GeneralizedSelfAdjointEigenSolver;
20 
21 namespace internal {
22 template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues;
23 
24 template<typename MatrixType, typename DiagType, typename SubDiagType>
25 EIGEN_DEVICE_FUNC
26 ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec);
27 }
28 
76 template<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;
90  typedef Eigen::Index Index;
91 
93 
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
173  explicit SelfAdjointEigenSolver(const EigenBase<InputType>& matrix, int options = ComputeEigenvectors)
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 
390 namespace internal {
411 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
412 EIGEN_DEVICE_FUNC
413 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
414 }
415 
416 template<typename MatrixType>
417 template<typename InputType>
418 EIGEN_DEVICE_FUNC
419 SelfAdjointEigenSolver<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 
470 template<typename MatrixType>
471 SelfAdjointEigenSolver<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 
490 namespace internal {
502 template<typename MatrixType, typename DiagType, typename SubDiagType>
503 EIGEN_DEVICE_FUNC
504 ComputationInfo 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 
575 template<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 
582 template<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
708  if(d0<=2*Eigen::NumTraits<Scalar>::epsilon()*d1)
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
740 template<typename SolverType>
741 struct 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 
825 template<typename MatrixType>
826 EIGEN_DEVICE_FUNC
827 SelfAdjointEigenSolver<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 
834 namespace internal {
835 
836 // Francis implicit QR step.
837 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
838 EIGEN_DEVICE_FUNC
839 static 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
Derived & setIdentity()
Definition: CwiseNullaryOp.h:873
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
const EigenvectorsType & eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: SelfAdjointEigenSolver.h:277
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
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
SelfAdjointEigenSolver(const EigenBase< InputType > &matrix, int options=ComputeEigenvectors)
Constructor; computes eigendecomposition of given matrix.
Definition: SelfAdjointEigenSolver.h:173
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
internal::plain_col_type< MatrixType, RealScalar >::type RealVectorType
Type for vector of eigenvalues as returned by eigenvalues().
Definition: SelfAdjointEigenSolver.h:109
SelfAdjointEigenSolver & compute(const EigenBase< InputType > &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix.
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
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_cos_op< typename Derived::Scalar >, const Derived > cos(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sin_op< typename Derived::Scalar >, const Derived > sin(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
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