Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.99 (git rev 10c77b0ff44d0b9cb0b252cfa0ccaaa39d3c5da4)
GeneralizedEigenSolver.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012-2016 Gael Guennebaud <gael.guennebaud@inria.fr>
5 // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
6 // Copyright (C) 2016 Tobias Wood <tobias@spinicist.org.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_GENERALIZEDEIGENSOLVER_H
13 #define EIGEN_GENERALIZEDEIGENSOLVER_H
14 
15 #include "./RealQZ.h"
16 
17 namespace Eigen {
18 
58 template<typename _MatrixType> class GeneralizedEigenSolver
59 {
60  public:
61 
63  typedef _MatrixType MatrixType;
64 
65  enum {
66  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
67  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
68  Options = MatrixType::Options,
69  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
70  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
71  };
72 
74  typedef typename MatrixType::Scalar Scalar;
75  typedef typename NumTraits<Scalar>::Real RealScalar;
76  typedef Eigen::Index Index;
77 
84  typedef std::complex<RealScalar> ComplexScalar;
85 
92 
99 
103 
110 
119  : m_eivec(),
120  m_alphas(),
121  m_betas(),
122  m_valuesOkay(false),
123  m_vectorsOkay(false),
124  m_realQZ()
125  {}
126 
134  : m_eivec(size, size),
135  m_alphas(size),
136  m_betas(size),
137  m_valuesOkay(false),
138  m_vectorsOkay(false),
139  m_realQZ(size),
140  m_tmp(size)
141  {}
142 
155  GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true)
156  : m_eivec(A.rows(), A.cols()),
157  m_alphas(A.cols()),
158  m_betas(A.cols()),
159  m_valuesOkay(false),
160  m_vectorsOkay(false),
161  m_realQZ(A.cols()),
162  m_tmp(A.cols())
163  {
164  compute(A, B, computeEigenvectors);
165  }
166 
167  /* \brief Returns the computed generalized eigenvectors.
168  *
169  * \returns %Matrix whose columns are the (possibly complex) right eigenvectors.
170  * i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues.
171  *
172  * \pre Either the constructor
173  * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function
174  * compute(const MatrixType&, const MatrixType& bool) has been called before, and
175  * \p computeEigenvectors was set to true (the default).
176  *
177  * \sa eigenvalues()
178  */
179  EigenvectorsType eigenvectors() const {
180  eigen_assert(m_vectorsOkay && "Eigenvectors for GeneralizedEigenSolver were not calculated.");
181  return m_eivec;
182  }
183 
203  {
204  eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
205  return EigenvalueType(m_alphas,m_betas);
206  }
207 
214  {
215  eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
216  return m_alphas;
217  }
218 
225  {
226  eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized.");
227  return m_betas;
228  }
229 
253  GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true);
254 
255  ComputationInfo info() const
256  {
257  eigen_assert(m_valuesOkay && "EigenSolver is not initialized.");
258  return m_realQZ.info();
259  }
260 
264  {
265  m_realQZ.setMaxIterations(maxIters);
266  return *this;
267  }
268 
269  protected:
270 
271  static void check_template_parameters()
272  {
273  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
274  EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
275  }
276 
277  EigenvectorsType m_eivec;
278  ComplexVectorType m_alphas;
279  VectorType m_betas;
280  bool m_valuesOkay, m_vectorsOkay;
281  RealQZ<MatrixType> m_realQZ;
282  ComplexVectorType m_tmp;
283 };
284 
285 template<typename MatrixType>
286 GeneralizedEigenSolver<MatrixType>&
287 GeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors)
288 {
289  check_template_parameters();
290 
291  using std::sqrt;
292  using std::abs;
293  eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
294  Index size = A.cols();
295  m_valuesOkay = false;
296  m_vectorsOkay = false;
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  if (m_realQZ.info() == Success)
301  {
302  // Resize storage
303  m_alphas.resize(size);
304  m_betas.resize(size);
305  if (computeEigenvectors)
306  {
307  m_eivec.resize(size,size);
308  m_tmp.resize(size);
309  }
310 
311  // Aliases:
312  Map<VectorType> v(reinterpret_cast<Scalar*>(m_tmp.data()), size);
313  ComplexVectorType &cv = m_tmp;
314  const MatrixType &mS = m_realQZ.matrixS();
315  const MatrixType &mT = m_realQZ.matrixT();
316 
317  Index i = 0;
318  while (i < size)
319  {
320  if (i == size - 1 || mS.coeff(i+1, i) == Scalar(0))
321  {
322  // Real eigenvalue
323  m_alphas.coeffRef(i) = mS.diagonal().coeff(i);
324  m_betas.coeffRef(i) = mT.diagonal().coeff(i);
325  if (computeEigenvectors)
326  {
327  v.setConstant(Scalar(0.0));
328  v.coeffRef(i) = Scalar(1.0);
329  // For singular eigenvalues do nothing more
330  if(abs(m_betas.coeffRef(i)) >= (std::numeric_limits<RealScalar>::min)())
331  {
332  // Non-singular eigenvalue
333  const Scalar alpha = real(m_alphas.coeffRef(i));
334  const Scalar beta = m_betas.coeffRef(i);
335  for (Index j = i-1; j >= 0; j--)
336  {
337  const Index st = j+1;
338  const Index sz = i-j;
339  if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
340  {
341  // 2x2 block
342  Matrix<Scalar, 2, 1> rhs = (alpha*mT.template block<2,Dynamic>(j-1,st,2,sz) - beta*mS.template block<2,Dynamic>(j-1,st,2,sz)) .lazyProduct( v.segment(st,sz) );
343  Matrix<Scalar, 2, 2> lhs = beta * mS.template block<2,2>(j-1,j-1) - alpha * mT.template block<2,2>(j-1,j-1);
344  v.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
345  j--;
346  }
347  else
348  {
349  v.coeffRef(j) = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() / (beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j));
350  }
351  }
352  }
353  m_eivec.col(i).real().noalias() = m_realQZ.matrixZ().transpose() * v;
354  m_eivec.col(i).real().normalize();
355  m_eivec.col(i).imag().setConstant(0);
356  }
357  ++i;
358  }
359  else
360  {
361  // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T
362  // 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):
363 
364  // T = [a 0]
365  // [0 b]
366  RealScalar a = mT.diagonal().coeff(i),
367  b = mT.diagonal().coeff(i+1);
368  const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b;
369 
370  // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug.
371  Matrix<RealScalar,2,2> S2 = mS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal();
372 
373  Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1));
374  Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1)));
375  const ComplexScalar alpha = ComplexScalar(S2.coeff(1,1) + p, (beta > 0) ? z : -z);
376  m_alphas.coeffRef(i) = conj(alpha);
377  m_alphas.coeffRef(i+1) = alpha;
378 
379  if (computeEigenvectors) {
380  // Compute eigenvector in position (i+1) and then position (i) is just the conjugate
381  cv.setZero();
382  cv.coeffRef(i+1) = Scalar(1.0);
383  // here, the "static_cast" workaound expression template issues.
384  cv.coeffRef(i) = -(static_cast<Scalar>(beta*mS.coeffRef(i,i+1)) - alpha*mT.coeffRef(i,i+1))
385  / (static_cast<Scalar>(beta*mS.coeffRef(i,i)) - alpha*mT.coeffRef(i,i));
386  for (Index j = i-1; j >= 0; j--)
387  {
388  const Index st = j+1;
389  const Index sz = i+1-j;
390  if (j > 0 && mS.coeff(j, j-1) != Scalar(0))
391  {
392  // 2x2 block
393  Matrix<ComplexScalar, 2, 1> rhs = (alpha*mT.template block<2,Dynamic>(j-1,st,2,sz) - beta*mS.template block<2,Dynamic>(j-1,st,2,sz)) .lazyProduct( cv.segment(st,sz) );
394  Matrix<ComplexScalar, 2, 2> lhs = beta * mS.template block<2,2>(j-1,j-1) - alpha * mT.template block<2,2>(j-1,j-1);
395  cv.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs);
396  j--;
397  } else {
398  cv.coeffRef(j) = cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum()
399  / (alpha*mT.coeffRef(j,j) - static_cast<Scalar>(beta*mS.coeffRef(j,j)));
400  }
401  }
402  m_eivec.col(i+1).noalias() = (m_realQZ.matrixZ().transpose() * cv);
403  m_eivec.col(i+1).normalize();
404  m_eivec.col(i) = m_eivec.col(i+1).conjugate();
405  }
406  i += 2;
407  }
408  }
409 
410  m_valuesOkay = true;
411  m_vectorsOkay = computeEigenvectors;
412  }
413  return *this;
414 }
415 
416 } // end namespace Eigen
417 
418 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H
Generic expression where a coefficient-wise binary operator is applied to two expressions.
Definition: CwiseBinaryOp.h:84
Computes the generalized eigenvalues and eigenvectors of a pair of general matrices.
Definition: GeneralizedEigenSolver.h:59
GeneralizedEigenSolver & compute(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Computes generalized eigendecomposition of given matrix.
Definition: GeneralizedEigenSolver.h:287
GeneralizedEigenSolver(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Constructor; computes the generalized eigendecomposition of given matrix pair.
Definition: GeneralizedEigenSolver.h:155
GeneralizedEigenSolver & setMaxIterations(Index maxIters)
Definition: GeneralizedEigenSolver.h:263
Eigen::Index Index
Definition: GeneralizedEigenSolver.h:76
_MatrixType MatrixType
Synonym for the template parameter _MatrixType.
Definition: GeneralizedEigenSolver.h:63
Matrix< Scalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > VectorType
Type for vector of real scalar values eigenvalues as returned by betas().
Definition: GeneralizedEigenSolver.h:91
EigenvalueType eigenvalues() const
Returns an expression of the computed generalized eigenvalues.
Definition: GeneralizedEigenSolver.h:202
ComplexVectorType alphas() const
Definition: GeneralizedEigenSolver.h:213
GeneralizedEigenSolver(Index size)
Default constructor with memory preallocation.
Definition: GeneralizedEigenSolver.h:133
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Definition: GeneralizedEigenSolver.h:84
VectorType betas() const
Definition: GeneralizedEigenSolver.h:224
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > ComplexVectorType
Type for vector of complex scalar values eigenvalues as returned by alphas().
Definition: GeneralizedEigenSolver.h:98
CwiseBinaryOp< internal::scalar_quotient_op< ComplexScalar, Scalar >, ComplexVectorType, VectorType > EigenvalueType
Expression type for the eigenvalues as returned by eigenvalues().
Definition: GeneralizedEigenSolver.h:102
GeneralizedEigenSolver()
Default constructor.
Definition: GeneralizedEigenSolver.h:118
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: GeneralizedEigenSolver.h:74
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorsType
Type for matrix of eigenvectors as returned by eigenvectors().
Definition: GeneralizedEigenSolver.h:109
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:96
const DiagonalWrapper< const Derived > asDiagonal() const
Definition: DiagonalMatrix.h:325
const PartialPivLU< PlainObject > partialPivLu() const
Definition: PartialPivLU.h:597
const Scalar & coeff(Index rowId, Index colId) const
Definition: PlainObjectBase.h:160
Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:183
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:562
RealQZ & setMaxIterations(Index maxIters)
Definition: RealQZ.h:186
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealQZ.h:169
ComputationInfo
Definition: Constants.h:440
@ Success
Definition: Constants.h:442
Namespace containing all symbols from the Eigen library.
Definition: Core:134
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_real_op< typename Derived::Scalar >, const Derived > real(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:213