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 
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()),
129  m_schur(matrix.rows()),
130  m_isInitialized(false),
131  m_eigenvectorsOk(false),
133  {
134  compute(matrix.derived(), computeEigenvectors);
135  }
136 
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 
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 
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 {
263  check_template_parameters();
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)
276  doComputeEigenvectors(m_schur.matrixT().norm());
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();
305  ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
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;
337  std::swap(m_eivalues[k],m_eivalues[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
Eigen::ComplexEigenSolver::ComplexEigenSolver
ComplexEigenSolver()
Default constructor.
Definition: ComplexEigenSolver.h:92
Eigen::ComplexEigenSolver::EigenvectorType
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorType
Type for matrix of eigenvectors as returned by eigenvectors().
Definition: ComplexEigenSolver.h:85
Eigen::ComplexEigenSolver::MaxColsAtCompileTime
@ MaxColsAtCompileTime
Definition: ComplexEigenSolver.h:57
Eigen::ComplexEigenSolver::ComplexEigenSolver
ComplexEigenSolver(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Constructor; computes eigendecomposition of given matrix.
Definition: ComplexEigenSolver.h:126
Eigen::ComplexEigenSolver::m_eivalues
EigenvalueType m_eivalues
Definition: ComplexEigenSolver.h:246
Eigen::ComplexEigenSolver::doComputeEigenvectors
void doComputeEigenvectors(RealScalar matrixnorm)
Definition: ComplexEigenSolver.h:287
Eigen
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Eigen::ComplexSchur::setMaxIterations
ComplexSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: ComplexSchur.h:228
Eigen::ComplexEigenSolver::m_isInitialized
bool m_isInitialized
Definition: ComplexEigenSolver.h:248
Eigen::ComplexEigenSolver::m_schur
ComplexSchur< MatrixType > m_schur
Definition: ComplexEigenSolver.h:247
Eigen::ComplexEigenSolver::ComplexScalar
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Definition: ComplexEigenSolver.h:71
Eigen::EigenBase
Definition: EigenBase.h:29
eigen_assert
#define eigen_assert(x)
Definition: Macros.h:1037
Eigen::ComplexEigenSolver::setMaxIterations
ComplexEigenSolver & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: ComplexEigenSolver.h:226
Eigen::ComplexEigenSolver::m_eigenvectorsOk
bool m_eigenvectorsOk
Definition: ComplexEigenSolver.h:249
Eigen::Success
@ Success
Definition: Constants.h:442
Eigen::ComplexEigenSolver::ColsAtCompileTime
@ ColsAtCompileTime
Definition: ComplexEigenSolver.h:54
Eigen::ComplexEigenSolver::MaxRowsAtCompileTime
@ MaxRowsAtCompileTime
Definition: ComplexEigenSolver.h:56
Eigen::RowMajor
@ RowMajor
Definition: Constants.h:321
Eigen::ComplexEigenSolver::check_template_parameters
static void check_template_parameters()
Definition: ComplexEigenSolver.h:240
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
size
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
Eigen::ComplexSchur::getMaxIterations
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: ComplexSchur.h:235
Eigen::ComplexEigenSolver::Options
@ Options
Definition: ComplexEigenSolver.h:55
Eigen::ComplexEigenSolver::eigenvalues
const EigenvalueType & eigenvalues() const
Returns the eigenvalues of given matrix.
Definition: ComplexEigenSolver.h:182
Eigen::ComplexSchur::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: ComplexSchur.h:217
n
int n
Definition: BiCGSTAB_simple.cpp:1
epsilon
static double epsilon
Definition: testRot3.cpp:37
Eigen::ComplexEigenSolver::RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: ComplexEigenSolver.h:62
Eigen::numext::real_ref
EIGEN_DEVICE_FUNC internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) >::type real_ref(const Scalar &x)
Definition: Eigen/src/Core/MathFunctions.h:1239
ComplexSchur.h
Eigen::ComplexEigenSolver::m_matX
EigenvectorType m_matX
Definition: ComplexEigenSolver.h:250
std::swap
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
Definition: NearestNeighbor.hpp:827
Eigen::ComplexEigenSolver::Index
Eigen::Index Index
Definition: ComplexEigenSolver.h:63
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
Eigen::ComplexEigenSolver::ComplexEigenSolver
ComplexEigenSolver(Index size)
Default Constructor with memory preallocation.
Definition: ComplexEigenSolver.h:107
matrix
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
Definition: gtsam/3rdparty/Eigen/blas/common.h:110
Eigen::ComplexEigenSolver::MatrixType
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
Definition: ComplexEigenSolver.h:50
Eigen::ComplexEigenSolver::eigenvectors
const EigenvectorType & eigenvectors() const
Returns the eigenvectors of given matrix.
Definition: ComplexEigenSolver.h:157
Eigen::ComplexEigenSolver::getMaxIterations
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: ComplexEigenSolver.h:233
Eigen::ComplexEigenSolver::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: ComplexEigenSolver.h:219
Eigen::ComplexEigenSolver::Scalar
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: ComplexEigenSolver.h:61
Eigen::ComplexEigenSolver
Computes eigenvalues and eigenvectors of general complex matrices.
Definition: ComplexEigenSolver.h:45
Eigen::ComplexEigenSolver::compute
ComplexEigenSolver & compute(const EigenBase< InputType > &matrix, bool computeEigenvectors=true)
Computes eigendecomposition of given matrix.
min
#define min(a, b)
Definition: datatypes.h:19
Eigen::ComplexEigenSolver::sortEigenvalues
void sortEigenvalues(bool computeEigenvectors)
Definition: ComplexEigenSolver.h:327
Eigen::ComplexEigenSolver::m_eivec
EigenvectorType m_eivec
Definition: ComplexEigenSolver.h:245
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition: 3rdparty/Eigen/Eigen/src/Core/Matrix.h:178
Eigen::numext::maxi
EIGEN_DEVICE_FUNC EIGEN_ALWAYS_INLINE T maxi(const T &x, const T &y)
Definition: Eigen/src/Core/MathFunctions.h:1093
cols
int cols
Definition: Tutorial_commainit_02.cpp:1
Eigen::ComputationInfo
ComputationInfo
Definition: Constants.h:440
Eigen::ComplexEigenSolver::RowsAtCompileTime
@ RowsAtCompileTime
Definition: ComplexEigenSolver.h:53
Eigen::ComplexSchur< MatrixType >
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:232
test_callbacks.value
value
Definition: test_callbacks.py:160
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
EIGEN_STATIC_ASSERT_NON_INTEGER
#define EIGEN_STATIC_ASSERT_NON_INTEGER(TYPE)
Definition: StaticAssert.h:187


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:02:00