123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226 |
- // This file is part of Eigen, a lightweight C++ template library
- // for linear algebra.
- //
- // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
- // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
- //
- // 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"
- 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<typename _MatrixType>
- class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType>
- {
- typedef SelfAdjointEigenSolver<_MatrixType> 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
- *
- * Accoring 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<typename MatrixType>
- GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
- 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<MatrixType> 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<Lower>();
- cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
- cholB.matrixU().template solveInPlace<OnTheRight>(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<Lower>();
- 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<Lower>();
- 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
|