157 typedef Teuchos::ScalarTraits<ScalarType> SCT;
158 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
159 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
210 const Teuchos::RCP<Teuchos::ParameterList> &pl );
216 Teuchos::RCP<SolverManager<ScalarType, MV, OP> >
clone ()
const override {
230 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters()
const override;
241 Teuchos::Array<Teuchos::RCP<Teuchos::Time> >
getTimers()
const {
242 return Teuchos::tuple(timerSolve_);
272 void setParameters(
const Teuchos::RCP<Teuchos::ParameterList> ¶ms )
override;
313 std::string description()
const override;
320 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
323 Teuchos::RCP<OutputManager<ScalarType> > printer_;
325 Teuchos::RCP<std::ostream> outputStream_;
331 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
334 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
337 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > convTest_;
340 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
343 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
346 Teuchos::RCP<Teuchos::ParameterList> params_;
351 static constexpr int maxIters_default_ = 1000;
352 static constexpr bool adaptiveBlockSize_default_ =
true;
353 static constexpr bool showMaxResNormOnly_default_ =
false;
354 static constexpr bool useSingleReduction_default_ =
false;
355 static constexpr int blockSize_default_ = 1;
358 static constexpr int outputFreq_default_ = -1;
359 static constexpr const char * resNorm_default_ =
"TwoNorm";
360 static constexpr bool foldConvergenceDetectionIntoAllreduce_default_ =
false;
361 static constexpr const char * resScale_default_ =
"Norm of Initial Residual";
362 static constexpr const char * label_default_ =
"Belos";
363 static constexpr const char * orthoType_default_ =
"ICGS";
364 static constexpr bool assertPositiveDefiniteness_default_ =
true;
371 MagnitudeType convtol_;
374 MagnitudeType orthoKappa_;
381 MagnitudeType achievedTol_;
390 int blockSize_, verbosity_, outputStyle_, outputFreq_;
391 bool adaptiveBlockSize_, showMaxResNormOnly_, useSingleReduction_;
392 std::string orthoType_, resScale_;
393 bool assertPositiveDefiniteness_;
394 bool foldConvergenceDetectionIntoAllreduce_;
400 Teuchos::RCP<Teuchos::Time> timerSolve_;
410 outputStream_(Teuchos::rcpFromRef(std::cout)),
413 achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
414 maxIters_(maxIters_default_),
416 blockSize_(blockSize_default_),
417 verbosity_(verbosity_default_),
418 outputStyle_(outputStyle_default_),
419 outputFreq_(outputFreq_default_),
420 adaptiveBlockSize_(adaptiveBlockSize_default_),
421 showMaxResNormOnly_(showMaxResNormOnly_default_),
422 useSingleReduction_(useSingleReduction_default_),
423 orthoType_(orthoType_default_),
424 resScale_(resScale_default_),
425 assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
426 foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
427 label_(label_default_),
436 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
438 outputStream_(Teuchos::rcpFromRef(std::cout)),
441 achievedTol_(Teuchos::ScalarTraits<MagnitudeType>::zero()),
442 maxIters_(maxIters_default_),
444 blockSize_(blockSize_default_),
445 verbosity_(verbosity_default_),
446 outputStyle_(outputStyle_default_),
447 outputFreq_(outputFreq_default_),
448 adaptiveBlockSize_(adaptiveBlockSize_default_),
449 showMaxResNormOnly_(showMaxResNormOnly_default_),
450 useSingleReduction_(useSingleReduction_default_),
451 orthoType_(orthoType_default_),
452 resScale_(resScale_default_),
453 assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
454 foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
455 label_(label_default_),
458 TEUCHOS_TEST_FOR_EXCEPTION(problem_.is_null(), std::invalid_argument,
459 "BlockCGSolMgr's constructor requires a nonnull LinearProblem instance.");
464 if (! pl.is_null()) {
472setParameters (
const Teuchos::RCP<Teuchos::ParameterList> ¶ms)
475 if (params_ == Teuchos::null) {
476 params_ = Teuchos::rcp(
new Teuchos::ParameterList(*getValidParameters()) );
479 params->validateParameters(*getValidParameters());
483 if (params->isParameter(
"Maximum Iterations")) {
484 maxIters_ = params->get(
"Maximum Iterations",maxIters_default_);
487 params_->set(
"Maximum Iterations", maxIters_);
488 if (maxIterTest_!=Teuchos::null)
489 maxIterTest_->setMaxIters( maxIters_ );
493 if (params->isParameter(
"Block Size")) {
494 blockSize_ = params->get(
"Block Size",blockSize_default_);
495 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
496 "Belos::BlockCGSolMgr: \"Block Size\" must be strictly positive.");
499 params_->set(
"Block Size", blockSize_);
503 if (params->isParameter(
"Adaptive Block Size")) {
504 adaptiveBlockSize_ = params->get(
"Adaptive Block Size",adaptiveBlockSize_default_);
507 params_->set(
"Adaptive Block Size", adaptiveBlockSize_);
511 if (params->isParameter(
"Use Single Reduction")) {
512 useSingleReduction_ = params->get(
"Use Single Reduction", useSingleReduction_default_);
515 if (params->isParameter(
"Fold Convergence Detection Into Allreduce")) {
516 foldConvergenceDetectionIntoAllreduce_ = params->get(
"Fold Convergence Detection Into Allreduce",
517 foldConvergenceDetectionIntoAllreduce_default_);
521 if (params->isParameter(
"Timer Label")) {
522 std::string tempLabel = params->get(
"Timer Label", label_default_);
525 if (tempLabel != label_) {
527 params_->set(
"Timer Label", label_);
528 std::string solveLabel = label_ +
": BlockCGSolMgr total solve time";
529#ifdef BELOS_TEUCHOS_TIME_MONITOR
530 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
532 if (ortho_ != Teuchos::null) {
533 ortho_->setLabel( label_ );
539 if (params->isParameter(
"Verbosity")) {
540 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
541 verbosity_ = params->get(
"Verbosity", verbosity_default_);
543 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
547 params_->set(
"Verbosity", verbosity_);
548 if (printer_ != Teuchos::null)
549 printer_->setVerbosity(verbosity_);
553 if (params->isParameter(
"Output Style")) {
554 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
555 outputStyle_ = params->get(
"Output Style", outputStyle_default_);
557 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
561 params_->set(
"Output Style", outputStyle_);
562 outputTest_ = Teuchos::null;
566 if (params->isParameter(
"Output Stream")) {
567 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
570 params_->set(
"Output Stream", outputStream_);
571 if (printer_ != Teuchos::null)
572 printer_->setOStream( outputStream_ );
577 if (params->isParameter(
"Output Frequency")) {
578 outputFreq_ = params->get(
"Output Frequency", outputFreq_default_);
582 params_->set(
"Output Frequency", outputFreq_);
583 if (outputTest_ != Teuchos::null)
584 outputTest_->setOutputFrequency( outputFreq_ );
588 if (printer_ == Teuchos::null) {
593 bool changedOrthoType =
false;
594 if (params->isParameter(
"Orthogonalization")) {
595 std::string tempOrthoType = params->get(
"Orthogonalization",orthoType_default_);
596 if (tempOrthoType != orthoType_) {
597 orthoType_ = tempOrthoType;
598 changedOrthoType =
true;
601 params_->set(
"Orthogonalization", orthoType_);
604 if (params->isParameter(
"Orthogonalization Constant")) {
605 if (params->isType<MagnitudeType> (
"Orthogonalization Constant")) {
606 orthoKappa_ = params->get (
"Orthogonalization Constant",
610 orthoKappa_ = params->get (
"Orthogonalization Constant",
615 params_->set(
"Orthogonalization Constant",orthoKappa_);
616 if (orthoType_==
"DGKS") {
617 if (orthoKappa_ > 0 && ortho_ != Teuchos::null && !changedOrthoType) {
618 Teuchos::rcp_dynamic_cast<DGKSOrthoManager<ScalarType,MV,OP> >(ortho_)->setDepTol( orthoKappa_ );
624 if (ortho_ == Teuchos::null || changedOrthoType) {
626 Teuchos::RCP<Teuchos::ParameterList> paramsOrtho;
627 if (orthoType_==
"DGKS" && orthoKappa_ > 0) {
628 paramsOrtho->set (
"depTol", orthoKappa_ );
631 ortho_ = factory.
makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, paramsOrtho);
639 if (params->isParameter(
"Convergence Tolerance")) {
640 if (params->isType<MagnitudeType> (
"Convergence Tolerance")) {
641 convtol_ = params->get (
"Convergence Tolerance",
649 params_->set(
"Convergence Tolerance", convtol_);
650 if (convTest_ != Teuchos::null)
654 if (params->isParameter(
"Show Maximum Residual Norm Only")) {
655 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,
"Show Maximum Residual Norm Only");
658 params_->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
659 if (convTest_ != Teuchos::null)
660 convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
664 bool newResTest =
false;
666 std::string tempResScale = resScale_;
667 if (params->isParameter (
"Implicit Residual Scaling")) {
668 tempResScale = params->get<std::string> (
"Implicit Residual Scaling");
672 if (resScale_ != tempResScale) {
675 resScale_ = tempResScale;
678 params_->set (
"Implicit Residual Scaling", resScale_);
680 if (! convTest_.is_null ()) {
683 if (params->isParameter(
"Residual Norm")) {
684 if (params->isType<std::string> (
"Residual Norm")) {
688 convTest_->defineResForm(StatusTestResNorm_t::Implicit, normType);
691 catch (std::exception& e) {
702 if (maxIterTest_ == Teuchos::null)
706 if (convTest_.is_null () || newResTest) {
709 if (params->isParameter(
"Residual Norm")) {
710 if (params->isType<std::string> (
"Residual Norm")) {
715 convTest_ = rcp (
new StatusTestResNorm_t (convtol_, 1, showMaxResNormOnly_));
716 convTest_->defineResForm(StatusTestResNorm_t::Implicit, normType);
720 if (sTest_.is_null () || newResTest)
721 sTest_ = Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
723 if (outputTest_.is_null () || newResTest) {
731 std::string solverDesc =
" Block CG ";
732 outputTest_->setSolverDesc( solverDesc );
737 if (params->isParameter(
"Assert Positive Definiteness")) {
738 assertPositiveDefiniteness_ = Teuchos::getParameter<bool>(*params,
"Assert Positive Definiteness");
739 params_->set(
"Assert Positive Definiteness", assertPositiveDefiniteness_);
743 if (timerSolve_ == Teuchos::null) {
744 std::string solveLabel = label_ +
": BlockCGSolMgr total solve time";
745#ifdef BELOS_TEUCHOS_TIME_MONITOR
746 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
819 using Teuchos::rcp_const_cast;
820 using Teuchos::rcp_dynamic_cast;
827 setParameters(Teuchos::parameterList(*getValidParameters()));
830 Teuchos::LAPACK<int,ScalarType> lapack;
832 TEUCHOS_TEST_FOR_EXCEPTION( !problem_->isProblemSet(),
834 "Belos::BlockCGSolMgr::solve(): Linear problem is not ready, setProblem() "
835 "has not been called.");
839 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
840 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
842 std::vector<int> currIdx, currIdx2;
847 if ( adaptiveBlockSize_ ) {
848 blockSize_ = numCurrRHS;
849 currIdx.resize( numCurrRHS );
850 currIdx2.resize( numCurrRHS );
851 for (
int i=0; i<numCurrRHS; ++i)
852 { currIdx[i] = startPtr+i; currIdx2[i]=i; }
856 currIdx.resize( blockSize_ );
857 currIdx2.resize( blockSize_ );
858 for (
int i=0; i<numCurrRHS; ++i)
859 { currIdx[i] = startPtr+i; currIdx2[i]=i; }
860 for (
int i=numCurrRHS; i<blockSize_; ++i)
861 { currIdx[i] = -1; currIdx2[i] = i; }
865 problem_->setLSIndex( currIdx );
869 Teuchos::ParameterList plist;
870 plist.set(
"Block Size",blockSize_);
873 outputTest_->reset();
877 bool isConverged =
true;
882 plist.set(
"Assert Positive Definiteness", assertPositiveDefiniteness_);
884 RCP<CGIteration<ScalarType,MV,OP> > block_cg_iter;
885 if (blockSize_ == 1) {
889 plist.set(
"Fold Convergence Detection Into Allreduce",
890 foldConvergenceDetectionIntoAllreduce_);
891 if (useSingleReduction_) {
894 outputTest_, convTest_, plist));
899 outputTest_, convTest_, plist));
910#ifdef BELOS_TEUCHOS_TIME_MONITOR
911 Teuchos::TimeMonitor slvtimer(*timerSolve_);
914 while ( numRHS2Solve > 0 ) {
917 std::vector<int> convRHSIdx;
918 std::vector<int> currRHSIdx( currIdx );
919 currRHSIdx.resize(numCurrRHS);
922 block_cg_iter->resetNumIters();
925 outputTest_->resetNumCalls();
928 RCP<MV> R_0 = MVT::CloneViewNonConst( *(rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
933 block_cg_iter->initializeCG(newstate);
939 block_cg_iter->iterate();
943 if (convTest_->getStatus() ==
Passed) {
948 std::vector<int> convIdx =
949 rcp_dynamic_cast<conv_test_type>(convTest_)->
convIndices();
954 if (convIdx.size() == currRHSIdx.size())
959 problem_->setCurrLS();
964 std::vector<int> unconvIdx(currRHSIdx.size());
965 for (
unsigned int i=0; i<currRHSIdx.size(); ++i) {
967 for (
unsigned int j=0; j<convIdx.size(); ++j) {
968 if (currRHSIdx[i] == convIdx[j]) {
974 currIdx2[have] = currIdx2[i];
975 currRHSIdx[have++] = currRHSIdx[i];
980 currRHSIdx.resize(have);
981 currIdx2.resize(have);
984 problem_->setLSIndex( currRHSIdx );
987 std::vector<MagnitudeType> norms;
988 R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
989 for (
int i=0; i<have; ++i) { currIdx2[i] = i; }
992 block_cg_iter->setBlockSize( have );
997 block_cg_iter->initializeCG(defstate);
1003 else if (maxIterTest_->getStatus() ==
Passed) {
1004 isConverged =
false;
1012 TEUCHOS_TEST_FOR_EXCEPTION(
true,std::logic_error,
1013 "Belos::BlockCGSolMgr::solve(): Neither the convergence test nor "
1014 "the maximum iteration count test passed. Please report this bug "
1015 "to the Belos developers.");
1018 catch (
const std::exception &e) {
1019 std::ostream& err = printer_->stream (
Errors);
1020 err <<
"Error! Caught std::exception in CGIteration::iterate() at "
1021 <<
"iteration " << block_cg_iter->getNumIters() << std::endl
1022 << e.what() << std::endl;
1029 problem_->setCurrLS();
1032 startPtr += numCurrRHS;
1033 numRHS2Solve -= numCurrRHS;
1034 if ( numRHS2Solve > 0 ) {
1035 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1037 if ( adaptiveBlockSize_ ) {
1038 blockSize_ = numCurrRHS;
1039 currIdx.resize( numCurrRHS );
1040 currIdx2.resize( numCurrRHS );
1041 for (
int i=0; i<numCurrRHS; ++i)
1042 { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1045 currIdx.resize( blockSize_ );
1046 currIdx2.resize( blockSize_ );
1047 for (
int i=0; i<numCurrRHS; ++i)
1048 { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1049 for (
int i=numCurrRHS; i<blockSize_; ++i)
1050 { currIdx[i] = -1; currIdx2[i] = i; }
1053 problem_->setLSIndex( currIdx );
1056 block_cg_iter->setBlockSize( blockSize_ );
1059 currIdx.resize( numRHS2Solve );
1070#ifdef BELOS_TEUCHOS_TIME_MONITOR
1076 Teuchos::TimeMonitor::summarize( printer_->stream(
TimingDetails) );
1081 numIters_ = maxIterTest_->getNumIters();
1087 const std::vector<MagnitudeType>* pTestValues =
1088 rcp_dynamic_cast<conv_test_type>(convTest_)->
getTestValue();
1090 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1091 "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
1092 "method returned NULL. Please report this bug to the Belos developers.");
1094 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1095 "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() "
1096 "method returned a vector of length zero. Please report this bug to the "
1097 "Belos developers.");
1102 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());