1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 // Copyright (C) 2010 Jitse Niesen <jitse (at) maths.leeds.ac.uk> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H 12 #define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H 13 14 #include "./Tridiagonalization.h" 15 16 namespace Eigen { 17 18 /** \eigenvalues_module \ingroup Eigenvalues_Module 19 * 20 * 21 * \class GeneralizedSelfAdjointEigenSolver 22 * 23 * \brief Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem 24 * 25 * \tparam _MatrixType the type of the matrix of which we are computing the 26 * eigendecomposition; this is expected to be an instantiation of the Matrix 27 * class template. 28 * 29 * This class solves the generalized eigenvalue problem 30 * \f$ Av = \lambda Bv \f$. In this case, the matrix \f$ A \f$ should be 31 * selfadjoint and the matrix \f$ B \f$ should be positive definite. 32 * 33 * Only the \b lower \b triangular \b part of the input matrix is referenced. 34 * 35 * Call the function compute() to compute the eigenvalues and eigenvectors of 36 * a given matrix. Alternatively, you can use the 37 * GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) 38 * constructor which computes the eigenvalues and eigenvectors at construction time. 39 * Once the eigenvalue and eigenvectors are computed, they can be retrieved with the eigenvalues() 40 * and eigenvectors() functions. 41 * 42 * The documentation for GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) 43 * contains an example of the typical use of this class. 44 * 45 * \sa class SelfAdjointEigenSolver, class EigenSolver, class ComplexEigenSolver 46 */ 47 template<typename _MatrixType> 48 class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType> 49 { 50 typedef SelfAdjointEigenSolver<_MatrixType> Base; 51 public: 52 53 typedef _MatrixType MatrixType; 54 55 /** \brief Default constructor for fixed-size matrices. 56 * 57 * The default constructor is useful in cases in which the user intends to 58 * perform decompositions via compute(). This constructor 59 * can only be used if \p _MatrixType is a fixed-size matrix; use 60 * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices. 61 */ 62 GeneralizedSelfAdjointEigenSolver() : Base() {} 63 64 /** \brief Constructor, pre-allocates memory for dynamic-size matrices. 65 * 66 * \param [in] size Positive integer, size of the matrix whose 67 * eigenvalues and eigenvectors will be computed. 68 * 69 * This constructor is useful for dynamic-size matrices, when the user 70 * intends to perform decompositions via compute(). The \p size 71 * parameter is only used as a hint. It is not an error to give a wrong 72 * \p size, but it may impair performance. 73 * 74 * \sa compute() for an example 75 */ 76 explicit GeneralizedSelfAdjointEigenSolver(Index size) 77 : Base(size) 78 {} 79 80 /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil. 81 * 82 * \param[in] matA Selfadjoint matrix in matrix pencil. 83 * Only the lower triangular part of the matrix is referenced. 84 * \param[in] matB Positive-definite matrix in matrix pencil. 85 * Only the lower triangular part of the matrix is referenced. 86 * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}. 87 * Default is #ComputeEigenvectors|#Ax_lBx. 88 * 89 * This constructor calls compute(const MatrixType&, const MatrixType&, int) 90 * to compute the eigenvalues and (if requested) the eigenvectors of the 91 * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the 92 * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix 93 * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property 94 * \f$ x^* B x = 1 \f$. The eigenvectors are computed if 95 * \a options contains ComputeEigenvectors. 96 * 97 * In addition, the two following variants can be solved via \p options: 98 * - \c ABx_lx: \f$ ABx = \lambda x \f$ 99 * - \c BAx_lx: \f$ BAx = \lambda x \f$ 100 * 101 * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp 102 * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out 103 * 104 * \sa compute(const MatrixType&, const MatrixType&, int) 105 */ 106 GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, 107 int options = ComputeEigenvectors|Ax_lBx) 108 : Base(matA.cols()) 109 { 110 compute(matA, matB, options); 111 } 112 113 /** \brief Computes generalized eigendecomposition of given matrix pencil. 114 * 115 * \param[in] matA Selfadjoint matrix in matrix pencil. 116 * Only the lower triangular part of the matrix is referenced. 117 * \param[in] matB Positive-definite matrix in matrix pencil. 118 * Only the lower triangular part of the matrix is referenced. 119 * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}. 120 * Default is #ComputeEigenvectors|#Ax_lBx. 121 * 122 * \returns Reference to \c *this 123 * 124 * Accoring to \p options, this function computes eigenvalues and (if requested) 125 * the eigenvectors of one of the following three generalized eigenproblems: 126 * - \c Ax_lBx: \f$ Ax = \lambda B x \f$ 127 * - \c ABx_lx: \f$ ABx = \lambda x \f$ 128 * - \c BAx_lx: \f$ BAx = \lambda x \f$ 129 * with \a matA the selfadjoint matrix \f$ A \f$ and \a matB the positive definite 130 * matrix \f$ B \f$. 131 * In addition, each eigenvector \f$ x \f$ satisfies the property \f$ x^* B x = 1 \f$. 132 * 133 * The eigenvalues() function can be used to retrieve 134 * the eigenvalues. If \p options contains ComputeEigenvectors, then the 135 * eigenvectors are also computed and can be retrieved by calling 136 * eigenvectors(). 137 * 138 * The implementation uses LLT to compute the Cholesky decomposition 139 * \f$ B = LL^* \f$ and computes the classical eigendecomposition 140 * of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx 141 * and of \f$ L^{*} A L \f$ otherwise. This solves the 142 * generalized eigenproblem, because any solution of the generalized 143 * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution 144 * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the 145 * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. Similar statements 146 * can be made for the two other variants. 147 * 148 * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp 149 * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out 150 * 151 * \sa GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) 152 */ 153 GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB, 154 int options = ComputeEigenvectors|Ax_lBx); 155 156 protected: 157 158 }; 159 160 161 template<typename MatrixType> 162 GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>:: 163 compute(const MatrixType& matA, const MatrixType& matB, int options) 164 { 165 eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); 166 eigen_assert((options&~(EigVecMask|GenEigMask))==0 167 && (options&EigVecMask)!=EigVecMask 168 && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx 169 || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx) 170 && "invalid option parameter"); 171 172 bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors); 173 174 // Compute the cholesky decomposition of matB = L L' = U'U 175 LLT<MatrixType> cholB(matB); 176 177 int type = (options&GenEigMask); 178 if(type==0) 179 type = Ax_lBx; 180 181 if(type==Ax_lBx) 182 { 183 // compute C = inv(L) A inv(L') 184 MatrixType matC = matA.template selfadjointView<Lower>(); 185 cholB.matrixL().template solveInPlace<OnTheLeft>(matC); 186 cholB.matrixU().template solveInPlace<OnTheRight>(matC); 187 188 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly ); 189 190 // transform back the eigen vectors: evecs = inv(U) * evecs 191 if(computeEigVecs) 192 cholB.matrixU().solveInPlace(Base::m_eivec); 193 } 194 else if(type==ABx_lx) 195 { 196 // compute C = L' A L 197 MatrixType matC = matA.template selfadjointView<Lower>(); 198 matC = matC * cholB.matrixL(); 199 matC = cholB.matrixU() * matC; 200 201 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); 202 203 // transform back the eigen vectors: evecs = inv(U) * evecs 204 if(computeEigVecs) 205 cholB.matrixU().solveInPlace(Base::m_eivec); 206 } 207 else if(type==BAx_lx) 208 { 209 // compute C = L' A L 210 MatrixType matC = matA.template selfadjointView<Lower>(); 211 matC = matC * cholB.matrixL(); 212 matC = cholB.matrixU() * matC; 213 214 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); 215 216 // transform back the eigen vectors: evecs = L * evecs 217 if(computeEigVecs) 218 Base::m_eivec = cholB.matrixL() * Base::m_eivec; 219 } 220 221 return *this; 222 } 223 224 } // end namespace Eigen 225 226 #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H 227
