Stratimikos Version of the Day
Loading...
Searching...
No Matches
Thyra_BelosLinearOpWithSolve_def.hpp
1/*
2// @HEADER
3// ***********************************************************************
4//
5// Stratimikos: Thyra-based strategies for linear solvers
6// Copyright (2006) Sandia Corporation
7//
8// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9// license for use of this work by or on behalf of the U.S. Government.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
39//
40// ***********************************************************************
41// @HEADER
42*/
43
44
45#ifndef THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP
46#define THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP
47
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 "Teuchos_TypeTraits.hpp"
55#include "Stratimikos_Config.h"
56#ifdef HAVE_STRATIMIKOS_THYRATPETRAADAPTERS
57# include "Thyra_TpetraThyraWrappers.hpp"
58# include <MatrixMarket_Tpetra.hpp>
59# include <TpetraExt_MatrixMatrix.hpp>
60#endif
61
62namespace {
63 // Set the Belos solver's parameter list to scale its residual norms
64 // in the specified way.
65 //
66 // We break this out in a separate function because the parameters
67 // to set depend on which parameters the Belos solver supports. Not
68 // all Belos solvers support both the "Implicit Residual Scaling"
69 // and "Explicit Residual Scaling" parameters, so we have to check
70 // the solver's list of valid parameters for the existence of these.
71 //
72 // Scaling options: Belos lets you decide whether the solver will
73 // scale residual norms by the (left-)preconditioned initial
74 // residual norms (residualScalingType = "Norm of Initial
75 // Residual"), or by the unpreconditioned initial residual norms
76 // (residualScalingType = "Norm of RHS"). Usually you want to scale
77 // by the unpreconditioned initial residual norms. This is because
78 // preconditioning is just an optimization, and you really want to
79 // make ||B - AX|| small, rather than ||M B - M (A X)||. If you're
80 // measuring ||B - AX|| and scaling by the initial residual, you
81 // should use the unpreconditioned initial residual to match it.
82 //
83 // Note, however, that the implicit residual test computes
84 // left-preconditioned residuals, if a left preconditioner was
85 // provided. That's OK because when Belos solvers (at least the
86 // GMRES variants) are given a left preconditioner, they first check
87 // the implicit residuals. If those converge, they then check the
88 // explicit residuals. The explicit residual test does _not_ apply
89 // the left preconditioner when computing the residual. The
90 // implicit residual test is just an optimization so that Belos
91 // doesn't have to compute explicit residuals B - A*X at every
92 // iteration. This is why we use the same scaling factor for both
93 // the implicit and explicit residuals.
94 //
95 // Arguments:
96 //
97 // solverParams [in/out] Parameters for the current solve.
98 //
99 // solverValidParams [in] Valid parameter list for the Belos solver.
100 // Result of calling the solver's getValidParameters() method.
101 //
102 // residualScalingType [in] String describing how the solver should
103 // scale residuals. Valid values include "Norm of RHS" and "Norm
104 // of Initial Residual" (these are the only two options this file
105 // currently uses, though Belos offers other options).
106 void
107 setResidualScalingType (const Teuchos::RCP<Teuchos::ParameterList>& solverParams,
108 const Teuchos::RCP<const Teuchos::ParameterList>& solverValidParams,
109 const std::string& residualScalingType)
110 {
111 // Many Belos solvers (especially the GMRES variants) define both
112 // "Implicit Residual Scaling" and "Explicit Residual Scaling"
113 // options.
114 //
115 // "Implicit" means "the left-preconditioned approximate
116 // a.k.a. 'recursive' residual as computed by the Krylov method."
117 //
118 // "Explicit" means ||B - A*X||, the unpreconditioned, "exact"
119 // residual.
120 //
121 // Belos' GMRES implementations chain these two tests in sequence.
122 // Implicit comes first, and explicit is not evaluated unless
123 // implicit passes. In some cases (e.g., no left preconditioner),
124 // GMRES _only_ uses the implicit tests. This means that only
125 // setting "Explicit Residual Scaling" won't change the solver's
126 // behavior. Stratimikos tends to prefer using a right
127 // preconditioner, in which case setting only the "Explicit
128 // Residual Scaling" argument has no effect. Furthermore, if
129 // "Explicit Residual Scaling" is set to something other than the
130 // default (initial residual norm), without "Implicit Residual
131 // Scaling" getting the same setting, then the implicit residual
132 // test will be using a radically different scaling factor than
133 // the user wanted.
134 //
135 // Not all Belos solvers support both options. We check the
136 // solver's valid parameter list first before attempting to set
137 // the option.
138 if (solverValidParams->isParameter ("Implicit Residual Scaling")) {
139 solverParams->set ("Implicit Residual Scaling", residualScalingType);
140 }
141 if (solverValidParams->isParameter ("Explicit Residual Scaling")) {
142 solverParams->set ("Explicit Residual Scaling", residualScalingType);
143 }
144 }
145
146} // namespace (anonymous)
147
148
149namespace Thyra {
150
151
152// Constructors/initializers/accessors
153
154
155template<class Scalar>
157 :convergenceTestFrequency_(-1),
158 isExternalPrec_(false),
159 supportSolveUse_(SUPPORT_SOLVE_UNSPECIFIED),
160 defaultTol_ (-1.0),
161 label_(""),
162 filenameLHS_(""),
163 filenameRHS_(""),
164 counter_(0)
165{}
166
167
168template<class Scalar>
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
179 )
180{
181 using Teuchos::as;
182 using Teuchos::TypeNameTraits;
183 using Teuchos::Exceptions::InvalidParameterType;
184 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
185
186 this->setLinePrefix("BELOS/T");
187 // ToDo: Validate input
188 lp_ = lp;
189 solverPL_ = solverPL;
190 iterativeSolver_ = iterativeSolver;
191 fwdOpSrc_ = fwdOpSrc;
192 prec_ = prec;
193 isExternalPrec_ = isExternalPrec_in;
194 approxFwdOpSrc_ = approxFwdOpSrc;
195 supportSolveUse_ = supportSolveUse_in;
196 convergenceTestFrequency_ = convergenceTestFrequency;
197 // Check if "Convergence Tolerance" is in the solver parameter list. If
198 // not, use the default from the solver.
199 if ( !is_null(solverPL_) ) {
200 if (solverPL_->isParameter("Convergence Tolerance")) {
201
202 // Stratimikos prefers tolerances as double, no matter the
203 // Scalar type. However, we also want it to accept the
204 // tolerance as magnitude_type, for example float if the Scalar
205 // type is float or std::complex<float>.
206 if (solverPL_->isType<double> ("Convergence Tolerance")) {
207 defaultTol_ =
208 as<magnitude_type> (solverPL_->get<double> ("Convergence Tolerance"));
209 }
210 else if (Teuchos::TypeTraits::is_same<double, magnitude_type>::value) {
211 // magnitude_type == double in this case, and we've already
212 // checked double above.
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).");
217 }
218 else if (solverPL_->isType<magnitude_type> ("Convergence Tolerance")) {
219 defaultTol_ = solverPL_->get<magnitude_type> ("Convergence Tolerance");
220 }
221 else {
222 // Throwing InvalidParameterType ensures that the exception's
223 // type is consistent both with what this method would have
224 // thrown before for an unrecognized type, and with what the
225 // user expects in general when the parameter doesn't have the
226 // right type.
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.");
234 }
235 }
236
237 if (solverPL_->isParameter("Timer Label") && solverPL_->isType<std::string>("Timer Label")) {
238 label_ = solverPL_->get<std::string>("Timer Label");
239 lp_->setLabel(label_);
240 }
241 if (solverPL_->isParameter("Filename LHS") && solverPL_->isType<std::string>("Filename LHS")) {
242 filenameLHS_ = solverPL_->get<std::string>("Filename LHS");
243 }
244 if (solverPL_->isParameter("Filename RHS") && solverPL_->isType<std::string>("Filename RHS")) {
245 filenameRHS_ = solverPL_->get<std::string>("Filename RHS");
246 }
247 }
248 else {
249 RCP<const Teuchos::ParameterList> defaultPL =
250 iterativeSolver->getValidParameters();
251
252 // Stratimikos prefers tolerances as double, no matter the
253 // Scalar type. However, we also want it to accept the
254 // tolerance as magnitude_type, for example float if the Scalar
255 // type is float or std::complex<float>.
256 if (defaultPL->isType<double> ("Convergence Tolerance")) {
257 defaultTol_ =
258 as<magnitude_type> (defaultPL->get<double> ("Convergence Tolerance"));
259 }
260 else if (Teuchos::TypeTraits::is_same<double, magnitude_type>::value) {
261 // magnitude_type == double in this case, and we've already
262 // checked double above.
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).");
267 }
268 else if (defaultPL->isType<magnitude_type> ("Convergence Tolerance")) {
269 defaultTol_ = defaultPL->get<magnitude_type> ("Convergence Tolerance");
270 }
271 else {
272 // Throwing InvalidParameterType ensures that the exception's
273 // type is consistent both with what this method would have
274 // thrown before for an unrecognized type, and with what the
275 // user expects in general when the parameter doesn't have the
276 // right type.
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.");
284 }
285 }
286}
287
288
289template<class Scalar>
290RCP<const LinearOpSourceBase<Scalar> >
292{
293 RCP<const LinearOpSourceBase<Scalar> >
294 _fwdOpSrc = fwdOpSrc_;
295 fwdOpSrc_ = Teuchos::null;
296 return _fwdOpSrc;
297}
298
299
300template<class Scalar>
301RCP<const PreconditionerBase<Scalar> >
303{
304 RCP<const PreconditionerBase<Scalar> >
305 _prec = prec_;
306 prec_ = Teuchos::null;
307 return _prec;
308}
309
310
311template<class Scalar>
313{
314 return isExternalPrec_;
315}
316
317
318template<class Scalar>
319RCP<const LinearOpSourceBase<Scalar> >
321{
322 RCP<const LinearOpSourceBase<Scalar> >
323 _approxFwdOpSrc = approxFwdOpSrc_;
324 approxFwdOpSrc_ = Teuchos::null;
325 return _approxFwdOpSrc;
326}
327
328
329template<class Scalar>
331{
332 return supportSolveUse_;
333}
334
335
336template<class Scalar>
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
346 )
347{
348 if (lp) *lp = lp_;
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_;
356
357 lp_ = Teuchos::null;
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;
365}
366
367
368// Overridden from LinearOpBase
369
370
371template<class Scalar>
372RCP< const VectorSpaceBase<Scalar> >
374{
375 if (!is_null(lp_))
376 return lp_->getOperator()->range();
377 return Teuchos::null;
378}
379
380
381template<class Scalar>
382RCP< const VectorSpaceBase<Scalar> >
384{
385 if (!is_null(lp_))
386 return lp_->getOperator()->domain();
387 return Teuchos::null;
388}
389
390
391template<class Scalar>
392RCP<const LinearOpBase<Scalar> >
394{
395 return Teuchos::null; // Not supported yet but could be
396}
397
398
399// Overridden from Teuchos::Describable
400
401
402template<class Scalar>
404{
405 std::ostringstream oss;
406 oss << Teuchos::Describable::description();
407 if ( !is_null(lp_) && !is_null(lp_->getOperator()) ) {
408 oss << "{";
409 oss << "iterativeSolver=\'"<<iterativeSolver_->description()<<"\'";
410 oss << ",fwdOp=\'"<<lp_->getOperator()->description()<<"\'";
411 if (lp_->getLeftPrec().get())
412 oss << ",leftPrecOp=\'"<<lp_->getLeftPrec()->description()<<"\'";
413 if (lp_->getRightPrec().get())
414 oss << ",rightPrecOp=\'"<<lp_->getRightPrec()->description()<<"\'";
415 oss << "}";
416 }
417 // ToDo: Make Belos::SolverManager derive from Teuchos::Describable so
418 // that we can get better information.
419 return oss.str();
420}
421
422
423template<class Scalar>
425 Teuchos::FancyOStream &out_arg,
426 const Teuchos::EVerbosityLevel verbLevel
427 ) const
428{
429 using Teuchos::FancyOStream;
430 using Teuchos::OSTab;
431 using Teuchos::describe;
432 RCP<FancyOStream> out = rcp(&out_arg,false);
433 OSTab tab(out);
434 switch (verbLevel) {
435 case Teuchos::VERB_LOW:
436 break;
437 case Teuchos::VERB_DEFAULT:
438 case Teuchos::VERB_MEDIUM:
439 *out << this->description() << std::endl;
440 break;
441 case Teuchos::VERB_HIGH:
442 case Teuchos::VERB_EXTREME:
443 {
444 *out
445 << Teuchos::Describable::description()<< "{"
446 << "rangeDim=" << this->range()->dim()
447 << ",domainDim=" << this->domain()->dim() << "}\n";
448 if (lp_->getOperator().get()) {
449 OSTab tab1(out);
450 *out
451 << "iterativeSolver = "<<describe(*iterativeSolver_,verbLevel)
452 << "fwdOp = " << describe(*lp_->getOperator(),verbLevel);
453 if (lp_->getLeftPrec().get())
454 *out << "leftPrecOp = "<<describe(*lp_->getLeftPrec(),verbLevel);
455 if (lp_->getRightPrec().get())
456 *out << "rightPrecOp = "<<describe(*lp_->getRightPrec(),verbLevel);
457 }
458 break;
459 }
460 default:
461 TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
462 }
463}
464
465
466// protected
467
468
469// Overridden from LinearOpBase
470
471
472template<class Scalar>
474{
475 return ::Thyra::opSupported(*lp_->getOperator(),M_trans);
476}
477
478
479template<class Scalar>
481 const EOpTransp M_trans,
482 const MultiVectorBase<Scalar> &X,
483 const Ptr<MultiVectorBase<Scalar> > &Y,
484 const Scalar alpha,
485 const Scalar beta
486 ) const
487{
488 ::Thyra::apply<Scalar>(*lp_->getOperator(), M_trans, X, Y, alpha, beta);
489}
490
491
492// Overridden from LinearOpWithSolveBase
493
494
495template<class Scalar>
496bool
498{
499 return solveSupportsNewImpl(M_trans, Teuchos::null);
500}
501
502
503template<class Scalar>
504bool
506 const Ptr<const SolveCriteria<Scalar> > /* solveCriteria */) const
507{
508 // Only support forward solve right now!
509 if (real_trans(transp)==NOTRANS) return true;
510 return false; // ToDo: Support adjoint solves!
511 // Otherwise, Thyra/Belos now supports every solve criteria type that exists
512 // because of the class Thyra::GeneralSolveCriteriaBelosStatusTest!
513/*
514 if (real_trans(M_trans)==NOTRANS) {
515 return (
516 solveMeasureType.useDefault()
517 ||
518 solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS)
519 ||
520 solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_INIT_RESIDUAL)
521 );
522 }
523*/
524}
525
526
527template<class Scalar>
528bool
530 EOpTransp M_trans, const SolveMeasureType& solveMeasureType) const
531{
532 SolveCriteria<Scalar> solveCriteria(solveMeasureType, SolveCriteria<Scalar>::unspecifiedTolerance());
533 return solveSupportsNewImpl(M_trans, Teuchos::constOptInArg(solveCriteria));
534}
535
536
537template<class Scalar>
538SolveStatus<Scalar>
540 const EOpTransp M_trans,
541 const MultiVectorBase<Scalar> &B,
542 const Ptr<MultiVectorBase<Scalar> > &X,
543 const Ptr<const SolveCriteria<Scalar> > solveCriteria
544 ) const
545{
546
547 THYRA_FUNC_TIME_MONITOR("Stratimikos: BelosLOWS");
548
549 using Teuchos::rcp;
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);
561
562 assertSolveSupports(*this, M_trans, solveCriteria);
563 // 2010/08/22: rabartl: Bug 4915 ToDo: Move the above into the NIV function
564 // solve(...).
565
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";
571 OSTab tab2(out);
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";
575 }
576
577#ifdef HAVE_STRATIMIKOS_THYRATPETRAADAPTERS
578 //
579 // Write RHS and LHS to file if desired
580 //
581 if (filenameLHS_ != "") {
582 try {
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";
587 }
588 }
589 if (filenameRHS_ != "") {
590 try {
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";
595 }
596 }
597 ++counter_;
598#endif
599
600 //
601 // Set RHS and LHS
602 //
603
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!"
608 );
609
610 //
611 // Set the solution criteria
612 //
613
614 // Parameter list for the current solve.
615 const RCP<ParameterList> tmpPL = Teuchos::parameterList();
616
617 // The solver's valid parameter list.
618 RCP<const ParameterList> validPL = iterativeSolver_->getValidParameters();
619
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_);
627 }
628 else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_RHS)) {
629 if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
630 tmpPL->set("Convergence Tolerance", requestedTol);
631 }
632 else {
633 tmpPL->set("Convergence Tolerance", defaultTol_);
634 }
635 setResidualScalingType (tmpPL, validPL, "Norm of RHS");
636 }
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);
640 }
641 else {
642 tmpPL->set("Convergence Tolerance", defaultTol_);
643 }
644 setResidualScalingType (tmpPL, validPL, "Norm of Initial Residual");
645 }
646 else {
647 // Set the most generic (and inefficient) solve criteria
648 generalSolveCriteriaBelosStatusTest = createGeneralSolveCriteriaBelosStatusTest(
649 *solveCriteria, convergenceTestFrequency_);
650 // Set the verbosity level (one level down)
651 generalSolveCriteriaBelosStatusTest->setOStream(out);
652 generalSolveCriteriaBelosStatusTest->setVerbLevel(incrVerbLevel(verbLevel, -1));
653 // Set the default convergence tolerance to always converged to allow
654 // the above status test to control things.
655 tmpPL->set("Convergence Tolerance", 1.0);
656 }
657 // maximum iterations
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"));
661 }
662 }
663 // If a preconditioner is on the left, then the implicit residual test
664 // scaling should be the preconditioned initial residual.
665 if (Teuchos::nonnull(lp_->getLeftPrec()) &&
666 validPL->isParameter ("Implicit Residual Scaling"))
667 tmpPL->set("Implicit Residual Scaling",
668 "Norm of Preconditioned Initial Residual");
669 }
670 else {
671 // No solveCriteria was even passed in!
672 tmpPL->set("Convergence Tolerance", defaultTol_);
673 }
674
675 //
676 // Solve the linear system
677 //
678
679 Belos::ReturnType belosSolveStatus;
680 {
681 // Write detailed convergence information if requested for levels >= VERB_LOW
682 RCP<std::ostream>
683 outUsed =
684 ( static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)
685 ? out
686 : rcp(new FancyOStream(rcp(new Teuchos::oblackholestream())))
687 );
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);
693 }
694 try {
695 belosSolveStatus = iterativeSolver_->solve();
696 }
697 catch (Belos::BelosError& ex) {
698 TEUCHOS_TEST_FOR_EXCEPTION( true,
699 CatastrophicSolveFailure,
700 ex.what() );
701 }
702 }
703
704 //
705 // Report the solve status
706 //
707
708 totalTimer.stop();
709
710 SolveStatus<Scalar> solveStatus;
711
712 switch (belosSolveStatus) {
713 case Belos::Unconverged: {
714 solveStatus.solveStatus = SOLVE_STATUS_UNCONVERGED;
715 // Set achievedTol even if the solver did not converge. This is
716 // helpful for things like nonlinear solvers, which might be
717 // able to use a partially converged result, and which would
718 // like to know the achieved convergence tolerance for use in
719 // computing bounds. It's also helpful for estimating whether a
720 // small increase in the maximum iteration count might be
721 // helpful next time.
722 try {
723 // Some solvers might not have implemented achievedTol().
724 // The default implementation throws std::runtime_error.
725 solveStatus.achievedTol = iterativeSolver_->achievedTol();
726 } catch (std::runtime_error&) {
727 // Do nothing; use the default value of achievedTol.
728 }
729 break;
730 }
731 case Belos::Converged: {
732 solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
733 if (nonnull(generalSolveCriteriaBelosStatusTest)) {
734 // The user set a custom status test. This means that we
735 // should ask the custom status test itself, rather than the
736 // Belos solver, what the final achieved convergence tolerance
737 // was.
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]);
743 }
744 }
745 else {
746 try {
747 // Some solvers might not have implemented achievedTol().
748 // The default implementation throws std::runtime_error.
749 solveStatus.achievedTol = iterativeSolver_->achievedTol();
750 } catch (std::runtime_error&) {
751 // Use the default convergence tolerance. This is a correct
752 // upper bound, since we did actually converge.
753 solveStatus.achievedTol = tmpPL->get("Convergence Tolerance", defaultTol_);
754 }
755 }
756 break;
757 }
758 TEUCHOS_SWITCH_DEFAULT_DEBUG_ASSERT();
759 }
760
761 std::ostringstream ossmessage;
762 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";
770
771 solveStatus.message = ossmessage.str();
772
773 // Dump the getNumIters() and the achieved convergence tolerance
774 // into solveStatus.extraParameters, as the "Belos/Iteration Count"
775 // resp. "Belos/Achieved Tolerance" parameters.
776 if (solveStatus.extraParameters.is_null()) {
777 solveStatus.extraParameters = parameterList ();
778 }
779 solveStatus.extraParameters->set ("Belos/Iteration Count",
780 iterativeSolver_->getNumIters());\
781 // package independent version of the same
782 solveStatus.extraParameters->set ("Iteration Count",
783 iterativeSolver_->getNumIters());\
784 // NOTE (mfh 13 Dec 2011) Though the most commonly used Belos
785 // solvers do implement achievedTol(), some Belos solvers currently
786 // do not. In the latter case, if the solver did not converge, the
787 // reported achievedTol() value may just be the default "invalid"
788 // value -1, and if the solver did converge, the reported value will
789 // just be the convergence tolerance (a correct upper bound).
790 solveStatus.extraParameters->set ("Belos/Achieved Tolerance",
791 solveStatus.achievedTol);
792
793// This information is in the previous line, which is printed anytime the verbosity
794// is not set to Teuchos::VERB_NONE, so I'm commenting this out for now.
795// if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE))
796// *out << "\nTotal solve time in Belos = "<<totalTimer.totalElapsedTime()<<" sec\n";
797
798 return solveStatus;
799
800}
801
802
803} // end namespace Thyra
804
805
806#endif // THYRA_BELOS_LINEAR_OP_WITH_SOLVE_HPP
Concrete LinearOpWithSolveBase subclass in terms of Belos.
RCP< const VectorSpaceBase< Scalar > > domain() const
RCP< const VectorSpaceBase< Scalar > > range() const
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.
NOTRANS

Generated for Stratimikos by doxygen 1.10.0