42 #ifndef BELOS_BLOCK_GMRES_SOLMGR_HPP 43 #define BELOS_BLOCK_GMRES_SOLMGR_HPP 66 #include "Teuchos_BLAS.hpp" 67 #ifdef BELOS_TEUCHOS_TIME_MONITOR 68 #include "Teuchos_TimeMonitor.hpp" 122 template<
class ScalarType,
class MV,
class OP>
128 typedef Teuchos::ScalarTraits<ScalarType> SCT;
129 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
130 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
165 const Teuchos::RCP<Teuchos::ParameterList> &pl );
171 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
198 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
199 return Teuchos::tuple(timerSolve_);
235 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms )
override;
285 describe (Teuchos::FancyOStream& out,
286 const Teuchos::EVerbosityLevel verbLevel =
287 Teuchos::Describable::verbLevel_default)
const override;
297 bool checkStatusTest();
300 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
303 Teuchos::RCP<OutputManager<ScalarType> > printer_;
304 Teuchos::RCP<std::ostream> outputStream_;
307 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugStatusTest_;
308 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
309 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
310 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
311 Teuchos::RCP<StatusTestResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
312 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
315 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
318 Teuchos::RCP<Teuchos::ParameterList> params_;
321 static constexpr
int maxRestarts_default_ = 20;
322 static constexpr
int maxIters_default_ = 1000;
323 static constexpr
bool adaptiveBlockSize_default_ =
true;
324 static constexpr
bool showMaxResNormOnly_default_ =
false;
325 static constexpr
bool flexibleGmres_default_ =
false;
326 static constexpr
bool expResTest_default_ =
false;
327 static constexpr
int blockSize_default_ = 1;
328 static constexpr
int numBlocks_default_ = 300;
331 static constexpr
int outputFreq_default_ = -1;
332 static constexpr
const char * impResScale_default_ =
"Norm of Preconditioned Initial Residual";
333 static constexpr
const char * expResScale_default_ =
"Norm of Initial Residual";
334 static constexpr
const char * label_default_ =
"Belos";
335 static constexpr
const char * orthoType_default_ =
"ICGS";
336 static constexpr std::ostream * outputStream_default_ = &std::cout;
339 MagnitudeType convtol_, orthoKappa_, achievedTol_;
340 int maxRestarts_, maxIters_, numIters_;
341 int blockSize_, numBlocks_, verbosity_, outputStyle_, outputFreq_;
342 bool adaptiveBlockSize_, showMaxResNormOnly_, isFlexible_, expResTest_;
343 std::string orthoType_;
344 std::string impResScale_, expResScale_;
348 Teuchos::RCP<Teuchos::Time> timerSolve_;
351 bool isSet_, isSTSet_;
357 template<
class ScalarType,
class MV,
class OP>
359 outputStream_(Teuchos::rcp(outputStream_default_,false)),
362 achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
363 maxRestarts_(maxRestarts_default_),
364 maxIters_(maxIters_default_),
366 blockSize_(blockSize_default_),
367 numBlocks_(numBlocks_default_),
368 verbosity_(verbosity_default_),
369 outputStyle_(outputStyle_default_),
370 outputFreq_(outputFreq_default_),
371 adaptiveBlockSize_(adaptiveBlockSize_default_),
372 showMaxResNormOnly_(showMaxResNormOnly_default_),
373 isFlexible_(flexibleGmres_default_),
374 expResTest_(expResTest_default_),
375 orthoType_(orthoType_default_),
376 impResScale_(impResScale_default_),
377 expResScale_(expResScale_default_),
378 label_(label_default_),
386 template<
class ScalarType,
class MV,
class OP>
389 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
391 outputStream_(Teuchos::rcp(outputStream_default_,false)),
394 achievedTol_(Teuchos::ScalarTraits<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType>::zero()),
395 maxRestarts_(maxRestarts_default_),
396 maxIters_(maxIters_default_),
398 blockSize_(blockSize_default_),
399 numBlocks_(numBlocks_default_),
400 verbosity_(verbosity_default_),
401 outputStyle_(outputStyle_default_),
402 outputFreq_(outputFreq_default_),
403 adaptiveBlockSize_(adaptiveBlockSize_default_),
404 showMaxResNormOnly_(showMaxResNormOnly_default_),
405 isFlexible_(flexibleGmres_default_),
406 expResTest_(expResTest_default_),
407 orthoType_(orthoType_default_),
408 impResScale_(impResScale_default_),
409 expResScale_(expResScale_default_),
410 label_(label_default_),
416 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument,
"Problem not given to solver manager.");
419 if ( !is_null(pl) ) {
426 template<
class ScalarType,
class MV,
class OP>
427 Teuchos::RCP<const Teuchos::ParameterList>
430 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
431 if (is_null(validPL)) {
432 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
437 "The relative residual tolerance that needs to be achieved by the\n" 438 "iterative solver in order for the linear system to be declared converged." );
439 pl->set(
"Maximum Restarts", static_cast<int>(maxRestarts_default_),
440 "The maximum number of restarts allowed for each\n" 441 "set of RHS solved.");
442 pl->set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
443 "The maximum number of block iterations allowed for each\n" 444 "set of RHS solved.");
445 pl->set(
"Num Blocks", static_cast<int>(numBlocks_default_),
446 "The maximum number of blocks allowed in the Krylov subspace\n" 447 "for each set of RHS solved.");
448 pl->set(
"Block Size", static_cast<int>(blockSize_default_),
449 "The number of vectors in each block. This number times the\n" 450 "number of blocks is the total Krylov subspace dimension.");
451 pl->set(
"Adaptive Block Size", static_cast<bool>(adaptiveBlockSize_default_),
452 "Whether the solver manager should adapt the block size\n" 453 "based on the number of RHS to solve.");
454 pl->set(
"Verbosity", static_cast<int>(verbosity_default_),
455 "What type(s) of solver information should be outputted\n" 456 "to the output stream.");
457 pl->set(
"Output Style", static_cast<int>(outputStyle_default_),
458 "What style is used for the solver information outputted\n" 459 "to the output stream.");
460 pl->set(
"Output Frequency", static_cast<int>(outputFreq_default_),
461 "How often convergence information should be outputted\n" 462 "to the output stream.");
463 pl->set(
"Output Stream", Teuchos::rcp(outputStream_default_,
false),
464 "A reference-counted pointer to the output stream where all\n" 465 "solver output is sent.");
466 pl->set(
"Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
467 "When convergence information is printed, only show the maximum\n" 468 "relative residual norm when the block size is greater than one.");
469 pl->set(
"Flexible Gmres", static_cast<bool>(flexibleGmres_default_),
470 "Whether the solver manager should use the flexible variant\n" 472 pl->set(
"Explicit Residual Test", static_cast<bool>(expResTest_default_),
473 "Whether the explicitly computed residual should be used in the convergence test.");
474 pl->set(
"Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
475 "The type of scaling used in the implicit residual convergence test.");
476 pl->set(
"Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
477 "The type of scaling used in the explicit residual convergence test.");
478 pl->set(
"Timer Label", static_cast<const char *>(label_default_),
479 "The string to use as a prefix for the timer labels.");
480 pl->set(
"Orthogonalization", static_cast<const char *>(orthoType_default_),
481 "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
483 "The constant used by DGKS orthogonalization to determine\n" 484 "whether another step of classical Gram-Schmidt is necessary.");
491 template<
class ScalarType,
class MV,
class OP>
496 if (params_ == Teuchos::null) {
497 params_ = Teuchos::rcp(
new Teuchos::ParameterList(*getValidParameters()) );
500 params->validateParameters(*getValidParameters());
504 if (params->isParameter(
"Maximum Restarts")) {
505 maxRestarts_ = params->get(
"Maximum Restarts",maxRestarts_default_);
508 params_->set(
"Maximum Restarts", maxRestarts_);
512 if (params->isParameter(
"Maximum Iterations")) {
513 maxIters_ = params->get(
"Maximum Iterations",maxIters_default_);
516 params_->set(
"Maximum Iterations", maxIters_);
517 if (maxIterTest_!=Teuchos::null)
518 maxIterTest_->setMaxIters( maxIters_ );
522 if (params->isParameter(
"Block Size")) {
523 blockSize_ = params->get(
"Block Size",blockSize_default_);
524 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
525 "Belos::BlockGmresSolMgr: \"Block Size\" must be strictly positive.");
528 params_->set(
"Block Size", blockSize_);
532 if (params->isParameter(
"Adaptive Block Size")) {
533 adaptiveBlockSize_ = params->get(
"Adaptive Block Size",adaptiveBlockSize_default_);
536 params_->set(
"Adaptive Block Size", adaptiveBlockSize_);
540 if (params->isParameter(
"Num Blocks")) {
541 numBlocks_ = params->get(
"Num Blocks",numBlocks_default_);
542 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
543 "Belos::BlockGmresSolMgr: \"Num Blocks\" must be strictly positive.");
546 params_->set(
"Num Blocks", numBlocks_);
550 if (params->isParameter(
"Timer Label")) {
551 std::string tempLabel = params->get(
"Timer Label", label_default_);
554 if (tempLabel != label_) {
556 params_->set(
"Timer Label", label_);
557 std::string solveLabel = label_ +
": BlockGmresSolMgr total solve time";
558 #ifdef BELOS_TEUCHOS_TIME_MONITOR 559 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
561 if (ortho_ != Teuchos::null) {
562 ortho_->setLabel( label_ );
568 if (params->isParameter(
"Flexible Gmres")) {
569 isFlexible_ = Teuchos::getParameter<bool>(*params,
"Flexible Gmres");
570 params_->set(
"Flexible Gmres", isFlexible_);
571 if (isFlexible_ && expResTest_) {
578 if (params->isParameter(
"Verbosity")) {
579 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
580 verbosity_ = params->get(
"Verbosity", verbosity_default_);
582 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
586 params_->set(
"Verbosity", verbosity_);
587 if (printer_ != Teuchos::null)
588 printer_->setVerbosity(verbosity_);
592 if (params->isParameter(
"Output Style")) {
593 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
594 outputStyle_ = params->get(
"Output Style", outputStyle_default_);
596 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
600 params_->set(
"Output Style", outputStyle_);
601 if (outputTest_ != Teuchos::null) {
607 if (params->isParameter(
"Output Stream")) {
608 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
611 params_->set(
"Output Stream", outputStream_);
612 if (printer_ != Teuchos::null)
613 printer_->setOStream( outputStream_ );
618 if (params->isParameter(
"Output Frequency")) {
619 outputFreq_ = params->get(
"Output Frequency", outputFreq_default_);
623 params_->set(
"Output Frequency", outputFreq_);
624 if (outputTest_ != Teuchos::null)
625 outputTest_->setOutputFrequency( outputFreq_ );
629 if (printer_ == Teuchos::null) {
634 bool changedOrthoType =
false;
635 if (params->isParameter(
"Orthogonalization")) {
636 std::string tempOrthoType = params->get(
"Orthogonalization",orthoType_default_);
637 if (tempOrthoType != orthoType_) {
638 orthoType_ = tempOrthoType;
639 changedOrthoType =
true;
642 params_->set(
"Orthogonalization", orthoType_);
645 if (params->isParameter(
"Orthogonalization Constant")) {
646 if (params->isType<MagnitudeType> (
"Orthogonalization Constant")) {
647 orthoKappa_ = params->get (
"Orthogonalization Constant",
651 orthoKappa_ = params->get (
"Orthogonalization Constant",
656 params_->set(
"Orthogonalization Constant",orthoKappa_);
657 if (orthoType_==
"DGKS") {
658 if (orthoKappa_ > 0 && ortho_ != Teuchos::null && !changedOrthoType) {
665 if (ortho_ == Teuchos::null || changedOrthoType) {
667 Teuchos::RCP<Teuchos::ParameterList> paramsOrtho;
668 if (orthoType_==
"DGKS" && orthoKappa_ > 0) {
669 paramsOrtho->set (
"depTol", orthoKappa_ );
672 ortho_ = factory.
makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
676 if (params->isParameter(
"Convergence Tolerance")) {
677 if (params->isType<MagnitudeType> (
"Convergence Tolerance")) {
678 convtol_ = params->get (
"Convergence Tolerance",
686 params_->set(
"Convergence Tolerance", convtol_);
687 if (impConvTest_ != Teuchos::null)
688 impConvTest_->setTolerance( convtol_ );
689 if (expConvTest_ != Teuchos::null)
690 expConvTest_->setTolerance( convtol_ );
694 if (params->isParameter(
"Implicit Residual Scaling")) {
695 std::string tempImpResScale = Teuchos::getParameter<std::string>( *params,
"Implicit Residual Scaling" );
698 if (impResScale_ != tempImpResScale) {
700 impResScale_ = tempImpResScale;
703 params_->set(
"Implicit Residual Scaling", impResScale_);
704 if (impConvTest_ != Teuchos::null) {
708 catch (std::exception& e) {
716 if (params->isParameter(
"Explicit Residual Scaling")) {
717 std::string tempExpResScale = Teuchos::getParameter<std::string>( *params,
"Explicit Residual Scaling" );
720 if (expResScale_ != tempExpResScale) {
722 expResScale_ = tempExpResScale;
725 params_->set(
"Explicit Residual Scaling", expResScale_);
726 if (expConvTest_ != Teuchos::null) {
730 catch (std::exception& e) {
738 if (params->isParameter(
"Explicit Residual Test")) {
739 expResTest_ = Teuchos::getParameter<bool>( *params,
"Explicit Residual Test" );
742 params_->set(
"Explicit Residual Test", expResTest_);
743 if (expConvTest_ == Teuchos::null) {
748 if (params->isParameter(
"Show Maximum Residual Norm Only")) {
749 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,
"Show Maximum Residual Norm Only");
752 params_->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
753 if (impConvTest_ != Teuchos::null)
754 impConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
755 if (expConvTest_ != Teuchos::null)
756 expConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
761 if (timerSolve_ == Teuchos::null) {
762 std::string solveLabel = label_ +
": BlockGmresSolMgr total solve time";
763 #ifdef BELOS_TEUCHOS_TIME_MONITOR 764 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
773 template<
class ScalarType,
class MV,
class OP>
786 if (isFlexible_ && Teuchos::is_null(problem_->getRightPrec())) {
788 params_->set(
"Flexible Gmres", isFlexible_);
792 "Belos::BlockGmresSolMgr::solve(): Linear problem has a left preconditioner, not a right preconditioner, which is required for flexible GMRES.");
797 if (!Teuchos::is_null(problem_->getLeftPrec()) && !isFlexible_) {
804 Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
805 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_ ) );
807 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
808 impConvTest_ = tmpImpConvTest;
811 Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
812 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_ ) );
813 tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit,
Belos::TwoNorm );
815 tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
816 expConvTest_ = tmpExpConvTest;
819 convTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
825 Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
826 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_ ) );
828 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
829 impConvTest_ = tmpImpConvTest;
834 Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
835 Teuchos::rcp(
new StatusTestImpResNorm_t( convtol_ ) );
837 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
838 impConvTest_ = tmpImpConvTest;
842 expConvTest_ = impConvTest_;
843 convTest_ = impConvTest_;
847 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
850 if (nonnull(debugStatusTest_) ) {
852 Teuchos::rcp_dynamic_cast<StatusTestCombo_t>(sTest_)->addStatusTest( debugStatusTest_ );
857 StatusTestOutputFactory<ScalarType,MV,OP> stoFactory( outputStyle_ );
861 std::string solverDesc =
" Block Gmres ";
863 solverDesc =
"Flexible" + solverDesc;
864 outputTest_->setSolverDesc( solverDesc );
872 template<
class ScalarType,
class MV,
class OP>
877 debugStatusTest_ = debugStatusTest;
882 template<
class ScalarType,
class MV,
class OP>
889 setParameters(Teuchos::parameterList(*getValidParameters()));
892 Teuchos::BLAS<int,ScalarType> blas;
895 "Belos::BlockGmresSolMgr::solve(): Linear problem is not a valid object.");
898 "Belos::BlockGmresSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
900 if (!isSTSet_ || (!expResTest_ && !Teuchos::is_null(problem_->getLeftPrec())) ) {
902 "Belos::BlockGmresSolMgr::solve(): Linear problem and requested status tests are incompatible.");
907 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
908 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
910 std::vector<int> currIdx;
913 if ( adaptiveBlockSize_ ) {
914 blockSize_ = numCurrRHS;
915 currIdx.resize( numCurrRHS );
916 for (
int i=0; i<numCurrRHS; ++i)
917 { currIdx[i] = startPtr+i; }
921 currIdx.resize( blockSize_ );
922 for (
int i=0; i<numCurrRHS; ++i)
923 { currIdx[i] = startPtr+i; }
924 for (
int i=numCurrRHS; i<blockSize_; ++i)
929 problem_->setLSIndex( currIdx );
933 Teuchos::ParameterList plist;
934 plist.set(
"Block Size",blockSize_);
936 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
937 if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {
938 int tmpNumBlocks = 0;
940 tmpNumBlocks = dim / blockSize_;
942 tmpNumBlocks = ( dim - blockSize_) / blockSize_;
944 "Belos::BlockGmresSolMgr::solve(): Warning! Requested Krylov subspace dimension is larger than operator dimension!" 945 << std::endl <<
" The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
946 plist.set(
"Num Blocks",tmpNumBlocks);
949 plist.set(
"Num Blocks",numBlocks_);
952 outputTest_->reset();
953 loaDetected_ =
false;
956 bool isConverged =
true;
961 Teuchos::RCP<GmresIteration<ScalarType,MV,OP> > block_gmres_iter;
970 #ifdef BELOS_TEUCHOS_TIME_MONITOR 971 Teuchos::TimeMonitor slvtimer(*timerSolve_);
974 while ( numRHS2Solve > 0 ) {
977 if (blockSize_*numBlocks_ > dim) {
978 int tmpNumBlocks = 0;
980 tmpNumBlocks = dim / blockSize_;
982 tmpNumBlocks = ( dim - blockSize_) / blockSize_;
983 block_gmres_iter->setSize( blockSize_, tmpNumBlocks );
986 block_gmres_iter->setSize( blockSize_, numBlocks_ );
989 block_gmres_iter->resetNumIters();
992 outputTest_->resetNumCalls();
995 Teuchos::RCP<MV> V_0;
998 if (currIdx[blockSize_-1] == -1) {
999 V_0 = MVT::Clone( *(problem_->getInitResVec()), blockSize_ );
1000 problem_->computeCurrResVec( &*V_0 );
1003 V_0 = MVT::CloneCopy( *(problem_->getInitResVec()), currIdx );
1008 if (currIdx[blockSize_-1] == -1) {
1009 V_0 = MVT::Clone( *(problem_->getInitPrecResVec()), blockSize_ );
1010 problem_->computeCurrPrecResVec( &*V_0 );
1013 V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1018 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z_0 =
1019 Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( blockSize_, blockSize_ ) );
1022 int rank = ortho_->normalize( *V_0, z_0 );
1024 "Belos::BlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors.");
1031 block_gmres_iter->initializeGmres(newstate);
1032 int numRestarts = 0;
1037 block_gmres_iter->iterate();
1044 if ( convTest_->getStatus() ==
Passed ) {
1045 if ( expConvTest_->getLOADetected() ) {
1047 loaDetected_ =
true;
1049 "Belos::BlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
1050 isConverged =
false;
1059 else if ( maxIterTest_->getStatus() ==
Passed ) {
1061 isConverged =
false;
1069 else if ( block_gmres_iter->getCurSubspaceDim() == block_gmres_iter->getMaxSubspaceDim() ) {
1071 if ( numRestarts >= maxRestarts_ ) {
1072 isConverged =
false;
1077 printer_->stream(
Debug) <<
" Performing restart number " << numRestarts <<
" of " << maxRestarts_ << std::endl << std::endl;
1080 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1083 Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1084 MVT::MvAddMv( 1.0, *curX, 1.0, *update, *curX );
1087 problem_->updateSolution( update,
true );
1094 V_0 = MVT::Clone( *(oldState.
V), blockSize_ );
1096 problem_->computeCurrResVec( &*V_0 );
1098 problem_->computeCurrPrecResVec( &*V_0 );
1101 z_0 = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( blockSize_, blockSize_ ) );
1104 rank = ortho_->normalize( *V_0, z_0 );
1106 "Belos::BlockGmresSolMgr::solve(): Failed to compute initial block of orthonormal vectors after restart.");
1112 block_gmres_iter->initializeGmres(newstate);
1124 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
1125 "Belos::BlockGmresSolMgr::solve(): Invalid return from BlockGmresIter::iterate().");
1130 if (blockSize_ != 1) {
1131 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 1132 << block_gmres_iter->getNumIters() << std::endl
1133 << e.what() << std::endl;
1134 if (convTest_->getStatus() !=
Passed)
1135 isConverged =
false;
1140 block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
1143 sTest_->checkStatus( &*block_gmres_iter );
1144 if (convTest_->getStatus() !=
Passed)
1145 isConverged =
false;
1149 catch (
const std::exception &e) {
1150 printer_->stream(
Errors) <<
"Error! Caught std::exception in BlockGmresIter::iterate() at iteration " 1151 << block_gmres_iter->getNumIters() << std::endl
1152 << e.what() << std::endl;
1161 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1162 Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1164 if (update != Teuchos::null)
1165 MVT::MvAddMv( 1.0, *curX, 1.0, *update, *curX );
1169 if ( !Teuchos::is_null(expConvTest_->getSolution()) ) {
1170 Teuchos::RCP<MV> newX = expConvTest_->getSolution();
1171 Teuchos::RCP<MV> curX = problem_->getCurrLHSVec();
1172 MVT::MvAddMv( 0.0, *newX, 1.0, *newX, *curX );
1175 Teuchos::RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1176 problem_->updateSolution( update,
true );
1181 problem_->setCurrLS();
1184 startPtr += numCurrRHS;
1185 numRHS2Solve -= numCurrRHS;
1186 if ( numRHS2Solve > 0 ) {
1187 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1189 if ( adaptiveBlockSize_ ) {
1190 blockSize_ = numCurrRHS;
1191 currIdx.resize( numCurrRHS );
1192 for (
int i=0; i<numCurrRHS; ++i)
1193 { currIdx[i] = startPtr+i; }
1196 currIdx.resize( blockSize_ );
1197 for (
int i=0; i<numCurrRHS; ++i)
1198 { currIdx[i] = startPtr+i; }
1199 for (
int i=numCurrRHS; i<blockSize_; ++i)
1200 { currIdx[i] = -1; }
1203 problem_->setLSIndex( currIdx );
1206 currIdx.resize( numRHS2Solve );
1217 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1222 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1226 numIters_ = maxIterTest_->getNumIters();
1240 const std::vector<MagnitudeType>* pTestValues = NULL;
1242 pTestValues = expConvTest_->getTestValue();
1243 if (pTestValues == NULL || pTestValues->size() < 1) {
1244 pTestValues = impConvTest_->getTestValue();
1249 pTestValues = impConvTest_->getTestValue();
1251 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1252 "Belos::BlockGmresSolMgr::solve(): The implicit convergence test's " 1253 "getTestValue() method returned NULL. Please report this bug to the " 1254 "Belos developers.");
1255 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1256 "Belos::BlockGmresSolMgr::solve(): The implicit convergence test's " 1257 "getTestValue() method returned a vector of length zero. Please report " 1258 "this bug to the Belos developers.");
1263 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1266 if (!isConverged || loaDetected_) {
1273 template<
class ScalarType,
class MV,
class OP>
1276 std::ostringstream out;
1277 out <<
"\"Belos::BlockGmresSolMgr\": {";
1278 if (this->getObjectLabel () !=
"") {
1279 out <<
"Label: " << this->getObjectLabel () <<
", ";
1281 out <<
"Flexible: " << (isFlexible_ ?
"true" :
"false")
1282 <<
", Num Blocks: " << numBlocks_
1283 <<
", Maximum Iterations: " << maxIters_
1284 <<
", Maximum Restarts: " << maxRestarts_
1285 <<
", Convergence Tolerance: " << convtol_
1291 template<
class ScalarType,
class MV,
class OP>
1295 const Teuchos::EVerbosityLevel verbLevel)
const 1297 using Teuchos::TypeNameTraits;
1298 using Teuchos::VERB_DEFAULT;
1299 using Teuchos::VERB_NONE;
1300 using Teuchos::VERB_LOW;
1307 const Teuchos::EVerbosityLevel vl =
1308 (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
1310 if (vl != VERB_NONE) {
1311 Teuchos::OSTab tab0 (out);
1313 out <<
"\"Belos::BlockGmresSolMgr\":" << endl;
1314 Teuchos::OSTab tab1 (out);
1315 out <<
"Template parameters:" << endl;
1317 Teuchos::OSTab tab2 (out);
1318 out <<
"ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1319 <<
"MV: " << TypeNameTraits<MV>::name () << endl
1320 <<
"OP: " << TypeNameTraits<OP>::name () << endl;
1322 if (this->getObjectLabel () !=
"") {
1323 out <<
"Label: " << this->getObjectLabel () << endl;
1325 out <<
"Flexible: " << (isFlexible_ ?
"true" :
"false") << endl
1326 <<
"Num Blocks: " << numBlocks_ << endl
1327 <<
"Maximum Iterations: " << maxIters_ << endl
1328 <<
"Maximum Restarts: " << maxRestarts_ << endl
1329 <<
"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.