Anasazi Version of the Day
AnasaziRTRSolMgr.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef ANASAZI_RTR_SOLMGR_HPP
00030 #define ANASAZI_RTR_SOLMGR_HPP
00031 
00037 #include "AnasaziConfigDefs.hpp"
00038 #include "AnasaziTypes.hpp"
00039 
00040 #include "AnasaziEigenproblem.hpp"
00041 #include "AnasaziSolverManager.hpp"
00042 #include "AnasaziSolverUtils.hpp"
00043 
00044 #include "AnasaziIRTR.hpp"
00045 #include "AnasaziSIRTR.hpp"
00046 #include "AnasaziBasicSort.hpp"
00047 #include "AnasaziICGSOrthoManager.hpp"
00048 #include "AnasaziStatusTestMaxIters.hpp"
00049 #include "AnasaziStatusTestResNorm.hpp"
00050 #include "AnasaziStatusTestWithOrdering.hpp"
00051 #include "AnasaziStatusTestCombo.hpp"
00052 #include "AnasaziStatusTestOutput.hpp"
00053 #include "AnasaziBasicOutputManager.hpp"
00054 
00055 #include <Teuchos_TimeMonitor.hpp>
00056 #include <Teuchos_FancyOStream.hpp>
00057 
00067 namespace Anasazi {
00068 
00069 template<class ScalarType, class MV, class OP>
00070 class RTRSolMgr : public SolverManager<ScalarType,MV,OP> {
00071 
00072   private:
00073     typedef MultiVecTraits<ScalarType,MV> MVT;
00074     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00075     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00076     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00077     typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00078 
00079   public:
00080 
00082 
00083 
00099   RTRSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00100              Teuchos::ParameterList &pl );
00101 
00103   virtual ~RTRSolMgr() {};
00105 
00107 
00108 
00110   const Eigenproblem<ScalarType,MV,OP>& getProblem() const {
00111     return *problem_;
00112   }
00113 
00119    Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00120      return Teuchos::tuple(_timerSolve);
00121    }
00122 
00124   int getNumIters() const {
00125     return numIters_;
00126   }
00127 
00128 
00130 
00132 
00133 
00142   ReturnType solve();
00144 
00145   private:
00146   Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
00147   std::string whch_;
00148   std::string ortho_;
00149 
00150   bool skinny_;
00151   MagnitudeType convtol_;
00152   int maxIters_;
00153   bool relconvtol_;
00154   enum StatusTestResNorm<ScalarType,MV,OP>::ResType convNorm_;
00155   int numIters_;
00156   int numICGS_;
00157 
00158   Teuchos::RCP<Teuchos::Time> _timerSolve;
00159   Teuchos::RCP<BasicOutputManager<ScalarType> > printer_;
00160   Teuchos::ParameterList pl_;
00161 };
00162 
00163 
00165 template<class ScalarType, class MV, class OP>
00166 RTRSolMgr<ScalarType,MV,OP>::RTRSolMgr(
00167         const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00168         Teuchos::ParameterList &pl ) :
00169   problem_(problem),
00170   whch_("SR"),
00171   skinny_(true),
00172   convtol_(MT::prec()),
00173   maxIters_(100),
00174   relconvtol_(true),
00175   numIters_(-1),
00176 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00177   _timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: RTRSolMgr::solve()")),
00178 #endif
00179   pl_(pl)
00180 {
00181   TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,              std::invalid_argument, "Problem not given to solver manager.");
00182   TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),              std::invalid_argument, "Problem not set.");
00183   TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(),               std::invalid_argument, "Problem not symmetric.");
00184   TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
00185 
00186   std::string strtmp;
00187 
00188   whch_ = pl_.get("Which","SR");
00189   TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SR" && whch_ != "LR",
00190       std::invalid_argument, "Anasazi::RTRSolMgr: Invalid sorting string. RTR solvers compute only LR or SR.");
00191 
00192   // convergence tolerance
00193   convtol_ = pl_.get("Convergence Tolerance",convtol_);
00194   relconvtol_ = pl_.get("Relative Convergence Tolerance",relconvtol_);
00195   strtmp = pl_.get("Convergence Norm",std::string("2"));
00196   if (strtmp == "2") {
00197     convNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_2NORM;
00198   }
00199   else if (strtmp == "M") {
00200     convNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_ORTH;
00201   }
00202   else {
00203     TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
00204         "Anasazi::RTRSolMgr: Invalid Convergence Norm.");
00205   }
00206 
00207 
00208   // maximum number of (outer) iterations
00209   maxIters_ = pl_.get("Maximum Iterations",maxIters_);
00210 
00211   // skinny solver or not
00212   skinny_ = pl_.get("Skinny Solver",skinny_);
00213 
00214   // number if ICGS iterations
00215   numICGS_ = pl_.get("Num ICGS",2);
00216 
00217   // output stream
00218   std::string fntemplate = "";
00219   bool allProcs = false;
00220   if (pl_.isParameter("Output on all processors")) {
00221     if (Teuchos::isParameterType<bool>(pl_,"Output on all processors")) {
00222       allProcs = pl_.get("Output on all processors",allProcs);
00223     } else {
00224       allProcs = ( Teuchos::getParameter<int>(pl_,"Output on all processors") != 0 );
00225     }
00226   }
00227   fntemplate = pl_.get("Output filename template",fntemplate);
00228   int MyPID;
00229 # ifdef HAVE_MPI
00230     // Initialize MPI
00231     int mpiStarted = 0;
00232     MPI_Initialized(&mpiStarted);
00233     if (mpiStarted) MPI_Comm_rank(MPI_COMM_WORLD, &MyPID);
00234     else MyPID=0;
00235 # else
00236     MyPID = 0;
00237 # endif
00238   if (fntemplate != "") {
00239     std::ostringstream MyPIDstr;
00240     MyPIDstr << MyPID;
00241     // replace %d in fntemplate with MyPID
00242     int pos, start=0;
00243     while ( (pos = fntemplate.find("%d",start)) != -1 ) {
00244       fntemplate.replace(pos,2,MyPIDstr.str());
00245       start = pos+2;
00246     }
00247   }
00248   Teuchos::RCP<ostream> osp;
00249   if (fntemplate != "") {
00250     osp = Teuchos::rcp( new std::ofstream(fntemplate.c_str(),std::ios::out | std::ios::app) );
00251     if (!*osp) {
00252       osp = Teuchos::rcpFromRef(std::cout);
00253       std::cout << "Anasazi::RTRSolMgr::constructor(): Could not open file for write: " << fntemplate << std::endl;
00254     }
00255   }
00256   else {
00257     osp = Teuchos::rcpFromRef(std::cout);
00258   }
00259   // Output manager
00260   int verbosity = Anasazi::Errors;
00261   if (pl_.isParameter("Verbosity")) {
00262     if (Teuchos::isParameterType<int>(pl_,"Verbosity")) {
00263       verbosity = pl_.get("Verbosity", verbosity);
00264     } else {
00265       verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl_,"Verbosity");
00266     }
00267   }
00268   if (allProcs) {
00269     // print on all procs
00270     printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,MyPID) );
00271   }
00272   else {
00273     // print only on proc 0
00274     printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,0) );
00275   }
00276 }
00277 
00278 
00279 // solve()
00280 template<class ScalarType, class MV, class OP>
00281 ReturnType
00282 RTRSolMgr<ScalarType,MV,OP>::solve() {
00283 
00284   using std::endl;
00285 
00286   // typedef SolverUtils<ScalarType,MV,OP> msutils; // unused
00287 
00288   const int nev = problem_->getNEV();
00289 
00290   // clear num iters
00291   numIters_ = -1;
00292 
00293 #ifdef TEUCHOS_DEBUG
00294     Teuchos::RCP<Teuchos::FancyOStream>
00295       out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(Debug)));
00296     out->setShowAllFrontMatter(false).setShowProcRank(true);
00297     *out << "Entering Anasazi::RTRSolMgr::solve()\n";
00298 #endif
00299 
00301   // Sort manager
00302   Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
00303 
00305   // Status tests
00306   //
00307   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > maxtest;
00308   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > ordertest;
00309   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > combotest;
00310   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
00311   // maximum iters
00312   if (maxIters_ > 0) {
00313     maxtest = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
00314   }
00315   else {
00316     maxtest = Teuchos::null;
00317   }
00318   //
00319   // convergence
00320   convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
00321   ordertest = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
00322   //
00323   // combo
00324   Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
00325   alltests.push_back(ordertest);
00326   if (maxtest != Teuchos::null) alltests.push_back(maxtest);
00327   combotest = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, alltests) );
00328   //
00329   // printing StatusTest
00330   Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
00331   if ( printer_->isVerbosity(Debug) ) {
00332     outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed+Failed+Undefined ) );
00333   }
00334   else {
00335     outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed ) );
00336   }
00337 
00339   // Orthomanager
00340   Teuchos::RCP<ICGSOrthoManager<ScalarType,MV,OP> > ortho
00341     = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM(),numICGS_) );
00342 
00344   // create an RTR solver
00345   // leftmost or rightmost?
00346   bool leftMost = true;
00347   if (whch_ == "LR" || whch_ == "LM") {
00348     leftMost = false;
00349   }
00350   pl_.set<bool>("Leftmost",leftMost);
00351   Teuchos::RCP<RTRBase<ScalarType,MV,OP> > rtr_solver;
00352   if (skinny_ == false) {
00353     // "hefty" IRTR
00354     rtr_solver = Teuchos::rcp( new  IRTR<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,pl_) );
00355   }
00356   else {
00357     // "skinny" IRTR
00358     rtr_solver = Teuchos::rcp( new SIRTR<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,pl_) );
00359   }
00360   // set any auxiliary vectors defined in the problem
00361   Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
00362   if (probauxvecs != Teuchos::null) {
00363     rtr_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
00364   }
00365 
00366   TEUCHOS_TEST_FOR_EXCEPTION(rtr_solver->getBlockSize() < problem_->getNEV(),std::logic_error,
00367             "Anasazi::RTRSolMgr requires block size >= requested number of eigenvalues.");
00368 
00369   int numfound = 0;
00370   Teuchos::RCP<MV> foundvecs;
00371   std::vector<MagnitudeType> foundvals;
00372 
00373   // tell the solver to iterate
00374   try {
00375 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00376     Teuchos::TimeMonitor slvtimer(*_timerSolve);
00377 #endif
00378     rtr_solver->iterate();
00379     numIters_ = rtr_solver->getNumIters();
00380   }
00381   catch (const std::exception &e) {
00382     // we are a simple solver manager. we don't catch exceptions. set solution empty, then rethrow.
00383     printer_->stream(Anasazi::Errors) << "Exception: " << e.what() << endl;
00384     Eigensolution<ScalarType,MV> sol;
00385     sol.numVecs = 0;
00386     problem_->setSolution(sol);
00387     throw;
00388   }
00389 
00390   // check the status tests
00391   if (convtest->getStatus() == Passed || (maxtest != Teuchos::null && maxtest->getStatus() == Passed))
00392   {
00393     int num = convtest->howMany();
00394     if (num > 0) {
00395       std::vector<int> ind = convtest->whichVecs();
00396       // copy the converged eigenvectors
00397       foundvecs = MVT::CloneCopy(*rtr_solver->getRitzVectors(),ind);
00398       // copy the converged eigenvalues
00399       foundvals.resize(num);
00400       std::vector<Value<ScalarType> > all = rtr_solver->getRitzValues();
00401       for (int i=0; i<num; i++) {
00402         foundvals[i] = all[ind[i]].realpart;
00403       }
00404       numfound = num;
00405     }
00406   }
00407   else {
00408     TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::RTRSolMgr::solve(): solver returned without satisfy status test.");
00409   }
00410 
00411   // create contiguous storage for all eigenvectors, eigenvalues
00412   Eigensolution<ScalarType,MV> sol;
00413   sol.numVecs = numfound;
00414   sol.Evecs = foundvecs;
00415   sol.Espace = sol.Evecs;
00416   sol.Evals.resize(sol.numVecs);
00417   for (int i=0; i<sol.numVecs; i++) {
00418     sol.Evals[i].realpart = foundvals[i];
00419   }
00420   // all real eigenvalues: set index vectors [0,...,numfound-1]
00421   sol.index.resize(numfound,0);
00422 
00423   // print final summary
00424   rtr_solver->currentStatus(printer_->stream(FinalSummary));
00425 
00426   // print timing information
00427 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
00428   if ( printer_->isVerbosity( TimingDetails ) ) {
00429     Teuchos::TimeMonitor::summarize( printer_->stream( TimingDetails ) );
00430   }
00431 #endif
00432 
00433   // send the solution to the eigenproblem
00434   problem_->setSolution(sol);
00435   printer_->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << endl;
00436 
00437   // return from SolMgr::solve()
00438   if (sol.numVecs < nev) return Unconverged;
00439   return Converged;
00440 }
00441 
00442 
00443 
00444 
00445 } // end Anasazi namespace
00446 
00447 #endif /* ANASAZI_RTR_SOLMGR_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends