![]() |
Eigen
3.3.3
|
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 00047 template<typename _MatrixType> 00048 class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType> 00049 { 00050 typedef SelfAdjointEigenSolver<_MatrixType> Base; 00051 public: 00052 00053 typedef _MatrixType MatrixType; 00054 00062 GeneralizedSelfAdjointEigenSolver() : Base() {} 00063 00076 explicit GeneralizedSelfAdjointEigenSolver(Index size) 00077 : Base(size) 00078 {} 00079 00106 GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, 00107 int options = ComputeEigenvectors|Ax_lBx) 00108 : Base(matA.cols()) 00109 { 00110 compute(matA, matB, options); 00111 } 00112 00153 GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB, 00154 int options = ComputeEigenvectors|Ax_lBx); 00155 00156 protected: 00157 00158 }; 00159 00160 00161 template<typename MatrixType> 00162 GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>:: 00163 compute(const MatrixType& matA, const MatrixType& matB, int options) 00164 { 00165 eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows()); 00166 eigen_assert((options&~(EigVecMask|GenEigMask))==0 00167 && (options&EigVecMask)!=EigVecMask 00168 && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx 00169 || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx) 00170 && "invalid option parameter"); 00171 00172 bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors); 00173 00174 // Compute the cholesky decomposition of matB = L L' = U'U 00175 LLT<MatrixType> cholB(matB); 00176 00177 int type = (options&GenEigMask); 00178 if(type==0) 00179 type = Ax_lBx; 00180 00181 if(type==Ax_lBx) 00182 { 00183 // compute C = inv(L) A inv(L') 00184 MatrixType matC = matA.template selfadjointView<Lower>(); 00185 cholB.matrixL().template solveInPlace<OnTheLeft>(matC); 00186 cholB.matrixU().template solveInPlace<OnTheRight>(matC); 00187 00188 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly ); 00189 00190 // transform back the eigen vectors: evecs = inv(U) * evecs 00191 if(computeEigVecs) 00192 cholB.matrixU().solveInPlace(Base::m_eivec); 00193 } 00194 else if(type==ABx_lx) 00195 { 00196 // compute C = L' A L 00197 MatrixType matC = matA.template selfadjointView<Lower>(); 00198 matC = matC * cholB.matrixL(); 00199 matC = cholB.matrixU() * matC; 00200 00201 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); 00202 00203 // transform back the eigen vectors: evecs = inv(U) * evecs 00204 if(computeEigVecs) 00205 cholB.matrixU().solveInPlace(Base::m_eivec); 00206 } 00207 else if(type==BAx_lx) 00208 { 00209 // compute C = L' A L 00210 MatrixType matC = matA.template selfadjointView<Lower>(); 00211 matC = matC * cholB.matrixL(); 00212 matC = cholB.matrixU() * matC; 00213 00214 Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly); 00215 00216 // transform back the eigen vectors: evecs = L * evecs 00217 if(computeEigVecs) 00218 Base::m_eivec = cholB.matrixL() * Base::m_eivec; 00219 } 00220 00221 return *this; 00222 } 00223 00224 } // end namespace Eigen 00225 00226 #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H