eigensolver/anasaziES.h

Go to the documentation of this file.

#pragma once

#include "anasaziHelper.h"
#include "eigensolver/eigens.hh"
#include "operator.hh"
#include "anasaziHelper.h"
#include "AnasaziBasicEigenproblem.hpp"
#include <AnasaziSolverManager.hpp>
#include <Teuchos_ParameterList.hpp>

namespace concepts {

template <class ScalarT>
class AnasaziES : public eigensolver::EigenSolver<ScalarT> {
public:
  enum State {INIT_PHASE, CONFIGURED_PHASE, SOLVED_PHASE};
  enum SolverType {BLOCK_KRYLOV_SCHUR, BLOCK_DAVIDSON, LOBPCG, RTR};

  typedef Anasazi::MultiVec<ScalarT> MV; // abstract base
  typedef Anasazi::Operator<ScalarT> OP; // abstract base
  typedef AnasaziMV<ScalarT> MV_C; // concrete implementation 
  typedef AnasaziOp<ScalarT> OP_C; // concrete implementation 
  typedef Anasazi::MultiVecTraits<ScalarT, MV> MVT;

  AnasaziES(SolverType solverType,
      concepts::Operator<ScalarT>& A
      );

  ~AnasaziES() {

  }

  void setM(concepts::Operator<ScalarT>& M) {
    assert(phase_ == INIT_PHASE);
    myProblem_->setM(Teuchos::rcp( new OP_C(Teuchos::rcp(&M, false) ) ));
  }

  void setHermitian(bool isHermite) {
    assert(phase_ == INIT_PHASE);
    myProblem_->setHermitian(isHermite);
  }

  void setBlockSize(int blocksize = 1) {
    myPL_.set( "Block Size", blocksize );
    int dim = A_->getOp().dimX();
     
    Teuchos::RCP<MV> ivec = Teuchos::rcp( new MV_C(dim,blocksize) );
    ivec->MvRandom();
    myProblem_->setInitVec(ivec);
  }

  void setNEV(int nev=3) {
    assert(phase_ == INIT_PHASE);
    myProblem_->setNEV(nev);
  }

  bool setProblem();

  int solve();

  /*
   Specify the verbosity level.  Options include:
   Anasazi::Errors 
     This option is always set
   Anasazi::Warnings 
     Warnings (less severe than errors)
   Anasazi::IterationDetails 
     Details at each iteration, such as the current eigenvalues
   Anasazi::OrthoDetails 
     Details about orthogonality
   Anasazi::TimingDetails
     A summary of the timing info for the solve() routine
   Anasazi::FinalSummary 
     A final summary 
   Anasazi::Debug 
     Debugging information
     */
  void setVerbosity(int verbosity = Anasazi::Errors) {
    assert(phase_ == INIT_PHASE);
    myPL_.set( "Verbosity", verbosity );
  }

  void setOrthogonalization(std::string mode = "SVQB") {
    assert(phase_ == INIT_PHASE);
    myPL_.set( "Orthogonalization", mode );
  }

  /*
   Choose which eigenvalues to compute.
   Choices are:
   LM - target the largest magnitude  [default]
   SM - target the smallest magnitude 
   LR - target the largest real 
   SR - target the smallest real 
   LI - target the largest imaginary
   SI - target the smallest imaginary
   */
  void setMode(std::string mode = "LM") {
    assert(phase_ == INIT_PHASE);
    myPL_.set( "Which", mode );
  }

  void setNumBlocks(int numBlocks) {
    assert(phase_ == INIT_PHASE);
    myPL_.set( "Num Blocks", numBlocks);
  }

  void setMaxRestarts(int maxRestarts = 100) {
    assert(phase_ == INIT_PHASE);
    myPL_.set( "Maximum Restarts", maxRestarts );
  }

  void setConvTol(double tolerance) {
    assert(phase_ == INIT_PHASE);
    myPL_.set( "Convergence Tolerance",  tolerance);
  }

  void setRelConvTol(bool rct) {
    assert(phase_ == INIT_PHASE);
    myPL_.set( "Relative Convergence Tolerance",  rct);
  }

  void setConvergenceNorm(std::string norm = "2") {
    assert(phase_ == INIT_PHASE);
    myPL_.set( "Convergence Norm", norm );
  }

  void setMaxLocked(int val) {
    assert(phase_ == INIT_PHASE);
    myPL_.set( "Max Locked", val);
  }

  virtual const concepts::Array<ScalarT>& getEV() {
    if(phase_ != SOLVED_PHASE) {
      int res = solve();
      assert(res == 0);
    }

    return evs_;
  }

  virtual const concepts::Array<concepts::Vector<ScalarT>*>& getEF() 
  {
    if(phase_ != SOLVED_PHASE) {
      int res = solve();
      assert(res == 0);
    }

    return efs_;
  }

  virtual uint converged() const 
  {
    if(phase_ != SOLVED_PHASE)
      return 0;

    const Anasazi::Eigensolution<ScalarT, MV>& sol = myProblem_->getSolution();
    return sol.numVecs;
  }

  virtual uint iterations() const {
     if(phase_ != SOLVED_PHASE)
      return -1;

    return solver_->getNumIters();
  }

  Teuchos::RCP< Anasazi::BasicEigenproblem<ScalarT, MV, OP> > getProblem()
  {
    return myProblem_;
  }

public: 
  virtual std::ostream& info(std::ostream& os) const {
     return os << "AnasaziES() in state " << (int)phase_;
  }

private:
  Teuchos::RCP<OP_C> A_;
  Teuchos::RCP< Anasazi::BasicEigenproblem<ScalarT, MV, OP> > myProblem_;

  Teuchos::ParameterList myPL_;

  Teuchos::RCP< Anasazi::SolverManager<ScalarT, MV, OP> > solver_;
  
  concepts::Array<ScalarT> evs_;
  concepts::Array<concepts::Vector<ScalarT>*> efs_;

  State phase_;
  SolverType solverType_;
};

}