Stratimikos  Version of the Day
Thyra_BelosLinearOpWithSolveFactory_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_FACTORY_HPP
46 #define THYRA_BELOS_LINEAR_OP_WITH_SOLVE_FACTORY_HPP
47 
48 
49 #include "Thyra_BelosLinearOpWithSolveFactory_decl.hpp"
50 #include "Thyra_BelosLinearOpWithSolve.hpp"
51 #include "Thyra_ScaledAdjointLinearOpBase.hpp"
52 
53 #include "BelosBlockGmresSolMgr.hpp"
54 #include "BelosPseudoBlockGmresSolMgr.hpp"
55 #include "BelosBlockCGSolMgr.hpp"
56 #include "BelosPseudoBlockCGSolMgr.hpp"
57 #include "BelosPseudoBlockStochasticCGSolMgr.hpp"
58 #include "BelosGCRODRSolMgr.hpp"
59 #include "BelosRCGSolMgr.hpp"
60 #include "BelosMinresSolMgr.hpp"
61 #include "BelosTFQMRSolMgr.hpp"
62 
63 #include "BelosThyraAdapter.hpp"
64 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
65 #include "Teuchos_StandardParameterEntryValidators.hpp"
66 #include "Teuchos_ParameterList.hpp"
67 #include "Teuchos_dyn_cast.hpp"
68 #include "Teuchos_ValidatorXMLConverterDB.hpp"
69 #include "Teuchos_StandardValidatorXMLConverters.hpp"
70 
71 
72 namespace Thyra {
73 
74 
75 // Parameter names for Parameter List
76 
77 template<class Scalar>
78 const std::string BelosLinearOpWithSolveFactory<Scalar>::SolverType_name = "Solver Type";
79 template<class Scalar>
80 const std::string BelosLinearOpWithSolveFactory<Scalar>::SolverType_default = "Pseudo Block GMRES";
81 template<class Scalar>
82 const std::string BelosLinearOpWithSolveFactory<Scalar>::SolverTypes_name = "Solver Types";
83 template<class Scalar>
84 const std::string BelosLinearOpWithSolveFactory<Scalar>::BlockGMRES_name = "Block GMRES";
85 template<class Scalar>
86 const std::string BelosLinearOpWithSolveFactory<Scalar>::PseudoBlockGMRES_name = "Pseudo Block GMRES";
87 template<class Scalar>
88 const std::string BelosLinearOpWithSolveFactory<Scalar>::BlockCG_name = "Block CG";
89 template<class Scalar>
90 const std::string BelosLinearOpWithSolveFactory<Scalar>::PseudoBlockCG_name = "Pseudo Block CG";
91 template<class Scalar>
92 const std::string BelosLinearOpWithSolveFactory<Scalar>::PseudoBlockStochasticCG_name = "Pseudo Block Stochastic CG";
93 template<class Scalar>
94 const std::string BelosLinearOpWithSolveFactory<Scalar>::GCRODR_name = "GCRODR";
95 template<class Scalar>
96 const std::string BelosLinearOpWithSolveFactory<Scalar>::RCG_name = "RCG";
97 template<class Scalar>
98 const std::string BelosLinearOpWithSolveFactory<Scalar>::MINRES_name = "MINRES";
99 template<class Scalar>
100 const std::string BelosLinearOpWithSolveFactory<Scalar>::TFQMR_name = "TFQMR";
101 template<class Scalar>
102 const std::string BelosLinearOpWithSolveFactory<Scalar>::ConvergenceTestFrequency_name = "Convergence Test Frequency";
103 
104 namespace {
105 const std::string LeftPreconditionerIfUnspecified_name = "Left Preconditioner If Unspecified";
106 }
107 
108 // Constructors/initializers/accessors
109 
110 
111 template<class Scalar>
113  :solverType_(SOLVER_TYPE_PSEUDO_BLOCK_GMRES),
114  convergenceTestFrequency_(1)
115 {
116  updateThisValidParamList();
117 }
118 
119 
120 template<class Scalar>
122  const RCP<PreconditionerFactoryBase<Scalar> > &precFactory
123  )
124  :solverType_(SOLVER_TYPE_PSEUDO_BLOCK_GMRES)
125 {
126  this->setPreconditionerFactory(precFactory, "");
127 }
128 
129 
130 // Overridden from LinearOpWithSolveFactoryBase
131 
132 
133 template<class Scalar>
135 {
136  return true;
137 }
138 
139 
140 template<class Scalar>
142  const RCP<PreconditionerFactoryBase<Scalar> > &precFactory,
143  const std::string &precFactoryName
144  )
145 {
146  TEUCHOS_TEST_FOR_EXCEPT(!precFactory.get());
147  RCP<const Teuchos::ParameterList>
148  precFactoryValidPL = precFactory->getValidParameters();
149  const std::string _precFactoryName =
150  ( precFactoryName != ""
151  ? precFactoryName
152  : ( precFactoryValidPL.get() ? precFactoryValidPL->name() : "GENERIC PRECONDITIONER FACTORY" )
153  );
154  precFactory_ = precFactory;
155  precFactoryName_ = _precFactoryName;
156  updateThisValidParamList();
157 }
158 
159 
160 template<class Scalar>
161 RCP<PreconditionerFactoryBase<Scalar> >
163 {
164  return precFactory_;
165 }
166 
167 
168 template<class Scalar>
170  RCP<PreconditionerFactoryBase<Scalar> > *precFactory,
171  std::string *precFactoryName
172  )
173 {
174  if(precFactory) *precFactory = precFactory_;
175  if(precFactoryName) *precFactoryName = precFactoryName_;
176  precFactory_ = Teuchos::null;
177  precFactoryName_ = "";
178  updateThisValidParamList();
179 }
180 
181 
182 template<class Scalar>
184  const LinearOpSourceBase<Scalar> &fwdOpSrc
185  ) const
186 {
187  if(precFactory_.get())
188  return precFactory_->isCompatible(fwdOpSrc);
189  return true; // Without a preconditioner, we are compatible with all linear operators!
190 }
191 
192 
193 template<class Scalar>
194 RCP<LinearOpWithSolveBase<Scalar> >
196 {
197  return Teuchos::rcp(new BelosLinearOpWithSolve<Scalar>());
198 }
199 
200 
201 template<class Scalar>
203  const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
204  LinearOpWithSolveBase<Scalar> *Op,
205  const ESupportSolveUse supportSolveUse
206  ) const
207 {
208  using Teuchos::null;
209  initializeOpImpl(fwdOpSrc,null,null,false,Op,supportSolveUse);
210 }
211 
212 
213 template<class Scalar>
215  const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
216  LinearOpWithSolveBase<Scalar> *Op
217  ) const
218 {
219  using Teuchos::null;
220  initializeOpImpl(fwdOpSrc,null,null,true,Op,SUPPORT_SOLVE_UNSPECIFIED);
221 }
222 
223 
224 template<class Scalar>
226  const EPreconditionerInputType precOpType
227  ) const
228 {
229  if(precFactory_.get())
230  return true;
231  return (precOpType==PRECONDITIONER_INPUT_TYPE_AS_OPERATOR);
232 }
233 
234 
235 template<class Scalar>
237  const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
238  const RCP<const PreconditionerBase<Scalar> > &prec,
239  LinearOpWithSolveBase<Scalar> *Op,
240  const ESupportSolveUse supportSolveUse
241  ) const
242 {
243  using Teuchos::null;
244  initializeOpImpl(fwdOpSrc,null,prec,false,Op,supportSolveUse);
245 }
246 
247 
248 template<class Scalar>
250  const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
251  const RCP<const LinearOpSourceBase<Scalar> > &approxFwdOpSrc,
252  LinearOpWithSolveBase<Scalar> *Op,
253  const ESupportSolveUse supportSolveUse
254  ) const
255 {
256  using Teuchos::null;
257  initializeOpImpl(fwdOpSrc,approxFwdOpSrc,null,false,Op,supportSolveUse);
258 }
259 
260 
261 template<class Scalar>
263  LinearOpWithSolveBase<Scalar> *Op,
264  RCP<const LinearOpSourceBase<Scalar> > *fwdOpSrc,
265  RCP<const PreconditionerBase<Scalar> > *prec,
266  RCP<const LinearOpSourceBase<Scalar> > *approxFwdOpSrc,
267  ESupportSolveUse *supportSolveUse
268  ) const
269 {
270 #ifdef TEUCHOS_DEBUG
271  TEUCHOS_TEST_FOR_EXCEPT(Op==NULL);
272 #endif
274  &belosOp = Teuchos::dyn_cast<BelosLinearOpWithSolve<Scalar> >(*Op);
275  RCP<const LinearOpSourceBase<Scalar> >
276  _fwdOpSrc = belosOp.extract_fwdOpSrc();
277  RCP<const PreconditionerBase<Scalar> >
278  _prec = ( belosOp.isExternalPrec() ? belosOp.extract_prec() : Teuchos::null );
279  // Note: above we only extract the preconditioner if it was passed in
280  // externally. Otherwise, we need to hold on to it so that we can reuse it
281  // in the next initialization.
282  RCP<const LinearOpSourceBase<Scalar> >
283  _approxFwdOpSrc = belosOp.extract_approxFwdOpSrc();
284  ESupportSolveUse
285  _supportSolveUse = belosOp.supportSolveUse();
286  if(fwdOpSrc) *fwdOpSrc = _fwdOpSrc;
287  if(prec) *prec = _prec;
288  if(approxFwdOpSrc) *approxFwdOpSrc = _approxFwdOpSrc;
289  if(supportSolveUse) *supportSolveUse = _supportSolveUse;
290 }
291 
292 
293 // Overridden from ParameterListAcceptor
294 
295 
296 template<class Scalar>
298  RCP<Teuchos::ParameterList> const& paramList
299  )
300 {
301  TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
302  paramList->validateParametersAndSetDefaults(*this->getValidParameters(), 1);
303  paramList_ = paramList;
304  solverType_ =
305  Teuchos::getIntegralValue<EBelosSolverType>(*paramList_, SolverType_name);
306  convergenceTestFrequency_ =
307  Teuchos::getParameter<int>(*paramList_, ConvergenceTestFrequency_name);
308  Teuchos::readVerboseObjectSublist(&*paramList_,this);
309 }
310 
311 
312 template<class Scalar>
313 RCP<Teuchos::ParameterList>
315 {
316  return paramList_;
317 }
318 
319 
320 template<class Scalar>
321 RCP<Teuchos::ParameterList>
323 {
324  RCP<Teuchos::ParameterList> _paramList = paramList_;
325  paramList_ = Teuchos::null;
326  return _paramList;
327 }
328 
329 
330 template<class Scalar>
331 RCP<const Teuchos::ParameterList>
333 {
334  return paramList_;
335 }
336 
337 
338 template<class Scalar>
339 RCP<const Teuchos::ParameterList>
341 {
342  return thisValidParamList_;
343 }
344 
345 
346 // Public functions overridden from Teuchos::Describable
347 
348 
349 template<class Scalar>
351 {
352  std::ostringstream oss;
353  oss << "Thyra::BelosLinearOpWithSolveFactory";
354  //oss << "{";
355  // ToDo: Fill this in some!
356  //oss << "}";
357  return oss.str();
358 }
359 
360 
361 // private
362 
363 
364 template<class Scalar>
365 RCP<const Teuchos::ParameterList>
367 {
368  using Teuchos::as;
369  using Teuchos::tuple;
370  using Teuchos::setStringToIntegralParameter;
371 Teuchos::ValidatorXMLConverterDB::addConverter(
372  Teuchos::DummyObjectGetter<
373  Teuchos::StringToIntegralParameterEntryValidator<EBelosSolverType>
374  >::getDummyObject(),
375  Teuchos::DummyObjectGetter<Teuchos::StringToIntegralValidatorXMLConverter<
376  EBelosSolverType> >::getDummyObject());
377 
378  typedef MultiVectorBase<Scalar> MV_t;
379  typedef LinearOpBase<Scalar> LO_t;
380  static RCP<Teuchos::ParameterList> validParamList;
381  if(validParamList.get()==NULL) {
382  validParamList = Teuchos::rcp(new Teuchos::ParameterList("BelosLinearOpWithSolveFactory"));
383  setStringToIntegralParameter<EBelosSolverType>(
384  SolverType_name, SolverType_default,
385  "Type of linear solver algorithm to use.",
386  tuple<std::string>(
387  "Block GMRES",
388  "Pseudo Block GMRES",
389  "Block CG",
390  "Pseudo Block CG",
391  "Pseudo Block Stochastic CG",
392  "GCRODR",
393  "RCG",
394  "MINRES",
395  "TFQMR"
396  ),
397  tuple<std::string>(
398  "Block GMRES solver for nonsymmetric linear systems. It can also solve "
399  "single right-hand side systems, and can also perform Flexible GMRES "
400  "(where the preconditioner may change at every iteration, for example "
401  "for inner-outer iterations), by setting options in the \"Block GMRES\" "
402  "sublist.",
403 
404  "GMRES solver for nonsymmetric linear systems, that performs single "
405  "right-hand side solves on multiple right-hand sides at once. It "
406  "exploits operator multivector multiplication in order to amortize "
407  "global communication costs. Individual linear systems are deflated "
408  "out as they are solved.",
409 
410  "Block CG solver for symmetric (Hermitian in complex arithmetic) "
411  "positive definite linear systems. It can also solve single "
412  "right-hand-side systems.",
413 
414  "CG solver that performs single right-hand side CG on multiple right-hand "
415  "sides at once. It exploits operator multivector multiplication in order "
416  "to amortize global communication costs. Individual linear systems are "
417  "deflated out as they are solved.",
418 
419  "stochastic CG solver that performs single right-hand side CG on multiple right-hand "
420  "sides at once. It exploits operator multivector multiplication in order "
421  "to amortize global communication costs. Individual linear systems are "
422  "deflated out as they are solved. [EXPERIMENTAL]",
423 
424  "Variant of GMRES that performs subspace recycling to accelerate "
425  "convergence for sequences of solves with related linear systems. "
426  "Individual linear systems are deflated out as they are solved. "
427  "The current implementation only supports real-valued Scalar types.",
428 
429  "CG solver for symmetric (Hermitian in complex arithmetic) positive "
430  "definite linear systems, that performs subspace recycling to "
431  "accelerate convergence for sequences of related linear systems.",
432 
433  "MINRES solver for symmetric indefinite linear systems. It performs "
434  "single-right-hand-side solves on multiple right-hand sides sequentially.",
435 
436  "TFQMR (Transpose-Free QMR) solver for nonsymmetric linear systems."
437  ),
438  tuple<EBelosSolverType>(
439  SOLVER_TYPE_BLOCK_GMRES,
440  SOLVER_TYPE_PSEUDO_BLOCK_GMRES,
441  SOLVER_TYPE_BLOCK_CG,
442  SOLVER_TYPE_PSEUDO_BLOCK_CG,
443  SOLVER_TYPE_PSEUDO_BLOCK_STOCHASTIC_CG,
444  SOLVER_TYPE_GCRODR,
445  SOLVER_TYPE_RCG,
446  SOLVER_TYPE_MINRES,
447  SOLVER_TYPE_TFQMR
448  ),
449  &*validParamList
450  );
451  validParamList->set(ConvergenceTestFrequency_name, as<int>(1),
452  "Number of linear solver iterations to skip between applying"
453  " user-defined convergence test.");
454  validParamList->set(
455  LeftPreconditionerIfUnspecified_name, false,
456  "If the preconditioner does not specify if it is left or right, and this\n"
457  "option is set to true, put the preconditioner on the left side.\n"
458  "Historically, preconditioning is on the right. Some solvers may not\n"
459  "support left preconditioning.");
460  Teuchos::ParameterList
461  &solverTypesSL = validParamList->sublist(SolverTypes_name);
462  {
463  Belos::BlockGmresSolMgr<Scalar,MV_t,LO_t> mgr;
464  solverTypesSL.sublist(BlockGMRES_name).setParameters(
465  *mgr.getValidParameters()
466  );
467  }
468  {
469  Belos::PseudoBlockGmresSolMgr<Scalar,MV_t,LO_t> mgr;
470  solverTypesSL.sublist(PseudoBlockGMRES_name).setParameters(
471  *mgr.getValidParameters()
472  );
473  }
474  {
475  Belos::BlockCGSolMgr<Scalar,MV_t,LO_t> mgr;
476  solverTypesSL.sublist(BlockCG_name).setParameters(
477  *mgr.getValidParameters()
478  );
479  }
480  {
481  Belos::PseudoBlockCGSolMgr<Scalar,MV_t,LO_t> mgr;
482  solverTypesSL.sublist(PseudoBlockCG_name).setParameters(
483  *mgr.getValidParameters()
484  );
485  }
486  {
487  Belos::PseudoBlockStochasticCGSolMgr<Scalar,MV_t,LO_t> mgr;
488  solverTypesSL.sublist(PseudoBlockStochasticCG_name).setParameters(
489  *mgr.getValidParameters()
490  );
491  }
492  {
493  Belos::GCRODRSolMgr<Scalar,MV_t,LO_t> mgr;
494  solverTypesSL.sublist(GCRODR_name).setParameters(
495  *mgr.getValidParameters()
496  );
497  }
498  {
499  Belos::RCGSolMgr<Scalar,MV_t,LO_t> mgr;
500  solverTypesSL.sublist(RCG_name).setParameters(
501  *mgr.getValidParameters()
502  );
503  }
504  {
505  Belos::MinresSolMgr<Scalar,MV_t,LO_t> mgr;
506  solverTypesSL.sublist(MINRES_name).setParameters(
507  *mgr.getValidParameters()
508  );
509  }
510  {
511  Belos::TFQMRSolMgr<Scalar,MV_t,LO_t> mgr;
512  solverTypesSL.sublist(TFQMR_name).setParameters(
513  *mgr.getValidParameters()
514  );
515  }
516  }
517  return validParamList;
518 }
519 
520 
521 template<class Scalar>
522 void BelosLinearOpWithSolveFactory<Scalar>::updateThisValidParamList()
523 {
524  thisValidParamList_ = Teuchos::rcp(
525  new Teuchos::ParameterList(*generateAndGetValidParameters())
526  );
527  Teuchos::setupVerboseObjectSublist(&*thisValidParamList_);
528 }
529 
530 
531 template<class Scalar>
532 void BelosLinearOpWithSolveFactory<Scalar>::initializeOpImpl(
533  const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
534  const RCP<const LinearOpSourceBase<Scalar> > &approxFwdOpSrc,
535  const RCP<const PreconditionerBase<Scalar> > &prec_in,
536  const bool reusePrec,
537  LinearOpWithSolveBase<Scalar> *Op,
538  const ESupportSolveUse supportSolveUse
539  ) const
540 {
541 
542  using Teuchos::rcp;
543  using Teuchos::set_extra_data;
544  typedef Teuchos::ScalarTraits<Scalar> ST;
545  typedef MultiVectorBase<Scalar> MV_t;
546  typedef LinearOpBase<Scalar> LO_t;
547 
548  const RCP<Teuchos::FancyOStream> out = this->getOStream();
549  const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
550  Teuchos::OSTab tab(out);
551  if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW))
552  *out << "\nEntering Thyra::BelosLinearOpWithSolveFactory<"<<ST::name()<<">::initializeOpImpl(...) ...\n";
553 
554  // These lines are changing the verbosity of the preconditioner, which has its own verbose object list,
555  // so I am commenting these out, as it is not the job of the linear solver to dictate preconditioner verbosity.
556  //typedef Teuchos::VerboseObjectTempState<PreconditionerFactoryBase<Scalar> > VOTSPF;
557  //VOTSPF precFactoryOutputTempState(precFactory_,out,verbLevel);
558 
559  TEUCHOS_TEST_FOR_EXCEPT(Op==NULL);
560  TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
561  TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc->getOp().get()==NULL);
562  RCP<const LinearOpBase<Scalar> >
563  fwdOp = fwdOpSrc->getOp(),
564  approxFwdOp = ( approxFwdOpSrc.get() ? approxFwdOpSrc->getOp() : Teuchos::null );
565 
566  //
567  // Get the BelosLinearOpWithSolve interface
568  //
569 
570  BelosLinearOpWithSolve<Scalar>
571  *belosOp = &Teuchos::dyn_cast<BelosLinearOpWithSolve<Scalar> >(*Op);
572 
573  //
574  // Get/Create the preconditioner
575  //
576 
577  RCP<PreconditionerBase<Scalar> > myPrec = Teuchos::null;
578  RCP<const PreconditionerBase<Scalar> > prec = Teuchos::null;
579  if(prec_in.get()) {
580  // Use an externally defined preconditioner
581  prec = prec_in;
582  }
583  else {
584  // Try and generate a preconditioner on our own
585  if(precFactory_.get()) {
586  myPrec =
587  ( !belosOp->isExternalPrec()
588  ? Teuchos::rcp_const_cast<PreconditionerBase<Scalar> >(belosOp->extract_prec())
589  : Teuchos::null
590  );
591  bool hasExistingPrec = false;
592  if(myPrec.get()) {
593  hasExistingPrec = true;
594  // ToDo: Get the forward operator and validate that it is the same
595  // operator that is used here!
596  }
597  else {
598  hasExistingPrec = false;
599  myPrec = precFactory_->createPrec();
600  }
601  if( hasExistingPrec && reusePrec ) {
602  // Just reuse the existing preconditioner again!
603  }
604  else {
605  // Update the preconditioner
606  if(approxFwdOp.get())
607  precFactory_->initializePrec(approxFwdOpSrc,&*myPrec);
608  else
609  precFactory_->initializePrec(fwdOpSrc,&*myPrec);
610  }
611  prec = myPrec;
612  }
613  }
614 
615  //
616  // Uninitialize the current solver object
617  //
618 
619  bool oldIsExternalPrec = false;
620  RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > oldLP = Teuchos::null;
621  RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > oldIterSolver = Teuchos::null;
622  RCP<const LinearOpSourceBase<Scalar> > oldFwdOpSrc = Teuchos::null;
623  RCP<const LinearOpSourceBase<Scalar> > oldApproxFwdOpSrc = Teuchos::null;
624  ESupportSolveUse oldSupportSolveUse = SUPPORT_SOLVE_UNSPECIFIED;
625 
626  belosOp->uninitialize( &oldLP, NULL, &oldIterSolver, &oldFwdOpSrc,
627  NULL, &oldIsExternalPrec, &oldApproxFwdOpSrc, &oldSupportSolveUse );
628 
629  //
630  // Create the Belos linear problem
631  // NOTE: If one exists already, reuse it.
632  //
633 
634  typedef Belos::LinearProblem<Scalar,MV_t,LO_t> LP_t;
635  RCP<LP_t> lp;
636  if (oldLP != Teuchos::null) {
637  lp = oldLP;
638  }
639  else {
640  lp = rcp(new LP_t());
641  }
642 
643  //
644  // Set the operator
645  //
646 
647  lp->setOperator(fwdOp);
648 
649  //
650  // Set the preconditioner
651  //
652 
653  if(prec.get()) {
654  RCP<const LinearOpBase<Scalar> > unspecified = prec->getUnspecifiedPrecOp();
655  RCP<const LinearOpBase<Scalar> > left = prec->getLeftPrecOp();
656  RCP<const LinearOpBase<Scalar> > right = prec->getRightPrecOp();
657  TEUCHOS_TEST_FOR_EXCEPTION(
658  !( left.get() || right.get() || unspecified.get() ), std::logic_error
659  ,"Error, at least one preconditoner linear operator objects must be set!"
660  );
661  if(unspecified.get()) {
662  if (paramList_->get<bool>(LeftPreconditionerIfUnspecified_name, false))
663  lp->setLeftPrec(unspecified);
664  else
665  lp->setRightPrec(unspecified);
666  }
667  else {
668  // Set a left, right or split preconditioner
669  TEUCHOS_TEST_FOR_EXCEPTION(
670  left.get(),std::logic_error
671  ,"Error, we can not currently handle a left preconditioner!"
672  );
673  lp->setRightPrec(right);
674  }
675  }
676  if(myPrec.get()) {
677  set_extra_data<RCP<PreconditionerBase<Scalar> > >(myPrec,"Belos::InternalPrec",
678  Teuchos::inOutArg(lp), Teuchos::POST_DESTROY, false);
679  }
680  else if(prec.get()) {
681  set_extra_data<RCP<const PreconditionerBase<Scalar> > >(prec,"Belos::ExternalPrec",
682  Teuchos::inOutArg(lp), Teuchos::POST_DESTROY, false);
683  }
684 
685  //
686  // Generate the parameter list.
687  //
688 
689  typedef Belos::SolverManager<Scalar,MV_t,LO_t> IterativeSolver_t;
690  RCP<IterativeSolver_t> iterativeSolver = Teuchos::null;
691  RCP<Teuchos::ParameterList> solverPL = Teuchos::rcp( new Teuchos::ParameterList() );
692 
693  switch(solverType_) {
694  case SOLVER_TYPE_BLOCK_GMRES:
695  {
696  // Set the PL
697  if(paramList_.get()) {
698  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
699  Teuchos::ParameterList &gmresPL = solverTypesPL.sublist(BlockGMRES_name);
700  solverPL = Teuchos::rcp( &gmresPL, false );
701  }
702  // Create the solver
703  if (oldIterSolver != Teuchos::null) {
704  iterativeSolver = oldIterSolver;
705  iterativeSolver->setProblem( lp );
706  iterativeSolver->setParameters( solverPL );
707  }
708  else {
709  iterativeSolver = rcp(new Belos::BlockGmresSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
710  }
711  break;
712  }
713  case SOLVER_TYPE_PSEUDO_BLOCK_GMRES:
714  {
715  // Set the PL
716  if(paramList_.get()) {
717  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
718  Teuchos::ParameterList &pbgmresPL = solverTypesPL.sublist(PseudoBlockGMRES_name);
719  solverPL = Teuchos::rcp( &pbgmresPL, false );
720  }
721  //
722  // Create the solver
723  //
724  if (oldIterSolver != Teuchos::null) {
725  iterativeSolver = oldIterSolver;
726  iterativeSolver->setProblem( lp );
727  iterativeSolver->setParameters( solverPL );
728  }
729  else {
730  iterativeSolver = rcp(new Belos::PseudoBlockGmresSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
731  }
732  break;
733  }
734  case SOLVER_TYPE_BLOCK_CG:
735  {
736  // Set the PL
737  if(paramList_.get()) {
738  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
739  Teuchos::ParameterList &cgPL = solverTypesPL.sublist(BlockCG_name);
740  solverPL = Teuchos::rcp( &cgPL, false );
741  }
742  // Create the solver
743  if (oldIterSolver != Teuchos::null) {
744  iterativeSolver = oldIterSolver;
745  iterativeSolver->setProblem( lp );
746  iterativeSolver->setParameters( solverPL );
747  }
748  else {
749  iterativeSolver = rcp(new Belos::BlockCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
750  }
751  break;
752  }
753  case SOLVER_TYPE_PSEUDO_BLOCK_CG:
754  {
755  // Set the PL
756  if(paramList_.get()) {
757  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
758  Teuchos::ParameterList &pbcgPL = solverTypesPL.sublist(PseudoBlockCG_name);
759  solverPL = Teuchos::rcp( &pbcgPL, false );
760  }
761  //
762  // Create the solver
763  //
764  if (oldIterSolver != Teuchos::null) {
765  iterativeSolver = oldIterSolver;
766  iterativeSolver->setProblem( lp );
767  iterativeSolver->setParameters( solverPL );
768  }
769  else {
770  iterativeSolver = rcp(new Belos::PseudoBlockCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
771  }
772  break;
773  }
774  case SOLVER_TYPE_PSEUDO_BLOCK_STOCHASTIC_CG:
775  {
776  // Set the PL
777  if(paramList_.get()) {
778  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
779  Teuchos::ParameterList &pbcgPL = solverTypesPL.sublist(PseudoBlockStochasticCG_name);
780  solverPL = Teuchos::rcp( &pbcgPL, false );
781  }
782  //
783  // Create the solver
784  //
785  if (oldIterSolver != Teuchos::null) {
786  iterativeSolver = oldIterSolver;
787  iterativeSolver->setProblem( lp );
788  iterativeSolver->setParameters( solverPL );
789  }
790  else {
791  iterativeSolver = rcp(new Belos::PseudoBlockStochasticCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
792  }
793  break;
794  }
795  case SOLVER_TYPE_GCRODR:
796  {
797  // Set the PL
798  if(paramList_.get()) {
799  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
800  Teuchos::ParameterList &gcrodrPL = solverTypesPL.sublist(GCRODR_name);
801  solverPL = Teuchos::rcp( &gcrodrPL, false );
802  }
803  // Create the solver
804  if (oldIterSolver != Teuchos::null) {
805  iterativeSolver = oldIterSolver;
806  iterativeSolver->setProblem( lp );
807  iterativeSolver->setParameters( solverPL );
808  }
809  else {
810  iterativeSolver = rcp(new Belos::GCRODRSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
811  }
812  break;
813  }
814  case SOLVER_TYPE_RCG:
815  {
816  // Set the PL
817  if(paramList_.get()) {
818  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
819  Teuchos::ParameterList &rcgPL = solverTypesPL.sublist(RCG_name);
820  solverPL = Teuchos::rcp( &rcgPL, false );
821  }
822  // Create the solver
823  if (oldIterSolver != Teuchos::null) {
824  iterativeSolver = oldIterSolver;
825  iterativeSolver->setProblem( lp );
826  iterativeSolver->setParameters( solverPL );
827  }
828  else {
829  iterativeSolver = rcp(new Belos::RCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
830  }
831  break;
832  }
833  case SOLVER_TYPE_MINRES:
834  {
835  // Set the PL
836  if(paramList_.get()) {
837  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
838  Teuchos::ParameterList &minresPL = solverTypesPL.sublist(MINRES_name);
839  solverPL = Teuchos::rcp( &minresPL, false );
840  }
841  // Create the solver
842  if (oldIterSolver != Teuchos::null) {
843  iterativeSolver = oldIterSolver;
844  iterativeSolver->setProblem( lp );
845  iterativeSolver->setParameters( solverPL );
846  }
847  else {
848  iterativeSolver = rcp(new Belos::MinresSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
849  }
850  break;
851  }
852  case SOLVER_TYPE_TFQMR:
853  {
854  // Set the PL
855  if(paramList_.get()) {
856  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
857  Teuchos::ParameterList &minresPL = solverTypesPL.sublist(TFQMR_name);
858  solverPL = Teuchos::rcp( &minresPL, false );
859  }
860  // Create the solver
861  if (oldIterSolver != Teuchos::null) {
862  iterativeSolver = oldIterSolver;
863  iterativeSolver->setProblem( lp );
864  iterativeSolver->setParameters( solverPL );
865  }
866  else {
867  iterativeSolver = rcp(new Belos::TFQMRSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
868  }
869  break;
870  }
871 
872  default:
873  {
874  TEUCHOS_TEST_FOR_EXCEPT(true);
875  }
876  }
877 
878  //
879  // Initialize the LOWS object
880  //
881 
882  belosOp->initialize(
883  lp, solverPL, iterativeSolver,
884  fwdOpSrc, prec, myPrec.get()==NULL, approxFwdOpSrc,
885  supportSolveUse, convergenceTestFrequency_
886  );
887  belosOp->setOStream(out);
888  belosOp->setVerbLevel(verbLevel);
889 #ifdef TEUCHOS_DEBUG
890  if(paramList_.get()) {
891  // Make sure we read the list correctly
892  paramList_->validateParameters(*this->getValidParameters(),1); // Validate 0th and 1st level deep
893  }
894 #endif
895  if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW))
896  *out << "\nLeaving Thyra::BelosLinearOpWithSolveFactory<"<<ST::name()<<">::initializeOpImpl(...) ...\n";
897 
898 }
899 
900 
901 } // namespace Thyra
902 
903 
904 #endif // THYRA_BELOS_LINEAR_OP_WITH_SOLVE_FACTORY_HPP
bool isCompatible(const LinearOpSourceBase< Scalar > &fwdOpSrc) const
RCP< const LinearOpSourceBase< Scalar > > extract_approxFwdOpSrc()
Concrete LinearOpWithSolveBase subclass in terms of Belos.
void initializeAndReuseOp(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, LinearOpWithSolveBase< Scalar > *Op) const
RCP< const PreconditionerBase< Scalar > > extract_prec()
bool supportsPreconditionerInputType(const EPreconditionerInputType precOpType) const
void unsetPreconditionerFactory(Teuchos::RCP< PreconditionerFactoryBase< Scalar > > *precFactory, std::string *precFactoryName)
Thyra specializations of MultiVecTraits and OperatorTraits.
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
void initializePreconditionedOp(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, const Teuchos::RCP< const PreconditionerBase< Scalar > > &prec, LinearOpWithSolveBase< Scalar > *Op, const ESupportSolveUse supportSolveUse) const
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
void setPreconditionerFactory(const Teuchos::RCP< PreconditionerFactoryBase< Scalar > > &precFactory, const std::string &precFactoryName)
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
Teuchos::RCP< PreconditionerFactoryBase< Scalar > > getPreconditionerFactory() const
LinearOpWithSolveFactoryBase subclass implemented in terms of Belos.
RCP< const LinearOpSourceBase< Scalar > > extract_fwdOpSrc()
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
BelosLinearOpWithSolveFactory()
Construct without preconditioner factory.
Teuchos::RCP< LinearOpWithSolveBase< Scalar > > createOp() const
void uninitializeOp(LinearOpWithSolveBase< Scalar > *Op, Teuchos::RCP< const LinearOpSourceBase< Scalar > > *fwdOpSrc, Teuchos::RCP< const PreconditionerBase< Scalar > > *prec, Teuchos::RCP< const LinearOpSourceBase< Scalar > > *approxFwdOpSrc, ESupportSolveUse *supportSolveUse) const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
void initializeApproxPreconditionedOp(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &approxFwdOpSrc, LinearOpWithSolveBase< Scalar > *Op, const ESupportSolveUse supportSolveUse) const
void initializeOp(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, LinearOpWithSolveBase< Scalar > *Op, const ESupportSolveUse supportSolveUse) const

Generated for Stratimikos by doxygen 1.8.14