Stratimikos  Version of the Day
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 
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 {}
156 
157 
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
169  )
170 {
171  using Teuchos::as;
172  using Teuchos::TypeNameTraits;
173  using Teuchos::Exceptions::InvalidParameterType;
174  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
175 
176  this->setLinePrefix("BELOS/T");
177  // ToDo: Validate input
178  lp_ = lp;
179  solverPL_ = solverPL;
180  iterativeSolver_ = iterativeSolver;
181  fwdOpSrc_ = fwdOpSrc;
182  prec_ = prec;
183  isExternalPrec_ = isExternalPrec_in;
184  approxFwdOpSrc_ = approxFwdOpSrc;
185  supportSolveUse_ = supportSolveUse_in;
186  convergenceTestFrequency_ = convergenceTestFrequency;
187  // Check if "Convergence Tolerance" is in the solver parameter list. If
188  // not, use the default from the solver.
189  if ( !is_null(solverPL_) ) {
190  if (solverPL_->isParameter("Convergence Tolerance")) {
191 
192  // Stratimikos prefers tolerances as double, no matter the
193  // Scalar type. However, we also want it to accept the
194  // tolerance as magnitude_type, for example float if the Scalar
195  // type is float or std::complex<float>.
196  if (solverPL_->isType<double> ("Convergence Tolerance")) {
197  defaultTol_ =
198  as<magnitude_type> (solverPL_->get<double> ("Convergence Tolerance"));
199  }
200  else if (Teuchos::TypeTraits::is_same<double, magnitude_type>::value) {
201  // magnitude_type == double in this case, and we've already
202  // checked double above.
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).");
207  }
208  else if (solverPL_->isType<magnitude_type> ("Convergence Tolerance")) {
209  defaultTol_ = solverPL_->get<magnitude_type> ("Convergence Tolerance");
210  }
211  else {
212  // Throwing InvalidParameterType ensures that the exception's
213  // type is consistent both with what this method would have
214  // thrown before for an unrecognized type, and with what the
215  // user expects in general when the parameter doesn't have the
216  // right type.
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.");
224  }
225  }
226  }
227  else {
228  RCP<const Teuchos::ParameterList> defaultPL =
229  iterativeSolver->getValidParameters();
230 
231  // Stratimikos prefers tolerances as double, no matter the
232  // Scalar type. However, we also want it to accept the
233  // tolerance as magnitude_type, for example float if the Scalar
234  // type is float or std::complex<float>.
235  if (defaultPL->isType<double> ("Convergence Tolerance")) {
236  defaultTol_ =
237  as<magnitude_type> (defaultPL->get<double> ("Convergence Tolerance"));
238  }
239  else if (Teuchos::TypeTraits::is_same<double, magnitude_type>::value) {
240  // magnitude_type == double in this case, and we've already
241  // checked double above.
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).");
246  }
247  else if (defaultPL->isType<magnitude_type> ("Convergence Tolerance")) {
248  defaultTol_ = defaultPL->get<magnitude_type> ("Convergence Tolerance");
249  }
250  else {
251  // Throwing InvalidParameterType ensures that the exception's
252  // type is consistent both with what this method would have
253  // thrown before for an unrecognized type, and with what the
254  // user expects in general when the parameter doesn't have the
255  // right type.
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.");
263  }
264  }
265 }
266 
267 
268 template<class Scalar>
269 RCP<const LinearOpSourceBase<Scalar> >
271 {
272  RCP<const LinearOpSourceBase<Scalar> >
273  _fwdOpSrc = fwdOpSrc_;
274  fwdOpSrc_ = Teuchos::null;
275  return _fwdOpSrc;
276 }
277 
278 
279 template<class Scalar>
280 RCP<const PreconditionerBase<Scalar> >
282 {
283  RCP<const PreconditionerBase<Scalar> >
284  _prec = prec_;
285  prec_ = Teuchos::null;
286  return _prec;
287 }
288 
289 
290 template<class Scalar>
292 {
293  return isExternalPrec_;
294 }
295 
296 
297 template<class Scalar>
298 RCP<const LinearOpSourceBase<Scalar> >
300 {
301  RCP<const LinearOpSourceBase<Scalar> >
302  _approxFwdOpSrc = approxFwdOpSrc_;
303  approxFwdOpSrc_ = Teuchos::null;
304  return _approxFwdOpSrc;
305 }
306 
307 
308 template<class Scalar>
310 {
311  return supportSolveUse_;
312 }
313 
314 
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
325  )
326 {
327  if (lp) *lp = lp_;
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_;
335 
336  lp_ = Teuchos::null;
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;
344 }
345 
346 
347 // Overridden from LinearOpBase
348 
349 
350 template<class Scalar>
351 RCP< const VectorSpaceBase<Scalar> >
353 {
354  if (!is_null(lp_))
355  return lp_->getOperator()->range();
356  return Teuchos::null;
357 }
358 
359 
360 template<class Scalar>
361 RCP< const VectorSpaceBase<Scalar> >
363 {
364  if (!is_null(lp_))
365  return lp_->getOperator()->domain();
366  return Teuchos::null;
367 }
368 
369 
370 template<class Scalar>
371 RCP<const LinearOpBase<Scalar> >
373 {
374  return Teuchos::null; // Not supported yet but could be
375 }
376 
377 
378 // Overridden from Teuchos::Describable
379 
380 
381 template<class Scalar>
383 {
384  std::ostringstream oss;
385  oss << Teuchos::Describable::description();
386  if ( !is_null(lp_) && !is_null(lp_->getOperator()) ) {
387  oss << "{";
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()<<"\'";
394  oss << "}";
395  }
396  // ToDo: Make Belos::SolverManager derive from Teuchos::Describable so
397  // that we can get better information.
398  return oss.str();
399 }
400 
401 
402 template<class Scalar>
404  Teuchos::FancyOStream &out_arg,
405  const Teuchos::EVerbosityLevel verbLevel
406  ) const
407 {
408  using Teuchos::FancyOStream;
409  using Teuchos::OSTab;
410  using Teuchos::describe;
411  RCP<FancyOStream> out = rcp(&out_arg,false);
412  OSTab tab(out);
413  switch (verbLevel) {
414  case Teuchos::VERB_LOW:
415  break;
416  case Teuchos::VERB_DEFAULT:
417  case Teuchos::VERB_MEDIUM:
418  *out << this->description() << std::endl;
419  break;
420  case Teuchos::VERB_HIGH:
421  case Teuchos::VERB_EXTREME:
422  {
423  *out
424  << Teuchos::Describable::description()<< "{"
425  << "rangeDim=" << this->range()->dim()
426  << ",domainDim=" << this->domain()->dim() << "}\n";
427  if (lp_->getOperator().get()) {
428  OSTab tab1(out);
429  *out
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);
436  }
437  break;
438  }
439  default:
440  TEUCHOS_TEST_FOR_EXCEPT(true); // Should never get here!
441  }
442 }
443 
444 
445 // protected
446 
447 
448 // Overridden from LinearOpBase
449 
450 
451 template<class Scalar>
453 {
454  return ::Thyra::opSupported(*lp_->getOperator(),M_trans);
455 }
456 
457 
458 template<class Scalar>
460  const EOpTransp M_trans,
461  const MultiVectorBase<Scalar> &X,
462  const Ptr<MultiVectorBase<Scalar> > &Y,
463  const Scalar alpha,
464  const Scalar beta
465  ) const
466 {
467  ::Thyra::apply<Scalar>(*lp_->getOperator(), M_trans, X, Y, alpha, beta);
468 }
469 
470 
471 // Overridden from LinearOpWithSolveBase
472 
473 
474 template<class Scalar>
475 bool
477 {
478  return solveSupportsNewImpl(M_trans, Teuchos::null);
479 }
480 
481 
482 template<class Scalar>
483 bool
485  const Ptr<const SolveCriteria<Scalar> > solveCriteria) const
486 {
487  // Only support forward solve right now!
488  if (real_trans(transp)==NOTRANS) return true;
489  return false; // ToDo: Support adjoint solves!
490  // Otherwise, Thyra/Belos now supports every solve criteria type that exists
491  // because of the class Thyra::GeneralSolveCriteriaBelosStatusTest!
492  return true;
493 /*
494  if (real_trans(M_trans)==NOTRANS) {
495  return (
496  solveMeasureType.useDefault()
497  ||
498  solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_RHS)
499  ||
500  solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL,SOLVE_MEASURE_NORM_INIT_RESIDUAL)
501  );
502  }
503 */
504 }
505 
506 
507 template<class Scalar>
508 bool
510  EOpTransp M_trans, const SolveMeasureType& solveMeasureType) const
511 {
512  SolveCriteria<Scalar> solveCriteria(solveMeasureType, SolveCriteria<Scalar>::unspecifiedTolerance());
513  return solveSupportsNewImpl(M_trans, Teuchos::constOptInArg(solveCriteria));
514 }
515 
516 
517 template<class Scalar>
518 SolveStatus<Scalar>
520  const EOpTransp M_trans,
521  const MultiVectorBase<Scalar> &B,
522  const Ptr<MultiVectorBase<Scalar> > &X,
523  const Ptr<const SolveCriteria<Scalar> > solveCriteria
524  ) const
525 {
526 
527  THYRA_FUNC_TIME_MONITOR("Stratimikos: BelosLOWS");
528 
529  using Teuchos::rcp;
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);
541 
542  assertSolveSupports(*this, M_trans, solveCriteria);
543  // 2010/08/22: rabartl: Bug 4915 ToDo: Move the above into the NIV function
544  // solve(...).
545 
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";
551  OSTab tab2(out);
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";
555  }
556 
557  //
558  // Set RHS and LHS
559  //
560 
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!"
565  );
566 
567  //
568  // Set the solution criteria
569  //
570 
571  // Parameter list for the current solve.
572  const RCP<ParameterList> tmpPL = Teuchos::parameterList();
573 
574  // The solver's valid parameter list.
575  RCP<const ParameterList> validPL = iterativeSolver_->getValidParameters();
576 
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_);
584  }
585  else if (solveMeasureType(SOLVE_MEASURE_NORM_RESIDUAL, SOLVE_MEASURE_NORM_RHS)) {
586  if (requestedTol != SolveCriteria<Scalar>::unspecifiedTolerance()) {
587  tmpPL->set("Convergence Tolerance", requestedTol);
588  }
589  else {
590  tmpPL->set("Convergence Tolerance", defaultTol_);
591  }
592  setResidualScalingType (tmpPL, validPL, "Norm of RHS");
593  }
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);
597  }
598  else {
599  tmpPL->set("Convergence Tolerance", defaultTol_);
600  }
601  setResidualScalingType (tmpPL, validPL, "Norm of Initial Residual");
602  }
603  else {
604  // Set the most generic (and inefficient) solve criteria
605  generalSolveCriteriaBelosStatusTest = createGeneralSolveCriteriaBelosStatusTest(
606  *solveCriteria, convergenceTestFrequency_);
607  // Set the verbosity level (one level down)
608  generalSolveCriteriaBelosStatusTest->setOStream(out);
609  generalSolveCriteriaBelosStatusTest->setVerbLevel(incrVerbLevel(verbLevel, -1));
610  // Set the default convergence tolerance to always converged to allow
611  // the above status test to control things.
612  tmpPL->set("Convergence Tolerance", 1.0);
613  }
614  // maximum iterations
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"));
618  }
619  }
620  // If a preconditioner is on the left, then the implicit residual test
621  // scaling should be the preconditioned initial residual.
622  if (Teuchos::nonnull(lp_->getLeftPrec()) &&
623  validPL->isParameter ("Implicit Residual Scaling"))
624  tmpPL->set("Implicit Residual Scaling",
625  "Norm of Preconditioned Initial Residual");
626  }
627  else {
628  // No solveCriteria was even passed in!
629  tmpPL->set("Convergence Tolerance", defaultTol_);
630  }
631 
632  //
633  // Solve the linear system
634  //
635 
636  Belos::ReturnType belosSolveStatus;
637  {
638  // Write detailed convergence information if requested for levels >= VERB_LOW
639  RCP<std::ostream>
640  outUsed =
641  ( static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW)
642  ? out
643  : rcp(new FancyOStream(rcp(new Teuchos::oblackholestream())))
644  );
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);
650  }
651  try {
652  belosSolveStatus = iterativeSolver_->solve();
653  }
654  catch (Belos::BelosError&) {
655  belosSolveStatus = Belos::Unconverged;
656  }
657  }
658 
659  //
660  // Report the solve status
661  //
662 
663  totalTimer.stop();
664 
665  SolveStatus<Scalar> solveStatus;
666 
667  switch (belosSolveStatus) {
668  case Belos::Unconverged: {
669  solveStatus.solveStatus = SOLVE_STATUS_UNCONVERGED;
670  // Set achievedTol even if the solver did not converge. This is
671  // helpful for things like nonlinear solvers, which might be
672  // able to use a partially converged result, and which would
673  // like to know the achieved convergence tolerance for use in
674  // computing bounds. It's also helpful for estimating whether a
675  // small increase in the maximum iteration count might be
676  // helpful next time.
677  try {
678  // Some solvers might not have implemented achievedTol().
679  // The default implementation throws std::runtime_error.
680  solveStatus.achievedTol = iterativeSolver_->achievedTol();
681  } catch (std::runtime_error&) {
682  // Do nothing; use the default value of achievedTol.
683  }
684  break;
685  }
686  case Belos::Converged: {
687  solveStatus.solveStatus = SOLVE_STATUS_CONVERGED;
688  if (nonnull(generalSolveCriteriaBelosStatusTest)) {
689  // The user set a custom status test. This means that we
690  // should ask the custom status test itself, rather than the
691  // Belos solver, what the final achieved convergence tolerance
692  // was.
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]);
698  }
699  }
700  else {
701  try {
702  // Some solvers might not have implemented achievedTol().
703  // The default implementation throws std::runtime_error.
704  solveStatus.achievedTol = iterativeSolver_->achievedTol();
705  } catch (std::runtime_error&) {
706  // Use the default convergence tolerance. This is a correct
707  // upper bound, since we did actually converge.
708  solveStatus.achievedTol = tmpPL->get("Convergence Tolerance", defaultTol_);
709  }
710  }
711  break;
712  }
713  TEUCHOS_SWITCH_DEFAULT_DEBUG_ASSERT();
714  }
715 
716  std::ostringstream ossmessage;
717  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";
724 
725  solveStatus.message = ossmessage.str();
726 
727  // Dump the getNumIters() and the achieved convergence tolerance
728  // into solveStatus.extraParameters, as the "Belos/Iteration Count"
729  // resp. "Belos/Achieved Tolerance" parameters.
730  if (solveStatus.extraParameters.is_null()) {
731  solveStatus.extraParameters = parameterList ();
732  }
733  solveStatus.extraParameters->set ("Belos/Iteration Count",
734  iterativeSolver_->getNumIters());\
735  // package independent version of the same
736  solveStatus.extraParameters->set ("Iteration Count",
737  iterativeSolver_->getNumIters());\
738  // NOTE (mfh 13 Dec 2011) Though the most commonly used Belos
739  // solvers do implement achievedTol(), some Belos solvers currently
740  // do not. In the latter case, if the solver did not converge, the
741  // reported achievedTol() value may just be the default "invalid"
742  // value -1, and if the solver did converge, the reported value will
743  // just be the convergence tolerance (a correct upper bound).
744  solveStatus.extraParameters->set ("Belos/Achieved Tolerance",
745  solveStatus.achievedTol);
746 
747 // This information is in the previous line, which is printed anytime the verbosity
748 // is not set to Teuchos::VERB_NONE, so I'm commenting this out for now.
749 // if (out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_NONE))
750 // *out << "\nTotal solve time in Belos = "<<totalTimer.totalElapsedTime()<<" sec\n";
751 
752  return solveStatus;
753 
754 }
755 
756 
757 } // end namespace Thyra
758 
759 
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
virtual bool solveSupportsSolveMeasureTypeImpl(EOpTransp M_trans, const SolveMeasureType &solveMeasureType) 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.
virtual bool opSupportedImpl(EOpTransp M_trans) const

Generated for Stratimikos by doxygen 1.8.14