170 const RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > &lp,
171 const RCP<Teuchos::ParameterList> &solverPL,
172 const RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > &iterativeSolver,
173 const RCP<
const LinearOpSourceBase<Scalar> > &fwdOpSrc,
174 const RCP<
const PreconditionerBase<Scalar> > &prec,
175 const bool isExternalPrec_in,
176 const RCP<
const LinearOpSourceBase<Scalar> > &approxFwdOpSrc,
177 const ESupportSolveUse &supportSolveUse_in,
178 const int convergenceTestFrequency
182 using Teuchos::TypeNameTraits;
183 using Teuchos::Exceptions::InvalidParameterType;
184 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
186 this->setLinePrefix(
"BELOS/T");
189 solverPL_ = solverPL;
190 iterativeSolver_ = iterativeSolver;
191 fwdOpSrc_ = fwdOpSrc;
193 isExternalPrec_ = isExternalPrec_in;
194 approxFwdOpSrc_ = approxFwdOpSrc;
195 supportSolveUse_ = supportSolveUse_in;
196 convergenceTestFrequency_ = convergenceTestFrequency;
199 if ( !is_null(solverPL_) ) {
200 if (solverPL_->isParameter(
"Convergence Tolerance")) {
206 if (solverPL_->isType<
double> (
"Convergence Tolerance")) {
208 as<magnitude_type> (solverPL_->get<
double> (
"Convergence Tolerance"));
210 else if (Teuchos::TypeTraits::is_same<double, magnitude_type>::value) {
213 TEUCHOS_TEST_FOR_EXCEPTION(
214 true, std::invalid_argument,
"BelosLinearOpWithSolve::initialize: "
215 "The \"Convergence Tolerance\" parameter, which you provided, must "
216 "have type double (the type of the magnitude of Scalar = double).");
218 else if (solverPL_->isType<magnitude_type> (
"Convergence Tolerance")) {
219 defaultTol_ = solverPL_->get<magnitude_type> (
"Convergence Tolerance");
227 TEUCHOS_TEST_FOR_EXCEPTION(
228 true, InvalidParameterType,
"BelosLinearOpWithSolve::initialize: "
229 "The \"Convergence Tolerance\" parameter, which you provided, must "
230 "have type double (preferred) or the type of the magnitude of Scalar "
231 "= " << TypeNameTraits<Scalar>::name () <<
", which is " <<
232 TypeNameTraits<magnitude_type>::name () <<
" in this case. You can "
233 "find that type using Teuchos::ScalarTraits<Scalar>::magnitudeType.");
237 if (solverPL_->isParameter(
"Timer Label") && solverPL_->isType<std::string>(
"Timer Label")) {
238 label_ = solverPL_->get<std::string>(
"Timer Label");
239 lp_->setLabel(label_);
241 if (solverPL_->isParameter(
"Filename LHS") && solverPL_->isType<std::string>(
"Filename LHS")) {
242 filenameLHS_ = solverPL_->get<std::string>(
"Filename LHS");
244 if (solverPL_->isParameter(
"Filename RHS") && solverPL_->isType<std::string>(
"Filename RHS")) {
245 filenameRHS_ = solverPL_->get<std::string>(
"Filename RHS");
249 RCP<const Teuchos::ParameterList> defaultPL =
250 iterativeSolver->getValidParameters();
256 if (defaultPL->isType<
double> (
"Convergence Tolerance")) {
258 as<magnitude_type> (defaultPL->get<
double> (
"Convergence Tolerance"));
260 else if (Teuchos::TypeTraits::is_same<double, magnitude_type>::value) {
263 TEUCHOS_TEST_FOR_EXCEPTION(
264 true, std::invalid_argument,
"BelosLinearOpWithSolve::initialize: "
265 "The \"Convergence Tolerance\" parameter, which you provided, must "
266 "have type double (the type of the magnitude of Scalar = double).");
268 else if (defaultPL->isType<magnitude_type> (
"Convergence Tolerance")) {
269 defaultTol_ = defaultPL->get<magnitude_type> (
"Convergence Tolerance");
277 TEUCHOS_TEST_FOR_EXCEPTION(
278 true, InvalidParameterType,
"BelosLinearOpWithSolve::initialize: "
279 "The \"Convergence Tolerance\" parameter, which you provided, must "
280 "have type double (preferred) or the type of the magnitude of Scalar "
281 "= " << TypeNameTraits<Scalar>::name () <<
", which is " <<
282 TypeNameTraits<magnitude_type>::name () <<
" in this case. You can "
283 "find that type using Teuchos::ScalarTraits<Scalar>::magnitudeType.");
338 RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > *lp,
339 RCP<Teuchos::ParameterList> *solverPL,
340 RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > *iterativeSolver,
341 RCP<
const LinearOpSourceBase<Scalar> > *fwdOpSrc,
342 RCP<
const PreconditionerBase<Scalar> > *prec,
343 bool *isExternalPrec_in,
344 RCP<
const LinearOpSourceBase<Scalar> > *approxFwdOpSrc,
345 ESupportSolveUse *supportSolveUse_in
349 if (solverPL) *solverPL = solverPL_;
350 if (iterativeSolver) *iterativeSolver = iterativeSolver_;
351 if (fwdOpSrc) *fwdOpSrc = fwdOpSrc_;
352 if (prec) *prec = prec_;
353 if (isExternalPrec_in) *isExternalPrec_in = isExternalPrec_;
354 if (approxFwdOpSrc) *approxFwdOpSrc = approxFwdOpSrc_;
355 if (supportSolveUse_in) *supportSolveUse_in = supportSolveUse_;
358 solverPL_ = Teuchos::null;
359 iterativeSolver_ = Teuchos::null;
360 fwdOpSrc_ = Teuchos::null;
361 prec_ = Teuchos::null;
362 isExternalPrec_ =
false;
363 approxFwdOpSrc_ = Teuchos::null;
364 supportSolveUse_ = SUPPORT_SOLVE_UNSPECIFIED;
540 const EOpTransp M_trans,
541 const MultiVectorBase<Scalar> &B,
542 const Ptr<MultiVectorBase<Scalar> > &X,
543 const Ptr<
const SolveCriteria<Scalar> > solveCriteria
547 THYRA_FUNC_TIME_MONITOR(
"Stratimikos: BelosLOWS");
550 using Teuchos::rcpFromRef;
551 using Teuchos::rcpFromPtr;
552 using Teuchos::FancyOStream;
553 using Teuchos::OSTab;
554 using Teuchos::ParameterList;
555 using Teuchos::parameterList;
556 using Teuchos::describe;
557 typedef Teuchos::ScalarTraits<Scalar>
ST;
558 typedef typename ST::magnitudeType ScalarMag;
559 Teuchos::Time totalTimer(
""), timer(
"");
560 totalTimer.start(
true);
562 assertSolveSupports(*
this, M_trans, solveCriteria);
566 const RCP<FancyOStream> out = this->getOStream();
567 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
568 OSTab tab = this->getOSTab();
569 if (out.get() &&
static_cast<int>(verbLevel) >
static_cast<int>(Teuchos::VERB_LOW)) {
570 *out <<
"\nStarting iterations with Belos:\n";
572 *out <<
"Using forward operator = " << describe(*fwdOpSrc_->getOp(),verbLevel);
573 *out <<
"Using iterative solver = " << describe(*iterativeSolver_,verbLevel);
574 *out <<
"With #Eqns="<<B.range()->dim()<<
", #RHSs="<<B.domain()->dim()<<
" ...\n";
577#ifdef HAVE_STRATIMIKOS_THYRATPETRAADAPTERS
581 if (filenameLHS_ !=
"") {
583 auto tmv = Thyra::TpetraOperatorVectorExtraction<Scalar>::getTpetraMultiVector(Teuchos::rcpFromPtr(X));
584 Tpetra::MatrixMarket::Writer<::Tpetra::CrsMatrix<Scalar> >::writeDenseFile(filenameLHS_+
"." + label_ +
"." + std::to_string(counter_), *tmv,
"",
"");
585 }
catch (
const std::logic_error&) {
586 *out <<
"Warning: Cannot write LHS multivector to file.\n";
589 if (filenameRHS_ !=
"") {
591 auto tmv = Thyra::TpetraOperatorVectorExtraction<Scalar>::getConstTpetraMultiVector(Teuchos::rcpFromRef(B));
592 Tpetra::MatrixMarket::Writer<::Tpetra::CrsMatrix<Scalar> >::writeDenseFile(filenameRHS_+
"." + label_ +
"." + std::to_string(counter_), *tmv,
"",
"");
593 }
catch (
const std::logic_error&) {
594 *out <<
"Warning: Cannot write RHS multivector to file.\n";
604 bool ret = lp_->setProblem( rcpFromPtr(X), rcpFromRef(B) );
605 TEUCHOS_TEST_FOR_EXCEPTION(
606 ret ==
false, CatastrophicSolveFailure
607 ,
"Error, the Belos::LinearProblem could not be set for the current solve!"
615 const RCP<ParameterList> tmpPL = Teuchos::parameterList();
618 RCP<const ParameterList> validPL = iterativeSolver_->getValidParameters();
620 SolveMeasureType solveMeasureType;
621 RCP<GeneralSolveCriteriaBelosStatusTest<Scalar> > generalSolveCriteriaBelosStatusTest;
622 if (nonnull(solveCriteria)) {
623 solveMeasureType = solveCriteria->solveMeasureType;
624 const ScalarMag requestedTol = solveCriteria->requestedTol;
625 if (solveMeasureType.useDefault()) {
626 tmpPL->set(
"Convergence Tolerance", defaultTol_);
628 else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_RHS)) {
629 if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
630 tmpPL->set(
"Convergence Tolerance", requestedTol);
633 tmpPL->set(
"Convergence Tolerance", defaultTol_);
635 setResidualScalingType (tmpPL, validPL,
"Norm of RHS");
637 else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_INIT_RESIDUAL)) {
638 if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
639 tmpPL->set(
"Convergence Tolerance", requestedTol);
642 tmpPL->set(
"Convergence Tolerance", defaultTol_);
644 setResidualScalingType (tmpPL, validPL,
"Norm of Initial Residual");
648 generalSolveCriteriaBelosStatusTest = createGeneralSolveCriteriaBelosStatusTest(
649 *solveCriteria, convergenceTestFrequency_);
651 generalSolveCriteriaBelosStatusTest->setOStream(out);
652 generalSolveCriteriaBelosStatusTest->setVerbLevel(incrVerbLevel(verbLevel, -1));
655 tmpPL->set(
"Convergence Tolerance", 1.0);
658 if (nonnull(solveCriteria->extraParameters)) {
659 if (Teuchos::isParameterType<int>(*solveCriteria->extraParameters,
"Maximum Iterations")) {
660 tmpPL->set(
"Maximum Iterations", Teuchos::get<int>(*solveCriteria->extraParameters,
"Maximum Iterations"));
665 if (Teuchos::nonnull(lp_->getLeftPrec()) &&
666 validPL->isParameter (
"Implicit Residual Scaling"))
667 tmpPL->set(
"Implicit Residual Scaling",
668 "Norm of Preconditioned Initial Residual");
672 tmpPL->set(
"Convergence Tolerance", defaultTol_);
679 Belos::ReturnType belosSolveStatus;
684 (
static_cast<int>(verbLevel) >=
static_cast<int>(Teuchos::VERB_LOW)
686 : rcp(
new FancyOStream(rcp(
new Teuchos::oblackholestream())))
688 Teuchos::OSTab tab1(outUsed,1,
"BELOS");
689 tmpPL->set(
"Output Stream", outUsed);
690 iterativeSolver_->setParameters(tmpPL);
691 if (nonnull(generalSolveCriteriaBelosStatusTest)) {
692 iterativeSolver_->setUserConvStatusTest(generalSolveCriteriaBelosStatusTest);
695 belosSolveStatus = iterativeSolver_->solve();
697 catch (Belos::BelosError& ex) {
698 TEUCHOS_TEST_FOR_EXCEPTION(
true,
699 CatastrophicSolveFailure,
710 SolveStatus<Scalar> solveStatus;
712 switch (belosSolveStatus) {
713 case Belos::Unconverged: {
714 solveStatus.solveStatus = SOLVE_STATUS_UNCONVERGED;
725 solveStatus.achievedTol = iterativeSolver_->achievedTol();
726 }
catch (std::runtime_error&) {
731 case Belos::Converged: {
732 solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
733 if (nonnull(generalSolveCriteriaBelosStatusTest)) {
738 const ArrayView<const ScalarMag> achievedTol =
739 generalSolveCriteriaBelosStatusTest->achievedTol();
740 solveStatus.achievedTol = Teuchos::ScalarTraits<ScalarMag>::zero();
741 for (Ordinal i = 0; i < achievedTol.size(); ++i) {
742 solveStatus.achievedTol = std::max(solveStatus.achievedTol, achievedTol[i]);
749 solveStatus.achievedTol = iterativeSolver_->achievedTol();
750 }
catch (std::runtime_error&) {
753 solveStatus.achievedTol = tmpPL->get(
"Convergence Tolerance", defaultTol_);
758 TEUCHOS_SWITCH_DEFAULT_DEBUG_ASSERT();
761 std::ostringstream ossmessage;
763 <<
"The Belos solver " << (label_ !=
"" ? (
"\"" + label_ +
"\" ") :
"")
764 <<
"of type \""<<iterativeSolver_->description()
765 <<
"\" returned a solve status of \""<< toString(solveStatus.solveStatus) <<
"\""
766 <<
" in " << iterativeSolver_->getNumIters() <<
" iterations"
767 <<
" with total CPU time of " << totalTimer.totalElapsedTime() <<
" sec" ;
768 if (out.get() &&
static_cast<int>(verbLevel) >
static_cast<int>(Teuchos::VERB_LOW))
769 *out <<
"\n" << ossmessage.str() <<
"\n";
771 solveStatus.message = ossmessage.str();
776 if (solveStatus.extraParameters.is_null()) {
777 solveStatus.extraParameters = parameterList ();
779 solveStatus.extraParameters->set (
"Belos/Iteration Count",
780 iterativeSolver_->getNumIters());\
782 solveStatus.extraParameters->set (
"Iteration Count",
783 iterativeSolver_->getNumIters());\
790 solveStatus.extraParameters->set (
"Belos/Achieved Tolerance",
791 solveStatus.achievedTol);