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_;
505 #if defined(_WIN32) && defined(__clang__) 507 __builtin_constant_p(reinterpret_cast<const std::ostream*>(&std::cout));
532 template<
class ScalarType,
class MV,
class OP>
534 outputStream_(
Teuchos::rcp(outputStream_default_,false)),
538 achievedTol_(
Teuchos::ScalarTraits<typename
Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
539 maxRestarts_(maxRestarts_default_),
540 maxIters_(maxIters_default_),
542 blockSize_(blockSize_default_),
543 numBlocks_(numBlocks_default_),
544 verbosity_(verbosity_default_),
545 outputStyle_(outputStyle_default_),
546 outputFreq_(outputFreq_default_),
547 defQuorum_(defQuorum_default_),
548 showMaxResNormOnly_(showMaxResNormOnly_default_),
549 orthoType_(orthoType_default_),
550 impResScale_(impResScale_default_),
551 expResScale_(expResScale_default_),
553 label_(label_default_),
561 template<
class ScalarType,
class MV,
class OP>
564 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
566 outputStream_(
Teuchos::rcp(outputStream_default_,false)),
570 achievedTol_(
Teuchos::ScalarTraits<typename
Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
571 maxRestarts_(maxRestarts_default_),
572 maxIters_(maxIters_default_),
574 blockSize_(blockSize_default_),
575 numBlocks_(numBlocks_default_),
576 verbosity_(verbosity_default_),
577 outputStyle_(outputStyle_default_),
578 outputFreq_(outputFreq_default_),
579 defQuorum_(defQuorum_default_),
580 showMaxResNormOnly_(showMaxResNormOnly_default_),
581 orthoType_(orthoType_default_),
582 impResScale_(impResScale_default_),
583 expResScale_(expResScale_default_),
585 label_(label_default_),
591 TEUCHOS_TEST_FOR_EXCEPTION(
problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
600 template<
class ScalarType,
class MV,
class OP>
605 using Teuchos::ParameterList;
606 using Teuchos::parameterList;
608 using Teuchos::rcp_dynamic_cast;
611 if (params_ == Teuchos::null) {
612 params_ = parameterList (*getValidParameters ());
619 params->validateParameters (*getValidParameters (), 0);
623 if (params->isParameter (
"Maximum Restarts")) {
624 maxRestarts_ = params->get (
"Maximum Restarts", maxRestarts_default_);
627 params_->set (
"Maximum Restarts", maxRestarts_);
631 if (params->isParameter (
"Maximum Iterations")) {
632 maxIters_ = params->get (
"Maximum Iterations", maxIters_default_);
635 params_->set (
"Maximum Iterations", maxIters_);
636 if (! maxIterTest_.is_null ()) {
637 maxIterTest_->setMaxIters (maxIters_);
642 if (params->isParameter (
"Block Size")) {
643 blockSize_ = params->get (
"Block Size", blockSize_default_);
644 TEUCHOS_TEST_FOR_EXCEPTION(
645 blockSize_ <= 0, std::invalid_argument,
646 "Belos::PseudoBlockGmresSolMgr::setParameters: " 647 "The \"Block Size\" parameter must be strictly positive, " 648 "but you specified a value of " << blockSize_ <<
".");
651 params_->set (
"Block Size", blockSize_);
655 if (params->isParameter (
"Num Blocks")) {
656 numBlocks_ = params->get (
"Num Blocks", numBlocks_default_);
657 TEUCHOS_TEST_FOR_EXCEPTION(
658 numBlocks_ <= 0, std::invalid_argument,
659 "Belos::PseudoBlockGmresSolMgr::setParameters: " 660 "The \"Num Blocks\" parameter must be strictly positive, " 661 "but you specified a value of " << numBlocks_ <<
".");
664 params_->set (
"Num Blocks", numBlocks_);
668 if (params->isParameter (
"Timer Label")) {
669 const std::string tempLabel = params->get (
"Timer Label", label_default_);
672 if (tempLabel != label_) {
674 params_->set (
"Timer Label", label_);
675 const std::string solveLabel =
676 label_ +
": PseudoBlockGmresSolMgr total solve time";
677 #ifdef BELOS_TEUCHOS_TIME_MONITOR 678 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
679 #endif // BELOS_TEUCHOS_TIME_MONITOR 680 if (ortho_ != Teuchos::null) {
681 ortho_->setLabel( label_ );
688 if (params->isParameter (
"Verbosity")) {
689 if (Teuchos::isParameterType<int> (*params,
"Verbosity")) {
690 verbosity_ = params->get (
"Verbosity", verbosity_default_);
692 verbosity_ = (int) Teuchos::getParameter<Belos::MsgType> (*params,
"Verbosity");
696 params_->set (
"Verbosity", verbosity_);
697 if (! printer_.is_null ()) {
698 printer_->setVerbosity (verbosity_);
703 if (params->isParameter (
"Output Style")) {
704 if (Teuchos::isParameterType<int> (*params,
"Output Style")) {
705 outputStyle_ = params->get (
"Output Style", outputStyle_default_);
707 outputStyle_ = (int) Teuchos::getParameter<Belos::OutputType> (*params,
"Output Style");
711 params_->set (
"Output Style", outputStyle_);
712 if (! outputTest_.is_null ()) {
719 if (params->isSublist (
"User Status Tests")) {
720 Teuchos::ParameterList userStatusTestsList = params->sublist(
"User Status Tests",
true);
721 if ( userStatusTestsList.numParams() > 0 ) {
722 std::string userCombo_string = params->get<std::string>(
"User Status Tests Combo Type",
"SEQ");
724 setUserConvStatusTest( testFactory->buildStatusTests(userStatusTestsList), testFactory->stringToComboType(userCombo_string) );
725 taggedTests_ = testFactory->getTaggedTests();
731 if (params->isParameter (
"Output Stream")) {
732 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> > (*params,
"Output Stream");
735 params_->set(
"Output Stream", outputStream_);
736 if (! printer_.is_null ()) {
737 printer_->setOStream (outputStream_);
743 if (params->isParameter (
"Output Frequency")) {
744 outputFreq_ = params->get (
"Output Frequency", outputFreq_default_);
748 params_->set (
"Output Frequency", outputFreq_);
749 if (! outputTest_.is_null ()) {
750 outputTest_->setOutputFrequency (outputFreq_);
755 if (printer_.is_null ()) {
760 bool changedOrthoType =
false;
761 if (params->isParameter (
"Orthogonalization")) {
762 std::string tempOrthoType = params->get (
"Orthogonalization", orthoType_default_);
763 if (tempOrthoType != orthoType_) {
764 orthoType_ = tempOrthoType;
765 changedOrthoType =
true;
768 params_->set(
"Orthogonalization", orthoType_);
771 if (params->isParameter (
"Orthogonalization Constant")) {
772 if (params->isType<
MagnitudeType> (
"Orthogonalization Constant")) {
773 orthoKappa_ = params->get (
"Orthogonalization Constant",
777 orthoKappa_ = params->get (
"Orthogonalization Constant",
782 params_->set (
"Orthogonalization Constant", orthoKappa_);
783 if (orthoType_ ==
"DGKS") {
784 if (orthoKappa_ > 0 && ! ortho_.is_null() && !changedOrthoType) {
786 rcp_dynamic_cast<ortho_type> (ortho_)->setDepTol (orthoKappa_);
792 if (ortho_.is_null() || changedOrthoType) {
794 Teuchos::RCP<Teuchos::ParameterList> paramsOrtho;
795 if (orthoType_==
"DGKS" && orthoKappa_ > 0) {
796 paramsOrtho->set (
"depTol", orthoKappa_ );
799 ortho_ = factory.
makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
805 if (params->isParameter (
"Convergence Tolerance")) {
806 if (params->isType<
MagnitudeType> (
"Convergence Tolerance")) {
807 convtol_ = params->get (
"Convergence Tolerance",
815 params_->set (
"Convergence Tolerance", convtol_);
816 if (! impConvTest_.is_null ()) {
817 impConvTest_->setTolerance (convtol_);
819 if (! expConvTest_.is_null ()) {
820 expConvTest_->setTolerance (convtol_);
825 bool userDefinedResidualScalingUpdated =
false;
826 if (params->isParameter (
"User Defined Residual Scaling")) {
828 if (params->isType<
MagnitudeType> (
"User Defined Residual Scaling")) {
829 tempResScaleFactor = params->get (
"User Defined Residual Scaling",
833 tempResScaleFactor = params->get (
"User Defined Residual Scaling",
838 if (resScaleFactor_ != tempResScaleFactor) {
839 resScaleFactor_ = tempResScaleFactor;
840 userDefinedResidualScalingUpdated =
true;
843 if(userDefinedResidualScalingUpdated)
845 if (! params->isParameter (
"Implicit Residual Scaling") && ! impConvTest_.is_null ()) {
847 if(impResScale_ ==
"User Provided")
850 catch (std::exception& e) {
855 if (! params->isParameter (
"Explicit Residual Scaling") && ! expConvTest_.is_null ()) {
857 if(expResScale_ ==
"User Provided")
860 catch (std::exception& e) {
869 if (params->isParameter (
"Implicit Residual Scaling")) {
870 const std::string tempImpResScale =
871 Teuchos::getParameter<std::string> (*params,
"Implicit Residual Scaling");
874 if (impResScale_ != tempImpResScale) {
876 impResScale_ = tempImpResScale;
879 params_->set (
"Implicit Residual Scaling", impResScale_);
880 if (! impConvTest_.is_null ()) {
882 if(impResScale_ ==
"User Provided")
883 impConvTest_->defineScaleForm (impResScaleType,
Belos::TwoNorm, resScaleFactor_);
887 catch (std::exception& e) {
893 else if (userDefinedResidualScalingUpdated) {
896 if (! impConvTest_.is_null ()) {
898 if(impResScale_ ==
"User Provided")
899 impConvTest_->defineScaleForm (impResScaleType,
Belos::TwoNorm, resScaleFactor_);
901 catch (std::exception& e) {
909 if (params->isParameter (
"Explicit Residual Scaling")) {
910 const std::string tempExpResScale =
911 Teuchos::getParameter<std::string> (*params,
"Explicit Residual Scaling");
914 if (expResScale_ != tempExpResScale) {
916 expResScale_ = tempExpResScale;
919 params_->set (
"Explicit Residual Scaling", expResScale_);
920 if (! expConvTest_.is_null ()) {
922 if(expResScale_ ==
"User Provided")
923 expConvTest_->defineScaleForm (expResScaleType,
Belos::TwoNorm, resScaleFactor_);
927 catch (std::exception& e) {
933 else if (userDefinedResidualScalingUpdated) {
936 if (! expConvTest_.is_null ()) {
938 if(expResScale_ ==
"User Provided")
939 expConvTest_->defineScaleForm (expResScaleType,
Belos::TwoNorm, resScaleFactor_);
941 catch (std::exception& e) {
950 if (params->isParameter (
"Show Maximum Residual Norm Only")) {
951 showMaxResNormOnly_ =
952 Teuchos::getParameter<bool> (*params,
"Show Maximum Residual Norm Only");
955 params_->set (
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
956 if (! impConvTest_.is_null ()) {
957 impConvTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
959 if (! expConvTest_.is_null ()) {
960 expConvTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
967 if (params->isParameter(
"Deflation Quorum")) {
968 defQuorum_ = params->get(
"Deflation Quorum", defQuorum_);
969 TEUCHOS_TEST_FOR_EXCEPTION(
970 defQuorum_ > blockSize_, std::invalid_argument,
971 "Belos::PseudoBlockGmresSolMgr::setParameters: " 972 "The \"Deflation Quorum\" parameter (= " << defQuorum_ <<
") must not be " 973 "larger than \"Block Size\" (= " << blockSize_ <<
").");
974 params_->set (
"Deflation Quorum", defQuorum_);
975 if (! impConvTest_.is_null ()) {
976 impConvTest_->setQuorum (defQuorum_);
978 if (! expConvTest_.is_null ()) {
979 expConvTest_->setQuorum (defQuorum_);
984 if (timerSolve_ == Teuchos::null) {
985 std::string solveLabel = label_ +
": PseudoBlockGmresSolMgr total solve time";
986 #ifdef BELOS_TEUCHOS_TIME_MONITOR 987 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
996 template<
class ScalarType,
class MV,
class OP>
1003 userConvStatusTest_ = userConvStatusTest;
1004 comboType_ = comboType;
1007 template<
class ScalarType,
class MV,
class OP>
1013 debugStatusTest_ = debugStatusTest;
1018 template<
class ScalarType,
class MV,
class OP>
1019 Teuchos::RCP<const Teuchos::ParameterList>
1022 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
1023 if (is_null(validPL)) {
1024 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
1029 pl= Teuchos::rcp(
new Teuchos::ParameterList() );
1031 "The relative residual tolerance that needs to be achieved by the\n" 1032 "iterative solver in order for the linear system to be declared converged.");
1033 pl->set(
"Maximum Restarts", static_cast<int>(maxRestarts_default_),
1034 "The maximum number of restarts allowed for each\n" 1035 "set of RHS solved.");
1036 pl->set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
1037 "The maximum number of block iterations allowed for each\n" 1038 "set of RHS solved.");
1039 pl->set(
"Num Blocks", static_cast<int>(numBlocks_default_),
1040 "The maximum number of vectors allowed in the Krylov subspace\n" 1041 "for each set of RHS solved.");
1042 pl->set(
"Block Size", static_cast<int>(blockSize_default_),
1043 "The number of RHS solved simultaneously.");
1044 pl->set(
"Verbosity", static_cast<int>(verbosity_default_),
1045 "What type(s) of solver information should be outputted\n" 1046 "to the output stream.");
1047 pl->set(
"Output Style", static_cast<int>(outputStyle_default_),
1048 "What style is used for the solver information outputted\n" 1049 "to the output stream.");
1050 pl->set(
"Output Frequency", static_cast<int>(outputFreq_default_),
1051 "How often convergence information should be outputted\n" 1052 "to the output stream.");
1053 pl->set(
"Deflation Quorum", static_cast<int>(defQuorum_default_),
1054 "The number of linear systems that need to converge before\n" 1055 "they are deflated. This number should be <= block size.");
1056 pl->set(
"Output Stream", Teuchos::rcp(outputStream_default_,
false),
1057 "A reference-counted pointer to the output stream where all\n" 1058 "solver output is sent.");
1059 pl->set(
"Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
1060 "When convergence information is printed, only show the maximum\n" 1061 "relative residual norm when the block size is greater than one.");
1062 pl->set(
"Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
1063 "The type of scaling used in the implicit residual convergence test.");
1064 pl->set(
"Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
1065 "The type of scaling used in the explicit residual convergence test.");
1066 pl->set(
"Timer Label", static_cast<const char *>(label_default_),
1067 "The string to use as a prefix for the timer labels.");
1068 pl->set(
"Orthogonalization", static_cast<const char *>(orthoType_default_),
1069 "The type of orthogonalization to use.");
1071 "The constant used by DGKS orthogonalization to determine\n" 1072 "whether another step of classical Gram-Schmidt is necessary.");
1073 pl->sublist(
"User Status Tests");
1074 pl->set(
"User Status Tests Combo Type",
"SEQ",
1075 "Type of logical combination operation of user-defined\n" 1076 "and/or solver-specific status tests.");
1083 template<
class ScalarType,
class MV,
class OP>
1095 if ( !Teuchos::is_null(problem_->getLeftPrec()) ) {
1102 Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
1103 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
1104 if(impResScale_ ==
"User Provided")
1109 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1110 impConvTest_ = tmpImpConvTest;
1113 Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
1114 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
1115 tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit,
Belos::TwoNorm );
1116 if(expResScale_ ==
"User Provided")
1120 tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1121 expConvTest_ = tmpExpConvTest;
1124 convTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
1130 Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
1131 Teuchos::rcp(
new StatusTestImpResNorm_t( convtol_, defQuorum_ ) );
1132 if(impResScale_ ==
"User Provided")
1136 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1137 impConvTest_ = tmpImpConvTest;
1140 expConvTest_ = impConvTest_;
1141 convTest_ = impConvTest_;
1144 if (nonnull(userConvStatusTest_) ) {
1146 Teuchos::RCP<StatusTestCombo_t> tmpComboTest = Teuchos::rcp_dynamic_cast<StatusTestCombo_t>(userConvStatusTest_);
1147 if (tmpComboTest != Teuchos::null) {
1148 std::vector<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > tmpVec = tmpComboTest->getStatusTests();
1150 const int numResTests =
static_cast<int>(tmpVec.size());
1151 auto newConvTest = Teuchos::rcp(
new StatusTestCombo_t(comboType_));
1152 newConvTest->addStatusTest(convTest_);
1153 for (
int j = 0; j < numResTests; ++j) {
1154 newConvTest->addStatusTest(tmpVec[j]);
1156 convTest_ = newConvTest;
1160 convTest_ = Teuchos::rcp(
1161 new StatusTestCombo_t( comboType_, convTest_, userConvStatusTest_ ) );
1169 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
1172 if (nonnull(debugStatusTest_) ) {
1174 Teuchos::rcp_dynamic_cast<StatusTestCombo_t>(sTest_)->addStatusTest( debugStatusTest_ );
1183 std::string solverDesc =
" Pseudo Block Gmres ";
1184 outputTest_->setSolverDesc( solverDesc );
1195 template<
class ScalarType,
class MV,
class OP>
1201 if (!isSet_) { setParameters( params_ ); }
1204 "Belos::PseudoBlockGmresSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1207 if (!isSTSet_ || (!expResTest_ && !Teuchos::is_null(problem_->getLeftPrec())) ) {
1209 "Belos::BlockGmresSolMgr::solve(): Linear problem and requested status tests are incompatible.");
1214 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1215 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1217 std::vector<int> currIdx( numCurrRHS );
1218 blockSize_ = numCurrRHS;
1219 for (
int i=0; i<numCurrRHS; ++i)
1220 { currIdx[i] = startPtr+i; }
1223 problem_->setLSIndex( currIdx );
1227 Teuchos::ParameterList plist;
1228 plist.set(
"Num Blocks",numBlocks_);
1231 outputTest_->reset();
1232 loaDetected_ =
false;
1235 bool isConverged =
true;
1240 Teuchos::RCP<PseudoBlockGmresIter<ScalarType,MV,OP> > block_gmres_iter
1245 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1246 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1249 while ( numRHS2Solve > 0 ) {
1252 std::vector<int> convRHSIdx;
1253 std::vector<int> currRHSIdx( currIdx );
1254 currRHSIdx.resize(numCurrRHS);
1257 block_gmres_iter->setNumBlocks( numBlocks_ );
1260 block_gmres_iter->resetNumIters();
1263 outputTest_->resetNumCalls();
1269 std::vector<int> index(1);
1270 Teuchos::RCP<MV> tmpV, R_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1271 newState.
V.resize( blockSize_ );
1272 newState.
Z.resize( blockSize_ );
1273 for (
int i=0; i<blockSize_; ++i) {
1275 tmpV = MVT::CloneViewNonConst( *R_0, index );
1278 Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > tmpZ
1279 = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>( 1 ));
1282 int rank = ortho_->normalize( *tmpV, tmpZ );
1284 "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
1286 newState.
V[i] = tmpV;
1287 newState.
Z[i] = tmpZ;
1291 block_gmres_iter->initialize(newState);
1292 int numRestarts = 0;
1298 block_gmres_iter->iterate();
1305 if ( convTest_->getStatus() ==
Passed ) {
1307 if ( expConvTest_->getLOADetected() ) {
1317 loaDetected_ =
true;
1319 "Belos::PseudoBlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
1320 isConverged =
false;
1324 std::vector<int> convIdx = expConvTest_->convIndices();
1328 if (convIdx.size() == currRHSIdx.size())
1335 problem_->setCurrLS();
1339 int curDim = oldState.
curDim;
1344 std::vector<int> oldRHSIdx( currRHSIdx );
1345 std::vector<int> defRHSIdx;
1346 for (
unsigned int i=0; i<currRHSIdx.size(); ++i) {
1348 for (
unsigned int j=0; j<convIdx.size(); ++j) {
1349 if (currRHSIdx[i] == convIdx[j]) {
1355 defRHSIdx.push_back( i );
1358 defState.
V.push_back( Teuchos::rcp_const_cast<MV>( oldState.
V[i] ) );
1359 defState.
Z.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.
Z[i] ) );
1360 defState.
H.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseMatrix<int,ScalarType> >( oldState.
H[i] ) );
1361 defState.
sn.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.
sn[i] ) );
1362 defState.
cs.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,MagnitudeType> >(oldState.
cs[i] ) );
1363 currRHSIdx[have] = currRHSIdx[i];
1367 defRHSIdx.resize(currRHSIdx.size()-have);
1368 currRHSIdx.resize(have);
1372 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1373 Teuchos::RCP<MV> defUpdate = MVT::CloneViewNonConst( *update, defRHSIdx );
1376 problem_->setLSIndex( convIdx );
1379 problem_->updateSolution( defUpdate,
true );
1383 problem_->setLSIndex( currRHSIdx );
1386 defState.
curDim = curDim;
1389 block_gmres_iter->initialize(defState);
1396 else if ( maxIterTest_->getStatus() ==
Passed ) {
1398 isConverged =
false;
1406 else if ( block_gmres_iter->getCurSubspaceDim() == block_gmres_iter->getMaxSubspaceDim() ) {
1408 if ( numRestarts >= maxRestarts_ ) {
1409 isConverged =
false;
1414 printer_->stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
1417 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1418 problem_->updateSolution( update,
true );
1425 newstate.
V.resize(currRHSIdx.size());
1426 newstate.
Z.resize(currRHSIdx.size());
1430 R_0 = MVT::Clone( *(problem_->getInitPrecResVec()), currRHSIdx.size() );
1431 problem_->computeCurrPrecResVec( &*R_0 );
1432 for (
unsigned int i=0; i<currRHSIdx.size(); ++i) {
1435 tmpV = MVT::CloneViewNonConst( *R_0, index );
1438 Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > tmpZ
1439 = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>( 1 ));
1442 int rank = ortho_->normalize( *tmpV, tmpZ );
1444 "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors after the restart.");
1446 newstate.
V[i] = tmpV;
1447 newstate.
Z[i] = tmpZ;
1452 block_gmres_iter->initialize(newstate);
1464 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
1465 "Belos::PseudoBlockGmresSolMgr::solve(): Invalid return from PseudoBlockGmresIter::iterate().");
1471 block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
1474 sTest_->checkStatus( &*block_gmres_iter );
1475 if (convTest_->getStatus() !=
Passed)
1476 isConverged =
false;
1479 catch (
const std::exception &e) {
1480 printer_->stream(
Errors) <<
"Error! Caught std::exception in PseudoBlockGmresIter::iterate() at iteration " 1481 << block_gmres_iter->getNumIters() << std::endl
1482 << e.what() << std::endl;
1489 if (nonnull(userConvStatusTest_)) {
1491 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1492 problem_->updateSolution( update,
true );
1494 else if (nonnull(expConvTest_->getSolution())) {
1496 Teuchos::RCP<MV> newX = expConvTest_->getSolution();
1497 Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1498 MVT::MvAddMv( 0.0, *newX, 1.0, *newX, *curX );
1502 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1503 problem_->updateSolution( update,
true );
1507 problem_->setCurrLS();
1510 startPtr += numCurrRHS;
1511 numRHS2Solve -= numCurrRHS;
1512 if ( numRHS2Solve > 0 ) {
1513 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1515 blockSize_ = numCurrRHS;
1516 currIdx.resize( numCurrRHS );
1517 for (
int i=0; i<numCurrRHS; ++i)
1518 { currIdx[i] = startPtr+i; }
1521 if (defQuorum_ > blockSize_) {
1522 if (impConvTest_ != Teuchos::null)
1523 impConvTest_->setQuorum( blockSize_ );
1524 if (expConvTest_ != Teuchos::null)
1525 expConvTest_->setQuorum( blockSize_ );
1529 problem_->setLSIndex( currIdx );
1532 currIdx.resize( numRHS2Solve );
1543 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1548 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1552 numIters_ = maxIterTest_->getNumIters();
1565 const std::vector<MagnitudeType>* pTestValues = NULL;
1567 pTestValues = expConvTest_->getTestValue();
1568 if (pTestValues == NULL || pTestValues->size() < 1) {
1569 pTestValues = impConvTest_->getTestValue();
1574 pTestValues = impConvTest_->getTestValue();
1576 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1577 "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's " 1578 "getTestValue() method returned NULL. Please report this bug to the " 1579 "Belos developers.");
1580 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1581 "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's " 1582 "getTestValue() method returned a vector of length zero. Please report " 1583 "this bug to the Belos developers.");
1588 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1591 if (!isConverged || loaDetected_) {
1598 template<
class ScalarType,
class MV,
class OP>
1601 std::ostringstream out;
1603 out <<
"\"Belos::PseudoBlockGmresSolMgr\": {";
1604 if (this->getObjectLabel () !=
"") {
1605 out <<
"Label: " << this->getObjectLabel () <<
", ";
1607 out <<
"Num Blocks: " << numBlocks_
1608 <<
", Maximum Iterations: " << maxIters_
1609 <<
", Maximum Restarts: " << maxRestarts_
1610 <<
", Convergence Tolerance: " << convtol_
1616 template<
class ScalarType,
class MV,
class OP>
1620 const Teuchos::EVerbosityLevel verbLevel)
const 1622 using Teuchos::TypeNameTraits;
1623 using Teuchos::VERB_DEFAULT;
1624 using Teuchos::VERB_NONE;
1625 using Teuchos::VERB_LOW;
1632 const Teuchos::EVerbosityLevel vl =
1633 (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
1635 if (vl != VERB_NONE) {
1636 Teuchos::OSTab tab0 (out);
1638 out <<
"\"Belos::PseudoBlockGmresSolMgr\":" << endl;
1639 Teuchos::OSTab tab1 (out);
1640 out <<
"Template parameters:" << endl;
1642 Teuchos::OSTab tab2 (out);
1643 out <<
"ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1644 <<
"MV: " << TypeNameTraits<MV>::name () << endl
1645 <<
"OP: " << TypeNameTraits<OP>::name () << endl;
1647 if (this->getObjectLabel () !=
"") {
1648 out <<
"Label: " << this->getObjectLabel () << endl;
1650 out <<
"Num Blocks: " << numBlocks_ << endl
1651 <<
"Maximum Iterations: " << maxIters_ << endl
1652 <<
"Maximum Restarts: " << maxRestarts_ << endl
1653 <<
"Convergence Tolerance: " << convtol_ << endl;
static const double orthoKappa
DGKS orthogonalization constant.
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > > H
The current Hessenberg matrix.
Collection of types and exceptions used within the Belos solvers.
ComboType
The test can be either the AND of all the component tests, or the OR of all the component tests...
virtual ~PseudoBlockGmresSolMgr()
Destructor.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶ms) override
Set the parameters the solver manager should use to solve the linear problem.
Class which manages the output and verbosity of the Belos solvers.
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
The current linear problem to solve.
Teuchos::RCP< Teuchos::ParameterList > params_
PseudoBlockGmresSolMgr()
Empty constructor.
PseudoBlockGmresSolMgrLinearProblemFailure(const std::string &what_arg)
MagnitudeType resScaleFactor_
ScaleType
The type of scaling to use on the residual norm value.
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Convergence test using the implicit residual norm(s), with an explicit residual norm(s) check for los...
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > debugStatusTest_
PseudoBlockGmresSolMgrOrthoFailure(const std::string &what_arg)
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
A factory class for generating StatusTestOutput objects.
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...
An abstract class of StatusTest for stopping criteria using residual norms.
std::vector< Teuchos::RCP< const MV > > V
The current Krylov basis.
An implementation of StatusTestResNorm using a family of residual norms.
static constexpr int outputFreq_default_
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
static const double convTol
Default convergence tolerance.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to PseudoBlockGmresIter state variables.
static constexpr const char * label_default_
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< std::ostream > outputStream_
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > convTest_
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
Traits class which defines basic operations on multivectors.
MagnitudeType achievedTol_
A Belos::StatusTest class for specifying a maximum number of iterations.
MultiVecTraits< ScalarType, MV > MVT
static constexpr int blockSize_default_
ResetType
How to reset the solver.
Pure virtual base class which describes the basic interface for a solver manager. ...
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.
static constexpr int numBlocks_default_
static constexpr int defQuorum_default_
static constexpr bool showMaxResNormOnly_default_
static constexpr int verbosity_default_
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
PseudoBlockGmresSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, MagnitudeType > > > cs
A linear system to solve, and its associated information.
Teuchos::RCP< std::map< std::string, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > > > taggedTests_
Class which describes the linear problem to be solved by the iterative solver.
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
Interface to standard and "pseudoblock" GMRES.
int getNumIters() const override
Iteration count for the most recent call to solve().
Teuchos::RCP< Teuchos::Time > timerSolve_
ReturnType
Whether the Belos solve converged for all linear systems.
bool checkStatusTest()
Check current status tests against current linear problem.
The Belos::SolverManager is a templated virtual base class that defines the basic interface that any ...
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
A list of valid default parameters for this solver.
static constexpr const char * orthoType_default_
static constexpr int maxRestarts_default_
static constexpr const char * impResScale_default_
bool isLOADetected() const override
Whether a "loss of accuracy" was detected during the last solve().
const StatusTestResNorm< ScalarType, MV, OP > * getResidualStatusTest() const
Return the residual status test.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > Z
The current right-hand side of the least squares system RY = Z.
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.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem to solve.
Belos concrete class for performing the pseudo-block GMRES iteration.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
ComboType getComboType() const
Return the type of combination (OR, AND, or SEQ).
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
A class for extending the status testing capabilities of Belos via logical combinations.
std::string description() const override
Return a one-line description of this object.
Class which defines basic traits for the operator type.
int curDim
The current dimension of the reduction.
Teuchos::ScalarTraits< MagnitudeType > MT
static constexpr const char * expResScale_default_
static constexpr int maxIters_default_
Parent class to all Belos exceptions.
ReturnType solve() override
This method performs possibly repeated calls to the underlying linear solver's iterate() routine unti...
Default parameters common to most Belos solvers.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > sn
The current Given's rotation coefficients.
static const double resScaleFactor
User-defined residual scaling factor.
static constexpr int outputStyle_default_
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< StatusTestResNorm< ScalarType, MV, OP > > expConvTest_
OperatorTraits< ScalarType, MV, OP > OPT
static constexpr std::ostream * outputStream_default_
This class implements the pseudo-block GMRES iteration, where a block Krylov subspace is constructed ...
Factory to build a set of status tests from a parameter list.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
The current parameters for this solver.
Teuchos::RCP< OutputManager< ScalarType > > printer_
PseudoBlockGmresIterOrthoFailure is thrown when the orthogonalization manager is unable to generate o...
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
Teuchos::ScalarTraits< ScalarType > SCT
StatusTestCombo< ScalarType, MV, OP >::ComboType comboType_
Teuchos::RCP< StatusTestResNorm< ScalarType, MV, OP > > impConvTest_
PseudoBlockGmresSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate...
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.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > userConvStatusTest_
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
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.
MagnitudeType orthoKappa_