42 #ifndef BELOS_BLOCK_GMRES_SOLMGR_HPP 43 #define BELOS_BLOCK_GMRES_SOLMGR_HPP 66 #ifdef BELOS_TEUCHOS_TIME_MONITOR 67 #include "Teuchos_TimeMonitor.hpp" 121 template<
class ScalarType,
class MV,
class OP>
127 typedef Teuchos::ScalarTraits<ScalarType> SCT;
128 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
129 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
164 const Teuchos::RCP<Teuchos::ParameterList> &pl );
170 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
197 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
198 return Teuchos::tuple(timerSolve_);
234 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms )
override;
284 describe (Teuchos::FancyOStream& out,
285 const Teuchos::EVerbosityLevel verbLevel =
286 Teuchos::Describable::verbLevel_default)
const override;
296 bool checkStatusTest();
299 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
302 Teuchos::RCP<OutputManager<ScalarType> > printer_;
303 Teuchos::RCP<std::ostream> outputStream_;
306 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugStatusTest_;
307 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
308 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
309 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
310 Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
311 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
314 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
317 Teuchos::RCP<Teuchos::ParameterList> params_;
320 static constexpr
int maxRestarts_default_ = 20;
321 static constexpr
int maxIters_default_ = 1000;
322 static constexpr
bool adaptiveBlockSize_default_ =
true;
323 static constexpr
bool showMaxResNormOnly_default_ =
false;
324 static constexpr
bool flexibleGmres_default_ =
false;
325 static constexpr
bool expResTest_default_ =
false;
326 static constexpr
int blockSize_default_ = 1;
327 static constexpr
int numBlocks_default_ = 300;
330 static constexpr
int outputFreq_default_ = -1;
331 static constexpr
const char * impResScale_default_ =
"Norm of Preconditioned Initial Residual";
332 static constexpr
const char * expResScale_default_ =
"Norm of Initial Residual";
333 static constexpr
const char * label_default_ =
"Belos";
334 static constexpr
const char * orthoType_default_ =
"ICGS";
336 #if defined(_WIN32) && defined(__clang__) 337 static constexpr std::ostream * outputStream_default_ =
338 __builtin_constant_p(reinterpret_cast<const std::ostream*>(&std::cout));
340 static constexpr std::ostream * outputStream_default_ = &std::cout;
344 MagnitudeType convtol_, orthoKappa_, achievedTol_;
345 int maxRestarts_, maxIters_, numIters_;
346 int blockSize_, numBlocks_, verbosity_, outputStyle_, outputFreq_;
347 bool adaptiveBlockSize_, showMaxResNormOnly_, isFlexible_, expResTest_;
348 std::string orthoType_;
349 std::string impResScale_, expResScale_;
353 Teuchos::RCP<Teuchos::Time> timerSolve_;
356 bool isSet_, isSTSet_;
362 template<
class ScalarType,
class MV,
class OP>
364 outputStream_(Teuchos::rcp(outputStream_default_,false)),
367 achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
368 maxRestarts_(maxRestarts_default_),
369 maxIters_(maxIters_default_),
371 blockSize_(blockSize_default_),
372 numBlocks_(numBlocks_default_),
373 verbosity_(verbosity_default_),
374 outputStyle_(outputStyle_default_),
375 outputFreq_(outputFreq_default_),
376 adaptiveBlockSize_(adaptiveBlockSize_default_),
377 showMaxResNormOnly_(showMaxResNormOnly_default_),
378 isFlexible_(flexibleGmres_default_),
379 expResTest_(expResTest_default_),
380 orthoType_(orthoType_default_),
381 impResScale_(impResScale_default_),
382 expResScale_(expResScale_default_),
383 label_(label_default_),
391 template<
class ScalarType,
class MV,
class OP>
394 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
396 outputStream_(Teuchos::rcp(outputStream_default_,false)),
399 achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
400 maxRestarts_(maxRestarts_default_),
401 maxIters_(maxIters_default_),
403 blockSize_(blockSize_default_),
404 numBlocks_(numBlocks_default_),
405 verbosity_(verbosity_default_),
406 outputStyle_(outputStyle_default_),
407 outputFreq_(outputFreq_default_),
408 adaptiveBlockSize_(adaptiveBlockSize_default_),
409 showMaxResNormOnly_(showMaxResNormOnly_default_),
410 isFlexible_(flexibleGmres_default_),
411 expResTest_(expResTest_default_),
412 orthoType_(orthoType_default_),
413 impResScale_(impResScale_default_),
414 expResScale_(expResScale_default_),
415 label_(label_default_),
421 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
424 if ( !is_null(pl) ) {
431 template<
class ScalarType,
class MV,
class OP>
432 Teuchos::RCP<const Teuchos::ParameterList>
435 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
436 if (is_null(validPL)) {
437 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
442 "The relative residual tolerance that needs to be achieved by the\n" 443 "iterative solver in order for the linear system to be declared converged." );
444 pl->set(
"Maximum Restarts", static_cast<int>(maxRestarts_default_),
445 "The maximum number of restarts allowed for each\n" 446 "set of RHS solved.");
447 pl->set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
448 "The maximum number of block iterations allowed for each\n" 449 "set of RHS solved.");
450 pl->set(
"Num Blocks", static_cast<int>(numBlocks_default_),
451 "The maximum number of blocks allowed in the Krylov subspace\n" 452 "for each set of RHS solved.");
453 pl->set(
"Block Size", static_cast<int>(blockSize_default_),
454 "The number of vectors in each block. This number times the\n" 455 "number of blocks is the total Krylov subspace dimension.");
456 pl->set(
"Adaptive Block Size", static_cast<bool>(adaptiveBlockSize_default_),
457 "Whether the solver manager should adapt the block size\n" 458 "based on the number of RHS to solve.");
459 pl->set(
"Verbosity", static_cast<int>(verbosity_default_),
460 "What type(s) of solver information should be outputted\n" 461 "to the output stream.");
462 pl->set(
"Output Style", static_cast<int>(outputStyle_default_),
463 "What style is used for the solver information outputted\n" 464 "to the output stream.");
465 pl->set(
"Output Frequency", static_cast<int>(outputFreq_default_),
466 "How often convergence information should be outputted\n" 467 "to the output stream.");
468 pl->set(
"Output Stream", Teuchos::rcp(outputStream_default_,
false),
469 "A reference-counted pointer to the output stream where all\n" 470 "solver output is sent.");
471 pl->set(
"Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
472 "When convergence information is printed, only show the maximum\n" 473 "relative residual norm when the block size is greater than one.");
474 pl->set(
"Flexible Gmres", static_cast<bool>(flexibleGmres_default_),
475 "Whether the solver manager should use the flexible variant\n" 477 pl->set(
"Explicit Residual Test", static_cast<bool>(expResTest_default_),
478 "Whether the explicitly computed residual should be used in the convergence test.");
479 pl->set(
"Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
480 "The type of scaling used in the implicit residual convergence test.");
481 pl->set(
"Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
482 "The type of scaling used in the explicit residual convergence test.");
483 pl->set(
"Timer Label", static_cast<const char *>(label_default_),
484 "The string to use as a prefix for the timer labels.");
485 pl->set(
"Orthogonalization", static_cast<const char *>(orthoType_default_),
486 "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
488 "The constant used by DGKS orthogonalization to determine\n" 489 "whether another step of classical Gram-Schmidt is necessary.");
496 template<
class ScalarType,
class MV,
class OP>
501 if (params_ == Teuchos::null) {
502 params_ = Teuchos::rcp(
new Teuchos::ParameterList(*getValidParameters()) );
505 params->validateParameters(*getValidParameters());
509 if (params->isParameter(
"Maximum Restarts")) {
510 maxRestarts_ = params->get(
"Maximum Restarts",maxRestarts_default_);
513 params_->set(
"Maximum Restarts", maxRestarts_);
517 if (params->isParameter(
"Maximum Iterations")) {
518 maxIters_ = params->get(
"Maximum Iterations",maxIters_default_);
521 params_->set(
"Maximum Iterations", maxIters_);
522 if (maxIterTest_!=Teuchos::null)
523 maxIterTest_->setMaxIters( maxIters_ );
527 if (params->isParameter(
"Block Size")) {
528 blockSize_ = params->get(
"Block Size",blockSize_default_);
529 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
530 "Belos::BlockGmresSolMgr: \"Block Size\" must be strictly positive.");
533 params_->set(
"Block Size", blockSize_);
537 if (params->isParameter(
"Adaptive Block Size")) {
538 adaptiveBlockSize_ = params->get(
"Adaptive Block Size",adaptiveBlockSize_default_);
541 params_->set(
"Adaptive Block Size", adaptiveBlockSize_);
545 if (params->isParameter(
"Num Blocks")) {
546 numBlocks_ = params->get(
"Num Blocks",numBlocks_default_);
547 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
548 "Belos::BlockGmresSolMgr: \"Num Blocks\" must be strictly positive.");
551 params_->set(
"Num Blocks", numBlocks_);
555 if (params->isParameter(
"Timer Label")) {
556 std::string tempLabel = params->get(
"Timer Label", label_default_);
559 if (tempLabel != label_) {
561 params_->set(
"Timer Label", label_);
562 std::string solveLabel = label_ +
": BlockGmresSolMgr total solve time";
563 #ifdef BELOS_TEUCHOS_TIME_MONITOR 564 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
566 if (ortho_ != Teuchos::null) {
567 ortho_->setLabel( label_ );
573 if (params->isParameter(
"Flexible Gmres")) {
574 isFlexible_ = Teuchos::getParameter<bool>(*params,
"Flexible Gmres");
575 params_->set(
"Flexible Gmres", isFlexible_);
576 if (isFlexible_ && expResTest_) {
583 if (params->isParameter(
"Verbosity")) {
584 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
585 verbosity_ = params->get(
"Verbosity", verbosity_default_);
587 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
591 params_->set(
"Verbosity", verbosity_);
592 if (printer_ != Teuchos::null)
593 printer_->setVerbosity(verbosity_);
597 if (params->isParameter(
"Output Style")) {
598 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
599 outputStyle_ = params->get(
"Output Style", outputStyle_default_);
601 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
605 params_->set(
"Output Style", outputStyle_);
606 if (outputTest_ != Teuchos::null) {
612 if (params->isParameter(
"Output Stream")) {
613 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
616 params_->set(
"Output Stream", outputStream_);
617 if (printer_ != Teuchos::null)
618 printer_->setOStream( outputStream_ );
623 if (params->isParameter(
"Output Frequency")) {
624 outputFreq_ = params->get(
"Output Frequency", outputFreq_default_);
628 params_->set(
"Output Frequency", outputFreq_);
629 if (outputTest_ != Teuchos::null)
630 outputTest_->setOutputFrequency( outputFreq_ );
634 if (printer_ == Teuchos::null) {
639 bool changedOrthoType =
false;
640 if (params->isParameter(
"Orthogonalization")) {
641 std::string tempOrthoType = params->get(
"Orthogonalization",orthoType_default_);
642 if (tempOrthoType != orthoType_) {
643 orthoType_ = tempOrthoType;
644 changedOrthoType =
true;
647 params_->set(
"Orthogonalization", orthoType_);
650 if (params->isParameter(
"Orthogonalization Constant")) {
651 if (params->isType<MagnitudeType> (
"Orthogonalization Constant")) {
652 orthoKappa_ = params->get (
"Orthogonalization Constant",
656 orthoKappa_ = params->get (
"Orthogonalization Constant",
661 params_->set(
"Orthogonalization Constant",orthoKappa_);
662 if (orthoType_==
"DGKS") {
663 if (orthoKappa_ > 0 && ortho_ != Teuchos::null && !changedOrthoType) {
670 if (ortho_ == Teuchos::null || changedOrthoType) {
672 Teuchos::RCP<Teuchos::ParameterList> paramsOrtho;
673 if (orthoType_==
"DGKS" && orthoKappa_ > 0) {
674 paramsOrtho->set (
"depTol", orthoKappa_ );
677 ortho_ = factory.
makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
681 if (params->isParameter(
"Convergence Tolerance")) {
682 if (params->isType<MagnitudeType> (
"Convergence Tolerance")) {
683 convtol_ = params->get (
"Convergence Tolerance",
691 params_->set(
"Convergence Tolerance", convtol_);
692 if (impConvTest_ != Teuchos::null)
693 impConvTest_->setTolerance( convtol_ );
694 if (expConvTest_ != Teuchos::null)
695 expConvTest_->setTolerance( convtol_ );
699 if (params->isParameter(
"Implicit Residual Scaling")) {
700 std::string tempImpResScale = Teuchos::getParameter<std::string>( *params,
"Implicit Residual Scaling" );
703 if (impResScale_ != tempImpResScale) {
705 impResScale_ = tempImpResScale;
708 params_->set(
"Implicit Residual Scaling", impResScale_);
709 if (impConvTest_ != Teuchos::null) {
713 catch (std::exception& e) {
721 if (params->isParameter(
"Explicit Residual Scaling")) {
722 std::string tempExpResScale = Teuchos::getParameter<std::string>( *params,
"Explicit Residual Scaling" );
725 if (expResScale_ != tempExpResScale) {
727 expResScale_ = tempExpResScale;
730 params_->set(
"Explicit Residual Scaling", expResScale_);
731 if (expConvTest_ != Teuchos::null) {
735 catch (std::exception& e) {
743 if (params->isParameter(
"Explicit Residual Test")) {
744 expResTest_ = Teuchos::getParameter<bool>( *params,
"Explicit Residual Test" );
747 params_->set(
"Explicit Residual Test", expResTest_);
748 if (expConvTest_ == Teuchos::null) {
753 if (params->isParameter(
"Show Maximum Residual Norm Only")) {
754 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,
"Show Maximum Residual Norm Only");
757 params_->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
758 if (impConvTest_ != Teuchos::null)
759 impConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
760 if (expConvTest_ != Teuchos::null)
761 expConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
766 if (timerSolve_ == Teuchos::null) {
767 std::string solveLabel = label_ +
": BlockGmresSolMgr total solve time";
768 #ifdef BELOS_TEUCHOS_TIME_MONITOR 769 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
778 template<
class ScalarType,
class MV,
class OP>
791 if (isFlexible_ && Teuchos::is_null(problem_->getRightPrec())) {
793 params_->set(
"Flexible Gmres", isFlexible_);
797 "Belos::BlockGmresSolMgr::solve(): Linear problem has a left preconditioner, not a right preconditioner, which is required for flexible GMRES.");
802 if (!Teuchos::is_null(problem_->getLeftPrec()) && !isFlexible_) {
809 Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
810 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_ ) );
812 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
813 impConvTest_ = tmpImpConvTest;
816 Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
817 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_ ) );
818 tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit,
Belos::TwoNorm );
820 tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
821 expConvTest_ = tmpExpConvTest;
824 convTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
830 Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
831 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_ ) );
833 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
834 impConvTest_ = tmpImpConvTest;
839 Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
840 Teuchos::rcp(
new StatusTestImpResNorm_t( convtol_ ) );
842 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
843 impConvTest_ = tmpImpConvTest;
847 expConvTest_ = impConvTest_;
848 convTest_ = impConvTest_;
852 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
855 if (nonnull(debugStatusTest_) ) {
857 Teuchos::rcp_dynamic_cast<StatusTestCombo_t>(sTest_)->addStatusTest( debugStatusTest_ );
862 StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
866 std::string solverDesc =
" Block Gmres ";
868 solverDesc =
"Flexible" + solverDesc;
869 outputTest_->setSolverDesc( solverDesc );
877 template<
class ScalarType,
class MV,
class OP>
882 debugStatusTest_ = debugStatusTest;
887 template<
class ScalarType,
class MV,
class OP>
894 setParameters(Teuchos::parameterList(*getValidParameters()));
898 "Belos::BlockGmresSolMgr::solve(): Linear problem is not a valid object.");
901 "Belos::BlockGmresSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
903 if (!isSTSet_ || (!expResTest_ && !Teuchos::is_null(problem_->getLeftPrec())) ) {
905 "Belos::BlockGmresSolMgr::solve(): Linear problem and requested status tests are incompatible.");
910 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
911 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
913 std::vector<int> currIdx;
916 if ( adaptiveBlockSize_ ) {
917 blockSize_ = numCurrRHS;
918 currIdx.resize( numCurrRHS );
919 for (
int i=0; i<numCurrRHS; ++i)
920 { currIdx[i] = startPtr+i; }
924 currIdx.resize( blockSize_ );
925 for (
int i=0; i<numCurrRHS; ++i)
926 { currIdx[i] = startPtr+i; }
927 for (
int i=numCurrRHS; i<blockSize_; ++i)
932 problem_->setLSIndex( currIdx );
936 Teuchos::ParameterList plist;
937 plist.set(
"Block Size",blockSize_);
939 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
940 if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {
941 int tmpNumBlocks = 0;
943 tmpNumBlocks = dim / blockSize_;
945 tmpNumBlocks = ( dim - blockSize_) / blockSize_;
947 "Belos::BlockGmresSolMgr::solve(): Warning! Requested Krylov subspace dimension is larger than operator dimension!" 948 << std::endl <<
" The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
949 plist.set(
"Num Blocks",tmpNumBlocks);
952 plist.set(
"Num Blocks",numBlocks_);
955 outputTest_->reset();
956 loaDetected_ =
false;
959 bool isConverged =
true;
964 Teuchos::RCP<GmresIteration<ScalarType,MV,OP> > block_gmres_iter;
973 #ifdef BELOS_TEUCHOS_TIME_MONITOR 974 Teuchos::TimeMonitor slvtimer(*timerSolve_);
977 while ( numRHS2Solve > 0 ) {
980 if (blockSize_*numBlocks_ > dim) {
981 int tmpNumBlocks = 0;
983 tmpNumBlocks = dim / blockSize_;
985 tmpNumBlocks = ( dim - blockSize_) / blockSize_;
986 block_gmres_iter->setSize( blockSize_, tmpNumBlocks );
989 block_gmres_iter->setSize( blockSize_, numBlocks_ );
992 block_gmres_iter->resetNumIters();
995 outputTest_->resetNumCalls();
998 Teuchos::RCP<MV> V_0;
1001 if (currIdx[blockSize_-1] == -1) {
1002 V_0 = MVT::Clone( *(problem_->getInitResVec()), blockSize_ );
1003 problem_->computeCurrResVec( &*V_0 );
1006 V_0 = MVT::CloneCopy( *(problem_->getInitResVec()), currIdx );
1011 if (currIdx[blockSize_-1] == -1) {
1012 V_0 = MVT::Clone( *(problem_->getInitPrecResVec()), blockSize_ );
1013 problem_->computeCurrPrecResVec( &*V_0 );
1016 V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1021 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z_0 =
1022 Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( blockSize_, blockSize_ ) );
1025 int rank = ortho_->normalize( *V_0, z_0 );
1027 "Belos::BlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
1034 block_gmres_iter->initializeGmres(newstate);
1035 int numRestarts = 0;
1040 block_gmres_iter->iterate();
1047 if ( convTest_->getStatus() ==
Passed ) {
1048 if ( expConvTest_->getLOADetected() ) {
1050 loaDetected_ =
true;
1052 "Belos::BlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
1053 isConverged =
false;
1062 else if ( maxIterTest_->getStatus() ==
Passed ) {
1064 isConverged =
false;
1072 else if ( block_gmres_iter->getCurSubspaceDim() == block_gmres_iter->getMaxSubspaceDim() ) {
1074 if ( numRestarts >= maxRestarts_ ) {
1075 isConverged =
false;
1080 printer_->stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
1083 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1086 Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1087 MVT::MvAddMv( 1.0, *curX, 1.0, *update, *curX );
1090 problem_->updateSolution( update,
true );
1097 V_0 = MVT::Clone( *(oldState.
V), blockSize_ );
1099 problem_->computeCurrResVec( &*V_0 );
1101 problem_->computeCurrPrecResVec( &*V_0 );
1104 z_0 = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( blockSize_, blockSize_ ) );
1107 rank = ortho_->normalize( *V_0, z_0 );
1109 "Belos::BlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors after restart.");
1115 block_gmres_iter->initializeGmres(newstate);
1127 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
1128 "Belos::BlockGmresSolMgr::solve(): Invalid return from BlockGmresIter::iterate().");
1133 if (blockSize_ != 1) {
1134 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 1135 << block_gmres_iter->getNumIters() << std::endl
1136 << e.what() << std::endl;
1137 if (convTest_->getStatus() !=
Passed)
1138 isConverged =
false;
1143 block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
1146 sTest_->checkStatus( &*block_gmres_iter );
1147 if (convTest_->getStatus() !=
Passed)
1148 isConverged =
false;
1152 catch (
const std::exception &e) {
1153 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 1154 << block_gmres_iter->getNumIters() << std::endl
1155 << e.what() << std::endl;
1164 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1165 Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1167 if (update != Teuchos::null)
1168 MVT::MvAddMv( 1.0, *curX, 1.0, *update, *curX );
1172 if ( !Teuchos::is_null(expConvTest_->getSolution()) ) {
1173 Teuchos::RCP<MV> newX = expConvTest_->getSolution();
1174 Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1175 MVT::MvAddMv( 0.0, *newX, 1.0, *newX, *curX );
1178 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1179 problem_->updateSolution( update,
true );
1184 problem_->setCurrLS();
1187 startPtr += numCurrRHS;
1188 numRHS2Solve -= numCurrRHS;
1189 if ( numRHS2Solve > 0 ) {
1190 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1192 if ( adaptiveBlockSize_ ) {
1193 blockSize_ = numCurrRHS;
1194 currIdx.resize( numCurrRHS );
1195 for (
int i=0; i<numCurrRHS; ++i)
1196 { currIdx[i] = startPtr+i; }
1199 currIdx.resize( blockSize_ );
1200 for (
int i=0; i<numCurrRHS; ++i)
1201 { currIdx[i] = startPtr+i; }
1202 for (
int i=numCurrRHS; i<blockSize_; ++i)
1203 { currIdx[i] = -1; }
1206 problem_->setLSIndex( currIdx );
1209 currIdx.resize( numRHS2Solve );
1220 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1225 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1229 numIters_ = maxIterTest_->getNumIters();
1243 const std::vector<MagnitudeType>* pTestValues = NULL;
1245 pTestValues = expConvTest_->getTestValue();
1246 if (pTestValues == NULL || pTestValues->size() < 1) {
1247 pTestValues = impConvTest_->getTestValue();
1252 pTestValues = impConvTest_->getTestValue();
1254 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1255 "Belos::BlockGmresSolMgr::solve(): The implicit convergence test's " 1256 "getTestValue() method returned NULL. Please report this bug to the " 1257 "Belos developers.");
1258 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1259 "Belos::BlockGmresSolMgr::solve(): The implicit convergence test's " 1260 "getTestValue() method returned a vector of length zero. Please report " 1261 "this bug to the Belos developers.");
1266 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1269 if (!isConverged || loaDetected_) {
1276 template<
class ScalarType,
class MV,
class OP>
1279 std::ostringstream out;
1280 out <<
"\"Belos::BlockGmresSolMgr\": {";
1281 if (this->getObjectLabel () !=
"") {
1282 out <<
"Label: " << this->getObjectLabel () <<
", ";
1284 out <<
"Flexible: " << (isFlexible_ ?
"true" :
"false")
1285 <<
", Num Blocks: " << numBlocks_
1286 <<
", Maximum Iterations: " << maxIters_
1287 <<
", Maximum Restarts: " << maxRestarts_
1288 <<
", Convergence Tolerance: " << convtol_
1294 template<
class ScalarType,
class MV,
class OP>
1298 const Teuchos::EVerbosityLevel verbLevel)
const 1300 using Teuchos::TypeNameTraits;
1301 using Teuchos::VERB_DEFAULT;
1302 using Teuchos::VERB_NONE;
1303 using Teuchos::VERB_LOW;
1310 const Teuchos::EVerbosityLevel vl =
1311 (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
1313 if (vl != VERB_NONE) {
1314 Teuchos::OSTab tab0 (out);
1316 out <<
"\"Belos::BlockGmresSolMgr\":" << endl;
1317 Teuchos::OSTab tab1 (out);
1318 out <<
"Template parameters:" << endl;
1320 Teuchos::OSTab tab2 (out);
1321 out <<
"ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1322 <<
"MV: " << TypeNameTraits<MV>::name () << endl
1323 <<
"OP: " << TypeNameTraits<OP>::name () << endl;
1325 if (this->getObjectLabel () !=
"") {
1326 out <<
"Label: " << this->getObjectLabel () << endl;
1328 out <<
"Flexible: " << (isFlexible_ ?
"true" :
"false") << endl
1329 <<
"Num Blocks: " << numBlocks_ << endl
1330 <<
"Maximum Iterations: " << maxIters_ << endl
1331 <<
"Maximum Restarts: " << maxRestarts_ << endl
1332 <<
"Convergence Tolerance: " << convtol_ << endl;
static const double orthoKappa
DGKS orthogonalization constant.
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
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.
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.
static const double convTol
Default convergence tolerance.
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.
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.
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. ...
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.
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.
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.
Teuchos::RCP< Belos::MatOrthoManager< Scalar, MV, OP > > makeMatOrthoManager(const std::string &ortho, const Teuchos::RCP< const OP > &M, const Teuchos::RCP< OutputManager< Scalar > > &, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Return an instance of the specified MatOrthoManager subclass.