AnasaziRTRSolMgr.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//                 Anasazi: Block Eigensolvers Package
//                 Copyright (2004) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library 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 2.1 of the
// License, or (at your option) any later version.
//
// This library 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 for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
// @HEADER

#ifndef ANASAZI_RTR_SOLMGR_HPP
#define ANASAZI_RTR_SOLMGR_HPP

#include "AnasaziConfigDefs.hpp"
#include "AnasaziTypes.hpp"

#include "AnasaziEigenproblem.hpp"
#include "AnasaziSolverManager.hpp"
#include "AnasaziSolverUtils.hpp"

#include "AnasaziIRTR.hpp"
#include "AnasaziSIRTR.hpp"
#include "AnasaziBasicSort.hpp"
#include "AnasaziICGSOrthoManager.hpp"
#include "AnasaziStatusTestMaxIters.hpp"
#include "AnasaziStatusTestResNorm.hpp"
#include "AnasaziStatusTestWithOrdering.hpp"
#include "AnasaziStatusTestCombo.hpp"
#include "AnasaziStatusTestOutput.hpp"
#include "AnasaziBasicOutputManager.hpp"

#include <Teuchos_TimeMonitor.hpp>
#include <Teuchos_FancyOStream.hpp>

namespace Anasazi {

template<class ScalarType, class MV, class OP>
class RTRSolMgr : public SolverManager<ScalarType,MV,OP> {

  private:
    typedef MultiVecTraits<ScalarType,MV> MVT;
    typedef OperatorTraits<ScalarType,MV,OP> OPT;
    typedef Teuchos::ScalarTraits<ScalarType> SCT;
    typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
    typedef Teuchos::ScalarTraits<MagnitudeType> MT;

  public:



  RTRSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
             Teuchos::ParameterList &pl );

  virtual ~RTRSolMgr() {};



  const Eigenproblem<ScalarType,MV,OP>& getProblem() const {
    return *problem_;
  }

   Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
     return Teuchos::tuple(_timerSolve);
   }

  int getNumIters() const {
    return numIters_;
  }





  ReturnType solve();

  private:
  Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
  std::string whch_;
  std::string ortho_;

  bool skinny_;
  MagnitudeType convtol_;
  int maxIters_;
  bool relconvtol_;
  enum ResType convNorm_;
  int numIters_;
  int numICGS_;

  Teuchos::RCP<Teuchos::Time> _timerSolve;
  Teuchos::RCP<BasicOutputManager<ScalarType> > printer_;
  Teuchos::ParameterList pl_;
};


template<class ScalarType, class MV, class OP>
RTRSolMgr<ScalarType,MV,OP>::RTRSolMgr(
        const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
        Teuchos::ParameterList &pl ) :
  problem_(problem),
  whch_("SR"),
  skinny_(true),
  convtol_(MT::prec()),
  maxIters_(100),
  relconvtol_(true),
  numIters_(-1),
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
  _timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: RTRSolMgr::solve()")),
#endif
  pl_(pl)
{
  TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,              std::invalid_argument, "Problem not given to solver manager.");
  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),              std::invalid_argument, "Problem not set.");
  TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(),               std::invalid_argument, "Problem not symmetric.");
  TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");

  std::string strtmp;

  whch_ = pl_.get("Which","SR");
  TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SR" && whch_ != "LR",
      std::invalid_argument, "Anasazi::RTRSolMgr: Invalid sorting string. RTR solvers compute only LR or SR.");

  // convergence tolerance
  convtol_ = pl_.get("Convergence Tolerance",convtol_);
  relconvtol_ = pl_.get("Relative Convergence Tolerance",relconvtol_);
  strtmp = pl_.get("Convergence Norm",std::string("2"));
  if (strtmp == "2") {
    convNorm_ = RES_2NORM;
  }
  else if (strtmp == "M") {
    convNorm_ = RES_ORTH;
  }
  else {
    TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
        "Anasazi::RTRSolMgr: Invalid Convergence Norm.");
  }


  // maximum number of (outer) iterations
  maxIters_ = pl_.get("Maximum Iterations",maxIters_);

  // skinny solver or not
  skinny_ = pl_.get("Skinny Solver",skinny_);

  // number if ICGS iterations
  numICGS_ = pl_.get("Num ICGS",2);

  // output stream
  std::string fntemplate = "";
  bool allProcs = false;
  if (pl_.isParameter("Output on all processors")) {
    if (Teuchos::isParameterType<bool>(pl_,"Output on all processors")) {
      allProcs = pl_.get("Output on all processors",allProcs);
    } else {
      allProcs = ( Teuchos::getParameter<int>(pl_,"Output on all processors") != 0 );
    }
  }
  fntemplate = pl_.get("Output filename template",fntemplate);
  int MyPID;
# ifdef HAVE_MPI
    // Initialize MPI
    int mpiStarted = 0;
    MPI_Initialized(&mpiStarted);
    if (mpiStarted) MPI_Comm_rank(MPI_COMM_WORLD, &MyPID);
    else MyPID=0;
# else
    MyPID = 0;
# endif
  if (fntemplate != "") {
    std::ostringstream MyPIDstr;
    MyPIDstr << MyPID;
    // replace %d in fntemplate with MyPID
    int pos, start=0;
    while ( (pos = fntemplate.find("%d",start)) != -1 ) {
      fntemplate.replace(pos,2,MyPIDstr.str());
      start = pos+2;
    }
  }
  Teuchos::RCP<ostream> osp;
  if (fntemplate != "") {
    osp = Teuchos::rcp( new std::ofstream(fntemplate.c_str(),std::ios::out | std::ios::app) );
    if (!*osp) {
      osp = Teuchos::rcpFromRef(std::cout);
      std::cout << "Anasazi::RTRSolMgr::constructor(): Could not open file for write: " << fntemplate << std::endl;
    }
  }
  else {
    osp = Teuchos::rcpFromRef(std::cout);
  }
  // Output manager
  int verbosity = Anasazi::Errors;
  if (pl_.isParameter("Verbosity")) {
    if (Teuchos::isParameterType<int>(pl_,"Verbosity")) {
      verbosity = pl_.get("Verbosity", verbosity);
    } else {
      verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl_,"Verbosity");
    }
  }
  if (allProcs) {
    // print on all procs
    printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,MyPID) );
  }
  else {
    // print only on proc 0
    printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,0) );
  }
}


// solve()
template<class ScalarType, class MV, class OP>
ReturnType
RTRSolMgr<ScalarType,MV,OP>::solve() {

  using std::endl;

  // typedef SolverUtils<ScalarType,MV,OP> msutils; // unused

  const int nev = problem_->getNEV();

  // clear num iters
  numIters_ = -1;

#ifdef TEUCHOS_DEBUG
    Teuchos::RCP<Teuchos::FancyOStream>
      out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(Debug)));
    out->setShowAllFrontMatter(false).setShowProcRank(true);
    *out << "Entering Anasazi::RTRSolMgr::solve()\n";
#endif

  // Sort manager
  Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );

  // Status tests
  //
  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > maxtest;
  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > ordertest;
  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > combotest;
  Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
  // maximum iters
  if (maxIters_ > 0) {
    maxtest = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
  }
  else {
    maxtest = Teuchos::null;
  }
  //
  // convergence
  convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
  ordertest = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
  //
  // combo
  Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
  alltests.push_back(ordertest);
  if (maxtest != Teuchos::null) alltests.push_back(maxtest);
  combotest = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, alltests) );
  //
  // printing StatusTest
  Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
  if ( printer_->isVerbosity(Debug) ) {
    outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed+Failed+Undefined ) );
  }
  else {
    outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed ) );
  }

  // Orthomanager
  Teuchos::RCP<ICGSOrthoManager<ScalarType,MV,OP> > ortho
    = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM(),numICGS_) );

  // create an RTR solver
  // leftmost or rightmost?
  bool leftMost = true;
  if (whch_ == "LR" || whch_ == "LM") {
    leftMost = false;
  }
  pl_.set<bool>("Leftmost",leftMost);
  Teuchos::RCP<RTRBase<ScalarType,MV,OP> > rtr_solver;
  if (skinny_ == false) {
    // "hefty" IRTR
    rtr_solver = Teuchos::rcp( new  IRTR<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,pl_) );
  }
  else {
    // "skinny" IRTR
    rtr_solver = Teuchos::rcp( new SIRTR<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,pl_) );
  }
  // set any auxiliary vectors defined in the problem
  Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
  if (probauxvecs != Teuchos::null) {
    rtr_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
  }

  TEUCHOS_TEST_FOR_EXCEPTION(rtr_solver->getBlockSize() < problem_->getNEV(),std::logic_error,
            "Anasazi::RTRSolMgr requires block size >= requested number of eigenvalues.");

  int numfound = 0;
  Teuchos::RCP<MV> foundvecs;
  std::vector<MagnitudeType> foundvals;

  // tell the solver to iterate
  try {
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
    Teuchos::TimeMonitor slvtimer(*_timerSolve);
#endif
    rtr_solver->iterate();
    numIters_ = rtr_solver->getNumIters();
  }
  catch (const std::exception &e) {
    // we are a simple solver manager. we don't catch exceptions. set solution empty, then rethrow.
    printer_->stream(Anasazi::Errors) << "Exception: " << e.what() << endl;
    Eigensolution<ScalarType,MV> sol;
    sol.numVecs = 0;
    problem_->setSolution(sol);
    throw;
  }

  // check the status tests
  if (convtest->getStatus() == Passed || (maxtest != Teuchos::null && maxtest->getStatus() == Passed))
  {
    int num = convtest->howMany();
    if (num > 0) {
      std::vector<int> ind = convtest->whichVecs();
      // copy the converged eigenvectors
      foundvecs = MVT::CloneCopy(*rtr_solver->getRitzVectors(),ind);
      // copy the converged eigenvalues
      foundvals.resize(num);
      std::vector<Value<ScalarType> > all = rtr_solver->getRitzValues();
      for (int i=0; i<num; i++) {
        foundvals[i] = all[ind[i]].realpart;
      }
      numfound = num;
    }
  }
  else {
    TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::RTRSolMgr::solve(): solver returned without satisfy status test.");
  }

  // create contiguous storage for all eigenvectors, eigenvalues
  Eigensolution<ScalarType,MV> sol;
  sol.numVecs = numfound;
  sol.Evecs = foundvecs;
  sol.Espace = sol.Evecs;
  sol.Evals.resize(sol.numVecs);
  for (int i=0; i<sol.numVecs; i++) {
    sol.Evals[i].realpart = foundvals[i];
  }
  // all real eigenvalues: set index vectors [0,...,numfound-1]
  sol.index.resize(numfound,0);

  // print final summary
  rtr_solver->currentStatus(printer_->stream(FinalSummary));

  // print timing information
#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
  if ( printer_->isVerbosity( TimingDetails ) ) {
    Teuchos::TimeMonitor::summarize( printer_->stream( TimingDetails ) );
  }
#endif

  // send the solution to the eigenproblem
  problem_->setSolution(sol);
  printer_->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << endl;

  // return from SolMgr::solve()
  if (sol.numVecs < nev) return Unconverged;
  return Converged;
}




} // end Anasazi namespace

#endif /* ANASAZI_RTR_SOLMGR_HPP */