![]() |
|
|||
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
[ Source navigation ] | [ Diff markup ] | [ Identifier search ] | [ general search ] |
This page was automatically generated by the 2.3.7 LXR engine. The LXR team |
![]() ![]() |