// This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2008-2010 Gael Guennebaud // Copyright (C) 2010 Jitse Niesen // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. #ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H #define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H #include "./Tridiagonalization.h" #include "./InternalHeaderCheck.h" namespace Eigen { /** \eigenvalues_module \ingroup Eigenvalues_Module * * * \class GeneralizedSelfAdjointEigenSolver * * \brief Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem * * \tparam MatrixType_ the type of the matrix of which we are computing the * eigendecomposition; this is expected to be an instantiation of the Matrix * class template. * * This class solves the generalized eigenvalue problem * \f$Av = \lambda Bv \f$. In this case, the matrix \f$A \f$ should be * selfadjoint and the matrix \f$B \f$ should be positive definite. * * Only the \b lower \b triangular \b part of the input matrix is referenced. * * Call the function compute() to compute the eigenvalues and eigenvectors of * a given matrix. Alternatively, you can use the * GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) * constructor which computes the eigenvalues and eigenvectors at construction time. * Once the eigenvalue and eigenvectors are computed, they can be retrieved with the eigenvalues() * and eigenvectors() functions. * * The documentation for GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) * contains an example of the typical use of this class. * * \sa class SelfAdjointEigenSolver, class EigenSolver, class ComplexEigenSolver */ template class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver { typedef SelfAdjointEigenSolver Base; public: typedef MatrixType_ MatrixType; /** \brief Default constructor for fixed-size matrices. * * The default constructor is useful in cases in which the user intends to * perform decompositions via compute(). This constructor * can only be used if \p MatrixType_ is a fixed-size matrix; use * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices. */ GeneralizedSelfAdjointEigenSolver() : Base() {} /** \brief Constructor, pre-allocates memory for dynamic-size matrices. * * \param [in] size Positive integer, size of the matrix whose * eigenvalues and eigenvectors will be computed. * * This constructor is useful for dynamic-size matrices, when the user * intends to perform decompositions via compute(). The \p size * parameter is only used as a hint. It is not an error to give a wrong * \p size, but it may impair performance. * * \sa compute() for an example */ explicit GeneralizedSelfAdjointEigenSolver(Index size) : Base(size) {} /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil. * * \param[in] matA Selfadjoint matrix in matrix pencil. * Only the lower triangular part of the matrix is referenced. * \param[in] matB Positive-definite matrix in matrix pencil. * Only the lower triangular part of the matrix is referenced. * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}. * Default is #ComputeEigenvectors|#Ax_lBx. * * This constructor calls compute(const MatrixType&, const MatrixType&, int) * to compute the eigenvalues and (if requested) the eigenvectors of the * generalized eigenproblem \f$Ax = \lambda B x \f$ with \a matA the * selfadjoint matrix \f$A \f$ and \a matB the positive definite matrix * \f$B \f$. Each eigenvector \f$x \f$ satisfies the property * \f$x^* B x = 1 \f$. The eigenvectors are computed if * \a options contains ComputeEigenvectors. * * In addition, the two following variants can be solved via \p options: * - \c ABx_lx: \f$ABx = \lambda x \f$ * - \c BAx_lx: \f$BAx = \lambda x \f$ * * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out * * \sa compute(const MatrixType&, const MatrixType&, int) */ GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, int options = ComputeEigenvectors|Ax_lBx) : Base(matA.cols()) { compute(matA, matB, options); } /** \brief Computes generalized eigendecomposition of given matrix pencil. * * \param[in] matA Selfadjoint matrix in matrix pencil. * Only the lower triangular part of the matrix is referenced. * \param[in] matB Positive-definite matrix in matrix pencil. * Only the lower triangular part of the matrix is referenced. * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}. * Default is #ComputeEigenvectors|#Ax_lBx. * * \returns Reference to \c *this * * According to \p options, this function computes eigenvalues and (if requested) * the eigenvectors of one of the following three generalized eigenproblems: * - \c Ax_lBx: \f$Ax = \lambda B x \f$ * - \c ABx_lx: \f$ABx = \lambda x \f$ * - \c BAx_lx: \f$BAx = \lambda x \f$ * with \a matA the selfadjoint matrix \f$A \f$ and \a matB the positive definite * matrix \f$B \f$. * In addition, each eigenvector \f$x \f$ satisfies the property \f$x^* B x = 1 \f$. * * The eigenvalues() function can be used to retrieve * the eigenvalues. If \p options contains ComputeEigenvectors, then the * eigenvectors are also computed and can be retrieved by calling * eigenvectors(). * * The implementation uses LLT to compute the Cholesky decomposition * \f$B = LL^* \f$ and computes the classical eigendecomposition * of the selfadjoint matrix \f$L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx * and of \f$L^{*} A L \f$ otherwise. This solves the * generalized eigenproblem, because any solution of the generalized * eigenproblem \f$Ax = \lambda B x \f$ corresponds to a solution * \f$L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the * eigenproblem for \f$L^{-1} A (L^*)^{-1} \f$. Similar statements * can be made for the two other variants. * * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out * * \sa GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) */ GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB, int options = ComputeEigenvectors|Ax_lBx); protected: }; template GeneralizedSelfAdjointEigenSolver& GeneralizedSelfAdjointEigenSolver:: compute(const MatrixType& matA, const MatrixType& matB, int options) { eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); eigen_assert((options&~(EigVecMask|GenEigMask))==0 && (options&EigVecMask)!=EigVecMask && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx) && "invalid option parameter"); bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors); // Compute the cholesky decomposition of matB = L L' = U'U LLT cholB(matB); int type = (options&GenEigMask); if(type==0) type = Ax_lBx; if(type==Ax_lBx) { // compute C = inv(L) A inv(L') MatrixType matC = matA.template selfadjointView(); cholB.matrixL().template solveInPlace(matC); cholB.matrixU().template solveInPlace(matC); Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly ); // transform back the eigen vectors: evecs = inv(U) * evecs if(computeEigVecs) cholB.matrixU().solveInPlace(Base::m_eivec); } else if(type==ABx_lx) { // compute C = L' A L MatrixType matC = matA.template selfadjointView(); matC = matC * cholB.matrixL(); matC = cholB.matrixU() * matC; Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); // transform back the eigen vectors: evecs = inv(U) * evecs if(computeEigVecs) cholB.matrixU().solveInPlace(Base::m_eivec); } else if(type==BAx_lx) { // compute C = L' A L MatrixType matC = matA.template selfadjointView(); matC = matC * cholB.matrixL(); matC = cholB.matrixU() * matC; Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); // transform back the eigen vectors: evecs = L * evecs if(computeEigVecs) Base::m_eivec = cholB.matrixL() * Base::m_eivec; } return *this; } } // end namespace Eigen #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H