EigenSolver.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2010,2012 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_EIGENSOLVER_H
12 #define EIGEN_EIGENSOLVER_H
13 
14 #include "./RealSchur.h"
15 
16 namespace Eigen {
17 
64 template<typename _MatrixType> class EigenSolver
65 {
66  public:
67 
69  typedef _MatrixType MatrixType;
70 
71  enum {
72  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
73  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
74  Options = MatrixType::Options,
75  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
76  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
77  };
78 
80  typedef typename MatrixType::Scalar Scalar;
82  typedef Eigen::Index Index;
83 
90  typedef std::complex<RealScalar> ComplexScalar;
91 
98 
105 
114 
121  explicit EigenSolver(Index size)
122  : m_eivec(size, size),
123  m_eivalues(size),
124  m_isInitialized(false),
125  m_eigenvectorsOk(false),
126  m_realSchur(size),
127  m_matT(size, size),
128  m_tmp(size)
129  {}
130 
146  template<typename InputType>
147  explicit EigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true)
148  : m_eivec(matrix.rows(), matrix.cols()),
149  m_eivalues(matrix.cols()),
150  m_isInitialized(false),
151  m_eigenvectorsOk(false),
152  m_realSchur(matrix.cols()),
153  m_matT(matrix.rows(), matrix.cols()),
154  m_tmp(matrix.cols())
155  {
156  compute(matrix.derived(), computeEigenvectors);
157  }
158 
179  EigenvectorsType eigenvectors() const;
180 
199  const MatrixType& pseudoEigenvectors() const
200  {
201  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
202  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
203  return m_eivec;
204  }
205 
224  MatrixType pseudoEigenvalueMatrix() const;
225 
244  const EigenvalueType& eigenvalues() const
245  {
246  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
247  return m_eivalues;
248  }
249 
277  template<typename InputType>
278  EigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true);
279 
282  {
283  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
284  return m_info;
285  }
286 
288  EigenSolver& setMaxIterations(Index maxIters)
289  {
290  m_realSchur.setMaxIterations(maxIters);
291  return *this;
292  }
293 
296  {
297  return m_realSchur.getMaxIterations();
298  }
299 
300  private:
301  void doComputeEigenvectors();
302 
303  protected:
304 
306  {
308  EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
309  }
310 
311  MatrixType m_eivec;
312  EigenvalueType m_eivalues;
317  MatrixType m_matT;
318 
320  ColumnVectorType m_tmp;
321 };
322 
323 template<typename MatrixType>
325 {
326  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
328  Index n = m_eivalues.rows();
329  MatrixType matD = MatrixType::Zero(n,n);
330  for (Index i=0; i<n; ++i)
331  {
333  matD.coeffRef(i,i) = numext::real(m_eivalues.coeff(i));
334  else
335  {
336  matD.template block<2,2>(i,i) << numext::real(m_eivalues.coeff(i)), numext::imag(m_eivalues.coeff(i)),
338  ++i;
339  }
340  }
341  return matD;
342 }
343 
344 template<typename MatrixType>
346 {
347  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
348  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
350  Index n = m_eivec.cols();
351  EigenvectorsType matV(n,n);
352  for (Index j=0; j<n; ++j)
353  {
355  {
356  // we have a real eigen value
357  matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
358  matV.col(j).normalize();
359  }
360  else
361  {
362  // we have a pair of complex eigen values
363  for (Index i=0; i<n; ++i)
364  {
365  matV.coeffRef(i,j) = ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1));
366  matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
367  }
368  matV.col(j).normalize();
369  matV.col(j+1).normalize();
370  ++j;
371  }
372  }
373  return matV;
374 }
375 
376 template<typename MatrixType>
377 template<typename InputType>
380 {
382 
383  using std::sqrt;
384  using std::abs;
385  using numext::isfinite;
386  eigen_assert(matrix.cols() == matrix.rows());
387 
388  // Reduce to real Schur form.
389  m_realSchur.compute(matrix.derived(), computeEigenvectors);
390 
391  m_info = m_realSchur.info();
392 
393  if (m_info == Success)
394  {
396  if (computeEigenvectors)
398 
399  // Compute eigenvalues from matT
400  m_eivalues.resize(matrix.cols());
401  Index i = 0;
402  while (i < matrix.cols())
403  {
404  if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0))
405  {
406  m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
407  if(!(isfinite)(m_eivalues.coeffRef(i)))
408  {
409  m_isInitialized = true;
410  m_eigenvectorsOk = false;
412  return *this;
413  }
414  ++i;
415  }
416  else
417  {
418  Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
419  Scalar z;
420  // Compute z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
421  // without overflow
422  {
423  Scalar t0 = m_matT.coeff(i+1, i);
424  Scalar t1 = m_matT.coeff(i, i+1);
425  Scalar maxval = numext::maxi<Scalar>(abs(p),numext::maxi<Scalar>(abs(t0),abs(t1)));
426  t0 /= maxval;
427  t1 /= maxval;
428  Scalar p0 = p/maxval;
429  z = maxval * sqrt(abs(p0 * p0 + t0 * t1));
430  }
431 
432  m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
433  m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
434  if(!((isfinite)(m_eivalues.coeffRef(i)) && (isfinite)(m_eivalues.coeffRef(i+1))))
435  {
436  m_isInitialized = true;
437  m_eigenvectorsOk = false;
439  return *this;
440  }
441  i += 2;
442  }
443  }
444 
445  // Compute eigenvectors.
446  if (computeEigenvectors)
448  }
449 
450  m_isInitialized = true;
451  m_eigenvectorsOk = computeEigenvectors;
452 
453  return *this;
454 }
455 
456 
457 template<typename MatrixType>
459 {
460  using std::abs;
461  const Index size = m_eivec.cols();
462  const Scalar eps = NumTraits<Scalar>::epsilon();
463 
464  // inefficient! this is already computed in RealSchur
465  Scalar norm(0);
466  for (Index j = 0; j < size; ++j)
467  {
468  norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
469  }
470 
471  // Backsubstitute to find vectors of upper triangular form
472  if (norm == Scalar(0))
473  {
474  return;
475  }
476 
477  for (Index n = size-1; n >= 0; n--)
478  {
479  Scalar p = m_eivalues.coeff(n).real();
480  Scalar q = m_eivalues.coeff(n).imag();
481 
482  // Scalar vector
483  if (q == Scalar(0))
484  {
485  Scalar lastr(0), lastw(0);
486  Index l = n;
487 
488  m_matT.coeffRef(n,n) = Scalar(1);
489  for (Index i = n-1; i >= 0; i--)
490  {
491  Scalar w = m_matT.coeff(i,i) - p;
492  Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
493 
494  if (m_eivalues.coeff(i).imag() < Scalar(0))
495  {
496  lastw = w;
497  lastr = r;
498  }
499  else
500  {
501  l = i;
502  if (m_eivalues.coeff(i).imag() == Scalar(0))
503  {
504  if (w != Scalar(0))
505  m_matT.coeffRef(i,n) = -r / w;
506  else
507  m_matT.coeffRef(i,n) = -r / (eps * norm);
508  }
509  else // Solve real equations
510  {
511  Scalar x = m_matT.coeff(i,i+1);
512  Scalar y = m_matT.coeff(i+1,i);
513  Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
514  Scalar t = (x * lastr - lastw * r) / denom;
515  m_matT.coeffRef(i,n) = t;
516  if (abs(x) > abs(lastw))
517  m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
518  else
519  m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
520  }
521 
522  // Overflow control
523  Scalar t = abs(m_matT.coeff(i,n));
524  if ((eps * t) * t > Scalar(1))
525  m_matT.col(n).tail(size-i) /= t;
526  }
527  }
528  }
529  else if (q < Scalar(0) && n > 0) // Complex vector
530  {
531  Scalar lastra(0), lastsa(0), lastw(0);
532  Index l = n-1;
533 
534  // Last vector component imaginary so matrix is triangular
535  if (abs(m_matT.coeff(n,n-1)) > abs(m_matT.coeff(n-1,n)))
536  {
537  m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
538  m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
539  }
540  else
541  {
542  ComplexScalar cc = ComplexScalar(Scalar(0),-m_matT.coeff(n-1,n)) / ComplexScalar(m_matT.coeff(n-1,n-1)-p,q);
543  m_matT.coeffRef(n-1,n-1) = numext::real(cc);
544  m_matT.coeffRef(n-1,n) = numext::imag(cc);
545  }
546  m_matT.coeffRef(n,n-1) = Scalar(0);
547  m_matT.coeffRef(n,n) = Scalar(1);
548  for (Index i = n-2; i >= 0; i--)
549  {
550  Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
551  Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
552  Scalar w = m_matT.coeff(i,i) - p;
553 
554  if (m_eivalues.coeff(i).imag() < Scalar(0))
555  {
556  lastw = w;
557  lastra = ra;
558  lastsa = sa;
559  }
560  else
561  {
562  l = i;
563  if (m_eivalues.coeff(i).imag() == RealScalar(0))
564  {
565  ComplexScalar cc = ComplexScalar(-ra,-sa) / ComplexScalar(w,q);
566  m_matT.coeffRef(i,n-1) = numext::real(cc);
567  m_matT.coeffRef(i,n) = numext::imag(cc);
568  }
569  else
570  {
571  // Solve complex equations
572  Scalar x = m_matT.coeff(i,i+1);
573  Scalar y = m_matT.coeff(i+1,i);
574  Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
575  Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
576  if ((vr == Scalar(0)) && (vi == Scalar(0)))
577  vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw));
578 
579  ComplexScalar cc = ComplexScalar(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra) / ComplexScalar(vr,vi);
580  m_matT.coeffRef(i,n-1) = numext::real(cc);
581  m_matT.coeffRef(i,n) = numext::imag(cc);
582  if (abs(x) > (abs(lastw) + abs(q)))
583  {
584  m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
585  m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
586  }
587  else
588  {
589  cc = ComplexScalar(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n)) / ComplexScalar(lastw,q);
590  m_matT.coeffRef(i+1,n-1) = numext::real(cc);
591  m_matT.coeffRef(i+1,n) = numext::imag(cc);
592  }
593  }
594 
595  // Overflow control
596  Scalar t = numext::maxi<Scalar>(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n)));
597  if ((eps * t) * t > Scalar(1))
598  m_matT.block(i, n-1, size-i, 2) /= t;
599 
600  }
601  }
602 
603  // We handled a pair of complex conjugate eigenvalues, so need to skip them both
604  n--;
605  }
606  else
607  {
608  eigen_assert(0 && "Internal bug in EigenSolver (INF or NaN has not been detected)"); // this should not happen
609  }
610  }
611 
612  // Back transformation to get eigenvectors of original matrix
613  for (Index j = size-1; j >= 0; j--)
614  {
615  m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
616  m_eivec.col(j) = m_tmp;
617  }
618 }
619 
620 } // end namespace Eigen
621 
622 #endif // EIGEN_EIGENSOLVER_H
MatrixType m_matT
Definition: EigenSolver.h:317
MatrixType pseudoEigenvalueMatrix() const
Returns the block-diagonal matrix in the pseudo-eigendecomposition.
Definition: EigenSolver.h:324
SCALAR Scalar
Definition: bench_gemm.cpp:46
const MatrixType & pseudoEigenvectors() const
Returns the pseudo-eigenvectors of given matrix.
Definition: EigenSolver.h:199
ColumnVectorType m_tmp
Definition: EigenSolver.h:320
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorsType
Type for matrix of eigenvectors as returned by eigenvectors().
Definition: EigenSolver.h:104
#define max(a, b)
Definition: datatypes.h:20
EIGEN_DEVICE_FUNC bool isMuchSmallerThan(const Scalar &x, const OtherScalar &y, const typename NumTraits< Scalar >::Real &precision=NumTraits< Scalar >::dummy_precision())
EigenSolver & compute(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Computes eigendecomposition of given matrix.
float real
Definition: datatypes.h:10
Scalar * y
EigenSolver(Index size)
Default constructor with memory preallocation.
Definition: EigenSolver.h:121
int n
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Vector3f p0
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
#define EIGEN_STATIC_ASSERT(CONDITION, MSG)
Definition: StaticAssert.h:127
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:63
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: EigenSolver.h:295
#define isfinite(X)
Definition: main.h:95
Eigen::Index Index
Definition: EigenSolver.h:82
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
RealSchur< MatrixType > m_realSchur
Definition: EigenSolver.h:316
static double epsilon
Definition: testRot3.cpp:37
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: RealSchur.h:213
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar & coeffRef(Index rowId, Index colId)
EigenSolver(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Constructor; computes eigendecomposition of given matrix.
Definition: EigenSolver.h:147
Matrix< Scalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > ColumnVectorType
Definition: EigenSolver.h:319
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
static const Line3 l(Rot3(), 1, 1)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
NumTraits< Scalar >::Real RealScalar
Definition: EigenSolver.h:81
ComputationInfo info() const
Definition: EigenSolver.h:281
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar & coeff(Index rowId, Index colId) const
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE bool() isfinite(const Eigen::bfloat16 &h)
Definition: BFloat16.h:671
#define eigen_assert(x)
Definition: Macros.h:1037
#define EIGEN_STATIC_ASSERT_NON_INTEGER(TYPE)
Definition: StaticAssert.h:187
EIGEN_DEVICE_FUNC const Scalar & q
RealSchur & compute(const EigenBase< InputType > &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
EigenSolver & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: EigenSolver.h:288
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Definition: EigenSolver.h:90
RowVector3d w
const MatrixType & matrixT() const
Returns the quasi-triangular matrix in the Schur decomposition.
Definition: RealSchur.h:144
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
Definition: EigenSolver.h:69
EigenvectorsType eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: EigenSolver.h:345
EigenSolver()
Default constructor.
Definition: EigenSolver.h:113
const MatrixType & matrixU() const
Returns the orthogonal matrix in the Schur decomposition.
Definition: RealSchur.h:127
RealSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: RealSchur.h:206
cout precision(2)
float * p
EIGEN_DEVICE_FUNC const ImagReturnType imag() const
void doComputeEigenvectors()
Definition: EigenSolver.h:458
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: EigenSolver.h:80
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
ComputationInfo m_info
Definition: EigenSolver.h:315
const EigenvalueType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: EigenSolver.h:244
Computes eigenvalues and eigenvectors of general matrices.
Definition: EigenSolver.h:64
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealSchur.h:195
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
#define abs(x)
Definition: datatypes.h:17
ComputationInfo
Definition: Constants.h:440
EigenvalueType m_eivalues
Definition: EigenSolver.h:312
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:46
std::ptrdiff_t j
Point2 t(10, 10)
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > EigenvalueType
Type for vector of eigenvalues as returned by eigenvalues().
Definition: EigenSolver.h:97
MatrixType m_eivec
Definition: EigenSolver.h:311
static void check_template_parameters()
Definition: EigenSolver.h:305


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:34:12