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 68 #include "Teuchos_BLAS.hpp" 69 #ifdef BELOS_TEUCHOS_TIME_MONITOR 70 #include "Teuchos_TimeMonitor.hpp" 128 template<
class ScalarType,
class MV,
class OP>
134 typedef Teuchos::ScalarTraits<ScalarType> SCT;
135 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
136 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
262 const Teuchos::RCP<Teuchos::ParameterList> &pl );
276 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const;
286 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
287 return Teuchos::tuple(timerSolve_);
377 void setParameters (
const Teuchos::RCP<Teuchos::ParameterList> ¶ms);
385 virtual void setUserConvStatusTest(
394 virtual void setDebugStatusTest(
443 describe (Teuchos::FancyOStream& out,
444 const Teuchos::EVerbosityLevel verbLevel =
445 Teuchos::Describable::verbLevel_default)
const;
448 std::string description ()
const;
468 bool checkStatusTest();
471 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
474 Teuchos::RCP<OutputManager<ScalarType> > printer_;
475 Teuchos::RCP<std::ostream> outputStream_;
478 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > userConvStatusTest_;
479 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugStatusTest_;
480 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
481 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
482 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
483 Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > impConvTest_, expConvTest_;
484 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
487 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
490 Teuchos::RCP<Teuchos::ParameterList> params_;
493 static const MagnitudeType convtol_default_;
494 static const MagnitudeType orthoKappa_default_;
495 static const int maxRestarts_default_;
496 static const int maxIters_default_;
497 static const bool showMaxResNormOnly_default_;
498 static const int blockSize_default_;
499 static const int numBlocks_default_;
500 static const int verbosity_default_;
501 static const int outputStyle_default_;
502 static const int outputFreq_default_;
503 static const int defQuorum_default_;
504 static const std::string impResScale_default_;
505 static const std::string expResScale_default_;
506 static const MagnitudeType resScaleFactor_default_;
507 static const std::string label_default_;
508 static const std::string orthoType_default_;
509 static const Teuchos::RCP<std::ostream> outputStream_default_;
512 MagnitudeType convtol_, orthoKappa_, achievedTol_;
513 int maxRestarts_, maxIters_, numIters_;
514 int blockSize_, numBlocks_, verbosity_, outputStyle_, outputFreq_, defQuorum_;
515 bool showMaxResNormOnly_;
516 std::string orthoType_;
517 std::string impResScale_, expResScale_;
518 MagnitudeType resScaleFactor_;
522 Teuchos::RCP<Teuchos::Time> timerSolve_;
525 bool isSet_, isSTSet_, expResTest_;
531 template<
class ScalarType,
class MV,
class OP>
534 template<
class ScalarType,
class MV,
class OP>
537 template<
class ScalarType,
class MV,
class OP>
540 template<
class ScalarType,
class MV,
class OP>
543 template<
class ScalarType,
class MV,
class OP>
546 template<
class ScalarType,
class MV,
class OP>
549 template<
class ScalarType,
class MV,
class OP>
552 template<
class ScalarType,
class MV,
class OP>
555 template<
class ScalarType,
class MV,
class OP>
558 template<
class ScalarType,
class MV,
class OP>
561 template<
class ScalarType,
class MV,
class OP>
564 template<
class ScalarType,
class MV,
class OP>
567 template<
class ScalarType,
class MV,
class OP>
570 template<
class ScalarType,
class MV,
class OP>
573 template<
class ScalarType,
class MV,
class OP>
576 template<
class ScalarType,
class MV,
class OP>
579 template<
class ScalarType,
class MV,
class OP>
584 template<
class ScalarType,
class MV,
class OP>
586 outputStream_(outputStream_default_),
587 convtol_(convtol_default_),
588 orthoKappa_(orthoKappa_default_),
589 achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
590 maxRestarts_(maxRestarts_default_),
591 maxIters_(maxIters_default_),
593 blockSize_(blockSize_default_),
594 numBlocks_(numBlocks_default_),
595 verbosity_(verbosity_default_),
596 outputStyle_(outputStyle_default_),
597 outputFreq_(outputFreq_default_),
598 defQuorum_(defQuorum_default_),
599 showMaxResNormOnly_(showMaxResNormOnly_default_),
600 orthoType_(orthoType_default_),
601 impResScale_(impResScale_default_),
602 expResScale_(expResScale_default_),
603 resScaleFactor_(resScaleFactor_default_),
604 label_(label_default_),
612 template<
class ScalarType,
class MV,
class OP>
615 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
617 outputStream_(outputStream_default_),
618 convtol_(convtol_default_),
619 orthoKappa_(orthoKappa_default_),
620 achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
621 maxRestarts_(maxRestarts_default_),
622 maxIters_(maxIters_default_),
624 blockSize_(blockSize_default_),
625 numBlocks_(numBlocks_default_),
626 verbosity_(verbosity_default_),
627 outputStyle_(outputStyle_default_),
628 outputFreq_(outputFreq_default_),
629 defQuorum_(defQuorum_default_),
630 showMaxResNormOnly_(showMaxResNormOnly_default_),
631 orthoType_(orthoType_default_),
632 impResScale_(impResScale_default_),
633 expResScale_(expResScale_default_),
634 resScaleFactor_(resScaleFactor_default_),
635 label_(label_default_),
641 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
650 template<
class ScalarType,
class MV,
class OP>
655 using Teuchos::ParameterList;
656 using Teuchos::parameterList;
658 using Teuchos::rcp_dynamic_cast;
661 if (params_ == Teuchos::null) {
668 if (params->isParameter (
"Maximum Restarts")) {
669 maxRestarts_ = params->get (
"Maximum Restarts", maxRestarts_default_);
672 params_->set (
"Maximum Restarts", maxRestarts_);
676 if (params->isParameter (
"Maximum Iterations")) {
677 maxIters_ = params->get (
"Maximum Iterations", maxIters_default_);
680 params_->set (
"Maximum Iterations", maxIters_);
681 if (! maxIterTest_.is_null ()) {
682 maxIterTest_->setMaxIters (maxIters_);
687 if (params->isParameter (
"Block Size")) {
688 blockSize_ = params->get (
"Block Size", blockSize_default_);
689 TEUCHOS_TEST_FOR_EXCEPTION(
690 blockSize_ <= 0, std::invalid_argument,
691 "Belos::PseudoBlockGmresSolMgr::setParameters: " 692 "The \"Block Size\" parameter must be strictly positive, " 693 "but you specified a value of " << blockSize_ <<
".");
696 params_->set (
"Block Size", blockSize_);
700 if (params->isParameter (
"Num Blocks")) {
701 numBlocks_ = params->get (
"Num Blocks", numBlocks_default_);
702 TEUCHOS_TEST_FOR_EXCEPTION(
703 numBlocks_ <= 0, std::invalid_argument,
704 "Belos::PseudoBlockGmresSolMgr::setParameters: " 705 "The \"Num Blocks\" parameter must be strictly positive, " 706 "but you specified a value of " << numBlocks_ <<
".");
709 params_->set (
"Num Blocks", numBlocks_);
713 if (params->isParameter (
"Timer Label")) {
714 const std::string tempLabel = params->get (
"Timer Label", label_default_);
717 if (tempLabel != label_) {
719 params_->set (
"Timer Label", label_);
720 const std::string solveLabel =
721 label_ +
": PseudoBlockGmresSolMgr total solve time";
722 #ifdef BELOS_TEUCHOS_TIME_MONITOR 723 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
724 #endif // BELOS_TEUCHOS_TIME_MONITOR 725 if (ortho_ != Teuchos::null) {
726 ortho_->setLabel( label_ );
732 if (params->isParameter (
"Orthogonalization")) {
733 std::string tempOrthoType = params->get (
"Orthogonalization", orthoType_default_);
734 #ifdef HAVE_BELOS_TSQR 735 TEUCHOS_TEST_FOR_EXCEPTION(
736 tempOrthoType !=
"DGKS" && tempOrthoType !=
"ICGS" &&
737 tempOrthoType !=
"IMGS" && tempOrthoType !=
"TSQR",
738 std::invalid_argument,
739 "Belos::PseudoBlockGmresSolMgr::setParameters: " 740 "The \"Orthogonalization\" parameter must be one of \"DGKS\", \"ICGS\", " 741 "\"IMGS\", or \"TSQR\".");
743 TEUCHOS_TEST_FOR_EXCEPTION(
744 tempOrthoType !=
"DGKS" && tempOrthoType !=
"ICGS" &&
745 tempOrthoType !=
"IMGS",
746 std::invalid_argument,
747 "Belos::PseudoBlockGmresSolMgr::setParameters: " 748 "The \"Orthogonalization\" parameter must be one of \"DGKS\", \"ICGS\", " 750 #endif // HAVE_BELOS_TSQR 752 if (tempOrthoType != orthoType_) {
753 orthoType_ = tempOrthoType;
755 if (orthoType_ ==
"DGKS") {
757 if (orthoKappa_ <= 0) {
758 ortho_ = rcp (
new ortho_type (label_));
761 ortho_ = rcp (
new ortho_type (label_));
762 rcp_dynamic_cast<ortho_type> (ortho_)->setDepTol (orthoKappa_);
765 else if (orthoType_ ==
"ICGS") {
767 ortho_ = rcp (
new ortho_type (label_));
769 else if (orthoType_ ==
"IMGS") {
771 ortho_ = rcp (
new ortho_type (label_));
773 #ifdef HAVE_BELOS_TSQR 774 else if (orthoType_ ==
"TSQR") {
776 ortho_ = rcp (
new ortho_type (label_));
778 #endif // HAVE_BELOS_TSQR 783 if (params->isParameter (
"Orthogonalization Constant")) {
784 orthoKappa_ = params->get (
"Orthogonalization Constant", orthoKappa_default_);
787 params_->set (
"Orthogonalization Constant", orthoKappa_);
788 if (orthoType_ ==
"DGKS") {
789 if (orthoKappa_ > 0 && ! ortho_.is_null ()) {
791 rcp_dynamic_cast<ortho_type> (ortho_)->setDepTol (orthoKappa_);
797 if (params->isParameter (
"Verbosity")) {
798 if (Teuchos::isParameterType<int> (*params,
"Verbosity")) {
799 verbosity_ = params->get (
"Verbosity", verbosity_default_);
801 verbosity_ = (int) Teuchos::getParameter<Belos::MsgType> (*params,
"Verbosity");
805 params_->set (
"Verbosity", verbosity_);
806 if (! printer_.is_null ()) {
807 printer_->setVerbosity (verbosity_);
812 if (params->isParameter (
"Output Style")) {
813 if (Teuchos::isParameterType<int> (*params,
"Output Style")) {
814 outputStyle_ = params->get (
"Output Style", outputStyle_default_);
816 outputStyle_ = (int) Teuchos::getParameter<Belos::OutputType> (*params,
"Output Style");
820 params_->set (
"Output Style", verbosity_);
821 if (! outputTest_.is_null ()) {
827 if (params->isParameter (
"Output Stream")) {
828 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> > (*params,
"Output Stream");
831 params_->set(
"Output Stream", outputStream_);
832 if (! printer_.is_null ()) {
833 printer_->setOStream (outputStream_);
839 if (params->isParameter (
"Output Frequency")) {
840 outputFreq_ = params->get (
"Output Frequency", outputFreq_default_);
844 params_->set (
"Output Frequency", outputFreq_);
845 if (! outputTest_.is_null ()) {
846 outputTest_->setOutputFrequency (outputFreq_);
851 if (printer_.is_null ()) {
858 if (params->isParameter (
"Convergence Tolerance")) {
859 convtol_ = params->get (
"Convergence Tolerance", convtol_default_);
862 params_->set (
"Convergence Tolerance", convtol_);
863 if (! impConvTest_.is_null ()) {
864 impConvTest_->setTolerance (convtol_);
866 if (! expConvTest_.is_null ()) {
867 expConvTest_->setTolerance (convtol_);
872 bool userDefinedResidualScalingUpdated =
false;
873 if (params->isParameter (
"User Defined Residual Scaling")) {
874 const MagnitudeType tempResScaleFactor =
875 Teuchos::getParameter<MagnitudeType> (*params,
"User Defined Residual Scaling");
878 if (resScaleFactor_ != tempResScaleFactor) {
879 resScaleFactor_ = tempResScaleFactor;
880 userDefinedResidualScalingUpdated =
true;
883 if(userDefinedResidualScalingUpdated)
885 if (! params->isParameter (
"Implicit Residual Scaling") && ! impConvTest_.is_null ()) {
887 if(impResScale_ ==
"User Provided")
890 catch (std::exception& e) {
895 if (! params->isParameter (
"Explicit Residual Scaling") && ! expConvTest_.is_null ()) {
897 if(expResScale_ ==
"User Provided")
900 catch (std::exception& e) {
909 if (params->isParameter (
"Implicit Residual Scaling")) {
910 const std::string tempImpResScale =
911 Teuchos::getParameter<std::string> (*params,
"Implicit Residual Scaling");
914 if (impResScale_ != tempImpResScale) {
916 impResScale_ = tempImpResScale;
919 params_->set (
"Implicit Residual Scaling", impResScale_);
920 if (! impConvTest_.is_null ()) {
922 if(impResScale_ ==
"User Provided")
923 impConvTest_->defineScaleForm (impResScaleType,
Belos::TwoNorm, resScaleFactor_);
927 catch (std::exception& e) {
933 else if (userDefinedResidualScalingUpdated) {
936 if (! impConvTest_.is_null ()) {
938 if(impResScale_ ==
"User Provided")
939 impConvTest_->defineScaleForm (impResScaleType,
Belos::TwoNorm, resScaleFactor_);
941 catch (std::exception& e) {
949 if (params->isParameter (
"Explicit Residual Scaling")) {
950 const std::string tempExpResScale =
951 Teuchos::getParameter<std::string> (*params,
"Explicit Residual Scaling");
954 if (expResScale_ != tempExpResScale) {
956 expResScale_ = tempExpResScale;
959 params_->set (
"Explicit Residual Scaling", expResScale_);
960 if (! expConvTest_.is_null ()) {
962 if(expResScale_ ==
"User Provided")
963 expConvTest_->defineScaleForm (expResScaleType,
Belos::TwoNorm, resScaleFactor_);
967 catch (std::exception& e) {
973 else if (userDefinedResidualScalingUpdated) {
976 if (! expConvTest_.is_null ()) {
978 if(expResScale_ ==
"User Provided")
979 expConvTest_->defineScaleForm (expResScaleType,
Belos::TwoNorm, resScaleFactor_);
981 catch (std::exception& e) {
990 if (params->isParameter (
"Show Maximum Residual Norm Only")) {
991 showMaxResNormOnly_ =
992 Teuchos::getParameter<bool> (*params,
"Show Maximum Residual Norm Only");
995 params_->set (
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
996 if (! impConvTest_.is_null ()) {
997 impConvTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
999 if (! expConvTest_.is_null ()) {
1000 expConvTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
1007 if (params->isParameter(
"Deflation Quorum")) {
1008 defQuorum_ = params->get(
"Deflation Quorum", defQuorum_);
1009 TEUCHOS_TEST_FOR_EXCEPTION(
1010 defQuorum_ > blockSize_, std::invalid_argument,
1011 "Belos::PseudoBlockGmresSolMgr::setParameters: " 1012 "The \"Deflation Quorum\" parameter (= " << defQuorum_ <<
") must not be " 1013 "larger than \"Block Size\" (= " << blockSize_ <<
").");
1014 params_->set (
"Deflation Quorum", defQuorum_);
1015 if (! impConvTest_.is_null ()) {
1016 impConvTest_->setQuorum (defQuorum_);
1018 if (! expConvTest_.is_null ()) {
1019 expConvTest_->setQuorum (defQuorum_);
1024 if (ortho_.is_null ()) {
1025 if (orthoType_ ==
"DGKS") {
1027 if (orthoKappa_ <= 0) {
1028 ortho_ = rcp (
new ortho_type (label_));
1031 ortho_ = rcp (
new ortho_type (label_));
1032 rcp_dynamic_cast<ortho_type> (ortho_)->setDepTol (orthoKappa_);
1035 else if (orthoType_ ==
"ICGS") {
1037 ortho_ = rcp (
new ortho_type (label_));
1039 else if (orthoType_ ==
"IMGS") {
1041 ortho_ = rcp (
new ortho_type (label_));
1043 #ifdef HAVE_BELOS_TSQR 1044 else if (orthoType_ ==
"TSQR") {
1046 ortho_ = rcp (
new ortho_type (label_));
1048 #endif // HAVE_BELOS_TSQR 1050 #ifdef HAVE_BELOS_TSQR 1051 TEUCHOS_TEST_FOR_EXCEPTION(
1052 orthoType_ !=
"ICGS" && orthoType_ !=
"DGKS" &&
1053 orthoType_ !=
"IMGS" && orthoType_ !=
"TSQR",
1055 "Belos::PseudoBlockGmresSolMgr::setParameters(): " 1056 "Invalid orthogonalization type \"" << orthoType_ <<
"\".");
1058 TEUCHOS_TEST_FOR_EXCEPTION(
1059 orthoType_ !=
"ICGS" && orthoType_ !=
"DGKS" &&
1060 orthoType_ !=
"IMGS",
1062 "Belos::PseudoBlockGmresSolMgr::setParameters(): " 1063 "Invalid orthogonalization type \"" << orthoType_ <<
"\".");
1064 #endif // HAVE_BELOS_TSQR 1069 if (timerSolve_ == Teuchos::null) {
1070 std::string solveLabel = label_ +
": PseudoBlockGmresSolMgr total solve time";
1071 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1072 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
1081 template<
class ScalarType,
class MV,
class OP>
1087 userConvStatusTest_ = userConvStatusTest;
1090 template<
class ScalarType,
class MV,
class OP>
1096 debugStatusTest_ = debugStatusTest;
1101 template<
class ScalarType,
class MV,
class OP>
1102 Teuchos::RCP<const Teuchos::ParameterList>
1105 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
1106 if (is_null(validPL)) {
1107 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
1109 pl= Teuchos::rcp(
new Teuchos::ParameterList() );
1110 pl->set(
"Convergence Tolerance", convtol_default_,
1111 "The relative residual tolerance that needs to be achieved by the\n" 1112 "iterative solver in order for the linear system to be declared converged.");
1113 pl->set(
"Maximum Restarts", maxRestarts_default_,
1114 "The maximum number of restarts allowed for each\n" 1115 "set of RHS solved.");
1116 pl->set(
"Maximum Iterations", maxIters_default_,
1117 "The maximum number of block iterations allowed for each\n" 1118 "set of RHS solved.");
1119 pl->set(
"Num Blocks", numBlocks_default_,
1120 "The maximum number of vectors allowed in the Krylov subspace\n" 1121 "for each set of RHS solved.");
1122 pl->set(
"Block Size", blockSize_default_,
1123 "The number of RHS solved simultaneously.");
1124 pl->set(
"Verbosity", verbosity_default_,
1125 "What type(s) of solver information should be outputted\n" 1126 "to the output stream.");
1127 pl->set(
"Output Style", outputStyle_default_,
1128 "What style is used for the solver information outputted\n" 1129 "to the output stream.");
1130 pl->set(
"Output Frequency", outputFreq_default_,
1131 "How often convergence information should be outputted\n" 1132 "to the output stream.");
1133 pl->set(
"Deflation Quorum", defQuorum_default_,
1134 "The number of linear systems that need to converge before\n" 1135 "they are deflated. This number should be <= block size.");
1136 pl->set(
"Output Stream", outputStream_default_,
1137 "A reference-counted pointer to the output stream where all\n" 1138 "solver output is sent.");
1139 pl->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
1140 "When convergence information is printed, only show the maximum\n" 1141 "relative residual norm when the block size is greater than one.");
1142 pl->set(
"Implicit Residual Scaling", impResScale_default_,
1143 "The type of scaling used in the implicit residual convergence test.");
1144 pl->set(
"Explicit Residual Scaling", expResScale_default_,
1145 "The type of scaling used in the explicit residual convergence test.");
1146 pl->set(
"Timer Label", label_default_,
1147 "The string to use as a prefix for the timer labels.");
1149 pl->set(
"Orthogonalization", orthoType_default_,
1150 "The type of orthogonalization to use: DGKS, ICGS, IMGS.");
1151 pl->set(
"Orthogonalization Constant",orthoKappa_default_,
1152 "The constant used by DGKS orthogonalization to determine\n" 1153 "whether another step of classical Gram-Schmidt is necessary.");
1160 template<
class ScalarType,
class MV,
class OP>
1172 if ( !Teuchos::is_null(problem_->getLeftPrec()) ) {
1179 Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
1180 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
1181 if(impResScale_ ==
"User Provided")
1186 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1187 impConvTest_ = tmpImpConvTest;
1190 Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
1191 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_, defQuorum_ ) );
1192 tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit,
Belos::TwoNorm );
1193 if(expResScale_ ==
"User Provided")
1197 tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1198 expConvTest_ = tmpExpConvTest;
1201 convTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
1207 Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
1208 Teuchos::rcp(
new StatusTestImpResNorm_t( convtol_, defQuorum_ ) );
1209 if(impResScale_ ==
"User Provided")
1213 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
1214 impConvTest_ = tmpImpConvTest;
1217 expConvTest_ = impConvTest_;
1218 convTest_ = impConvTest_;
1221 if (nonnull(debugStatusTest_) ) {
1223 convTest_ = Teuchos::rcp(
1224 new StatusTestCombo_t( StatusTestCombo_t::AND, convTest_, debugStatusTest_ ) );
1227 if (nonnull(userConvStatusTest_) ) {
1229 convTest_ = Teuchos::rcp(
1230 new StatusTestCombo_t( StatusTestCombo_t::SEQ, convTest_, userConvStatusTest_ ) );
1235 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
1243 std::string solverDesc =
" Pseudo Block Gmres ";
1244 outputTest_->setSolverDesc( solverDesc );
1255 template<
class ScalarType,
class MV,
class OP>
1263 Teuchos::BLAS<int,ScalarType> blas;
1266 "Belos::PseudoBlockGmresSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1269 if (!isSTSet_ || (!expResTest_ && !Teuchos::is_null(problem_->getLeftPrec())) ) {
1271 "Belos::BlockGmresSolMgr::solve(): Linear problem and requested status tests are incompatible.");
1277 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1279 std::vector<int> currIdx( numCurrRHS );
1280 blockSize_ = numCurrRHS;
1281 for (
int i=0; i<numCurrRHS; ++i)
1282 { currIdx[i] = startPtr+i; }
1285 problem_->setLSIndex( currIdx );
1289 Teuchos::ParameterList plist;
1290 plist.set(
"Num Blocks",numBlocks_);
1293 outputTest_->reset();
1294 loaDetected_ =
false;
1297 bool isConverged =
true;
1302 Teuchos::RCP<PseudoBlockGmresIter<ScalarType,MV,OP> > block_gmres_iter
1307 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1308 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1311 while ( numRHS2Solve > 0 ) {
1314 std::vector<int> convRHSIdx;
1315 std::vector<int> currRHSIdx( currIdx );
1316 currRHSIdx.resize(numCurrRHS);
1319 block_gmres_iter->setNumBlocks( numBlocks_ );
1322 block_gmres_iter->resetNumIters();
1325 outputTest_->resetNumCalls();
1331 std::vector<int> index(1);
1332 Teuchos::RCP<MV> tmpV, R_0 =
MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1333 newState.
V.resize( blockSize_ );
1334 newState.
Z.resize( blockSize_ );
1335 for (
int i=0; i<blockSize_; ++i) {
1340 Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > tmpZ
1341 = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>( 1 ));
1344 int rank = ortho_->normalize( *tmpV, tmpZ );
1346 "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
1348 newState.
V[i] = tmpV;
1349 newState.
Z[i] = tmpZ;
1353 block_gmres_iter->initialize(newState);
1354 int numRestarts = 0;
1360 block_gmres_iter->iterate();
1367 if ( convTest_->getStatus() ==
Passed ) {
1369 if ( expConvTest_->getLOADetected() ) {
1379 loaDetected_ =
true;
1381 "Belos::PseudoBlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
1382 isConverged =
false;
1386 std::vector<int> convIdx = expConvTest_->convIndices();
1390 if (convIdx.size() == currRHSIdx.size())
1397 problem_->setCurrLS();
1401 int curDim = oldState.
curDim;
1406 std::vector<int> oldRHSIdx( currRHSIdx );
1407 std::vector<int> defRHSIdx;
1408 for (
unsigned int i=0; i<currRHSIdx.size(); ++i) {
1410 for (
unsigned int j=0; j<convIdx.size(); ++j) {
1411 if (currRHSIdx[i] == convIdx[j]) {
1417 defRHSIdx.push_back( i );
1420 defState.
V.push_back( Teuchos::rcp_const_cast<MV>( oldState.
V[i] ) );
1421 defState.
Z.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.
Z[i] ) );
1422 defState.
H.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseMatrix<int,ScalarType> >( oldState.
H[i] ) );
1423 defState.
sn.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,ScalarType> >( oldState.
sn[i] ) );
1424 defState.
cs.push_back( Teuchos::rcp_const_cast<Teuchos::SerialDenseVector<int,MagnitudeType> >(oldState.
cs[i] ) );
1425 currRHSIdx[have] = currRHSIdx[i];
1429 defRHSIdx.resize(currRHSIdx.size()-have);
1430 currRHSIdx.resize(have);
1434 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1438 problem_->setLSIndex( convIdx );
1441 problem_->updateSolution( defUpdate,
true );
1445 problem_->setLSIndex( currRHSIdx );
1448 defState.
curDim = curDim;
1451 block_gmres_iter->initialize(defState);
1458 else if ( maxIterTest_->getStatus() ==
Passed ) {
1460 isConverged =
false;
1468 else if ( block_gmres_iter->getCurSubspaceDim() == block_gmres_iter->getMaxSubspaceDim() ) {
1470 if ( numRestarts >= maxRestarts_ ) {
1471 isConverged =
false;
1476 printer_->stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
1479 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1480 problem_->updateSolution( update,
true );
1487 newstate.
V.resize(currRHSIdx.size());
1488 newstate.
Z.resize(currRHSIdx.size());
1492 R_0 =
MVT::Clone( *(problem_->getInitPrecResVec()), currRHSIdx.size() );
1493 problem_->computeCurrPrecResVec( &*R_0 );
1494 for (
unsigned int i=0; i<currRHSIdx.size(); ++i) {
1500 Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > tmpZ
1501 = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>( 1 ));
1504 int rank = ortho_->normalize( *tmpV, tmpZ );
1506 "Belos::PseudoBlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors after the restart.");
1508 newstate.
V[i] = tmpV;
1509 newstate.
Z[i] = tmpZ;
1514 block_gmres_iter->initialize(newstate);
1526 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
1527 "Belos::PseudoBlockGmresSolMgr::solve(): Invalid return from PseudoBlockGmresIter::iterate().");
1533 block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
1536 sTest_->checkStatus( &*block_gmres_iter );
1537 if (convTest_->getStatus() !=
Passed)
1538 isConverged =
false;
1541 catch (
const std::exception &e) {
1542 printer_->stream(
Errors) <<
"Error! Caught std::exception in PseudoBlockGmresIter::iterate() at iteration " 1543 << block_gmres_iter->getNumIters() << std::endl
1544 << e.what() << std::endl;
1551 if (nonnull(userConvStatusTest_)) {
1553 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1554 problem_->updateSolution( update,
true );
1556 else if (nonnull(expConvTest_->getSolution())) {
1558 Teuchos::RCP<MV> newX = expConvTest_->getSolution();
1559 Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1564 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1565 problem_->updateSolution( update,
true );
1569 problem_->setCurrLS();
1572 startPtr += numCurrRHS;
1573 numRHS2Solve -= numCurrRHS;
1574 if ( numRHS2Solve > 0 ) {
1575 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1577 blockSize_ = numCurrRHS;
1578 currIdx.resize( numCurrRHS );
1579 for (
int i=0; i<numCurrRHS; ++i)
1580 { currIdx[i] = startPtr+i; }
1583 if (defQuorum_ > blockSize_) {
1584 if (impConvTest_ != Teuchos::null)
1585 impConvTest_->setQuorum( blockSize_ );
1586 if (expConvTest_ != Teuchos::null)
1587 expConvTest_->setQuorum( blockSize_ );
1591 problem_->setLSIndex( currIdx );
1594 currIdx.resize( numRHS2Solve );
1605 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1610 Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1614 numIters_ = maxIterTest_->getNumIters();
1627 const std::vector<MagnitudeType>* pTestValues = NULL;
1629 pTestValues = expConvTest_->getTestValue();
1630 if (pTestValues == NULL || pTestValues->size() < 1) {
1631 pTestValues = impConvTest_->getTestValue();
1636 pTestValues = impConvTest_->getTestValue();
1638 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1639 "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's " 1640 "getTestValue() method returned NULL. Please report this bug to the " 1641 "Belos developers.");
1642 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1643 "Belos::PseudoBlockGmresSolMgr::solve(): The implicit convergence test's " 1644 "getTestValue() method returned a vector of length zero. Please report " 1645 "this bug to the Belos developers.");
1650 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1653 if (!isConverged || loaDetected_) {
1660 template<
class ScalarType,
class MV,
class OP>
1663 std::ostringstream out;
1665 out <<
"\"Belos::PseudoBlockGmresSolMgr\": {";
1666 if (this->getObjectLabel () !=
"") {
1667 out <<
"Label: " << this->getObjectLabel () <<
", ";
1669 out <<
"Num Blocks: " << numBlocks_
1670 <<
", Maximum Iterations: " << maxIters_
1671 <<
", Maximum Restarts: " << maxRestarts_
1672 <<
", Convergence Tolerance: " << convtol_
1678 template<
class ScalarType,
class MV,
class OP>
1682 const Teuchos::EVerbosityLevel verbLevel)
const 1684 using Teuchos::TypeNameTraits;
1685 using Teuchos::VERB_DEFAULT;
1686 using Teuchos::VERB_NONE;
1687 using Teuchos::VERB_LOW;
1694 const Teuchos::EVerbosityLevel vl =
1695 (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
1697 if (vl != VERB_NONE) {
1698 Teuchos::OSTab tab0 (out);
1700 out <<
"\"Belos::PseudoBlockGmresSolMgr\":" << endl;
1701 Teuchos::OSTab tab1 (out);
1702 out <<
"Template parameters:" << endl;
1704 Teuchos::OSTab tab2 (out);
1705 out <<
"ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1706 <<
"MV: " << TypeNameTraits<MV>::name () << endl
1707 <<
"OP: " << TypeNameTraits<OP>::name () << endl;
1709 if (this->getObjectLabel () !=
"") {
1710 out <<
"Label: " << this->getObjectLabel () << endl;
1712 out <<
"Num Blocks: " << numBlocks_ << endl
1713 <<
"Maximum Iterations: " << maxIters_ << endl
1714 <<
"Maximum Restarts: " << maxRestarts_ << endl
1715 <<
"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.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
The current parameters for this solver.
virtual ~PseudoBlockGmresSolMgr()
Destructor.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Class which manages the output and verbosity of the Belos solvers.
virtual void setUserConvStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &userConvStatusTest)
Set a custom status test.
PseudoBlockGmresSolMgr()
Empty constructor.
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Set the parameters the solver manager should use to solve the linear problem.
PseudoBlockGmresSolMgrLinearProblemFailure(const std::string &what_arg)
ScaleType
The type of scaling to use on the residual norm value.
Convergence test using the implicit residual norm(s), with an explicit residual norm(s) check for los...
PseudoBlockGmresSolMgrOrthoFailure(const std::string &what_arg)
A factory class for generating StatusTestOutput objects.
std::vector< Teuchos::RCP< const MV > > V
The current Krylov basis.
An implementation of StatusTestResNorm using a family of residual norms.
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Return a reference to the linear problem being solved by this solver manager.
Belos::StatusTest class for specifying a maximum number of iterations.
A pure virtual class for defining the status tests for the Belos iterative solvers.
Structure to contain pointers to PseudoBlockGmresIter state variables.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
A factory class for generating StatusTestOutput objects.
Iterated Modified Gram-Schmidt (IMGS) implementation of the Belos::OrthoManager class.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
A Belos::StatusTest class for specifying a maximum number of iterations.
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
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
Print the object with the given verbosity level to a FancyOStream.
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.
Class which describes the linear problem to be solved by the iterative solver.
Interface to standard and "pseudoblock" GMRES.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
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 ...
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
bool isLOADetected() const
Whether a "loss of accuracy" was detected during the last solve().
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
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.
Belos concrete class for performing the pseudo-block GMRES iteration.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest for specifying an implicit residual norm stopping criteria that checks for loss of ...
virtual void setDebugStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &debugStatusTest)
Set a debug status test.
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.
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.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem to solve.
Parent class to all Belos exceptions.
std::vector< Teuchos::RCP< const Teuchos::SerialDenseVector< int, ScalarType > > > sn
The current Given's rotation coefficients.
ReturnType solve()
This method performs possibly repeated calls to the underlying linear solver's iterate() routine unti...
MatOrthoManager subclass using TSQR or DGKS.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
A list of valid default parameters for this solver.
int getNumIters() const
Iteration count for the most recent call to solve().
Belos header file which uses auto-configuration information to include necessary C++ headers...
std::string description() const
Return a one-line description of this object.
This class implements the pseudo-block GMRES iteration, where a block Krylov subspace is constructed ...
PseudoBlockGmresIterOrthoFailure is thrown when the orthogonalization manager is unable to generate o...
MagnitudeType achievedTol() const
Tolerance achieved by the last solve() invocation.
PseudoBlockGmresSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate...