ComplexEigenSolver.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) 2009 Claire Maurice
5 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
6 // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
7 //
8 // This Source Code Form is subject to the terms of the Mozilla
9 // Public License v. 2.0. If a copy of the MPL was not distributed
10 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
11 
12 #ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H
13 #define EIGEN_COMPLEX_EIGEN_SOLVER_H
14 
15 #include "./ComplexSchur.h"
16 
17 namespace Eigen {
18 
45 template<typename _MatrixType> class ComplexEigenSolver
46 {
47  public:
48 
50  typedef _MatrixType MatrixType;
51 
52  enum {
53  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
54  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
55  Options = MatrixType::Options,
56  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
57  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
58  };
59 
61  typedef typename MatrixType::Scalar Scalar;
63  typedef Eigen::Index Index;
64 
71  typedef std::complex<RealScalar> ComplexScalar;
72 
79 
86 
93  : m_eivec(),
94  m_eivalues(),
95  m_schur(),
96  m_isInitialized(false),
97  m_eigenvectorsOk(false),
98  m_matX()
99  {}
100 
107  explicit ComplexEigenSolver(Index size)
108  : m_eivec(size, size),
109  m_eivalues(size),
110  m_schur(size),
111  m_isInitialized(false),
112  m_eigenvectorsOk(false),
113  m_matX(size, size)
114  {}
115 
125  template<typename InputType>
126  explicit ComplexEigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true)
127  : m_eivec(matrix.rows(),matrix.cols()),
128  m_eivalues(matrix.cols()),
129  m_schur(matrix.rows()),
130  m_isInitialized(false),
131  m_eigenvectorsOk(false),
132  m_matX(matrix.rows(),matrix.cols())
133  {
134  compute(matrix.derived(), computeEigenvectors);
135  }
136 
157  const EigenvectorType& eigenvectors() const
158  {
159  eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
160  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
161  return m_eivec;
162  }
163 
182  const EigenvalueType& eigenvalues() const
183  {
184  eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
185  return m_eivalues;
186  }
187 
212  template<typename InputType>
213  ComplexEigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true);
214 
220  {
221  eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
222  return m_schur.info();
223  }
224 
227  {
228  m_schur.setMaxIterations(maxIters);
229  return *this;
230  }
231 
234  {
235  return m_schur.getMaxIterations();
236  }
237 
238  protected:
239 
241  {
243  }
244 
245  EigenvectorType m_eivec;
246  EigenvalueType m_eivalues;
250  EigenvectorType m_matX;
251 
252  private:
253  void doComputeEigenvectors(RealScalar matrixnorm);
254  void sortEigenvalues(bool computeEigenvectors);
255 };
256 
257 
258 template<typename MatrixType>
259 template<typename InputType>
262 {
264 
265  // this code is inspired from Jampack
266  eigen_assert(matrix.cols() == matrix.rows());
267 
268  // Do a complex Schur decomposition, A = U T U^*
269  // The eigenvalues are on the diagonal of T.
270  m_schur.compute(matrix.derived(), computeEigenvectors);
271 
272  if(m_schur.info() == Success)
273  {
274  m_eivalues = m_schur.matrixT().diagonal();
275  if(computeEigenvectors)
277  sortEigenvalues(computeEigenvectors);
278  }
279 
280  m_isInitialized = true;
281  m_eigenvectorsOk = computeEigenvectors;
282  return *this;
283 }
284 
285 
286 template<typename MatrixType>
288 {
289  const Index n = m_eivalues.size();
290 
291  matrixnorm = numext::maxi(matrixnorm,(std::numeric_limits<RealScalar>::min)());
292 
293  // Compute X such that T = X D X^(-1), where D is the diagonal of T.
294  // The matrix X is unit triangular.
295  m_matX = EigenvectorType::Zero(n, n);
296  for(Index k=n-1 ; k>=0 ; k--)
297  {
298  m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
299  // Compute X(i,k) using the (i,k) entry of the equation X T = D X
300  for(Index i=k-1 ; i>=0 ; i--)
301  {
302  m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
303  if(k-i-1>0)
304  m_matX.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * m_matX.col(k).segment(i+1,k-i-1)).value();
306  if(z==ComplexScalar(0))
307  {
308  // If the i-th and k-th eigenvalue are equal, then z equals 0.
309  // Use a small value instead, to prevent division by zero.
311  }
312  m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
313  }
314  }
315 
316  // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
317  m_eivec.noalias() = m_schur.matrixU() * m_matX;
318  // .. and normalize the eigenvectors
319  for(Index k=0 ; k<n ; k++)
320  {
321  m_eivec.col(k).normalize();
322  }
323 }
324 
325 
326 template<typename MatrixType>
328 {
329  const Index n = m_eivalues.size();
330  for (Index i=0; i<n; i++)
331  {
332  Index k;
333  m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
334  if (k != 0)
335  {
336  k += i;
338  if(computeEigenvectors)
339  m_eivec.col(i).swap(m_eivec.col(k));
340  }
341  }
342 }
343 
344 } // end namespace Eigen
345 
346 #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H
SCALAR Scalar
Definition: bench_gemm.cpp:46
EIGEN_DEVICE_FUNC internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) >::type real_ref(const Scalar &x)
static void check_template_parameters()
ComplexSchur< MatrixType > m_schur
ComplexSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: ComplexSchur.h:228
ComplexEigenSolver(Index size)
Default Constructor with memory preallocation.
#define min(a, b)
Definition: datatypes.h:19
ComplexEigenSolver & compute(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Computes eigendecomposition of given matrix.
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &(~RowMajor), MaxColsAtCompileTime, 1 > EigenvalueType
Type for vector of eigenvalues as returned by eigenvalues().
int n
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:63
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void swap(DenseBase< OtherDerived > &other)
void doComputeEigenvectors(RealScalar matrixnorm)
ComplexSchur & compute(const EigenBase< InputType > &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: ComplexSchur.h:217
static double epsilon
Definition: testRot3.cpp:37
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar & coeffRef(Index rowId, Index colId)
const EigenvalueType & eigenvalues() const
Returns the eigenvalues of given matrix.
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorType
Type for matrix of eigenvectors as returned by eigenvectors().
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
NumTraits< Scalar >::Real RealScalar
ComputationInfo info() const
Reports whether previous computation was successful.
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
#define eigen_assert(x)
Definition: Macros.h:1037
ComplexEigenSolver(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Constructor; computes eigendecomposition of given matrix.
#define EIGEN_STATIC_ASSERT_NON_INTEGER(TYPE)
Definition: StaticAssert.h:187
const EigenvectorType & eigenvectors() const
Returns the eigenvectors of given matrix.
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: ComplexSchur.h:235
EIGEN_DEVICE_FUNC EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
ComplexEigenSolver()
Default constructor.
Index getMaxIterations()
Returns the maximum number of iterations.
ComplexEigenSolver & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
const ComplexMatrixType & matrixU() const
Returns the unitary matrix in the Schur decomposition.
Definition: ComplexSchur.h:138
void sortEigenvalues(bool computeEigenvectors)
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
Computes eigenvalues and eigenvectors of general complex matrices.
ComputationInfo
Definition: Constants.h:440
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:46
const ComplexMatrixType & matrixT() const
Returns the triangular matrix in the Schur decomposition.
Definition: ComplexSchur.h:162


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