GeneralizedSelfAdjointEigenSolver.h
Go to the documentation of this file.
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


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 11:58:22