42 #ifndef BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP 43 #define BELOS_PSEUDO_BLOCK_CG_SOLMGR_HPP 63 #ifdef BELOS_TEUCHOS_TIME_MONITOR 111 template<
class ScalarType,
class MV,
class OP,
112 const bool supportsScalarType =
116 Belos::Details::LapackSupportsScalar<ScalarType>::value>
118 static const bool scalarTypeIsSupported =
138 template<
class ScalarType,
class MV,
class OP>
204 return Teuchos::tuple(timerSolve_);
291 std::string description()
const;
298 ScalarType & lambda_min,
299 ScalarType & lambda_max,
300 ScalarType & ConditionNumber );
329 static const MagnitudeType convtol_default_;
330 static const int maxIters_default_;
331 static const bool assertPositiveDefiniteness_default_;
332 static const bool showMaxResNormOnly_default_;
333 static const int verbosity_default_;
334 static const int outputStyle_default_;
335 static const int outputFreq_default_;
336 static const int defQuorum_default_;
337 static const std::string resScale_default_;
338 static const std::string label_default_;
340 static const bool genCondEst_default_;
343 MagnitudeType convtol_,achievedTol_;
344 int maxIters_, numIters_;
345 int verbosity_, outputStyle_, outputFreq_, defQuorum_;
346 bool assertPositiveDefiniteness_, showMaxResNormOnly_;
347 std::string resScale_;
349 ScalarType condEstimate_;
361 template<
class ScalarType,
class MV,
class OP>
362 const typename PseudoBlockCGSolMgr<ScalarType,MV,OP,true>::MagnitudeType PseudoBlockCGSolMgr<ScalarType,MV,OP,true>::convtol_default_ = 1e-8;
364 template<
class ScalarType,
class MV,
class OP>
365 const int PseudoBlockCGSolMgr<ScalarType,MV,OP,true>::maxIters_default_ = 1000;
367 template<
class ScalarType,
class MV,
class OP>
368 const bool PseudoBlockCGSolMgr<ScalarType,MV,OP,true>::assertPositiveDefiniteness_default_ =
true;
370 template<
class ScalarType,
class MV,
class OP>
371 const bool PseudoBlockCGSolMgr<ScalarType,MV,OP,true>::showMaxResNormOnly_default_ =
false;
373 template<
class ScalarType,
class MV,
class OP>
374 const int PseudoBlockCGSolMgr<ScalarType,MV,OP,true>::verbosity_default_ =
Belos::Errors;
376 template<
class ScalarType,
class MV,
class OP>
377 const int PseudoBlockCGSolMgr<ScalarType,MV,OP,true>::outputStyle_default_ =
Belos::General;
379 template<
class ScalarType,
class MV,
class OP>
380 const int PseudoBlockCGSolMgr<ScalarType,MV,OP,true>::outputFreq_default_ = -1;
382 template<
class ScalarType,
class MV,
class OP>
383 const int PseudoBlockCGSolMgr<ScalarType,MV,OP,true>::defQuorum_default_ = 1;
385 template<
class ScalarType,
class MV,
class OP>
386 const std::string PseudoBlockCGSolMgr<ScalarType,MV,OP,true>::resScale_default_ =
"Norm of Initial Residual";
388 template<
class ScalarType,
class MV,
class OP>
389 const std::string PseudoBlockCGSolMgr<ScalarType,MV,OP,true>::label_default_ =
"Belos";
391 template<
class ScalarType,
class MV,
class OP>
394 template<
class ScalarType,
class MV,
class OP>
395 const bool PseudoBlockCGSolMgr<ScalarType,MV,OP,true>::genCondEst_default_ =
false;
398 template<
class ScalarType,
class MV,
class OP>
400 outputStream_(outputStream_default_),
401 convtol_(convtol_default_),
402 maxIters_(maxIters_default_),
404 verbosity_(verbosity_default_),
405 outputStyle_(outputStyle_default_),
406 outputFreq_(outputFreq_default_),
407 defQuorum_(defQuorum_default_),
408 assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
409 showMaxResNormOnly_(showMaxResNormOnly_default_),
410 resScale_(resScale_default_),
411 genCondEst_(genCondEst_default_),
412 condEstimate_(-
Teuchos::ScalarTraits<ScalarType>::one()),
413 label_(label_default_),
418 template<
class ScalarType,
class MV,
class OP>
423 outputStream_(outputStream_default_),
424 convtol_(convtol_default_),
425 maxIters_(maxIters_default_),
427 verbosity_(verbosity_default_),
428 outputStyle_(outputStyle_default_),
429 outputFreq_(outputFreq_default_),
430 defQuorum_(defQuorum_default_),
431 assertPositiveDefiniteness_(assertPositiveDefiniteness_default_),
432 showMaxResNormOnly_(showMaxResNormOnly_default_),
433 resScale_(resScale_default_),
434 genCondEst_(genCondEst_default_),
435 condEstimate_(-
Teuchos::ScalarTraits<ScalarType>::one()),
436 label_(label_default_),
440 problem_.is_null (), std::invalid_argument,
441 "Belos::PseudoBlockCGSolMgr two-argument constructor: " 442 "'problem' is null. You must supply a non-null Belos::LinearProblem " 443 "instance when calling this constructor.");
451 template<
class ScalarType,
class MV,
class OP>
457 using Teuchos::parameterList;
461 RCP<const ParameterList> defaultParams = this->getValidParameters ();
480 if (params_.is_null ()) {
482 params_ = parameterList (defaultParams->name ());
489 maxIters_ = params->
get (
"Maximum Iterations", maxIters_default_);
492 params_->set (
"Maximum Iterations", maxIters_);
493 if (! maxIterTest_.is_null ()) {
494 maxIterTest_->setMaxIters (maxIters_);
499 if (params->
isParameter (
"Assert Positive Definiteness")) {
500 assertPositiveDefiniteness_ =
501 params->
get (
"Assert Positive Definiteness",
502 assertPositiveDefiniteness_default_);
505 params_->set (
"Assert Positive Definiteness", assertPositiveDefiniteness_);
510 const std::string tempLabel = params->
get (
"Timer Label", label_default_);
513 if (tempLabel != label_) {
515 params_->set (
"Timer Label", label_);
516 const std::string solveLabel =
517 label_ +
": PseudoBlockCGSolMgr total solve time";
518 #ifdef BELOS_TEUCHOS_TIME_MONITOR 521 if (! ortho_.is_null ()) {
522 ortho_->setLabel (label_);
529 if (Teuchos::isParameterType<int> (*params,
"Verbosity")) {
530 verbosity_ = params->
get (
"Verbosity", verbosity_default_);
532 verbosity_ = (int) Teuchos::getParameter<Belos::MsgType> (*params,
"Verbosity");
536 params_->set (
"Verbosity", verbosity_);
537 if (! printer_.is_null ()) {
538 printer_->setVerbosity (verbosity_);
544 if (Teuchos::isParameterType<int> (*params,
"Output Style")) {
545 outputStyle_ = params->
get (
"Output Style", outputStyle_default_);
548 outputStyle_ = (int) Teuchos::getParameter<Belos::OutputType> (*params,
"Output Style");
553 params_->set (
"Output Style", outputStyle_);
554 outputTest_ = Teuchos::null;
559 outputStream_ = params->
get<RCP<std::ostream> > (
"Output Stream");
562 params_->set (
"Output Stream", outputStream_);
563 if (! printer_.is_null ()) {
564 printer_->setOStream (outputStream_);
571 outputFreq_ = params->
get (
"Output Frequency", outputFreq_default_);
575 params_->set (
"Output Frequency", outputFreq_);
576 if (! outputTest_.is_null ()) {
577 outputTest_->setOutputFrequency (outputFreq_);
582 if (params->
isParameter (
"Estimate Condition Number")) {
583 genCondEst_ = params->
get (
"Estimate Condition Number", genCondEst_default_);
587 if (printer_.is_null ()) {
596 if (params->
isParameter (
"Convergence Tolerance")) {
597 convtol_ = params->
get (
"Convergence Tolerance", convtol_default_);
600 params_->set (
"Convergence Tolerance", convtol_);
601 if (! convTest_.is_null ()) {
602 convTest_->setTolerance (convtol_);
606 if (params->
isParameter (
"Show Maximum Residual Norm Only")) {
607 showMaxResNormOnly_ = params->
get<
bool> (
"Show Maximum Residual Norm Only");
610 params_->set (
"Show Maximum Residual Norm Only", showMaxResNormOnly_);
611 if (! convTest_.is_null ()) {
612 convTest_->setShowMaxResNormOnly (showMaxResNormOnly_);
617 bool newResTest =
false;
622 std::string tempResScale = resScale_;
623 bool implicitResidualScalingName =
false;
625 tempResScale = params->
get<std::string> (
"Residual Scaling");
627 else if (params->
isParameter (
"Implicit Residual Scaling")) {
628 tempResScale = params->
get<std::string> (
"Implicit Residual Scaling");
629 implicitResidualScalingName =
true;
633 if (resScale_ != tempResScale) {
636 resScale_ = tempResScale;
640 if (implicitResidualScalingName) {
641 params_->set (
"Implicit Residual Scaling", resScale_);
644 params_->set (
"Residual Scaling", resScale_);
647 if (! convTest_.is_null ()) {
651 catch (std::exception& e) {
661 defQuorum_ = params->
get (
"Deflation Quorum", defQuorum_);
662 params_->set (
"Deflation Quorum", defQuorum_);
663 if (! convTest_.is_null ()) {
664 convTest_->setQuorum( defQuorum_ );
671 if (maxIterTest_.is_null ()) {
676 if (convTest_.is_null () || newResTest) {
677 convTest_ =
rcp (
new StatusTestResNorm_t (convtol_, defQuorum_, showMaxResNormOnly_));
681 if (sTest_.is_null () || newResTest) {
682 sTest_ =
rcp (
new StatusTestCombo_t (StatusTestCombo_t::OR, maxIterTest_, convTest_));
685 if (outputTest_.is_null () || newResTest) {
689 outputTest_ = stoFactory.
create (printer_, sTest_, outputFreq_,
693 const std::string solverDesc =
" Pseudo Block CG ";
694 outputTest_->setSolverDesc (solverDesc);
698 if (timerSolve_.is_null ()) {
699 const std::string solveLabel =
700 label_ +
": PseudoBlockCGSolMgr total solve time";
701 #ifdef BELOS_TEUCHOS_TIME_MONITOR 711 template<
class ScalarType,
class MV,
class OP>
716 using Teuchos::parameterList;
719 if (validParams_.is_null()) {
721 RCP<ParameterList> pl = parameterList ();
722 pl->set(
"Convergence Tolerance", convtol_default_,
723 "The relative residual tolerance that needs to be achieved by the\n" 724 "iterative solver in order for the linera system to be declared converged.");
725 pl->set(
"Maximum Iterations", maxIters_default_,
726 "The maximum number of block iterations allowed for each\n" 727 "set of RHS solved.");
728 pl->set(
"Assert Positive Definiteness", assertPositiveDefiniteness_default_,
729 "Whether or not to assert that the linear operator\n" 730 "and the preconditioner are indeed positive definite.");
731 pl->set(
"Verbosity", verbosity_default_,
732 "What type(s) of solver information should be outputted\n" 733 "to the output stream.");
734 pl->set(
"Output Style", outputStyle_default_,
735 "What style is used for the solver information outputted\n" 736 "to the output stream.");
737 pl->set(
"Output Frequency", outputFreq_default_,
738 "How often convergence information should be outputted\n" 739 "to the output stream.");
740 pl->set(
"Deflation Quorum", defQuorum_default_,
741 "The number of linear systems that need to converge before\n" 742 "they are deflated. This number should be <= block size.");
743 pl->set(
"Output Stream", outputStream_default_,
744 "A reference-counted pointer to the output stream where all\n" 745 "solver output is sent.");
746 pl->set(
"Show Maximum Residual Norm Only", showMaxResNormOnly_default_,
747 "When convergence information is printed, only show the maximum\n" 748 "relative residual norm when the block size is greater than one.");
749 pl->set(
"Implicit Residual Scaling", resScale_default_,
750 "The type of scaling used in the residual convergence test.");
751 pl->set(
"Estimate Condition Number", genCondEst_default_,
752 "Whether or not to estimate the condition number of the preconditioned system.");
758 pl->set(
"Residual Scaling", resScale_default_,
759 "The type of scaling used in the residual convergence test. This " 760 "name is deprecated; the new name is \"Implicit Residual Scaling\".");
761 pl->set(
"Timer Label", label_default_,
762 "The string to use as a prefix for the timer labels.");
771 template<
class ScalarType,
class MV,
class OP>
774 const char prefix[] =
"Belos::PseudoBlockCGSolMgr::solve: ";
779 if (!isSet_) { setParameters( params_ ); }
785 prefix <<
"The linear problem to solve is not ready. You must call " 786 "setProblem() on the Belos::LinearProblem instance before telling the " 787 "Belos solver to solve it.");
791 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
792 int numCurrRHS = numRHS2Solve;
794 std::vector<int> currIdx( numRHS2Solve ), currIdx2( numRHS2Solve );
795 for (
int i=0; i<numRHS2Solve; ++i) {
796 currIdx[i] = startPtr+i;
801 problem_->setLSIndex( currIdx );
807 plist.
set(
"Assert Positive Definiteness",assertPositiveDefiniteness_);
808 if(genCondEst_) plist.
set(
"Max Size For Condest",maxIters_);
811 outputTest_->reset();
814 bool isConverged =
true;
824 #ifdef BELOS_TEUCHOS_TIME_MONITOR 828 bool first_time=
true;
829 while ( numRHS2Solve > 0 ) {
830 if(genCondEst_ && first_time) block_cg_iter->setDoCondEst(
true);
831 else block_cg_iter->setDoCondEst(
false);
834 std::vector<int> convRHSIdx;
835 std::vector<int> currRHSIdx( currIdx );
836 currRHSIdx.resize(numCurrRHS);
839 block_cg_iter->resetNumIters();
842 outputTest_->resetNumCalls();
845 Teuchos::RCP<MV> R_0 = MVT::CloneViewNonConst( *(Teuchos::rcp_const_cast<MV>(problem_->getInitResVec())), currIdx );
850 block_cg_iter->initializeCG(newState);
857 block_cg_iter->iterate();
864 if ( convTest_->getStatus() ==
Passed ) {
871 if (convIdx.size() == currRHSIdx.size())
875 problem_->setCurrLS();
879 std::vector<int> unconvIdx(currRHSIdx.size());
880 for (
unsigned int i=0; i<currRHSIdx.size(); ++i) {
882 for (
unsigned int j=0; j<convIdx.size(); ++j) {
883 if (currRHSIdx[i] == convIdx[j]) {
889 currIdx2[have] = currIdx2[i];
890 currRHSIdx[have++] = currRHSIdx[i];
893 currRHSIdx.resize(have);
894 currIdx2.resize(have);
897 problem_->setLSIndex( currRHSIdx );
900 std::vector<MagnitudeType> norms;
901 R_0 = MVT::CloneCopy( *(block_cg_iter->getNativeResiduals(&norms)),currIdx2 );
902 for (
int i=0; i<have; ++i) { currIdx2[i] = i; }
907 block_cg_iter->initializeCG(defstate);
915 else if ( maxIterTest_->getStatus() ==
Passed ) {
930 "Belos::PseudoBlockCGSolMgr::solve(): Invalid return from PseudoBlockCGIter::iterate().");
933 catch (
const std::exception &e) {
934 printer_->stream(
Errors) <<
"Error! Caught std::exception in PseudoBlockCGIter::iterate() at iteration " 935 << block_cg_iter->getNumIters() << std::endl
936 << e.what() << std::endl;
942 problem_->setCurrLS();
945 startPtr += numCurrRHS;
946 numRHS2Solve -= numCurrRHS;
948 if ( numRHS2Solve > 0 ) {
950 numCurrRHS = numRHS2Solve;
951 currIdx.resize( numCurrRHS );
952 currIdx2.resize( numCurrRHS );
953 for (
int i=0; i<numCurrRHS; ++i)
954 { currIdx[i] = startPtr+i; currIdx2[i] = i; }
957 problem_->setLSIndex( currIdx );
960 currIdx.resize( numRHS2Solve );
972 #ifdef BELOS_TEUCHOS_TIME_MONITOR 981 numIters_ = maxIterTest_->getNumIters();
985 const std::vector<MagnitudeType>* pTestValues = convTest_->getTestValue();
986 if (pTestValues != NULL && pTestValues->size () > 0) {
987 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
992 ScalarType l_min, l_max;
995 compute_condnum_tridiag_sym(diag,offdiag,l_min,l_max,condEstimate_);
1005 template<
class ScalarType,
class MV,
class OP>
1008 std::ostringstream oss;
1016 template<
class ScalarType,
class MV,
class OP>
1021 ScalarType & lambda_min,
1022 ScalarType & lambda_max,
1023 ScalarType & ConditionNumber )
1034 ScalarType scalar_dummy;
1035 MagnitudeType mag_dummy;
1038 const int N = diag.
size ();
1040 lambda_min = STS::one ();
1041 lambda_max = STS::one ();
1044 &scalar_dummy, 1, &mag_dummy, &info);
1046 (info < 0, std::logic_error,
"Belos::PseudoBlockCGSolMgr::" 1047 "compute_condnum_tridiag_sym: LAPACK's _STEQR failed with info = " 1048 << info <<
" < 0. This suggests there might be a bug in the way Belos " 1049 "is calling LAPACK. Please report this to the Belos developers.");
1050 lambda_min = Teuchos::as<ScalarType> (diag[0]);
1051 lambda_max = Teuchos::as<ScalarType> (diag[N-1]);
1058 if (STS::real (lambda_max) < STS::real (lambda_min)) {
1059 ConditionNumber = STS::one ();
1063 ConditionNumber = lambda_max / lambda_min;
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
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.
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem that needs to be solved.
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)
Structure to contain pointers to CGIteration state variables.
T & get(const std::string &name, T def_value)
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)
A factory class for generating StatusTestOutput objects.
PseudoBlockCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i...
bool isLOADetected() const
Return whether a loss of accuracy was detected by this solver during the most current solve...
An implementation of StatusTestResNorm using a family of residual norms.
Belos::StatusTest class for specifying a maximum number of iterations.
static std::string name()
A factory class for generating StatusTestOutput objects.
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
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Return a reference to the linear problem being solved by this solver manager.
This class implements the pseudo-block CG iteration, where the basic CG algorithm is performed on all...
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 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)
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
MagnitudeType achievedTol() const
Tolerance achieved by the last solve() invocation.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
A linear system to solve, and its associated information.
int getNumIters() const
Get the iteration count for the most recent call to solve().
Class which describes the linear problem to be solved by the iterative solver.
PseudoBlockCGSolMgrOrthoFailure(const std::string &what_arg)
void STEQR(const char COMPZ, const OrdinalType n, ScalarType *D, ScalarType *E, ScalarType *Z, const OrdinalType ldz, ScalarType *WORK, OrdinalType *info) const
Type traits class that says whether Teuchos::LAPACK has a valid implementation for the given ScalarTy...
ReturnType
Whether the Belos solve converged for all linear systems.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > getResidualStatusTest() const
Return the residual status test.
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.
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.
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...