42 #ifndef BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP 43 #define BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP 59 #ifdef HAVE_BELOS_TSQR 61 #endif // HAVE_BELOS_TSQR 65 #include "Teuchos_BLAS.hpp" 66 #ifdef BELOS_TEUCHOS_TIME_MONITOR 67 #include "Teuchos_TimeMonitor.hpp" 125 template<
class ScalarType,
class MV,
class OP>
131 typedef Teuchos::ScalarTraits<ScalarType> SCT;
132 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
133 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
259 const Teuchos::RCP<Teuchos::ParameterList> &pl );
265 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
288 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
289 return Teuchos::tuple(timerSolve_);
379 void setParameters (
const Teuchos::RCP<Teuchos::ParameterList> ¶ms)
override;
441 describe (Teuchos::FancyOStream& out,
442 const Teuchos::EVerbosityLevel verbLevel =
443 Teuchos::Describable::verbLevel_default)
const override;
466 bool checkStatusTest();
469 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
472 Teuchos::RCP<OutputManager<ScalarType> > printer_;
473 Teuchos::RCP<std::ostream> outputStream_;
476 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > userConvStatusTest_;
477 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugStatusTest_;
478 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
479 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
480 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
481 Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > impConvTest_, expConvTest_;
482 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
484 Teuchos::RCP<std::map<std::string, Teuchos::RCP<StatusTest<ScalarType, MV, OP> > > > taggedTests_;
487 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
490 Teuchos::RCP<Teuchos::ParameterList> params_;
493 static constexpr
int maxRestarts_default_ = 20;
494 static constexpr
int maxIters_default_ = 1000;
495 static constexpr
bool showMaxResNormOnly_default_ =
false;
496 static constexpr
int blockSize_default_ = 1;
497 static constexpr
int numBlocks_default_ = 300;
500 static constexpr
int outputFreq_default_ = -1;
501 static constexpr
int defQuorum_default_ = 1;
502 static constexpr
const char * impResScale_default_ =
"Norm of Preconditioned Initial Residual";
503 static constexpr
const char * expResScale_default_ =
"Norm of Initial Residual";
504 static constexpr
const char * label_default_ =
"Belos";
505 static constexpr
const char * orthoType_default_ =
"DGKS";
506 static constexpr std::ostream * outputStream_default_ = &std::cout;
509 MagnitudeType convtol_, orthoKappa_, achievedTol_;
510 int maxRestarts_, maxIters_, numIters_;
511 int blockSize_, numBlocks_, verbosity_, outputStyle_, outputFreq_, defQuorum_;
512 bool showMaxResNormOnly_;
513 std::string orthoType_;
514 std::string impResScale_, expResScale_;
515 MagnitudeType resScaleFactor_;
519 Teuchos::RCP<Teuchos::Time> timerSolve_;
522 bool isSet_, isSTSet_, expResTest_;
528 template<
class ScalarType,
class MV,
class OP>
530 outputStream_(Teuchos::rcp(outputStream_default_,false)),
531 taggedTests_(Teuchos::null),
534 achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
535 maxRestarts_(maxRestarts_default_),
536 maxIters_(maxIters_default_),
538 blockSize_(blockSize_default_),
539 numBlocks_(numBlocks_default_),
540 verbosity_(verbosity_default_),
541 outputStyle_(outputStyle_default_),
542 outputFreq_(outputFreq_default_),
543 defQuorum_(defQuorum_default_),
544 showMaxResNormOnly_(showMaxResNormOnly_default_),
545 orthoType_(orthoType_default_),
546 impResScale_(impResScale_default_),
547 expResScale_(expResScale_default_),
549 label_(label_default_),
557 template<
class ScalarType,
class MV,
class OP>
560 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
562 outputStream_(Teuchos::rcp(outputStream_default_,false)),
563 taggedTests_(Teuchos::null),
566 achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
567 maxRestarts_(maxRestarts_default_),
568 maxIters_(maxIters_default_),
570 blockSize_(blockSize_default_),
571 numBlocks_(numBlocks_default_),
572 verbosity_(verbosity_default_),
573 outputStyle_(outputStyle_default_),
574 outputFreq_(outputFreq_default_),
575 defQuorum_(defQuorum_default_),
576 showMaxResNormOnly_(showMaxResNormOnly_default_),
577 orthoType_(orthoType_default_),
578 impResScale_(impResScale_default_),
579 expResScale_(expResScale_default_),
581 label_(label_default_),
587 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
596 template<
class ScalarType,
class MV,
class OP>
601 using Teuchos::ParameterList;
602 using Teuchos::parameterList;
604 using Teuchos::rcp_dynamic_cast;
607 if (params_ == Teuchos::null) {
608 params_ = parameterList (*getValidParameters ());
615 params->validateParameters (*getValidParameters (), 0);
619 if (params->isParameter (
"Maximum Restarts")) {
620 maxRestarts_ = params->get (
"Maximum Restarts", maxRestarts_default_);
623 params_->set (
"Maximum Restarts", maxRestarts_);
627 if (params->isParameter (
"Maximum Iterations")) {
628 maxIters_ = params->get (
"Maximum Iterations", maxIters_default_);
631 params_->set (
"Maximum Iterations", maxIters_);
632 if (! maxIterTest_.is_null ()) {
633 maxIterTest_->setMaxIters (maxIters_);
638 if (params->isParameter (
"Block Size")) {
639 blockSize_ = params->get (
"Block Size", blockSize_default_);
640 TEUCHOS_TEST_FOR_EXCEPTION(
641 blockSize_ <= 0, std::invalid_argument,
642 "Belos::PseudoBlockGmresSolMgr::setParameters: " 643 "The \"Block Size\" parameter must be strictly positive, " 644 "but you specified a value of " << blockSize_ <<
".");
647 params_->set (
"Block Size", blockSize_);
651 if (params->isParameter (
"Num Blocks")) {
652 numBlocks_ = params->get (
"Num Blocks", numBlocks_default_);
653 TEUCHOS_TEST_FOR_EXCEPTION(
654 numBlocks_ <= 0, std::invalid_argument,
655 "Belos::PseudoBlockGmresSolMgr::setParameters: " 656 "The \"Num Blocks\" parameter must be strictly positive, " 657 "but you specified a value of " << numBlocks_ <<
".");
660 params_->set (
"Num Blocks", numBlocks_);
664 if (params->isParameter (
"Timer Label")) {
665 const std::string tempLabel = params->get (
"Timer Label", label_default_);
668 if (tempLabel != label_) {
670 params_->set (
"Timer Label", label_);
671 const std::string solveLabel =
672 label_ +
": PseudoBlockGmresSolMgr total solve time";
673 #ifdef BELOS_TEUCHOS_TIME_MONITOR 674 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
675 #endif // BELOS_TEUCHOS_TIME_MONITOR 676 if (ortho_ != Teuchos::null) {
677 ortho_->setLabel( label_ );
683 if (params->isParameter (
"Orthogonalization")) {
684 std::string tempOrthoType = params->get (
"Orthogonalization", orthoType_default_);
685 #ifdef HAVE_BELOS_TSQR 686 TEUCHOS_TEST_FOR_EXCEPTION(
687 tempOrthoType !=
"DGKS" && tempOrthoType !=
"ICGS" &&
688 tempOrthoType !=
"IMGS" && tempOrthoType !=
"TSQR",
689 std::invalid_argument,
690 "Belos::PseudoBlockGmresSolMgr::setParameters: " 691 "The \"Orthogonalization\" parameter must be one of \"DGKS\", \"ICGS\", " 692 "\"IMGS\", or \"TSQR\".");
694 TEUCHOS_TEST_FOR_EXCEPTION(
695 tempOrthoType !=
"DGKS" && tempOrthoType !=
"ICGS" &&
696 tempOrthoType !=
"IMGS",
697 std::invalid_argument,
698 "Belos::PseudoBlockGmresSolMgr::setParameters: " 699 "The \"Orthogonalization\" parameter must be one of \"DGKS\", \"ICGS\", " 701 #endif // HAVE_BELOS_TSQR 703 if (tempOrthoType != orthoType_) {
704 orthoType_ = tempOrthoType;
705 params_->set(
"Orthogonalization", orthoType_);
707 if (orthoType_ ==
"DGKS") {
709 if (orthoKappa_ <= 0) {
710 ortho_ = rcp (
new ortho_type (label_));
713 ortho_ = rcp (
new ortho_type (label_));
714 rcp_dynamic_cast<ortho_type> (ortho_)->setDepTol (orthoKappa_);
717 else if (orthoType_ ==
"ICGS") {
719 ortho_ = rcp (
new ortho_type (label_));
721 else if (orthoType_ ==
"IMGS") {
723 ortho_ = rcp (
new ortho_type (label_));
725 #ifdef HAVE_BELOS_TSQR 726 else if (orthoType_ ==
"TSQR") {
728 ortho_ = rcp (
new ortho_type (label_));
730 #endif // HAVE_BELOS_TSQR 735 if (params->isParameter (
"Orthogonalization Constant")) {
736 if (params->isType<MagnitudeType> (
"Orthogonalization Constant")) {
737 orthoKappa_ = params->get (
"Orthogonalization Constant",
741 orthoKappa_ = params->get (
"Orthogonalization Constant",
746 params_->set (
"Orthogonalization Constant", orthoKappa_);
747 if (orthoType_ ==
"DGKS") {
748 if (orthoKappa_ > 0 && ! ortho_.is_null ()) {
750 rcp_dynamic_cast<ortho_type> (ortho_)->setDepTol (orthoKappa_);
756 if (params->isParameter (
"Verbosity")) {
757 if (Teuchos::isParameterType<int> (*params,
"Verbosity")) {
758 verbosity_ = params->get (
"Verbosity", verbosity_default_);
760 verbosity_ = (int) Teuchos::getParameter<Belos::MsgType> (*params,
"Verbosity");
764 params_->set (
"Verbosity", verbosity_);
765 if (! printer_.is_null ()) {
766 printer_->setVerbosity (verbosity_);
771 if (params->isParameter (
"Output Style")) {
772 if (Teuchos::isParameterType<int> (*params,
"Output Style")) {
773 outputStyle_ = params->get (
"Output Style", outputStyle_default_);
775 outputStyle_ = (int) Teuchos::getParameter<Belos::OutputType> (*params,
"Output Style");
779 params_->set (
"Output Style", verbosity_);
780 if (! outputTest_.is_null ()) {
787 if (params->isSublist (
"User Status Tests")) {
788 Teuchos::ParameterList userStatusTestsList = params->sublist(
"User Status Tests",
true);
789 if ( userStatusTestsList.numParams() > 0 ) {
790 std::string userCombo_string = params->get<std::string>(
"User Status Tests Combo Type",
"SEQ");
792 setUserConvStatusTest( testFactory->buildStatusTests(userStatusTestsList), testFactory->stringToComboType(userCombo_string) );
793 taggedTests_ = testFactory->getTaggedTests();
799 if (params->isParameter (
"Output Stream")) {
800 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> > (*params,
"Output Stream");
803 params_->set(
"Output Stream", outputStream_);
804 if (! printer_.is_null ()) {
805 printer_->setOStream (outputStream_);
811 if (params->isParameter (
"Output Frequency")) {
812 outputFreq_ = params->get (
"Output Frequency", outputFreq_default_);
816 params_->set (
"Output Frequency", outputFreq_);
817 if (! outputTest_.is_null ()) {
818 outputTest_->setOutputFrequency (outputFreq_);
823 if (printer_.is_null ()) {
830 if (params->isParameter (
"Convergence Tolerance")) {
831 if (params->isType<MagnitudeType> (
"Convergence Tolerance")) {
832 convtol_ = params->get (
"Convergence Tolerance",
840 params_->set (
"Convergence Tolerance", convtol_);
841 if (! impConvTest_.is_null ()) {
842 impConvTest_->setTolerance (convtol_);
844 if (! expConvTest_.is_null ()) {
845 expConvTest_->setTolerance (convtol_);
850 bool userDefinedResidualScalingUpdated =
false;
851 if (params->isParameter (
"User Defined Residual Scaling")) {
853 if (params->isType<MagnitudeType> (
"User Defined Residual Scaling")) {
854 tempResScaleFactor = params->get (
"User Defined Residual Scaling",
858 tempResScaleFactor = params->get (
"User Defined Residual Scaling",
863 if (resScaleFactor_ != tempResScaleFactor) {
864 resScaleFactor_ = tempResScaleFactor;
865 userDefinedResidualScalingUpdated =
true;
868 if(userDefinedResidualScalingUpdated)
870 if (! params->isParameter (
"Implicit Residual Scaling") && ! impConvTest_.is_null ()) {
872 if(impResScale_ ==
"User Provided")
875 catch (std::exception& e) {
880 if (! params->isParameter (
"Explicit Residual Scaling") && ! expConvTest_.is_null ()) {
882 if(expResScale_ ==
"User Provided")
885 catch (std::exception& e) {
894 if (params->isParameter (
"Implicit Residual Scaling")) {
895 const std::string tempImpResScale =
896 Teuchos::getParameter<std::string> (*params,
"Implicit Residual Scaling");
899 if (impResScale_ != tempImpResScale) {
901 impResScale_ = tempImpResScale;
904 params_->set (
"Implicit Residual Scaling", impResScale_);
905 if (! impConvTest_.is_null ()) {
907 if(impResScale_ ==
"User Provided")
908 impConvTest_->defineScaleForm (impResScaleType,
Belos::TwoNorm, resScaleFactor_);
912 catch (std::exception& e) {
918 else if (userDefinedResidualScalingUpdated) {
921 if (! impConvTest_.is_null ()) {
923 if(impResScale_ ==
"User Provided")
924 impConvTest_->defineScaleForm (impResScaleType,
Belos::TwoNorm, resScaleFactor_);
926 catch (std::exception& e) {
934 if (params->isParameter (
"Explicit Residual Scaling")) {
935 const std::string tempExpResScale =
936 Teuchos::getParameter<std::string> (*params,
"Explicit Residual Scaling");
939 if (expResScale_ != tempExpResScale) {
941 expResScale_ = tempExpResScale;
944 params_->set (
"Explicit Residual Scaling", expResScale_);
945 if (! expConvTest_.is_null ()) {
947 if(expResScale_ ==
"User Provided")
948 expConvTest_->defineScaleForm (expResScaleType,
Belos::TwoNorm, resScaleFactor_);
952 catch (std::exception& e) {
958 else if (userDefinedResidualScalingUpdated) {
961 if (! expConvTest_.is_null ()) {
963 if(expResScale_ ==
"User Provided")
964 expConvTest_->defineScaleForm (expResScaleType,
Belos::TwoNorm, resScaleFactor_);
966 catch (std::exception& e) {
975 if (params->isParameter (
"Show Maximum Residual Norm Only")) {
976 showMaxResNormOnly_ =
977 Teuchos::getParameter<bool> (*params,
"Show Maximum Residual Norm Only");
980 params_->set (
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
981 if (! impConvTest_.is_null ()) {
982 impConvTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
984 if (! expConvTest_.is_null ()) {
985 expConvTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
992 if (params->isParameter(
"Deflation Quorum")) {
993 defQuorum_ = params->get(
"Deflation Quorum", defQuorum_);
994 TEUCHOS_TEST_FOR_EXCEPTION(
995 defQuorum_ > blockSize_, std::invalid_argument,
996 "Belos::PseudoBlockGmresSolMgr::setParameters: " 997 "The \"Deflation Quorum\" parameter (= " << defQuorum_ <<
") must not be " 998 "larger than \"Block Size\" (= " << blockSize_ <<
").");
999 params_->set (
"Deflation Quorum", defQuorum_);
1000 if (! impConvTest_.is_null ()) {
1001 impConvTest_->setQuorum (defQuorum_);
1003 if (! expConvTest_.is_null ()) {
1004 expConvTest_->setQuorum (defQuorum_);
1009 if (ortho_.is_null ()) {
1010 params_->set(
"Orthogonalization", orthoType_);
1011 if (orthoType_ ==
"DGKS") {
1013 if (orthoKappa_ <= 0) {
1014 ortho_ = rcp (
new ortho_type (label_));
1017 ortho_ = rcp (
new ortho_type (label_));
1018 rcp_dynamic_cast<ortho_type> (ortho_)->setDepTol (orthoKappa_);
1021 else if (orthoType_ ==
"ICGS") {
1023 ortho_ = rcp (
new ortho_type (label_));
1025 else if (orthoType_ ==
"IMGS") {
1027 ortho_ = rcp (
new ortho_type (label_));
1029 #ifdef HAVE_BELOS_TSQR 1030 else if (orthoType_ ==
"TSQR") {
1032 ortho_ = rcp (
new ortho_type (label_));
1034 #endif // HAVE_BELOS_TSQR 1036 #ifdef HAVE_BELOS_TSQR 1037 TEUCHOS_TEST_FOR_EXCEPTION(
1038 orthoType_ !=
"ICGS" && orthoType_ !=
"DGKS" &&
1039 orthoType_ !=
"IMGS" && orthoType_ !=
"TSQR",
1041 "Belos::PseudoBlockGmresSolMgr::setParameters(): " 1042 "Invalid orthogonalization type \"" << orthoType_ <<
"\".");
1044 TEUCHOS_TEST_FOR_EXCEPTION(
1045 orthoType_ !=
"ICGS" && orthoType_ !=
"DGKS" &&
1046 orthoType_ !=
"IMGS",
1048 "Belos::PseudoBlockGmresSolMgr::setParameters(): " 1049 "Invalid orthogonalization type \"" << orthoType_ <<
"\".");
1050 #endif // HAVE_BELOS_TSQR 1055 if (timerSolve_ == Teuchos::null) {
1056 std::string solveLabel = label_ +
": PseudoBlockGmresSolMgr total solve time";
1057 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1058 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
1067 template<
class ScalarType,
class MV,
class OP>
1074 userConvStatusTest_ = userConvStatusTest;
1075 comboType_ = comboType;
1078 template<
class ScalarType,
class MV,
class OP>
1084 debugStatusTest_ = debugStatusTest;
1089 template<
class ScalarType,
class MV,
class OP>
1090 Teuchos::RCP<const Teuchos::ParameterList>
1093 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
1094 if (is_null(validPL)) {
1095 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
1100 pl= Teuchos::rcp(
new Teuchos::ParameterList() );
1102 "The relative residual tolerance that needs to be achieved by the\n" 1103 "iterative solver in order for the linear system to be declared converged.");
1104 pl->set(
"Maximum Restarts", static_cast<int>(maxRestarts_default_),
1105 "The maximum number of restarts allowed for each\n" 1106 "set of RHS solved.");
1107 pl->set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
1108 "The maximum number of block iterations allowed for each\n" 1109 "set of RHS solved.");
1110 pl->set(
"Num Blocks", static_cast<int>(numBlocks_default_),
1111 "The maximum number of vectors allowed in the Krylov subspace\n" 1112 "for each set of RHS solved.");
1113 pl->set(
"Block Size", static_cast<int>(blockSize_default_),
1114 "The number of RHS solved simultaneously.");
1115 pl->set(
"Verbosity", static_cast<int>(verbosity_default_),
1116 "What type(s) of solver information should be outputted\n" 1117 "to the output stream.");
1118 pl->set(
"Output Style", static_cast<int>(outputStyle_default_),
1119 "What style is used for the solver information outputted\n" 1120 "to the output stream.");
1121 pl->set(
"Output Frequency", static_cast<int>(outputFreq_default_),
1122 "How often convergence information should be outputted\n" 1123 "to the output stream.");
1124 pl->set(
"Deflation Quorum", static_cast<int>(defQuorum_default_),
1125 "The number of linear systems that need to converge before\n" 1126 "they are deflated. This number should be <= block size.");
1127 pl->set(
"Output Stream", Teuchos::rcp(outputStream_default_,
false),
1128 "A reference-counted pointer to the output stream where all\n" 1129 "solver output is sent.");
1130 pl->set(
"Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
1131 "When convergence information is printed, only show the maximum\n" 1132 "relative residual norm when the block size is greater than one.");
1133 pl->set(
"Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
1134 "The type of scaling used in the implicit residual convergence test.");
1135 pl->set(
"Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
1136 "The type of scaling used in the explicit residual convergence test.");
1137 pl->set(
"Timer Label", static_cast<const char *>(label_default_),
1138 "The string to use as a prefix for the timer labels.");
1139 pl->set(
"Orthogonalization", static_cast<const char *>(orthoType_default_),
1140 "The type of orthogonalization to use: DGKS, ICGS, IMGS.");
1142 "The constant used by DGKS orthogonalization to determine\n" 1143 "whether another step of classical Gram-Schmidt is necessary.");
1144 pl->sublist(
"User Status Tests");
1145 pl->set(
"User Status Tests Combo Type",
"SEQ",
1146 "Type of logical combination operation of user-defined\n" 1147 "and/or solver-specific status tests.");
1154 template<
class ScalarType,
class MV,
class OP>
1166 if ( !Teuchos::is_null(problem_->getLeftPrec()) ) {
1173 Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
1174 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
1175 if(impResScale_ ==
"User Provided")
1180 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1181 impConvTest_ = tmpImpConvTest;
1184 Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
1185 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
1186 tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit,
Belos::TwoNorm );
1187 if(expResScale_ ==
"User Provided")
1191 tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1192 expConvTest_ = tmpExpConvTest;
1195 convTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
1201 Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
1202 Teuchos::rcp(
new StatusTestImpResNorm_t( convtol_, defQuorum_ ) );
1203 if(impResScale_ ==
"User Provided")
1207 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1208 impConvTest_ = tmpImpConvTest;
1211 expConvTest_ = impConvTest_;
1212 convTest_ = impConvTest_;
1215 if (nonnull(userConvStatusTest_) ) {
1217 convTest_ = Teuchos::rcp(
1218 new StatusTestCombo_t( comboType_, convTest_, userConvStatusTest_ ) );
1225 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
1228 if (nonnull(debugStatusTest_) ) {
1230 Teuchos::rcp_dynamic_cast<StatusTestCombo_t>(sTest_)->addStatusTest( debugStatusTest_ );
1235 StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_, taggedTests_ );
1239 std::string solverDesc =
" Pseudo Block Gmres ";
1240 outputTest_->setSolverDesc( solverDesc );
1251 template<
class ScalarType,
class MV,
class OP>
1257 if (!isSet_) { setParameters( params_ ); }
1259 Teuchos::BLAS<int,ScalarType> blas;
1262 "Belos::PseudoBlockGmresSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1265 if (!isSTSet_ || (!expResTest_ && !Teuchos::is_null(problem_->getLeftPrec())) ) {
1267 "Belos::BlockGmresSolMgr::solve(): Linear problem and requested status tests are incompatible.");
1272 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1273 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1275 std::vector<int> currIdx( numCurrRHS );
1276 blockSize_ = numCurrRHS;
1277 for (
int i=0; i<numCurrRHS; ++i)
1278 { currIdx[i] = startPtr+i; }
1281 problem_->setLSIndex( currIdx );
1285 Teuchos::ParameterList plist;
1286 plist.set(
"Num Blocks",numBlocks_);
1289 outputTest_->reset();
1290 loaDetected_ =
false;
1293 bool isConverged =
true;
1298 Teuchos::RCP<PseudoBlockGmresIter<ScalarType,MV,OP> > block_gmres_iter
1303 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1304 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1307 while ( numRHS2Solve > 0 ) {
1310 std::vector<int> convRHSIdx;
1311 std::vector<int> currRHSIdx( currIdx );
1312 currRHSIdx.resize(numCurrRHS);
1315 block_gmres_iter->setNumBlocks( numBlocks_ );
1318 block_gmres_iter->resetNumIters();
1321 outputTest_->resetNumCalls();
1327 std::vector<int> index(1);
1328 Teuchos::RCP<MV> tmpV, R_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1329 newState.
V.resize( blockSize_ );
1330 newState.
Z.resize( blockSize_ );
1331 for (
int i=0; i<blockSize_; ++i) {
1333 tmpV = MVT::CloneViewNonConst( *R_0, index );
1336 Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > tmpZ
1337 = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>( 1 ));
1340 int rank = ortho_->normalize( *tmpV, tmpZ );
1342 "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
1344 newState.
V[i] = tmpV;
1345 newState.
Z[i] = tmpZ;
1349 block_gmres_iter->initialize(newState);
1350 int numRestarts = 0;
1356 block_gmres_iter->iterate();
1363 if ( convTest_->getStatus() ==
Passed ) {
1365 if ( expConvTest_->getLOADetected() ) {
1375 loaDetected_ =
true;
1377 "Belos::PseudoBlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
1378 isConverged =
false;
1382 std::vector<int> convIdx = expConvTest_->convIndices();
1386 if (convIdx.size() == currRHSIdx.size())
1393 problem_->setCurrLS();
1397 int curDim = oldState.
curDim;
1402 std::vector<int> oldRHSIdx( currRHSIdx );
1403 std::vector<int> defRHSIdx;
1404 for (
unsigned int i=0; i<currRHSIdx.size(); ++i) {
1406 for (
unsigned int j=0; j<convIdx.size(); ++j) {
1407 if (currRHSIdx[i] == convIdx[j]) {
1413 defRHSIdx.push_back( i );
1416 defState.
V.push_back( Teuchos::rcp_const_cast<MV>( oldState.
V[i] ) );
1417 defState.
Z.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.
Z[i] ) );
1418 defState.
H.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseMatrix<int,ScalarType> >( oldState.
H[i] ) );
1419 defState.
sn.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.
sn[i] ) );
1420 defState.
cs.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,MagnitudeType> >(oldState.
cs[i] ) );
1421 currRHSIdx[have] = currRHSIdx[i];
1425 defRHSIdx.resize(currRHSIdx.size()-have);
1426 currRHSIdx.resize(have);
1430 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1431 Teuchos::RCP<MV> defUpdate = MVT::CloneViewNonConst( *update, defRHSIdx );
1434 problem_->setLSIndex( convIdx );
1437 problem_->updateSolution( defUpdate,
true );
1441 problem_->setLSIndex( currRHSIdx );
1444 defState.
curDim = curDim;
1447 block_gmres_iter->initialize(defState);
1454 else if ( maxIterTest_->getStatus() ==
Passed ) {
1456 isConverged =
false;
1464 else if ( block_gmres_iter->getCurSubspaceDim() == block_gmres_iter->getMaxSubspaceDim() ) {
1466 if ( numRestarts >= maxRestarts_ ) {
1467 isConverged =
false;
1472 printer_->stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
1475 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1476 problem_->updateSolution( update,
true );
1483 newstate.
V.resize(currRHSIdx.size());
1484 newstate.
Z.resize(currRHSIdx.size());
1488 R_0 = MVT::Clone( *(problem_->getInitPrecResVec()), currRHSIdx.size() );
1489 problem_->computeCurrPrecResVec( &*R_0 );
1490 for (
unsigned int i=0; i<currRHSIdx.size(); ++i) {
1493 tmpV = MVT::CloneViewNonConst( *R_0, index );
1496 Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > tmpZ
1497 = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>( 1 ));
1500 int rank = ortho_->normalize( *tmpV, tmpZ );
1502 "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors after the restart.");
1504 newstate.
V[i] = tmpV;
1505 newstate.
Z[i] = tmpZ;
1510 block_gmres_iter->initialize(newstate);
1522 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
1523 "Belos::PseudoBlockGmresSolMgr::solve(): Invalid return from PseudoBlockGmresIter::iterate().");
1529 block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
1532 sTest_->checkStatus( &*block_gmres_iter );
1533 if (convTest_->getStatus() !=
Passed)
1534 isConverged =
false;
1537 catch (
const std::exception &e) {
1538 printer_->stream(
Errors) <<
"Error! Caught std::exception in PseudoBlockGmresIter::iterate() at iteration " 1539 << block_gmres_iter->getNumIters() << std::endl
1540 << e.what() << std::endl;
1547 if (nonnull(userConvStatusTest_)) {
1549 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1550 problem_->updateSolution( update,
true );
1552 else if (nonnull(expConvTest_->getSolution())) {
1554 Teuchos::RCP<MV> newX = expConvTest_->getSolution();
1555 Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1556 MVT::MvAddMv( 0.0, *newX, 1.0, *newX, *curX );
1560 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1561 problem_->updateSolution( update,
true );
1565 problem_->setCurrLS();
1568 startPtr += numCurrRHS;
1569 numRHS2Solve -= numCurrRHS;
1570 if ( numRHS2Solve > 0 ) {
1571 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1573 blockSize_ = numCurrRHS;
1574 currIdx.resize( numCurrRHS );
1575 for (
int i=0; i<numCurrRHS; ++i)
1576 { currIdx[i] = startPtr+i; }
1579 if (defQuorum_ > blockSize_) {
1580 if (impConvTest_ != Teuchos::null)
1581 impConvTest_->setQuorum( blockSize_ );
1582 if (expConvTest_ != Teuchos::null)
1583 expConvTest_->setQuorum( blockSize_ );
1587 problem_->setLSIndex( currIdx );
1590 currIdx.resize( numRHS2Solve );
1601 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1606 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1610 numIters_ = maxIterTest_->getNumIters();
1623 const std::vector<MagnitudeType>* pTestValues = NULL;
1625 pTestValues = expConvTest_->getTestValue();
1626 if (pTestValues == NULL || pTestValues->size() < 1) {
1627 pTestValues = impConvTest_->getTestValue();
1632 pTestValues = impConvTest_->getTestValue();
1634 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1635 "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's " 1636 "getTestValue() method returned NULL. Please report this bug to the " 1637 "Belos developers.");
1638 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1639 "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's " 1640 "getTestValue() method returned a vector of length zero. Please report " 1641 "this bug to the Belos developers.");
1646 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1649 if (!isConverged || loaDetected_) {
1656 template<
class ScalarType,
class MV,
class OP>
1659 std::ostringstream out;
1661 out <<
"\"Belos::PseudoBlockGmresSolMgr\": {";
1662 if (this->getObjectLabel () !=
"") {
1663 out <<
"Label: " << this->getObjectLabel () <<
", ";
1665 out <<
"Num Blocks: " << numBlocks_
1666 <<
", Maximum Iterations: " << maxIters_
1667 <<
", Maximum Restarts: " << maxRestarts_
1668 <<
", Convergence Tolerance: " << convtol_
1674 template<
class ScalarType,
class MV,
class OP>
1678 const Teuchos::EVerbosityLevel verbLevel)
const 1680 using Teuchos::TypeNameTraits;
1681 using Teuchos::VERB_DEFAULT;
1682 using Teuchos::VERB_NONE;
1683 using Teuchos::VERB_LOW;
1690 const Teuchos::EVerbosityLevel vl =
1691 (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
1693 if (vl != VERB_NONE) {
1694 Teuchos::OSTab tab0 (out);
1696 out <<
"\"Belos::PseudoBlockGmresSolMgr\":" << endl;
1697 Teuchos::OSTab tab1 (out);
1698 out <<
"Template parameters:" << endl;
1700 Teuchos::OSTab tab2 (out);
1701 out <<
"ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1702 <<
"MV: " << TypeNameTraits<MV>::name () << endl
1703 <<
"OP: " << TypeNameTraits<OP>::name () << endl;
1705 if (this->getObjectLabel () !=
"") {
1706 out <<
"Label: " << this->getObjectLabel () << endl;
1708 out <<
"Num Blocks: " << numBlocks_ << endl
1709 <<
"Maximum Iterations: " << maxIters_ << endl
1710 <<
"Maximum Restarts: " << maxRestarts_ << endl
1711 <<
"Convergence Tolerance: " << convtol_ << endl;
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.
PseudoBlockGmresSolMgr()
Empty constructor.
PseudoBlockGmresSolMgrLinearProblemFailure(const std::string &what_arg)
ScaleType
The type of scaling to use on the residual norm value.
static constexpr double orthoKappa
DGKS orthogonalization constant.
Convergence test using the implicit residual norm(s), with an explicit residual norm(s) check for los...
PseudoBlockGmresSolMgrOrthoFailure(const std::string &what_arg)
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
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...
std::vector< Teuchos::RCP< const MV > > V
The current Krylov basis.
An implementation of StatusTestResNorm using a family of residual norms.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to PseudoBlockGmresIter state variables.
A factory class for generating StatusTestOutput objects.
Iterated Modified Gram-Schmidt (IMGS) implementation of the Belos::OrthoManager class.
Traits class which defines basic operations on multivectors.
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
Pure virtual base class which describes the basic interface for a solver manager. ...
static constexpr double convTol
Default convergence tolerance.
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.
PseudoBlockGmresSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
static constexpr double resScaleFactor
User-defined residual scaling factor.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, MagnitudeType > > > cs
A linear system to solve, and its associated information.
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().
ReturnType
Whether the Belos solve converged for all linear systems.
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.
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
bool isLOADetected() const override
Whether a "loss of accuracy" was detected during the last solve().
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > Z
The current right-hand side of the least squares system RY = Z.
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< 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.
Orthogonalization manager based on Tall Skinny QR (TSQR)
Class which defines basic traits for the operator type.
int curDim
The current dimension of the reduction.
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.
MatOrthoManager subclass using TSQR or DGKS.
Belos header file which uses auto-configuration information to include necessary C++ headers...
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.
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.
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.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.