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