42 #ifndef BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP 43 #define BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP 65 #ifdef BELOS_TEUCHOS_TIME_MONITOR 113 template<
class ScalarType,
class MV,
class OP,
114 const bool supportsScalarType =
118 Belos::Details::LapackSupportsScalar<ScalarType>::value>
140 template<
class ScalarType,
class MV,
class OP>
211 return Teuchos::tuple(timerSolve_);
298 std::string description()
const override;
305 ScalarType & lambda_min,
306 ScalarType & lambda_max,
307 ScalarType & ConditionNumber );
333 static constexpr
int maxIters_default_ = 1000;
334 static constexpr
bool assertPositiveDefiniteness_default_ =
true;
335 static constexpr
bool showMaxResNormOnly_default_ =
false;
338 static constexpr
int outputFreq_default_ = -1;
339 static constexpr
int defQuorum_default_ = 1;
340 static constexpr
const char * resScale_default_ =
"Norm of Initial Residual";
341 static constexpr
const char * label_default_ =
"Belos";
342 static constexpr std::ostream * outputStream_default_ = &std::cout;
343 static constexpr
bool genCondEst_default_ =
false;
364 template<
class ScalarType,
class MV,
class OP>
366 outputStream_(
Teuchos::
rcp(outputStream_default_,false)),
368 maxIters_(maxIters_default_),
370 verbosity_(verbosity_default_),
371 outputStyle_(outputStyle_default_),
372 outputFreq_(outputFreq_default_),
373 defQuorum_(defQuorum_default_),
374 assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
375 showMaxResNormOnly_(showMaxResNormOnly_default_),
376 resScale_(resScale_default_),
377 genCondEst_(genCondEst_default_),
378 condEstimate_(-
Teuchos::ScalarTraits<ScalarType>::one()),
379 label_(label_default_),
384 template<
class ScalarType,
class MV,
class OP>
389 outputStream_(
Teuchos::
rcp(outputStream_default_,false)),
391 maxIters_(maxIters_default_),
393 verbosity_(verbosity_default_),
394 outputStyle_(outputStyle_default_),
395 outputFreq_(outputFreq_default_),
396 defQuorum_(defQuorum_default_),
397 assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
398 showMaxResNormOnly_(showMaxResNormOnly_default_),
399 resScale_(resScale_default_),
400 genCondEst_(genCondEst_default_),
401 condEstimate_(-
Teuchos::ScalarTraits<ScalarType>::one()),
402 label_(label_default_),
406 problem_.is_null (), std::invalid_argument,
407 "Belos::PseudoBlockCGSolMgr two-argument constructor: " 408 "'problem' is null. You must supply a non-null Belos::LinearProblem " 409 "instance when calling this constructor.");
417 template<
class ScalarType,
class MV,
class OP>
423 using Teuchos::parameterList;
427 RCP<const ParameterList> defaultParams = this->getValidParameters ();
446 if (params_.is_null ()) {
448 params_ = parameterList (defaultParams->name ());
455 maxIters_ = params->
get (
"Maximum Iterations", maxIters_default_);
458 params_->set (
"Maximum Iterations", maxIters_);
459 if (! maxIterTest_.is_null ()) {
460 maxIterTest_->setMaxIters (maxIters_);
465 if (params->
isParameter (
"Assert Positive Definiteness")) {
466 assertPositiveDefiniteness_ =
467 params->
get (
"Assert Positive Definiteness",
468 assertPositiveDefiniteness_default_);
471 params_->set (
"Assert Positive Definiteness", assertPositiveDefiniteness_);
476 const std::string tempLabel = params->
get (
"Timer Label", label_default_);
479 if (tempLabel != label_) {
481 params_->set (
"Timer Label", label_);
482 const std::string solveLabel =
483 label_ +
": PseudoBlockCGSolMgr total solve time";
484 #ifdef BELOS_TEUCHOS_TIME_MONITOR 492 if (Teuchos::isParameterType<int> (*params,
"Verbosity")) {
493 verbosity_ = params->
get (
"Verbosity", verbosity_default_);
495 verbosity_ = (int) Teuchos::getParameter<Belos::MsgType> (*params,
"Verbosity");
499 params_->set (
"Verbosity", verbosity_);
500 if (! printer_.is_null ()) {
501 printer_->setVerbosity (verbosity_);
507 if (Teuchos::isParameterType<int> (*params,
"Output Style")) {
508 outputStyle_ = params->
get (
"Output Style", outputStyle_default_);
511 outputStyle_ = (int) Teuchos::getParameter<Belos::OutputType> (*params,
"Output Style");
516 params_->set (
"Output Style", outputStyle_);
517 outputTest_ = Teuchos::null;
522 outputStream_ = params->
get<RCP<std::ostream> > (
"Output Stream");
525 params_->set (
"Output Stream", outputStream_);
526 if (! printer_.is_null ()) {
527 printer_->setOStream (outputStream_);
534 outputFreq_ = params->
get (
"Output Frequency", outputFreq_default_);
538 params_->set (
"Output Frequency", outputFreq_);
539 if (! outputTest_.is_null ()) {
540 outputTest_->setOutputFrequency (outputFreq_);
545 if (params->
isParameter (
"Estimate Condition Number")) {
546 genCondEst_ = params->
get (
"Estimate Condition Number", genCondEst_default_);
550 if (printer_.is_null ()) {
559 if (params->
isParameter (
"Convergence Tolerance")) {
561 convtol_ = params->
get (
"Convergence Tolerance",
569 params_->set (
"Convergence Tolerance", convtol_);
570 if (! convTest_.is_null ()) {
571 convTest_->setTolerance (convtol_);
575 if (params->
isParameter (
"Show Maximum Residual Norm Only")) {
576 showMaxResNormOnly_ = params->
get<
bool> (
"Show Maximum Residual Norm Only");
579 params_->set (
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
580 if (! convTest_.is_null ()) {
581 convTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
586 bool newResTest =
false;
591 std::string tempResScale = resScale_;
592 bool implicitResidualScalingName =
false;
594 tempResScale = params->
get<std::string> (
"Residual Scaling");
596 else if (params->
isParameter (
"Implicit Residual Scaling")) {
597 tempResScale = params->
get<std::string> (
"Implicit Residual Scaling");
598 implicitResidualScalingName =
true;
602 if (resScale_ != tempResScale) {
605 resScale_ = tempResScale;
609 if (implicitResidualScalingName) {
610 params_->set (
"Implicit Residual Scaling", resScale_);
613 params_->set (
"Residual Scaling", resScale_);
616 if (! convTest_.is_null ()) {
620 catch (std::exception& e) {
630 defQuorum_ = params->
get (
"Deflation Quorum", defQuorum_);
631 params_->set (
"Deflation Quorum", defQuorum_);
632 if (! convTest_.is_null ()) {
633 convTest_->setQuorum( defQuorum_ );
640 if (maxIterTest_.is_null ()) {
645 if (convTest_.is_null () || newResTest) {
646 convTest_ =
rcp (
new StatusTestResNorm_t (convtol_, defQuorum_, showMaxResNormOnly_));
650 if (sTest_.is_null () || newResTest) {
651 sTest_ =
rcp (
new StatusTestCombo_t (StatusTestCombo_t::OR, maxIterTest_, convTest_));
654 if (outputTest_.is_null () || newResTest) {
658 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
662 const std::string solverDesc =
" Pseudo Block CG ";
663 outputTest_->setSolverDesc (solverDesc);
667 if (timerSolve_.is_null ()) {
668 const std::string solveLabel =
669 label_ +
": PseudoBlockCGSolMgr total solve time";
670 #ifdef BELOS_TEUCHOS_TIME_MONITOR 680 template<
class ScalarType,
class MV,
class OP>
685 using Teuchos::parameterList;
688 if (validParams_.is_null()) {
690 RCP<ParameterList> pl = parameterList ();
692 "The relative residual tolerance that needs to be achieved by the\n" 693 "iterative solver in order for the linear system to be declared converged.");
694 pl->set(
"Maximum Iterations", static_cast<int>(maxIters_default_),
695 "The maximum number of block iterations allowed for each\n" 696 "set of RHS solved.");
697 pl->set(
"Assert Positive Definiteness", static_cast<bool>(assertPositiveDefiniteness_default_),
698 "Whether or not to assert that the linear operator\n" 699 "and the preconditioner are indeed positive definite.");
700 pl->set(
"Verbosity", static_cast<int>(verbosity_default_),
701 "What type(s) of solver information should be outputted\n" 702 "to the output stream.");
703 pl->set(
"Output Style", static_cast<int>(outputStyle_default_),
704 "What style is used for the solver information outputted\n" 705 "to the output stream.");
706 pl->set(
"Output Frequency", static_cast<int>(outputFreq_default_),
707 "How often convergence information should be outputted\n" 708 "to the output stream.");
709 pl->set(
"Deflation Quorum", static_cast<int>(defQuorum_default_),
710 "The number of linear systems that need to converge before\n" 711 "they are deflated. This number should be <= block size.");
712 pl->set(
"Output Stream",
Teuchos::rcp(outputStream_default_,
false),
713 "A reference-counted pointer to the output stream where all\n" 714 "solver output is sent.");
715 pl->set(
"Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
716 "When convergence information is printed, only show the maximum\n" 717 "relative residual norm when the block size is greater than one.");
718 pl->set(
"Implicit Residual Scaling", resScale_default_,
719 "The type of scaling used in the residual convergence test.");
720 pl->set(
"Estimate Condition Number", static_cast<bool>(genCondEst_default_),
721 "Whether or not to estimate the condition number of the preconditioned system.");
727 pl->set(
"Residual Scaling", resScale_default_,
728 "The type of scaling used in the residual convergence test. This " 729 "name is deprecated; the new name is \"Implicit Residual Scaling\".");
730 pl->set(
"Timer Label", static_cast<const char *>(label_default_),
731 "The string to use as a prefix for the timer labels.");
739 template<
class ScalarType,
class MV,
class OP>
742 const char prefix[] =
"Belos::PseudoBlockCGSolMgr::solve: ";
747 if (!isSet_) { setParameters( params_ ); }
753 prefix <<
"The linear problem to solve is not ready. You must call " 754 "setProblem() on the Belos::LinearProblem instance before telling the " 755 "Belos solver to solve it.");
759 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
760 int numCurrRHS = numRHS2Solve;
762 std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
763 for (
int i=0; i<numRHS2Solve; ++i) {
764 currIdx[i] = startPtr+i;
769 problem_->setLSIndex( currIdx );
775 plist.
set(
"Assert Positive Definiteness",assertPositiveDefiniteness_);
776 if(genCondEst_) plist.
set(
"Max Size For Condest",maxIters_);
779 outputTest_->reset();
782 bool isConverged =
true;
787 if (numRHS2Solve == 1) {
798 #ifdef BELOS_TEUCHOS_TIME_MONITOR 802 bool first_time=
true;
803 while ( numRHS2Solve > 0 ) {
804 if(genCondEst_ && first_time) block_cg_iter->setDoCondEst(
true);
805 else block_cg_iter->setDoCondEst(
false);
808 std::vector<int> convRHSIdx;
809 std::vector<int> currRHSIdx( currIdx );
810 currRHSIdx.resize(numCurrRHS);
813 block_cg_iter->resetNumIters();
816 outputTest_->resetNumCalls();
819 Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
824 block_cg_iter->initializeCG(newState);
831 block_cg_iter->iterate();
838 if ( convTest_->getStatus() ==
Passed ) {
845 if (convIdx.size() == currRHSIdx.size())
849 problem_->setCurrLS();
853 std::vector<int> unconvIdx(currRHSIdx.size());
854 for (
unsigned int i=0; i<currRHSIdx.size(); ++i) {
856 for (
unsigned int j=0; j<convIdx.size(); ++j) {
857 if (currRHSIdx[i] == convIdx[j]) {
863 currIdx2[have] = currIdx2[i];
864 currRHSIdx[have++] = currRHSIdx[i];
867 currRHSIdx.resize(have);
868 currIdx2.resize(have);
871 problem_->setLSIndex( currRHSIdx );
874 std::vector<MagnitudeType> norms;
875 R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
876 for (
int i=0; i<have; ++i) { currIdx2[i] = i; }
881 block_cg_iter->initializeCG(defstate);
889 else if ( maxIterTest_->getStatus() ==
Passed ) {
904 "Belos::PseudoBlockCGSolMgr::solve(): Invalid return from PseudoBlockCGIter::iterate().");
907 catch (
const std::exception &e) {
908 printer_->stream(
Errors) <<
"Error! Caught std::exception in PseudoBlockCGIter::iterate() at iteration " 909 << block_cg_iter->getNumIters() << std::endl
910 << e.what() << std::endl;
916 problem_->setCurrLS();
919 startPtr += numCurrRHS;
920 numRHS2Solve -= numCurrRHS;
922 if ( numRHS2Solve > 0 ) {
924 numCurrRHS = numRHS2Solve;
925 currIdx.resize( numCurrRHS );
926 currIdx2.resize( numCurrRHS );
927 for (
int i=0; i<numCurrRHS; ++i)
928 { currIdx[i] = startPtr+i; currIdx2[i] = i; }
931 problem_->setLSIndex( currIdx );
934 currIdx.resize( numRHS2Solve );
946 #ifdef BELOS_TEUCHOS_TIME_MONITOR 955 numIters_ = maxIterTest_->getNumIters();
959 const std::vector<MagnitudeType>* pTestValues = convTest_->getTestValue();
960 if (pTestValues != NULL && pTestValues->size () > 0) {
961 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
966 ScalarType l_min, l_max;
969 compute_condnum_tridiag_sym(diag,offdiag,l_min,l_max,condEstimate_);
979 template<
class ScalarType,
class MV,
class OP>
982 std::ostringstream oss;
990 template<
class ScalarType,
class MV,
class OP>
995 ScalarType & lambda_min,
996 ScalarType & lambda_max,
997 ScalarType & ConditionNumber )
1008 ScalarType scalar_dummy;
1012 const int N = diag.
size ();
1014 lambda_min = STS::one ();
1015 lambda_max = STS::one ();
1018 &scalar_dummy, 1, &mag_dummy, &info);
1020 (info < 0, std::logic_error,
"Belos::PseudoBlockCGSolMgr::" 1021 "compute_condnum_tridiag_sym: LAPACK's _STEQR failed with info = " 1022 << info <<
" < 0. This suggests there might be a bug in the way Belos " 1023 "is calling LAPACK. Please report this to the Belos developers.");
1024 lambda_min = Teuchos::as<ScalarType> (diag[0]);
1025 lambda_max = Teuchos::as<ScalarType> (diag[N-1]);
1032 if (STS::real (lambda_max) < STS::real (lambda_min)) {
1033 ConditionNumber = STS::one ();
1037 ConditionNumber = lambda_max / lambda_min;
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Teuchos::RCP< const MV > R
The current residual.
Collection of types and exceptions used within the Belos solvers.
ScalarType getConditionEstimate() const
Gets the estimated condition number.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
virtual ~PseudoBlockCGSolMgr()
Destructor.
Class which manages the output and verbosity of the Belos solvers.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
ScaleType
The type of scaling to use on the residual norm value.
PseudoBlockCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
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)
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > getResidualStatusTest() const
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void STEQR(const char &COMPZ, const OrdinalType &n, ScalarType *D, ScalarType *E, ScalarType *Z, const OrdinalType &ldz, ScalarType *WORK, OrdinalType *info) const
A factory class for generating StatusTestOutput objects.
PseudoBlockCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
An implementation of StatusTestResNorm using a family of residual norms.
Teuchos::RCP< Teuchos::ParameterList > params_
Belos::StatusTest class for specifying a maximum number of iterations.
Details::SolverManagerRequiresLapack< ScalarType, MV, OP, scalarTypeIsSupported > base_type
Teuchos::RCP< std::ostream > outputStream_
static std::string name()
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
A factory class for generating StatusTestOutput objects.
Teuchos::ScalarTraits< ScalarType > SCT
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
Traits class which defines basic operations on multivectors.
The Belos::PseudoBlockCGSolMgr provides a powerful and fully-featured solver manager over the pseudo-...
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
This class implements the pseudo-block CG iteration, where the basic CG algorithm is performed on all...
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > convTest_
A Belos::StatusTest class for specifying a maximum number of iterations.
ResetType
How to reset the solver.
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.
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
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)
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Teuchos::RCP< const Teuchos::ParameterList > validParams_
List of valid parameters and their default values.
A linear system to solve, and its associated information.
Class which describes the linear problem to be solved by the iterative solver.
PseudoBlockCGSolMgrOrthoFailure(const std::string &what_arg)
bool isType(const std::string &name) const
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve...
Teuchos::RCP< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
ReturnType
Whether the Belos solve converged for all linear systems.
Teuchos::ScalarTraits< MagnitudeType > MT
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > getResidualStatusTest() const
Return the residual status test.
Teuchos::RCP< Teuchos::Time > timerSolve_
virtual ~PseudoBlockCGSolMgr()
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.
PseudoBlockCGSolMgrLinearProblemFailure(const std::string &what_arg)
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Teuchos::RCP< OutputManager< ScalarType > > printer_
A class for extending the status testing capabilities of Belos via logical combinations.
bool isParameter(const std::string &name) const
PseudoBlockCGSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate or...
Class which defines basic traits for the operator type.
Belos concrete class for performing the pseudo-block CG iteration.
Parent class to all Belos exceptions.
Default parameters common to most Belos solvers.
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...
MultiVecTraits< ScalarType, MV > MVT
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
static const bool scalarTypeIsSupported
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.