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