42 #ifndef BELOS_BLOCK_GMRES_SOLMGR_HPP 43 #define BELOS_BLOCK_GMRES_SOLMGR_HPP 68 #include "Teuchos_BLAS.hpp" 69 #ifdef BELOS_TEUCHOS_TIME_MONITOR 70 #include "Teuchos_TimeMonitor.hpp" 124 template<
class ScalarType,
class MV,
class OP>
130 typedef Teuchos::ScalarTraits<ScalarType> SCT;
131 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
132 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
167 const Teuchos::RCP<Teuchos::ParameterList> &pl );
173 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
200 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
201 return Teuchos::tuple(timerSolve_);
237 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms )
override;
287 describe (Teuchos::FancyOStream& out,
288 const Teuchos::EVerbosityLevel verbLevel =
289 Teuchos::Describable::verbLevel_default)
const override;
299 bool checkStatusTest();
302 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
305 Teuchos::RCP<OutputManager<ScalarType> > printer_;
306 Teuchos::RCP<std::ostream> outputStream_;
309 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugStatusTest_;
310 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
311 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
312 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
313 Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
314 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
317 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
320 Teuchos::RCP<Teuchos::ParameterList> params_;
323 static constexpr
int maxRestarts_default_ = 20;
324 static constexpr
int maxIters_default_ = 1000;
325 static constexpr
bool adaptiveBlockSize_default_ =
true;
326 static constexpr
bool showMaxResNormOnly_default_ =
false;
327 static constexpr
bool flexibleGmres_default_ =
false;
328 static constexpr
bool expResTest_default_ =
false;
329 static constexpr
int blockSize_default_ = 1;
330 static constexpr
int numBlocks_default_ = 300;
333 static constexpr
int outputFreq_default_ = -1;
334 static constexpr
const char * impResScale_default_ =
"Norm of Preconditioned Initial Residual";
335 static constexpr
const char * expResScale_default_ =
"Norm of Initial Residual";
336 static constexpr
const char * label_default_ =
"Belos";
337 static constexpr
const char * orthoType_default_ =
"DGKS";
338 static constexpr std::ostream * outputStream_default_ = &std::cout;
341 MagnitudeType convtol_, orthoKappa_, achievedTol_;
342 int maxRestarts_, maxIters_, numIters_;
343 int blockSize_, numBlocks_, verbosity_, outputStyle_, outputFreq_;
344 bool adaptiveBlockSize_, showMaxResNormOnly_, isFlexible_, expResTest_;
345 std::string orthoType_;
346 std::string impResScale_, expResScale_;
350 Teuchos::RCP<Teuchos::Time> timerSolve_;
353 bool isSet_, isSTSet_;
359 template<
class ScalarType,
class MV,
class OP>
361 outputStream_(Teuchos::rcp(outputStream_default_,false)),
364 achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
365 maxRestarts_(maxRestarts_default_),
366 maxIters_(maxIters_default_),
368 blockSize_(blockSize_default_),
369 numBlocks_(numBlocks_default_),
370 verbosity_(verbosity_default_),
371 outputStyle_(outputStyle_default_),
372 outputFreq_(outputFreq_default_),
373 adaptiveBlockSize_(adaptiveBlockSize_default_),
374 showMaxResNormOnly_(showMaxResNormOnly_default_),
375 isFlexible_(flexibleGmres_default_),
376 expResTest_(expResTest_default_),
377 orthoType_(orthoType_default_),
378 impResScale_(impResScale_default_),
379 expResScale_(expResScale_default_),
380 label_(label_default_),
388 template<
class ScalarType,
class MV,
class OP>
391 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
393 outputStream_(Teuchos::rcp(outputStream_default_,false)),
396 achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
397 maxRestarts_(maxRestarts_default_),
398 maxIters_(maxIters_default_),
400 blockSize_(blockSize_default_),
401 numBlocks_(numBlocks_default_),
402 verbosity_(verbosity_default_),
403 outputStyle_(outputStyle_default_),
404 outputFreq_(outputFreq_default_),
405 adaptiveBlockSize_(adaptiveBlockSize_default_),
406 showMaxResNormOnly_(showMaxResNormOnly_default_),
407 isFlexible_(flexibleGmres_default_),
408 expResTest_(expResTest_default_),
409 orthoType_(orthoType_default_),
410 impResScale_(impResScale_default_),
411 expResScale_(expResScale_default_),
412 label_(label_default_),
418 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
421 if ( !is_null(pl) ) {
428 template<
class ScalarType,
class MV,
class OP>
429 Teuchos::RCP<const Teuchos::ParameterList>
432 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
433 if (is_null(validPL)) {
434 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
439 "The relative residual tolerance that needs to be achieved by the\n" 440 "iterative solver in order for the linear system to be declared converged." );
441 pl->set(
"Maximum Restarts", static_cast<int>(maxRestarts_default_),
442 "The maximum number of restarts allowed for each\n" 443 "set of RHS solved.");
444 pl->set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
445 "The maximum number of block iterations allowed for each\n" 446 "set of RHS solved.");
447 pl->set(
"Num Blocks", static_cast<int>(numBlocks_default_),
448 "The maximum number of blocks allowed in the Krylov subspace\n" 449 "for each set of RHS solved.");
450 pl->set(
"Block Size", static_cast<int>(blockSize_default_),
451 "The number of vectors in each block. This number times the\n" 452 "number of blocks is the total Krylov subspace dimension.");
453 pl->set(
"Adaptive Block Size", static_cast<bool>(adaptiveBlockSize_default_),
454 "Whether the solver manager should adapt the block size\n" 455 "based on the number of RHS to solve.");
456 pl->set(
"Verbosity", static_cast<int>(verbosity_default_),
457 "What type(s) of solver information should be outputted\n" 458 "to the output stream.");
459 pl->set(
"Output Style", static_cast<int>(outputStyle_default_),
460 "What style is used for the solver information outputted\n" 461 "to the output stream.");
462 pl->set(
"Output Frequency", static_cast<int>(outputFreq_default_),
463 "How often convergence information should be outputted\n" 464 "to the output stream.");
465 pl->set(
"Output Stream", Teuchos::rcp(outputStream_default_,
false),
466 "A reference-counted pointer to the output stream where all\n" 467 "solver output is sent.");
468 pl->set(
"Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
469 "When convergence information is printed, only show the maximum\n" 470 "relative residual norm when the block size is greater than one.");
471 pl->set(
"Flexible Gmres", static_cast<bool>(flexibleGmres_default_),
472 "Whether the solver manager should use the flexible variant\n" 474 pl->set(
"Explicit Residual Test", static_cast<bool>(expResTest_default_),
475 "Whether the explicitly computed residual should be used in the convergence test.");
476 pl->set(
"Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
477 "The type of scaling used in the implicit residual convergence test.");
478 pl->set(
"Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
479 "The type of scaling used in the explicit residual convergence test.");
480 pl->set(
"Timer Label", static_cast<const char *>(label_default_),
481 "The string to use as a prefix for the timer labels.");
482 pl->set(
"Orthogonalization", static_cast<const char *>(orthoType_default_),
483 "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
485 "The constant used by DGKS orthogonalization to determine\n" 486 "whether another step of classical Gram-Schmidt is necessary.");
493 template<
class ScalarType,
class MV,
class OP>
498 if (params_ == Teuchos::null) {
499 params_ = Teuchos::rcp(
new Teuchos::ParameterList(*getValidParameters()) );
502 params->validateParameters(*getValidParameters());
506 if (params->isParameter(
"Maximum Restarts")) {
507 maxRestarts_ = params->get(
"Maximum Restarts",maxRestarts_default_);
510 params_->set(
"Maximum Restarts", maxRestarts_);
514 if (params->isParameter(
"Maximum Iterations")) {
515 maxIters_ = params->get(
"Maximum Iterations",maxIters_default_);
518 params_->set(
"Maximum Iterations", maxIters_);
519 if (maxIterTest_!=Teuchos::null)
520 maxIterTest_->setMaxIters( maxIters_ );
524 if (params->isParameter(
"Block Size")) {
525 blockSize_ = params->get(
"Block Size",blockSize_default_);
526 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
527 "Belos::BlockGmresSolMgr: \"Block Size\" must be strictly positive.");
530 params_->set(
"Block Size", blockSize_);
534 if (params->isParameter(
"Adaptive Block Size")) {
535 adaptiveBlockSize_ = params->get(
"Adaptive Block Size",adaptiveBlockSize_default_);
538 params_->set(
"Adaptive Block Size", adaptiveBlockSize_);
542 if (params->isParameter(
"Num Blocks")) {
543 numBlocks_ = params->get(
"Num Blocks",numBlocks_default_);
544 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
545 "Belos::BlockGmresSolMgr: \"Num Blocks\" must be strictly positive.");
548 params_->set(
"Num Blocks", numBlocks_);
552 if (params->isParameter(
"Timer Label")) {
553 std::string tempLabel = params->get(
"Timer Label", label_default_);
556 if (tempLabel != label_) {
558 params_->set(
"Timer Label", label_);
559 std::string solveLabel = label_ +
": BlockGmresSolMgr total solve time";
560 #ifdef BELOS_TEUCHOS_TIME_MONITOR 561 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
563 if (ortho_ != Teuchos::null) {
564 ortho_->setLabel( label_ );
570 if (params->isParameter(
"Flexible Gmres")) {
571 isFlexible_ = Teuchos::getParameter<bool>(*params,
"Flexible Gmres");
572 params_->set(
"Flexible Gmres", isFlexible_);
573 if (isFlexible_ && expResTest_) {
581 if (params->isParameter(
"Orthogonalization")) {
582 std::string tempOrthoType = params->get(
"Orthogonalization",orthoType_default_);
583 TEUCHOS_TEST_FOR_EXCEPTION( tempOrthoType !=
"DGKS" && tempOrthoType !=
"ICGS" && tempOrthoType !=
"IMGS",
584 std::invalid_argument,
585 "Belos::BlockGmresSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
586 if (tempOrthoType != orthoType_) {
587 orthoType_ = tempOrthoType;
588 params_->set(
"Orthogonalization", orthoType_);
590 if (orthoType_==
"DGKS") {
591 if (orthoKappa_ <= 0) {
599 else if (orthoType_==
"ICGS") {
602 else if (orthoType_==
"IMGS") {
609 if (params->isParameter(
"Orthogonalization Constant")) {
610 if (params->isType<MagnitudeType> (
"Orthogonalization Constant")) {
611 orthoKappa_ = params->get (
"Orthogonalization Constant",
615 orthoKappa_ = params->get (
"Orthogonalization Constant",
620 params_->set(
"Orthogonalization Constant",orthoKappa_);
621 if (orthoType_==
"DGKS") {
622 if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
629 if (params->isParameter(
"Verbosity")) {
630 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
631 verbosity_ = params->get(
"Verbosity", verbosity_default_);
633 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
637 params_->set(
"Verbosity", verbosity_);
638 if (printer_ != Teuchos::null)
639 printer_->setVerbosity(verbosity_);
643 if (params->isParameter(
"Output Style")) {
644 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
645 outputStyle_ = params->get(
"Output Style", outputStyle_default_);
647 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
651 params_->set(
"Output Style", outputStyle_);
652 if (outputTest_ != Teuchos::null) {
658 if (params->isParameter(
"Output Stream")) {
659 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
662 params_->set(
"Output Stream", outputStream_);
663 if (printer_ != Teuchos::null)
664 printer_->setOStream( outputStream_ );
669 if (params->isParameter(
"Output Frequency")) {
670 outputFreq_ = params->get(
"Output Frequency", outputFreq_default_);
674 params_->set(
"Output Frequency", outputFreq_);
675 if (outputTest_ != Teuchos::null)
676 outputTest_->setOutputFrequency( outputFreq_ );
680 if (printer_ == Teuchos::null) {
685 if (params->isParameter(
"Convergence Tolerance")) {
686 if (params->isType<MagnitudeType> (
"Convergence Tolerance")) {
687 convtol_ = params->get (
"Convergence Tolerance",
695 params_->set(
"Convergence Tolerance", convtol_);
696 if (impConvTest_ != Teuchos::null)
697 impConvTest_->setTolerance( convtol_ );
698 if (expConvTest_ != Teuchos::null)
699 expConvTest_->setTolerance( convtol_ );
703 if (params->isParameter(
"Implicit Residual Scaling")) {
704 std::string tempImpResScale = Teuchos::getParameter<std::string>( *params,
"Implicit Residual Scaling" );
707 if (impResScale_ != tempImpResScale) {
709 impResScale_ = tempImpResScale;
712 params_->set(
"Implicit Residual Scaling", impResScale_);
713 if (impConvTest_ != Teuchos::null) {
717 catch (std::exception& e) {
725 if (params->isParameter(
"Explicit Residual Scaling")) {
726 std::string tempExpResScale = Teuchos::getParameter<std::string>( *params,
"Explicit Residual Scaling" );
729 if (expResScale_ != tempExpResScale) {
731 expResScale_ = tempExpResScale;
734 params_->set(
"Explicit Residual Scaling", expResScale_);
735 if (expConvTest_ != Teuchos::null) {
739 catch (std::exception& e) {
747 if (params->isParameter(
"Explicit Residual Test")) {
748 expResTest_ = Teuchos::getParameter<bool>( *params,
"Explicit Residual Test" );
751 params_->set(
"Explicit Residual Test", expResTest_);
752 if (expConvTest_ == Teuchos::null) {
758 if (params->isParameter(
"Show Maximum Residual Norm Only")) {
759 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,
"Show Maximum Residual Norm Only");
762 params_->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
763 if (impConvTest_ != Teuchos::null)
764 impConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
765 if (expConvTest_ != Teuchos::null)
766 expConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
770 if (ortho_ == Teuchos::null) {
771 params_->set(
"Orthogonalization", orthoType_);
772 if (orthoType_==
"DGKS") {
773 if (orthoKappa_ <= 0) {
781 else if (orthoType_==
"ICGS") {
784 else if (orthoType_==
"IMGS") {
788 TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!=
"ICGS"&&orthoType_!=
"DGKS"&&orthoType_!=
"IMGS",std::logic_error,
789 "Belos::BlockGmresSolMgr(): Invalid orthogonalization type.");
794 if (timerSolve_ == Teuchos::null) {
795 std::string solveLabel = label_ +
": BlockGmresSolMgr total solve time";
796 #ifdef BELOS_TEUCHOS_TIME_MONITOR 797 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
806 template<
class ScalarType,
class MV,
class OP>
818 if (!Teuchos::is_null(problem_->getLeftPrec()) && !isFlexible_) {
825 Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
826 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_ ) );
828 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
829 impConvTest_ = tmpImpConvTest;
832 Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
833 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_ ) );
834 tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit,
Belos::TwoNorm );
836 tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
837 expConvTest_ = tmpExpConvTest;
840 convTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
846 Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
847 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_ ) );
849 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
850 impConvTest_ = tmpImpConvTest;
855 Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
856 Teuchos::rcp(
new StatusTestImpResNorm_t( convtol_ ) );
858 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
859 impConvTest_ = tmpImpConvTest;
863 expConvTest_ = impConvTest_;
864 convTest_ = impConvTest_;
868 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
871 if (nonnull(debugStatusTest_) ) {
873 Teuchos::rcp_dynamic_cast<StatusTestCombo_t>(sTest_)->addStatusTest( debugStatusTest_ );
878 StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
882 std::string solverDesc =
" Block Gmres ";
884 solverDesc =
"Flexible" + solverDesc;
885 outputTest_->setSolverDesc( solverDesc );
893 template<
class ScalarType,
class MV,
class OP>
898 debugStatusTest_ = debugStatusTest;
903 template<
class ScalarType,
class MV,
class OP>
910 setParameters(Teuchos::parameterList(*getValidParameters()));
913 Teuchos::BLAS<int,ScalarType> blas;
916 "Belos::BlockGmresSolMgr::solve(): Linear problem is not a valid object.");
919 "Belos::BlockGmresSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
923 "Belos::BlockGmresSolMgr::solve(): Linear problem does not have a preconditioner required for flexible GMRES, call setRightPrec().");
926 if (!isSTSet_ || (!expResTest_ && !Teuchos::is_null(problem_->getLeftPrec())) ) {
928 "Belos::BlockGmresSolMgr::solve(): Linear problem and requested status tests are incompatible.");
933 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
934 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
936 std::vector<int> currIdx;
939 if ( adaptiveBlockSize_ ) {
940 blockSize_ = numCurrRHS;
941 currIdx.resize( numCurrRHS );
942 for (
int i=0; i<numCurrRHS; ++i)
943 { currIdx[i] = startPtr+i; }
947 currIdx.resize( blockSize_ );
948 for (
int i=0; i<numCurrRHS; ++i)
949 { currIdx[i] = startPtr+i; }
950 for (
int i=numCurrRHS; i<blockSize_; ++i)
955 problem_->setLSIndex( currIdx );
959 Teuchos::ParameterList plist;
960 plist.set(
"Block Size",blockSize_);
962 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
963 if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {
964 int tmpNumBlocks = 0;
966 tmpNumBlocks = dim / blockSize_;
968 tmpNumBlocks = ( dim - blockSize_) / blockSize_;
970 "Belos::BlockGmresSolMgr::solve(): Warning! Requested Krylov subspace dimension is larger than operator dimension!" 971 << std::endl <<
" The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
972 plist.set(
"Num Blocks",tmpNumBlocks);
975 plist.set(
"Num Blocks",numBlocks_);
978 outputTest_->reset();
979 loaDetected_ =
false;
982 bool isConverged =
true;
987 Teuchos::RCP<GmresIteration<ScalarType,MV,OP> > block_gmres_iter;
996 #ifdef BELOS_TEUCHOS_TIME_MONITOR 997 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1000 while ( numRHS2Solve > 0 ) {
1003 if (blockSize_*numBlocks_ > dim) {
1004 int tmpNumBlocks = 0;
1005 if (blockSize_ == 1)
1006 tmpNumBlocks = dim / blockSize_;
1008 tmpNumBlocks = ( dim - blockSize_) / blockSize_;
1009 block_gmres_iter->setSize( blockSize_, tmpNumBlocks );
1012 block_gmres_iter->setSize( blockSize_, numBlocks_ );
1015 block_gmres_iter->resetNumIters();
1018 outputTest_->resetNumCalls();
1021 Teuchos::RCP<MV> V_0;
1024 if (currIdx[blockSize_-1] == -1) {
1025 V_0 = MVT::Clone( *(problem_->getInitResVec()), blockSize_ );
1026 problem_->computeCurrResVec( &*V_0 );
1029 V_0 = MVT::CloneCopy( *(problem_->getInitResVec()), currIdx );
1034 if (currIdx[blockSize_-1] == -1) {
1035 V_0 = MVT::Clone( *(problem_->getInitPrecResVec()), blockSize_ );
1036 problem_->computeCurrPrecResVec( &*V_0 );
1039 V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1044 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z_0 =
1045 Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( blockSize_, blockSize_ ) );
1048 int rank = ortho_->normalize( *V_0, z_0 );
1050 "Belos::BlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
1057 block_gmres_iter->initializeGmres(newstate);
1058 int numRestarts = 0;
1063 block_gmres_iter->iterate();
1070 if ( convTest_->getStatus() ==
Passed ) {
1071 if ( expConvTest_->getLOADetected() ) {
1073 loaDetected_ =
true;
1075 "Belos::BlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
1076 isConverged =
false;
1085 else if ( maxIterTest_->getStatus() ==
Passed ) {
1087 isConverged =
false;
1095 else if ( block_gmres_iter->getCurSubspaceDim() == block_gmres_iter->getMaxSubspaceDim() ) {
1097 if ( numRestarts >= maxRestarts_ ) {
1098 isConverged =
false;
1103 printer_->stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
1106 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1109 Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1110 MVT::MvAddMv( 1.0, *curX, 1.0, *update, *curX );
1113 problem_->updateSolution( update,
true );
1120 V_0 = MVT::Clone( *(oldState.
V), blockSize_ );
1122 problem_->computeCurrResVec( &*V_0 );
1124 problem_->computeCurrPrecResVec( &*V_0 );
1127 z_0 = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( blockSize_, blockSize_ ) );
1130 rank = ortho_->normalize( *V_0, z_0 );
1132 "Belos::BlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors after restart.");
1138 block_gmres_iter->initializeGmres(newstate);
1150 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
1151 "Belos::BlockGmresSolMgr::solve(): Invalid return from BlockGmresIter::iterate().");
1156 if (blockSize_ != 1) {
1157 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 1158 << block_gmres_iter->getNumIters() << std::endl
1159 << e.what() << std::endl;
1160 if (convTest_->getStatus() !=
Passed)
1161 isConverged =
false;
1166 block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
1169 sTest_->checkStatus( &*block_gmres_iter );
1170 if (convTest_->getStatus() !=
Passed)
1171 isConverged =
false;
1175 catch (
const std::exception &e) {
1176 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 1177 << block_gmres_iter->getNumIters() << std::endl
1178 << e.what() << std::endl;
1187 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1188 Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1190 if (update != Teuchos::null)
1191 MVT::MvAddMv( 1.0, *curX, 1.0, *update, *curX );
1195 if ( !Teuchos::is_null(expConvTest_->getSolution()) ) {
1196 Teuchos::RCP<MV> newX = expConvTest_->getSolution();
1197 Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1198 MVT::MvAddMv( 0.0, *newX, 1.0, *newX, *curX );
1201 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1202 problem_->updateSolution( update,
true );
1207 problem_->setCurrLS();
1210 startPtr += numCurrRHS;
1211 numRHS2Solve -= numCurrRHS;
1212 if ( numRHS2Solve > 0 ) {
1213 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1215 if ( adaptiveBlockSize_ ) {
1216 blockSize_ = numCurrRHS;
1217 currIdx.resize( numCurrRHS );
1218 for (
int i=0; i<numCurrRHS; ++i)
1219 { currIdx[i] = startPtr+i; }
1222 currIdx.resize( blockSize_ );
1223 for (
int i=0; i<numCurrRHS; ++i)
1224 { currIdx[i] = startPtr+i; }
1225 for (
int i=numCurrRHS; i<blockSize_; ++i)
1226 { currIdx[i] = -1; }
1229 problem_->setLSIndex( currIdx );
1232 currIdx.resize( numRHS2Solve );
1243 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1248 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1252 numIters_ = maxIterTest_->getNumIters();
1266 const std::vector<MagnitudeType>* pTestValues = NULL;
1268 pTestValues = expConvTest_->getTestValue();
1269 if (pTestValues == NULL || pTestValues->size() < 1) {
1270 pTestValues = impConvTest_->getTestValue();
1275 pTestValues = impConvTest_->getTestValue();
1277 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1278 "Belos::BlockGmresSolMgr::solve(): The implicit convergence test's " 1279 "getTestValue() method returned NULL. Please report this bug to the " 1280 "Belos developers.");
1281 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1282 "Belos::BlockGmresSolMgr::solve(): The implicit convergence test's " 1283 "getTestValue() method returned a vector of length zero. Please report " 1284 "this bug to the Belos developers.");
1289 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1292 if (!isConverged || loaDetected_) {
1299 template<
class ScalarType,
class MV,
class OP>
1302 std::ostringstream out;
1303 out <<
"\"Belos::BlockGmresSolMgr\": {";
1304 if (this->getObjectLabel () !=
"") {
1305 out <<
"Label: " << this->getObjectLabel () <<
", ";
1307 out <<
"Flexible: " << (isFlexible_ ?
"true" :
"false")
1308 <<
", Num Blocks: " << numBlocks_
1309 <<
", Maximum Iterations: " << maxIters_
1310 <<
", Maximum Restarts: " << maxRestarts_
1311 <<
", Convergence Tolerance: " << convtol_
1317 template<
class ScalarType,
class MV,
class OP>
1321 const Teuchos::EVerbosityLevel verbLevel)
const 1323 using Teuchos::TypeNameTraits;
1324 using Teuchos::VERB_DEFAULT;
1325 using Teuchos::VERB_NONE;
1326 using Teuchos::VERB_LOW;
1333 const Teuchos::EVerbosityLevel vl =
1334 (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
1336 if (vl != VERB_NONE) {
1337 Teuchos::OSTab tab0 (out);
1339 out <<
"\"Belos::BlockGmresSolMgr\":" << endl;
1340 Teuchos::OSTab tab1 (out);
1341 out <<
"Template parameters:" << endl;
1343 Teuchos::OSTab tab2 (out);
1344 out <<
"ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1345 <<
"MV: " << TypeNameTraits<MV>::name () << endl
1346 <<
"OP: " << TypeNameTraits<OP>::name () << endl;
1348 if (this->getObjectLabel () !=
"") {
1349 out <<
"Label: " << this->getObjectLabel () << endl;
1351 out <<
"Flexible: " << (isFlexible_ ?
"true" :
"false") << endl
1352 <<
"Num Blocks: " << numBlocks_ << endl
1353 <<
"Maximum Iterations: " << maxIters_ << endl
1354 <<
"Maximum Restarts: " << maxRestarts_ << endl
1355 <<
"Convergence Tolerance: " << convtol_ << endl;
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::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...
ScaleType
The type of scaling to use on the residual norm value.
static constexpr double orthoKappa
DGKS orthogonalization constant.
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.
Virtual base class for StatusTest that printing status tests.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
Get a parameter list containing the valid parameters for this object.
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.
BlockGmresSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
Structure to contain pointers to GmresIteration state variables.
Belos::StatusTest class for specifying a maximum number of iterations.
Interface to Block GMRES and Flexible GMRES.
A pure virtual class for defining the status tests for the Belos iterative solvers.
virtual ~BlockGmresSolMgr()
Destructor.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
A factory class for generating StatusTestOutput objects.
Iterated Modified Gram-Schmidt (IMGS) implementation of the Belos::OrthoManager class.
ReturnType solve() override
This method performs possibly repeated calls to the underlying linear solver's iterate() routine unti...
Traits class which defines basic operations on multivectors.
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶ms) override
Set the parameters the solver manager should use to solve the linear problem.
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.
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.
int curDim
The current dimension of the reduction.
BlockGmresSolMgrLinearProblemFailure(const std::string &what_arg)
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
This class implements the block flexible GMRES iteration, where a block Krylov subspace is constructe...
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
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 concrete class for performing the block GMRES iteration.
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
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.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest for specifying an implicit residual norm stopping criteria that checks for loss of ...
BlockGmresSolMgr()
Empty constructor for BlockGmresSolMgr. This constructor takes no arguments and sets the default valu...
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
BlockGmresSolMgrOrthoFailure(const std::string &what_arg)
A class for extending the status testing capabilities of Belos via logical combinations.
std::string description() const override
Return a one-line description of this object.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
BlockGmresSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate ortho...
Default parameters common to most Belos solvers.
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 ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Get current linear problem being solved for in this object.
Belos concrete class for performing the block, flexible GMRES iteration.