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);
179 EigenvectorsType eigenvectors()
const {
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);
277 EigenvectorsType m_eivec;
278 ComplexVectorType m_alphas;
280 bool m_valuesOkay, m_vectorsOkay;
282 ComplexVectorType m_tmp;
285 template<
typename 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);
305 if (computeEigenvectors)
307 m_eivec.
resize(size,size);
320 if (i == size - 1 || mS.coeff(i+1, i) ==
Scalar(0))
323 m_alphas.
coeffRef(i) = mS.diagonal().coeff(i);
324 m_betas.
coeffRef(i) = mT.diagonal().coeff(i);
325 if (computeEigenvectors)
327 v.setConstant(
Scalar(0.0));
328 v.coeffRef(i) =
Scalar(1.0);
330 if(
abs(m_betas.
coeffRef(i)) >= (std::numeric_limits<RealScalar>::min)())
335 for (
Index j = i-1; j >= 0; j--)
337 const Index st = j+1;
338 const Index sz = i-j;
339 if (j > 0 && mS.coeff(j, j-1) !=
Scalar(0))
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);
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));
353 m_eivec.col(i).real().noalias() = m_realQZ.
matrixZ().transpose() * v;
354 m_eivec.col(i).real().normalize();
366 RealScalar a = mT.diagonal().coeff(i),
367 b = mT.diagonal().coeff(i+1);
379 if (computeEigenvectors) {
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--)
388 const Index st = j+1;
389 const Index sz = i+1-j;
390 if (j > 0 && mS.coeff(j, j-1) !=
Scalar(0))
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) );
395 cv.template segment<2>(j-1) = lhs.
partialPivLu().solve(rhs);
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)));
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();
411 m_vectorsOkay = computeEigenvectors;
418 #endif // EIGEN_GENERALIZEDEIGENSOLVER_H const MatrixType & matrixS() const
Returns matrix S in the QZ decomposition.
Definition: RealQZ.h:142
EigenvalueType eigenvalues() const
Returns an expression of the computed generalized eigenvalues.
Definition: GeneralizedEigenSolver.h:202
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
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)
Namespace containing all symbols from the Eigen library.
Definition: Core:141
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
std::complex< RealScalar > ComplexScalar
Complex scalar type for MatrixType.
Definition: GeneralizedEigenSolver.h:84
RealQZ & setMaxIterations(Index maxIters)
Definition: RealQZ.h:186
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
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
GeneralizedEigenSolver & setMaxIterations(Index maxIters)
Definition: GeneralizedEigenSolver.h:263
Scalar & coeffRef(Index rowId, Index colId)
Definition: PlainObjectBase.h:175
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:562
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:74
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 DiagonalWrapper< const Derived > asDiagonal() const
Definition: DiagonalMatrix.h:325
Derived & setConstant(Index size, const Scalar &val)
Definition: CwiseNullaryOp.h:361
const MatrixType & matrixZ() const
Returns matrix Z in the QZ decomposition.
Definition: RealQZ.h:132
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealQZ.h:169
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_real_op< typename Derived::Scalar >, const Derived > real(const Eigen::ArrayBase< Derived > &x)
const PartialPivLU< PlainObject > partialPivLu() const
Definition: PartialPivLU.h:602
Matrix< ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > EigenvectorsType
Type for matrix of eigenvectors as returned by eigenvectors().
Definition: GeneralizedEigenSolver.h:109
Definition: Constants.h:442
const Scalar * data() const
Definition: PlainObjectBase.h:247
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:152
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
const MatrixType & matrixT() const
Returns matrix S in the QZ decomposition.
Definition: RealQZ.h:151
RealQZ & compute(const MatrixType &A, const MatrixType &B, bool computeQZ=true)
Computes QZ decomposition of given matrix.
Definition: RealQZ.h:559
GeneralizedEigenSolver(const MatrixType &A, const MatrixType &B, bool computeEigenvectors=true)
Constructor; computes the generalized eigendecomposition of given matrix pair.
Definition: GeneralizedEigenSolver.h:155
Computes the generalized eigenvalues and eigenvectors of a pair of general matrices.
Definition: GeneralizedEigenSolver.h:58
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:440
VectorType betas() const
Definition: GeneralizedEigenSolver.h:224
ComplexVectorType alphas() const
Definition: GeneralizedEigenSolver.h:213
GeneralizedEigenSolver()
Default constructor.
Definition: GeneralizedEigenSolver.h:118
Eigen::Index Index
Definition: GeneralizedEigenSolver.h:76