GeneralizedEigenSolver.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) 2012 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_GENERALIZEDEIGENSOLVER_H
12 #define EIGEN_GENERALIZEDEIGENSOLVER_H
13 
14 #include "./RealQZ.h"
15 
16 namespace Eigen {
17 
57 template<typename _MatrixType> class GeneralizedEigenSolver
58 {
59  public:
60 
62  typedef _MatrixType MatrixType;
63 
64  enum {
65  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
66  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
67  Options = MatrixType::Options,
68  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
69  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
70  };
71 
73  typedef typename MatrixType::Scalar Scalar;
75  typedef typename MatrixType::Index Index;
76 
83  typedef std::complex<RealScalar> ComplexScalar;
84 
91 
98 
102 
109 
118 
126  : m_eivec(size, size),
127  m_alphas(size),
128  m_betas(size),
129  m_isInitialized(false),
130  m_eigenvectorsOk(false),
131  m_realQZ(size),
132  m_matS(size, size),
133  m_tmp(size)
134  {}
135 
148  GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true)
149  : m_eivec(A.rows(), A.cols()),
150  m_alphas(A.cols()),
151  m_betas(A.cols()),
152  m_isInitialized(false),
153  m_eigenvectorsOk(false),
154  m_realQZ(A.cols()),
155  m_matS(A.rows(), A.cols()),
156  m_tmp(A.cols())
157  {
158  compute(A, B, computeEigenvectors);
159  }
160 
161  /* \brief Returns the computed generalized eigenvectors.
162  *
163  * \returns %Matrix whose columns are the (possibly complex) eigenvectors.
164  *
165  * \pre Either the constructor
166  * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function
167  * compute(const MatrixType&, const MatrixType& bool) has been called before, and
168  * \p computeEigenvectors was set to true (the default).
169  *
170  * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
171  * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The
172  * eigenvectors are normalized to have (Euclidean) norm equal to one. The
173  * matrix returned by this function is the matrix \f$ V \f$ in the
174  * generalized eigendecomposition \f$ A = B V D V^{-1} \f$, if it exists.
175  *
176  * \sa eigenvalues()
177  */
178 // EigenvectorsType eigenvectors() const;
179 
198  EigenvalueType eigenvalues() const
199  {
200  eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
202  }
203 
209  ComplexVectorType alphas() const
210  {
211  eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
212  return m_alphas;
213  }
214 
220  VectorType betas() const
221  {
222  eigen_assert(m_isInitialized && "GeneralizedEigenSolver is not initialized.");
223  return m_betas;
224  }
225 
249  GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true);
250 
252  {
253  eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
254  return m_realQZ.info();
255  }
256 
260  {
261  m_realQZ.setMaxIterations(maxIters);
262  return *this;
263  }
264 
265  protected:
266  MatrixType m_eivec;
267  ComplexVectorType m_alphas;
268  VectorType m_betas;
272  MatrixType m_matS;
273 
275  ColumnVectorType m_tmp;
276 };
277 
278 //template<typename MatrixType>
279 //typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const
280 //{
281 // eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
282 // eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
283 // Index n = m_eivec.cols();
284 // EigenvectorsType matV(n,n);
285 // // TODO
286 // return matV;
287 //}
288 
289 template<typename MatrixType>
291 GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
292 {
293  using std::sqrt;
294  using std::abs;
295  eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
296 
297  // Reduce to generalized real Schur form:
298  // A = Q S Z and B = Q T Z
299  m_realQZ.compute(A, B, computeEigenvectors);
300 
301  if (m_realQZ.info() == Success)
302  {
303  m_matS = m_realQZ.matrixS();
304  if (computeEigenvectors)
305  m_eivec = m_realQZ.matrixZ().transpose();
306 
307  // Compute eigenvalues from matS
308  m_alphas.resize(A.cols());
309  m_betas.resize(A.cols());
310  Index i = 0;
311  while (i < A.cols())
312  {
313  if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0))
314  {
315  m_alphas.coeffRef(i) = m_matS.coeff(i, i);
316  m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i);
317  ++i;
318  }
319  else
320  {
321  Scalar p = Scalar(0.5) * (m_matS.coeff(i, i) - m_matS.coeff(i+1, i+1));
322  Scalar z = sqrt(abs(p * p + m_matS.coeff(i+1, i) * m_matS.coeff(i, i+1)));
323  m_alphas.coeffRef(i) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, z);
324  m_alphas.coeffRef(i+1) = ComplexScalar(m_matS.coeff(i+1, i+1) + p, -z);
325 
326  m_betas.coeffRef(i) = m_realQZ.matrixT().coeff(i,i);
327  m_betas.coeffRef(i+1) = m_realQZ.matrixT().coeff(i,i);
328  i += 2;
329  }
330  }
331  }
332 
333  m_isInitialized = true;
334  m_eigenvectorsOk = false;//computeEigenvectors;
335 
336  return *this;
337 }
338 
339 } // end namespace Eigen
340 
341 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H
IntermediateState sqrt(const Expression &arg)
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealQZ.h:166
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
NumTraits< Scalar >::Real RealScalar
const MatrixType & matrixT() const
Returns matrix S in the QZ decomposition.
Definition: RealQZ.h:148
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
ComplexVectorType alphas() const
RealQZ & setMaxIterations(Index maxIters)
Definition: RealQZ.h:183
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
GeneralizedEigenSolver & compute(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Computes generalized eigendecomposition of given matrix.
GeneralizedEigenSolver & setMaxIterations(Index maxIters)
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > abs() const
Generic expression where a coefficient-wise binary operator is applied to two expressions.
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
EIGEN_STRONG_INLINE Scalar & coeffRef(Index rowId, Index colId)
EIGEN_STRONG_INLINE void resize(Index nbRows, Index nbCols)
CwiseBinaryOp< internal::scalar_quotient_op< ComplexScalar, Scalar >, ComplexVectorType, VectorType > EigenvalueType
Expression type for the eigenvalues as returned by eigenvalues().
GeneralizedEigenSolver(Index size)
Default constructor with memory preallocation.
Provides a generic way to set and pass user-specified options.
Definition: options.hpp:65
const MatrixType & matrixS() const
Returns matrix S in the QZ decomposition.
Definition: RealQZ.h:139
const MatrixType & matrixZ() const
Returns matrix Z in the QZ decomposition.
Definition: RealQZ.h:129
EigenvalueType eigenvalues() const
Returns an expression of the computed generalized eigenvalues.
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorsType
Type for matrix of eigenvectors as returned by eigenvectors().
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > ComplexVectorType
Type for vector of complex scalar values eigenvalues as returned by betas().
RealQZ & compute(const MatrixType &A, const MatrixType &B, bool computeQZ=true)
Computes QZ decomposition of given matrix.
Definition: RealQZ.h:557
GeneralizedEigenSolver(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Constructor; computes the generalized eigendecomposition of given matrix pair.
Computes the generalized eigenvalues and eigenvectors of a pair of general matrices.
#define eigen_assert(x)
Matrix< Scalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > VectorType
Type for vector of real scalar values eigenvalues as returned by betas().
GeneralizedEigenSolver()
Default constructor.
Matrix< Scalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > ColumnVectorType


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