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