42 #ifndef BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP
43 #define BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP
64 #ifdef BELOS_TEUCHOS_TIME_MONITOR
102 template<
class ScalarType,
class MV,
class OP,
103 const bool supportsScalarType =
107 Belos::Details::LapackSupportsScalar<ScalarType>::value>
129 template<
class ScalarType,
class MV,
class OP>
200 return Teuchos::tuple(timerSolve_);
288 std::string description()
const override;
296 ScalarType & lambda_min,
297 ScalarType & lambda_max,
298 ScalarType & ConditionNumber );
324 static constexpr
int maxIters_default_ = 1000;
325 static constexpr
bool assertPositiveDefiniteness_default_ =
true;
326 static constexpr
bool showMaxResNormOnly_default_ =
false;
329 static constexpr
int outputFreq_default_ = -1;
330 static constexpr
int defQuorum_default_ = 1;
331 static constexpr
bool foldConvergenceDetectionIntoAllreduce_default_ =
false;
332 static constexpr
const char * resScale_default_ =
"Norm of Initial Residual";
333 static constexpr
const char * label_default_ =
"Belos";
334 static constexpr
bool genCondEst_default_ =
false;
357 template<
class ScalarType,
class MV,
class OP>
359 outputStream_(Teuchos::rcpFromRef(std::cout)),
361 maxIters_(maxIters_default_),
363 verbosity_(verbosity_default_),
364 outputStyle_(outputStyle_default_),
365 outputFreq_(outputFreq_default_),
366 defQuorum_(defQuorum_default_),
367 assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
368 showMaxResNormOnly_(showMaxResNormOnly_default_),
369 foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
370 resScale_(resScale_default_),
371 genCondEst_(genCondEst_default_),
372 condEstimate_(-Teuchos::ScalarTraits<ScalarType>::one()),
373 label_(label_default_),
378 template<
class ScalarType,
class MV,
class OP>
383 outputStream_(Teuchos::rcpFromRef(std::cout)),
385 maxIters_(maxIters_default_),
387 verbosity_(verbosity_default_),
388 outputStyle_(outputStyle_default_),
389 outputFreq_(outputFreq_default_),
390 defQuorum_(defQuorum_default_),
391 assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
392 showMaxResNormOnly_(showMaxResNormOnly_default_),
393 foldConvergenceDetectionIntoAllreduce_(foldConvergenceDetectionIntoAllreduce_default_),
394 resScale_(resScale_default_),
395 genCondEst_(genCondEst_default_),
396 condEstimate_(-Teuchos::ScalarTraits<ScalarType>::one()),
397 label_(label_default_),
401 problem_.is_null (), std::invalid_argument,
402 "Belos::PseudoBlockCGSolMgr two-argument constructor: "
403 "'problem' is null. You must supply a non-null Belos::LinearProblem "
404 "instance when calling this constructor.");
412 template<
class ScalarType,
class MV,
class OP>
418 using Teuchos::parameterList;
422 RCP<const ParameterList> defaultParams = this->getValidParameters ();
441 if (params_.is_null ()) {
443 params_ = parameterList (defaultParams->name ());
450 maxIters_ = params->
get (
"Maximum Iterations", maxIters_default_);
453 params_->set (
"Maximum Iterations", maxIters_);
454 if (! maxIterTest_.is_null ()) {
455 maxIterTest_->setMaxIters (maxIters_);
460 if (params->
isParameter (
"Assert Positive Definiteness")) {
461 assertPositiveDefiniteness_ =
462 params->
get (
"Assert Positive Definiteness",
463 assertPositiveDefiniteness_default_);
466 params_->set (
"Assert Positive Definiteness", assertPositiveDefiniteness_);
469 if (params->
isParameter(
"Fold Convergence Detection Into Allreduce")) {
470 foldConvergenceDetectionIntoAllreduce_ = params->
get(
"Fold Convergence Detection Into Allreduce",
471 foldConvergenceDetectionIntoAllreduce_default_);
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_);
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::rcpFromRef(std::cout),
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.");
732 pl->set(
"Fold Convergence Detection Into Allreduce",static_cast<bool>(foldConvergenceDetectionIntoAllreduce_default_),
733 "Merge the allreduce for convergence detection with the one for CG.\n"
734 "This saves one all-reduce, but incurs more computation.");
742 template<
class ScalarType,
class MV,
class OP>
745 const char prefix[] =
"Belos::PseudoBlockCGSolMgr::solve: ";
750 if (!isSet_) { setParameters( params_ ); }
754 prefix <<
"The linear problem to solve is not ready. You must call "
755 "setProblem() on the Belos::LinearProblem instance before telling the "
756 "Belos solver to solve it.");
760 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
761 int numCurrRHS = numRHS2Solve;
763 std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
764 for (
int i=0; i<numRHS2Solve; ++i) {
765 currIdx[i] = startPtr+i;
770 problem_->setLSIndex( currIdx );
776 plist.
set(
"Assert Positive Definiteness",assertPositiveDefiniteness_);
777 if(genCondEst_) plist.
set(
"Max Size For Condest",maxIters_);
780 outputTest_->reset();
783 bool isConverged =
true;
788 if (numRHS2Solve == 1) {
789 plist.
set(
"Fold Convergence Detection Into Allreduce",
790 foldConvergenceDetectionIntoAllreduce_);
799 block_cg_iter->setDoCondEst(genCondEst_);
800 bool condEstPerf =
false;
804 #ifdef BELOS_TEUCHOS_TIME_MONITOR
808 while ( numRHS2Solve > 0 ) {
811 std::vector<int> convRHSIdx;
812 std::vector<int> currRHSIdx( currIdx );
813 currRHSIdx.resize(numCurrRHS);
816 block_cg_iter->resetNumIters();
819 outputTest_->resetNumCalls();
822 Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
827 block_cg_iter->initializeCG(newState);
834 block_cg_iter->iterate();
841 if ( convTest_->getStatus() ==
Passed ) {
848 if (convIdx.size() == currRHSIdx.size())
852 problem_->setCurrLS();
856 for (
unsigned int i=0; i<currRHSIdx.size(); ++i) {
858 for (
unsigned int j=0; j<convIdx.size(); ++j) {
859 if (currRHSIdx[i] == convIdx[j]) {
865 currIdx2[have] = currIdx2[i];
866 currRHSIdx[have++] = currRHSIdx[i];
869 currRHSIdx.resize(have);
870 currIdx2.resize(have);
873 if (currRHSIdx[0] != 0 && genCondEst_ && !condEstPerf)
876 ScalarType l_min, l_max;
879 compute_condnum_tridiag_sym(diag,offdiag,eigenEstimates_,l_min,l_max,condEstimate_);
882 block_cg_iter->setDoCondEst(
false);
887 problem_->setLSIndex( currRHSIdx );
890 std::vector<MagnitudeType> norms;
891 R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
892 for (
int i=0; i<have; ++i) { currIdx2[i] = i; }
897 block_cg_iter->initializeCG(defstate);
905 else if ( maxIterTest_->getStatus() ==
Passed ) {
920 "Belos::PseudoBlockCGSolMgr::solve(): Invalid return from PseudoBlockCGIter::iterate().");
925 achievedTol_ = MT::one();
927 MVT::MvInit( *X, SCT::zero() );
928 printer_->stream(
Warnings) <<
"Belos::PseudoBlockCGSolMgr::solve(): Warning! NaN has been detected!"
932 catch (
const std::exception &e) {
933 printer_->stream(
Errors) <<
"Error! Caught std::exception in PseudoBlockCGIter::iterate() at iteration "
934 << block_cg_iter->getNumIters() << std::endl
935 << e.what() << std::endl;
941 problem_->setCurrLS();
944 startPtr += numCurrRHS;
945 numRHS2Solve -= numCurrRHS;
947 if ( numRHS2Solve > 0 ) {
949 numCurrRHS = numRHS2Solve;
950 currIdx.resize( numCurrRHS );
951 currIdx2.resize( numCurrRHS );
952 for (
int i=0; i<numCurrRHS; ++i)
953 { currIdx[i] = startPtr+i; currIdx2[i] = i; }
956 problem_->setLSIndex( currIdx );
959 currIdx.resize( numRHS2Solve );
970 #ifdef BELOS_TEUCHOS_TIME_MONITOR
979 numIters_ = maxIterTest_->getNumIters();
983 const std::vector<MagnitudeType>* pTestValues = convTest_->getTestValue();
984 if (pTestValues != NULL && pTestValues->size () > 0) {
985 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
989 if (genCondEst_ && !condEstPerf) {
990 ScalarType l_min, l_max;
993 compute_condnum_tridiag_sym(diag,offdiag,eigenEstimates_,l_min,l_max,condEstimate_);
1004 template<
class ScalarType,
class MV,
class OP>
1007 std::ostringstream oss;
1015 template<
class ScalarType,
class MV,
class OP>
1021 ScalarType & lambda_min,
1022 ScalarType & lambda_max,
1023 ScalarType & ConditionNumber )
1034 const int N = diag.
size ();
1035 ScalarType scalar_dummy;
1036 std::vector<MagnitudeType> mag_dummy(4*N);
1041 lambda_min = STS::one ();
1042 lambda_max = STS::one ();
1045 &scalar_dummy, 1, &mag_dummy[0], &info);
1047 (info < 0, std::logic_error,
"Belos::PseudoBlockCGSolMgr::"
1048 "compute_condnum_tridiag_sym: LAPACK's _PTEQR failed with info = "
1049 << info <<
" < 0. This suggests there might be a bug in the way Belos "
1050 "is calling LAPACK. Please report this to the Belos developers.");
1051 for (
int k = 0; k < N; k++) {
1052 lambdas[k] = diag[N - 1 - k];
1054 lambda_min = Teuchos::as<ScalarType> (diag[N-1]);
1055 lambda_max = Teuchos::as<ScalarType> (diag[0]);
1062 if (STS::real (lambda_max) < STS::real (lambda_min)) {
1063 ConditionNumber = STS::one ();
1067 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.
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::ArrayRCP< MagnitudeType > eigenEstimates_
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)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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.
ScalarType getConditionEstimate() const
Gets the estimated condition number.
void PTEQR(const char &COMPZ, const OrdinalType &n, MagnitudeType *D, MagnitudeType *E, ScalarType *Z, const OrdinalType &ldz, MagnitudeType *WORK, OrdinalType *info) const
Teuchos::RCP< Teuchos::ParameterList > params_
static const double convTol
Default convergence tolerance.
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 resize(const size_type n, const T &val=T())
bool isParameter(const std::string &name) 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. ...
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)
bool foldConvergenceDetectionIntoAllreduce_
Teuchos::ArrayRCP< MagnitudeType > getEigenEstimates() const
Teuchos::RCP< const Teuchos::ParameterList > validParams_
List of valid parameters and their default values.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > getResidualStatusTest() const
Return the residual status test.
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.
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
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
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...
bool isType(const std::string &name) const
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.
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.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > getResidualStatusTest() const