easyArpackpp.hh

Go to the documentation of this file.

 
#ifndef EASYARPACKPP_HH_
#define EASYARPACKPP_HH_
 
#include "arpackpp.hh"
#include "basics/exceptions.hh"
 
namespace eigensolver {
 
  /*** struct to compute the type of the operator for the EasyArpack classes ***/
  template<typename F, typename G>
  struct OperatorType {
    typedef concepts::Cmplx type;
  };
 
  template<>
  struct OperatorType<concepts::Real, concepts::Real> {
    typedef concepts::Real type;
  };
 
  template<class T>
  class EasyArPackpp {
 
  public:
    EasyArPackpp() {
    }
    ;
 
    virtual ~EasyArPackpp() {
    }
    ;
 
    virtual EigenSolver<T>* getSolver() = 0;
 
  protected:
    int dim_;
  };
 
  template<class F = Cmplx, class G = Real,
      class H = typename eigensolver::OperatorType<F, G>::type>
  class EasyArPackppStd: public EasyArPackpp<concepts::Cmplx> {
  public:
 
    EasyArPackppStd(concepts::SparseMatrix<F> &A, int kmax, G shift) {
 
      //set dimension
      dim_ = A.dimX();
 
      //calculate shifted matrix: shifted = A - shift*I
      concepts::SparseMatrix<H> shifted(dim_, dim_);
      A.addInto(shifted, 1.);
      concepts::DiagonalMatrix<concepts::Real> eye(dim_);
      for (int i = 0; i < dim_; ++i)
        eye(i, i) = 1;
      eye.addInto(shifted, -shift);
 
      //set operator
#ifdef HAS_MUMPS      
      OP_ = new concepts::Mumps<H>(shifted);
#else
#ifdef HAS_SuperLU
      OP_ = new concepts::SuperLU<H>(shifted);
#else
    throw conceptsException(concepts::ExceptionBase("No valid solver was defined"));
#endif
#endif
      
      //set solver
      solver_ = new ArPackppStd<H> (*OP_, std::min(std::max(kmax, 1), dim_ - 2), shift);
    }
    ;
 
    EasyArPackppStd(concepts::SparseMatrix<F> &A, int kmax = 1, char* which =
        (char*) "LM") {
 
      //ensure we're not calling with SM
      if ( strcmp(which,"SM")==0 )
        throw(concepts::MissingFeature("SM mode is not allowed. Use the shift and invert method instead."));
      
      //set dimension
      dim_ = A.dimX();
      
      //set empty operator
      OP_ = 0;
      
      //set solver
      solver_ = new ArPackppStd<H> (A, std::min(std::max(kmax, 1), dim_ - 2), which);
      
    }
    ;
 
    virtual ~EasyArPackppStd() {
      if (OP_ != 0)
        delete OP_;
      delete solver_;
    }
    ;
 
    virtual ArPackppStd<H>* getSolver() {
      return solver_;
    }
    ;
 
  private:
    concepts::Operator<H>* OP_;
    ArPackppStd<H>* solver_;
  };
 
  template<class F, class G = concepts::Real,
      class H = typename eigensolver::OperatorType<F, G>::type>
  class EasyArPackppGen: public EasyArPackpp<concepts::Cmplx> {
  public:
 
    EasyArPackppGen(concepts::SparseMatrix<F> &A, concepts::SparseMatrix<
        concepts::Real> &B, int kmax, G shift) {
      //set dimension
      dim_ = A.dimX();
 
      //calculate shifted matrix: shifted = A - shift*B
      concepts::SparseMatrix<H> shifted(dim_, dim_);
      A.addInto(shifted, 1.);
      B.addInto(shifted, -shift);
      
      //set operator
#ifdef HAS_MUMPS      
      OP_ = new concepts::Mumps<H>(shifted);
#else
#ifdef HAS_SuperLU
      OP_ = new concepts::SuperLU<H>(shifted);
#else
    throw conceptsException(concepts::ExceptionBase("No valid solver was defined"));
#endif
#endif
      
      //set solver
      solver_ = new ArPackppGen<H, F, concepts::Real> (*OP_, A, B,
          std::min(std::max(kmax, 1), dim_ - 2), (concepts::Cmplx) shift);
    }
 
    EasyArPackppGen(concepts::SparseMatrix<F> &A, concepts::SparseMatrix<
        concepts::Real> &B, int kmax = 1, char* which = (char*) "LM") {
 
      //ensure we're not calling with SM
      if ( strcmp(which,"SM")==0 )
        throw(concepts::MissingFeature("SM mode is not allowed. Use the shift and invert method instead."));
      
      //set dimension
      dim_ = A.dimX();
 
      /********************UNSECURE CAST**************************
       *
       * In The case H = concepts::Real the cast is well defined.
       * If H = concepts::Cmplx the cast is not well defined and it can occur many problems.
       * But if H = concepts::Cmplx the only method we use from concepts::SuperLU<concepts::Real> is
       * operator()(vector<concepts::Cmplx> v1, vector<concepts::Cmplx> v2) and that is well defined for
       * an instance of concepts::SuperLU<concepts::Real>. Never Use
       * operator()(vector<concepts::Real> v1, vector<concepts::Cmplx> v2) cause this method is not implemented in
       * instances of Operators<concepts::Real>.
       * The class member OP_ is private, so nobody can kid around with it (exept in the class).
       *
       */
      
      //set operator
#ifdef HAS_MUMPS      
      OP_ = (concepts::Operator<H>*) new concepts::Mumps<concepts::Real>(B);
#else
#ifdef HAS_SuperLU
      OP_ = (concepts::Operator<H>*) new concepts::SuperLU<concepts::Real>(B);
#else
    throw conceptsException(concepts::ExceptionBase("No valid solver was defined"));
#endif
#endif
    
      //set solver
      solver_ = new ArPackppGen<H, F, concepts::Real> (*OP_, A, B,
          std::min(std::max(kmax, 1), dim_ - 2), which);
    }
    ;
 
    virtual ~EasyArPackppGen() {
      if (OP_ != 0)
        delete OP_;
      delete solver_;
 
    }
    ;
 
    virtual ArPackppGen<H, F, concepts::Real>* getSolver() {
      return solver_;
    }
    ;
 
  private:
    concepts::Operator<H>* OP_;
    ArPackppGen<H, F, concepts::Real>* solver_;
  };
 
  class EasyArPackppSymGen: public EasyArPackpp<concepts::Real> {
  public:
 
    EasyArPackppSymGen(concepts::SparseMatrix<concepts::Real> &A,
        concepts::SparseMatrix<concepts::Real> &B, int kmax,
        concepts::Real shift);
 
    EasyArPackppSymGen(char invertType,
        concepts::SparseMatrix<concepts::Real> &A, concepts::SparseMatrix<
            concepts::Real> &B, int kmax, concepts::Real shift);
 
    /******************* Not working yet ******************************/
    //    /**Regular mode constructor. Builds an eigenvalue solver
    //     * that solves the given eigenvalue problem using arpack++ with the
    //     * regular mode.
    //     *@param A symmetric matrix of the left hand side
    //     *@param B symmetric matrix of the right hand side; it has to be symmetric positive definite
    //     *@param kmax number of eigenvalues that should be calculated
    //     *@param which defines which eigenvalue should be calculated (default "LM")
    //     */
    //    EasyArPackppSymGen(concepts::SparseMatrix<concepts::Real> &A,
    //        concepts::SparseMatrix<concepts::Real> &B, int kmax,
    //        char* which = (char*)("LM"));
 
    virtual ~EasyArPackppSymGen();
 
    virtual ArPackppSymGen* getSolver();
 
  private:
    concepts::Operator<concepts::Real>* OP_;
    ArPackppSymGen* solver_;
  };
 
}//end namespace eigensolver
 
#endif /* EASYARPACKPP_HH_ */