45 #ifndef THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP
46 #define THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP
49 #include "Thyra_GeneralSolveCriteriaBelosStatusTest.hpp"
50 #include "Thyra_LinearOpWithSolveHelpers.hpp"
54 #include "Stratimikos_Config.h"
55 #ifdef HAVE_STRATIMIKOS_THYRATPETRAADAPTERS
57 # include <MatrixMarket_Tpetra.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),
158 supportSolveUse_(SUPPORT_SOLVE_UNSPECIFIED),
167 template<
class Scalar>
169 const RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > &lp,
171 const RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > &iterativeSolver,
172 const RCP<
const LinearOpSourceBase<Scalar> > &fwdOpSrc,
173 const RCP<
const PreconditionerBase<Scalar> > &prec,
174 const bool isExternalPrec_in,
175 const RCP<
const LinearOpSourceBase<Scalar> > &approxFwdOpSrc,
176 const ESupportSolveUse &supportSolveUse_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_;
299 template<
class Scalar>
310 template<
class Scalar>
313 return isExternalPrec_;
317 template<
class Scalar>
322 _approxFwdOpSrc = approxFwdOpSrc_;
324 return _approxFwdOpSrc;
328 template<
class Scalar>
331 return supportSolveUse_;
335 template<
class Scalar>
337 RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > *lp,
339 RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > *iterativeSolver,
340 RCP<
const LinearOpSourceBase<Scalar> > *fwdOpSrc,
341 RCP<
const PreconditionerBase<Scalar> > *prec,
342 bool *isExternalPrec_in,
343 RCP<
const LinearOpSourceBase<Scalar> > *approxFwdOpSrc,
344 ESupportSolveUse *supportSolveUse_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_;
361 isExternalPrec_ =
false;
363 supportSolveUse_ = SUPPORT_SOLVE_UNSPECIFIED;
370 template<
class Scalar>
375 return lp_->getOperator()->range();
380 template<
class Scalar>
385 return lp_->getOperator()->domain();
390 template<
class Scalar>
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>
479 const EOpTransp M_trans,
480 const MultiVectorBase<Scalar> &X,
481 const Ptr<MultiVectorBase<Scalar> > &Y,
486 ::Thyra::apply<Scalar>(*lp_->getOperator(), M_trans, X, Y, alpha, beta);
493 template<
class Scalar>
501 template<
class Scalar>
504 const Ptr<
const SolveCriteria<Scalar> > )
const
507 if (real_trans(transp)==
NOTRANS)
return true;
525 template<
class Scalar>
528 EOpTransp M_trans,
const SolveMeasureType& solveMeasureType)
const
530 SolveCriteria<Scalar> solveCriteria(solveMeasureType, SolveCriteria<Scalar>::unspecifiedTolerance());
531 return solveSupportsNewImpl(M_trans, Teuchos::constOptInArg(solveCriteria));
535 template<
class Scalar>
538 const EOpTransp M_trans,
539 const MultiVectorBase<Scalar> &
B,
540 const Ptr<MultiVectorBase<Scalar> > &X,
541 const Ptr<
const SolveCriteria<Scalar> > solveCriteria
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_ !=
"") {
582 Tpetra::MatrixMarket::Writer<::Tpetra::CrsMatrix<Scalar> >::writeDenseFile(filenameLHS_+
"." + label_ +
"." + std::to_string(counter_), *tmv,
"",
"");
583 }
catch (
const std::logic_error&) {
584 *out <<
"Warning: Cannot write LHS multivector to file.\n";
587 if (filenameRHS_ !=
"") {
590 Tpetra::MatrixMarket::Writer<::Tpetra::CrsMatrix<Scalar> >::writeDenseFile(filenameRHS_+
"." + label_ +
"." + std::to_string(counter_), *tmv,
"",
"");
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) );
604 ret ==
false, CatastrophicSolveFailure
605 ,
"Error, the Belos::LinearProblem could not be set for the current solve!"
618 SolveMeasureType solveMeasureType;
621 solveMeasureType = solveCriteria->solveMeasureType;
622 const ScalarMag requestedTol = solveCriteria->requestedTol;
623 if (solveMeasureType.useDefault()) {
624 tmpPL->
set(
"Convergence Tolerance", defaultTol_);
626 else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_RHS)) {
627 if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
628 tmpPL->
set(
"Convergence Tolerance", requestedTol);
631 tmpPL->
set(
"Convergence Tolerance", defaultTol_);
633 setResidualScalingType (tmpPL, validPL,
"Norm of RHS");
635 else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_INIT_RESIDUAL)) {
636 if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
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_);
677 Belos::ReturnType belosSolveStatus;
687 tmpPL->
set(
"Output Stream", outUsed);
689 if (
nonnull(generalSolveCriteriaBelosStatusTest)) {
690 iterativeSolver_->setUserConvStatusTest(generalSolveCriteriaBelosStatusTest);
693 belosSolveStatus = iterativeSolver_->solve();
695 catch (Belos::BelosError& ex) {
697 CatastrophicSolveFailure,
708 SolveStatus<Scalar> solveStatus;
710 switch (belosSolveStatus) {
711 case Belos::Unconverged: {
712 solveStatus.solveStatus = SOLVE_STATUS_UNCONVERGED;
723 solveStatus.achievedTol = iterativeSolver_->achievedTol();
724 }
catch (std::runtime_error&) {
729 case Belos::Converged: {
730 solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
731 if (
nonnull(generalSolveCriteriaBelosStatusTest)) {
736 const ArrayView<const ScalarMag> achievedTol =
737 generalSolveCriteriaBelosStatusTest->achievedTol();
739 for (
Ordinal i = 0; i < achievedTol.size(); ++i) {
740 solveStatus.achievedTol = std::max(solveStatus.achievedTol, achievedTol[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()
763 <<
"\" returned a solve status of \""<<
toString(solveStatus.solveStatus) <<
"\""
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();
774 if (solveStatus.extraParameters.is_null()) {
775 solveStatus.extraParameters = parameterList ();
777 solveStatus.extraParameters->set (
"Belos/Iteration Count",
778 iterativeSolver_->getNumIters());\
780 solveStatus.extraParameters->set (
"Iteration Count",
781 iterativeSolver_->getNumIters());\
788 solveStatus.extraParameters->set (
"Belos/Achieved Tolerance",
789 solveStatus.achievedTol);
804 #endif // THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP
virtual bool solveSupportsImpl(EOpTransp M_trans) const
RCP< const VectorSpaceBase< Scalar > > domain() const
const std::string toString(const EAdjointEpetraOp adjointEpetraOp)
bool is_null(const boost::shared_ptr< T > &p)
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(ParameterList &l, const std::string &name)
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)
#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
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
basic_OSTab< char > OSTab
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)
#define TEUCHOS_SWITCH_DEFAULT_DEBUG_ASSERT()
virtual std::string description() const
std::string description() const
TypeTo as(const TypeFrom &t)
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)
bool nonnull(const boost::shared_ptr< T > &p)
basic_FancyOStream< char > FancyOStream
bool isType(const std::string &name) const
BelosLinearOpWithSolve()
Construct to unintialize.
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)