GeneralizedSelfAdjointEigenSolver.h

Go to the documentation of this file.

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.

#ifndef EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H
#define EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H

#include "./EigenvaluesCommon.h"
#include "./Tridiagonalization.h"

template<typename _MatrixType>
class GeneralizedSelfAdjointEigenSolver : public SelfAdjointEigenSolver<_MatrixType>
{
    typedef SelfAdjointEigenSolver<_MatrixType> Base;
  public:

    typedef typename Base::Index Index;
    typedef _MatrixType MatrixType;

    GeneralizedSelfAdjointEigenSolver() : Base() {}

    GeneralizedSelfAdjointEigenSolver(Index size)
        : Base(size)
    {}

    GeneralizedSelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB,
                                      int options = ComputeEigenvectors|Ax_lBx)
      : Base(matA.cols())
    {
      compute(matA, matB, options);
    }

    GeneralizedSelfAdjointEigenSolver& compute(const MatrixType& matA, const MatrixType& matB,
                                               int options = ComputeEigenvectors|Ax_lBx);

  protected:

};


template<typename MatrixType>
GeneralizedSelfAdjointEigenSolver<MatrixType>& GeneralizedSelfAdjointEigenSolver<MatrixType>::
compute(const MatrixType& matA, const MatrixType& matB, int options)
{
  eigen_assert(matA.cols()==matA.rows() && matB.rows()==matA.rows() && matB.cols()==matB.rows());
  eigen_assert((options&~(EigVecMask|GenEigMask))==0
          && (options&EigVecMask)!=EigVecMask
          && ((options&GenEigMask)==0 || (options&GenEigMask)==Ax_lBx
           || (options&GenEigMask)==ABx_lx || (options&GenEigMask)==BAx_lx)
          && "invalid option parameter");

  bool computeEigVecs = ((options&EigVecMask)==0) || ((options&EigVecMask)==ComputeEigenvectors);

  // Compute the cholesky decomposition of matB = L L' = U'U
  LLT<MatrixType> cholB(matB);

  int type = (options&GenEigMask);
  if(type==0)
    type = Ax_lBx;

  if(type==Ax_lBx)
  {
    // compute C = inv(L) A inv(L')
    MatrixType matC = matA.template selfadjointView<Lower>();
    cholB.matrixL().template solveInPlace<OnTheLeft>(matC);
    cholB.matrixU().template solveInPlace<OnTheRight>(matC);

    Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly );

    // transform back the eigen vectors: evecs = inv(U) * evecs
    if(computeEigVecs)
      cholB.matrixU().solveInPlace(Base::m_eivec);
  }
  else if(type==ABx_lx)
  {
    // compute C = L' A L
    MatrixType matC = matA.template selfadjointView<Lower>();
    matC = matC * cholB.matrixL();
    matC = cholB.matrixU() * matC;

    Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);

    // transform back the eigen vectors: evecs = inv(U) * evecs
    if(computeEigVecs)
      cholB.matrixU().solveInPlace(Base::m_eivec);
  }
  else if(type==BAx_lx)
  {
    // compute C = L' A L
    MatrixType matC = matA.template selfadjointView<Lower>();
    matC = matC * cholB.matrixL();
    matC = cholB.matrixU() * matC;

    Base::compute(matC, computeEigVecs ? ComputeEigenvectors : EigenvaluesOnly);

    // transform back the eigen vectors: evecs = L * evecs
    if(computeEigVecs)
      Base::m_eivec = cholB.matrixL() * Base::m_eivec;
  }

  return *this;
}

#endif // EIGEN_GENERALIZEDSELFADJOINTEIGENSOLVER_H