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 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   _timerSolve(Teuchos::TimeMonitor::getNewTimer("RTRSolMgr::solve()")),
00177   pl_(pl)
00178 {
00179   TEST_FOR_EXCEPTION(problem_ == Teuchos::null,              std::invalid_argument, "Problem not given to solver manager.");
00180   TEST_FOR_EXCEPTION(!problem_->isProblemSet(),              std::invalid_argument, "Problem not set.");
00181   TEST_FOR_EXCEPTION(!problem_->isHermitian(),               std::invalid_argument, "Problem not symmetric.");
00182   TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
00183 
00184   std::string strtmp;
00185 
00186   whch_ = pl_.get("Which","SR");
00187   TEST_FOR_EXCEPTION(whch_ != "SR" && whch_ != "LR",
00188       std::invalid_argument, "Anasazi::RTRSolMgr: Invalid sorting string. RTR solvers compute only LR or SR.");
00189 
00190   // convergence tolerance
00191   convtol_ = pl_.get("Convergence Tolerance",convtol_);
00192   relconvtol_ = pl_.get("Relative Convergence Tolerance",relconvtol_);
00193   strtmp = pl_.get("Convergence Norm",std::string("2"));
00194   if (strtmp == "2") {
00195     convNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_2NORM;
00196   }
00197   else if (strtmp == "M") {
00198     convNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_ORTH;
00199   }
00200   else {
00201     TEST_FOR_EXCEPTION(true, std::invalid_argument, 
00202         "Anasazi::RTRSolMgr: Invalid Convergence Norm.");
00203   }
00204   
00205 
00206   // maximum number of (outer) iterations
00207   maxIters_ = pl_.get("Maximum Iterations",maxIters_);
00208 
00209   // skinny solver or not
00210   skinny_ = pl_.get("Skinny Solver",skinny_);
00211 
00212   // number if ICGS iterations
00213   numICGS_ = pl_.get("Num ICGS",2);
00214 
00215   // output stream
00216   std::string fntemplate = "";
00217   bool allProcs = false;
00218   if (pl_.isParameter("Output on all processors")) {
00219     if (Teuchos::isParameterType<bool>(pl_,"Output on all processors")) {
00220       allProcs = pl_.get("Output on all processors",allProcs);
00221     } else {
00222       allProcs = ( Teuchos::getParameter<int>(pl_,"Output on all processors") != 0 );
00223     }
00224   }
00225   fntemplate = pl_.get("Output filename template",fntemplate);
00226   int MyPID;
00227 # ifdef HAVE_MPI
00228     // Initialize MPI
00229     int mpiStarted = 0;
00230     MPI_Initialized(&mpiStarted);
00231     if (mpiStarted) MPI_Comm_rank(MPI_COMM_WORLD, &MyPID);
00232     else MyPID=0;
00233 # else 
00234     MyPID = 0;
00235 # endif
00236   if (fntemplate != "") {
00237     std::ostringstream MyPIDstr;
00238     MyPIDstr << MyPID;
00239     // replace %d in fntemplate with MyPID
00240     int pos, start=0;
00241     while ( (pos = fntemplate.find("%d",start)) != -1 ) {
00242       fntemplate.replace(pos,2,MyPIDstr.str());
00243       start = pos+2;
00244     }
00245   }
00246   Teuchos::RCP<ostream> osp;
00247   if (fntemplate != "") {
00248     osp = Teuchos::rcp( new std::ofstream(fntemplate.c_str(),std::ios::out | std::ios::app) );
00249     if (!*osp) {
00250       osp = Teuchos::rcpFromRef(std::cout);
00251       std::cout << "Anasazi::RTRSolMgr::constructor(): Could not open file for write: " << fntemplate << std::endl;
00252     }
00253   }
00254   else {
00255     osp = Teuchos::rcpFromRef(std::cout);
00256   }
00257   // Output manager
00258   int verbosity = Anasazi::Errors;
00259   if (pl_.isParameter("Verbosity")) {
00260     if (Teuchos::isParameterType<int>(pl_,"Verbosity")) {
00261       verbosity = pl_.get("Verbosity", verbosity);
00262     } else {
00263       verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl_,"Verbosity");
00264     }
00265   }
00266   if (allProcs) {
00267     // print on all procs
00268     printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,MyPID) );
00269   }
00270   else {
00271     // print only on proc 0
00272     printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,0) );
00273   }
00274 }
00275 
00276 
00277 // solve()
00278 template<class ScalarType, class MV, class OP>
00279 ReturnType 
00280 RTRSolMgr<ScalarType,MV,OP>::solve() {
00281 
00282   using std::endl;
00283 
00284   typedef SolverUtils<ScalarType,MV,OP> msutils;
00285 
00286   const int nev = problem_->getNEV();
00287 
00288   // clear num iters
00289   numIters_ = -1;
00290 
00291 #ifdef TEUCHOS_DEBUG
00292     Teuchos::RCP<Teuchos::FancyOStream>
00293       out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(Debug)));
00294     out->setShowAllFrontMatter(false).setShowProcRank(true);
00295     *out << "Entering Anasazi::RTRSolMgr::solve()\n";
00296 #endif
00297 
00299   // Sort manager
00300   Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
00301 
00303   // Status tests
00304   //
00305   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > maxtest;
00306   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > ordertest;
00307   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > combotest;
00308   Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
00309   // maximum iters
00310   if (maxIters_ > 0) {
00311     maxtest = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
00312   }
00313   else {
00314     maxtest = Teuchos::null;
00315   }
00316   //
00317   // convergence
00318   convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
00319   ordertest = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
00320   //
00321   // combo 
00322   Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
00323   alltests.push_back(ordertest);
00324   if (maxtest != Teuchos::null) alltests.push_back(maxtest);
00325   combotest = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, alltests) );
00326   //
00327   // printing StatusTest
00328   Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
00329   if ( printer_->isVerbosity(Debug) ) {
00330     outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed+Failed+Undefined ) );
00331   }
00332   else {
00333     outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed ) );
00334   }
00335 
00337   // Orthomanager
00338   Teuchos::RCP<ICGSOrthoManager<ScalarType,MV,OP> > ortho 
00339     = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM(),numICGS_) );
00340 
00342   // create an RTR solver
00343   // leftmost or rightmost?
00344   bool leftMost = true;
00345   if (whch_ == "LR" || whch_ == "LM") {
00346     leftMost = false;
00347   }
00348   pl_.set<bool>("Leftmost",leftMost);
00349   Teuchos::RCP<RTRBase<ScalarType,MV,OP> > rtr_solver;
00350   if (skinny_ == false) {
00351     // "hefty" IRTR
00352     rtr_solver = Teuchos::rcp( new  IRTR<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,pl_) );
00353   }
00354   else {
00355     // "skinny" IRTR
00356     rtr_solver = Teuchos::rcp( new SIRTR<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,pl_) );
00357   }
00358   // set any auxiliary vectors defined in the problem
00359   Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
00360   if (probauxvecs != Teuchos::null) {
00361     rtr_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
00362   }
00363 
00364   TEST_FOR_EXCEPTION(rtr_solver->getBlockSize() < problem_->getNEV(),std::logic_error,
00365             "Anasazi::RTRSolMgr requires block size >= requested number of eigenvalues.");
00366 
00367   int numfound = 0;
00368   Teuchos::RCP<MV> foundvecs;
00369   std::vector<MagnitudeType> foundvals;
00370 
00371   // tell the solver to iterate
00372   try {
00373     Teuchos::TimeMonitor slvtimer(*_timerSolve);
00374     rtr_solver->iterate();
00375     numIters_ = rtr_solver->getNumIters();
00376   }
00377   catch (const std::exception &e) {
00378     // we are a simple solver manager. we don't catch exceptions. set solution empty, then rethrow.
00379     printer_->stream(Anasazi::Errors) << "Exception: " << e.what() << endl;
00380     Eigensolution<ScalarType,MV> sol;
00381     sol.numVecs = 0;
00382     problem_->setSolution(sol);
00383     throw;
00384   }
00385 
00386   // check the status tests
00387   if (convtest->getStatus() == Passed || (maxtest != Teuchos::null && maxtest->getStatus() == Passed)) 
00388   {
00389     int num = convtest->howMany();
00390     if (num > 0) {
00391       std::vector<int> ind = convtest->whichVecs();
00392       // copy the converged eigenvectors
00393       foundvecs = MVT::CloneCopy(*rtr_solver->getRitzVectors(),ind);
00394       // copy the converged eigenvalues
00395       foundvals.resize(num);
00396       std::vector<Value<ScalarType> > all = rtr_solver->getRitzValues();
00397       for (int i=0; i<num; i++) {
00398         foundvals[i] = all[ind[i]].realpart;
00399       }
00400       numfound = num;
00401     }
00402   }
00403   else {
00404     TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::RTRSolMgr::solve(): solver returned without satisfy status test.");
00405   }
00406 
00407   // create contiguous storage for all eigenvectors, eigenvalues
00408   Eigensolution<ScalarType,MV> sol;
00409   sol.numVecs = numfound;
00410   sol.Evecs = foundvecs;
00411   sol.Espace = sol.Evecs;
00412   sol.Evals.resize(sol.numVecs);
00413   for (int i=0; i<sol.numVecs; i++) {
00414     sol.Evals[i].realpart = foundvals[i];
00415   }
00416   // all real eigenvalues: set index vectors [0,...,numfound-1]
00417   sol.index.resize(numfound,0);
00418 
00419   // print final summary
00420   rtr_solver->currentStatus(printer_->stream(FinalSummary));
00421 
00422   // print timing information
00423   Teuchos::TimeMonitor::summarize(printer_->stream(TimingDetails));
00424 
00425   // send the solution to the eigenproblem
00426   problem_->setSolution(sol);
00427   printer_->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << endl;
00428 
00429   // return from SolMgr::solve()
00430   if (sol.numVecs < nev) return Unconverged;
00431   return Converged;
00432 }
00433 
00434 
00435 
00436 
00437 } // end Anasazi namespace
00438 
00439 #endif /* ANASAZI_RTR_SOLMGR_HPP */

Generated on Wed May 12 21:24:34 2010 for Anasazi by  doxygen 1.4.7