$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 typename Base::Index Index; 00054 typedef _MatrixType MatrixType; 00055 00063 GeneralizedSelfAdjointEigenSolver() : Base() {} 00064 00077 GeneralizedSelfAdjointEigenSolver(Index size) 00078 : Base(size) 00079 {} 00080 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 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