12#ifndef EIGEN_GENERALIZEDEIGENSOLVER_H
13#define EIGEN_GENERALIZEDEIGENSOLVER_H
17#include "./InternalHeaderCheck.h"
68 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
69 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
70 Options = MatrixType::Options,
71 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
72 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
76 typedef typename MatrixType::Scalar
Scalar;
125 m_vectorsOkay(false),
136 : m_eivec(size, size),
140 m_vectorsOkay(false),
158 : m_eivec(A.rows(), A.cols()),
162 m_vectorsOkay(false),
166 compute(A, B, computeEigenvectors);
182 eigen_assert(m_vectorsOkay &&
"Eigenvectors for GeneralizedEigenSolver were not calculated.");
206 eigen_assert(m_valuesOkay &&
"GeneralizedEigenSolver is not initialized.");
217 eigen_assert(m_valuesOkay &&
"GeneralizedEigenSolver is not initialized.");
228 eigen_assert(m_valuesOkay &&
"GeneralizedEigenSolver is not initialized.");
259 eigen_assert(m_valuesOkay &&
"EigenSolver is not initialized.");
260 return m_realQZ.info();
267 m_realQZ.setMaxIterations(maxIters);
273 EIGEN_STATIC_ASSERT_NON_INTEGER(
Scalar)
279 bool m_valuesOkay, m_vectorsOkay;
284template<
typename MatrixType>
285GeneralizedEigenSolver<MatrixType>&
290 eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows());
291 Index size = A.cols();
292 m_valuesOkay =
false;
293 m_vectorsOkay =
false;
296 m_realQZ.compute(A, B, computeEigenvectors);
297 if (m_realQZ.info() ==
Success)
300 m_alphas.resize(size);
301 m_betas.resize(size);
302 if (computeEigenvectors)
304 m_eivec.resize(size,size);
317 if (i == size - 1 || mS.coeff(i+1, i) ==
Scalar(0))
320 m_alphas.
coeffRef(i) = mS.diagonal().coeff(i);
321 m_betas.coeffRef(i) = mT.diagonal().coeff(i);
322 if (computeEigenvectors)
324 v.setConstant(
Scalar(0.0));
325 v.coeffRef(i) =
Scalar(1.0);
327 if(
abs(m_betas.coeffRef(i)) >= (std::numeric_limits<RealScalar>::min)())
330 const Scalar alpha =
real(m_alphas.coeffRef(i));
331 const Scalar beta = m_betas.coeffRef(i);
332 for (
Index j = i-1; j >= 0; j--)
334 const Index st = j+1;
335 const Index sz = i-j;
336 if (j > 0 && mS.coeff(j, j-1) !=
Scalar(0))
339 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) );
340 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);
341 v.template segment<2>(j-1) = lhs.
partialPivLu().solve(rhs);
346 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 m_eivec.col(i).real().noalias() = m_realQZ.matrixZ().transpose() * v;
351 m_eivec.col(i).real().normalize();
352 m_eivec.col(i).imag().setConstant(0);
363 RealScalar a = mT.diagonal().coeff(i),
364 b = mT.diagonal().coeff(i+1);
365 const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b;
373 m_alphas.coeffRef(i) =
conj(alpha);
374 m_alphas.coeffRef(i+1) = alpha;
376 if (computeEigenvectors) {
381 cv.
coeffRef(i) = -(
static_cast<Scalar>(beta*mS.coeffRef(i,i+1)) - alpha*mT.coeffRef(i,i+1))
382 / (
static_cast<Scalar>(beta*mS.coeffRef(i,i)) - alpha*mT.coeffRef(i,i));
383 for (
Index j = i-1; j >= 0; j--)
385 const Index st = j+1;
386 const Index sz = i+1-j;
387 if (j > 0 && mS.coeff(j, j-1) !=
Scalar(0))
390 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) );
392 cv.template segment<2>(j-1) = lhs.
partialPivLu().solve(rhs);
395 cv.
coeffRef(j) = cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum()
396 / (alpha*mT.coeffRef(j,j) -
static_cast<Scalar>(beta*mS.coeffRef(j,j)));
399 m_eivec.col(i+1).noalias() = (m_realQZ.matrixZ().transpose() * cv);
400 m_eivec.col(i+1).normalize();
401 m_eivec.col(i) = m_eivec.col(i+1).conjugate();
408 m_vectorsOkay = computeEigenvectors;
Generic expression where a coefficient-wise binary operator is applied to two expressions.
Definition: CwiseBinaryOp.h:86
Computes the generalized eigenvalues and eigenvectors of a pair of general matrices.
Definition: GeneralizedEigenSolver.h:61
GeneralizedEigenSolver & compute(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Computes generalized eigendecomposition of given matrix.
Definition: GeneralizedEigenSolver.h:286
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Definition: GeneralizedEigenSolver.h:86
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > ComplexVectorType
Type for vector of complex scalar values eigenvalues as returned by alphas().
Definition: GeneralizedEigenSolver.h:100
Matrix< Scalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > VectorType
Type for vector of real scalar values eigenvalues as returned by betas().
Definition: GeneralizedEigenSolver.h:93
CwiseBinaryOp< internal::scalar_quotient_op< ComplexScalar, Scalar >, ComplexVectorType, VectorType > EigenvalueType
Expression type for the eigenvalues as returned by eigenvalues().
Definition: GeneralizedEigenSolver.h:104
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: GeneralizedEigenSolver.h:76
ComplexVectorType alphas() const
Definition: GeneralizedEigenSolver.h:215
GeneralizedEigenSolver()
Default constructor.
Definition: GeneralizedEigenSolver.h:120
VectorType betas() const
Definition: GeneralizedEigenSolver.h:226
Eigen::Index Index
Definition: GeneralizedEigenSolver.h:78
GeneralizedEigenSolver(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Constructor; computes the generalized eigendecomposition of given matrix pair.
Definition: GeneralizedEigenSolver.h:157
EigenvalueType eigenvalues() const
Returns an expression of the computed generalized eigenvalues.
Definition: GeneralizedEigenSolver.h:204
GeneralizedEigenSolver(Index size)
Default constructor with memory preallocation.
Definition: GeneralizedEigenSolver.h:135
GeneralizedEigenSolver & setMaxIterations(Index maxIters)
Definition: GeneralizedEigenSolver.h:265
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorsType
Type for matrix of eigenvectors as returned by eigenvectors().
Definition: GeneralizedEigenSolver.h:111
MatrixType_ MatrixType
Synonym for the template parameter MatrixType_.
Definition: GeneralizedEigenSolver.h:65
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:98
const DiagonalWrapper< const Derived > asDiagonal() const
Definition: DiagonalMatrix.h:325
const PartialPivLU< PlainObject > partialPivLu() const
Definition: PartialPivLU.h:599
Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:187
const Scalar & coeff(Index rowId, Index colId) const
Definition: PlainObjectBase.h:164
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:564
Performs a real QZ decomposition of a pair of square matrices.
Definition: RealQZ.h:60
ComputationInfo
Definition: Constants.h:442
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: B01_Experimental.dox:1
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(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_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:235