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;