00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk> 00006 // 00007 // Eigen is free software; you can redistribute it and/or 00008 // modify it under the terms of the GNU Lesser General Public 00009 // License as published by the Free Software Foundation; either 00010 // version 3 of the License, or (at your option) any later version. 00011 // 00012 // Alternatively, you can redistribute it and/or 00013 // modify it under the terms of the GNU General Public License as 00014 // published by the Free Software Foundation; either version 2 of 00015 // the License, or (at your option) any later version. 00016 // 00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY 00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the 00020 // GNU General Public License for more details. 00021 // 00022 // You should have received a copy of the GNU Lesser General Public 00023 // License and a copy of the GNU General Public License along with 00024 // Eigen. If not, see <http://www.gnu.org/licenses/>. 00025 00026 #ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H 00027 #define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H 00028 00029 #include "./EigenvaluesCommon.h" 00030 #include "./Tridiagonalization.h" 00031 00032 /** \eigenvalues_module \ingroup Eigenvalues_Module 00033 * 00034 * 00035 * \class GeneralizedSelfAdjointEigenSolver 00036 * 00037 * \brief Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem 00038 * 00039 * \tparam _MatrixType the type of the matrix of which we are computing the 00040 * eigendecomposition; this is expected to be an instantiation of the Matrix 00041 * class template. 00042 * 00043 * This class solves the generalized eigenvalue problem 00044 * \f$ Av = \lambda Bv \f$. In this case, the matrix \f$ A \f$ should be 00045 * selfadjoint and the matrix \f$ B \f$ should be positive definite. 00046 * 00047 * Only the \b lower \b triangular \b part of the input matrix is referenced. 00048 * 00049 * Call the function compute() to compute the eigenvalues and eigenvectors of 00050 * a given matrix. Alternatively, you can use the 00051 * GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) 00052 * constructor which computes the eigenvalues and eigenvectors at construction time. 00053 * Once the eigenvalue and eigenvectors are computed, they can be retrieved with the eigenvalues() 00054 * and eigenvectors() functions. 00055 * 00056 * The documentation for GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) 00057 * contains an example of the typical use of this class. 00058 * 00059 * \sa class SelfAdjointEigenSolver, class EigenSolver, class ComplexEigenSolver 00060 */ 00061 template<typename _MatrixType> 00062 class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType> 00063 { 00064 typedef SelfAdjointEigenSolver<_MatrixType> Base; 00065 public: 00066 00067 typedef typename Base::Index Index; 00068 typedef _MatrixType MatrixType; 00069 00070 /** \brief Default constructor for fixed-size matrices. 00071 * 00072 * The default constructor is useful in cases in which the user intends to 00073 * perform decompositions via compute(const MatrixType&, bool) or 00074 * compute(const MatrixType&, const MatrixType&, bool). This constructor 00075 * can only be used if \p _MatrixType is a fixed-size matrix; use 00076 * SelfAdjointEigenSolver(Index) for dynamic-size matrices. 00077 * 00078 * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp 00079 * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out 00080 */ 00081 GeneralizedSelfAdjointEigenSolver() : Base() {} 00082 00083 /** \brief Constructor, pre-allocates memory for dynamic-size matrices. 00084 * 00085 * \param [in] size Positive integer, size of the matrix whose 00086 * eigenvalues and eigenvectors will be computed. 00087 * 00088 * This constructor is useful for dynamic-size matrices, when the user 00089 * intends to perform decompositions via compute(const MatrixType&, bool) 00090 * or compute(const MatrixType&, const MatrixType&, bool). The \p size 00091 * parameter is only used as a hint. It is not an error to give a wrong 00092 * \p size, but it may impair performance. 00093 * 00094 * \sa compute(const MatrixType&, bool) for an example 00095 */ 00096 GeneralizedSelfAdjointEigenSolver(Index size) 00097 : Base(size) 00098 {} 00099 00100 /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil. 00101 * 00102 * \param[in] matA Selfadjoint matrix in matrix pencil. 00103 * Only the lower triangular part of the matrix is referenced. 00104 * \param[in] matB Positive-definite matrix in matrix pencil. 00105 * Only the lower triangular part of the matrix is referenced. 00106 * \param[in] options A or-ed set of flags {ComputeEigenvectors,EigenvaluesOnly} | {Ax_lBx,ABx_lx,BAx_lx}. 00107 * Default is ComputeEigenvectors|Ax_lBx. 00108 * 00109 * This constructor calls compute(const MatrixType&, const MatrixType&, int) 00110 * to compute the eigenvalues and (if requested) the eigenvectors of the 00111 * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the 00112 * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix 00113 * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property 00114 * \f$ x^* B x = 1 \f$. The eigenvectors are computed if 00115 * \a options contains ComputeEigenvectors. 00116 * 00117 * In addition, the two following variants can be solved via \p options: 00118 * - \c ABx_lx: \f$ ABx = \lambda x \f$ 00119 * - \c BAx_lx: \f$ BAx = \lambda x \f$ 00120 * 00121 * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp 00122 * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out 00123 * 00124 * \sa compute(const MatrixType&, const MatrixType&, int) 00125 */ 00126 GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, 00127 int options = ComputeEigenvectors|Ax_lBx) 00128 : Base(matA.cols()) 00129 { 00130 compute(matA, matB, options); 00131 } 00132 00133 /** \brief Computes generalized eigendecomposition of given matrix pencil. 00134 * 00135 * \param[in] matA Selfadjoint matrix in matrix pencil. 00136 * Only the lower triangular part of the matrix is referenced. 00137 * \param[in] matB Positive-definite matrix in matrix pencil. 00138 * Only the lower triangular part of the matrix is referenced. 00139 * \param[in] options A or-ed set of flags {ComputeEigenvectors,EigenvaluesOnly} | {Ax_lBx,ABx_lx,BAx_lx}. 00140 * Default is ComputeEigenvectors|Ax_lBx. 00141 * 00142 * \returns Reference to \c *this 00143 * 00144 * Accoring to \p options, this function computes eigenvalues and (if requested) 00145 * the eigenvectors of one of the following three generalized eigenproblems: 00146 * - \c Ax_lBx: \f$ Ax = \lambda B x \f$ 00147 * - \c ABx_lx: \f$ ABx = \lambda x \f$ 00148 * - \c BAx_lx: \f$ BAx = \lambda x \f$ 00149 * with \a matA the selfadjoint matrix \f$ A \f$ and \a matB the positive definite 00150 * matrix \f$ B \f$. 00151 * In addition, each eigenvector \f$ x \f$ satisfies the property \f$ x^* B x = 1 \f$. 00152 * 00153 * The eigenvalues() function can be used to retrieve 00154 * the eigenvalues. If \p options contains ComputeEigenvectors, then the 00155 * eigenvectors are also computed and can be retrieved by calling 00156 * eigenvectors(). 00157 * 00158 * The implementation uses LLT to compute the Cholesky decomposition 00159 * \f$ B = LL^* \f$ and computes the classical eigendecomposition 00160 * of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx 00161 * and of \f$ L^{*} A L \f$ otherwise. This solves the 00162 * generalized eigenproblem, because any solution of the generalized 00163 * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution 00164 * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the 00165 * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. Similar statements 00166 * can be made for the two other variants. 00167 * 00168 * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp 00169 * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out 00170 * 00171 * \sa GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) 00172 */ 00173 GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB, 00174 int options = ComputeEigenvectors|Ax_lBx); 00175 00176 protected: 00177 00178 }; 00179 00180 00181 template<typename MatrixType> 00182 GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>:: 00183 compute(const MatrixType& matA, const MatrixType& matB, int options) 00184 { 00185 eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); 00186 eigen_assert((options&~(EigVecMask|GenEigMask))==0 00187 && (options&EigVecMask)!=EigVecMask 00188 && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx 00189 || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx) 00190 && "invalid option parameter"); 00191 00192 bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors); 00193 00194 // Compute the cholesky decomposition of matB = L L' = U'U 00195 LLT<MatrixType> cholB(matB); 00196 00197 int type = (options&GenEigMask); 00198 if(type==0) 00199 type = Ax_lBx; 00200 00201 if(type==Ax_lBx) 00202 { 00203 // compute C = inv(L) A inv(L') 00204 MatrixType matC = matA.template selfadjointView<Lower>(); 00205 cholB.matrixL().template solveInPlace<OnTheLeft>(matC); 00206 cholB.matrixU().template solveInPlace<OnTheRight>(matC); 00207 00208 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly ); 00209 00210 // transform back the eigen vectors: evecs = inv(U) * evecs 00211 if(computeEigVecs) 00212 cholB.matrixU().solveInPlace(Base::m_eivec); 00213 } 00214 else if(type==ABx_lx) 00215 { 00216 // compute C = L' A L 00217 MatrixType matC = matA.template selfadjointView<Lower>(); 00218 matC = matC * cholB.matrixL(); 00219 matC = cholB.matrixU() * matC; 00220 00221 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); 00222 00223 // transform back the eigen vectors: evecs = inv(U) * evecs 00224 if(computeEigVecs) 00225 cholB.matrixU().solveInPlace(Base::m_eivec); 00226 } 00227 else if(type==BAx_lx) 00228 { 00229 // compute C = L' A L 00230 MatrixType matC = matA.template selfadjointView<Lower>(); 00231 matC = matC * cholB.matrixL(); 00232 matC = cholB.matrixU() * matC; 00233 00234 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); 00235 00236 // transform back the eigen vectors: evecs = L * evecs 00237 if(computeEigVecs) 00238 Base::m_eivec = cholB.matrixL() * Base::m_eivec; 00239 } 00240 00241 return *this; 00242 } 00243 00244 #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN:exported at Tue Jan 25 21:56:31 UTC 2011 |