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 typename Base::Index Index; 54 typedef _MatrixType MatrixType; 55 56 /** \brief Default constructor for fixed-size matrices. 57 * 58 * The default constructor is useful in cases in which the user intends to 59 * perform decompositions via compute(). This constructor 60 * can only be used if \p _MatrixType is a fixed-size matrix; use 61 * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices. 62 */ 63 GeneralizedSelfAdjointEigenSolver() : Base() {} 64 65 /** \brief Constructor, pre-allocates memory for dynamic-size matrices. 66 * 67 * \param [in] size Positive integer, size of the matrix whose 68 * eigenvalues and eigenvectors will be computed. 69 * 70 * This constructor is useful for dynamic-size matrices, when the user 71 * intends to perform decompositions via compute(). The \p size 72 * parameter is only used as a hint. It is not an error to give a wrong 73 * \p size, but it may impair performance. 74 * 75 * \sa compute() for an example 76 */ 77 GeneralizedSelfAdjointEigenSolver(Index size) 78 : Base(size) 79 {} 80 81 /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil. 82 * 83 * \param[in] matA Selfadjoint matrix in matrix pencil. 84 * Only the lower triangular part of the matrix is referenced. 85 * \param[in] matB Positive-definite matrix in matrix pencil. 86 * Only the lower triangular part of the matrix is referenced. 87 * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}. 88 * Default is #ComputeEigenvectors|#Ax_lBx. 89 * 90 * This constructor calls compute(const MatrixType&, const MatrixType&, int) 91 * to compute the eigenvalues and (if requested) the eigenvectors of the 92 * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the 93 * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix 94 * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property 95 * \f$ x^* B x = 1 \f$. The eigenvectors are computed if 96 * \a options contains ComputeEigenvectors. 97 * 98 * In addition, the two following variants can be solved via \p options: 99 * - \c ABx_lx: \f$ ABx = \lambda x \f$ 100 * - \c BAx_lx: \f$ BAx = \lambda x \f$ 101 * 102 * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp 103 * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out 104 * 105 * \sa compute(const MatrixType&, const MatrixType&, int) 106 */ 107 GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, 108 int options = ComputeEigenvectors|Ax_lBx) 109 : Base(matA.cols()) 110 { 111 compute(matA, matB, options); 112 } 113 114 /** \brief Computes generalized eigendecomposition of given matrix pencil. 115 * 116 * \param[in] matA Selfadjoint matrix in matrix pencil. 117 * Only the lower triangular part of the matrix is referenced. 118 * \param[in] matB Positive-definite matrix in matrix pencil. 119 * Only the lower triangular part of the matrix is referenced. 120 * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}. 121 * Default is #ComputeEigenvectors|#Ax_lBx. 122 * 123 * \returns Reference to \c *this 124 * 125 * Accoring to \p options, this function computes eigenvalues and (if requested) 126 * the eigenvectors of one of the following three generalized eigenproblems: 127 * - \c Ax_lBx: \f$ Ax = \lambda B x \f$ 128 * - \c ABx_lx: \f$ ABx = \lambda x \f$ 129 * - \c BAx_lx: \f$ BAx = \lambda x \f$ 130 * with \a matA the selfadjoint matrix \f$ A \f$ and \a matB the positive definite 131 * matrix \f$ B \f$. 132 * In addition, each eigenvector \f$ x \f$ satisfies the property \f$ x^* B x = 1 \f$. 133 * 134 * The eigenvalues() function can be used to retrieve 135 * the eigenvalues. If \p options contains ComputeEigenvectors, then the 136 * eigenvectors are also computed and can be retrieved by calling 137 * eigenvectors(). 138 * 139 * The implementation uses LLT to compute the Cholesky decomposition 140 * \f$ B = LL^* \f$ and computes the classical eigendecomposition 141 * of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx 142 * and of \f$ L^{*} A L \f$ otherwise. This solves the 143 * generalized eigenproblem, because any solution of the generalized 144 * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution 145 * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the 146 * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. Similar statements 147 * can be made for the two other variants. 148 * 149 * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp 150 * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out 151 * 152 * \sa GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) 153 */ 154 GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB, 155 int options = ComputeEigenvectors|Ax_lBx); 156 157 protected: 158 159 }; 160 161 162 template<typename MatrixType> 163 GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>:: 164 compute(const MatrixType& matA, const MatrixType& matB, int options) 165 { 166 eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); 167 eigen_assert((options&~(EigVecMask|GenEigMask))==0 168 && (options&EigVecMask)!=EigVecMask 169 && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx 170 || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx) 171 && "invalid option parameter"); 172 173 bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors); 174 175 // Compute the cholesky decomposition of matB = L L' = U'U 176 LLT<MatrixType> cholB(matB); 177 178 int type = (options&GenEigMask); 179 if(type==0) 180 type = Ax_lBx; 181 182 if(type==Ax_lBx) 183 { 184 // compute C = inv(L) A inv(L') 185 MatrixType matC = matA.template selfadjointView<Lower>(); 186 cholB.matrixL().template solveInPlace<OnTheLeft>(matC); 187 cholB.matrixU().template solveInPlace<OnTheRight>(matC); 188 189 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly ); 190 191 // transform back the eigen vectors: evecs = inv(U) * evecs 192 if(computeEigVecs) 193 cholB.matrixU().solveInPlace(Base::m_eivec); 194 } 195 else if(type==ABx_lx) 196 { 197 // compute C = L' A L 198 MatrixType matC = matA.template selfadjointView<Lower>(); 199 matC = matC * cholB.matrixL(); 200 matC = cholB.matrixU() * matC; 201 202 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); 203 204 // transform back the eigen vectors: evecs = inv(U) * evecs 205 if(computeEigVecs) 206 cholB.matrixU().solveInPlace(Base::m_eivec); 207 } 208 else if(type==BAx_lx) 209 { 210 // compute C = L' A L 211 MatrixType matC = matA.template selfadjointView<Lower>(); 212 matC = matC * cholB.matrixL(); 213 matC = cholB.matrixU() * matC; 214 215 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); 216 217 // transform back the eigen vectors: evecs = L * evecs 218 if(computeEigVecs) 219 Base::m_eivec = cholB.matrixL() * Base::m_eivec; 220 } 221 222 return *this; 223 } 224 225 } // end namespace Eigen 226 227 #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H 228
