42#ifndef BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP
43#define BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP
60#ifdef BELOS_TEUCHOS_TIME_MONITOR
61#include "Teuchos_TimeMonitor.hpp"
119 template<
class ScalarType,
class MV,
class OP>
125 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
126 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
MagnitudeType;
127 typedef Teuchos::ScalarTraits<MagnitudeType>
MT;
254 const Teuchos::RCP<Teuchos::ParameterList> &pl );
260 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
283 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
377 void setParameters (
const Teuchos::RCP<Teuchos::ParameterList> ¶ms)
override;
439 describe (Teuchos::FancyOStream& out,
440 const Teuchos::EVerbosityLevel verbLevel =
441 Teuchos::Describable::verbLevel_default)
const override;
467 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
problem_;
476 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
sTest_;
482 Teuchos::RCP<std::map<std::string, Teuchos::RCP<StatusTest<ScalarType, MV, OP> > > >
taggedTests_;
485 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> >
ortho_;
525template<
class ScalarType,
class MV,
class OP>
527 outputStream_(Teuchos::rcpFromRef(std::cout)),
528 taggedTests_(Teuchos::null),
531 achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
532 maxRestarts_(maxRestarts_default_),
533 maxIters_(maxIters_default_),
535 blockSize_(blockSize_default_),
536 numBlocks_(numBlocks_default_),
537 verbosity_(verbosity_default_),
538 outputStyle_(outputStyle_default_),
539 outputFreq_(outputFreq_default_),
540 defQuorum_(defQuorum_default_),
541 showMaxResNormOnly_(showMaxResNormOnly_default_),
542 orthoType_(orthoType_default_),
543 impResScale_(impResScale_default_),
544 expResScale_(expResScale_default_),
546 label_(label_default_),
554template<
class ScalarType,
class MV,
class OP>
557 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
559 outputStream_(Teuchos::rcpFromRef(std::cout)),
560 taggedTests_(Teuchos::null),
563 achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
564 maxRestarts_(maxRestarts_default_),
565 maxIters_(maxIters_default_),
567 blockSize_(blockSize_default_),
568 numBlocks_(numBlocks_default_),
569 verbosity_(verbosity_default_),
570 outputStyle_(outputStyle_default_),
571 outputFreq_(outputFreq_default_),
572 defQuorum_(defQuorum_default_),
573 showMaxResNormOnly_(showMaxResNormOnly_default_),
574 orthoType_(orthoType_default_),
575 impResScale_(impResScale_default_),
576 expResScale_(expResScale_default_),
578 label_(label_default_),
584 TEUCHOS_TEST_FOR_EXCEPTION(
problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
593template<
class ScalarType,
class MV,
class OP>
596setParameters (
const Teuchos::RCP<Teuchos::ParameterList>& params)
598 using Teuchos::ParameterList;
599 using Teuchos::parameterList;
601 using Teuchos::rcp_dynamic_cast;
604 if (params_ == Teuchos::null) {
605 params_ = parameterList (*getValidParameters ());
612 params->validateParameters (*getValidParameters (), 0);
616 if (params->isParameter (
"Maximum Restarts")) {
617 maxRestarts_ = params->get (
"Maximum Restarts", maxRestarts_default_);
620 params_->set (
"Maximum Restarts", maxRestarts_);
624 if (params->isParameter (
"Maximum Iterations")) {
625 maxIters_ = params->get (
"Maximum Iterations", maxIters_default_);
628 params_->set (
"Maximum Iterations", maxIters_);
629 if (! maxIterTest_.is_null ()) {
630 maxIterTest_->setMaxIters (maxIters_);
635 if (params->isParameter (
"Block Size")) {
636 blockSize_ = params->get (
"Block Size", blockSize_default_);
637 TEUCHOS_TEST_FOR_EXCEPTION(
638 blockSize_ <= 0, std::invalid_argument,
639 "Belos::PseudoBlockGmresSolMgr::setParameters: "
640 "The \"Block Size\" parameter must be strictly positive, "
641 "but you specified a value of " << blockSize_ <<
".");
644 params_->set (
"Block Size", blockSize_);
648 if (params->isParameter (
"Num Blocks")) {
649 numBlocks_ = params->get (
"Num Blocks", numBlocks_default_);
650 TEUCHOS_TEST_FOR_EXCEPTION(
651 numBlocks_ <= 0, std::invalid_argument,
652 "Belos::PseudoBlockGmresSolMgr::setParameters: "
653 "The \"Num Blocks\" parameter must be strictly positive, "
654 "but you specified a value of " << numBlocks_ <<
".");
657 params_->set (
"Num Blocks", numBlocks_);
661 if (params->isParameter (
"Timer Label")) {
662 const std::string tempLabel = params->get (
"Timer Label", label_default_);
665 if (tempLabel != label_) {
667 params_->set (
"Timer Label", label_);
668 const std::string solveLabel =
669 label_ +
": PseudoBlockGmresSolMgr total solve time";
670#ifdef BELOS_TEUCHOS_TIME_MONITOR
671 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
673 if (ortho_ != Teuchos::null) {
674 ortho_->setLabel( label_ );
681 if (params->isParameter (
"Verbosity")) {
682 if (Teuchos::isParameterType<int> (*params,
"Verbosity")) {
683 verbosity_ = params->get (
"Verbosity", verbosity_default_);
685 verbosity_ = (int) Teuchos::getParameter<Belos::MsgType> (*params,
"Verbosity");
689 params_->set (
"Verbosity", verbosity_);
690 if (! printer_.is_null ()) {
691 printer_->setVerbosity (verbosity_);
696 if (params->isParameter (
"Output Style")) {
697 if (Teuchos::isParameterType<int> (*params,
"Output Style")) {
698 outputStyle_ = params->get (
"Output Style", outputStyle_default_);
700 outputStyle_ = (int) Teuchos::getParameter<Belos::OutputType> (*params,
"Output Style");
704 params_->set (
"Output Style", outputStyle_);
705 if (! outputTest_.is_null ()) {
712 if (params->isSublist (
"User Status Tests")) {
713 Teuchos::ParameterList userStatusTestsList = params->sublist(
"User Status Tests",
true);
714 if ( userStatusTestsList.numParams() > 0 ) {
715 std::string userCombo_string = params->get<std::string>(
"User Status Tests Combo Type",
"SEQ");
717 setUserConvStatusTest( testFactory->buildStatusTests(userStatusTestsList), testFactory->stringToComboType(userCombo_string) );
718 taggedTests_ = testFactory->getTaggedTests();
724 if (params->isParameter (
"Output Stream")) {
725 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> > (*params,
"Output Stream");
728 params_->set(
"Output Stream", outputStream_);
729 if (! printer_.is_null ()) {
730 printer_->setOStream (outputStream_);
736 if (params->isParameter (
"Output Frequency")) {
737 outputFreq_ = params->get (
"Output Frequency", outputFreq_default_);
741 params_->set (
"Output Frequency", outputFreq_);
742 if (! outputTest_.is_null ()) {
743 outputTest_->setOutputFrequency (outputFreq_);
748 if (printer_.is_null ()) {
753 bool changedOrthoType =
false;
754 if (params->isParameter (
"Orthogonalization")) {
755 std::string tempOrthoType = params->get (
"Orthogonalization", orthoType_default_);
756 if (tempOrthoType != orthoType_) {
757 orthoType_ = tempOrthoType;
758 changedOrthoType =
true;
761 params_->set(
"Orthogonalization", orthoType_);
764 if (params->isParameter (
"Orthogonalization Constant")) {
765 if (params->isType<
MagnitudeType> (
"Orthogonalization Constant")) {
766 orthoKappa_ = params->get (
"Orthogonalization Constant",
770 orthoKappa_ = params->get (
"Orthogonalization Constant",
775 params_->set (
"Orthogonalization Constant", orthoKappa_);
776 if (orthoType_ ==
"DGKS") {
777 if (orthoKappa_ > 0 && ! ortho_.is_null() && !changedOrthoType) {
779 rcp_dynamic_cast<ortho_type> (ortho_)->
setDepTol (orthoKappa_);
785 if (ortho_.is_null() || changedOrthoType) {
787 Teuchos::RCP<Teuchos::ParameterList> paramsOrtho;
788 if (orthoType_==
"DGKS" && orthoKappa_ > 0) {
789 paramsOrtho->set (
"depTol", orthoKappa_ );
792 ortho_ = factory.
makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
798 if (params->isParameter (
"Convergence Tolerance")) {
799 if (params->isType<
MagnitudeType> (
"Convergence Tolerance")) {
800 convtol_ = params->get (
"Convergence Tolerance",
808 params_->set (
"Convergence Tolerance", convtol_);
809 if (! impConvTest_.is_null ()) {
810 impConvTest_->setTolerance (convtol_);
812 if (! expConvTest_.is_null ()) {
813 expConvTest_->setTolerance (convtol_);
818 bool userDefinedResidualScalingUpdated =
false;
819 if (params->isParameter (
"User Defined Residual Scaling")) {
821 if (params->isType<
MagnitudeType> (
"User Defined Residual Scaling")) {
822 tempResScaleFactor = params->get (
"User Defined Residual Scaling",
826 tempResScaleFactor = params->get (
"User Defined Residual Scaling",
831 if (resScaleFactor_ != tempResScaleFactor) {
832 resScaleFactor_ = tempResScaleFactor;
833 userDefinedResidualScalingUpdated =
true;
836 if(userDefinedResidualScalingUpdated)
838 if (! params->isParameter (
"Implicit Residual Scaling") && ! impConvTest_.is_null ()) {
840 if(impResScale_ ==
"User Provided")
843 catch (std::exception& e) {
848 if (! params->isParameter (
"Explicit Residual Scaling") && ! expConvTest_.is_null ()) {
850 if(expResScale_ ==
"User Provided")
853 catch (std::exception& e) {
862 if (params->isParameter (
"Implicit Residual Scaling")) {
863 const std::string tempImpResScale =
864 Teuchos::getParameter<std::string> (*params,
"Implicit Residual Scaling");
867 if (impResScale_ != tempImpResScale) {
869 impResScale_ = tempImpResScale;
872 params_->set (
"Implicit Residual Scaling", impResScale_);
873 if (! impConvTest_.is_null ()) {
875 if(impResScale_ ==
"User Provided")
876 impConvTest_->defineScaleForm (impResScaleType,
Belos::TwoNorm, resScaleFactor_);
880 catch (std::exception& e) {
886 else if (userDefinedResidualScalingUpdated) {
889 if (! impConvTest_.is_null ()) {
891 if(impResScale_ ==
"User Provided")
892 impConvTest_->defineScaleForm (impResScaleType,
Belos::TwoNorm, resScaleFactor_);
894 catch (std::exception& e) {
902 if (params->isParameter (
"Explicit Residual Scaling")) {
903 const std::string tempExpResScale =
904 Teuchos::getParameter<std::string> (*params,
"Explicit Residual Scaling");
907 if (expResScale_ != tempExpResScale) {
909 expResScale_ = tempExpResScale;
912 params_->set (
"Explicit Residual Scaling", expResScale_);
913 if (! expConvTest_.is_null ()) {
915 if(expResScale_ ==
"User Provided")
916 expConvTest_->defineScaleForm (expResScaleType,
Belos::TwoNorm, resScaleFactor_);
920 catch (std::exception& e) {
926 else if (userDefinedResidualScalingUpdated) {
929 if (! expConvTest_.is_null ()) {
931 if(expResScale_ ==
"User Provided")
932 expConvTest_->defineScaleForm (expResScaleType,
Belos::TwoNorm, resScaleFactor_);
934 catch (std::exception& e) {
943 if (params->isParameter (
"Show Maximum Residual Norm Only")) {
944 showMaxResNormOnly_ =
945 Teuchos::getParameter<bool> (*params,
"Show Maximum Residual Norm Only");
948 params_->set (
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
949 if (! impConvTest_.is_null ()) {
950 impConvTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
952 if (! expConvTest_.is_null ()) {
953 expConvTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
960 if (params->isParameter(
"Deflation Quorum")) {
961 defQuorum_ = params->get(
"Deflation Quorum", defQuorum_);
962 TEUCHOS_TEST_FOR_EXCEPTION(
963 defQuorum_ > blockSize_, std::invalid_argument,
964 "Belos::PseudoBlockGmresSolMgr::setParameters: "
965 "The \"Deflation Quorum\" parameter (= " << defQuorum_ <<
") must not be "
966 "larger than \"Block Size\" (= " << blockSize_ <<
").");
967 params_->set (
"Deflation Quorum", defQuorum_);
968 if (! impConvTest_.is_null ()) {
969 impConvTest_->setQuorum (defQuorum_);
971 if (! expConvTest_.is_null ()) {
972 expConvTest_->setQuorum (defQuorum_);
977 if (timerSolve_ == Teuchos::null) {
978 std::string solveLabel = label_ +
": PseudoBlockGmresSolMgr total solve time";
979#ifdef BELOS_TEUCHOS_TIME_MONITOR
980 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
989template<
class ScalarType,
class MV,
class OP>
996 userConvStatusTest_ = userConvStatusTest;
997 comboType_ = comboType;
1000template<
class ScalarType,
class MV,
class OP>
1006 debugStatusTest_ = debugStatusTest;
1011template<
class ScalarType,
class MV,
class OP>
1012Teuchos::RCP<const Teuchos::ParameterList>
1015 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
1016 if (is_null(validPL)) {
1017 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
1022 pl= Teuchos::rcp(
new Teuchos::ParameterList() );
1024 "The relative residual tolerance that needs to be achieved by the\n"
1025 "iterative solver in order for the linear system to be declared converged.");
1026 pl->set(
"Maximum Restarts",
static_cast<int>(maxRestarts_default_),
1027 "The maximum number of restarts allowed for each\n"
1028 "set of RHS solved.");
1029 pl->set(
"Maximum Iterations",
static_cast<int>(maxIters_default_),
1030 "The maximum number of block iterations allowed for each\n"
1031 "set of RHS solved.");
1032 pl->set(
"Num Blocks",
static_cast<int>(numBlocks_default_),
1033 "The maximum number of vectors allowed in the Krylov subspace\n"
1034 "for each set of RHS solved.");
1035 pl->set(
"Block Size",
static_cast<int>(blockSize_default_),
1036 "The number of RHS solved simultaneously.");
1037 pl->set(
"Verbosity",
static_cast<int>(verbosity_default_),
1038 "What type(s) of solver information should be outputted\n"
1039 "to the output stream.");
1040 pl->set(
"Output Style",
static_cast<int>(outputStyle_default_),
1041 "What style is used for the solver information outputted\n"
1042 "to the output stream.");
1043 pl->set(
"Output Frequency",
static_cast<int>(outputFreq_default_),
1044 "How often convergence information should be outputted\n"
1045 "to the output stream.");
1046 pl->set(
"Deflation Quorum",
static_cast<int>(defQuorum_default_),
1047 "The number of linear systems that need to converge before\n"
1048 "they are deflated. This number should be <= block size.");
1049 pl->set(
"Output Stream", Teuchos::rcpFromRef(std::cout),
1050 "A reference-counted pointer to the output stream where all\n"
1051 "solver output is sent.");
1052 pl->set(
"Show Maximum Residual Norm Only",
static_cast<bool>(showMaxResNormOnly_default_),
1053 "When convergence information is printed, only show the maximum\n"
1054 "relative residual norm when the block size is greater than one.");
1055 pl->set(
"Implicit Residual Scaling",
static_cast<const char *
>(impResScale_default_),
1056 "The type of scaling used in the implicit residual convergence test.");
1057 pl->set(
"Explicit Residual Scaling",
static_cast<const char *
>(expResScale_default_),
1058 "The type of scaling used in the explicit residual convergence test.");
1059 pl->set(
"Timer Label",
static_cast<const char *
>(label_default_),
1060 "The string to use as a prefix for the timer labels.");
1061 pl->set(
"Orthogonalization",
static_cast<const char *
>(orthoType_default_),
1062 "The type of orthogonalization to use.");
1064 "The constant used by DGKS orthogonalization to determine\n"
1065 "whether another step of classical Gram-Schmidt is necessary.");
1066 pl->sublist(
"User Status Tests");
1067 pl->set(
"User Status Tests Combo Type",
"SEQ",
1068 "Type of logical combination operation of user-defined\n"
1069 "and/or solver-specific status tests.");
1076template<
class ScalarType,
class MV,
class OP>
1088 if ( !Teuchos::is_null(problem_->getLeftPrec()) ) {
1095 Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
1096 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
1097 if(impResScale_ ==
"User Provided")
1102 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1103 impConvTest_ = tmpImpConvTest;
1106 Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
1107 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
1108 tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit,
Belos::TwoNorm );
1109 if(expResScale_ ==
"User Provided")
1113 tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1114 expConvTest_ = tmpExpConvTest;
1117 convTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
1123 Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
1124 Teuchos::rcp(
new StatusTestImpResNorm_t( convtol_, defQuorum_ ) );
1125 if(impResScale_ ==
"User Provided")
1129 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1130 impConvTest_ = tmpImpConvTest;
1133 expConvTest_ = impConvTest_;
1134 convTest_ = impConvTest_;
1137 if (nonnull(userConvStatusTest_) ) {
1139 Teuchos::RCP<StatusTestCombo_t> tmpComboTest = Teuchos::rcp_dynamic_cast<StatusTestCombo_t>(userConvStatusTest_);
1140 if (tmpComboTest != Teuchos::null) {
1141 std::vector<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > tmpVec = tmpComboTest->getStatusTests();
1142 comboType_ = tmpComboTest->getComboType();
1143 const int numResTests =
static_cast<int>(tmpVec.size());
1144 auto newConvTest = Teuchos::rcp(
new StatusTestCombo_t(comboType_));
1145 newConvTest->addStatusTest(convTest_);
1146 for (
int j = 0; j < numResTests; ++j) {
1147 newConvTest->addStatusTest(tmpVec[j]);
1149 convTest_ = newConvTest;
1153 convTest_ = Teuchos::rcp(
1154 new StatusTestCombo_t( comboType_, convTest_, userConvStatusTest_ ) );
1162 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
1165 if (nonnull(debugStatusTest_) ) {
1167 Teuchos::rcp_dynamic_cast<StatusTestCombo_t>(sTest_)->addStatusTest( debugStatusTest_ );
1176 std::string solverDesc =
" Pseudo Block Gmres ";
1177 outputTest_->setSolverDesc( solverDesc );
1188template<
class ScalarType,
class MV,
class OP>
1194 if (!isSet_) { setParameters( params_ ); }
1197 "Belos::PseudoBlockGmresSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1200 if (!isSTSet_ || (!expResTest_ && !Teuchos::is_null(problem_->getLeftPrec())) ) {
1202 "Belos::BlockGmresSolMgr::solve(): Linear problem and requested status tests are incompatible.");
1207 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1208 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1210 std::vector<int> currIdx( numCurrRHS );
1211 blockSize_ = numCurrRHS;
1212 for (
int i=0; i<numCurrRHS; ++i)
1213 { currIdx[i] = startPtr+i; }
1216 problem_->setLSIndex( currIdx );
1220 Teuchos::ParameterList plist;
1221 plist.set(
"Num Blocks",numBlocks_);
1224 outputTest_->reset();
1225 loaDetected_ =
false;
1228 bool isConverged =
true;
1233 Teuchos::RCP<PseudoBlockGmresIter<ScalarType,MV,OP> > block_gmres_iter
1238#ifdef BELOS_TEUCHOS_TIME_MONITOR
1239 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1242 while ( numRHS2Solve > 0 ) {
1245 std::vector<int> convRHSIdx;
1246 std::vector<int> currRHSIdx( currIdx );
1247 currRHSIdx.resize(numCurrRHS);
1250 block_gmres_iter->setNumBlocks( numBlocks_ );
1253 block_gmres_iter->resetNumIters();
1256 outputTest_->resetNumCalls();
1262 std::vector<int> index(1);
1263 Teuchos::RCP<MV> tmpV, R_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1264 newState.
V.resize( blockSize_ );
1265 newState.
Z.resize( blockSize_ );
1266 for (
int i=0; i<blockSize_; ++i) {
1268 tmpV = MVT::CloneViewNonConst( *R_0, index );
1271 Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > tmpZ
1272 = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>( 1 ));
1275 int rank = ortho_->normalize( *tmpV, tmpZ );
1277 "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
1279 newState.
V[i] = tmpV;
1280 newState.
Z[i] = tmpZ;
1284 block_gmres_iter->initialize(newState);
1285 int numRestarts = 0;
1291 block_gmres_iter->iterate();
1298 if ( convTest_->getStatus() ==
Passed ) {
1300 if ( expConvTest_->getLOADetected() ) {
1310 loaDetected_ =
true;
1312 "Belos::PseudoBlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
1313 isConverged =
false;
1317 std::vector<int> convIdx = expConvTest_->convIndices();
1321 if (convIdx.size() == currRHSIdx.size())
1328 problem_->setCurrLS();
1332 int curDim = oldState.
curDim;
1337 std::vector<int> oldRHSIdx( currRHSIdx );
1338 std::vector<int> defRHSIdx;
1339 for (
unsigned int i=0; i<currRHSIdx.size(); ++i) {
1341 for (
unsigned int j=0; j<convIdx.size(); ++j) {
1342 if (currRHSIdx[i] == convIdx[j]) {
1348 defRHSIdx.push_back( i );
1351 defState.
V.push_back( Teuchos::rcp_const_cast<MV>( oldState.
V[i] ) );
1352 defState.
Z.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.
Z[i] ) );
1353 defState.
H.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseMatrix<int,ScalarType> >( oldState.
H[i] ) );
1354 defState.
sn.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.
sn[i] ) );
1355 defState.
cs.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,MagnitudeType> >(oldState.
cs[i] ) );
1356 currRHSIdx[have] = currRHSIdx[i];
1360 defRHSIdx.resize(currRHSIdx.size()-have);
1361 currRHSIdx.resize(have);
1365 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1366 Teuchos::RCP<MV> defUpdate = MVT::CloneViewNonConst( *update, defRHSIdx );
1369 problem_->setLSIndex( convIdx );
1372 problem_->updateSolution( defUpdate,
true );
1376 problem_->setLSIndex( currRHSIdx );
1379 defState.
curDim = curDim;
1382 block_gmres_iter->initialize(defState);
1389 else if ( maxIterTest_->getStatus() ==
Passed ) {
1391 isConverged =
false;
1399 else if ( block_gmres_iter->getCurSubspaceDim() == block_gmres_iter->getMaxSubspaceDim() ) {
1401 if ( numRestarts >= maxRestarts_ ) {
1402 isConverged =
false;
1407 printer_->stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
1410 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1411 problem_->updateSolution( update,
true );
1418 newstate.
V.resize(currRHSIdx.size());
1419 newstate.
Z.resize(currRHSIdx.size());
1423 R_0 = MVT::Clone( *(problem_->getInitPrecResVec()), currRHSIdx.size() );
1424 problem_->computeCurrPrecResVec( &*R_0 );
1425 for (
unsigned int i=0; i<currRHSIdx.size(); ++i) {
1428 tmpV = MVT::CloneViewNonConst( *R_0, index );
1431 Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > tmpZ
1432 = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>( 1 ));
1435 int rank = ortho_->normalize( *tmpV, tmpZ );
1437 "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors after the restart.");
1439 newstate.
V[i] = tmpV;
1440 newstate.
Z[i] = tmpZ;
1445 block_gmres_iter->initialize(newstate);
1457 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
1458 "Belos::PseudoBlockGmresSolMgr::solve(): Invalid return from PseudoBlockGmresIter::iterate().");
1464 block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
1467 sTest_->checkStatus( &*block_gmres_iter );
1468 if (convTest_->getStatus() !=
Passed)
1469 isConverged =
false;
1472 catch (
const std::exception &e) {
1473 printer_->stream(
Errors) <<
"Error! Caught std::exception in PseudoBlockGmresIter::iterate() at iteration "
1474 << block_gmres_iter->getNumIters() << std::endl
1475 << e.what() << std::endl;
1482 if (nonnull(userConvStatusTest_)) {
1484 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1485 problem_->updateSolution( update,
true );
1487 else if (nonnull(expConvTest_->getSolution())) {
1489 Teuchos::RCP<MV> newX = expConvTest_->getSolution();
1490 Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1491 MVT::MvAddMv( 0.0, *newX, 1.0, *newX, *curX );
1495 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1496 problem_->updateSolution( update,
true );
1500 problem_->setCurrLS();
1503 startPtr += numCurrRHS;
1504 numRHS2Solve -= numCurrRHS;
1505 if ( numRHS2Solve > 0 ) {
1506 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1508 blockSize_ = numCurrRHS;
1509 currIdx.resize( numCurrRHS );
1510 for (
int i=0; i<numCurrRHS; ++i)
1511 { currIdx[i] = startPtr+i; }
1514 if (defQuorum_ > blockSize_) {
1515 if (impConvTest_ != Teuchos::null)
1516 impConvTest_->setQuorum( blockSize_ );
1517 if (expConvTest_ != Teuchos::null)
1518 expConvTest_->setQuorum( blockSize_ );
1522 problem_->setLSIndex( currIdx );
1525 currIdx.resize( numRHS2Solve );
1536#ifdef BELOS_TEUCHOS_TIME_MONITOR
1541 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1545 numIters_ = maxIterTest_->getNumIters();
1558 const std::vector<MagnitudeType>* pTestValues = NULL;
1560 pTestValues = expConvTest_->getTestValue();
1561 if (pTestValues == NULL || pTestValues->size() < 1) {
1562 pTestValues = impConvTest_->getTestValue();
1567 pTestValues = impConvTest_->getTestValue();
1569 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1570 "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's "
1571 "getTestValue() method returned NULL. Please report this bug to the "
1572 "Belos developers.");
1573 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1574 "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's "
1575 "getTestValue() method returned a vector of length zero. Please report "
1576 "this bug to the Belos developers.");
1581 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1584 if (!isConverged || loaDetected_) {
1591template<
class ScalarType,
class MV,
class OP>
1594 std::ostringstream out;
1596 out <<
"\"Belos::PseudoBlockGmresSolMgr\": {";
1597 if (this->getObjectLabel () !=
"") {
1598 out <<
"Label: " << this->getObjectLabel () <<
", ";
1600 out <<
"Num Blocks: " << numBlocks_
1601 <<
", Maximum Iterations: " << maxIters_
1602 <<
", Maximum Restarts: " << maxRestarts_
1603 <<
", Convergence Tolerance: " << convtol_
1609template<
class ScalarType,
class MV,
class OP>
1612describe (Teuchos::FancyOStream &out,
1613 const Teuchos::EVerbosityLevel verbLevel)
const
1615 using Teuchos::TypeNameTraits;
1616 using Teuchos::VERB_DEFAULT;
1617 using Teuchos::VERB_NONE;
1618 using Teuchos::VERB_LOW;
1625 const Teuchos::EVerbosityLevel vl =
1626 (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
1628 if (vl != VERB_NONE) {
1629 Teuchos::OSTab tab0 (out);
1631 out <<
"\"Belos::PseudoBlockGmresSolMgr\":" << endl;
1632 Teuchos::OSTab tab1 (out);
1633 out <<
"Template parameters:" << endl;
1635 Teuchos::OSTab tab2 (out);
1636 out <<
"ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1637 <<
"MV: " << TypeNameTraits<MV>::name () << endl
1638 <<
"OP: " << TypeNameTraits<OP>::name () << endl;
1640 if (this->getObjectLabel () !=
"") {
1641 out <<
"Label: " << this->getObjectLabel () << endl;
1643 out <<
"Num Blocks: " << numBlocks_ << endl
1644 <<
"Maximum Iterations: " << maxIters_ << endl
1645 <<
"Maximum Restarts: " << maxRestarts_ << endl
1646 <<
"Convergence Tolerance: " << convtol_ << endl;
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
Class which manages the output and verbosity of the Belos solvers.
Belos concrete class for performing the pseudo-block GMRES iteration.
Pure virtual base class which describes the basic interface for a solver manager.
A factory class for generating StatusTestOutput objects.
Collection of types and exceptions used within the Belos solvers.
Parent class to all Belos exceptions.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
void setDepTol(const MagnitudeType dep_tol)
Set parameter for re-orthogonalization threshhold.
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Enumeration of all valid Belos (Mat)OrthoManager classes.
Teuchos::RCP< Belos::MatOrthoManager< Scalar, MV, OP > > makeMatOrthoManager(const std::string &ortho, const Teuchos::RCP< const OP > &M, const Teuchos::RCP< OutputManager< Scalar > > &, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Return an instance of the specified MatOrthoManager subclass.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
PseudoBlockGmresIterOrthoFailure is thrown when the orthogonalization manager is unable to generate o...
This class implements the pseudo-block GMRES iteration, where a block Krylov subspace is constructed ...
PseudoBlockGmresSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i....
PseudoBlockGmresSolMgrLinearProblemFailure(const std::string &what_arg)
PseudoBlockGmresSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate...
PseudoBlockGmresSolMgrOrthoFailure(const std::string &what_arg)
Interface to standard and "pseudoblock" GMRES.
static constexpr int defQuorum_default_
static constexpr int outputFreq_default_
Teuchos::RCP< std::ostream > outputStream_
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
MagnitudeType achievedTol_
MagnitudeType orthoKappa_
Teuchos::RCP< std::map< std::string, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > > > taggedTests_
bool isLOADetected() const override
Whether a "loss of accuracy" was detected during the last solve().
OperatorTraits< ScalarType, MV, OP > OPT
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const override
Print the object with the given verbosity level to a FancyOStream.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
A list of valid default parameters for this solver.
static constexpr int numBlocks_default_
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > debugStatusTest_
Teuchos::ScalarTraits< MagnitudeType > MT
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > convTest_
PseudoBlockGmresSolMgr()
Empty constructor.
static constexpr const char * orthoType_default_
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
static constexpr const char * label_default_
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶ms) override
Set the parameters the solver manager should use to solve the linear problem.
ReturnType solve() override
This method performs possibly repeated calls to the underlying linear solver's iterate() routine unti...
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
static constexpr const char * expResScale_default_
static constexpr int outputStyle_default_
MultiVecTraits< ScalarType, MV > MVT
const StatusTestResNorm< ScalarType, MV, OP > * getResidualStatusTest() const
Return the residual status test.
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
The current linear problem to solve.
void setDebugStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &debugStatusTest) override
Set a debug status test that will be checked at the same time as the top-level status test.
MagnitudeType resScaleFactor_
Teuchos::RCP< Teuchos::Time > timerSolve_
static constexpr int blockSize_default_
std::string description() const override
Return a one-line description of this object.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Teuchos::RCP< StatusTestResNorm< ScalarType, MV, OP > > impConvTest_
Teuchos::RCP< StatusTestResNorm< ScalarType, MV, OP > > expConvTest_
Teuchos::ScalarTraits< ScalarType > SCT
virtual void setUserConvStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &userConvStatusTest, const typename StatusTestCombo< ScalarType, MV, OP >::ComboType &comboType=StatusTestCombo< ScalarType, MV, OP >::SEQ) override
Set a custom status test.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
bool checkStatusTest()
Check current status tests against current linear problem.
Teuchos::RCP< OutputManager< ScalarType > > printer_
StatusTestCombo< ScalarType, MV, OP >::ComboType comboType_
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
virtual ~PseudoBlockGmresSolMgr()
Destructor.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
static constexpr int verbosity_default_
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
The current parameters for this solver.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem to solve.
static constexpr int maxIters_default_
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
static constexpr int maxRestarts_default_
static constexpr const char * impResScale_default_
Teuchos::RCP< Teuchos::ParameterList > params_
static constexpr bool showMaxResNormOnly_default_
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > userConvStatusTest_
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
int getNumIters() const override
Iteration count for the most recent call to solve().
The Belos::SolverManager is a templated virtual base class that defines the basic interface that any ...
A class for extending the status testing capabilities of Belos via logical combinations.
ComboType
The test can be either the AND of all the component tests, or the OR of all the component tests,...
Factory to build a set of status tests from a parameter list.
An implementation of StatusTestResNorm using a family of residual norms.
Convergence test using the implicit residual norm(s), with an explicit residual norm(s) check for los...
A Belos::StatusTest class for specifying a maximum number of iterations.
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
An abstract class of StatusTest for stopping criteria using residual norms.
A pure virtual class for defining the status tests for the Belos iterative solvers.
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
ReturnType
Whether the Belos solve converged for all linear systems.
ScaleType
The type of scaling to use on the residual norm value.
ResetType
How to reset the solver.
Default parameters common to most Belos solvers.
static const double resScaleFactor
User-defined residual scaling factor.
static const double convTol
Default convergence tolerance.
static const double orthoKappa
DGKS orthogonalization constant.
Structure to contain pointers to PseudoBlockGmresIter state variables.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > sn
The current Given's rotation coefficients.
std::vector< Teuchos::RCP< const MV > > V
The current Krylov basis.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > Z
The current right-hand side of the least squares system RY = Z.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > > H
The current Hessenberg matrix.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, MagnitudeType > > > cs
int curDim
The current dimension of the reduction.