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 "./EigenvaluesCommon.h"
00030 #include "./Tridiagonalization.h"
00031 
00061 template<typename _MatrixType>
00062 class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType>
00063 {
00064     typedef SelfAdjointEigenSolver<_MatrixType> Base;
00065   public:
00066 
00067     typedef typename Base::Index Index;
00068     typedef _MatrixType MatrixType;
00069 
00077     GeneralizedSelfAdjointEigenSolver() : Base() {}
00078 
00091     GeneralizedSelfAdjointEigenSolver(Index size)
00092         : Base(size)
00093     {}
00094 
00121     GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
00122                                       int options = ComputeEigenvectors|Ax_lBx)
00123       : Base(matA.cols())
00124     {
00125       compute(matA, matB, options);
00126     }
00127 
00168     GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB,
00169                                                int options = ComputeEigenvectors|Ax_lBx);
00170 
00171   protected:
00172 
00173 };
00174 
00175 
00176 template<typename MatrixType>
00177 GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
00178 compute(const MatrixType& matA, const MatrixType& matB, int options)
00179 {
00180   eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
00181   eigen_assert((options&~(EigVecMask|GenEigMask))==0
00182           && (options&EigVecMask)!=EigVecMask
00183           && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
00184            || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
00185           && "invalid option parameter");
00186 
00187   bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);
00188 
00189   // Compute the cholesky decomposition of matB = L L' = U'U
00190   LLT<MatrixType> cholB(matB);
00191 
00192   int type = (options&GenEigMask);
00193   if(type==0)
00194     type = Ax_lBx;
00195 
00196   if(type==Ax_lBx)
00197   {
00198     // compute C = inv(L) A inv(L')
00199     MatrixType matC = matA.template selfadjointView<Lower>();
00200     cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
00201     cholB.matrixU().template solveInPlace<OnTheRight>(matC);
00202 
00203     Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );
00204 
00205     // transform back the eigen vectors: evecs = inv(U) * evecs
00206     if(computeEigVecs)
00207       cholB.matrixU().solveInPlace(Base::m_eivec);
00208   }
00209   else if(type==ABx_lx)
00210   {
00211     // compute C = L' A L
00212     MatrixType matC = matA.template selfadjointView<Lower>();
00213     matC = matC * cholB.matrixL();
00214     matC = cholB.matrixU() * matC;
00215 
00216     Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
00217 
00218     // transform back the eigen vectors: evecs = inv(U) * evecs
00219     if(computeEigVecs)
00220       cholB.matrixU().solveInPlace(Base::m_eivec);
00221   }
00222   else if(type==BAx_lx)
00223   {
00224     // compute C = L' A L
00225     MatrixType matC = matA.template selfadjointView<Lower>();
00226     matC = matC * cholB.matrixL();
00227     matC = cholB.matrixU() * matC;
00228 
00229     Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);
00230 
00231     // transform back the eigen vectors: evecs = L * evecs
00232     if(computeEigVecs)
00233       Base::m_eivec = cholB.matrixL() * Base::m_eivec;
00234   }
00235 
00236   return *this;
00237 }
00238 
00239 #endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:44