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

Generated on Wed May 12 21:40:23 2010 for Anasazi by  doxygen 1.4.7