Eigen  3.3.3
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 
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
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends