Stratimikos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_BelosLinearOpWithSolveFactory_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_FACTORY_HPP
46 #define THYRA_BELOS_LINEAR_OP_WITH_SOLVE_FACTORY_HPP
47 
48 
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"
67 #include "Teuchos_dyn_cast.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 {
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());
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>
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>
196 {
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
276  _fwdOpSrc = belosOp.extract_fwdOpSrc();
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.
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>
315 {
316  return paramList_;
317 }
318 
319 
320 template<class Scalar>
323 {
324  RCP<Teuchos::ParameterList> _paramList = paramList_;
325  paramList_ = Teuchos::null;
326  return _paramList;
327 }
328 
329 
330 template<class Scalar>
333 {
334  return paramList_;
335 }
336 
337 
338 template<class Scalar>
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>
367 {
368  using Teuchos::as;
369  using Teuchos::tuple;
370  using Teuchos::setStringToIntegralParameter;
374  >::getDummyObject(),
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.");
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>
523 {
524  thisValidParamList_ = Teuchos::rcp(
525  new Teuchos::ParameterList(*generateAndGetValidParameters())
526  );
527  Teuchos::setupVerboseObjectSublist(&*thisValidParamList_);
528 }
529 
530 
531 template<class Scalar>
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;
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);
563  fwdOp = fwdOpSrc->getOp(),
564  approxFwdOp = ( approxFwdOpSrc.get() ? approxFwdOpSrc->getOp() : Teuchos::null );
565 
566  //
567  // Get the BelosLinearOpWithSolve interface
568  //
569 
572 
573  //
574  // Get/Create the preconditioner
575  //
576 
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;
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();
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 (nonnull(unspecified)) {
662  if (paramList_->get<bool>(LeftPreconditionerIfUnspecified_name, false))
663  lp->setLeftPrec(unspecified);
664  else
665  lp->setRightPrec(unspecified);
666  }
667  else if (nonnull(left)) {
668  lp->setLeftPrec(left);
669  }
670  else if (nonnull(right)) {
671  lp->setRightPrec(right);
672  }
673  else {
674  // Set a left, right or split preconditioner
676  nonnull(left) && nonnull(right),std::logic_error
677  ,"Error, we can not currently handle split preconditioners!"
678  );
679  }
680  }
681  if(myPrec.get()) {
682  set_extra_data<RCP<PreconditionerBase<Scalar> > >(myPrec,"Belos::InternalPrec",
683  Teuchos::inOutArg(lp), Teuchos::POST_DESTROY, false);
684  }
685  else if(prec.get()) {
686  set_extra_data<RCP<const PreconditionerBase<Scalar> > >(prec,"Belos::ExternalPrec",
687  Teuchos::inOutArg(lp), Teuchos::POST_DESTROY, false);
688  }
689 
690  //
691  // Generate the parameter list.
692  //
693 
694  typedef Belos::SolverManager<Scalar,MV_t,LO_t> IterativeSolver_t;
695  RCP<IterativeSolver_t> iterativeSolver = Teuchos::null;
697 
698  switch(solverType_) {
699  case SOLVER_TYPE_BLOCK_GMRES:
700  {
701  // Set the PL
702  if(paramList_.get()) {
703  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
704  Teuchos::ParameterList &gmresPL = solverTypesPL.sublist(BlockGMRES_name);
705  solverPL = Teuchos::rcp( &gmresPL, false );
706  }
707  // Create the solver
708  if (oldIterSolver != Teuchos::null) {
709  iterativeSolver = oldIterSolver;
710  iterativeSolver->setProblem( lp );
711  iterativeSolver->setParameters( solverPL );
712  }
713  else {
714  iterativeSolver = rcp(new Belos::BlockGmresSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
715  }
716  break;
717  }
718  case SOLVER_TYPE_PSEUDO_BLOCK_GMRES:
719  {
720  // Set the PL
721  if(paramList_.get()) {
722  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
723  Teuchos::ParameterList &pbgmresPL = solverTypesPL.sublist(PseudoBlockGMRES_name);
724  solverPL = Teuchos::rcp( &pbgmresPL, false );
725  }
726  //
727  // Create the solver
728  //
729  if (oldIterSolver != Teuchos::null) {
730  iterativeSolver = oldIterSolver;
731  iterativeSolver->setProblem( lp );
732  iterativeSolver->setParameters( solverPL );
733  }
734  else {
735  iterativeSolver = rcp(new Belos::PseudoBlockGmresSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
736  }
737  break;
738  }
739  case SOLVER_TYPE_BLOCK_CG:
740  {
741  // Set the PL
742  if(paramList_.get()) {
743  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
744  Teuchos::ParameterList &cgPL = solverTypesPL.sublist(BlockCG_name);
745  solverPL = Teuchos::rcp( &cgPL, false );
746  }
747  // Create the solver
748  if (oldIterSolver != Teuchos::null) {
749  iterativeSolver = oldIterSolver;
750  iterativeSolver->setProblem( lp );
751  iterativeSolver->setParameters( solverPL );
752  }
753  else {
754  iterativeSolver = rcp(new Belos::BlockCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
755  }
756  break;
757  }
758  case SOLVER_TYPE_PSEUDO_BLOCK_CG:
759  {
760  // Set the PL
761  if(paramList_.get()) {
762  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
763  Teuchos::ParameterList &pbcgPL = solverTypesPL.sublist(PseudoBlockCG_name);
764  solverPL = Teuchos::rcp( &pbcgPL, false );
765  }
766  //
767  // Create the solver
768  //
769  if (oldIterSolver != Teuchos::null) {
770  iterativeSolver = oldIterSolver;
771  iterativeSolver->setProblem( lp );
772  iterativeSolver->setParameters( solverPL );
773  }
774  else {
775  iterativeSolver = rcp(new Belos::PseudoBlockCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
776  }
777  break;
778  }
779  case SOLVER_TYPE_PSEUDO_BLOCK_STOCHASTIC_CG:
780  {
781  // Set the PL
782  if(paramList_.get()) {
783  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
784  Teuchos::ParameterList &pbcgPL = solverTypesPL.sublist(PseudoBlockStochasticCG_name);
785  solverPL = Teuchos::rcp( &pbcgPL, false );
786  }
787  //
788  // Create the solver
789  //
790  if (oldIterSolver != Teuchos::null) {
791  iterativeSolver = oldIterSolver;
792  iterativeSolver->setProblem( lp );
793  iterativeSolver->setParameters( solverPL );
794  }
795  else {
796  iterativeSolver = rcp(new Belos::PseudoBlockStochasticCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
797  }
798  break;
799  }
800  case SOLVER_TYPE_GCRODR:
801  {
802  // Set the PL
803  if(paramList_.get()) {
804  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
805  Teuchos::ParameterList &gcrodrPL = solverTypesPL.sublist(GCRODR_name);
806  solverPL = Teuchos::rcp( &gcrodrPL, false );
807  }
808  // Create the solver
809  if (oldIterSolver != Teuchos::null) {
810  iterativeSolver = oldIterSolver;
811  iterativeSolver->setProblem( lp );
812  iterativeSolver->setParameters( solverPL );
813  }
814  else {
815  iterativeSolver = rcp(new Belos::GCRODRSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
816  }
817  break;
818  }
819  case SOLVER_TYPE_RCG:
820  {
821  // Set the PL
822  if(paramList_.get()) {
823  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
824  Teuchos::ParameterList &rcgPL = solverTypesPL.sublist(RCG_name);
825  solverPL = Teuchos::rcp( &rcgPL, false );
826  }
827  // Create the solver
828  if (oldIterSolver != Teuchos::null) {
829  iterativeSolver = oldIterSolver;
830  iterativeSolver->setProblem( lp );
831  iterativeSolver->setParameters( solverPL );
832  }
833  else {
834  iterativeSolver = rcp(new Belos::RCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
835  }
836  break;
837  }
838  case SOLVER_TYPE_MINRES:
839  {
840  // Set the PL
841  if(paramList_.get()) {
842  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
843  Teuchos::ParameterList &minresPL = solverTypesPL.sublist(MINRES_name);
844  solverPL = Teuchos::rcp( &minresPL, false );
845  }
846  // Create the solver
847  if (oldIterSolver != Teuchos::null) {
848  iterativeSolver = oldIterSolver;
849  iterativeSolver->setProblem( lp );
850  iterativeSolver->setParameters( solverPL );
851  }
852  else {
853  iterativeSolver = rcp(new Belos::MinresSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
854  }
855  break;
856  }
857  case SOLVER_TYPE_TFQMR:
858  {
859  // Set the PL
860  if(paramList_.get()) {
861  Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
862  Teuchos::ParameterList &minresPL = solverTypesPL.sublist(TFQMR_name);
863  solverPL = Teuchos::rcp( &minresPL, false );
864  }
865  // Create the solver
866  if (oldIterSolver != Teuchos::null) {
867  iterativeSolver = oldIterSolver;
868  iterativeSolver->setProblem( lp );
869  iterativeSolver->setParameters( solverPL );
870  }
871  else {
872  iterativeSolver = rcp(new Belos::TFQMRSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
873  }
874  break;
875  }
876 
877  default:
878  {
880  }
881  }
882 
883  //
884  // Initialize the LOWS object
885  //
886 
887  belosOp->initialize(
888  lp, solverPL, iterativeSolver,
889  fwdOpSrc, prec, myPrec.get()==NULL, approxFwdOpSrc,
890  supportSolveUse, convergenceTestFrequency_
891  );
892  belosOp->setOStream(out);
893  belosOp->setVerbLevel(verbLevel);
894 #ifdef TEUCHOS_DEBUG
895  if(paramList_.get()) {
896  // Make sure we read the list correctly
897  paramList_->validateParameters(*this->getValidParameters(),1); // Validate 0th and 1st level deep
898  }
899 #endif
900  if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW))
901  *out << "\nLeaving Thyra::BelosLinearOpWithSolveFactory<"<<ST::name()<<">::initializeOpImpl(...) ...\n";
902 
903 }
904 
905 
906 } // namespace Thyra
907 
908 
909 #endif // THYRA_BELOS_LINEAR_OP_WITH_SOLVE_FACTORY_HPP
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() 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
RCP< const LinearOpSourceBase< Scalar > > extract_approxFwdOpSrc()
static void addConverter(RCP< const ParameterEntryValidator > validator, RCP< ValidatorXMLConverter > converterToAdd)
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.
Concrete LinearOpWithSolveBase subclass in terms of Belos.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
RCP< const PreconditionerBase< Scalar > > extract_prec()
double ST
void initializeAndReuseOp(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, LinearOpWithSolveBase< Scalar > *Op) const
RCP< ParameterList > sublist(const RCP< ParameterList > &paramList, const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
T * get() const
static Teuchos::RCP< const Teuchos::ParameterList > generateAndGetValidParameters()
void unsetPreconditionerFactory(Teuchos::RCP< PreconditionerFactoryBase< Scalar > > *precFactory, std::string *precFactoryName)
Thyra specializations of MultiVecTraits and OperatorTraits.
void setParameterList(Teuchos::RCP< Teuchos::ParameterList > const &paramList)
bool isCompatible(const LinearOpSourceBase< Scalar > &fwdOpSrc) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
T_To & dyn_cast(T_From &from)
Teuchos::RCP< Teuchos::ParameterList > getNonconstParameterList()
void initializeOp(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, LinearOpWithSolveBase< Scalar > *Op, const ESupportSolveUse supportSolveUse) const
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
void setPreconditionerFactory(const Teuchos::RCP< PreconditionerFactoryBase< Scalar > > &precFactory, const std::string &precFactoryName)
TypeTo as(const TypeFrom &t)
bool supportsPreconditionerInputType(const EPreconditionerInputType precOpType) const
Teuchos::RCP< const Teuchos::ParameterList > getParameterList() const
RCP< const LinearOpSourceBase< Scalar > > extract_fwdOpSrc()
bool nonnull(const boost::shared_ptr< T > &p)
Teuchos::RCP< PreconditionerFactoryBase< Scalar > > getPreconditionerFactory() const
void initializeApproxPreconditionedOp(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &approxFwdOpSrc, LinearOpWithSolveBase< Scalar > *Op, const ESupportSolveUse supportSolveUse) const
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< LinearOpWithSolveBase< Scalar > > createOp() const
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.
BelosLinearOpWithSolveFactory()
Construct without preconditioner factory.
void initializeOpImpl(const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &fwdOpSrc, const Teuchos::RCP< const LinearOpSourceBase< Scalar > > &approxFwdOpSrc, const Teuchos::RCP< const PreconditionerBase< Scalar > > &prec, const bool reusePrec, LinearOpWithSolveBase< Scalar > *Op, const ESupportSolveUse supportSolveUse) const
Teuchos::RCP< Teuchos::ParameterList > unsetParameterList()
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)