12 #ifndef EIGEN_GENERALIZEDEIGENSOLVER_H
13 #define EIGEN_GENERALIZEDEIGENSOLVER_H
66 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
67 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
68 Options = MatrixType::Options,
69 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
70 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
74 typedef typename MatrixType::Scalar
Scalar;
123 m_vectorsOkay(false),
134 : m_eivec(size, size),
138 m_vectorsOkay(false),
156 : m_eivec(A.rows(), A.cols()),
160 m_vectorsOkay(false),
164 compute(A, B, computeEigenvectors);
180 eigen_assert(m_vectorsOkay &&
"Eigenvectors for GeneralizedEigenSolver were not calculated.");
204 eigen_assert(m_valuesOkay &&
"GeneralizedEigenSolver is not initialized.");
215 eigen_assert(m_valuesOkay &&
"GeneralizedEigenSolver is not initialized.");
226 eigen_assert(m_valuesOkay &&
"GeneralizedEigenSolver is not initialized.");
257 eigen_assert(m_valuesOkay &&
"EigenSolver is not initialized.");
258 return m_realQZ.
info();
271 static void check_template_parameters()
273 EIGEN_STATIC_ASSERT_NON_INTEGER(
Scalar);
280 bool m_valuesOkay, m_vectorsOkay;
281 RealQZ<MatrixType> m_realQZ;
285 template<
typename MatrixType>
286 GeneralizedEigenSolver<MatrixType>&
289 check_template_parameters();
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;
299 m_realQZ.compute(A, B, computeEigenvectors);
300 if (m_realQZ.info() ==
Success)
303 m_alphas.resize(size);
304 m_betas.resize(size);
305 if (computeEigenvectors)
307 m_eivec.resize(size,size);
321 if (i == size - 1 || mS.coeff(i+1, i) ==
Scalar(0))
324 m_alphas.
coeffRef(i) = mS.diagonal().coeff(i);
325 m_betas.coeffRef(i) = mT.diagonal().coeff(i);
326 if (computeEigenvectors)
328 v.setConstant(
Scalar(0.0));
329 v.coeffRef(i) =
Scalar(1.0);
331 if(
abs(m_betas.coeffRef(i)) >= (std::numeric_limits<RealScalar>::min)())
334 const Scalar alpha =
real(m_alphas.coeffRef(i));
335 const Scalar beta = m_betas.coeffRef(i);
336 for (
Index j = i-1; j >= 0; j--)
338 const Index st = j+1;
339 const Index sz = i-j;
340 if (j > 0 && mS.coeff(j, j-1) !=
Scalar(0))
343 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) );
344 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);
345 v.template segment<2>(j-1) = lhs.
partialPivLu().solve(rhs);
350 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));
354 m_eivec.col(i).real().noalias() = mZ.transpose() * v;
355 m_eivec.col(i).real().normalize();
356 m_eivec.col(i).imag().setConstant(0);
367 RealScalar a = mT.diagonal().coeff(i),
368 b = mT.diagonal().coeff(i+1);
369 const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b;
377 m_alphas.coeffRef(i) =
conj(alpha);
378 m_alphas.coeffRef(i+1) = alpha;
380 if (computeEigenvectors) {
385 cv.
coeffRef(i) = -(
static_cast<Scalar>(beta*mS.coeffRef(i,i+1)) - alpha*mT.coeffRef(i,i+1))
386 / (static_cast<Scalar>(beta*mS.coeffRef(i,i)) - alpha*mT.coeffRef(i,i));
387 for (
Index j = i-1; j >= 0; j--)
389 const Index st = j+1;
390 const Index sz = i+1-j;
391 if (j > 0 && mS.coeff(j, j-1) !=
Scalar(0))
394 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) );
396 cv.template segment<2>(j-1) = lhs.
partialPivLu().solve(rhs);
399 cv.
coeffRef(j) = cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum()
400 / (alpha*mT.coeffRef(j,j) -
static_cast<Scalar>(beta*mS.coeffRef(j,j)));
403 m_eivec.col(i+1).noalias() = (mZ.transpose() * cv);
404 m_eivec.col(i+1).normalize();
405 m_eivec.col(i) = m_eivec.col(i+1).conjugate();
412 m_vectorsOkay = computeEigenvectors;
419 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:88
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_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealQZ.h:166
VectorType betas() const
Definition: GeneralizedEigenSolver.h:224
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Definition: GeneralizedEigenSolver.h:84
ComplexVectorType alphas() const
Definition: GeneralizedEigenSolver.h:213
RealQZ & setMaxIterations(Index maxIters)
Definition: RealQZ.h:183
MatrixType::Scalar Scalar
Scalar type for matrices of type MatrixType.
Definition: GeneralizedEigenSolver.h:74
GeneralizedEigenSolver & compute(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Computes generalized eigendecomposition of given matrix.
Definition: GeneralizedEigenSolver.h:287
const PartialPivLU< PlainObject > partialPivLu() const
Definition: PartialPivLU.h:589
GeneralizedEigenSolver & setMaxIterations(Index maxIters)
Definition: GeneralizedEigenSolver.h:263
Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:183
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:515
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:63
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
CwiseBinaryOp< internal::scalar_quotient_op< ComplexScalar, Scalar >, ComplexVectorType, VectorType > EigenvalueType
Expression type for the eigenvalues as returned by eigenvalues().
Definition: GeneralizedEigenSolver.h:102
GeneralizedEigenSolver(Index size)
Default constructor with memory preallocation.
Definition: GeneralizedEigenSolver.h:133
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_real_op< typename Derived::Scalar >, const Derived > real(const Eigen::ArrayBase< Derived > &x)
EigenvalueType eigenvalues() const
Returns an expression of the computed generalized eigenvalues.
Definition: GeneralizedEigenSolver.h:202
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorsType
Type for matrix of eigenvectors as returned by eigenvectors().
Definition: GeneralizedEigenSolver.h:109
Definition: Constants.h:432
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
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
GeneralizedEigenSolver(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Constructor; computes the generalized eigendecomposition of given matrix pair.
Definition: GeneralizedEigenSolver.h:155
const DiagonalWrapper< const Derived > asDiagonal() const
Definition: DiagonalMatrix.h:277
Computes the generalized eigenvalues and eigenvectors of a pair of general matrices.
Definition: GeneralizedEigenSolver.h:58
const Scalar & coeff(Index rowId, Index colId) const
Definition: PlainObjectBase.h:160
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
ComputationInfo
Definition: Constants.h:430
GeneralizedEigenSolver()
Default constructor.
Definition: GeneralizedEigenSolver.h:118
Eigen::Index Index
Definition: GeneralizedEigenSolver.h:76