gevp_solver.hh

Go to the documentation of this file.

 
#ifndef EIGENSOLVER_GEVP_SOLVER_HH_
#define EIGENSOLVER_GEVP_SOLVER_HH_
 
#include "arpackpp.hh"
#include "operator.hh"
#include "basics/typedefs.hh"
 
namespace concepts
{
 
  template<class F>
  class ShiftAndInvertOperatorForGEVPs: public VecOperator<F>
  {
 
  public:
    ShiftAndInvertOperatorForGEVPs(SparseMatrix<F> &A, SparseMatrix<F> &B,
           F shift) :
      VecOperator<F>(A.nofRows(), A.nofCols()), _B(B), _shift(shift)
    {
      SparseMatrix<F> shifted(A);
      B.addInto(shifted, -shift); // shifted_matrix = A - shift*B
#ifdef HAS_MUMPS
      _shift_and_invert_pointer.reset(new Mumps<F>(shifted));
#else      
      _shift_and_invert_pointer.reset(new SuperLU<F>(shifted));
#endif      
    }
 
    virtual ~ShiftAndInvertOperatorForGEVPs()
    {
    }
 
    F getShift();
 
  protected:
    virtual std::ostream&
    info(std::ostream& os) const
    {
      // TODO ergaenzen
      return os << concepts::typeOf(*this);
    }
 
    virtual void apply_(const Vector<F>& fncY, Vector<F>& fncX)
    {
      Vector<F> B_times_y(fncY.n());
      _B(fncY, B_times_y);
      (*_shift_and_invert_pointer)(B_times_y, fncX);
    }
 
    virtual void apply_()
    {
    }
 
  private:
    SparseMatrix<F>& _B;
#ifdef HAS_MUMPS
    std::unique_ptr<Mumps<F> > _shift_and_invert_pointer;
#else
    std::unique_ptr<SuperLU<F> > _shift_and_invert_pointer;
#endif      
    F _shift;
 
  };
 
}
 
using concepts::SparseMatrix;
 
namespace eigensolver
{
 
  template<typename F>
  class GEVPSolver: public EigenSolver<Cmplx>
  {
 
  public:
    GEVPSolver(SparseMatrix<F>& A, SparseMatrix<F>& B, F shift, int kmax = 6,
         Real tol = EPS,int maxiter=300):
      standard_evp_solver(*(new concepts::ShiftAndInvertOperatorForGEVPs<F>(A,B,shift)),kmax,(char*)"LM",tol,maxiter),
      _shift(shift)
    {
    }
 
    GEVPSolver(concepts::ShiftAndInvertOperatorForGEVPs<F> OP, int kmax = 6,
         Real tol = EPS,int maxiter=300):
      standard_evp_solver(OP,kmax,(char*)"LM",tol,maxiter), _shift(OP.getShift())
    {
    }
 
    virtual ~GEVPSolver()
    {
    }
 
    virtual const concepts::Array<Cmplx>& getEV();
 
    virtual const concepts::Array<concepts::Vector<Cmplx>*>& getEF();
 
    virtual uint converged() const;
 
    virtual uint iterations() const;
 
  protected:
    virtual std::ostream& info(std::ostream& os) const
    {
      return os << concepts::typeOf(*this);
    }
 
  private:
    ArPackppStd<F> standard_evp_solver;
    F _shift;
    concepts::Array<Cmplx> _unshifted_eigenvalues;
 
  };
 
}
 
#endif /* EIGENSOLVER_GEVP_SOLVER_HH_ */