42 #ifndef BELOS_GMRES_POLY_SOLMGR_HPP 43 #define BELOS_GMRES_POLY_SOLMGR_HPP 67 #include "Teuchos_BLAS.hpp" 68 #include "Teuchos_as.hpp" 69 #ifdef BELOS_TEUCHOS_TIME_MONITOR 70 #include "Teuchos_TimeMonitor.hpp" 154 template<
class ScalarType,
class MV,
class OP>
159 typedef Teuchos::ScalarTraits<ScalarType>
STS;
160 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
MagnitudeType;
161 typedef Teuchos::ScalarTraits<MagnitudeType>
MT;
189 const Teuchos::RCP<Teuchos::ParameterList> &pl );
195 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
222 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
245 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms )
override;
308 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
problem_;
315 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
sTest_;
322 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> >
ortho_;
355 Teuchos::RCP<Teuchos::SerialDenseVector<int, ScalarType> >
poly_r0_;
356 Teuchos::RCP<Belos::GmresPolyOp<ScalarType, MV, OP> >
poly_Op_;
368 mutable Teuchos::RCP<const Teuchos::ParameterList>
validPL_;
372 template<
class ScalarType,
class MV,
class OP>
374 outputStream_ (outputStream_default_),
378 maxDegree_ (maxDegree_default_),
379 maxRestarts_ (maxRestarts_default_),
380 maxIters_ (maxIters_default_),
382 blockSize_ (blockSize_default_),
383 numBlocks_ (numBlocks_default_),
384 verbosity_ (verbosity_default_),
385 outputStyle_ (outputStyle_default_),
386 outputFreq_ (outputFreq_default_),
387 strictConvTol_ (strictConvTol_default_),
388 showMaxResNormOnly_ (showMaxResNormOnly_default_),
389 orthoType_ (orthoType_default_),
390 impResScale_ (impResScale_default_),
391 expResScale_ (expResScale_default_),
393 label_ (label_default_),
394 isPolyBuilt_ (false),
402 template<
class ScalarType,
class MV,
class OP>
405 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
407 outputStream_ (outputStream_default_),
411 maxDegree_ (maxDegree_default_),
412 maxRestarts_ (maxRestarts_default_),
413 maxIters_ (maxIters_default_),
415 blockSize_ (blockSize_default_),
416 numBlocks_ (numBlocks_default_),
417 verbosity_ (verbosity_default_),
418 outputStyle_ (outputStyle_default_),
419 outputFreq_ (outputFreq_default_),
420 strictConvTol_ (strictConvTol_default_),
421 showMaxResNormOnly_ (showMaxResNormOnly_default_),
422 orthoType_ (orthoType_default_),
423 impResScale_ (impResScale_default_),
424 expResScale_ (expResScale_default_),
426 label_ (label_default_),
427 isPolyBuilt_ (false),
433 TEUCHOS_TEST_FOR_EXCEPTION(
434 problem_.is_null (), std::invalid_argument,
435 "Belos::GmresPolySolMgr: The given linear problem is null. " 436 "Please call this constructor with a nonnull LinearProblem argument, " 437 "or call the constructor that does not take a LinearProblem.");
441 if (! pl.is_null ()) {
447 template<
class ScalarType,
class MV,
class OP>
448 Teuchos::RCP<const Teuchos::ParameterList>
451 if (validPL_.is_null ()) {
452 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList ();
457 "The relative residual tolerance that used to construct the GMRES polynomial.");
458 pl->set(
"Maximum Degree", static_cast<int>(maxDegree_default_),
459 "The maximum degree allowed for any GMRES polynomial.");
461 "The relative residual tolerance that needs to be achieved by the\n" 462 "iterative solver in order for the linear system to be declared converged." );
463 pl->set(
"Maximum Restarts", static_cast<int>(maxRestarts_default_),
464 "The maximum number of restarts allowed for each\n" 465 "set of RHS solved.");
466 pl->set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
467 "The maximum number of block iterations allowed for each\n" 468 "set of RHS solved.");
469 pl->set(
"Num Blocks", static_cast<int>(numBlocks_default_),
470 "The maximum number of blocks allowed in the Krylov subspace\n" 471 "for each set of RHS solved.");
472 pl->set(
"Block Size", static_cast<int>(blockSize_default_),
473 "The number of vectors in each block. This number times the\n" 474 "number of blocks is the total Krylov subspace dimension.");
475 pl->set(
"Verbosity", static_cast<int>(verbosity_default_),
476 "What type(s) of solver information should be outputted\n" 477 "to the output stream.");
478 pl->set(
"Output Style", static_cast<int>(outputStyle_default_),
479 "What style is used for the solver information outputted\n" 480 "to the output stream.");
481 pl->set(
"Output Frequency", static_cast<int>(outputFreq_default_),
482 "How often convergence information should be outputted\n" 483 "to the output stream.");
484 pl->set(
"Output Stream", Teuchos::rcp(outputStream_default_,
false),
485 "A reference-counted pointer to the output stream where all\n" 486 "solver output is sent.");
487 pl->set(
"Strict Convergence", static_cast<bool>(strictConvTol_default_),
488 "After polynomial is applied, whether solver should try to achieve\n" 489 "the relative residual tolerance.");
490 pl->set(
"Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
491 "When convergence information is printed, only show the maximum\n" 492 "relative residual norm when the block size is greater than one.");
493 pl->set(
"Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
494 "The type of scaling used in the implicit residual convergence test.");
495 pl->set(
"Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
496 "The type of scaling used in the explicit residual convergence test.");
497 pl->set(
"Timer Label", static_cast<const char *>(label_default_),
498 "The string to use as a prefix for the timer labels.");
499 pl->set(
"Orthogonalization", static_cast<const char *>(orthoType_default_),
500 "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
502 "The constant used by DGKS orthogonalization to determine\n" 503 "whether another step of classical Gram-Schmidt is necessary.");
510 template<
class ScalarType,
class MV,
class OP>
515 if (params_.is_null ()) {
516 params_ = Teuchos::parameterList (*getValidParameters ());
519 params->validateParameters (*getValidParameters ());
523 if (params->isParameter(
"Maximum Degree")) {
524 maxDegree_ = params->get(
"Maximum Degree",maxDegree_default_);
527 params_->set(
"Maximum Degree", maxDegree_);
531 if (params->isParameter(
"Maximum Restarts")) {
532 maxRestarts_ = params->get(
"Maximum Restarts",maxRestarts_default_);
535 params_->set(
"Maximum Restarts", maxRestarts_);
539 if (params->isParameter(
"Maximum Iterations")) {
540 maxIters_ = params->get(
"Maximum Iterations",maxIters_default_);
543 params_->set(
"Maximum Iterations", maxIters_);
544 if (maxIterTest_!=Teuchos::null)
545 maxIterTest_->setMaxIters( maxIters_ );
549 if (params->isParameter(
"Block Size")) {
550 blockSize_ = params->get(
"Block Size",blockSize_default_);
551 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
552 "Belos::GmresPolySolMgr: \"Block Size\" must be strictly positive.");
555 params_->set(
"Block Size", blockSize_);
559 if (params->isParameter(
"Num Blocks")) {
560 numBlocks_ = params->get(
"Num Blocks",numBlocks_default_);
561 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
562 "Belos::GmresPolySolMgr: \"Num Blocks\" must be strictly positive.");
565 params_->set(
"Num Blocks", numBlocks_);
569 if (params->isParameter(
"Timer Label")) {
570 std::string tempLabel = params->get(
"Timer Label", label_default_);
573 if (tempLabel != label_) {
575 params_->set(
"Timer Label", label_);
576 std::string solveLabel = label_ +
": GmresPolySolMgr total solve time";
577 #ifdef BELOS_TEUCHOS_TIME_MONITOR 578 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
580 std::string polyLabel = label_ +
": GmresPolySolMgr polynomial creation time";
581 #ifdef BELOS_TEUCHOS_TIME_MONITOR 582 timerPoly_ = Teuchos::TimeMonitor::getNewCounter(polyLabel);
584 if (ortho_ != Teuchos::null) {
585 ortho_->setLabel( label_ );
591 if (params->isParameter(
"Orthogonalization")) {
592 std::string tempOrthoType = params->get(
"Orthogonalization",orthoType_default_);
593 TEUCHOS_TEST_FOR_EXCEPTION( tempOrthoType !=
"DGKS" && tempOrthoType !=
"ICGS" && tempOrthoType !=
"IMGS",
594 std::invalid_argument,
595 "Belos::GmresPolySolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
596 if (tempOrthoType != orthoType_) {
597 params_->set(
"Orthogonalization", orthoType_);
598 orthoType_ = tempOrthoType;
600 if (orthoType_==
"DGKS") {
601 if (orthoKappa_ <= 0) {
609 else if (orthoType_==
"ICGS") {
612 else if (orthoType_==
"IMGS") {
619 if (params->isParameter(
"Orthogonalization Constant")) {
620 if (params->isType<
MagnitudeType> (
"Orthogonalization Constant")) {
621 orthoKappa_ = params->get (
"Orthogonalization Constant",
625 orthoKappa_ = params->get (
"Orthogonalization Constant",
630 params_->set(
"Orthogonalization Constant",orthoKappa_);
631 if (orthoType_==
"DGKS") {
632 if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
639 if (params->isParameter(
"Verbosity")) {
640 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
641 verbosity_ = params->get(
"Verbosity", verbosity_default_);
643 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
647 params_->set(
"Verbosity", verbosity_);
648 if (printer_ != Teuchos::null)
649 printer_->setVerbosity(verbosity_);
653 if (params->isParameter(
"Output Style")) {
654 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
655 outputStyle_ = params->get(
"Output Style", outputStyle_default_);
657 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
661 params_->set(
"Output Style", outputStyle_);
662 if (outputTest_ != Teuchos::null) {
668 if (params->isParameter(
"Output Stream")) {
669 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
672 params_->set(
"Output Stream", outputStream_);
673 if (printer_ != Teuchos::null)
674 printer_->setOStream( outputStream_ );
679 if (params->isParameter(
"Output Frequency")) {
680 outputFreq_ = params->get(
"Output Frequency", outputFreq_default_);
684 params_->set(
"Output Frequency", outputFreq_);
685 if (outputTest_ != Teuchos::null)
686 outputTest_->setOutputFrequency( outputFreq_ );
690 if (printer_ == Teuchos::null) {
699 if (params->isParameter(
"Polynomial Tolerance")) {
700 if (params->isType<
MagnitudeType> (
"Convergence Tolerance")) {
701 polytol_ = params->get (
"Polynomial Tolerance",
709 params_->set(
"Polynomial Tolerance", polytol_);
713 if (params->isParameter(
"Convergence Tolerance")) {
714 if (params->isType<
MagnitudeType> (
"Convergence Tolerance")) {
715 convtol_ = params->get (
"Convergence Tolerance",
723 params_->set(
"Convergence Tolerance", convtol_);
724 if (impConvTest_ != Teuchos::null)
725 impConvTest_->setTolerance( convtol_ );
726 if (expConvTest_ != Teuchos::null)
727 expConvTest_->setTolerance( convtol_ );
731 if (params->isParameter(
"Strict Convergence")) {
732 strictConvTol_ = params->get(
"Strict Convergence",strictConvTol_default_);
735 params_->set(
"Strict Convergence", strictConvTol_);
739 if (params->isParameter(
"Implicit Residual Scaling")) {
740 std::string tempImpResScale = Teuchos::getParameter<std::string>( *params,
"Implicit Residual Scaling" );
743 if (impResScale_ != tempImpResScale) {
745 impResScale_ = tempImpResScale;
748 params_->set(
"Implicit Residual Scaling", impResScale_);
749 if (impConvTest_ != Teuchos::null) {
753 catch (std::exception& e) {
761 if (params->isParameter(
"Explicit Residual Scaling")) {
762 std::string tempExpResScale = Teuchos::getParameter<std::string>( *params,
"Explicit Residual Scaling" );
765 if (expResScale_ != tempExpResScale) {
767 expResScale_ = tempExpResScale;
770 params_->set(
"Explicit Residual Scaling", expResScale_);
771 if (expConvTest_ != Teuchos::null) {
775 catch (std::exception& e) {
784 if (params->isParameter(
"Show Maximum Residual Norm Only")) {
785 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,
"Show Maximum Residual Norm Only");
788 params_->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
789 if (impConvTest_ != Teuchos::null)
790 impConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
791 if (expConvTest_ != Teuchos::null)
792 expConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
796 if (ortho_ == Teuchos::null) {
797 params_->set(
"Orthogonalization", orthoType_);
798 if (orthoType_==
"DGKS") {
799 if (orthoKappa_ <= 0) {
807 else if (orthoType_==
"ICGS") {
810 else if (orthoType_==
"IMGS") {
814 TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!=
"ICGS"&&orthoType_!=
"DGKS"&&orthoType_!=
"IMGS",std::logic_error,
815 "Belos::GmresPolySolMgr(): Invalid orthogonalization type.");
820 if (timerSolve_ == Teuchos::null) {
821 std::string solveLabel = label_ +
": GmresPolySolMgr total solve time";
822 #ifdef BELOS_TEUCHOS_TIME_MONITOR 823 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
827 if (timerPoly_ == Teuchos::null) {
828 std::string polyLabel = label_ +
": GmresPolySolMgr polynomial creation time";
829 #ifdef BELOS_TEUCHOS_TIME_MONITOR 830 timerPoly_ = Teuchos::TimeMonitor::getNewCounter(polyLabel);
839 template<
class ScalarType,
class MV,
class OP>
851 if (!Teuchos::is_null(problem_->getLeftPrec())) {
858 Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
859 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_ ) );
861 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
862 impConvTest_ = tmpImpConvTest;
865 Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
866 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_ ) );
867 tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit,
Belos::TwoNorm );
869 tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
870 expConvTest_ = tmpExpConvTest;
873 convTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
879 Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
880 Teuchos::rcp(
new StatusTestImpResNorm_t( convtol_ ) );
882 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
883 impConvTest_ = tmpImpConvTest;
886 expConvTest_ = impConvTest_;
887 convTest_ = impConvTest_;
890 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
898 std::string solverDesc =
" Gmres Polynomial ";
899 outputTest_->setSolverDesc( solverDesc );
908 template<
class ScalarType,
class MV,
class OP>
912 Teuchos::RCP<MV> newX = MVT::Clone( *(problem_->getLHS()), 1 );
913 Teuchos::RCP<MV> newB = MVT::Clone( *(problem_->getRHS()), 1 );
914 MVT::MvInit( *newX, STS::zero() );
915 MVT::MvRandom( *newB );
916 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > newProblem =
918 newProblem->setLeftPrec( problem_->getLeftPrec() );
919 newProblem->setRightPrec( problem_->getRightPrec() );
920 newProblem->setLabel(
"Belos GMRES Poly Generation");
921 newProblem->setProblem();
922 std::vector<int> idx(1,0);
923 newProblem->setLSIndex( idx );
926 Teuchos::ParameterList polyList;
929 polyList.set(
"Num Blocks",maxDegree_);
930 polyList.set(
"Block Size",1);
931 polyList.set(
"Keep Hessenberg",
true);
934 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxItrTst =
938 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTst =
943 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > polyTest =
947 Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > gmres_iter;
951 Teuchos::RCP<MV> V_0 = MVT::Clone( *(newProblem->getRHS()), 1 );
952 newProblem->computeCurrPrecResVec( &*V_0 );
955 poly_r0_ = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>(1) );
958 int rank = ortho_->normalize( *V_0, poly_r0_ );
960 "Belos::GmresPolySolMgr::generatePoly(): Failed to compute initial block of orthonormal vectors for polynomial generation.");
965 newstate.
z = poly_r0_;
967 gmres_iter->initializeGmres(newstate);
972 gmres_iter->iterate();
975 if ( convTst->getStatus() ==
Passed ) {
982 gmres_iter->updateLSQR( gmres_iter->getCurSubspaceDim() );
985 polyTest->checkStatus( &*gmres_iter );
986 if (convTst->getStatus() ==
Passed) {
990 catch (std::exception e) {
992 printer_->stream(
Errors) <<
"Error! Caught exception in " 993 "BlockGmresIter::iterate() at iteration " << gmres_iter->getNumIters()
994 << endl << e.what () << endl;
1004 Teuchos::RCP<MV> currX = gmres_iter->getCurrentUpdate();
1010 poly_dim_ = gmresState.
curDim;
1011 if (poly_dim_ == 0) {
1017 poly_y_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *gmresState.
z, poly_dim_, 1 ) );
1018 poly_H_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( *gmresState.
H ) );
1022 const ScalarType one = STS::one ();
1023 Teuchos::BLAS<int,ScalarType> blas;
1024 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
1025 Teuchos::NON_UNIT_DIAG, poly_dim_, 1, one,
1026 gmresState.
R->values(), gmresState.
R->stride(),
1027 poly_y_->values(), poly_y_->stride() );
1031 poly_Op_ = Teuchos::rcp(
1038 template<
class ScalarType,
class MV,
class OP>
1043 typedef Teuchos::SerialDenseMatrix<int, ScalarType> SDM;
1050 setParameters (Teuchos::parameterList (*getValidParameters ()));
1053 TEUCHOS_TEST_FOR_EXCEPTION(
1055 "Belos::GmresPolySolMgr::solve: The linear problem has not been set yet, " 1056 "or was set to null. Please call setProblem() with a nonnull input before " 1057 "calling solve().");
1059 TEUCHOS_TEST_FOR_EXCEPTION(
1061 "Belos::GmresPolySolMgr::solve: The linear problem is not ready. Please " 1062 "call setProblem() on the LinearProblem object before calling solve().");
1064 if (! isSTSet_ || (! expResTest_ && ! problem_->getLeftPrec ().is_null ())) {
1069 const bool check = checkStatusTest ();
1070 TEUCHOS_TEST_FOR_EXCEPTION(
1072 "Belos::GmresPolySolMgr::solve: Linear problem and requested status " 1073 "tests are incompatible.");
1078 if (! isPolyBuilt_) {
1079 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1080 Teuchos::TimeMonitor slvtimer(*timerPoly_);
1082 isPolyBuilt_ = generatePoly();
1084 "Belos::GmresPolySolMgr::generatePoly(): Failed to generate a non-trivial polynomial, reduce polynomial tolerance.");
1086 "Belos::GmresPolySolMgr::generatePoly(): Failed to generate polynomial that satisfied requirements.");
1090 bool isConverged =
true;
1094 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1095 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1099 poly_Op_->Apply( *problem_->getRHS(), *problem_->getLHS() );
1102 problem_->setProblem ();
1105 if (strictConvTol_) {
1109 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1110 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1118 std::vector<int> currIdx (blockSize_);
1119 for (
int i = 0; i < numCurrRHS; ++i) {
1120 currIdx[i] = startPtr+i;
1122 for (
int i = numCurrRHS; i < blockSize_; ++i) {
1127 problem_->setLSIndex (currIdx);
1131 Teuchos::ParameterList plist;
1132 plist.set (
"Block Size", blockSize_);
1134 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1135 if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {
1136 ptrdiff_t tmpNumBlocks = 0;
1137 if (blockSize_ == 1) {
1138 tmpNumBlocks = dim / blockSize_;
1140 tmpNumBlocks = ( dim - blockSize_) / blockSize_;
1143 <<
"Warning! Requested Krylov subspace dimension is larger than " 1144 <<
"operator dimension!" << std::endl <<
"The maximum number of " 1145 <<
"blocks allowed for the Krylov subspace will be adjusted to " 1146 << tmpNumBlocks << std::endl;
1147 plist.set (
"Num Blocks", Teuchos::asSafe<int> (tmpNumBlocks));
1150 plist.set (
"Num Blocks", numBlocks_);
1154 outputTest_->reset ();
1155 loaDetected_ =
false;
1165 RCP<GmresIteration<ScalarType,MV,OP> > block_gmres_iter =
1167 outputTest_, ortho_, plist));
1170 while (numRHS2Solve > 0) {
1173 if (blockSize_*numBlocks_ > dim) {
1174 int tmpNumBlocks = 0;
1175 if (blockSize_ == 1) {
1177 tmpNumBlocks = dim / blockSize_;
1180 tmpNumBlocks = (dim - blockSize_) / blockSize_;
1182 block_gmres_iter->setSize (blockSize_, tmpNumBlocks);
1185 block_gmres_iter->setSize (blockSize_, numBlocks_);
1189 block_gmres_iter->resetNumIters ();
1192 outputTest_->resetNumCalls ();
1195 RCP<MV> V_0 = MVT::CloneCopy (*(problem_->getInitPrecResVec ()), currIdx);
1198 RCP<SDM> z_0 = rcp (
new SDM (blockSize_, blockSize_));
1201 int rank = ortho_->normalize (*V_0, z_0);
1202 TEUCHOS_TEST_FOR_EXCEPTION(
1204 "Belos::GmresPolySolMgr::solve: Failed to compute initial block of " 1205 "orthonormal vectors.");
1212 block_gmres_iter->initializeGmres(newstate);
1213 int numRestarts = 0;
1217 block_gmres_iter->iterate();
1220 if ( convTest_->getStatus() ==
Passed ) {
1221 if ( expConvTest_->getLOADetected() ) {
1223 loaDetected_ =
true;
1224 isConverged =
false;
1230 else if ( maxIterTest_->getStatus() ==
Passed ) {
1232 isConverged =
false;
1236 else if (block_gmres_iter->getCurSubspaceDim () ==
1237 block_gmres_iter->getMaxSubspaceDim ()) {
1238 if (numRestarts >= maxRestarts_) {
1239 isConverged =
false;
1244 printer_->stream(
Debug)
1245 <<
" Performing restart number " << numRestarts <<
" of " 1246 << maxRestarts_ << std::endl << std::endl;
1249 RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1250 problem_->updateSolution (update,
true);
1260 RCP<MV> VV_0 = MVT::Clone (*(oldState.
V), blockSize_);
1261 problem_->computeCurrPrecResVec (&*VV_0);
1267 RCP<SDM> zz_0 = rcp (
new SDM (blockSize_, blockSize_));
1270 const int theRank = ortho_->normalize( *VV_0, zz_0 );
1271 TEUCHOS_TEST_FOR_EXCEPTION(
1273 "Belos::GmresPolySolMgr::solve: Failed to compute initial " 1274 "block of orthonormal vectors after restart.");
1278 theNewState.
V = VV_0;
1279 theNewState.
z = zz_0;
1281 block_gmres_iter->initializeGmres (theNewState);
1289 TEUCHOS_TEST_FOR_EXCEPTION(
1290 true, std::logic_error,
1291 "Belos::GmresPolySolMgr::solve: Invalid return from " 1292 "BlockGmresIter::iterate(). Please report this bug " 1293 "to the Belos developers.");
1298 if (blockSize_ != 1) {
1300 <<
"Error! Caught std::exception in BlockGmresIter::iterate() " 1301 <<
"at iteration " << block_gmres_iter->getNumIters()
1302 << std::endl << e.what() << std::endl;
1303 if (convTest_->getStatus() !=
Passed) {
1304 isConverged =
false;
1311 block_gmres_iter->updateLSQR (block_gmres_iter->getCurSubspaceDim ());
1315 sTest_->checkStatus (&*block_gmres_iter);
1316 if (convTest_->getStatus() !=
Passed) {
1317 isConverged =
false;
1322 catch (
const std::exception &e) {
1324 <<
"Error! Caught std::exception in BlockGmresIter::iterate() " 1325 <<
"at iteration " << block_gmres_iter->getNumIters() << std::endl
1326 << e.what() << std::endl;
1334 if (! Teuchos::is_null (expConvTest_->getSolution ()) ) {
1335 RCP<MV> newX = expConvTest_->getSolution ();
1336 RCP<MV> curX = problem_->getCurrLHSVec ();
1337 MVT::MvAddMv (STS::zero (), *newX, STS::one (), *newX, *curX);
1340 RCP<MV> update = block_gmres_iter->getCurrentUpdate ();
1341 problem_->updateSolution (update,
true);
1345 problem_->setCurrLS ();
1348 startPtr += numCurrRHS;
1349 numRHS2Solve -= numCurrRHS;
1350 if (numRHS2Solve > 0) {
1351 numCurrRHS = (numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1353 currIdx.resize (blockSize_);
1354 for (
int i=0; i<numCurrRHS; ++i) {
1355 currIdx[i] = startPtr+i;
1357 for (
int i=numCurrRHS; i<blockSize_; ++i) {
1362 problem_->setLSIndex( currIdx );
1365 currIdx.resize( numRHS2Solve );
1377 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1382 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1385 if (!isConverged || loaDetected_) {
1392 template<
class ScalarType,
class MV,
class OP>
1395 std::ostringstream out;
1397 out <<
"\"Belos::GmresPolySolMgr\": {" 1398 <<
"ScalarType: " << Teuchos::TypeNameTraits<ScalarType>::name ()
1399 <<
", Ortho Type: " << orthoType_
1400 <<
", Block Size: " << blockSize_
1401 <<
", Num Blocks: " << numBlocks_
1402 <<
", Max Restarts: " << maxRestarts_;
1409 #endif // BELOS_GMRES_POLY_SOLMGR_HPP Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Collection of types and exceptions used within the Belos solvers.
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.
Teuchos::RCP< const Teuchos::ParameterList > validPL_
Cached default (valid) parameters.
Teuchos::RCP< OutputManager< ScalarType > > printer_
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
ReturnType solve() override
This method performs possibly repeated calls to the underlying linear solver's iterate() routine unti...
virtual ~GmresPolySolMgr()
Destructor.
Defines the GMRES polynomial operator hybrid-GMRES iterative linear solver.
Teuchos::RCP< Teuchos::SerialDenseVector< int, ScalarType > > poly_r0_
static constexpr int outputStyle_default_
ScaleType
The type of scaling to use on the residual norm value.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
static constexpr double orthoKappa
DGKS orthogonalization constant.
MagnitudeType orthoKappa_
static constexpr const char * label_default_
This class implements the block GMRES iteration, where a block Krylov subspace is constructed...
Convergence test using the implicit residual norm(s), with an explicit residual norm(s) check for los...
Teuchos::RCP< const MV > V
The current Krylov basis.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > convTest_
static constexpr int maxDegree_default_
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< Belos::GmresPolyOp< ScalarType, MV, OP > > poly_Op_
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
An implementation of StatusTestResNorm using a family of residual norms.
Structure to contain pointers to GmresIteration state variables.
Belos::StatusTest class for specifying a maximum number of iterations.
Teuchos::RCP< Teuchos::ParameterList > params_
static constexpr int maxRestarts_default_
static constexpr int maxIters_default_
A factory class for generating StatusTestOutput objects.
static constexpr int blockSize_default_
Iterated Modified Gram-Schmidt (IMGS) implementation of the Belos::OrthoManager class.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
Teuchos::RCP< StatusTestResNorm< ScalarType, MV, OP > > impConvTest_
static constexpr const char * impResScale_default_
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
OperatorTraits< ScalarType, MV, OP > OPT
static constexpr std::ostream * outputStream_default_
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
GmresPolySolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
Pure virtual base class which describes the basic interface for a solver manager. ...
static constexpr double convTol
Default convergence tolerance.
Teuchos::RCP< Teuchos::Time > timerSolve_
Teuchos::ScalarTraits< MagnitudeType > MT
int curDim
The current dimension of the reduction.
Teuchos::ScalarTraits< ScalarType > STS
static constexpr bool strictConvTol_default_
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > R
The current upper-triangular matrix from the QR reduction of H.
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Get current linear problem being solved for in this object.
A linear system to solve, and its associated information.
static constexpr double polyTol
Relative residual tolerance for matrix polynomial construction.
static constexpr bool showMaxResNormOnly_default_
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Class which describes the linear problem to be solved by the iterative solver.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
MultiVecTraits< ScalarType, MV > MVT
Hybrid block GMRES iterative linear solver.
static constexpr int numBlocks_default_
static constexpr const char * orthoType_default_
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Teuchos::RCP< StatusTestResNorm< ScalarType, MV, OP > > expConvTest_
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 ...
Belos's class for applying the GMRES polynomial operator that is used by the hybrid-GMRES linear solv...
void reset(const ResetType type) override
Reset the solver.
Belos concrete class for performing the block GMRES iteration.
Teuchos::RCP< Teuchos::Time > timerPoly_
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
std::string description() const override
Method to return description of the hybrid block GMRES solver manager.
static constexpr const char * expResScale_default_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > poly_H_
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
GmresPolySolMgrPolynomialFailure is thrown when their is a problem generating the GMRES polynomial fo...
GmresPolySolMgrPolynomialFailure(const std::string &what_arg)
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
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 setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶ms) override
Set the parameters the solver manager should use to solve the linear problem.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest for specifying an implicit residual norm stopping criteria that checks for loss of ...
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
A class for extending the status testing capabilities of Belos via logical combinations.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
Get a parameter list containing the valid parameters for this object.
Class which defines basic traits for the operator type.
Teuchos::RCP< std::ostream > outputStream_
Parent class to all Belos exceptions.
Default parameters common to most Belos solvers.
GmresPolySolMgr()
Empty constructor for GmresPolySolMgr. This constructor takes no arguments and sets the default value...
Pure virtual base class which augments the basic interface for a Gmres linear solver iteration...
GmresIterationOrthoFailure is thrown when the GmresIteration object is unable to compute independent ...
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > poly_y_
Belos header file which uses auto-configuration information to include necessary C++ headers...
static constexpr int verbosity_default_
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
GmresPolySolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthon...
GmresPolySolMgrOrthoFailure(const std::string &what_arg)
GmresPolySolMgrLinearProblemFailure(const std::string &what_arg)
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
static constexpr int outputFreq_default_
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.