00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
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
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
00206 maxIters_ = pl_.get("Maximum Iterations",maxIters_);
00207
00208
00209 skinny_ = pl_.get("Skinny Solver",skinny_);
00210
00211
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
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
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
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
00264 printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,MyPID) );
00265 }
00266 else {
00267
00268 printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,0) );
00269 }
00270 }
00271
00272
00273
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
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
00296 Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
00297
00299
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
00306 if (maxIters_ > 0) {
00307 maxtest = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>(maxIters_) );
00308 }
00309 else {
00310 maxtest = Teuchos::null;
00311 }
00312
00313
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
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
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
00334 Teuchos::RCP<ICGSOrthoManager<ScalarType,MV,OP> > ortho
00335 = Teuchos::rcp( new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
00336
00338
00339
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
00348 rtr_solver = Teuchos::rcp( new IRTR<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,pl_) );
00349 }
00350 else {
00351
00352 rtr_solver = Teuchos::rcp( new SIRTR<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,pl_) );
00353 }
00354
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
00368 try {
00369 Teuchos::TimeMonitor slvtimer(*_timerSolve);
00370 rtr_solver->iterate();
00371 numIters_ = rtr_solver->getNumIters();
00372 }
00373 catch (std::exception e) {
00374
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
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
00389 foundvecs = MVT::CloneCopy(*rtr_solver->getRitzVectors(),ind);
00390
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
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
00413 sol.index.resize(numfound,0);
00414
00415
00416 rtr_solver->currentStatus(printer_->stream(FinalSummary));
00417
00418
00419 Teuchos::TimeMonitor::summarize(printer_->stream(TimingDetails));
00420
00421
00422 problem_->setSolution(sol);
00423 printer_->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << endl;
00424
00425
00426 if (sol.numVecs < nev) return Unconverged;
00427 return Converged;
00428 }
00429
00430
00431
00432
00433 }
00434
00435 #endif