42 #ifndef BELOS_BLOCK_CG_SOLMGR_HPP 43 #define BELOS_BLOCK_CG_SOLMGR_HPP 68 #ifdef BELOS_TEUCHOS_TIME_MONITOR 71 #if defined(HAVE_TEUCHOSCORE_CXX11) 72 # include <type_traits> 73 #endif // defined(HAVE_TEUCHOSCORE_CXX11) 118 template<
class ScalarType,
class MV,
class OP,
119 const bool lapackSupportsScalarType =
144 template<
class ScalarType,
class MV,
class OP>
256 return Teuchos::tuple(timerSolve_);
327 std::string description()
const override;
365 static constexpr
int maxIters_default_ = 1000;
366 static constexpr
bool adaptiveBlockSize_default_ =
true;
367 static constexpr
bool showMaxResNormOnly_default_ =
false;
368 static constexpr
bool useSingleReduction_default_ =
false;
369 static constexpr
int blockSize_default_ = 1;
372 static constexpr
int outputFreq_default_ = -1;
373 static constexpr
const char * resScale_default_ =
"Norm of Initial Residual";
374 static constexpr
const char * label_default_ =
"Belos";
375 static constexpr
const char * orthoType_default_ =
"DGKS";
376 static constexpr std::ostream * outputStream_default_ = &std::cout;
418 template<
class ScalarType,
class MV,
class OP>
420 outputStream_(
Teuchos::
rcp(outputStream_default_,false)),
424 maxIters_(maxIters_default_),
426 blockSize_(blockSize_default_),
427 verbosity_(verbosity_default_),
428 outputStyle_(outputStyle_default_),
429 outputFreq_(outputFreq_default_),
430 adaptiveBlockSize_(adaptiveBlockSize_default_),
431 showMaxResNormOnly_(showMaxResNormOnly_default_),
432 useSingleReduction_(useSingleReduction_default_),
433 orthoType_(orthoType_default_),
434 resScale_(resScale_default_),
435 label_(label_default_),
441 template<
class ScalarType,
class MV,
class OP>
446 outputStream_(
Teuchos::
rcp(outputStream_default_,false)),
450 maxIters_(maxIters_default_),
452 blockSize_(blockSize_default_),
453 verbosity_(verbosity_default_),
454 outputStyle_(outputStyle_default_),
455 outputFreq_(outputFreq_default_),
456 adaptiveBlockSize_(adaptiveBlockSize_default_),
457 showMaxResNormOnly_(showMaxResNormOnly_default_),
458 useSingleReduction_(useSingleReduction_default_),
459 orthoType_(orthoType_default_),
460 resScale_(resScale_default_),
461 label_(label_default_),
465 "BlockCGSolMgr's constructor requires a nonnull LinearProblem instance.");
475 template<
class ScalarType,
class MV,
class OP>
481 if (params_ == Teuchos::null) {
490 maxIters_ = params->
get(
"Maximum Iterations",maxIters_default_);
493 params_->set(
"Maximum Iterations", maxIters_);
494 if (maxIterTest_!=Teuchos::null)
495 maxIterTest_->setMaxIters( maxIters_ );
500 blockSize_ = params->
get(
"Block Size",blockSize_default_);
502 "Belos::BlockCGSolMgr: \"Block Size\" must be strictly positive.");
505 params_->set(
"Block Size", blockSize_);
510 adaptiveBlockSize_ = params->
get(
"Adaptive Block Size",adaptiveBlockSize_default_);
513 params_->set(
"Adaptive Block Size", adaptiveBlockSize_);
518 useSingleReduction_ = params->
get(
"Use Single Reduction", useSingleReduction_default_);
523 std::string tempLabel = params->
get(
"Timer Label", label_default_);
526 if (tempLabel != label_) {
528 params_->set(
"Timer Label", label_);
529 std::string solveLabel = label_ +
": BlockCGSolMgr total solve time";
530 #ifdef BELOS_TEUCHOS_TIME_MONITOR 533 if (ortho_ != Teuchos::null) {
534 ortho_->setLabel( label_ );
541 std::string tempOrthoType = params->
get(
"Orthogonalization",orthoType_default_);
543 std::invalid_argument,
544 "Belos::BlockCGSolMgr: \"Orthogonalization\" must be either \"DGKS\", \"ICGS\", or \"IMGS\".");
545 if (tempOrthoType != orthoType_) {
546 orthoType_ = tempOrthoType;
547 params_->set(
"Orthogonalization", orthoType_);
549 if (orthoType_==
"DGKS") {
550 if (orthoKappa_ <= 0) {
558 else if (orthoType_==
"ICGS") {
561 else if (orthoType_==
"IMGS") {
568 if (params->
isParameter(
"Orthogonalization Constant")) {
570 orthoKappa_ = params->
get (
"Orthogonalization Constant",
574 orthoKappa_ = params->
get (
"Orthogonalization Constant",
579 params_->set(
"Orthogonalization Constant",orthoKappa_);
580 if (orthoType_==
"DGKS") {
581 if (orthoKappa_ > 0 && ortho_ != Teuchos::null) {
589 if (Teuchos::isParameterType<int>(*params,
"Verbosity")) {
590 verbosity_ = params->
get(
"Verbosity", verbosity_default_);
592 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,
"Verbosity");
596 params_->set(
"Verbosity", verbosity_);
597 if (printer_ != Teuchos::null)
598 printer_->setVerbosity(verbosity_);
603 if (Teuchos::isParameterType<int>(*params,
"Output Style")) {
604 outputStyle_ = params->
get(
"Output Style", outputStyle_default_);
606 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,
"Output Style");
610 params_->set(
"Output Style", outputStyle_);
611 outputTest_ = Teuchos::null;
616 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,
"Output Stream");
619 params_->set(
"Output Stream", outputStream_);
620 if (printer_ != Teuchos::null)
621 printer_->setOStream( outputStream_ );
627 outputFreq_ = params->
get(
"Output Frequency", outputFreq_default_);
631 params_->set(
"Output Frequency", outputFreq_);
632 if (outputTest_ != Teuchos::null)
633 outputTest_->setOutputFrequency( outputFreq_ );
637 if (printer_ == Teuchos::null) {
646 if (params->
isParameter(
"Convergence Tolerance")) {
648 convtol_ = params->
get (
"Convergence Tolerance",
656 params_->set(
"Convergence Tolerance", convtol_);
657 if (convTest_ != Teuchos::null)
658 convTest_->setTolerance( convtol_ );
661 if (params->
isParameter(
"Show Maximum Residual Norm Only")) {
662 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,
"Show Maximum Residual Norm Only");
665 params_->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
666 if (convTest_ != Teuchos::null)
667 convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
671 bool newResTest =
false;
673 std::string tempResScale = resScale_;
674 if (params->
isParameter (
"Implicit Residual Scaling")) {
675 tempResScale = params->
get<std::string> (
"Implicit Residual Scaling");
679 if (resScale_ != tempResScale) {
682 resScale_ = tempResScale;
685 params_->set (
"Implicit Residual Scaling", resScale_);
687 if (! convTest_.is_null ()) {
691 catch (std::exception& e) {
702 if (maxIterTest_ == Teuchos::null)
706 if (convTest_.is_null () || newResTest) {
707 convTest_ =
rcp (
new StatusTestResNorm_t (convtol_, 1, showMaxResNormOnly_));
711 if (sTest_.is_null () || newResTest)
712 sTest_ =
Teuchos::rcp(
new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
714 if (outputTest_.is_null () || newResTest) {
722 std::string solverDesc =
" Block CG ";
723 outputTest_->setSolverDesc( solverDesc );
728 if (ortho_ == Teuchos::null) {
729 params_->set(
"Orthogonalization", orthoType_);
730 if (orthoType_==
"DGKS") {
731 if (orthoKappa_ <= 0) {
739 else if (orthoType_==
"ICGS") {
742 else if (orthoType_==
"IMGS") {
747 "Belos::BlockCGSolMgr(): Invalid orthogonalization type.");
752 if (timerSolve_ == Teuchos::null) {
753 std::string solveLabel = label_ +
": BlockCGSolMgr total solve time";
754 #ifdef BELOS_TEUCHOS_TIME_MONITOR 764 template<
class ScalarType,
class MV,
class OP>
774 "The relative residual tolerance that needs to be achieved by the\n" 775 "iterative solver in order for the linear system to be declared converged.");
776 pl->
set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
777 "The maximum number of block iterations allowed for each\n" 778 "set of RHS solved.");
779 pl->
set(
"Block Size", static_cast<int>(blockSize_default_),
780 "The number of vectors in each block.");
781 pl->
set(
"Adaptive Block Size", static_cast<bool>(adaptiveBlockSize_default_),
782 "Whether the solver manager should adapt to the block size\n" 783 "based on the number of RHS to solve.");
784 pl->
set(
"Verbosity", static_cast<int>(verbosity_default_),
785 "What type(s) of solver information should be outputted\n" 786 "to the output stream.");
787 pl->
set(
"Output Style", static_cast<int>(outputStyle_default_),
788 "What style is used for the solver information outputted\n" 789 "to the output stream.");
790 pl->
set(
"Output Frequency", static_cast<int>(outputFreq_default_),
791 "How often convergence information should be outputted\n" 792 "to the output stream.");
794 "A reference-counted pointer to the output stream where all\n" 795 "solver output is sent.");
796 pl->
set(
"Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
797 "When convergence information is printed, only show the maximum\n" 798 "relative residual norm when the block size is greater than one.");
799 pl->
set(
"Use Single Reduction", static_cast<bool>(useSingleReduction_default_),
800 "Use single reduction iteration when the block size is one.");
801 pl->
set(
"Implicit Residual Scaling", resScale_default_,
802 "The type of scaling used in the residual convergence test.");
803 pl->
set(
"Timer Label", static_cast<const char *>(label_default_),
804 "The string to use as a prefix for the timer labels.");
805 pl->
set(
"Orthogonalization", static_cast<const char *>(orthoType_default_),
806 "The type of orthogonalization to use: DGKS, ICGS, or IMGS.");
808 "The constant used by DGKS orthogonalization to determine\n" 809 "whether another step of classical Gram-Schmidt is necessary.");
817 template<
class ScalarType,
class MV,
class OP>
821 using Teuchos::rcp_const_cast;
822 using Teuchos::rcp_dynamic_cast;
829 setParameters(Teuchos::parameterList(*getValidParameters()));
837 "Belos::BlockCGSolMgr::solve(): Linear problem is not ready, setProblem() " 838 "has not been called.");
842 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
843 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
845 std::vector<int> currIdx, currIdx2;
850 if ( adaptiveBlockSize_ ) {
851 blockSize_ = numCurrRHS;
852 currIdx.resize( numCurrRHS );
853 currIdx2.resize( numCurrRHS );
854 for (
int i=0; i<numCurrRHS; ++i)
855 { currIdx[i] = startPtr+i; currIdx2[i]=i; }
859 currIdx.resize( blockSize_ );
860 currIdx2.resize( blockSize_ );
861 for (
int i=0; i<numCurrRHS; ++i)
862 { currIdx[i] = startPtr+i; currIdx2[i]=i; }
863 for (
int i=numCurrRHS; i<blockSize_; ++i)
864 { currIdx[i] = -1; currIdx2[i] = i; }
868 problem_->setLSIndex( currIdx );
873 plist.
set(
"Block Size",blockSize_);
876 outputTest_->reset();
880 bool isConverged =
true;
885 RCP<CGIteration<ScalarType,MV,OP> > block_cg_iter;
886 if (blockSize_ == 1) {
890 if (useSingleReduction_) {
893 outputTest_, plist));
898 outputTest_, plist));
909 #ifdef BELOS_TEUCHOS_TIME_MONITOR 913 while ( numRHS2Solve > 0 ) {
916 std::vector<int> convRHSIdx;
917 std::vector<int> currRHSIdx( currIdx );
918 currRHSIdx.resize(numCurrRHS);
921 block_cg_iter->resetNumIters();
924 outputTest_->resetNumCalls();
927 RCP<MV> R_0 = MVT::CloneViewNonConst( *(rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
932 block_cg_iter->initializeCG(newstate);
938 block_cg_iter->iterate();
942 if (convTest_->getStatus() ==
Passed) {
947 std::vector<int> convIdx =
948 rcp_dynamic_cast<conv_test_type>(convTest_)->convIndices();
953 if (convIdx.size() == currRHSIdx.size())
958 problem_->setCurrLS();
963 std::vector<int> unconvIdx(currRHSIdx.size());
964 for (
unsigned int i=0; i<currRHSIdx.size(); ++i) {
966 for (
unsigned int j=0; j<convIdx.size(); ++j) {
967 if (currRHSIdx[i] == convIdx[j]) {
973 currIdx2[have] = currIdx2[i];
974 currRHSIdx[have++] = currRHSIdx[i];
979 currRHSIdx.resize(have);
980 currIdx2.resize(have);
983 problem_->setLSIndex( currRHSIdx );
986 std::vector<MagnitudeType> norms;
987 R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
988 for (
int i=0; i<have; ++i) { currIdx2[i] = i; }
991 block_cg_iter->setBlockSize( have );
996 block_cg_iter->initializeCG(defstate);
1002 else if (maxIterTest_->getStatus() ==
Passed) {
1003 isConverged =
false;
1012 "Belos::BlockCGSolMgr::solve(): Neither the convergence test nor " 1013 "the maximum iteration count test passed. Please report this bug " 1014 "to the Belos developers.");
1017 catch (
const std::exception &e) {
1018 std::ostream& err = printer_->stream (
Errors);
1019 err <<
"Error! Caught std::exception in CGIteration::iterate() at " 1020 <<
"iteration " << block_cg_iter->getNumIters() << std::endl
1021 << e.what() << std::endl;
1028 problem_->setCurrLS();
1031 startPtr += numCurrRHS;
1032 numRHS2Solve -= numCurrRHS;
1033 if ( numRHS2Solve > 0 ) {
1034 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1036 if ( adaptiveBlockSize_ ) {
1037 blockSize_ = numCurrRHS;
1038 currIdx.resize( numCurrRHS );
1039 currIdx2.resize( numCurrRHS );
1040 for (
int i=0; i<numCurrRHS; ++i)
1041 { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1044 currIdx.resize( blockSize_ );
1045 currIdx2.resize( blockSize_ );
1046 for (
int i=0; i<numCurrRHS; ++i)
1047 { currIdx[i] = startPtr+i; currIdx2[i] = i; }
1048 for (
int i=numCurrRHS; i<blockSize_; ++i)
1049 { currIdx[i] = -1; currIdx2[i] = i; }
1052 problem_->setLSIndex( currIdx );
1055 block_cg_iter->setBlockSize( blockSize_ );
1058 currIdx.resize( numRHS2Solve );
1069 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1080 numIters_ = maxIterTest_->getNumIters();
1086 const std::vector<MagnitudeType>* pTestValues =
1087 rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1090 "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() " 1091 "method returned NULL. Please report this bug to the Belos developers.");
1094 "Belos::BlockCGSolMgr::solve(): The convergence test's getTestValue() " 1095 "method returned a vector of length zero. Please report this bug to the " 1096 "Belos developers.");
1101 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1111 template<
class ScalarType,
class MV,
class OP>
1114 std::ostringstream oss;
1117 oss <<
"Ortho Type='"<<orthoType_<<
"\', Block Size=" << blockSize_;
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Teuchos::RCP< const MV > R
The current residual.
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
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...
int maxIters_
Maximum iteration count (read from parameter list).
Class which manages the output and verbosity of the Belos solvers.
bool is_null(const boost::shared_ptr< T > &p)
bool isSet_
Whether or not the parameters have been set (via setParameters()).
Teuchos::RCP< Teuchos::ParameterList > params_
Current parameter list.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Output "status test" that controls all the other status tests.
ScaleType
The type of scaling to use on the residual norm value.
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
Teuchos::ScalarTraits< ScalarType > SCT
static constexpr double orthoKappa
DGKS orthogonalization constant.
This class implements the preconditioned Conjugate Gradient (CG) iteration.
T & get(ParameterList &l, const std::string &name)
Structure to contain pointers to CGIteration state variables.
Belos concrete class for performing the conjugate-gradient (CG) iteration.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
A factory class for generating StatusTestOutput objects.
Belos concrete class for performing the block conjugate-gradient (CG) iteration.
An implementation of StatusTestResNorm using a family of residual norms.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
std::string label_
Prefix label for all the timers.
The Belos::BlockCGSolMgr provides a powerful and fully-featured solver manager over the CG and BlockC...
Belos::StatusTest class for specifying a maximum number of iterations.
static std::string name()
BlockCGSolMgrLinearProblemFailure(const std::string &what_arg)
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
BlockCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
A factory class for generating StatusTestOutput objects.
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
Orthogonalization manager.
Iterated Modified Gram-Schmidt (IMGS) implementation of the Belos::OrthoManager class.
BlockCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Maximum iteration count stopping criterion.
Traits class which defines basic operations on multivectors.
Belos::StatusTest for logically combining several status tests.
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
Aggregate stopping criterion.
Details::SolverManagerRequiresLapack< ScalarType, MV, OP, requiresLapack > base_type
virtual ~BlockCGSolMgr()
Destructor.
Classical Gram-Schmidt (with DGKS correction) implementation of the Belos::OrthoManager class...
OperatorTraits< ScalarType, MV, OP > OPT
MultiVecTraits< ScalarType, MV > MVT
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Pure virtual base class which describes the basic interface for a solver manager. ...
static constexpr double convTol
Default convergence tolerance.
Belos concrete class for performing a single-reduction conjugate-gradient (CG) iteration.
static void summarize(Ptr< const Comm< int > > comm, std::ostream &out=std::cout, const bool alwaysWriteLocal=false, const bool writeGlobalStats=true, const bool writeZeroTimers=true, const ECounterSetOp setOp=Intersection, const std::string &filter="", const bool ignoreZeroTimers=false)
bool is_null(const RCP< T > &p)
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
MagnitudeType achievedTol_
Tolerance achieved by the last solve() invocation.
BlockCGSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonor...
bool isType(const std::string &name) const
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
ReturnType
Whether the Belos solve converged for all linear systems.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > convTest_
Convergence stopping criterion.
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
The linear problem to solve.
Iterated Classical Gram-Schmidt (ICGS) implementation of the Belos::OrthoManager class.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
BlockCGSolMgrOrthoFailure(const std::string &what_arg)
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
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.
MagnitudeType convtol_
Convergence tolerance (read from parameter list).
Teuchos::RCP< Teuchos::Time > timerSolve_
Solve timer.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
A class for extending the status testing capabilities of Belos via logical combinations.
bool isParameter(const std::string &name) const
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
Class which defines basic traits for the operator type.
Parent class to all Belos exceptions.
Default parameters common to most Belos solvers.
This class implements the preconditioned single-reduction Conjugate Gradient (CG) iteration...
Teuchos::RCP< OutputManager< ScalarType > > printer_
Output manager, that handles printing of different kinds of messages.
Base class for Belos::SolverManager subclasses which normally can only compile with ScalarType types ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
int numIters_
Number of iterations taken by the last solve() invocation.
Teuchos::ScalarTraits< MagnitudeType > MT
static const bool requiresLapack
Teuchos::RCP< std::ostream > outputStream_
Output stream to which the output manager prints.
MagnitudeType orthoKappa_
Orthogonalization parameter (read from parameter list).
Stub implementation of BlockCGIter, for ScalarType types for which Teuchos::LAPACK does NOT have a va...