Back to home page

EIC code displayed by LXR

 
 

    


Warning, file /include/eigen3/Eigen/src/Eigenvalues/GeneralizedSelfAdjointEigenSolver.h was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010 
0011 #ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
0012 #define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
0013 
0014 #include "./Tridiagonalization.h"
0015 
0016 namespace Eigen { 
0017 
0018 /** \eigenvalues_module \ingroup Eigenvalues_Module
0019   *
0020   *
0021   * \class GeneralizedSelfAdjointEigenSolver
0022   *
0023   * \brief Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem
0024   *
0025   * \tparam _MatrixType the type of the matrix of which we are computing the
0026   * eigendecomposition; this is expected to be an instantiation of the Matrix
0027   * class template.
0028   *
0029   * This class solves the generalized eigenvalue problem
0030   * \f$ Av = \lambda Bv \f$. In this case, the matrix \f$ A \f$ should be
0031   * selfadjoint and the matrix \f$ B \f$ should be positive definite.
0032   *
0033   * Only the \b lower \b triangular \b part of the input matrix is referenced.
0034   *
0035   * Call the function compute() to compute the eigenvalues and eigenvectors of
0036   * a given matrix. Alternatively, you can use the
0037   * GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
0038   * constructor which computes the eigenvalues and eigenvectors at construction time.
0039   * Once the eigenvalue and eigenvectors are computed, they can be retrieved with the eigenvalues()
0040   * and eigenvectors() functions.
0041   *
0042   * The documentation for GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
0043   * contains an example of the typical use of this class.
0044   *
0045   * \sa class SelfAdjointEigenSolver, class EigenSolver, class ComplexEigenSolver
0046   */
0047 template<typename _MatrixType>
0048 class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType>
0049 {
0050     typedef SelfAdjointEigenSolver<_MatrixType> Base;
0051   public:
0052 
0053     typedef _MatrixType MatrixType;
0054 
0055     /** \brief Default constructor for fixed-size matrices.
0056       *
0057       * The default constructor is useful in cases in which the user intends to
0058       * perform decompositions via compute(). This constructor
0059       * can only be used if \p _MatrixType is a fixed-size matrix; use
0060       * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices.
0061       */
0062     GeneralizedSelfAdjointEigenSolver() : Base() {}
0063 
0064     /** \brief Constructor, pre-allocates memory for dynamic-size matrices.
0065       *
0066       * \param [in]  size  Positive integer, size of the matrix whose
0067       * eigenvalues and eigenvectors will be computed.
0068       *
0069       * This constructor is useful for dynamic-size matrices, when the user
0070       * intends to perform decompositions via compute(). The \p size
0071       * parameter is only used as a hint. It is not an error to give a wrong
0072       * \p size, but it may impair performance.
0073       *
0074       * \sa compute() for an example
0075       */
0076     explicit GeneralizedSelfAdjointEigenSolver(Index size)
0077         : Base(size)
0078     {}
0079 
0080     /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil.
0081       *
0082       * \param[in]  matA  Selfadjoint matrix in matrix pencil.
0083       *                   Only the lower triangular part of the matrix is referenced.
0084       * \param[in]  matB  Positive-definite matrix in matrix pencil.
0085       *                   Only the lower triangular part of the matrix is referenced.
0086       * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
0087       *                     Default is #ComputeEigenvectors|#Ax_lBx.
0088       *
0089       * This constructor calls compute(const MatrixType&, const MatrixType&, int)
0090       * to compute the eigenvalues and (if requested) the eigenvectors of the
0091       * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the
0092       * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix
0093       * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property
0094       * \f$ x^* B x = 1 \f$. The eigenvectors are computed if
0095       * \a options contains ComputeEigenvectors.
0096       *
0097       * In addition, the two following variants can be solved via \p options:
0098       * - \c ABx_lx: \f$ ABx = \lambda x \f$
0099       * - \c BAx_lx: \f$ BAx = \lambda x \f$
0100       *
0101       * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp
0102       * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out
0103       *
0104       * \sa compute(const MatrixType&, const MatrixType&, int)
0105       */
0106     GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
0107                                       int options = ComputeEigenvectors|Ax_lBx)
0108       : Base(matA.cols())
0109     {
0110       compute(matA, matB, options);
0111     }
0112 
0113     /** \brief Computes generalized eigendecomposition of given matrix pencil.
0114       *
0115       * \param[in]  matA  Selfadjoint matrix in matrix pencil.
0116       *                   Only the lower triangular part of the matrix is referenced.
0117       * \param[in]  matB  Positive-definite matrix in matrix pencil.
0118       *                   Only the lower triangular part of the matrix is referenced.
0119       * \param[in]  options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}.
0120       *                     Default is #ComputeEigenvectors|#Ax_lBx.
0121       *
0122       * \returns    Reference to \c *this
0123       *
0124       * According to \p options, this function computes eigenvalues and (if requested)
0125       * the eigenvectors of one of the following three generalized eigenproblems:
0126       * - \c Ax_lBx: \f$ Ax = \lambda B x \f$
0127       * - \c ABx_lx: \f$ ABx = \lambda x \f$
0128       * - \c BAx_lx: \f$ BAx = \lambda x \f$
0129       * with \a matA the selfadjoint matrix \f$ A \f$ and \a matB the positive definite
0130       * matrix \f$ B \f$.
0131       * In addition, each eigenvector \f$ x \f$ satisfies the property \f$ x^* B x = 1 \f$.
0132       *
0133       * The eigenvalues() function can be used to retrieve
0134       * the eigenvalues. If \p options contains ComputeEigenvectors, then the
0135       * eigenvectors are also computed and can be retrieved by calling
0136       * eigenvectors().
0137       *
0138       * The implementation uses LLT to compute the Cholesky decomposition
0139       * \f$ B = LL^* \f$ and computes the classical eigendecomposition
0140       * of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx
0141       * and of \f$ L^{*} A L \f$ otherwise. This solves the
0142       * generalized eigenproblem, because any solution of the generalized
0143       * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution
0144       * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the
0145       * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. Similar statements
0146       * can be made for the two other variants.
0147       *
0148       * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp
0149       * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out
0150       *
0151       * \sa GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int)
0152       */
0153     GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB,
0154                                                int options = ComputeEigenvectors|Ax_lBx);
0155 
0156   protected:
0157 
0158 };
0159 
0160 
0161 template<typename MatrixType>
0162 GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
0163 compute(const MatrixType& matA, const MatrixType& matB, int options)
0164 {
0165   eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
0166   eigen_assert((options&~(EigVecMask|GenEigMask))==0
0167           && (options&EigVecMask)!=EigVecMask
0168           && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
0169            || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
0170           && "invalid option parameter");
0171 
0172   bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);
0173 
0174   // Compute the cholesky decomposition of matB = L L' = U'U
0175   LLT<MatrixType> cholB(matB);
0176 
0177   int type = (options&GenEigMask);
0178   if(type==0)
0179     type = Ax_lBx;
0180 
0181   if(type==Ax_lBx)
0182   {
0183     // compute C = inv(L) A inv(L')
0184     MatrixType matC = matA.template selfadjointView<Lower>();
0185     cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
0186     cholB.matrixU().template solveInPlace<OnTheRight>(matC);
0187 
0188     Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );
0189 
0190     // transform back the eigen vectors: evecs = inv(U) * evecs
0191     if(computeEigVecs)
0192       cholB.matrixU().solveInPlace(Base::m_eivec);
0193   }
0194   else if(type==ABx_lx)
0195   {
0196     // compute C = L' A L
0197     MatrixType matC = matA.template selfadjointView<Lower>();
0198     matC = matC * cholB.matrixL();
0199     matC = cholB.matrixU() * matC;
0200 
0201     Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
0202 
0203     // transform back the eigen vectors: evecs = inv(U) * evecs
0204     if(computeEigVecs)
0205       cholB.matrixU().solveInPlace(Base::m_eivec);
0206   }
0207   else if(type==BAx_lx)
0208   {
0209     // compute C = L' A L
0210     MatrixType matC = matA.template selfadjointView<Lower>();
0211     matC = matC * cholB.matrixL();
0212     matC = cholB.matrixU() * matC;
0213 
0214     Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
0215 
0216     // transform back the eigen vectors: evecs = L * evecs
0217     if(computeEigVecs)
0218       Base::m_eivec = cholB.matrixL() * Base::m_eivec;
0219   }
0220 
0221   return *this;
0222 }
0223 
0224 } // end namespace Eigen
0225 
0226 #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H