Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.3.8
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 template<typename MatrixType, typename DiagType, typename SubDiagType>
24 ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec);
25 }
26 
70 template<typename _MatrixType> class SelfAdjointEigenSolver
71 {
72  public:
73 
74  typedef _MatrixType MatrixType;
75  enum {
76  Size = MatrixType::RowsAtCompileTime,
77  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
78  Options = MatrixType::Options,
79  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
80  };
81 
83  typedef typename MatrixType::Scalar Scalar;
84  typedef Eigen::Index Index;
85 
87 
95 
96  friend struct internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>;
97 
103  typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
104  typedef Tridiagonalization<MatrixType> TridiagonalizationType;
106 
117  EIGEN_DEVICE_FUNC
119  : m_eivec(),
120  m_eivalues(),
121  m_subdiag(),
122  m_isInitialized(false)
123  { }
124 
137  EIGEN_DEVICE_FUNC
139  : m_eivec(size, size),
140  m_eivalues(size),
141  m_subdiag(size > 1 ? size - 1 : 1),
142  m_isInitialized(false)
143  {}
144 
160  template<typename InputType>
161  EIGEN_DEVICE_FUNC
162  explicit SelfAdjointEigenSolver(const EigenBase<InputType>& matrix, int options = ComputeEigenvectors)
163  : m_eivec(matrix.rows(), matrix.cols()),
164  m_eivalues(matrix.cols()),
165  m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
166  m_isInitialized(false)
167  {
168  compute(matrix.derived(), options);
169  }
170 
201  template<typename InputType>
202  EIGEN_DEVICE_FUNC
204 
223  EIGEN_DEVICE_FUNC
224  SelfAdjointEigenSolver& computeDirect(const MatrixType& matrix, int options = ComputeEigenvectors);
225 
239 
258  EIGEN_DEVICE_FUNC
260  {
261  eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
262  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
263  return m_eivec;
264  }
265 
281  EIGEN_DEVICE_FUNC
283  {
284  eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
285  return m_eivalues;
286  }
287 
305  EIGEN_DEVICE_FUNC
306  MatrixType operatorSqrt() const
307  {
308  eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
309  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
310  return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
311  }
312 
330  EIGEN_DEVICE_FUNC
331  MatrixType operatorInverseSqrt() const
332  {
333  eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
334  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
335  return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
336  }
337 
342  EIGEN_DEVICE_FUNC
344  {
345  eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
346  return m_info;
347  }
348 
354  static const int m_maxIterations = 30;
355 
356  protected:
357  static void check_template_parameters()
358  {
359  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
360  }
361 
362  EigenvectorsType m_eivec;
363  RealVectorType m_eivalues;
364  typename TridiagonalizationType::SubDiagonalType m_subdiag;
365  ComputationInfo m_info;
366  bool m_isInitialized;
367  bool m_eigenvectorsOk;
368 };
369 
370 namespace internal {
391 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
392 EIGEN_DEVICE_FUNC
393 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
394 }
395 
396 template<typename MatrixType>
397 template<typename InputType>
398 EIGEN_DEVICE_FUNC
399 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
400 ::compute(const EigenBase<InputType>& a_matrix, int options)
401 {
402  check_template_parameters();
403 
404  const InputType &matrix(a_matrix.derived());
405 
406  using std::abs;
407  eigen_assert(matrix.cols() == matrix.rows());
408  eigen_assert((options&~(EigVecMask|GenEigMask))==0
409  && (options&EigVecMask)!=EigVecMask
410  && "invalid option parameter");
411  bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
412  Index n = matrix.cols();
413  m_eivalues.resize(n,1);
414 
415  if(n==1)
416  {
417  m_eivec = matrix;
418  m_eivalues.coeffRef(0,0) = numext::real(m_eivec.coeff(0,0));
419  if(computeEigenvectors)
420  m_eivec.setOnes(n,n);
421  m_info = Success;
422  m_isInitialized = true;
423  m_eigenvectorsOk = computeEigenvectors;
424  return *this;
425  }
426 
427  // declare some aliases
428  RealVectorType& diag = m_eivalues;
429  EigenvectorsType& mat = m_eivec;
430 
431  // map the matrix coefficients to [-1:1] to avoid over- and underflow.
432  mat = matrix.template triangularView<Lower>();
433  RealScalar scale = mat.cwiseAbs().maxCoeff();
434  if(scale==RealScalar(0)) scale = RealScalar(1);
435  mat.template triangularView<Lower>() /= scale;
436  m_subdiag.resize(n-1);
437  internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors);
438 
439  m_info = internal::computeFromTridiagonal_impl(diag, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
440 
441  // scale back the eigen values
442  m_eivalues *= scale;
443 
444  m_isInitialized = true;
445  m_eigenvectorsOk = computeEigenvectors;
446  return *this;
447 }
448 
449 template<typename MatrixType>
450 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
451 ::computeFromTridiagonal(const RealVectorType& diag, const SubDiagonalType& subdiag , int options)
452 {
453  //TODO : Add an option to scale the values beforehand
454  bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
455 
456  m_eivalues = diag;
457  m_subdiag = subdiag;
458  if (computeEigenvectors)
459  {
460  m_eivec.setIdentity(diag.size(), diag.size());
461  }
462  m_info = internal::computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
463 
464  m_isInitialized = true;
465  m_eigenvectorsOk = computeEigenvectors;
466  return *this;
467 }
468 
469 namespace internal {
481 template<typename MatrixType, typename DiagType, typename SubDiagType>
482 ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec)
483 {
484  using std::abs;
485 
486  ComputationInfo info;
487  typedef typename MatrixType::Scalar Scalar;
488 
489  Index n = diag.size();
490  Index end = n-1;
491  Index start = 0;
492  Index iter = 0; // total number of iterations
493 
494  typedef typename DiagType::RealScalar RealScalar;
495  const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
496  const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon();
497 
498  while (end>0)
499  {
500  for (Index i = start; i<end; ++i)
501  if (internal::isMuchSmallerThan(abs(subdiag[i]),(abs(diag[i])+abs(diag[i+1])),precision) || abs(subdiag[i]) <= considerAsZero)
502  subdiag[i] = 0;
503 
504  // find the largest unreduced block
505  while (end>0 && subdiag[end-1]==RealScalar(0))
506  {
507  end--;
508  }
509  if (end<=0)
510  break;
511 
512  // if we spent too many iterations, we give up
513  iter++;
514  if(iter > maxIterations * n) break;
515 
516  start = end - 1;
517  while (start>0 && subdiag[start-1]!=0)
518  start--;
519 
520  internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), subdiag.data(), start, end, computeEigenvectors ? eivec.data() : (Scalar*)0, n);
521  }
522  if (iter <= maxIterations * n)
523  info = Success;
524  else
525  info = NoConvergence;
526 
527  // Sort eigenvalues and corresponding vectors.
528  // TODO make the sort optional ?
529  // TODO use a better sort algorithm !!
530  if (info == Success)
531  {
532  for (Index i = 0; i < n-1; ++i)
533  {
534  Index k;
535  diag.segment(i,n-i).minCoeff(&k);
536  if (k > 0)
537  {
538  std::swap(diag[i], diag[k+i]);
539  if(computeEigenvectors)
540  eivec.col(i).swap(eivec.col(k+i));
541  }
542  }
543  }
544  return info;
545 }
546 
547 template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues
548 {
549  EIGEN_DEVICE_FUNC
550  static inline void run(SolverType& eig, const typename SolverType::MatrixType& A, int options)
551  { eig.compute(A,options); }
552 };
553 
554 template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3,false>
555 {
556  typedef typename SolverType::MatrixType MatrixType;
557  typedef typename SolverType::RealVectorType VectorType;
558  typedef typename SolverType::Scalar Scalar;
559  typedef typename SolverType::EigenvectorsType EigenvectorsType;
560 
561 
566  EIGEN_DEVICE_FUNC
567  static inline void computeRoots(const MatrixType& m, VectorType& roots)
568  {
569  EIGEN_USING_STD_MATH(sqrt)
570  EIGEN_USING_STD_MATH(atan2)
571  EIGEN_USING_STD_MATH(cos)
572  EIGEN_USING_STD_MATH(sin)
573  const Scalar s_inv3 = Scalar(1)/Scalar(3);
574  const Scalar s_sqrt3 = sqrt(Scalar(3));
575 
576  // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0. The
577  // eigenvalues are the roots to this equation, all guaranteed to be
578  // real-valued, because the matrix is symmetric.
579  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);
580  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);
581  Scalar c2 = m(0,0) + m(1,1) + m(2,2);
582 
583  // Construct the parameters used in classifying the roots of the equation
584  // and in solving the equation for the roots in closed form.
585  Scalar c2_over_3 = c2*s_inv3;
586  Scalar a_over_3 = (c2*c2_over_3 - c1)*s_inv3;
587  a_over_3 = numext::maxi(a_over_3, Scalar(0));
588 
589  Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1));
590 
591  Scalar q = a_over_3*a_over_3*a_over_3 - half_b*half_b;
592  q = numext::maxi(q, Scalar(0));
593 
594  // Compute the eigenvalues by solving for the roots of the polynomial.
595  Scalar rho = sqrt(a_over_3);
596  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]
597  Scalar cos_theta = cos(theta);
598  Scalar sin_theta = sin(theta);
599  // roots are already sorted, since cos is monotonically decreasing on [0, pi]
600  roots(0) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta); // == 2*rho*cos(theta+2pi/3)
601  roots(1) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta); // == 2*rho*cos(theta+ pi/3)
602  roots(2) = c2_over_3 + Scalar(2)*rho*cos_theta;
603  }
604 
605  EIGEN_DEVICE_FUNC
606  static inline bool extract_kernel(MatrixType& mat, Ref<VectorType> res, Ref<VectorType> representative)
607  {
608  EIGEN_USING_STD_MATH(sqrt)
609  EIGEN_USING_STD_MATH(abs)
610  Index i0;
611  // Find non-zero column i0 (by construction, there must exist a non zero coefficient on the diagonal):
612  mat.diagonal().cwiseAbs().maxCoeff(&i0);
613  // mat.col(i0) is a good candidate for an orthogonal vector to the current eigenvector,
614  // so let's save it:
615  representative = mat.col(i0);
616  Scalar n0, n1;
617  VectorType c0, c1;
618  n0 = (c0 = representative.cross(mat.col((i0+1)%3))).squaredNorm();
619  n1 = (c1 = representative.cross(mat.col((i0+2)%3))).squaredNorm();
620  if(n0>n1) res = c0/sqrt(n0);
621  else res = c1/sqrt(n1);
622 
623  return true;
624  }
625 
626  EIGEN_DEVICE_FUNC
627  static inline void run(SolverType& solver, const MatrixType& mat, int options)
628  {
629  eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows());
630  eigen_assert((options&~(EigVecMask|GenEigMask))==0
631  && (options&EigVecMask)!=EigVecMask
632  && "invalid option parameter");
633  bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
634 
635  EigenvectorsType& eivecs = solver.m_eivec;
636  VectorType& eivals = solver.m_eivalues;
637 
638  // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow.
639  Scalar shift = mat.trace() / Scalar(3);
640  // TODO Avoid this copy. Currently it is necessary to suppress bogus values when determining maxCoeff and for computing the eigenvectors later
641  MatrixType scaledMat = mat.template selfadjointView<Lower>();
642  scaledMat.diagonal().array() -= shift;
643  Scalar scale = scaledMat.cwiseAbs().maxCoeff();
644  if(scale > 0) scaledMat /= scale; // TODO for scale==0 we could save the remaining operations
645 
646  // compute the eigenvalues
647  computeRoots(scaledMat,eivals);
648 
649  // compute the eigenvectors
650  if(computeEigenvectors)
651  {
652  if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::epsilon())
653  {
654  // All three eigenvalues are numerically the same
655  eivecs.setIdentity();
656  }
657  else
658  {
659  MatrixType tmp;
660  tmp = scaledMat;
661 
662  // Compute the eigenvector of the most distinct eigenvalue
663  Scalar d0 = eivals(2) - eivals(1);
664  Scalar d1 = eivals(1) - eivals(0);
665  Index k(0), l(2);
666  if(d0 > d1)
667  {
668  numext::swap(k,l);
669  d0 = d1;
670  }
671 
672  // Compute the eigenvector of index k
673  {
674  tmp.diagonal().array () -= eivals(k);
675  // By construction, 'tmp' is of rank 2, and its kernel corresponds to the respective eigenvector.
676  extract_kernel(tmp, eivecs.col(k), eivecs.col(l));
677  }
678 
679  // Compute eigenvector of index l
680  if(d0<=2*Eigen::NumTraits<Scalar>::epsilon()*d1)
681  {
682  // If d0 is too small, then the two other eigenvalues are numerically the same,
683  // and thus we only have to ortho-normalize the near orthogonal vector we saved above.
684  eivecs.col(l) -= eivecs.col(k).dot(eivecs.col(l))*eivecs.col(l);
685  eivecs.col(l).normalize();
686  }
687  else
688  {
689  tmp = scaledMat;
690  tmp.diagonal().array () -= eivals(l);
691 
692  VectorType dummy;
693  extract_kernel(tmp, eivecs.col(l), dummy);
694  }
695 
696  // Compute last eigenvector from the other two
697  eivecs.col(1) = eivecs.col(2).cross(eivecs.col(0)).normalized();
698  }
699  }
700 
701  // Rescale back to the original size.
702  eivals *= scale;
703  eivals.array() += shift;
704 
705  solver.m_info = Success;
706  solver.m_isInitialized = true;
707  solver.m_eigenvectorsOk = computeEigenvectors;
708  }
709 };
710 
711 // 2x2 direct eigenvalues decomposition, code from Hauke Heibel
712 template<typename SolverType>
713 struct direct_selfadjoint_eigenvalues<SolverType,2,false>
714 {
715  typedef typename SolverType::MatrixType MatrixType;
716  typedef typename SolverType::RealVectorType VectorType;
717  typedef typename SolverType::Scalar Scalar;
718  typedef typename SolverType::EigenvectorsType EigenvectorsType;
719 
720  EIGEN_DEVICE_FUNC
721  static inline void computeRoots(const MatrixType& m, VectorType& roots)
722  {
723  using std::sqrt;
724  const Scalar t0 = Scalar(0.5) * sqrt( numext::abs2(m(0,0)-m(1,1)) + Scalar(4)*numext::abs2(m(1,0)));
725  const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1));
726  roots(0) = t1 - t0;
727  roots(1) = t1 + t0;
728  }
729 
730  EIGEN_DEVICE_FUNC
731  static inline void run(SolverType& solver, const MatrixType& mat, int options)
732  {
733  EIGEN_USING_STD_MATH(sqrt);
734  EIGEN_USING_STD_MATH(abs);
735 
736  eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
737  eigen_assert((options&~(EigVecMask|GenEigMask))==0
738  && (options&EigVecMask)!=EigVecMask
739  && "invalid option parameter");
740  bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
741 
742  EigenvectorsType& eivecs = solver.m_eivec;
743  VectorType& eivals = solver.m_eivalues;
744 
745  // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow.
746  Scalar shift = mat.trace() / Scalar(2);
747  MatrixType scaledMat = mat;
748  scaledMat.coeffRef(0,1) = mat.coeff(1,0);
749  scaledMat.diagonal().array() -= shift;
750  Scalar scale = scaledMat.cwiseAbs().maxCoeff();
751  if(scale > Scalar(0))
752  scaledMat /= scale;
753 
754  // Compute the eigenvalues
755  computeRoots(scaledMat,eivals);
756 
757  // compute the eigen vectors
758  if(computeEigenvectors)
759  {
760  if((eivals(1)-eivals(0))<=abs(eivals(1))*Eigen::NumTraits<Scalar>::epsilon())
761  {
762  eivecs.setIdentity();
763  }
764  else
765  {
766  scaledMat.diagonal().array () -= eivals(1);
767  Scalar a2 = numext::abs2(scaledMat(0,0));
768  Scalar c2 = numext::abs2(scaledMat(1,1));
769  Scalar b2 = numext::abs2(scaledMat(1,0));
770  if(a2>c2)
771  {
772  eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0);
773  eivecs.col(1) /= sqrt(a2+b2);
774  }
775  else
776  {
777  eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0);
778  eivecs.col(1) /= sqrt(c2+b2);
779  }
780 
781  eivecs.col(0) << eivecs.col(1).unitOrthogonal();
782  }
783  }
784 
785  // Rescale back to the original size.
786  eivals *= scale;
787  eivals.array() += shift;
788 
789  solver.m_info = Success;
790  solver.m_isInitialized = true;
791  solver.m_eigenvectorsOk = computeEigenvectors;
792  }
793 };
794 
795 }
796 
797 template<typename MatrixType>
798 EIGEN_DEVICE_FUNC
799 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
800 ::computeDirect(const MatrixType& matrix, int options)
801 {
802  internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>::run(*this,matrix,options);
803  return *this;
804 }
805 
806 namespace internal {
807 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
808 EIGEN_DEVICE_FUNC
809 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
810 {
811  using std::abs;
812  RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
813  RealScalar e = subdiag[end-1];
814  // Note that thanks to scaling, e^2 or td^2 cannot overflow, however they can still
815  // underflow thus leading to inf/NaN values when using the following commented code:
816 // RealScalar e2 = numext::abs2(subdiag[end-1]);
817 // RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
818  // This explain the following, somewhat more complicated, version:
819  RealScalar mu = diag[end];
820  if(td==RealScalar(0))
821  mu -= abs(e);
822  else
823  {
824  RealScalar e2 = numext::abs2(subdiag[end-1]);
825  RealScalar h = numext::hypot(td,e);
826  if(e2==RealScalar(0)) mu -= (e / (td + (td>RealScalar(0) ? RealScalar(1) : RealScalar(-1)))) * (e / h);
827  else mu -= e2 / (td + (td>RealScalar(0) ? h : -h));
828  }
829 
830  RealScalar x = diag[start] - mu;
831  RealScalar z = subdiag[start];
832  for (Index k = start; k < end; ++k)
833  {
834  JacobiRotation<RealScalar> rot;
835  rot.makeGivens(x, z);
836 
837  // do T = G' T G
838  RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
839  RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
840 
841  diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
842  diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
843  subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
844 
845 
846  if (k > start)
847  subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
848 
849  x = subdiag[k];
850 
851  if (k < end - 1)
852  {
853  z = -rot.s() * subdiag[k+1];
854  subdiag[k + 1] = rot.c() * subdiag[k+1];
855  }
856 
857  // apply the givens rotation to the unit matrix Q = Q * G
858  if (matrixQ)
859  {
860  // FIXME if StorageOrder == RowMajor this operation is not very efficient
861  Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n);
862  q.applyOnTheRight(k,k+1,rot);
863  }
864  }
865 }
866 
867 } // end namespace internal
868 
869 } // end namespace Eigen
870 
871 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H
Eigen
Namespace containing all symbols from the Eigen library.
Definition: Core:309
Eigen::SelfAdjointEigenSolver::SelfAdjointEigenSolver
SelfAdjointEigenSolver(const EigenBase< InputType > &matrix, int options=ComputeEigenvectors)
Constructor; computes eigendecomposition of given matrix.
Definition: SelfAdjointEigenSolver.h:162
Eigen::Tridiagonalization
Tridiagonal decomposition of a selfadjoint matrix.
Definition: Tridiagonalization.h:64
Eigen::sqrt
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Eigen::SelfAdjointEigenSolver::computeFromTridiagonal
SelfAdjointEigenSolver & computeFromTridiagonal(const RealVectorType &diag, const SubDiagonalType &subdiag, int options=ComputeEigenvectors)
Computes the eigen decomposition from a tridiagonal symmetric matrix.
Definition: SelfAdjointEigenSolver.h:451
Eigen::EigenBase
Definition: EigenBase.h:30
Eigen::SelfAdjointEigenSolver::eigenvalues
const RealVectorType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: SelfAdjointEigenSolver.h:282
Eigen::sin
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sin_op< typename Derived::Scalar >, const Derived > sin(const Eigen::ArrayBase< Derived > &x)
Eigen::SelfAdjointEigenSolver::RealScalar
NumTraits< Scalar >::Real RealScalar
Real scalar type for _MatrixType.
Definition: SelfAdjointEigenSolver.h:94
Eigen::Success
@ Success
Definition: Constants.h:432
Eigen::ComputeEigenvectors
@ ComputeEigenvectors
Definition: Constants.h:395
Eigen::SelfAdjointEigenSolver::RealVectorType
internal::plain_col_type< MatrixType, RealScalar >::type RealVectorType
Type for vector of eigenvalues as returned by eigenvalues().
Definition: SelfAdjointEigenSolver.h:103
Eigen::SelfAdjointEigenSolver::operatorInverseSqrt
MatrixType operatorInverseSqrt() const
Computes the inverse square root of the matrix.
Definition: SelfAdjointEigenSolver.h:331
Eigen::NoConvergence
@ NoConvergence
Definition: Constants.h:436
Eigen::SelfAdjointEigenSolver::SelfAdjointEigenSolver
SelfAdjointEigenSolver(Index size)
Constructor, pre-allocates memory for dynamic-size matrices.
Definition: SelfAdjointEigenSolver.h:138
Eigen::SelfAdjointEigenSolver
Computes eigenvalues and eigenvectors of selfadjoint matrices.
Definition: SelfAdjointEigenSolver.h:71
Eigen::SelfAdjointEigenSolver::computeDirect
SelfAdjointEigenSolver & computeDirect(const MatrixType &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix using a closed-form algorithm.
Definition: SelfAdjointEigenSolver.h:800
Eigen::SelfAdjointEigenSolver::Index
Eigen::Index Index
Definition: SelfAdjointEigenSolver.h:84
Eigen::abs
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Eigen::cos
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_cos_op< typename Derived::Scalar >, const Derived > cos(const Eigen::ArrayBase< Derived > &x)
Eigen::SelfAdjointEigenSolver::eigenvectors
const EigenvectorsType & eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: SelfAdjointEigenSolver.h:259
Eigen::EigenBase::derived
Derived & derived()
Definition: EigenBase.h:45
Eigen::SelfAdjointEigenSolver::operatorSqrt
MatrixType operatorSqrt() const
Computes the positive-definite square root of the matrix.
Definition: SelfAdjointEigenSolver.h:306
Eigen::Matrix< Scalar, Size, Size, ColMajor, MaxColsAtCompileTime, MaxColsAtCompileTime >
Eigen::SelfAdjointEigenSolver::compute
SelfAdjointEigenSolver & compute(const EigenBase< InputType > &matrix, int options=ComputeEigenvectors)
Computes eigendecomposition of given matrix.
Eigen::SelfAdjointEigenSolver::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SelfAdjointEigenSolver.h:343
Eigen::SelfAdjointEigenSolver::m_maxIterations
static const int m_maxIterations
Maximum number of iterations.
Definition: SelfAdjointEigenSolver.h:354
Eigen::ComputationInfo
ComputationInfo
Definition: Constants.h:430
Eigen::MatrixBase::setIdentity
Derived & setIdentity()
Definition: CwiseNullaryOp.h:774
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:151
Eigen::SelfAdjointEigenSolver::Scalar
MatrixType::Scalar Scalar
Scalar type for matrices of type _MatrixType.
Definition: SelfAdjointEigenSolver.h:83
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33