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" 51 #include "Teuchos_DebugDefaultAsserts.hpp" 52 #include "Teuchos_Assert.hpp" 53 #include "Teuchos_TimeMonitor.hpp" 54 #include "Teuchos_TypeTraits.hpp" 101 setResidualScalingType (
const Teuchos::RCP<Teuchos::ParameterList>& solverParams,
102 const Teuchos::RCP<const Teuchos::ParameterList>& solverValidParams,
103 const std::string& residualScalingType)
132 if (solverValidParams->isParameter (
"Implicit Residual Scaling")) {
133 solverParams->set (
"Implicit Residual Scaling", residualScalingType);
135 if (solverValidParams->isParameter (
"Explicit Residual Scaling")) {
136 solverParams->set (
"Explicit Residual Scaling", residualScalingType);
149 template<
class Scalar>
151 :convergenceTestFrequency_(-1),
152 isExternalPrec_(false),
153 supportSolveUse_(SUPPORT_SOLVE_UNSPECIFIED),
158 template<
class Scalar>
160 const RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > &lp,
161 const RCP<Teuchos::ParameterList> &solverPL,
162 const RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > &iterativeSolver,
163 const RCP<
const LinearOpSourceBase<Scalar> > &fwdOpSrc,
164 const RCP<
const PreconditionerBase<Scalar> > &prec,
165 const bool isExternalPrec_in,
166 const RCP<
const LinearOpSourceBase<Scalar> > &approxFwdOpSrc,
167 const ESupportSolveUse &supportSolveUse_in,
168 const int convergenceTestFrequency
172 using Teuchos::TypeNameTraits;
173 using Teuchos::Exceptions::InvalidParameterType;
174 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
176 this->setLinePrefix(
"BELOS/T");
179 solverPL_ = solverPL;
180 iterativeSolver_ = iterativeSolver;
181 fwdOpSrc_ = fwdOpSrc;
183 isExternalPrec_ = isExternalPrec_in;
184 approxFwdOpSrc_ = approxFwdOpSrc;
185 supportSolveUse_ = supportSolveUse_in;
186 convergenceTestFrequency_ = convergenceTestFrequency;
189 if ( !is_null(solverPL_) ) {
190 if (solverPL_->isParameter(
"Convergence Tolerance")) {
196 if (solverPL_->isType<
double> (
"Convergence Tolerance")) {
198 as<magnitude_type> (solverPL_->get<
double> (
"Convergence Tolerance"));
200 else if (Teuchos::TypeTraits::is_same<double, magnitude_type>::value) {
203 TEUCHOS_TEST_FOR_EXCEPTION(
204 true, std::invalid_argument,
"BelosLinearOpWithSolve::initialize: " 205 "The \"Convergence Tolerance\" parameter, which you provided, must " 206 "have type double (the type of the magnitude of Scalar = double).");
208 else if (solverPL_->isType<magnitude_type> (
"Convergence Tolerance")) {
209 defaultTol_ = solverPL_->get<magnitude_type> (
"Convergence Tolerance");
217 TEUCHOS_TEST_FOR_EXCEPTION(
218 true, InvalidParameterType,
"BelosLinearOpWithSolve::initialize: " 219 "The \"Convergence Tolerance\" parameter, which you provided, must " 220 "have type double (preferred) or the type of the magnitude of Scalar " 221 "= " << TypeNameTraits<Scalar>::name () <<
", which is " <<
222 TypeNameTraits<magnitude_type>::name () <<
" in this case. You can " 223 "find that type using Teuchos::ScalarTraits<Scalar>::magnitudeType.");
228 RCP<const Teuchos::ParameterList> defaultPL =
229 iterativeSolver->getValidParameters();
235 if (defaultPL->isType<
double> (
"Convergence Tolerance")) {
237 as<magnitude_type> (defaultPL->get<
double> (
"Convergence Tolerance"));
239 else if (Teuchos::TypeTraits::is_same<double, magnitude_type>::value) {
242 TEUCHOS_TEST_FOR_EXCEPTION(
243 true, std::invalid_argument,
"BelosLinearOpWithSolve::initialize: " 244 "The \"Convergence Tolerance\" parameter, which you provided, must " 245 "have type double (the type of the magnitude of Scalar = double).");
247 else if (defaultPL->isType<magnitude_type> (
"Convergence Tolerance")) {
248 defaultTol_ = defaultPL->get<magnitude_type> (
"Convergence Tolerance");
256 TEUCHOS_TEST_FOR_EXCEPTION(
257 true, InvalidParameterType,
"BelosLinearOpWithSolve::initialize: " 258 "The \"Convergence Tolerance\" parameter, which you provided, must " 259 "have type double (preferred) or the type of the magnitude of Scalar " 260 "= " << TypeNameTraits<Scalar>::name () <<
", which is " <<
261 TypeNameTraits<magnitude_type>::name () <<
" in this case. You can " 262 "find that type using Teuchos::ScalarTraits<Scalar>::magnitudeType.");
268 template<
class Scalar>
269 RCP<const LinearOpSourceBase<Scalar> >
272 RCP<const LinearOpSourceBase<Scalar> >
273 _fwdOpSrc = fwdOpSrc_;
274 fwdOpSrc_ = Teuchos::null;
279 template<
class Scalar>
280 RCP<const PreconditionerBase<Scalar> >
283 RCP<const PreconditionerBase<Scalar> >
285 prec_ = Teuchos::null;
290 template<
class Scalar>
293 return isExternalPrec_;
297 template<
class Scalar>
298 RCP<const LinearOpSourceBase<Scalar> >
301 RCP<const LinearOpSourceBase<Scalar> >
302 _approxFwdOpSrc = approxFwdOpSrc_;
303 approxFwdOpSrc_ = Teuchos::null;
304 return _approxFwdOpSrc;
308 template<
class Scalar>
311 return supportSolveUse_;
315 template<
class Scalar>
317 RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > *lp,
318 RCP<Teuchos::ParameterList> *solverPL,
319 RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > *iterativeSolver,
320 RCP<
const LinearOpSourceBase<Scalar> > *fwdOpSrc,
321 RCP<
const PreconditionerBase<Scalar> > *prec,
322 bool *isExternalPrec_in,
323 RCP<
const LinearOpSourceBase<Scalar> > *approxFwdOpSrc,
324 ESupportSolveUse *supportSolveUse_in
328 if (solverPL) *solverPL = solverPL_;
329 if (iterativeSolver) *iterativeSolver = iterativeSolver_;
330 if (fwdOpSrc) *fwdOpSrc = fwdOpSrc_;
331 if (prec) *prec = prec_;
332 if (isExternalPrec_in) *isExternalPrec_in = isExternalPrec_;
333 if (approxFwdOpSrc) *approxFwdOpSrc = approxFwdOpSrc_;
334 if (supportSolveUse_in) *supportSolveUse_in = supportSolveUse_;
337 solverPL_ = Teuchos::null;
338 iterativeSolver_ = Teuchos::null;
339 fwdOpSrc_ = Teuchos::null;
340 prec_ = Teuchos::null;
341 isExternalPrec_ =
false;
342 approxFwdOpSrc_ = Teuchos::null;
343 supportSolveUse_ = SUPPORT_SOLVE_UNSPECIFIED;
350 template<
class Scalar>
351 RCP< const VectorSpaceBase<Scalar> >
355 return lp_->getOperator()->range();
356 return Teuchos::null;
360 template<
class Scalar>
361 RCP< const VectorSpaceBase<Scalar> >
365 return lp_->getOperator()->domain();
366 return Teuchos::null;
370 template<
class Scalar>
371 RCP<const LinearOpBase<Scalar> >
374 return Teuchos::null;
381 template<
class Scalar>
384 std::ostringstream oss;
385 oss << Teuchos::Describable::description();
386 if ( !is_null(lp_) && !is_null(lp_->getOperator()) ) {
388 oss <<
"iterativeSolver=\'"<<iterativeSolver_->description()<<
"\'";
389 oss <<
",fwdOp=\'"<<lp_->getOperator()->description()<<
"\'";
390 if (lp_->getLeftPrec().get())
391 oss <<
",leftPrecOp=\'"<<lp_->getLeftPrec()->description()<<
"\'";
392 if (lp_->getRightPrec().get())
393 oss <<
",rightPrecOp=\'"<<lp_->getRightPrec()->description()<<
"\'";
402 template<
class Scalar>
404 Teuchos::FancyOStream &out_arg,
405 const Teuchos::EVerbosityLevel verbLevel
408 using Teuchos::FancyOStream;
409 using Teuchos::OSTab;
410 using Teuchos::describe;
411 RCP<FancyOStream> out = rcp(&out_arg,
false);
414 case Teuchos::VERB_LOW:
416 case Teuchos::VERB_DEFAULT:
417 case Teuchos::VERB_MEDIUM:
418 *out << this->description() << std::endl;
420 case Teuchos::VERB_HIGH:
421 case Teuchos::VERB_EXTREME:
424 << Teuchos::Describable::description()<<
"{" 425 <<
"rangeDim=" << this->range()->dim()
426 <<
",domainDim=" << this->domain()->dim() <<
"}\n";
427 if (lp_->getOperator().get()) {
430 <<
"iterativeSolver = "<<describe(*iterativeSolver_,verbLevel)
431 <<
"fwdOp = " << describe(*lp_->getOperator(),verbLevel);
432 if (lp_->getLeftPrec().get())
433 *out <<
"leftPrecOp = "<<describe(*lp_->getLeftPrec(),verbLevel);
434 if (lp_->getRightPrec().get())
435 *out <<
"rightPrecOp = "<<describe(*lp_->getRightPrec(),verbLevel);
440 TEUCHOS_TEST_FOR_EXCEPT(
true);
451 template<
class Scalar>
454 return ::Thyra::opSupported(*lp_->getOperator(),M_trans);
458 template<
class Scalar>
460 const EOpTransp M_trans,
461 const MultiVectorBase<Scalar> &X,
462 const Ptr<MultiVectorBase<Scalar> > &Y,
467 ::Thyra::apply<Scalar>(*lp_->getOperator(), M_trans, X, Y, alpha, beta);
474 template<
class Scalar>
478 return solveSupportsNewImpl(M_trans, Teuchos::null);
482 template<
class Scalar>
485 const Ptr<
const SolveCriteria<Scalar> > solveCriteria)
const 488 if (real_trans(transp)==
NOTRANS)
return true;
507 template<
class Scalar>
510 EOpTransp M_trans,
const SolveMeasureType& solveMeasureType)
const 512 SolveCriteria<Scalar> solveCriteria(solveMeasureType, SolveCriteria<Scalar>::unspecifiedTolerance());
513 return solveSupportsNewImpl(M_trans, Teuchos::constOptInArg(solveCriteria));
517 template<
class Scalar>
520 const EOpTransp M_trans,
521 const MultiVectorBase<Scalar> &B,
522 const Ptr<MultiVectorBase<Scalar> > &X,
523 const Ptr<
const SolveCriteria<Scalar> > solveCriteria
527 THYRA_FUNC_TIME_MONITOR(
"Stratimikos: BelosLOWS");
530 using Teuchos::rcpFromRef;
531 using Teuchos::rcpFromPtr;
532 using Teuchos::FancyOStream;
533 using Teuchos::OSTab;
534 using Teuchos::ParameterList;
535 using Teuchos::parameterList;
536 using Teuchos::describe;
537 typedef Teuchos::ScalarTraits<Scalar>
ST;
538 typedef typename ST::magnitudeType ScalarMag;
539 Teuchos::Time totalTimer(
""), timer(
"");
540 totalTimer.start(
true);
542 assertSolveSupports(*
this, M_trans, solveCriteria);
546 const RCP<FancyOStream> out = this->getOStream();
547 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
548 OSTab tab = this->getOSTab();
549 if (out.get() &&
static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW)) {
550 *out <<
"\nStarting iterations with Belos:\n";
552 *out <<
"Using forward operator = " << describe(*fwdOpSrc_->getOp(),verbLevel);
553 *out <<
"Using iterative solver = " << describe(*iterativeSolver_,verbLevel);
554 *out <<
"With #Eqns="<<B.range()->dim()<<
", #RHSs="<<B.domain()->dim()<<
" ...\n";
561 bool ret = lp_->setProblem( rcpFromPtr(X), rcpFromRef(B) );
562 TEUCHOS_TEST_FOR_EXCEPTION(
563 ret ==
false, CatastrophicSolveFailure
564 ,
"Error, the Belos::LinearProblem could not be set for the current solve!" 572 const RCP<ParameterList> tmpPL = Teuchos::parameterList();
575 RCP<const ParameterList> validPL = iterativeSolver_->getValidParameters();
577 SolveMeasureType solveMeasureType;
578 RCP<GeneralSolveCriteriaBelosStatusTest<Scalar> > generalSolveCriteriaBelosStatusTest;
579 if (nonnull(solveCriteria)) {
580 solveMeasureType = solveCriteria->solveMeasureType;
581 const ScalarMag requestedTol = solveCriteria->requestedTol;
582 if (solveMeasureType.useDefault()) {
583 tmpPL->set(
"Convergence Tolerance", defaultTol_);
585 else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_RHS)) {
586 if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
587 tmpPL->set(
"Convergence Tolerance", requestedTol);
590 tmpPL->set(
"Convergence Tolerance", defaultTol_);
592 setResidualScalingType (tmpPL, validPL,
"Norm of RHS");
594 else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_INIT_RESIDUAL)) {
595 if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
596 tmpPL->set(
"Convergence Tolerance", requestedTol);
599 tmpPL->set(
"Convergence Tolerance", defaultTol_);
601 setResidualScalingType (tmpPL, validPL,
"Norm of Initial Residual");
605 generalSolveCriteriaBelosStatusTest = createGeneralSolveCriteriaBelosStatusTest(
606 *solveCriteria, convergenceTestFrequency_);
608 generalSolveCriteriaBelosStatusTest->setOStream(out);
609 generalSolveCriteriaBelosStatusTest->setVerbLevel(incrVerbLevel(verbLevel, -1));
612 tmpPL->set(
"Convergence Tolerance", 1.0);
615 if (nonnull(solveCriteria->extraParameters)) {
616 if (Teuchos::isParameterType<int>(*solveCriteria->extraParameters,
"Maximum Iterations")) {
617 tmpPL->set(
"Maximum Iterations", Teuchos::get<int>(*solveCriteria->extraParameters,
"Maximum Iterations"));
622 if (Teuchos::nonnull(lp_->getLeftPrec()) &&
623 validPL->isParameter (
"Implicit Residual Scaling"))
624 tmpPL->set(
"Implicit Residual Scaling",
625 "Norm of Preconditioned Initial Residual");
629 tmpPL->set(
"Convergence Tolerance", defaultTol_);
636 Belos::ReturnType belosSolveStatus;
641 (
static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)
643 : rcp(
new FancyOStream(rcp(
new Teuchos::oblackholestream())))
645 Teuchos::OSTab tab1(outUsed,1,
"BELOS");
646 tmpPL->set(
"Output Stream", outUsed);
647 iterativeSolver_->setParameters(tmpPL);
648 if (nonnull(generalSolveCriteriaBelosStatusTest)) {
649 iterativeSolver_->setUserConvStatusTest(generalSolveCriteriaBelosStatusTest);
652 belosSolveStatus = iterativeSolver_->solve();
654 catch (Belos::BelosError&) {
655 belosSolveStatus = Belos::Unconverged;
665 SolveStatus<Scalar> solveStatus;
667 switch (belosSolveStatus) {
668 case Belos::Unconverged: {
669 solveStatus.solveStatus = SOLVE_STATUS_UNCONVERGED;
680 solveStatus.achievedTol = iterativeSolver_->achievedTol();
681 }
catch (std::runtime_error&) {
686 case Belos::Converged: {
687 solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
688 if (nonnull(generalSolveCriteriaBelosStatusTest)) {
693 const ArrayView<const ScalarMag> achievedTol =
694 generalSolveCriteriaBelosStatusTest->achievedTol();
695 solveStatus.achievedTol = ST::zero();
696 for (
Ordinal i = 0; i < achievedTol.size(); ++i) {
697 solveStatus.achievedTol = std::max(solveStatus.achievedTol, achievedTol[i]);
704 solveStatus.achievedTol = iterativeSolver_->achievedTol();
705 }
catch (std::runtime_error&) {
708 solveStatus.achievedTol = tmpPL->get(
"Convergence Tolerance", defaultTol_);
713 TEUCHOS_SWITCH_DEFAULT_DEBUG_ASSERT();
716 std::ostringstream ossmessage;
718 <<
"The Belos solver of type \""<<iterativeSolver_->description()
719 <<
"\" returned a solve status of \""<<
toString(solveStatus.solveStatus) <<
"\"" 720 <<
" in " << iterativeSolver_->getNumIters() <<
" iterations" 721 <<
" with total CPU time of " << totalTimer.totalElapsedTime() <<
" sec" ;
722 if (out.get() &&
static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW))
723 *out <<
"\n" << ossmessage.str() <<
"\n";
725 solveStatus.message = ossmessage.str();
730 if (solveStatus.extraParameters.is_null()) {
731 solveStatus.extraParameters = parameterList ();
733 solveStatus.extraParameters->set (
"Belos/Iteration Count",
734 iterativeSolver_->getNumIters());\
736 solveStatus.extraParameters->set (
"Iteration Count",
737 iterativeSolver_->getNumIters());\
744 solveStatus.extraParameters->set (
"Belos/Achieved Tolerance",
745 solveStatus.achievedTol);
760 #endif // THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP virtual bool solveSupportsNewImpl(EOpTransp transp, const Ptr< const SolveCriteria< Scalar > > solveCriteria) const
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel) const
RCP< const LinearOpSourceBase< Scalar > > extract_approxFwdOpSrc()
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.
RCP< const PreconditionerBase< Scalar > > extract_prec()
RCP< const VectorSpaceBase< Scalar > > domain() const
RCP< const LinearOpBase< Scalar > > clone() const
std::string description() const
virtual bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) const
bool isExternalPrec() const
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 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
RCP< const LinearOpSourceBase< Scalar > > extract_fwdOpSrc()
virtual bool solveSupportsImpl(EOpTransp M_trans) 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.
ESupportSolveUse supportSolveUse() const
const char * toString(const ESolverType solverType)
virtual bool opSupportedImpl(EOpTransp M_trans) const