Important changes to repositories hosted on mbed.com
Mbed hosted mercurial repositories are deprecated and are due to be permanently deleted in July 2026.
To keep a copy of this software download the repository Zip archive or clone locally using Mercurial.
It is also possible to export all your personal repositories from the account settings page.
Dependents: Eigen_test Odometry_test AttitudeEstimation_usingTicker MPU9250_Quaternion_Binary_Serial ... more
GeneralizedSelfAdjointEigenSolver.h
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 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H 00012 #define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H 00013 00014 #include "./Tridiagonalization.h" 00015 00016 namespace Eigen { 00017 00018 /** \eigenvalues_module \ingroup Eigenvalues_Module 00019 * 00020 * 00021 * \class GeneralizedSelfAdjointEigenSolver 00022 * 00023 * \brief Computes eigenvalues and eigenvectors of the generalized selfadjoint eigen problem 00024 * 00025 * \tparam _MatrixType the type of the matrix of which we are computing the 00026 * eigendecomposition; this is expected to be an instantiation of the Matrix 00027 * class template. 00028 * 00029 * This class solves the generalized eigenvalue problem 00030 * \f$ Av = \lambda Bv \f$. In this case, the matrix \f$ A \f$ should be 00031 * selfadjoint and the matrix \f$ B \f$ should be positive definite. 00032 * 00033 * Only the \b lower \b triangular \b part of the input matrix is referenced. 00034 * 00035 * Call the function compute() to compute the eigenvalues and eigenvectors of 00036 * a given matrix. Alternatively, you can use the 00037 * GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) 00038 * constructor which computes the eigenvalues and eigenvectors at construction time. 00039 * Once the eigenvalue and eigenvectors are computed, they can be retrieved with the eigenvalues() 00040 * and eigenvectors() functions. 00041 * 00042 * The documentation for GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) 00043 * contains an example of the typical use of this class. 00044 * 00045 * \sa class SelfAdjointEigenSolver, class EigenSolver, class ComplexEigenSolver 00046 */ 00047 template<typename _MatrixType> 00048 class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver <_MatrixType> 00049 { 00050 typedef SelfAdjointEigenSolver<_MatrixType> Base ; 00051 public: 00052 00053 typedef typename Base::Index Index; 00054 typedef _MatrixType MatrixType; 00055 00056 /** \brief Default constructor for fixed-size matrices. 00057 * 00058 * The default constructor is useful in cases in which the user intends to 00059 * perform decompositions via compute(). This constructor 00060 * can only be used if \p _MatrixType is a fixed-size matrix; use 00061 * GeneralizedSelfAdjointEigenSolver(Index) for dynamic-size matrices. 00062 */ 00063 GeneralizedSelfAdjointEigenSolver() : Base () {} 00064 00065 /** \brief Constructor, pre-allocates memory for dynamic-size matrices. 00066 * 00067 * \param [in] size Positive integer, size of the matrix whose 00068 * eigenvalues and eigenvectors will be computed. 00069 * 00070 * This constructor is useful for dynamic-size matrices, when the user 00071 * intends to perform decompositions via compute(). The \p size 00072 * parameter is only used as a hint. It is not an error to give a wrong 00073 * \p size, but it may impair performance. 00074 * 00075 * \sa compute() for an example 00076 */ 00077 GeneralizedSelfAdjointEigenSolver(Index size) 00078 : Base (size) 00079 {} 00080 00081 /** \brief Constructor; computes generalized eigendecomposition of given matrix pencil. 00082 * 00083 * \param[in] matA Selfadjoint matrix in matrix pencil. 00084 * Only the lower triangular part of the matrix is referenced. 00085 * \param[in] matB Positive-definite matrix in matrix pencil. 00086 * Only the lower triangular part of the matrix is referenced. 00087 * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}. 00088 * Default is #ComputeEigenvectors|#Ax_lBx. 00089 * 00090 * This constructor calls compute(const MatrixType&, const MatrixType&, int) 00091 * to compute the eigenvalues and (if requested) the eigenvectors of the 00092 * generalized eigenproblem \f$ Ax = \lambda B x \f$ with \a matA the 00093 * selfadjoint matrix \f$ A \f$ and \a matB the positive definite matrix 00094 * \f$ B \f$. Each eigenvector \f$ x \f$ satisfies the property 00095 * \f$ x^* B x = 1 \f$. The eigenvectors are computed if 00096 * \a options contains ComputeEigenvectors. 00097 * 00098 * In addition, the two following variants can be solved via \p options: 00099 * - \c ABx_lx: \f$ ABx = \lambda x \f$ 00100 * - \c BAx_lx: \f$ BAx = \lambda x \f$ 00101 * 00102 * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.cpp 00103 * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType2.out 00104 * 00105 * \sa compute(const MatrixType&, const MatrixType&, int) 00106 */ 00107 GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, 00108 int options = ComputeEigenvectors|Ax_lBx) 00109 : Base (matA.cols()) 00110 { 00111 compute(matA, matB, options); 00112 } 00113 00114 /** \brief Computes generalized eigendecomposition of given matrix pencil. 00115 * 00116 * \param[in] matA Selfadjoint matrix in matrix pencil. 00117 * Only the lower triangular part of the matrix is referenced. 00118 * \param[in] matB Positive-definite matrix in matrix pencil. 00119 * Only the lower triangular part of the matrix is referenced. 00120 * \param[in] options A or-ed set of flags {#ComputeEigenvectors,#EigenvaluesOnly} | {#Ax_lBx,#ABx_lx,#BAx_lx}. 00121 * Default is #ComputeEigenvectors|#Ax_lBx. 00122 * 00123 * \returns Reference to \c *this 00124 * 00125 * Accoring to \p options, this function computes eigenvalues and (if requested) 00126 * the eigenvectors of one of the following three generalized eigenproblems: 00127 * - \c Ax_lBx: \f$ Ax = \lambda B x \f$ 00128 * - \c ABx_lx: \f$ ABx = \lambda x \f$ 00129 * - \c BAx_lx: \f$ BAx = \lambda x \f$ 00130 * with \a matA the selfadjoint matrix \f$ A \f$ and \a matB the positive definite 00131 * matrix \f$ B \f$. 00132 * In addition, each eigenvector \f$ x \f$ satisfies the property \f$ x^* B x = 1 \f$. 00133 * 00134 * The eigenvalues() function can be used to retrieve 00135 * the eigenvalues. If \p options contains ComputeEigenvectors, then the 00136 * eigenvectors are also computed and can be retrieved by calling 00137 * eigenvectors(). 00138 * 00139 * The implementation uses LLT to compute the Cholesky decomposition 00140 * \f$ B = LL^* \f$ and computes the classical eigendecomposition 00141 * of the selfadjoint matrix \f$ L^{-1} A (L^*)^{-1} \f$ if \p options contains Ax_lBx 00142 * and of \f$ L^{*} A L \f$ otherwise. This solves the 00143 * generalized eigenproblem, because any solution of the generalized 00144 * eigenproblem \f$ Ax = \lambda B x \f$ corresponds to a solution 00145 * \f$ L^{-1} A (L^*)^{-1} (L^* x) = \lambda (L^* x) \f$ of the 00146 * eigenproblem for \f$ L^{-1} A (L^*)^{-1} \f$. Similar statements 00147 * can be made for the two other variants. 00148 * 00149 * Example: \include SelfAdjointEigenSolver_compute_MatrixType2.cpp 00150 * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType2.out 00151 * 00152 * \sa GeneralizedSelfAdjointEigenSolver(const MatrixType&, const MatrixType&, int) 00153 */ 00154 GeneralizedSelfAdjointEigenSolver & compute(const MatrixType& matA, const MatrixType& matB, 00155 int options = ComputeEigenvectors|Ax_lBx); 00156 00157 protected: 00158 00159 }; 00160 00161 00162 template<typename MatrixType> 00163 GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>:: 00164 compute(const MatrixType& matA, const MatrixType& matB, int options) 00165 { 00166 eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); 00167 eigen_assert((options&~(EigVecMask|GenEigMask))==0 00168 && (options&EigVecMask)!=EigVecMask 00169 && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx 00170 || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx) 00171 && "invalid option parameter"); 00172 00173 bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors); 00174 00175 // Compute the cholesky decomposition of matB = L L' = U'U 00176 LLT<MatrixType> cholB(matB); 00177 00178 int type = (options&GenEigMask); 00179 if(type==0) 00180 type = Ax_lBx; 00181 00182 if(type==Ax_lBx) 00183 { 00184 // compute C = inv(L) A inv(L') 00185 MatrixType matC = matA.template selfadjointView<Lower>(); 00186 cholB.matrixL ().template solveInPlace<OnTheLeft>(matC); 00187 cholB.matrixU ().template solveInPlace<OnTheRight>(matC); 00188 00189 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly ); 00190 00191 // transform back the eigen vectors: evecs = inv(U) * evecs 00192 if(computeEigVecs) 00193 cholB.matrixU ().solveInPlace(Base::m_eivec); 00194 } 00195 else if(type==ABx_lx) 00196 { 00197 // compute C = L' A L 00198 MatrixType matC = matA.template selfadjointView<Lower>(); 00199 matC = matC * cholB.matrixL (); 00200 matC = cholB.matrixU () * matC; 00201 00202 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); 00203 00204 // transform back the eigen vectors: evecs = inv(U) * evecs 00205 if(computeEigVecs) 00206 cholB.matrixU ().solveInPlace(Base::m_eivec); 00207 } 00208 else if(type==BAx_lx) 00209 { 00210 // compute C = L' A L 00211 MatrixType matC = matA.template selfadjointView<Lower>(); 00212 matC = matC * cholB.matrixL (); 00213 matC = cholB.matrixU () * matC; 00214 00215 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); 00216 00217 // transform back the eigen vectors: evecs = L * evecs 00218 if(computeEigVecs) 00219 Base::m_eivec = cholB.matrixL () * Base::m_eivec; 00220 } 00221 00222 return *this; 00223 } 00224 00225 } // end namespace Eigen 00226 00227 #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
Generated on Tue Jul 12 2022 17:46:54 by
