42 #ifndef BELOS_GMRES_POLY_SOLMGR_HPP 43 #define BELOS_GMRES_POLY_SOLMGR_HPP 67 #include "Teuchos_BLAS.hpp" 68 #include "Teuchos_as.hpp" 69 #ifdef BELOS_TEUCHOS_TIME_MONITOR 70 #include "Teuchos_TimeMonitor.hpp" 154 template<
class ScalarType,
class MV,
class OP>
159 typedef Teuchos::ScalarTraits<ScalarType>
STS;
160 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
MagnitudeType;
161 typedef Teuchos::ScalarTraits<MagnitudeType>
MT;
189 const Teuchos::RCP<Teuchos::ParameterList> &pl );
217 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
240 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms );
303 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >
problem_;
310 Teuchos::RCP<StatusTest<ScalarType,MV,OP> >
sTest_;
317 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> >
ortho_;
353 Teuchos::RCP<Teuchos::SerialDenseVector<int, ScalarType> >
poly_r0_;
354 Teuchos::RCP<Belos::GmresPolyOp<ScalarType, MV, OP> >
poly_Op_;
366 mutable Teuchos::RCP<const Teuchos::ParameterList>
validPL_;
371 template<
class ScalarType,
class MV,
class OP>
375 template<
class ScalarType,
class MV,
class OP>
379 template<
class ScalarType,
class MV,
class OP>
382 -Teuchos::ScalarTraits<MagnitudeType>::one();
384 template<
class ScalarType,
class MV,
class OP>
387 template<
class ScalarType,
class MV,
class OP>
390 template<
class ScalarType,
class MV,
class OP>
393 template<
class ScalarType,
class MV,
class OP>
396 template<
class ScalarType,
class MV,
class OP>
399 template<
class ScalarType,
class MV,
class OP>
402 template<
class ScalarType,
class MV,
class OP>
405 template<
class ScalarType,
class MV,
class OP>
408 template<
class ScalarType,
class MV,
class OP>
411 template<
class ScalarType,
class MV,
class OP>
414 template<
class ScalarType,
class MV,
class OP>
417 template<
class ScalarType,
class MV,
class OP>
420 template<
class ScalarType,
class MV,
class OP>
423 template<
class ScalarType,
class MV,
class OP>
426 template<
class ScalarType,
class MV,
class OP>
427 const Teuchos::RCP<std::ostream>
431 template<
class ScalarType,
class MV,
class OP>
433 outputStream_ (outputStream_default_),
434 polytol_ (polytol_default_),
435 convtol_ (convtol_default_),
436 orthoKappa_ (orthoKappa_default_),
437 maxDegree_ (maxDegree_default_),
438 maxRestarts_ (maxRestarts_default_),
439 maxIters_ (maxIters_default_),
441 blockSize_ (blockSize_default_),
442 numBlocks_ (numBlocks_default_),
443 verbosity_ (verbosity_default_),
444 outputStyle_ (outputStyle_default_),
445 outputFreq_ (outputFreq_default_),
446 strictConvTol_ (strictConvTol_default_),
447 showMaxResNormOnly_ (showMaxResNormOnly_default_),
448 orthoType_ (orthoType_default_),
449 impResScale_ (impResScale_default_),
450 expResScale_ (expResScale_default_),
452 label_ (label_default_),
453 isPolyBuilt_ (false),
461 template<
class ScalarType,
class MV,
class OP>
464 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
466 outputStream_ (outputStream_default_),
467 polytol_ (polytol_default_),
468 convtol_ (convtol_default_),
469 orthoKappa_ (orthoKappa_default_),
470 maxDegree_ (maxDegree_default_),
471 maxRestarts_ (maxRestarts_default_),
472 maxIters_ (maxIters_default_),
474 blockSize_ (blockSize_default_),
475 numBlocks_ (numBlocks_default_),
476 verbosity_ (verbosity_default_),
477 outputStyle_ (outputStyle_default_),
478 outputFreq_ (outputFreq_default_),
479 strictConvTol_ (strictConvTol_default_),
480 showMaxResNormOnly_ (showMaxResNormOnly_default_),
481 orthoType_ (orthoType_default_),
482 impResScale_ (impResScale_default_),
483 expResScale_ (expResScale_default_),
485 label_ (label_default_),
486 isPolyBuilt_ (false),
492 TEUCHOS_TEST_FOR_EXCEPTION(
493 problem_.is_null (), std::invalid_argument,
494 "Belos::GmresPolySolMgr: The given linear problem is null. " 495 "Please call this constructor with a nonnull LinearProblem argument, " 496 "or call the constructor that does not take a LinearProblem.");
500 if (! pl.is_null ()) {
506 template<
class ScalarType,
class MV,
class OP>
507 Teuchos::RCP<const Teuchos::ParameterList>
510 if (validPL_.is_null ()) {
511 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList ();
512 pl->set(
"Polynomial Tolerance", polytol_default_,
513 "The relative residual tolerance that used to construct the GMRES polynomial.");
514 pl->set(
"Maximum Degree", maxDegree_default_,
515 "The maximum degree allowed for any GMRES polynomial.");
516 pl->set(
"Convergence Tolerance", convtol_default_,
517 "The relative residual tolerance that needs to be achieved by the\n" 518 "iterative solver in order for the linear system to be declared converged." );
519 pl->set(
"Maximum Restarts", maxRestarts_default_,
520 "The maximum number of restarts allowed for each\n" 521 "set of RHS solved.");
522 pl->set(
"Maximum Iterations", maxIters_default_,
523 "The maximum number of block iterations allowed for each\n" 524 "set of RHS solved.");
525 pl->set(
"Num Blocks", numBlocks_default_,
526 "The maximum number of blocks allowed in the Krylov subspace\n" 527 "for each set of RHS solved.");
528 pl->set(
"Block Size", blockSize_default_,
529 "The number of vectors in each block. This number times the\n" 530 "number of blocks is the total Krylov subspace dimension.");
531 pl->set(
"Verbosity", verbosity_default_,
532 "What type(s) of solver information should be outputted\n" 533 "to the output stream.");
534 pl->set(
"Output Style", outputStyle_default_,
535 "What style is used for the solver information outputted\n" 536 "to the output stream.");
537 pl->set(
"Output Frequency", outputFreq_default_,
538 "How often convergence information should be outputted\n" 539 "to the output stream.");
540 pl->set(
"Output Stream", outputStream_default_,
541 "A reference-counted pointer to the output stream where all\n" 542 "solver output is sent.");
543 pl->set(
"Strict Convergence", strictConvTol_default_,
544 "After polynomial is applied, whether solver should try to achieve\n" 545 "the relative residual tolerance.");
546 pl->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
547 "When convergence information is printed, only show the maximum\n" 548 "relative residual norm when the block size is greater than one.");
549 pl->set(
"Implicit Residual Scaling", impResScale_default_,
550 "The type of scaling used in the implicit residual convergence test.");
551 pl->set(
"Explicit Residual Scaling", expResScale_default_,
552 "The type of scaling used in the explicit residual convergence test.");
553 pl->set(
"Timer Label", label_default_,
554 "The string to use as a prefix for the timer labels.");
556 pl->set(
"Orthogonalization", orthoType_default_,
557 "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
558 pl->set(
"Orthogonalization Constant",orthoKappa_default_,
559 "The constant used by DGKS orthogonalization to determine\n" 560 "whether another step of classical Gram-Schmidt is necessary.");
567 template<
class ScalarType,
class MV,
class OP>
572 if (params_.is_null ()) {
573 params_ = Teuchos::parameterList (*getValidParameters ());
576 params->validateParameters (*getValidParameters ());
580 if (params->isParameter(
"Maximum Degree")) {
581 maxDegree_ = params->get(
"Maximum Degree",maxDegree_default_);
584 params_->set(
"Maximum Degree", maxDegree_);
588 if (params->isParameter(
"Maximum Restarts")) {
589 maxRestarts_ = params->get(
"Maximum Restarts",maxRestarts_default_);
592 params_->set(
"Maximum Restarts", maxRestarts_);
596 if (params->isParameter(
"Maximum Iterations")) {
597 maxIters_ = params->get(
"Maximum Iterations",maxIters_default_);
600 params_->set(
"Maximum Iterations", maxIters_);
601 if (maxIterTest_!=Teuchos::null)
602 maxIterTest_->setMaxIters( maxIters_ );
606 if (params->isParameter(
"Block Size")) {
607 blockSize_ = params->get(
"Block Size",blockSize_default_);
608 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
609 "Belos::GmresPolySolMgr: \"Block Size\" must be strictly positive.");
612 params_->set(
"Block Size", blockSize_);
616 if (params->isParameter(
"Num Blocks")) {
617 numBlocks_ = params->get(
"Num Blocks",numBlocks_default_);
618 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
619 "Belos::GmresPolySolMgr: \"Num Blocks\" must be strictly positive.");
622 params_->set(
"Num Blocks", numBlocks_);
626 if (params->isParameter(
"Timer Label")) {
627 std::string tempLabel = params->get(
"Timer Label", label_default_);
630 if (tempLabel != label_) {
632 params_->set(
"Timer Label", label_);
633 std::string solveLabel = label_ +
": GmresPolySolMgr total solve time";
634 #ifdef BELOS_TEUCHOS_TIME_MONITOR 635 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
637 std::string polyLabel = label_ +
": GmresPolySolMgr polynomial creation time";
638 #ifdef BELOS_TEUCHOS_TIME_MONITOR 639 timerPoly_ = Teuchos::TimeMonitor::getNewCounter(polyLabel);
641 if (ortho_ != Teuchos::null) {
642 ortho_->setLabel( label_ );
648 if (params->isParameter(
"Orthogonalization")) {
649 std::string tempOrthoType = params->get(
"Orthogonalization",orthoType_default_);
650 TEUCHOS_TEST_FOR_EXCEPTION( tempOrthoType !=
"DGKS" && tempOrthoType !=
"ICGS" && tempOrthoType !=
"IMGS",
651 std::invalid_argument,
652 "Belos::GmresPolySolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
653 if (tempOrthoType != orthoType_) {
654 orthoType_ = tempOrthoType;
656 if (orthoType_==
"DGKS") {
657 if (orthoKappa_ <= 0) {
665 else if (orthoType_==
"ICGS") {
668 else if (orthoType_==
"IMGS") {
675 if (params->isParameter(
"Orthogonalization Constant")) {
676 orthoKappa_ = params->get(
"Orthogonalization Constant",orthoKappa_default_);
679 params_->set(
"Orthogonalization Constant",orthoKappa_);
680 if (orthoType_==
"DGKS") {
681 if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
688 if (params->isParameter(
"Verbosity")) {
689 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
690 verbosity_ = params->get(
"Verbosity", verbosity_default_);
692 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
696 params_->set(
"Verbosity", verbosity_);
697 if (printer_ != Teuchos::null)
698 printer_->setVerbosity(verbosity_);
702 if (params->isParameter(
"Output Style")) {
703 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
704 outputStyle_ = params->get(
"Output Style", outputStyle_default_);
706 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
710 params_->set(
"Output Style", outputStyle_);
711 if (outputTest_ != Teuchos::null) {
717 if (params->isParameter(
"Output Stream")) {
718 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
721 params_->set(
"Output Stream", outputStream_);
722 if (printer_ != Teuchos::null)
723 printer_->setOStream( outputStream_ );
728 if (params->isParameter(
"Output Frequency")) {
729 outputFreq_ = params->get(
"Output Frequency", outputFreq_default_);
733 params_->set(
"Output Frequency", outputFreq_);
734 if (outputTest_ != Teuchos::null)
735 outputTest_->setOutputFrequency( outputFreq_ );
739 if (printer_ == Teuchos::null) {
748 if (params->isParameter(
"Polynomial Tolerance")) {
749 polytol_ = params->get(
"Polynomial Tolerance",polytol_default_);
752 params_->set(
"Polynomial Tolerance", polytol_);
756 if (params->isParameter(
"Convergence Tolerance")) {
757 convtol_ = params->get(
"Convergence Tolerance",convtol_default_);
760 params_->set(
"Convergence Tolerance", convtol_);
761 if (impConvTest_ != Teuchos::null)
762 impConvTest_->setTolerance( convtol_ );
763 if (expConvTest_ != Teuchos::null)
764 expConvTest_->setTolerance( convtol_ );
768 if (params->isParameter(
"Strict Convergence")) {
769 strictConvTol_ = params->get(
"Strict Convergence",strictConvTol_default_);
772 params_->set(
"Strict Convergence", strictConvTol_);
776 if (params->isParameter(
"Implicit Residual Scaling")) {
777 std::string tempImpResScale = Teuchos::getParameter<std::string>( *params,
"Implicit Residual Scaling" );
780 if (impResScale_ != tempImpResScale) {
782 impResScale_ = tempImpResScale;
785 params_->set(
"Implicit Residual Scaling", impResScale_);
786 if (impConvTest_ != Teuchos::null) {
790 catch (std::exception& e) {
798 if (params->isParameter(
"Explicit Residual Scaling")) {
799 std::string tempExpResScale = Teuchos::getParameter<std::string>( *params,
"Explicit Residual Scaling" );
802 if (expResScale_ != tempExpResScale) {
804 expResScale_ = tempExpResScale;
807 params_->set(
"Explicit Residual Scaling", expResScale_);
808 if (expConvTest_ != Teuchos::null) {
812 catch (std::exception& e) {
821 if (params->isParameter(
"Show Maximum Residual Norm Only")) {
822 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,
"Show Maximum Residual Norm Only");
825 params_->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
826 if (impConvTest_ != Teuchos::null)
827 impConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
828 if (expConvTest_ != Teuchos::null)
829 expConvTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
833 if (ortho_ == Teuchos::null) {
834 if (orthoType_==
"DGKS") {
835 if (orthoKappa_ <= 0) {
843 else if (orthoType_==
"ICGS") {
846 else if (orthoType_==
"IMGS") {
850 TEUCHOS_TEST_FOR_EXCEPTION(orthoType_!=
"ICGS"&&orthoType_!=
"DGKS"&&orthoType_!=
"IMGS",std::logic_error,
851 "Belos::GmresPolySolMgr(): Invalid orthogonalization type.");
856 if (timerSolve_ == Teuchos::null) {
857 std::string solveLabel = label_ +
": GmresPolySolMgr total solve time";
858 #ifdef BELOS_TEUCHOS_TIME_MONITOR 859 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
863 if (timerPoly_ == Teuchos::null) {
864 std::string polyLabel = label_ +
": GmresPolySolMgr polynomial creation time";
865 #ifdef BELOS_TEUCHOS_TIME_MONITOR 866 timerPoly_ = Teuchos::TimeMonitor::getNewCounter(polyLabel);
875 template<
class ScalarType,
class MV,
class OP>
887 if (!Teuchos::is_null(problem_->getLeftPrec())) {
894 Teuchos::RCP<StatusTestGenResNorm_t> tmpImpConvTest =
895 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_ ) );
897 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
898 impConvTest_ = tmpImpConvTest;
901 Teuchos::RCP<StatusTestGenResNorm_t> tmpExpConvTest =
902 Teuchos::rcp(
new StatusTestGenResNorm_t( convtol_ ) );
903 tmpExpConvTest->defineResForm( StatusTestGenResNorm_t::Explicit,
Belos::TwoNorm );
905 tmpExpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
906 expConvTest_ = tmpExpConvTest;
909 convTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::SEQ, impConvTest_, expConvTest_ ) );
915 Teuchos::RCP<StatusTestImpResNorm_t> tmpImpConvTest =
916 Teuchos::rcp(
new StatusTestImpResNorm_t( convtol_ ) );
918 tmpImpConvTest->setShowMaxResNormOnly( showMaxResNormOnly_ );
919 impConvTest_ = tmpImpConvTest;
922 expConvTest_ = impConvTest_;
923 convTest_ = impConvTest_;
926 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
934 std::string solverDesc =
" Gmres Polynomial ";
935 outputTest_->setSolverDesc( solverDesc );
944 template<
class ScalarType,
class MV,
class OP>
948 Teuchos::RCP<MV> newX = MVT::Clone( *(problem_->getLHS()), 1 );
949 Teuchos::RCP<MV> newB = MVT::Clone( *(problem_->getRHS()), 1 );
950 MVT::MvInit( *newX, STS::zero() );
951 MVT::MvRandom( *newB );
952 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > newProblem =
954 newProblem->setLeftPrec( problem_->getLeftPrec() );
955 newProblem->setRightPrec( problem_->getRightPrec() );
956 newProblem->setLabel(
"Belos GMRES Poly Generation");
957 newProblem->setProblem();
958 std::vector<int> idx(1,0);
959 newProblem->setLSIndex( idx );
962 Teuchos::ParameterList polyList;
965 polyList.set(
"Num Blocks",maxDegree_);
966 polyList.set(
"Block Size",1);
967 polyList.set(
"Keep Hessenberg",
true);
970 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxItrTst =
974 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTst =
979 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > polyTest =
983 Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > gmres_iter;
987 Teuchos::RCP<MV> V_0 = MVT::Clone( *(newProblem->getRHS()), 1 );
988 newProblem->computeCurrPrecResVec( &*V_0 );
991 poly_r0_ = Teuchos::rcp(
new Teuchos::SerialDenseVector<int,ScalarType>(1) );
994 int rank = ortho_->normalize( *V_0, poly_r0_ );
996 "Belos::GmresPolySolMgr::generatePoly(): Failed to compute initial block of orthonormal vectors for polynomial generation.");
1001 newstate.
z = poly_r0_;
1003 gmres_iter->initializeGmres(newstate);
1006 bool polyConverged =
false;
1008 gmres_iter->iterate();
1011 if ( convTst->getStatus() ==
Passed ) {
1013 polyConverged =
true;
1018 gmres_iter->updateLSQR( gmres_iter->getCurSubspaceDim() );
1021 polyTest->checkStatus( &*gmres_iter );
1022 if (convTst->getStatus() ==
Passed)
1023 polyConverged =
true;
1025 catch (std::exception e) {
1026 printer_->stream(
Errors) <<
"Error! Caught exception in BlockGmresIter::iterate() at iteration " 1027 << gmres_iter->getNumIters() << std::endl
1028 << e.what() << std::endl;
1035 (void) polyConverged;
1038 Teuchos::RCP<MV> currX = gmres_iter->getCurrentUpdate();
1044 poly_dim_ = gmresState.
curDim;
1045 if (poly_dim_ == 0) {
1051 poly_y_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::Copy, *gmresState.
z, poly_dim_, 1 ) );
1052 poly_H_ = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>( *gmresState.
H ) );
1056 const ScalarType one = STS::one ();
1057 Teuchos::BLAS<int,ScalarType> blas;
1058 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
1059 Teuchos::NON_UNIT_DIAG, poly_dim_, 1, one,
1060 gmresState.
R->values(), gmresState.
R->stride(),
1061 poly_y_->values(), poly_y_->stride() );
1065 poly_Op_ = Teuchos::rcp(
1072 template<
class ScalarType,
class MV,
class OP>
1077 typedef Teuchos::SerialDenseMatrix<int, ScalarType> SDM;
1084 setParameters (Teuchos::parameterList (*getValidParameters ()));
1087 TEUCHOS_TEST_FOR_EXCEPTION(
1089 "Belos::GmresPolySolMgr::solve: The linear problem has not been set yet, " 1090 "or was set to null. Please call setProblem() with a nonnull input before " 1091 "calling solve().");
1093 TEUCHOS_TEST_FOR_EXCEPTION(
1095 "Belos::GmresPolySolMgr::solve: The linear problem is not ready. Please " 1096 "call setProblem() on the LinearProblem object before calling solve().");
1098 if (! isSTSet_ || (! expResTest_ && ! problem_->getLeftPrec ().is_null ())) {
1103 const bool check = checkStatusTest ();
1104 TEUCHOS_TEST_FOR_EXCEPTION(
1106 "Belos::GmresPolySolMgr::solve: Linear problem and requested status " 1107 "tests are incompatible.");
1112 if (! isPolyBuilt_) {
1113 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1114 Teuchos::TimeMonitor slvtimer(*timerPoly_);
1116 isPolyBuilt_ = generatePoly();
1118 "Belos::GmresPolySolMgr::generatePoly(): Failed to generate a non-trivial polynomial, reduce polynomial tolerance.");
1120 "Belos::GmresPolySolMgr::generatePoly(): Failed to generate polynomial that satisfied requirements.");
1124 bool isConverged =
true;
1128 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1129 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1133 poly_Op_->Apply( *problem_->getRHS(), *problem_->getLHS() );
1136 problem_->setProblem ();
1139 if (strictConvTol_) {
1143 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1144 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1152 std::vector<int> currIdx (blockSize_);
1153 for (
int i = 0; i < numCurrRHS; ++i) {
1154 currIdx[i] = startPtr+i;
1156 for (
int i = numCurrRHS; i < blockSize_; ++i) {
1161 problem_->setLSIndex (currIdx);
1165 Teuchos::ParameterList plist;
1166 plist.set (
"Block Size", blockSize_);
1168 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1169 if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {
1170 ptrdiff_t tmpNumBlocks = 0;
1171 if (blockSize_ == 1) {
1172 tmpNumBlocks = dim / blockSize_;
1174 tmpNumBlocks = ( dim - blockSize_) / blockSize_;
1177 <<
"Warning! Requested Krylov subspace dimension is larger than " 1178 <<
"operator dimension!" << std::endl <<
"The maximum number of " 1179 <<
"blocks allowed for the Krylov subspace will be adjusted to " 1180 << tmpNumBlocks << std::endl;
1181 plist.set (
"Num Blocks", Teuchos::asSafe<int> (tmpNumBlocks));
1184 plist.set (
"Num Blocks", numBlocks_);
1188 outputTest_->reset ();
1189 loaDetected_ =
false;
1199 RCP<GmresIteration<ScalarType,MV,OP> > block_gmres_iter =
1201 outputTest_, ortho_, plist));
1204 while (numRHS2Solve > 0) {
1207 if (blockSize_*numBlocks_ > dim) {
1208 int tmpNumBlocks = 0;
1209 if (blockSize_ == 1) {
1211 tmpNumBlocks = dim / blockSize_;
1214 tmpNumBlocks = (dim - blockSize_) / blockSize_;
1216 block_gmres_iter->setSize (blockSize_, tmpNumBlocks);
1219 block_gmres_iter->setSize (blockSize_, numBlocks_);
1223 block_gmres_iter->resetNumIters ();
1226 outputTest_->resetNumCalls ();
1229 RCP<MV> V_0 = MVT::CloneCopy (*(problem_->getInitPrecResVec ()), currIdx);
1232 RCP<SDM> z_0 = rcp (
new SDM (blockSize_, blockSize_));
1235 int rank = ortho_->normalize (*V_0, z_0);
1236 TEUCHOS_TEST_FOR_EXCEPTION(
1238 "Belos::GmresPolySolMgr::solve: Failed to compute initial block of " 1239 "orthonormal vectors.");
1246 block_gmres_iter->initializeGmres(newstate);
1247 int numRestarts = 0;
1251 block_gmres_iter->iterate();
1254 if ( convTest_->getStatus() ==
Passed ) {
1255 if ( expConvTest_->getLOADetected() ) {
1257 loaDetected_ =
true;
1258 isConverged =
false;
1264 else if ( maxIterTest_->getStatus() ==
Passed ) {
1266 isConverged =
false;
1270 else if (block_gmres_iter->getCurSubspaceDim () ==
1271 block_gmres_iter->getMaxSubspaceDim ()) {
1272 if (numRestarts >= maxRestarts_) {
1273 isConverged =
false;
1278 printer_->stream(
Debug)
1279 <<
" Performing restart number " << numRestarts <<
" of " 1280 << maxRestarts_ << std::endl << std::endl;
1283 RCP<MV> update = block_gmres_iter->getCurrentUpdate();
1284 problem_->updateSolution (update,
true);
1294 RCP<MV> VV_0 = MVT::Clone (*(oldState.
V), blockSize_);
1295 problem_->computeCurrPrecResVec (&*VV_0);
1301 RCP<SDM> zz_0 = rcp (
new SDM (blockSize_, blockSize_));
1304 const int theRank = ortho_->normalize( *VV_0, zz_0 );
1305 TEUCHOS_TEST_FOR_EXCEPTION(
1307 "Belos::GmresPolySolMgr::solve: Failed to compute initial " 1308 "block of orthonormal vectors after restart.");
1312 theNewState.
V = VV_0;
1313 theNewState.
z = zz_0;
1315 block_gmres_iter->initializeGmres (theNewState);
1323 TEUCHOS_TEST_FOR_EXCEPTION(
1324 true, std::logic_error,
1325 "Belos::GmresPolySolMgr::solve: Invalid return from " 1326 "BlockGmresIter::iterate(). Please report this bug " 1327 "to the Belos developers.");
1332 if (blockSize_ != 1) {
1334 <<
"Error! Caught std::exception in BlockGmresIter::iterate() " 1335 <<
"at iteration " << block_gmres_iter->getNumIters()
1336 << std::endl << e.what() << std::endl;
1337 if (convTest_->getStatus() !=
Passed) {
1338 isConverged =
false;
1345 block_gmres_iter->updateLSQR (block_gmres_iter->getCurSubspaceDim ());
1349 sTest_->checkStatus (&*block_gmres_iter);
1350 if (convTest_->getStatus() !=
Passed) {
1351 isConverged =
false;
1356 catch (
const std::exception &e) {
1358 <<
"Error! Caught std::exception in BlockGmresIter::iterate() " 1359 <<
"at iteration " << block_gmres_iter->getNumIters() << std::endl
1360 << e.what() << std::endl;
1368 if (! Teuchos::is_null (expConvTest_->getSolution ()) ) {
1369 RCP<MV> newX = expConvTest_->getSolution ();
1370 RCP<MV> curX = problem_->getCurrLHSVec ();
1371 MVT::MvAddMv (STS::zero (), *newX, STS::one (), *newX, *curX);
1374 RCP<MV> update = block_gmres_iter->getCurrentUpdate ();
1375 problem_->updateSolution (update,
true);
1379 problem_->setCurrLS ();
1382 startPtr += numCurrRHS;
1383 numRHS2Solve -= numCurrRHS;
1384 if (numRHS2Solve > 0) {
1385 numCurrRHS = (numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1387 currIdx.resize (blockSize_);
1388 for (
int i=0; i<numCurrRHS; ++i) {
1389 currIdx[i] = startPtr+i;
1391 for (
int i=numCurrRHS; i<blockSize_; ++i) {
1396 problem_->setLSIndex( currIdx );
1399 currIdx.resize( numRHS2Solve );
1411 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1416 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1419 if (!isConverged || loaDetected_) {
1426 template<
class ScalarType,
class MV,
class OP>
1429 std::ostringstream out;
1431 out <<
"\"Belos::GmresPolySolMgr\": {" 1432 <<
"ScalarType: " << Teuchos::TypeNameTraits<ScalarType>::name ()
1433 <<
", Ortho Type: " << orthoType_
1434 <<
", Block Size: " << blockSize_
1435 <<
", Num Blocks: " << numBlocks_
1436 <<
", Max Restarts: " << maxRestarts_;
1443 #endif // BELOS_GMRES_POLY_SOLMGR_HPP ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Collection of types and exceptions used within the Belos solvers.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Class which manages the output and verbosity of the Belos solvers.
Teuchos::RCP< const Teuchos::ParameterList > validPL_
Cached default (valid) parameters.
Teuchos::RCP< OutputManager< ScalarType > > printer_
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
virtual ~GmresPolySolMgr()
Destructor.
Defines the GMRES polynomial operator hybrid-GMRES iterative linear solver.
Teuchos::RCP< Teuchos::SerialDenseVector< int, ScalarType > > poly_r0_
ScaleType
The type of scaling to use on the residual norm value.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
MagnitudeType orthoKappa_
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...
static const std::string orthoType_default_
Teuchos::RCP< const MV > V
The current Krylov basis.
static const int numBlocks_default_
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > convTest_
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
static const Teuchos::RCP< std::ostream > outputStream_default_
static const int verbosity_default_
A factory class for generating StatusTestOutput objects.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
Teuchos::RCP< Belos::GmresPolyOp< ScalarType, MV, OP > > poly_Op_
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
An implementation of StatusTestResNorm using a family of residual norms.
Structure to contain pointers to GmresIteration state variables.
static const int maxRestarts_default_
static const int maxIters_default_
Belos::StatusTest class for specifying a maximum number of iterations.
Teuchos::RCP< Teuchos::ParameterList > params_
int check(Epetra_CrsGraph &A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1, int NumGlobalNonzeros1, int *MyGlobalElements, bool verbose)
static const bool showMaxResNormOnly_default_
A factory class for generating StatusTestOutput objects.
Iterated Modified Gram-Schmidt (IMGS) implementation of the Belos::OrthoManager class.
void reset(const ResetType type)
Reset the solver.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
Teuchos::RCP< StatusTestResNorm< ScalarType, MV, OP > > impConvTest_
std::string description() const
Method to return description of the hybrid block GMRES solver manager.
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
OperatorTraits< ScalarType, MV, OP > OPT
A Belos::StatusTest class for specifying a maximum number of iterations.
static const int blockSize_default_
ResetType
How to reset the solver.
GmresPolySolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
Pure virtual base class which describes the basic interface for a solver manager. ...
static const int outputFreq_default_
Teuchos::RCP< Teuchos::Time > timerSolve_
Teuchos::ScalarTraits< MagnitudeType > MT
int curDim
The current dimension of the reduction.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get a parameter list containing the valid parameters for this object.
Teuchos::ScalarTraits< ScalarType > STS
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > R
The current upper-triangular matrix from the QR reduction of H.
static const int maxDegree_default_
A linear system to solve, and its associated information.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Class which describes the linear problem to be solved by the iterative solver.
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Set the parameters the solver manager should use to solve the linear problem.
static const std::string expResScale_default_
MultiVecTraits< ScalarType, MV > MVT
static const MagnitudeType orthoKappa_default_
static const bool strictConvTol_default_
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
Hybrid block GMRES iterative linear solver.
Teuchos::RCP< StatusTestResNorm< ScalarType, MV, OP > > expConvTest_
ReturnType
Whether the Belos solve converged for all linear systems.
The Belos::SolverManager is a templated virtual base class that defines the basic interface that any ...
ReturnType solve()
This method performs possibly repeated calls to the underlying linear solver's iterate() routine unti...
Belos's class for applying the GMRES polynomial operator that is used by the hybrid-GMRES linear solv...
Belos concrete class for performing the block GMRES iteration.
Teuchos::RCP< Teuchos::Time > timerPoly_
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
static const MagnitudeType convtol_default_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > poly_H_
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
GmresPolySolMgrPolynomialFailure is thrown when their is a problem generating the GMRES polynomial fo...
GmresPolySolMgrPolynomialFailure(const std::string &what_arg)
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
Belos::StatusTest for specifying an implicit residual norm stopping criteria that checks for loss of ...
static const std::string label_default_
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
A class for extending the status testing capabilities of Belos via logical combinations.
Class which defines basic traits for the operator type.
Teuchos::RCP< std::ostream > outputStream_
Parent class to all Belos exceptions.
GmresPolySolMgr()
Empty constructor for GmresPolySolMgr. This constructor takes no arguments and sets the default value...
Pure virtual base class which augments the basic interface for a Gmres linear solver iteration...
GmresIterationOrthoFailure is thrown when the GmresIteration object is unable to compute independent ...
int getNumIters() const
Get the iteration count for the most recent call to solve().
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > poly_y_
Belos header file which uses auto-configuration information to include necessary C++ headers...
static const std::string impResScale_default_
static const MagnitudeType polytol_default_
bool isLOADetected() const
Return whether a loss of accuracy was detected by this solver during the most current solve...
GmresPolySolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthon...
static const int outputStyle_default_
GmresPolySolMgrOrthoFailure(const std::string &what_arg)
GmresPolySolMgrLinearProblemFailure(const std::string &what_arg)
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.