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 typename MatrixType::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  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  ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
126  : m_eivec(matrix.rows(),matrix.cols()),
127  m_eivalues(matrix.cols()),
128  m_schur(matrix.rows()),
129  m_isInitialized(false),
130  m_eigenvectorsOk(false),
131  m_matX(matrix.rows(),matrix.cols())
132  {
133  compute(matrix, computeEigenvectors);
134  }
135 
156  const EigenvectorType& eigenvectors() const
157  {
158  eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
159  eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
160  return m_eivec;
161  }
162 
181  const EigenvalueType& eigenvalues() const
182  {
183  eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
184  return m_eivalues;
185  }
186 
211  ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true);
212 
218  {
219  eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized.");
220  return m_schur.info();
221  }
222 
225  {
226  m_schur.setMaxIterations(maxIters);
227  return *this;
228  }
229 
232  {
233  return m_schur.getMaxIterations();
234  }
235 
236  protected:
237  EigenvectorType m_eivec;
238  EigenvalueType m_eivalues;
242  EigenvectorType m_matX;
243 
244  private:
245  void doComputeEigenvectors(const RealScalar& matrixnorm);
246  void sortEigenvalues(bool computeEigenvectors);
247 };
248 
249 
250 template<typename MatrixType>
252 ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
253 {
254  // this code is inspired from Jampack
255  eigen_assert(matrix.cols() == matrix.rows());
256 
257  // Do a complex Schur decomposition, A = U T U^*
258  // The eigenvalues are on the diagonal of T.
259  m_schur.compute(matrix, computeEigenvectors);
260 
261  if(m_schur.info() == Success)
262  {
263  m_eivalues = m_schur.matrixT().diagonal();
264  if(computeEigenvectors)
265  doComputeEigenvectors(matrix.norm());
266  sortEigenvalues(computeEigenvectors);
267  }
268 
269  m_isInitialized = true;
270  m_eigenvectorsOk = computeEigenvectors;
271  return *this;
272 }
273 
274 
275 template<typename MatrixType>
277 {
278  const Index n = m_eivalues.size();
279 
280  // Compute X such that T = X D X^(-1), where D is the diagonal of T.
281  // The matrix X is unit triangular.
282  m_matX = EigenvectorType::Zero(n, n);
283  for(Index k=n-1 ; k>=0 ; k--)
284  {
285  m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0);
286  // Compute X(i,k) using the (i,k) entry of the equation X T = D X
287  for(Index i=k-1 ; i>=0 ; i--)
288  {
289  m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k);
290  if(k-i-1>0)
291  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();
292  ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k);
293  if(z==ComplexScalar(0))
294  {
295  // If the i-th and k-th eigenvalue are equal, then z equals 0.
296  // Use a small value instead, to prevent division by zero.
298  }
299  m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z;
300  }
301  }
302 
303  // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1)
304  m_eivec.noalias() = m_schur.matrixU() * m_matX;
305  // .. and normalize the eigenvectors
306  for(Index k=0 ; k<n ; k++)
307  {
308  m_eivec.col(k).normalize();
309  }
310 }
311 
312 
313 template<typename MatrixType>
315 {
316  const Index n = m_eivalues.size();
317  for (Index i=0; i<n; i++)
318  {
319  Index k;
320  m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k);
321  if (k != 0)
322  {
323  k += i;
324  std::swap(m_eivalues[k],m_eivalues[i]);
325  if(computeEigenvectors)
326  m_eivec.col(i).swap(m_eivec.col(k));
327  }
328  }
329 }
330 
331 } // end namespace Eigen
332 
333 #endif // EIGEN_COMPLEX_EIGEN_SOLVER_H
ComplexSchur< MatrixType > m_schur
void swap(MatrixBase< OtherDerived > const &other)
Definition: Matrix.h:318
ComplexSchur & setMaxIterations(Index maxIters)
Sets the maximum number of iterations allowed.
Definition: ComplexSchur.h:226
ComplexEigenSolver(Index size)
Default Constructor with memory preallocation.
ComplexSchur & compute(const MatrixType &matrix, bool computeU=true)
Computes Schur decomposition of given matrix.
Definition: ComplexSchur.h:316
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &(~RowMajor), MaxColsAtCompileTime, 1 > EigenvalueType
Type for vector of eigenvalues as returned by eigenvalues().
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorType
Type for matrix of eigenvectors as returned by eigenvectors().
ComputationInfo info() const
Reports whether previous computation was successful.
ComplexEigenSolver(const MatrixType &matrix, bool computeEigenvectors=true)
Constructor; computes eigendecomposition of given matrix.
NumTraits< Scalar >::Real RealScalar
EIGEN_STRONG_INLINE const Scalar & coeff(Index rowId, Index colId) const
EIGEN_STRONG_INLINE Scalar & coeffRef(Index rowId, Index colId)
Provides a generic way to set and pass user-specified options.
Definition: options.hpp:65
Index getMaxIterations()
Returns the maximum number of iterations.
Definition: ComplexSchur.h:233
void doComputeEigenvectors(const RealScalar &matrixnorm)
const EigenvectorType & eigenvectors() const
Returns the eigenvectors of given matrix.
const ComplexMatrixType & matrixT() const
Returns the triangular matrix in the Schur decomposition.
Definition: ComplexSchur.h:161
_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.
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: ComplexSchur.h:215
void sortEigenvalues(bool computeEigenvectors)
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) >::type real_ref(const Scalar &x)
const EigenvalueType & eigenvalues() const
Returns the eigenvalues of given matrix.
#define eigen_assert(x)
Computes eigenvalues and eigenvectors of general complex matrices.
ComplexEigenSolver & compute(const MatrixType &matrix, bool computeEigenvectors=true)
Computes eigendecomposition of given matrix.
const ComplexMatrixType & matrixU() const
Returns the unitary matrix in the Schur decomposition.
Definition: ComplexSchur.h:137


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:31