Stratimikos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_BelosLinearOpWithSolve_def.hpp
Go to the documentation of this file.
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 
49 #include "Thyra_GeneralSolveCriteriaBelosStatusTest.hpp"
50 #include "Thyra_LinearOpWithSolveHelpers.hpp"
52 #include "Teuchos_Assert.hpp"
53 #include "Teuchos_TimeMonitor.hpp"
54 #include "Teuchos_TypeTraits.hpp"
55 
56 namespace {
57  // Set the Belos solver's parameter list to scale its residual norms
58  // in the specified way.
59  //
60  // We break this out in a separate function because the parameters
61  // to set depend on which parameters the Belos solver supports. Not
62  // all Belos solvers support both the "Implicit Residual Scaling"
63  // and "Explicit Residual Scaling" parameters, so we have to check
64  // the solver's list of valid parameters for the existence of these.
65  //
66  // Scaling options: Belos lets you decide whether the solver will
67  // scale residual norms by the (left-)preconditioned initial
68  // residual norms (residualScalingType = "Norm of Initial
69  // Residual"), or by the unpreconditioned initial residual norms
70  // (residualScalingType = "Norm of RHS"). Usually you want to scale
71  // by the unpreconditioned initial residual norms. This is because
72  // preconditioning is just an optimization, and you really want to
73  // make ||B - AX|| small, rather than ||M B - M (A X)||. If you're
74  // measuring ||B - AX|| and scaling by the initial residual, you
75  // should use the unpreconditioned initial residual to match it.
76  //
77  // Note, however, that the implicit residual test computes
78  // left-preconditioned residuals, if a left preconditioner was
79  // provided. That's OK because when Belos solvers (at least the
80  // GMRES variants) are given a left preconditioner, they first check
81  // the implicit residuals. If those converge, they then check the
82  // explicit residuals. The explicit residual test does _not_ apply
83  // the left preconditioner when computing the residual. The
84  // implicit residual test is just an optimization so that Belos
85  // doesn't have to compute explicit residuals B - A*X at every
86  // iteration. This is why we use the same scaling factor for both
87  // the implicit and explicit residuals.
88  //
89  // Arguments:
90  //
91  // solverParams [in/out] Parameters for the current solve.
92  //
93  // solverValidParams [in] Valid parameter list for the Belos solver.
94  // Result of calling the solver's getValidParameters() method.
95  //
96  // residualScalingType [in] String describing how the solver should
97  // scale residuals. Valid values include "Norm of RHS" and "Norm
98  // of Initial Residual" (these are the only two options this file
99  // currently uses, though Belos offers other options).
100  void
101  setResidualScalingType (const Teuchos::RCP<Teuchos::ParameterList>& solverParams,
102  const Teuchos::RCP<const Teuchos::ParameterList>& solverValidParams,
103  const std::string& residualScalingType)
104  {
105  // Many Belos solvers (especially the GMRES variants) define both
106  // "Implicit Residual Scaling" and "Explicit Residual Scaling"
107  // options.
108  //
109  // "Implicit" means "the left-preconditioned approximate
110  // a.k.a. 'recursive' residual as computed by the Krylov method."
111  //
112  // "Explicit" means ||B - A*X||, the unpreconditioned, "exact"
113  // residual.
114  //
115  // Belos' GMRES implementations chain these two tests in sequence.
116  // Implicit comes first, and explicit is not evaluated unless
117  // implicit passes. In some cases (e.g., no left preconditioner),
118  // GMRES _only_ uses the implicit tests. This means that only
119  // setting "Explicit Residual Scaling" won't change the solver's
120  // behavior. Stratimikos tends to prefer using a right
121  // preconditioner, in which case setting only the "Explicit
122  // Residual Scaling" argument has no effect. Furthermore, if
123  // "Explicit Residual Scaling" is set to something other than the
124  // default (initial residual norm), without "Implicit Residual
125  // Scaling" getting the same setting, then the implicit residual
126  // test will be using a radically different scaling factor than
127  // the user wanted.
128  //
129  // Not all Belos solvers support both options. We check the
130  // solver's valid parameter list first before attempting to set
131  // the option.
132  if (solverValidParams->isParameter ("Implicit Residual Scaling")) {
133  solverParams->set ("Implicit Residual Scaling", residualScalingType);
134  }
135  if (solverValidParams->isParameter ("Explicit Residual Scaling")) {
136  solverParams->set ("Explicit Residual Scaling", residualScalingType);
137  }
138  }
139 
140 } // namespace (anonymous)
141 
142 
143 namespace Thyra {
144 
145 
146 // Constructors/initializers/accessors
147 
148 
149 template<class Scalar>
151  :convergenceTestFrequency_(-1),
152  isExternalPrec_(false),
153  supportSolveUse_(SUPPORT_SOLVE_UNSPECIFIED),
154  defaultTol_ (-1.0),
155  label_("")
156 {}
157 
158 
159 template<class Scalar>
161  const RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > &lp,
162  const RCP<Teuchos::ParameterList> &solverPL,
163  const RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > &iterativeSolver,
164  const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
165  const RCP<const PreconditionerBase<Scalar> > &prec,
166  const bool isExternalPrec_in,
167  const RCP<const LinearOpSourceBase<Scalar> > &approxFwdOpSrc,
168  const ESupportSolveUse &supportSolveUse_in,
169  const int convergenceTestFrequency
170  )
171 {
172  using Teuchos::as;
175  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
176 
177  this->setLinePrefix("BELOS/T");
178  // ToDo: Validate input
179  lp_ = lp;
180  solverPL_ = solverPL;
181  iterativeSolver_ = iterativeSolver;
182  fwdOpSrc_ = fwdOpSrc;
183  prec_ = prec;
184  isExternalPrec_ = isExternalPrec_in;
185  approxFwdOpSrc_ = approxFwdOpSrc;
186  supportSolveUse_ = supportSolveUse_in;
187  convergenceTestFrequency_ = convergenceTestFrequency;
188  // Check if "Convergence Tolerance" is in the solver parameter list. If
189  // not, use the default from the solver.
190  if ( !is_null(solverPL_) ) {
191  if (solverPL_->isParameter("Convergence Tolerance")) {
192 
193  // Stratimikos prefers tolerances as double, no matter the
194  // Scalar type. However, we also want it to accept the
195  // tolerance as magnitude_type, for example float if the Scalar
196  // type is float or std::complex<float>.
197  if (solverPL_->isType<double> ("Convergence Tolerance")) {
198  defaultTol_ =
199  as<magnitude_type> (solverPL_->get<double> ("Convergence Tolerance"));
200  }
202  // magnitude_type == double in this case, and we've already
203  // checked double above.
205  true, std::invalid_argument, "BelosLinearOpWithSolve::initialize: "
206  "The \"Convergence Tolerance\" parameter, which you provided, must "
207  "have type double (the type of the magnitude of Scalar = double).");
208  }
209  else if (solverPL_->isType<magnitude_type> ("Convergence Tolerance")) {
210  defaultTol_ = solverPL_->get<magnitude_type> ("Convergence Tolerance");
211  }
212  else {
213  // Throwing InvalidParameterType ensures that the exception's
214  // type is consistent both with what this method would have
215  // thrown before for an unrecognized type, and with what the
216  // user expects in general when the parameter doesn't have the
217  // right type.
219  true, InvalidParameterType, "BelosLinearOpWithSolve::initialize: "
220  "The \"Convergence Tolerance\" parameter, which you provided, must "
221  "have type double (preferred) or the type of the magnitude of Scalar "
222  "= " << TypeNameTraits<Scalar>::name () << ", which is " <<
223  TypeNameTraits<magnitude_type>::name () << " in this case. You can "
224  "find that type using Teuchos::ScalarTraits<Scalar>::magnitudeType.");
225  }
226  }
227 
228  if (solverPL_->isParameter("Timer Label") && solverPL_->isType<std::string>("Timer Label")) {
229  label_ = solverPL_->get<std::string>("Timer Label");
230  lp_->setLabel(label_);
231  }
232  }
233  else {
235  iterativeSolver->getValidParameters();
236 
237  // Stratimikos prefers tolerances as double, no matter the
238  // Scalar type. However, we also want it to accept the
239  // tolerance as magnitude_type, for example float if the Scalar
240  // type is float or std::complex<float>.
241  if (defaultPL->isType<double> ("Convergence Tolerance")) {
242  defaultTol_ =
243  as<magnitude_type> (defaultPL->get<double> ("Convergence Tolerance"));
244  }
246  // magnitude_type == double in this case, and we've already
247  // checked double above.
249  true, std::invalid_argument, "BelosLinearOpWithSolve::initialize: "
250  "The \"Convergence Tolerance\" parameter, which you provided, must "
251  "have type double (the type of the magnitude of Scalar = double).");
252  }
253  else if (defaultPL->isType<magnitude_type> ("Convergence Tolerance")) {
254  defaultTol_ = defaultPL->get<magnitude_type> ("Convergence Tolerance");
255  }
256  else {
257  // Throwing InvalidParameterType ensures that the exception's
258  // type is consistent both with what this method would have
259  // thrown before for an unrecognized type, and with what the
260  // user expects in general when the parameter doesn't have the
261  // right type.
263  true, InvalidParameterType, "BelosLinearOpWithSolve::initialize: "
264  "The \"Convergence Tolerance\" parameter, which you provided, must "
265  "have type double (preferred) or the type of the magnitude of Scalar "
266  "= " << TypeNameTraits<Scalar>::name () << ", which is " <<
267  TypeNameTraits<magnitude_type>::name () << " in this case. You can "
268  "find that type using Teuchos::ScalarTraits<Scalar>::magnitudeType.");
269  }
270  }
271 }
272 
273 
274 template<class Scalar>
277 {
279  _fwdOpSrc = fwdOpSrc_;
280  fwdOpSrc_ = Teuchos::null;
281  return _fwdOpSrc;
282 }
283 
284 
285 template<class Scalar>
288 {
290  _prec = prec_;
291  prec_ = Teuchos::null;
292  return _prec;
293 }
294 
295 
296 template<class Scalar>
298 {
299  return isExternalPrec_;
300 }
301 
302 
303 template<class Scalar>
306 {
308  _approxFwdOpSrc = approxFwdOpSrc_;
309  approxFwdOpSrc_ = Teuchos::null;
310  return _approxFwdOpSrc;
311 }
312 
313 
314 template<class Scalar>
316 {
317  return supportSolveUse_;
318 }
319 
320 
321 template<class Scalar>
323  RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > *lp,
324  RCP<Teuchos::ParameterList> *solverPL,
325  RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > *iterativeSolver,
326  RCP<const LinearOpSourceBase<Scalar> > *fwdOpSrc,
327  RCP<const PreconditionerBase<Scalar> > *prec,
328  bool *isExternalPrec_in,
329  RCP<const LinearOpSourceBase<Scalar> > *approxFwdOpSrc,
330  ESupportSolveUse *supportSolveUse_in
331  )
332 {
333  if (lp) *lp = lp_;
334  if (solverPL) *solverPL = solverPL_;
335  if (iterativeSolver) *iterativeSolver = iterativeSolver_;
336  if (fwdOpSrc) *fwdOpSrc = fwdOpSrc_;
337  if (prec) *prec = prec_;
338  if (isExternalPrec_in) *isExternalPrec_in = isExternalPrec_;
339  if (approxFwdOpSrc) *approxFwdOpSrc = approxFwdOpSrc_;
340  if (supportSolveUse_in) *supportSolveUse_in = supportSolveUse_;
341 
342  lp_ = Teuchos::null;
343  solverPL_ = Teuchos::null;
344  iterativeSolver_ = Teuchos::null;
345  fwdOpSrc_ = Teuchos::null;
346  prec_ = Teuchos::null;
347  isExternalPrec_ = false;
348  approxFwdOpSrc_ = Teuchos::null;
349  supportSolveUse_ = SUPPORT_SOLVE_UNSPECIFIED;
350 }
351 
352 
353 // Overridden from LinearOpBase
354 
355 
356 template<class Scalar>
359 {
360  if (!is_null(lp_))
361  return lp_->getOperator()->range();
362  return Teuchos::null;
363 }
364 
365 
366 template<class Scalar>
369 {
370  if (!is_null(lp_))
371  return lp_->getOperator()->domain();
372  return Teuchos::null;
373 }
374 
375 
376 template<class Scalar>
379 {
380  return Teuchos::null; // Not supported yet but could be
381 }
382 
383 
384 // Overridden from Teuchos::Describable
385 
386 
387 template<class Scalar>
389 {
390  std::ostringstream oss;
392  if ( !is_null(lp_) && !is_null(lp_->getOperator()) ) {
393  oss << "{";
394  oss << "iterativeSolver=\'"<<iterativeSolver_->description()<<"\'";
395  oss << ",fwdOp=\'"<<lp_->getOperator()->description()<<"\'";
396  if (lp_->getLeftPrec().get())
397  oss << ",leftPrecOp=\'"<<lp_->getLeftPrec()->description()<<"\'";
398  if (lp_->getRightPrec().get())
399  oss << ",rightPrecOp=\'"<<lp_->getRightPrec()->description()<<"\'";
400  oss << "}";
401  }
402  // ToDo: Make Belos::SolverManager derive from Teuchos::Describable so
403  // that we can get better information.
404  return oss.str();
405 }
406 
407 
408 template<class Scalar>
410  Teuchos::FancyOStream &out_arg,
411  const Teuchos::EVerbosityLevel verbLevel
412  ) const
413 {
414  using Teuchos::FancyOStream;
415  using Teuchos::OSTab;
416  using Teuchos::describe;
417  RCP<FancyOStream> out = rcp(&out_arg,false);
418  OSTab tab(out);
419  switch (verbLevel) {
420  case Teuchos::VERB_LOW:
421  break;
424  *out << this->description() << std::endl;
425  break;
426  case Teuchos::VERB_HIGH:
428  {
429  *out
431  << "rangeDim=" << this->range()->dim()
432  << ",domainDim=" << this->domain()->dim() << "}\n";
433  if (lp_->getOperator().get()) {
434  OSTab tab1(out);
435  *out
436  << "iterativeSolver = "<<describe(*iterativeSolver_,verbLevel)
437  << "fwdOp = " << describe(*lp_->getOperator(),verbLevel);
438  if (lp_->getLeftPrec().get())
439  *out << "leftPrecOp = "<<describe(*lp_->getLeftPrec(),verbLevel);
440  if (lp_->getRightPrec().get())
441  *out << "rightPrecOp = "<<describe(*lp_->getRightPrec(),verbLevel);
442  }
443  break;
444  }
445  default:
446  TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
447  }
448 }
449 
450 
451 // protected
452 
453 
454 // Overridden from LinearOpBase
455 
456 
457 template<class Scalar>
459 {
460  return ::Thyra::opSupported(*lp_->getOperator(),M_trans);
461 }
462 
463 
464 template<class Scalar>
466  const EOpTransp M_trans,
467  const MultiVectorBase<Scalar> &X,
468  const Ptr<MultiVectorBase<Scalar> > &Y,
469  const Scalar alpha,
470  const Scalar beta
471  ) const
472 {
473  ::Thyra::apply<Scalar>(*lp_->getOperator(), M_trans, X, Y, alpha, beta);
474 }
475 
476 
477 // Overridden from LinearOpWithSolveBase
478 
479 
480 template<class Scalar>
481 bool
483 {
484  return solveSupportsNewImpl(M_trans, Teuchos::null);
485 }
486 
487 
488 template<class Scalar>
489 bool
491  const Ptr<const SolveCriteria<Scalar> > /* solveCriteria */) const
492 {
493  // Only support forward solve right now!
494  if (real_trans(transp)==NOTRANS) return true;
495  return false; // ToDo: Support adjoint solves!
496  // Otherwise, Thyra/Belos now supports every solve criteria type that exists
497  // because of the class Thyra::GeneralSolveCriteriaBelosStatusTest!
498 /*
499  if (real_trans(M_trans)==NOTRANS) {
500  return (
501  solveMeasureType.useDefault()
502  ||
503  solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS)
504  ||
505  solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_INIT_RESIDUAL)
506  );
507  }
508 */
509 }
510 
511 
512 template<class Scalar>
513 bool
515  EOpTransp M_trans, const SolveMeasureType& solveMeasureType) const
516 {
517  SolveCriteria<Scalar> solveCriteria(solveMeasureType, SolveCriteria<Scalar>::unspecifiedTolerance());
518  return solveSupportsNewImpl(M_trans, Teuchos::constOptInArg(solveCriteria));
519 }
520 
521 
522 template<class Scalar>
523 SolveStatus<Scalar>
525  const EOpTransp M_trans,
526  const MultiVectorBase<Scalar> &B,
527  const Ptr<MultiVectorBase<Scalar> > &X,
528  const Ptr<const SolveCriteria<Scalar> > solveCriteria
529  ) const
530 {
531 
532  THYRA_FUNC_TIME_MONITOR("Stratimikos: BelosLOWS");
533 
534  using Teuchos::rcp;
535  using Teuchos::rcpFromRef;
536  using Teuchos::rcpFromPtr;
537  using Teuchos::FancyOStream;
538  using Teuchos::OSTab;
540  using Teuchos::parameterList;
541  using Teuchos::describe;
543  typedef typename ST::magnitudeType ScalarMag;
544  Teuchos::Time totalTimer(""), timer("");
545  totalTimer.start(true);
546 
547  assertSolveSupports(*this, M_trans, solveCriteria);
548  // 2010/08/22: rabartl: Bug 4915 ToDo: Move the above into the NIV function
549  // solve(...).
550 
551  const RCP<FancyOStream> out = this->getOStream();
552  const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
553  OSTab tab = this->getOSTab();
554  if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW)) {
555  *out << "\nStarting iterations with Belos:\n";
556  OSTab tab2(out);
557  *out << "Using forward operator = " << describe(*fwdOpSrc_->getOp(),verbLevel);
558  *out << "Using iterative solver = " << describe(*iterativeSolver_,verbLevel);
559  *out << "With #Eqns="<<B.range()->dim()<<", #RHSs="<<B.domain()->dim()<<" ...\n";
560  }
561 
562  //
563  // Set RHS and LHS
564  //
565 
566  bool ret = lp_->setProblem( rcpFromPtr(X), rcpFromRef(B) );
568  ret == false, CatastrophicSolveFailure
569  ,"Error, the Belos::LinearProblem could not be set for the current solve!"
570  );
571 
572  //
573  // Set the solution criteria
574  //
575 
576  // Parameter list for the current solve.
577  const RCP<ParameterList> tmpPL = Teuchos::parameterList();
578 
579  // The solver's valid parameter list.
580  RCP<const ParameterList> validPL = iterativeSolver_->getValidParameters();
581 
582  SolveMeasureType solveMeasureType;
583  RCP<GeneralSolveCriteriaBelosStatusTest<Scalar> > generalSolveCriteriaBelosStatusTest;
584  if (nonnull(solveCriteria)) {
585  solveMeasureType = solveCriteria->solveMeasureType;
586  const ScalarMag requestedTol = solveCriteria->requestedTol;
587  if (solveMeasureType.useDefault()) {
588  tmpPL->set("Convergence Tolerance", defaultTol_);
589  }
590  else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_RHS)) {
591  if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
592  tmpPL->set("Convergence Tolerance", requestedTol);
593  }
594  else {
595  tmpPL->set("Convergence Tolerance", defaultTol_);
596  }
597  setResidualScalingType (tmpPL, validPL, "Norm of RHS");
598  }
599  else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_INIT_RESIDUAL)) {
600  if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
601  tmpPL->set("Convergence Tolerance", requestedTol);
602  }
603  else {
604  tmpPL->set("Convergence Tolerance", defaultTol_);
605  }
606  setResidualScalingType (tmpPL, validPL, "Norm of Initial Residual");
607  }
608  else {
609  // Set the most generic (and inefficient) solve criteria
610  generalSolveCriteriaBelosStatusTest = createGeneralSolveCriteriaBelosStatusTest(
611  *solveCriteria, convergenceTestFrequency_);
612  // Set the verbosity level (one level down)
613  generalSolveCriteriaBelosStatusTest->setOStream(out);
614  generalSolveCriteriaBelosStatusTest->setVerbLevel(incrVerbLevel(verbLevel, -1));
615  // Set the default convergence tolerance to always converged to allow
616  // the above status test to control things.
617  tmpPL->set("Convergence Tolerance", 1.0);
618  }
619  // maximum iterations
620  if (nonnull(solveCriteria->extraParameters)) {
621  if (Teuchos::isParameterType<int>(*solveCriteria->extraParameters,"Maximum Iterations")) {
622  tmpPL->set("Maximum Iterations", Teuchos::get<int>(*solveCriteria->extraParameters,"Maximum Iterations"));
623  }
624  }
625  // If a preconditioner is on the left, then the implicit residual test
626  // scaling should be the preconditioned initial residual.
627  if (Teuchos::nonnull(lp_->getLeftPrec()) &&
628  validPL->isParameter ("Implicit Residual Scaling"))
629  tmpPL->set("Implicit Residual Scaling",
630  "Norm of Preconditioned Initial Residual");
631  }
632  else {
633  // No solveCriteria was even passed in!
634  tmpPL->set("Convergence Tolerance", defaultTol_);
635  }
636 
637  //
638  // Solve the linear system
639  //
640 
641  Belos::ReturnType belosSolveStatus;
642  {
643  // Write detailed convergence information if requested for levels >= VERB_LOW
645  outUsed =
646  ( static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)
647  ? out
649  );
650  Teuchos::OSTab tab1(outUsed,1,"BELOS");
651  tmpPL->set("Output Stream", outUsed);
652  iterativeSolver_->setParameters(tmpPL);
653  if (nonnull(generalSolveCriteriaBelosStatusTest)) {
654  iterativeSolver_->setUserConvStatusTest(generalSolveCriteriaBelosStatusTest);
655  }
656  try {
657  belosSolveStatus = iterativeSolver_->solve();
658  }
659  catch (Belos::BelosError& ex) {
661  CatastrophicSolveFailure,
662  ex.what() );
663  }
664  }
665 
666  //
667  // Report the solve status
668  //
669 
670  totalTimer.stop();
671 
672  SolveStatus<Scalar> solveStatus;
673 
674  switch (belosSolveStatus) {
675  case Belos::Unconverged: {
676  solveStatus.solveStatus = SOLVE_STATUS_UNCONVERGED;
677  // Set achievedTol even if the solver did not converge. This is
678  // helpful for things like nonlinear solvers, which might be
679  // able to use a partially converged result, and which would
680  // like to know the achieved convergence tolerance for use in
681  // computing bounds. It's also helpful for estimating whether a
682  // small increase in the maximum iteration count might be
683  // helpful next time.
684  try {
685  // Some solvers might not have implemented achievedTol().
686  // The default implementation throws std::runtime_error.
687  solveStatus.achievedTol = iterativeSolver_->achievedTol();
688  } catch (std::runtime_error&) {
689  // Do nothing; use the default value of achievedTol.
690  }
691  break;
692  }
693  case Belos::Converged: {
694  solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
695  if (nonnull(generalSolveCriteriaBelosStatusTest)) {
696  // The user set a custom status test. This means that we
697  // should ask the custom status test itself, rather than the
698  // Belos solver, what the final achieved convergence tolerance
699  // was.
700  const ArrayView<const ScalarMag> achievedTol =
701  generalSolveCriteriaBelosStatusTest->achievedTol();
702  solveStatus.achievedTol = ST::zero();
703  for (Ordinal i = 0; i < achievedTol.size(); ++i) {
704  solveStatus.achievedTol = std::max(solveStatus.achievedTol, achievedTol[i]);
705  }
706  }
707  else {
708  try {
709  // Some solvers might not have implemented achievedTol().
710  // The default implementation throws std::runtime_error.
711  solveStatus.achievedTol = iterativeSolver_->achievedTol();
712  } catch (std::runtime_error&) {
713  // Use the default convergence tolerance. This is a correct
714  // upper bound, since we did actually converge.
715  solveStatus.achievedTol = tmpPL->get("Convergence Tolerance", defaultTol_);
716  }
717  }
718  break;
719  }
721  }
722 
723  std::ostringstream ossmessage;
724  ossmessage
725  << "The Belos solver " << (label_ != "" ? ("\"" + label_ + "\" ") : "")
726  << "of type \""<<iterativeSolver_->description()
727  <<"\" returned a solve status of \""<< toString(solveStatus.solveStatus) << "\""
728  << " in " << iterativeSolver_->getNumIters() << " iterations"
729  << " with total CPU time of " << totalTimer.totalElapsedTime() << " sec" ;
730  if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW))
731  *out << "\n" << ossmessage.str() << "\n";
732 
733  solveStatus.message = ossmessage.str();
734 
735  // Dump the getNumIters() and the achieved convergence tolerance
736  // into solveStatus.extraParameters, as the "Belos/Iteration Count"
737  // resp. "Belos/Achieved Tolerance" parameters.
738  if (solveStatus.extraParameters.is_null()) {
739  solveStatus.extraParameters = parameterList ();
740  }
741  solveStatus.extraParameters->set ("Belos/Iteration Count",
742  iterativeSolver_->getNumIters());\
743  // package independent version of the same
744  solveStatus.extraParameters->set ("Iteration Count",
745  iterativeSolver_->getNumIters());\
746  // NOTE (mfh 13 Dec 2011) Though the most commonly used Belos
747  // solvers do implement achievedTol(), some Belos solvers currently
748  // do not. In the latter case, if the solver did not converge, the
749  // reported achievedTol() value may just be the default "invalid"
750  // value -1, and if the solver did converge, the reported value will
751  // just be the convergence tolerance (a correct upper bound).
752  solveStatus.extraParameters->set ("Belos/Achieved Tolerance",
753  solveStatus.achievedTol);
754 
755 // This information is in the previous line, which is printed anytime the verbosity
756 // is not set to Teuchos::VERB_NONE, so I'm commenting this out for now.
757 // if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE))
758 // *out << "\nTotal solve time in Belos = "<<totalTimer.totalElapsedTime()<<" sec\n";
759 
760  return solveStatus;
761 
762 }
763 
764 
765 } // end namespace Thyra
766 
767 
768 #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
double ST
T * get() 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
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
TypeTo as(const TypeFrom &t)
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
Teuchos_Ordinal Ordinal
virtual bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)