45 #ifndef THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP
46 #define THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP
48 #include "Thyra_BelosLinearOpWithSolve_decl.hpp"
49 #include "Thyra_GeneralSolveCriteriaBelosStatusTest.hpp"
50 #include "Thyra_LinearOpWithSolveHelpers.hpp"
51 #include "Teuchos_DebugDefaultAsserts.hpp"
52 #include "Teuchos_Assert.hpp"
53 #include "Teuchos_TimeMonitor.hpp"
54 #include "Stratimikos_Config.h"
55 #ifdef HAVE_STRATIMIKOS_THYRATPETRAADAPTERS
56 # include "Thyra_TpetraThyraWrappers.hpp"
58 # include <TpetraExt_MatrixMatrix.hpp>
108 const std::string& residualScalingType)
137 if (solverValidParams->
isParameter (
"Implicit Residual Scaling")) {
138 solverParams->
set (
"Implicit Residual Scaling", residualScalingType);
140 if (solverValidParams->
isParameter (
"Explicit Residual Scaling")) {
141 solverParams->
set (
"Explicit Residual Scaling", residualScalingType);
154 template<
class Scalar>
156 :convergenceTestFrequency_(-1),
157 isExternalPrec_(false),
167 template<
class Scalar>
174 const bool isExternalPrec_in,
177 const int convergenceTestFrequency
185 this->setLinePrefix(
"BELOS/T");
188 solverPL_ = solverPL;
189 iterativeSolver_ = iterativeSolver;
190 fwdOpSrc_ = fwdOpSrc;
192 isExternalPrec_ = isExternalPrec_in;
193 approxFwdOpSrc_ = approxFwdOpSrc;
194 supportSolveUse_ = supportSolveUse_in;
195 convergenceTestFrequency_ = convergenceTestFrequency;
199 if (solverPL_->isParameter(
"Convergence Tolerance")) {
205 if (solverPL_->isType<
double> (
"Convergence Tolerance")) {
207 as<magnitude_type> (solverPL_->get<
double> (
"Convergence Tolerance"));
209 else if (std::is_same_v<double, magnitude_type>) {
213 true, std::invalid_argument,
"BelosLinearOpWithSolve::initialize: "
214 "The \"Convergence Tolerance\" parameter, which you provided, must "
215 "have type double (the type of the magnitude of Scalar = double).");
217 else if (solverPL_->isType<magnitude_type> (
"Convergence Tolerance")) {
218 defaultTol_ = solverPL_->get<magnitude_type> (
"Convergence Tolerance");
227 true, InvalidParameterType,
"BelosLinearOpWithSolve::initialize: "
228 "The \"Convergence Tolerance\" parameter, which you provided, must "
229 "have type double (preferred) or the type of the magnitude of Scalar "
230 "= " << TypeNameTraits<Scalar>::name () <<
", which is " <<
231 TypeNameTraits<magnitude_type>::name () <<
" in this case. You can "
232 "find that type using Teuchos::ScalarTraits<Scalar>::magnitudeType.");
236 if (solverPL_->isParameter(
"Timer Label") && solverPL_->isType<std::string>(
"Timer Label")) {
237 label_ = solverPL_->get<std::string>(
"Timer Label");
238 lp_->setLabel(label_);
240 if (solverPL_->isParameter(
"Filename LHS") && solverPL_->isType<std::string>(
"Filename LHS")) {
241 filenameLHS_ = solverPL_->get<std::string>(
"Filename LHS");
243 if (solverPL_->isParameter(
"Filename RHS") && solverPL_->isType<std::string>(
"Filename RHS")) {
244 filenameRHS_ = solverPL_->get<std::string>(
"Filename RHS");
249 iterativeSolver->getValidParameters();
255 if (defaultPL->
isType<
double> (
"Convergence Tolerance")) {
257 as<magnitude_type> (defaultPL->
get<
double> (
"Convergence Tolerance"));
259 else if (std::is_same_v<double, magnitude_type>) {
263 true, std::invalid_argument,
"BelosLinearOpWithSolve::initialize: "
264 "The \"Convergence Tolerance\" parameter, which you provided, must "
265 "have type double (the type of the magnitude of Scalar = double).");
267 else if (defaultPL->
isType<magnitude_type> (
"Convergence Tolerance")) {
268 defaultTol_ = defaultPL->
get<magnitude_type> (
"Convergence Tolerance");
277 true, InvalidParameterType,
"BelosLinearOpWithSolve::initialize: "
278 "The \"Convergence Tolerance\" parameter, which you provided, must "
279 "have type double (preferred) or the type of the magnitude of Scalar "
280 "= " << TypeNameTraits<Scalar>::name () <<
", which is " <<
281 TypeNameTraits<magnitude_type>::name () <<
" in this case. You can "
282 "find that type using Teuchos::ScalarTraits<Scalar>::magnitudeType.");
288 template<
class Scalar>
293 _fwdOpSrc = fwdOpSrc_;
294 fwdOpSrc_ = Teuchos::null;
299 template<
class Scalar>
305 prec_ = Teuchos::null;
310 template<
class Scalar>
313 return isExternalPrec_;
317 template<
class Scalar>
322 _approxFwdOpSrc = approxFwdOpSrc_;
323 approxFwdOpSrc_ = Teuchos::null;
324 return _approxFwdOpSrc;
328 template<
class Scalar>
331 return supportSolveUse_;
335 template<
class Scalar>
342 bool *isExternalPrec_in,
348 if (solverPL) *solverPL = solverPL_;
349 if (iterativeSolver) *iterativeSolver = iterativeSolver_;
350 if (fwdOpSrc) *fwdOpSrc = fwdOpSrc_;
351 if (prec) *prec = prec_;
352 if (isExternalPrec_in) *isExternalPrec_in = isExternalPrec_;
353 if (approxFwdOpSrc) *approxFwdOpSrc = approxFwdOpSrc_;
354 if (supportSolveUse_in) *supportSolveUse_in = supportSolveUse_;
357 solverPL_ = Teuchos::null;
358 iterativeSolver_ = Teuchos::null;
359 fwdOpSrc_ = Teuchos::null;
360 prec_ = Teuchos::null;
361 isExternalPrec_ =
false;
362 approxFwdOpSrc_ = Teuchos::null;
370 template<
class Scalar>
375 return lp_->getOperator()->range();
376 return Teuchos::null;
380 template<
class Scalar>
385 return lp_->getOperator()->domain();
386 return Teuchos::null;
390 template<
class Scalar>
394 return Teuchos::null;
401 template<
class Scalar>
404 std::ostringstream oss;
408 oss <<
"iterativeSolver=\'"<<iterativeSolver_->description()<<
"\'";
409 oss <<
",fwdOp=\'"<<lp_->getOperator()->description()<<
"\'";
410 if (lp_->getLeftPrec().get())
411 oss <<
",leftPrecOp=\'"<<lp_->getLeftPrec()->description()<<
"\'";
412 if (lp_->getRightPrec().get())
413 oss <<
",rightPrecOp=\'"<<lp_->getRightPrec()->description()<<
"\'";
422 template<
class Scalar>
430 using Teuchos::describe;
437 *out << this->description() << std::endl;
444 <<
"rangeDim=" << this->range()->dim()
445 <<
",domainDim=" << this->domain()->dim() <<
"}\n";
446 if (lp_->getOperator().get()) {
449 <<
"iterativeSolver = "<<describe(*iterativeSolver_,verbLevel)
450 <<
"fwdOp = " << describe(*lp_->getOperator(),verbLevel);
451 if (lp_->getLeftPrec().get())
452 *out <<
"leftPrecOp = "<<describe(*lp_->getLeftPrec(),verbLevel);
453 if (lp_->getRightPrec().get())
454 *out <<
"rightPrecOp = "<<describe(*lp_->getRightPrec(),verbLevel);
470 template<
class Scalar>
473 return ::Thyra::opSupported(*lp_->getOperator(),M_trans);
477 template<
class Scalar>
486 ::Thyra::apply<Scalar>(*lp_->getOperator(), M_trans, X, Y, alpha, beta);
493 template<
class Scalar>
497 return solveSupportsNewImpl(M_trans, Teuchos::null);
501 template<
class Scalar>
525 template<
class Scalar>
531 return solveSupportsNewImpl(M_trans, Teuchos::constOptInArg(solveCriteria));
535 template<
class Scalar>
545 THYRA_FUNC_TIME_MONITOR(
"Stratimikos: BelosLOWS");
548 using Teuchos::rcpFromRef;
549 using Teuchos::rcpFromPtr;
553 using Teuchos::parameterList;
554 using Teuchos::describe;
556 typedef typename ST::magnitudeType ScalarMag;
558 totalTimer.start(
true);
560 assertSolveSupports(*
this, M_trans, solveCriteria);
566 OSTab tab = this->getOSTab();
568 *out <<
"\nStarting iterations with Belos:\n";
570 *out <<
"Using forward operator = " << describe(*fwdOpSrc_->getOp(),verbLevel);
571 *out <<
"Using iterative solver = " << describe(*iterativeSolver_,verbLevel);
572 *out <<
"With #Eqns="<<B.
range()->dim()<<
", #RHSs="<<B.
domain()->dim()<<
" ...\n";
575 #ifdef HAVE_STRATIMIKOS_THYRATPETRAADAPTERS
579 if (filenameLHS_ !=
"") {
583 }
catch (
const std::logic_error&) {
584 *out <<
"Warning: Cannot write LHS multivector to file.\n";
587 if (filenameRHS_ !=
"") {
591 }
catch (
const std::logic_error&) {
592 *out <<
"Warning: Cannot write RHS multivector to file.\n";
602 bool ret = lp_->setProblem( rcpFromPtr(X), rcpFromRef(B) );
605 ,
"Error, the Belos::LinearProblem could not be set for the current solve!"
621 solveMeasureType = solveCriteria->solveMeasureType;
622 const ScalarMag requestedTol = solveCriteria->requestedTol;
624 tmpPL->
set(
"Convergence Tolerance", defaultTol_);
628 tmpPL->
set(
"Convergence Tolerance", requestedTol);
631 tmpPL->
set(
"Convergence Tolerance", defaultTol_);
633 setResidualScalingType (tmpPL, validPL,
"Norm of RHS");
637 tmpPL->
set(
"Convergence Tolerance", requestedTol);
640 tmpPL->
set(
"Convergence Tolerance", defaultTol_);
642 setResidualScalingType (tmpPL, validPL,
"Norm of Initial Residual");
646 generalSolveCriteriaBelosStatusTest = createGeneralSolveCriteriaBelosStatusTest(
647 *solveCriteria, convergenceTestFrequency_);
649 generalSolveCriteriaBelosStatusTest->setOStream(out);
650 generalSolveCriteriaBelosStatusTest->setVerbLevel(
incrVerbLevel(verbLevel, -1));
653 tmpPL->
set(
"Convergence Tolerance", 1.0);
656 if (
nonnull(solveCriteria->extraParameters)) {
657 if (Teuchos::isParameterType<int>(*solveCriteria->extraParameters,
"Maximum Iterations")) {
658 tmpPL->
set(
"Maximum Iterations", Teuchos::get<int>(*solveCriteria->extraParameters,
"Maximum Iterations"));
664 validPL->
isParameter (
"Implicit Residual Scaling"))
665 tmpPL->
set(
"Implicit Residual Scaling",
666 "Norm of Preconditioned Initial Residual");
670 tmpPL->
set(
"Convergence Tolerance", defaultTol_);
687 tmpPL->
set(
"Output Stream", outUsed);
689 if (
nonnull(generalSolveCriteriaBelosStatusTest)) {
690 iterativeSolver_->setUserConvStatusTest(generalSolveCriteriaBelosStatusTest);
693 belosSolveStatus = iterativeSolver_->solve();
710 switch (belosSolveStatus) {
723 solveStatus.
achievedTol = iterativeSolver_->achievedTol();
724 }
catch (std::runtime_error&) {
731 if (
nonnull(generalSolveCriteriaBelosStatusTest)) {
736 const ArrayView<const ScalarMag> achievedTol =
737 generalSolveCriteriaBelosStatusTest->achievedTol();
739 for (
Ordinal i = 0; i < achievedTol.size(); ++i) {
747 solveStatus.
achievedTol = iterativeSolver_->achievedTol();
748 }
catch (std::runtime_error&) {
751 solveStatus.
achievedTol = tmpPL->
get(
"Convergence Tolerance", defaultTol_);
759 std::ostringstream ossmessage;
761 <<
"The Belos solver " << (label_ !=
"" ? (
"\"" + label_ +
"\" ") :
"")
762 <<
"of type \""<<iterativeSolver_->description()
764 <<
" in " << iterativeSolver_->getNumIters() <<
" iterations"
765 <<
" with total CPU time of " << totalTimer.totalElapsedTime() <<
" sec" ;
767 *out <<
"\n" << ossmessage.str() <<
"\n";
769 solveStatus.
message = ossmessage.str();
778 iterativeSolver_->getNumIters());\
781 iterativeSolver_->getNumIters());\
804 #endif // THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP
virtual bool solveSupportsImpl(EOpTransp M_trans) const
RCP< const VectorSpaceBase< Scalar > > domain() const
bool is_null(const boost::shared_ptr< T > &p)
virtual RCP< const VectorSpaceBase< Scalar > > range() const =0
#define TEUCHOS_SWITCH_DEFAULT_DEBUG_ASSERT()
basic_OSTab< char > OSTab
virtual SolveStatus< Scalar > solveImpl(const EOpTransp transp, const MultiVectorBase< Scalar > &B, const Ptr< MultiVectorBase< Scalar > > &X, const Ptr< const SolveCriteria< Scalar > > solveCriteria) const
RCP< const LinearOpSourceBase< Scalar > > extract_approxFwdOpSrc()
T & get(const std::string &name, T def_value)
void initialize(const RCP< Belos::LinearProblem< Scalar, MV_t, LO_t > > &lp, const RCP< Teuchos::ParameterList > &solverPL, const RCP< Belos::SolverManager< Scalar, MV_t, LO_t > > &iterativeSolver, const RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, const RCP< const PreconditionerBase< Scalar > > &prec, const bool isExternalPrec, const RCP< const LinearOpSourceBase< Scalar > > &approxFwdOpSrc, const ESupportSolveUse &supportSolveUse, const int convergenceTestFrequency)
Initializes given precreated solver objects.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
bool nonnull(const std::shared_ptr< T > &p)
basic_FancyOStream< char > FancyOStream
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
RCP< const PreconditionerBase< Scalar > > extract_prec()
virtual bool solveSupportsNewImpl(EOpTransp transp, const Ptr< const SolveCriteria< Scalar > > solveCriteria) const
EOpTransp real_trans(EOpTransp transp)
RCP< const VectorSpaceBase< Scalar > > range() const
virtual void applyImpl(const EOpTransp M_trans, const MultiVectorBase< Scalar > &X, const Ptr< MultiVectorBase< Scalar > > &Y, const Scalar alpha, const Scalar beta) const
bool isExternalPrec() const
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
bool isParameter(const std::string &name) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual std::string description() const
TEUCHOSCORE_LIB_DLL_EXPORT std::string toString(const EVerbosityLevel verbLevel)
std::string description() const
ESupportSolveUse supportSolveUse() const
ParameterList & setParameters(const ParameterList &source)
RCP< const LinearOpSourceBase< Scalar > > extract_fwdOpSrc()
TEUCHOSCORE_LIB_DLL_EXPORT EVerbosityLevel incrVerbLevel(const EVerbosityLevel inputVerbLevel, const int numLevels)
virtual RCP< const VectorSpaceBase< Scalar > > domain() const =0
bool nonnull(const boost::shared_ptr< T > &p)
TypeTo as(const TypeFrom &t)
bool isType(const std::string &name) const
BelosLinearOpWithSolve()
Construct to unintialize.
RCP< ParameterList > extraParameters
void uninitialize(RCP< Belos::LinearProblem< Scalar, MV_t, LO_t > > *lp=NULL, RCP< Teuchos::ParameterList > *solverPL=NULL, RCP< Belos::SolverManager< Scalar, MV_t, LO_t > > *iterativeSolver=NULL, RCP< const LinearOpSourceBase< Scalar > > *fwdOpSrc=NULL, RCP< const PreconditionerBase< Scalar > > *prec=NULL, bool *isExternalPrec=NULL, RCP< const LinearOpSourceBase< Scalar > > *approxFwdOpSrc=NULL, ESupportSolveUse *supportSolveUse=NULL)
Uninitializes and returns stored quantities.
RCP< const LinearOpBase< Scalar > > clone() const
virtual bool opSupportedImpl(EOpTransp M_trans) const
virtual bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)