Eigen  3.2.93
GeneralizedEigenSolver.h
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;
74  typedef typename NumTraits<Scalar>::Real RealScalar;
75  typedef Eigen::Index Index;
76 
83  typedef std::complex<RealScalar> ComplexScalar;
84 
91 
98 
102 
109 
117  GeneralizedEigenSolver() : m_eivec(), m_alphas(), m_betas(), m_isInitialized(false), m_realQZ(), m_matS(), m_tmp() {}
118 
125  explicit GeneralizedEigenSolver(Index size)
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.");
201  return EigenvalueType(m_alphas,m_betas);
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 
251  ComputationInfo info() const
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 
267  static void check_template_parameters()
268  {
269  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
270  EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
271  }
272 
273  MatrixType m_eivec;
274  ComplexVectorType m_alphas;
275  VectorType m_betas;
276  bool m_isInitialized;
277  bool m_eigenvectorsOk;
278  RealQZ<MatrixType> m_realQZ;
279  MatrixType m_matS;
280 
282  ColumnVectorType m_tmp;
283 };
284 
285 //template<typename MatrixType>
286 //typename GeneralizedEigenSolver<MatrixType>::EigenvectorsType GeneralizedEigenSolver<MatrixType>::eigenvectors() const
287 //{
288 // eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
289 // eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
290 // Index n = m_eivec.cols();
291 // EigenvectorsType matV(n,n);
292 // // TODO
293 // return matV;
294 //}
295 
296 template<typename MatrixType>
298 GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
299 {
300  check_template_parameters();
301 
302  using std::sqrt;
303  using std::abs;
304  eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
305 
306  // Reduce to generalized real Schur form:
307  // A = Q S Z and B = Q T Z
308  m_realQZ.compute(A, B, computeEigenvectors);
309 
310  if (m_realQZ.info() == Success)
311  {
312  m_matS = m_realQZ.matrixS();
313  const MatrixType &matT = m_realQZ.matrixT();
314  if (computeEigenvectors)
315  m_eivec = m_realQZ.matrixZ().transpose();
316 
317  // Compute eigenvalues from matS
318  m_alphas.resize(A.cols());
319  m_betas.resize(A.cols());
320  Index i = 0;
321  while (i < A.cols())
322  {
323  if (i == A.cols() - 1 || m_matS.coeff(i+1, i) == Scalar(0))
324  {
325  m_alphas.coeffRef(i) = m_matS.coeff(i, i);
326  m_betas.coeffRef(i) = matT.coeff(i,i);
327  ++i;
328  }
329  else
330  {
331  // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T
332  // Then taking beta=T_00*T_11, we can avoid any division, and alpha is the eigenvalues of A = (U^-1 * S * U) * diag(T_11,T_00):
333 
334  // T = [a 0]
335  // [0 b]
336  RealScalar a = matT.diagonal().coeff(i),
337  b = matT.diagonal().coeff(i+1);
338  // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug.
339  Matrix<RealScalar,2,2> S2 = m_matS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal();
340 
341  Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1));
342  Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1)));
343  m_alphas.coeffRef(i) = ComplexScalar(S2.coeff(1,1) + p, z);
344  m_alphas.coeffRef(i+1) = ComplexScalar(S2.coeff(1,1) + p, -z);
345 
346  m_betas.coeffRef(i) =
347  m_betas.coeffRef(i+1) = a*b;
348 
349  i += 2;
350  }
351  }
352  }
353 
354  m_isInitialized = true;
355  m_eigenvectorsOk = false;//computeEigenvectors;
356 
357  return *this;
358 }
359 
360 } // end namespace Eigen
361 
362 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H
const MatrixType & matrixS() const
Returns matrix S in the QZ decomposition.
Definition: RealQZ.h:139
EigenvalueType eigenvalues() const
Returns an expression of the computed generalized eigenvalues.
Definition: GeneralizedEigenSolver.h:198
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Namespace containing all symbols from the Eigen library.
Definition: Core:271
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:167
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Definition: GeneralizedEigenSolver.h:83
RealQZ & setMaxIterations(Index maxIters)
Definition: RealQZ.h:183
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:273
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: GeneralizedEigenSolver.h:73
GeneralizedEigenSolver & compute(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Computes generalized eigendecomposition of given matrix.
Definition: GeneralizedEigenSolver.h:298
GeneralizedEigenSolver & setMaxIterations(Index maxIters)
Definition: GeneralizedEigenSolver.h:259
Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:177
Generic expression where a coefficient-wise binary operator is applied to two expressions.
Definition: CwiseBinaryOp.h:77
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
Definition: GeneralizedEigenSolver.h:62
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: XprHelper.h:35
CwiseBinaryOp< internal::scalar_quotient_op< ComplexScalar, Scalar >, ComplexVectorType, VectorType > EigenvalueType
Expression type for the eigenvalues as returned by eigenvalues().
Definition: GeneralizedEigenSolver.h:101
GeneralizedEigenSolver(Index size)
Default constructor with memory preallocation.
Definition: GeneralizedEigenSolver.h:125
const DiagonalWrapper< const Derived > asDiagonal() const
Definition: DiagonalMatrix.h:277
const MatrixType & matrixZ() const
Returns matrix Z in the QZ decomposition.
Definition: RealQZ.h:129
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealQZ.h:166
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorsType
Type for matrix of eigenvectors as returned by eigenvectors().
Definition: GeneralizedEigenSolver.h:108
Definition: Constants.h:432
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
const Scalar & coeff(Index rowId, Index colId) const
Definition: PlainObjectBase.h:154
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > ComplexVectorType
Type for vector of complex scalar values eigenvalues as returned by betas().
Definition: GeneralizedEigenSolver.h:97
const MatrixType & matrixT() const
Returns matrix S in the QZ decomposition.
Definition: RealQZ.h:148
RealQZ & compute(const MatrixType &A, const MatrixType &B, bool computeQZ=true)
Computes QZ decomposition of given matrix.
Definition: RealQZ.h:556
GeneralizedEigenSolver(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Constructor; computes the generalized eigendecomposition of given matrix pair.
Definition: GeneralizedEigenSolver.h:148
Computes the generalized eigenvalues and eigenvectors of a pair of general matrices.
Definition: GeneralizedEigenSolver.h:57
Matrix< Scalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > VectorType
Type for vector of real scalar values eigenvalues as returned by betas().
Definition: GeneralizedEigenSolver.h:90
ComputationInfo
Definition: Constants.h:430
VectorType betas() const
Definition: GeneralizedEigenSolver.h:220
ComplexVectorType alphas() const
Definition: GeneralizedEigenSolver.h:209
GeneralizedEigenSolver()
Default constructor.
Definition: GeneralizedEigenSolver.h:117
Eigen::Index Index
Definition: GeneralizedEigenSolver.h:75