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 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
00027 #define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
00028 
00029 #include "./Tridiagonalization.h"
00030 
00031 namespace Eigen { 
00032 
00062 template<typename _MatrixType>
00063 class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType>
00064 {
00065     typedef SelfAdjointEigenSolver<_MatrixType> Base;
00066   public:
00067 
00068     typedef typename Base::Index Index;
00069     typedef _MatrixType MatrixType;
00070 
00078     GeneralizedSelfAdjointEigenSolver() : Base() {}
00079 
00092     GeneralizedSelfAdjointEigenSolver(Index size)
00093         : Base(size)
00094     {}
00095 
00122     GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
00123                                       int options = ComputeEigenvectors|Ax_lBx)
00124       : Base(matA.cols())
00125     {
00126       compute(matA, matB, options);
00127     }
00128 
00169     GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB,
00170                                                int options = ComputeEigenvectors|Ax_lBx);
00171 
00172   protected:
00173 
00174 };
00175 
00176 
00177 template<typename MatrixType>
00178 GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
00179 compute(const MatrixType& matA, const MatrixType& matB, int options)
00180 {
00181   eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
00182   eigen_assert((options&~(EigVecMask|GenEigMask))==0
00183           && (options&EigVecMask)!=EigVecMask
00184           && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
00185            || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
00186           && "invalid option parameter");
00187 
00188   bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);
00189 
00190   // Compute the cholesky decomposition of matB = L L' = U'U
00191   LLT<MatrixType> cholB(matB);
00192 
00193   int type = (options&GenEigMask);
00194   if(type==0)
00195     type = Ax_lBx;
00196 
00197   if(type==Ax_lBx)
00198   {
00199     // compute C = inv(L) A inv(L')
00200     MatrixType matC = matA.template selfadjointView<Lower>();
00201     cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
00202     cholB.matrixU().template solveInPlace<OnTheRight>(matC);
00203 
00204     Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );
00205 
00206     // transform back the eigen vectors: evecs = inv(U) * evecs
00207     if(computeEigVecs)
00208       cholB.matrixU().solveInPlace(Base::m_eivec);
00209   }
00210   else if(type==ABx_lx)
00211   {
00212     // compute C = L' A L
00213     MatrixType matC = matA.template selfadjointView<Lower>();
00214     matC = matC * cholB.matrixL();
00215     matC = cholB.matrixU() * matC;
00216 
00217     Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
00218 
00219     // transform back the eigen vectors: evecs = inv(U) * evecs
00220     if(computeEigVecs)
00221       cholB.matrixU().solveInPlace(Base::m_eivec);
00222   }
00223   else if(type==BAx_lx)
00224   {
00225     // compute C = L' A L
00226     MatrixType matC = matA.template selfadjointView<Lower>();
00227     matC = matC * cholB.matrixL();
00228     matC = cholB.matrixU() * matC;
00229 
00230     Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
00231 
00232     // transform back the eigen vectors: evecs = L * evecs
00233     if(computeEigVecs)
00234       Base::m_eivec = cholB.matrixL() * Base::m_eivec;
00235   }
00236 
00237   return *this;
00238 }
00239 
00240 } // end namespace Eigen
00241 
00242 #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H