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 BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP
00030 #define BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP
00031
00036 #include "BelosConfigDefs.hpp"
00037 #include "BelosTypes.hpp"
00038
00039 #include "BelosLinearProblem.hpp"
00040 #include "BelosSolverManager.hpp"
00041
00042 #include "BelosPseudoBlockCGIter.hpp"
00043 #include "BelosStatusTestMaxIters.hpp"
00044 #include "BelosStatusTestGenResNorm.hpp"
00045 #include "BelosStatusTestCombo.hpp"
00046 #include "BelosStatusTestOutput.hpp"
00047 #include "BelosOutputManager.hpp"
00048 #include "Teuchos_BLAS.hpp"
00049 #include "Teuchos_TimeMonitor.hpp"
00050
00066 namespace Belos {
00067
00069
00070
00077 class PseudoBlockCGSolMgrLinearProblemFailure : public BelosError {public:
00078 PseudoBlockCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
00079 {}};
00080
00087 class PseudoBlockCGSolMgrOrthoFailure : public BelosError {public:
00088 PseudoBlockCGSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
00089 {}};
00090
00091 template<class ScalarType, class MV, class OP>
00092 class PseudoBlockCGSolMgr : public SolverManager<ScalarType,MV,OP> {
00093
00094 private:
00095 typedef MultiVecTraits<ScalarType,MV> MVT;
00096 typedef OperatorTraits<ScalarType,MV,OP> OPT;
00097 typedef Teuchos::ScalarTraits<ScalarType> SCT;
00098 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00099 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00100
00101 public:
00102
00104
00105
00111 PseudoBlockCGSolMgr();
00112
00121 PseudoBlockCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00122 const Teuchos::RCP<Teuchos::ParameterList> &pl );
00123
00125 virtual ~PseudoBlockCGSolMgr() {};
00127
00129
00130
00131 const LinearProblem<ScalarType,MV,OP>& getProblem() const {
00132 return *problem_;
00133 }
00134
00137 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
00138
00141 Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
00142
00148 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00149 return tuple(timerSolve_);
00150 }
00151
00153 int getNumIters() const {
00154 return numIters_;
00155 }
00156
00160 bool isLOADetected() const { return false; }
00161
00163
00165
00166
00168 void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) { problem_ = problem; }
00169
00171 void setParameters( const Teuchos::RCP<Teuchos::ParameterList> ¶ms );
00172
00174
00176
00177
00181 void reset( const ResetType type ) { if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem(); }
00183
00185
00186
00204 ReturnType solve();
00205
00207
00210
00212 std::string description() const;
00213
00215
00216 private:
00217
00218
00219 Belos::ScaleType convertStringToScaleType( const std::string& scaleType ) {
00220 if (scaleType == "Norm of Initial Residual") {
00221 return Belos::NormOfInitRes;
00222 } else if (scaleType == "Norm of Preconditioned Initial Residual") {
00223 return Belos::NormOfPrecInitRes;
00224 } else if (scaleType == "Norm of RHS") {
00225 return Belos::NormOfRHS;
00226 } else if (scaleType == "None") {
00227 return Belos::None;
00228 } else
00229 TEST_FOR_EXCEPTION( true ,std::logic_error,
00230 "Belos::PseudoBlockCGSolMgr(): Invalid residual scaling type.");
00231 }
00232
00233
00234 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
00235
00236
00237 Teuchos::RCP<OutputManager<ScalarType> > printer_;
00238 Teuchos::RCP<std::ostream> outputStream_;
00239
00240
00241 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
00242 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
00243 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
00244 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
00245
00246
00247 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
00248
00249
00250 Teuchos::RCP<ParameterList> params_;
00251
00252
00253 static const MagnitudeType convtol_default_;
00254 static const int maxIters_default_;
00255 static const bool showMaxResNormOnly_default_;
00256 static const int verbosity_default_;
00257 static const int outputFreq_default_;
00258 static const int defQuorum_default_;
00259 static const std::string resScale_default_;
00260 static const std::string label_default_;
00261 static const Teuchos::RCP<std::ostream> outputStream_default_;
00262
00263
00264 MagnitudeType convtol_;
00265 int maxIters_, numIters_;
00266 int verbosity_, outputFreq_, defQuorum_;
00267 bool showMaxResNormOnly_;
00268 std::string resScale_;
00269
00270
00271 std::string label_;
00272 Teuchos::RCP<Teuchos::Time> timerSolve_;
00273
00274
00275 bool isSet_, isSTSet_;
00276 };
00277
00278
00279
00280 template<class ScalarType, class MV, class OP>
00281 const typename Teuchos::ScalarTraits<ScalarType>::magnitudeType PseudoBlockCGSolMgr<ScalarType,MV,OP>::convtol_default_ = 1e-8;
00282
00283 template<class ScalarType, class MV, class OP>
00284 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::maxIters_default_ = 1000;
00285
00286 template<class ScalarType, class MV, class OP>
00287 const bool PseudoBlockCGSolMgr<ScalarType,MV,OP>::showMaxResNormOnly_default_ = false;
00288
00289 template<class ScalarType, class MV, class OP>
00290 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::verbosity_default_ = Belos::Errors;
00291
00292 template<class ScalarType, class MV, class OP>
00293 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::outputFreq_default_ = -1;
00294
00295 template<class ScalarType, class MV, class OP>
00296 const int PseudoBlockCGSolMgr<ScalarType,MV,OP>::defQuorum_default_ = 1;
00297
00298 template<class ScalarType, class MV, class OP>
00299 const std::string PseudoBlockCGSolMgr<ScalarType,MV,OP>::resScale_default_ = "Norm of Initial Residual";
00300
00301 template<class ScalarType, class MV, class OP>
00302 const std::string PseudoBlockCGSolMgr<ScalarType,MV,OP>::label_default_ = "Belos";
00303
00304 template<class ScalarType, class MV, class OP>
00305 const Teuchos::RCP<std::ostream> PseudoBlockCGSolMgr<ScalarType,MV,OP>::outputStream_default_ = Teuchos::rcp(&std::cout,false);
00306
00307
00308
00309 template<class ScalarType, class MV, class OP>
00310 PseudoBlockCGSolMgr<ScalarType,MV,OP>::PseudoBlockCGSolMgr() :
00311 outputStream_(outputStream_default_),
00312 convtol_(convtol_default_),
00313 maxIters_(maxIters_default_),
00314 verbosity_(verbosity_default_),
00315 outputFreq_(outputFreq_default_),
00316 defQuorum_(defQuorum_default_),
00317 showMaxResNormOnly_(showMaxResNormOnly_default_),
00318 resScale_(resScale_default_),
00319 label_(label_default_),
00320 isSet_(false),
00321 isSTSet_(false)
00322 {}
00323
00324
00325 template<class ScalarType, class MV, class OP>
00326 PseudoBlockCGSolMgr<ScalarType,MV,OP>::PseudoBlockCGSolMgr(
00327 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00328 const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
00329 problem_(problem),
00330 outputStream_(outputStream_default_),
00331 convtol_(convtol_default_),
00332 maxIters_(maxIters_default_),
00333 verbosity_(verbosity_default_),
00334 outputFreq_(outputFreq_default_),
00335 defQuorum_(defQuorum_default_),
00336 showMaxResNormOnly_(showMaxResNormOnly_default_),
00337 resScale_(resScale_default_),
00338 label_(label_default_),
00339 isSet_(false),
00340 isSTSet_(false)
00341 {
00342 TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00343
00344
00345 if (!is_null(pl)) {
00346
00347 setParameters( pl );
00348 }
00349 }
00350
00351 template<class ScalarType, class MV, class OP>
00352 void PseudoBlockCGSolMgr<ScalarType,MV,OP>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> ¶ms )
00353 {
00354
00355 if (params_ == Teuchos::null) {
00356 params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
00357 }
00358 else {
00359 params->validateParameters(*getValidParameters());
00360 }
00361
00362
00363 if (params->isParameter("Maximum Iterations")) {
00364 maxIters_ = params->get("Maximum Iterations",maxIters_default_);
00365
00366
00367 params_->set("Maximum Iterations", maxIters_);
00368 if (maxIterTest_!=Teuchos::null)
00369 maxIterTest_->setMaxIters( maxIters_ );
00370 }
00371
00372
00373 if (params->isParameter("Timer Label")) {
00374 std::string tempLabel = params->get("Timer Label", label_default_);
00375
00376
00377 if (tempLabel != label_) {
00378 label_ = tempLabel;
00379 params_->set("Timer Label", label_);
00380 std::string solveLabel = label_ + ": PseudoBlockCGSolMgr total solve time";
00381 timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00382 }
00383 }
00384
00385
00386 if (params->isParameter("Verbosity")) {
00387 if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
00388 verbosity_ = params->get("Verbosity", verbosity_default_);
00389 } else {
00390 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
00391 }
00392
00393
00394 params_->set("Verbosity", verbosity_);
00395 if (printer_ != Teuchos::null)
00396 printer_->setVerbosity(verbosity_);
00397 }
00398
00399
00400 if (params->isParameter("Output Stream")) {
00401 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
00402
00403
00404 params_->set("Output Stream", outputStream_);
00405 if (printer_ != Teuchos::null)
00406 printer_->setOStream( outputStream_ );
00407 }
00408
00409
00410 if (verbosity_ & Belos::StatusTestDetails) {
00411 if (params->isParameter("Output Frequency")) {
00412 outputFreq_ = params->get("Output Frequency", outputFreq_default_);
00413 }
00414
00415
00416 params_->set("Output Frequency", outputFreq_);
00417 if (outputTest_ != Teuchos::null)
00418 outputTest_->setOutputFrequency( outputFreq_ );
00419 }
00420
00421
00422 if (printer_ == Teuchos::null) {
00423 printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
00424 }
00425
00426
00427 typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
00428 typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
00429
00430
00431 if (params->isParameter("Convergence Tolerance")) {
00432 convtol_ = params->get("Convergence Tolerance",convtol_default_);
00433
00434
00435 params_->set("Convergence Tolerance", convtol_);
00436 if (convTest_ != Teuchos::null)
00437 convTest_->setTolerance( convtol_ );
00438 }
00439
00440 if (params->isParameter("Show Maximum Residual Norm Only")) {
00441 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
00442
00443
00444 params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
00445 if (convTest_ != Teuchos::null)
00446 convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
00447 }
00448
00449
00450 bool newResTest = false;
00451 if (params->isParameter("Residual Scaling")) {
00452 std::string tempResScale = Teuchos::getParameter<std::string>( *params, "Residual Scaling" );
00453
00454
00455 if (resScale_ != tempResScale) {
00456 Belos::ScaleType resScaleType = convertStringToScaleType( tempResScale );
00457 resScale_ = tempResScale;
00458
00459
00460 params_->set("Residual Scaling", resScale_);
00461 if (convTest_ != Teuchos::null) {
00462 try {
00463 convTest_->defineScaleForm( resScaleType, Belos::TwoNorm );
00464 }
00465 catch (std::exception& e) {
00466
00467 newResTest = true;
00468 }
00469 }
00470 }
00471 }
00472
00473
00474 if (params->isParameter("Deflation Quorum")) {
00475 defQuorum_ = params->get("Deflation Quorum", defQuorum_);
00476 params_->set("Deflation Quorum", defQuorum_);
00477 if (convTest_ != Teuchos::null)
00478 convTest_->setQuorum( defQuorum_ );
00479 }
00480
00481
00482
00483
00484 if (maxIterTest_ == Teuchos::null)
00485 maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP>( maxIters_ ) );
00486
00487
00488 if (convTest_ == Teuchos::null || newResTest) {
00489 convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, defQuorum_, showMaxResNormOnly_ ) );
00490 convTest_->defineScaleForm( convertStringToScaleType( resScale_ ), Belos::TwoNorm );
00491 }
00492
00493 if (sTest_ == Teuchos::null || newResTest)
00494 sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
00495
00496 if (outputTest_ == Teuchos::null || newResTest) {
00497 if (outputFreq_ > 0) {
00498 outputTest_ = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,
00499 sTest_,
00500 outputFreq_,
00501 Passed+Failed+Undefined ) );
00502 }
00503 else {
00504 outputTest_ = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,
00505 sTest_, 1 ) );
00506 }
00507 }
00508
00509
00510 if (timerSolve_ == Teuchos::null) {
00511 std::string solveLabel = label_ + ": PseudoBlockCGSolMgr total solve time";
00512 timerSolve_ = Teuchos::TimeMonitor::getNewTimer(solveLabel);
00513 }
00514
00515
00516 isSet_ = true;
00517 }
00518
00519
00520 template<class ScalarType, class MV, class OP>
00521 Teuchos::RCP<const Teuchos::ParameterList>
00522 PseudoBlockCGSolMgr<ScalarType,MV,OP>::getValidParameters() const
00523 {
00524 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
00525 if (is_null(validPL)) {
00526 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
00527
00528 pl= Teuchos::rcp( new Teuchos::ParameterList() );
00529 pl->set("Convergence Tolerance", convtol_default_,
00530 "The relative residual tolerance that needs to be achieved by the\n"
00531 "iterative solver in order for the linera system to be declared converged.");
00532 pl->set("Maximum Iterations", maxIters_default_,
00533 "The maximum number of block iterations allowed for each\n"
00534 "set of RHS solved.");
00535 pl->set("Verbosity", verbosity_default_,
00536 "What type(s) of solver information should be outputted\n"
00537 "to the output stream.");
00538 pl->set("Output Frequency", outputFreq_default_,
00539 "How often convergence information should be outputted\n"
00540 "to the output stream.");
00541 pl->set("Deflation Quorum", defQuorum_default_,
00542 "The number of linear systems that need to converge before\n"
00543 "they are deflated. This number should be <= block size.");
00544 pl->set("Output Stream", outputStream_default_,
00545 "A reference-counted pointer to the output stream where all\n"
00546 "solver output is sent.");
00547 pl->set("Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
00548 "When convergence information is printed, only show the maximum\n"
00549 "relative residual norm when the block size is greater than one.");
00550 pl->set("Residual Scaling", resScale_default_,
00551 "The type of scaling used in the residual convergence test.");
00552 pl->set("Timer Label", label_default_,
00553 "The string to use as a prefix for the timer labels.");
00554
00555 validPL = pl;
00556 }
00557 return validPL;
00558 }
00559
00560
00561
00562 template<class ScalarType, class MV, class OP>
00563 ReturnType PseudoBlockCGSolMgr<ScalarType,MV,OP>::solve() {
00564
00565
00566
00567
00568 if (!isSet_) { setParameters( params_ ); }
00569
00570 Teuchos::BLAS<int,ScalarType> blas;
00571
00572 TEST_FOR_EXCEPTION(!problem_->isProblemSet(),PseudoBlockCGSolMgrLinearProblemFailure,
00573 "Belos::PseudoBlockCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
00574
00575
00576 int startPtr = 0;
00577 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
00578 int numCurrRHS = numRHS2Solve;
00579
00580 std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
00581 for (int i=0; i<numRHS2Solve; ++i) {
00582 currIdx[i] = startPtr+i;
00583 currIdx2[i]=i;
00584 }
00585
00586
00587 problem_->setLSIndex( currIdx );
00588
00590
00591 Teuchos::ParameterList plist;
00592
00593
00594 outputTest_->reset();
00595
00596
00597 bool isConverged = true;
00598
00600
00601
00602 Teuchos::RCP<PseudoBlockCGIter<ScalarType,MV,OP> > block_cg_iter
00603 = Teuchos::rcp( new PseudoBlockCGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );
00604
00605
00606 {
00607 Teuchos::TimeMonitor slvtimer(*timerSolve_);
00608
00609 while ( numRHS2Solve > 0 ) {
00610
00611
00612 std::vector<int> convRHSIdx;
00613 std::vector<int> currRHSIdx( currIdx );
00614 currRHSIdx.resize(numCurrRHS);
00615
00616
00617 block_cg_iter->resetNumIters();
00618
00619
00620 outputTest_->resetNumCalls();
00621
00622
00623 Teuchos::RCP<MV> R_0 = MVT::CloneView( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
00624
00625
00626 CGIterationState<ScalarType,MV> newState;
00627 newState.R = R_0;
00628 block_cg_iter->initializeCG(newState);
00629
00630 while(1) {
00631
00632
00633 try {
00634 block_cg_iter->iterate();
00635
00637
00638
00639
00641 if ( convTest_->getStatus() == Passed ) {
00642
00643
00644 std::vector<int> convIdx = Teuchos::rcp_dynamic_cast<StatusTestGenResNorm<ScalarType,MV,OP> >(convTest_)->convIndices();
00645
00646
00647
00648 if (convIdx.size() == currRHSIdx.size())
00649 break;
00650
00651
00652 problem_->setCurrLS();
00653
00654
00655 int have = 0;
00656 std::vector<int> unconvIdx(currRHSIdx.size());
00657 for (unsigned int i=0; i<currRHSIdx.size(); ++i) {
00658 bool found = false;
00659 for (unsigned int j=0; j<convIdx.size(); ++j) {
00660 if (currRHSIdx[i] == convIdx[j]) {
00661 found = true;
00662 break;
00663 }
00664 }
00665 if (!found) {
00666 currIdx2[have] = currIdx2[i];
00667 currRHSIdx[have++] = currRHSIdx[i];
00668 }
00669 }
00670 currRHSIdx.resize(have);
00671 currIdx2.resize(have);
00672
00673
00674 problem_->setLSIndex( currRHSIdx );
00675
00676
00677 std::vector<MagnitudeType> norms;
00678 R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
00679 for (int i=0; i<have; ++i) { currIdx2[i] = i; }
00680
00681
00682 CGIterationState<ScalarType,MV> defstate;
00683 defstate.R = R_0;
00684 block_cg_iter->initializeCG(defstate);
00685 }
00686
00688
00689
00690
00692 else if ( maxIterTest_->getStatus() == Passed ) {
00693
00694 isConverged = false;
00695 break;
00696 }
00697
00699
00700
00701
00702
00704
00705 else {
00706 TEST_FOR_EXCEPTION(true,std::logic_error,
00707 "Belos::PseudoBlockCGSolMgr::solve(): Invalid return from PseudoBlockCGIter::iterate().");
00708 }
00709 }
00710 catch (const std::exception &e) {
00711 printer_->stream(Errors) << "Error! Caught std::exception in PseudoBlockCGIter::iterate() at iteration "
00712 << block_cg_iter->getNumIters() << std::endl
00713 << e.what() << std::endl;
00714 throw;
00715 }
00716 }
00717
00718
00719 problem_->setCurrLS();
00720
00721
00722 startPtr += numCurrRHS;
00723 numRHS2Solve -= numCurrRHS;
00724
00725 if ( numRHS2Solve > 0 ) {
00726
00727 numCurrRHS = numRHS2Solve;
00728 currIdx.resize( numCurrRHS );
00729 currIdx2.resize( numCurrRHS );
00730 for (int i=0; i<numCurrRHS; ++i)
00731 { currIdx[i] = startPtr+i; currIdx2[i] = i; }
00732
00733
00734 problem_->setLSIndex( currIdx );
00735 }
00736 else {
00737 currIdx.resize( numRHS2Solve );
00738 }
00739
00740 }
00741
00742 }
00743
00744
00745 sTest_->print( printer_->stream(FinalSummary) );
00746
00747
00748 Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
00749
00750
00751 numIters_ = maxIterTest_->getNumIters();
00752
00753 if (!isConverged ) {
00754 return Unconverged;
00755 }
00756 return Converged;
00757 }
00758
00759
00760 template<class ScalarType, class MV, class OP>
00761 std::string PseudoBlockCGSolMgr<ScalarType,MV,OP>::description() const
00762 {
00763 std::ostringstream oss;
00764 oss << "Belos::PseudoBlockCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
00765 oss << "{";
00766 oss << "}";
00767 return oss.str();
00768 }
00769
00770 }
00771
00772 #endif