Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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#include "BelosBiCGStabSolMgr.hpp"
63#include "BelosFixedPointSolMgr.hpp"
64#include "BelosThyraAdapter.hpp"
65
67
68#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
69#include "Teuchos_StandardParameterEntryValidators.hpp"
70#include "Teuchos_ParameterList.hpp"
71#include "Teuchos_dyn_cast.hpp"
72#include "Teuchos_ValidatorXMLConverterDB.hpp"
73#include "Teuchos_StandardValidatorXMLConverters.hpp"
74
75
76namespace Thyra {
77
78
79// Parameter names for Parameter List
80
81template<class Scalar>
82const std::string BelosLinearOpWithSolveFactory<Scalar>::SolverType_name = "Solver Type";
83template<class Scalar>
84const std::string BelosLinearOpWithSolveFactory<Scalar>::SolverType_default = "Pseudo Block GMRES";
85template<class Scalar>
86const std::string BelosLinearOpWithSolveFactory<Scalar>::SolverTypes_name = "Solver Types";
87template<class Scalar>
88const std::string BelosLinearOpWithSolveFactory<Scalar>::BlockGMRES_name = "Block GMRES";
89template<class Scalar>
90const std::string BelosLinearOpWithSolveFactory<Scalar>::PseudoBlockGMRES_name = "Pseudo Block GMRES";
91template<class Scalar>
92const std::string BelosLinearOpWithSolveFactory<Scalar>::BlockCG_name = "Block CG";
93template<class Scalar>
94const std::string BelosLinearOpWithSolveFactory<Scalar>::PseudoBlockCG_name = "Pseudo Block CG";
95template<class Scalar>
96const std::string BelosLinearOpWithSolveFactory<Scalar>::PseudoBlockStochasticCG_name = "Pseudo Block Stochastic CG";
97template<class Scalar>
99template<class Scalar>
101template<class Scalar>
102const std::string BelosLinearOpWithSolveFactory<Scalar>::MINRES_name = "MINRES";
103template<class Scalar>
104const std::string BelosLinearOpWithSolveFactory<Scalar>::TFQMR_name = "TFQMR";
105template<class Scalar>
106const std::string BelosLinearOpWithSolveFactory<Scalar>::BiCGStab_name = "BiCGStab";
107template<class Scalar>
108const std::string BelosLinearOpWithSolveFactory<Scalar>::FixedPoint_name = "Fixed Point";
109template<class Scalar>
110const std::string BelosLinearOpWithSolveFactory<Scalar>::TpetraGmres_name = "TPETRA GMRES";
111template<class Scalar>
112const std::string BelosLinearOpWithSolveFactory<Scalar>::TpetraGmresPipeline_name = "TPETRA GMRES PIPELINE";
113template<class Scalar>
114const std::string BelosLinearOpWithSolveFactory<Scalar>::TpetraGmresSingleReduce_name = "TPETRA GMRES SINGLE REDUCE";
115template<class Scalar>
116const std::string BelosLinearOpWithSolveFactory<Scalar>::TpetraGmresSstep_name = "TPETRA GMRES S-STEP";
117
118template<class Scalar>
119const std::string BelosLinearOpWithSolveFactory<Scalar>::ConvergenceTestFrequency_name = "Convergence Test Frequency";
120
121namespace {
122const std::string LeftPreconditionerIfUnspecified_name = "Left Preconditioner If Unspecified";
123}
124
125// Constructors/initializers/accessors
126
127
128template<class Scalar>
135
136
137template<class Scalar>
139 const RCP<PreconditionerFactoryBase<Scalar> > &precFactory
140 )
142{
143 this->setPreconditionerFactory(precFactory, "");
144}
145
146
147// Overridden from LinearOpWithSolveFactoryBase
148
149
150template<class Scalar>
155
156
157template<class Scalar>
159 const RCP<PreconditionerFactoryBase<Scalar> > &precFactory,
160 const std::string &precFactoryName
161 )
162{
163 TEUCHOS_TEST_FOR_EXCEPT(!precFactory.get());
164 RCP<const Teuchos::ParameterList>
165 precFactoryValidPL = precFactory->getValidParameters();
166 const std::string _precFactoryName =
167 ( precFactoryName != ""
168 ? precFactoryName
169 : ( precFactoryValidPL.get() ? precFactoryValidPL->name() : "GENERIC PRECONDITIONER FACTORY" )
170 );
171 precFactory_ = precFactory;
172 precFactoryName_ = _precFactoryName;
173 updateThisValidParamList();
174}
175
176
177template<class Scalar>
178RCP<PreconditionerFactoryBase<Scalar> >
183
184
185template<class Scalar>
187 RCP<PreconditionerFactoryBase<Scalar> > *precFactory,
188 std::string *precFactoryName
189 )
190{
191 if(precFactory) *precFactory = precFactory_;
192 if(precFactoryName) *precFactoryName = precFactoryName_;
193 precFactory_ = Teuchos::null;
194 precFactoryName_ = "";
195 updateThisValidParamList();
196}
197
198
199template<class Scalar>
201 const LinearOpSourceBase<Scalar> &fwdOpSrc
202 ) const
203{
204 if(precFactory_.get())
205 return precFactory_->isCompatible(fwdOpSrc);
206 return true; // Without a preconditioner, we are compatible with all linear operators!
207}
208
209
210template<class Scalar>
211RCP<LinearOpWithSolveBase<Scalar> >
213{
214 return Teuchos::rcp(new BelosLinearOpWithSolve<Scalar>());
215}
216
217
218template<class Scalar>
220 const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
221 LinearOpWithSolveBase<Scalar> *Op,
222 const ESupportSolveUse supportSolveUse
223 ) const
224{
225 using Teuchos::null;
226 initializeOpImpl(fwdOpSrc,null,null,false,Op,supportSolveUse);
227}
228
229
230template<class Scalar>
232 const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
233 LinearOpWithSolveBase<Scalar> *Op
234 ) const
235{
236 using Teuchos::null;
237 initializeOpImpl(fwdOpSrc,null,null,true,Op,SUPPORT_SOLVE_UNSPECIFIED);
238}
239
240
241template<class Scalar>
243 const EPreconditionerInputType precOpType
244 ) const
245{
246 if(precFactory_.get())
247 return true;
248 return (precOpType==PRECONDITIONER_INPUT_TYPE_AS_OPERATOR);
249}
250
251
252template<class Scalar>
254 const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
255 const RCP<const PreconditionerBase<Scalar> > &prec,
256 LinearOpWithSolveBase<Scalar> *Op,
257 const ESupportSolveUse supportSolveUse
258 ) const
259{
260 using Teuchos::null;
261 initializeOpImpl(fwdOpSrc,null,prec,false,Op,supportSolveUse);
262}
263
264
265template<class Scalar>
267 const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
268 const RCP<const LinearOpSourceBase<Scalar> > &approxFwdOpSrc,
269 LinearOpWithSolveBase<Scalar> *Op,
270 const ESupportSolveUse supportSolveUse
271 ) const
272{
273 using Teuchos::null;
274 initializeOpImpl(fwdOpSrc,approxFwdOpSrc,null,false,Op,supportSolveUse);
275}
276
277
278template<class Scalar>
280 LinearOpWithSolveBase<Scalar> *Op,
281 RCP<const LinearOpSourceBase<Scalar> > *fwdOpSrc,
282 RCP<const PreconditionerBase<Scalar> > *prec,
283 RCP<const LinearOpSourceBase<Scalar> > *approxFwdOpSrc,
284 ESupportSolveUse *supportSolveUse
285 ) const
286{
287#ifdef TEUCHOS_DEBUG
288 TEUCHOS_TEST_FOR_EXCEPT(Op==NULL);
289#endif
290 BelosLinearOpWithSolve<Scalar>
291 &belosOp = Teuchos::dyn_cast<BelosLinearOpWithSolve<Scalar> >(*Op);
292 RCP<const LinearOpSourceBase<Scalar> >
293 _fwdOpSrc = belosOp.extract_fwdOpSrc();
294 RCP<const PreconditionerBase<Scalar> >
295 _prec = ( belosOp.isExternalPrec() ? belosOp.extract_prec() : Teuchos::null );
296 // Note: above we only extract the preconditioner if it was passed in
297 // externally. Otherwise, we need to hold on to it so that we can reuse it
298 // in the next initialization.
299 RCP<const LinearOpSourceBase<Scalar> >
300 _approxFwdOpSrc = belosOp.extract_approxFwdOpSrc();
301 ESupportSolveUse
302 _supportSolveUse = belosOp.supportSolveUse();
303 if(fwdOpSrc) *fwdOpSrc = _fwdOpSrc;
304 if(prec) *prec = _prec;
305 if(approxFwdOpSrc) *approxFwdOpSrc = _approxFwdOpSrc;
306 if(supportSolveUse) *supportSolveUse = _supportSolveUse;
307}
308
309
310// Overridden from ParameterListAcceptor
311
312
313template<class Scalar>
315 RCP<Teuchos::ParameterList> const& paramList
316 )
317{
318 TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
319 paramList->validateParametersAndSetDefaults(*this->getValidParameters(), 1);
320 paramList_ = paramList;
321 solverType_ =
322 Teuchos::getIntegralValue<EBelosSolverType>(*paramList_, SolverType_name);
323 convergenceTestFrequency_ =
324 Teuchos::getParameter<int>(*paramList_, ConvergenceTestFrequency_name);
325 Teuchos::readVerboseObjectSublist(&*paramList_,this);
326}
327
328
329template<class Scalar>
330RCP<Teuchos::ParameterList>
335
336
337template<class Scalar>
338RCP<Teuchos::ParameterList>
340{
341 RCP<Teuchos::ParameterList> _paramList = paramList_;
342 paramList_ = Teuchos::null;
343 return _paramList;
344}
345
346
347template<class Scalar>
348RCP<const Teuchos::ParameterList>
350{
351 return paramList_;
352}
353
354
355template<class Scalar>
356RCP<const Teuchos::ParameterList>
358{
359 return thisValidParamList_;
360}
361
362
363// Public functions overridden from Teuchos::Describable
364
365
366template<class Scalar>
368{
369 std::ostringstream oss;
370 oss << "Thyra::BelosLinearOpWithSolveFactory";
371 //oss << "{";
372 // ToDo: Fill this in some!
373 //oss << "}";
374 return oss.str();
375}
376
377
378// private
379
380
381template<class Scalar>
382RCP<const Teuchos::ParameterList>
384{
385 using Teuchos::as;
386 using Teuchos::tuple;
387 using Teuchos::setStringToIntegralParameter;
388Teuchos::ValidatorXMLConverterDB::addConverter(
389 Teuchos::DummyObjectGetter<
390 Teuchos::StringToIntegralParameterEntryValidator<EBelosSolverType>
391 >::getDummyObject(),
392 Teuchos::DummyObjectGetter<Teuchos::StringToIntegralValidatorXMLConverter<
393 EBelosSolverType> >::getDummyObject());
394
395 typedef MultiVectorBase<Scalar> MV_t;
396 typedef LinearOpBase<Scalar> LO_t;
397 static RCP<Teuchos::ParameterList> validParamList;
398 if(validParamList.get()==NULL) {
399 validParamList = Teuchos::rcp(new Teuchos::ParameterList("BelosLinearOpWithSolveFactory"));
400 setStringToIntegralParameter<EBelosSolverType>(
401 SolverType_name, SolverType_default,
402 "Type of linear solver algorithm to use.",
403 tuple<std::string>(
404 "Block GMRES",
405 "Pseudo Block GMRES",
406 "Block CG",
407 "Pseudo Block CG",
408 "Pseudo Block Stochastic CG",
409 "GCRODR",
410 "RCG",
411 "MINRES",
412 "TFQMR",
413 "BiCGStab",
414 "Fixed Point",
415 "TPETRA GMRES",
416 "TPETRA GMRES PIPELINE",
417 "TPETRA GMRES SINGLE REDUCE",
418 "TPETRA GMRES S-STEP"
419 ),
420 tuple<std::string>(
421 "Block GMRES solver for nonsymmetric linear systems. It can also solve "
422 "single right-hand side systems, and can also perform Flexible GMRES "
423 "(where the preconditioner may change at every iteration, for example "
424 "for inner-outer iterations), by setting options in the \"Block GMRES\" "
425 "sublist.",
426
427 "GMRES solver for nonsymmetric linear systems, that performs single "
428 "right-hand side solves on multiple right-hand sides at once. It "
429 "exploits operator multivector multiplication in order to amortize "
430 "global communication costs. Individual linear systems are deflated "
431 "out as they are solved.",
432
433 "Block CG solver for symmetric (Hermitian in complex arithmetic) "
434 "positive definite linear systems. It can also solve single "
435 "right-hand-side systems.",
436
437 "CG solver that performs single right-hand side CG on multiple right-hand "
438 "sides at once. It exploits operator multivector multiplication in order "
439 "to amortize global communication costs. Individual linear systems are "
440 "deflated out as they are solved.",
441
442 "stochastic CG solver that performs single right-hand side CG on multiple right-hand "
443 "sides at once. It exploits operator multivector multiplication in order "
444 "to amortize global communication costs. Individual linear systems are "
445 "deflated out as they are solved. [EXPERIMENTAL]",
446
447 "Variant of GMRES that performs subspace recycling to accelerate "
448 "convergence for sequences of solves with related linear systems. "
449 "Individual linear systems are deflated out as they are solved. "
450 "The current implementation only supports real-valued Scalar types.",
451
452 "CG solver for symmetric (Hermitian in complex arithmetic) positive "
453 "definite linear systems, that performs subspace recycling to "
454 "accelerate convergence for sequences of related linear systems.",
455
456 "MINRES solver for symmetric indefinite linear systems. It performs "
457 "single-right-hand-side solves on multiple right-hand sides sequentially.",
458
459 "TFQMR (Transpose-Free QMR) solver for nonsymmetric linear systems.",
460
461 "BiCGStab solver for nonsymmetric linear systems.",
462
463 "Fixed point iteration",
464
465 "Native Tpetra implementation of GMRES",
466
467 "Native Tpetra implementation of pipeline GMRES",
468
469 "Native Tpetra implementation of single-reduce GMRES",
470
471 "Native Tpetra implementation of s-step GMRES"
472 ),
473 tuple<EBelosSolverType>(
489 ),
490 &*validParamList
491 );
492 validParamList->set(ConvergenceTestFrequency_name, as<int>(1),
493 "Number of linear solver iterations to skip between applying"
494 " user-defined convergence test.");
495 validParamList->set(
496 LeftPreconditionerIfUnspecified_name, false,
497 "If the preconditioner does not specify if it is left or right, and this\n"
498 "option is set to true, put the preconditioner on the left side.\n"
499 "Historically, preconditioning is on the right. Some solvers may not\n"
500 "support left preconditioning.");
501 Teuchos::ParameterList
502 &solverTypesSL = validParamList->sublist(SolverTypes_name);
503
504 const bool lapackSupportsScalar = Belos::Details::LapackSupportsScalar<Scalar>::value;
505 const bool scalarIsReal = !Teuchos::ScalarTraits<Scalar>::isComplex;
506
507 {
508 Belos::BlockGmresSolMgr<Scalar,MV_t,LO_t> mgr;
509 solverTypesSL.sublist(BlockGMRES_name).setParameters(
510 *mgr.getValidParameters()
511 );
512 }
513 {
514 Belos::PseudoBlockGmresSolMgr<Scalar,MV_t,LO_t> mgr;
515 solverTypesSL.sublist(PseudoBlockGMRES_name).setParameters(
516 *mgr.getValidParameters()
517 );
518 }
519 if (lapackSupportsScalar) {
520 Belos::BlockCGSolMgr<Scalar,MV_t,LO_t> mgr;
521 solverTypesSL.sublist(BlockCG_name).setParameters(
522 *mgr.getValidParameters()
523 );
524 }
525 if (lapackSupportsScalar) {
526 Belos::PseudoBlockCGSolMgr<Scalar,MV_t,LO_t> mgr;
527 solverTypesSL.sublist(PseudoBlockCG_name).setParameters(
528 *mgr.getValidParameters()
529 );
530 }
531 {
532 Belos::PseudoBlockStochasticCGSolMgr<Scalar,MV_t,LO_t> mgr;
533 solverTypesSL.sublist(PseudoBlockStochasticCG_name).setParameters(
534 *mgr.getValidParameters()
535 );
536 }
537 if (lapackSupportsScalar) {
538 Belos::GCRODRSolMgr<Scalar,MV_t,LO_t> mgr;
539 solverTypesSL.sublist(GCRODR_name).setParameters(
540 *mgr.getValidParameters()
541 );
542 }
543 if (lapackSupportsScalar && scalarIsReal) {
544 Belos::RCGSolMgr<Scalar,MV_t,LO_t> mgr;
545 solverTypesSL.sublist(RCG_name).setParameters(
546 *mgr.getValidParameters()
547 );
548 }
549 {
550 Belos::MinresSolMgr<Scalar,MV_t,LO_t> mgr;
551 solverTypesSL.sublist(MINRES_name).setParameters(
552 *mgr.getValidParameters()
553 );
554 }
555 {
556 Belos::TFQMRSolMgr<Scalar,MV_t,LO_t> mgr;
557 solverTypesSL.sublist(TFQMR_name).setParameters(
558 *mgr.getValidParameters()
559 );
560 }
561 {
562 Belos::BiCGStabSolMgr<Scalar,MV_t,LO_t> mgr;
563 solverTypesSL.sublist(BiCGStab_name).setParameters(
564 *mgr.getValidParameters()
565 );
566 }
567 {
568 Belos::FixedPointSolMgr<Scalar,MV_t,LO_t> mgr;
569 solverTypesSL.sublist(FixedPoint_name).setParameters(
570 *mgr.getValidParameters()
571 );
572 }
573 {
574 Thyra::BelosTpetraGmres<Scalar,MV_t,LO_t> mgr;
575 solverTypesSL.sublist(TpetraGmres_name).setParameters(
576 *mgr.getValidParameters()
577 );
578 }
579 {
580 Thyra::BelosTpetraGmresPipeline<Scalar,MV_t,LO_t> mgr;
581 solverTypesSL.sublist(TpetraGmresPipeline_name).setParameters(
582 *mgr.getValidParameters()
583 );
584 }
585 {
586 Thyra::BelosTpetraGmresSingleReduce<Scalar,MV_t,LO_t> mgr;
587 solverTypesSL.sublist(TpetraGmresSingleReduce_name).setParameters(
588 *mgr.getValidParameters()
589 );
590 }
591 {
592 Thyra::BelosTpetraGmresSstep<Scalar,MV_t,LO_t> mgr;
593 solverTypesSL.sublist(TpetraGmresSstep_name).setParameters(
594 *mgr.getValidParameters()
595 );
596 }
597 }
598 return validParamList;
599}
600
601
602template<class Scalar>
604{
605 thisValidParamList_ = Teuchos::rcp(
606 new Teuchos::ParameterList(*generateAndGetValidParameters())
607 );
608 Teuchos::setupVerboseObjectSublist(&*thisValidParamList_);
609}
610
611
612template<class Scalar>
614 const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
615 const RCP<const LinearOpSourceBase<Scalar> > &approxFwdOpSrc,
616 const RCP<const PreconditionerBase<Scalar> > &prec_in,
617 const bool reusePrec,
618 LinearOpWithSolveBase<Scalar> *Op,
619 const ESupportSolveUse supportSolveUse
620 ) const
621{
622
623 using Teuchos::rcp;
624 using Teuchos::set_extra_data;
625 typedef Teuchos::ScalarTraits<Scalar> ST;
626 typedef MultiVectorBase<Scalar> MV_t;
627 typedef LinearOpBase<Scalar> LO_t;
628
629 const RCP<Teuchos::FancyOStream> out = this->getOStream();
630 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
631 Teuchos::OSTab tab(out);
632 if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW))
633 *out << "\nEntering Thyra::BelosLinearOpWithSolveFactory<"<<ST::name()<<">::initializeOpImpl(...) ...\n";
634
635 // These lines are changing the verbosity of the preconditioner, which has its own verbose object list,
636 // so I am commenting these out, as it is not the job of the linear solver to dictate preconditioner verbosity.
637 //typedef Teuchos::VerboseObjectTempState<PreconditionerFactoryBase<Scalar> > VOTSPF;
638 //VOTSPF precFactoryOutputTempState(precFactory_,out,verbLevel);
639
640 TEUCHOS_TEST_FOR_EXCEPT(Op==NULL);
641 TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
642 TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc->getOp().get()==NULL);
643 RCP<const LinearOpBase<Scalar> >
644 fwdOp = fwdOpSrc->getOp(),
645 approxFwdOp = ( approxFwdOpSrc.get() ? approxFwdOpSrc->getOp() : Teuchos::null );
646
647 //
648 // Get the BelosLinearOpWithSolve interface
649 //
650
651 BelosLinearOpWithSolve<Scalar>
652 *belosOp = &Teuchos::dyn_cast<BelosLinearOpWithSolve<Scalar> >(*Op);
653
654 //
655 // Get/Create the preconditioner
656 //
657
658 RCP<PreconditionerBase<Scalar> > myPrec = Teuchos::null;
659 RCP<const PreconditionerBase<Scalar> > prec = Teuchos::null;
660 if(prec_in.get()) {
661 // Use an externally defined preconditioner
662 prec = prec_in;
663 }
664 else {
665 // Try and generate a preconditioner on our own
666 if(precFactory_.get()) {
667 myPrec =
668 ( !belosOp->isExternalPrec()
669 ? Teuchos::rcp_const_cast<PreconditionerBase<Scalar> >(belosOp->extract_prec())
670 : Teuchos::null
671 );
672 bool hasExistingPrec = false;
673 if(myPrec.get()) {
674 hasExistingPrec = true;
675 // ToDo: Get the forward operator and validate that it is the same
676 // operator that is used here!
677 }
678 else {
679 hasExistingPrec = false;
680 myPrec = precFactory_->createPrec();
681 }
682 if( hasExistingPrec && reusePrec ) {
683 // Just reuse the existing preconditioner again!
684 }
685 else {
686 // Update the preconditioner
687 if(approxFwdOp.get())
688 precFactory_->initializePrec(approxFwdOpSrc,&*myPrec);
689 else
690 precFactory_->initializePrec(fwdOpSrc,&*myPrec);
691 }
692 prec = myPrec;
693 }
694 }
695
696 //
697 // Uninitialize the current solver object
698 //
699
700 bool oldIsExternalPrec = false;
701 RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > oldLP = Teuchos::null;
702 RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > oldIterSolver = Teuchos::null;
703 RCP<const LinearOpSourceBase<Scalar> > oldFwdOpSrc = Teuchos::null;
704 RCP<const LinearOpSourceBase<Scalar> > oldApproxFwdOpSrc = Teuchos::null;
705 ESupportSolveUse oldSupportSolveUse = SUPPORT_SOLVE_UNSPECIFIED;
706
707 belosOp->uninitialize( &oldLP, NULL, &oldIterSolver, &oldFwdOpSrc,
708 NULL, &oldIsExternalPrec, &oldApproxFwdOpSrc, &oldSupportSolveUse );
709
710 //
711 // Create the Belos linear problem
712 // NOTE: If one exists already, reuse it.
713 //
714
715 typedef Belos::LinearProblem<Scalar,MV_t,LO_t> LP_t;
716 RCP<LP_t> lp;
717 if (oldLP != Teuchos::null) {
718 lp = oldLP;
719 }
720 else {
721 lp = rcp(new LP_t());
722 }
723
724 //
725 // Set the operator
726 //
727
728 lp->setOperator(fwdOp);
729
730 //
731 // Set the preconditioner
732 //
733
734 if(prec.get()) {
735 RCP<const LinearOpBase<Scalar> > unspecified = prec->getUnspecifiedPrecOp();
736 RCP<const LinearOpBase<Scalar> > left = prec->getLeftPrecOp();
737 RCP<const LinearOpBase<Scalar> > right = prec->getRightPrecOp();
738 TEUCHOS_TEST_FOR_EXCEPTION(
739 !( left.get() || right.get() || unspecified.get() ), std::logic_error
740 ,"Error, at least one preconditoner linear operator objects must be set!"
741 );
742 if (nonnull(unspecified)) {
743 if (paramList_->get<bool>(LeftPreconditionerIfUnspecified_name, false))
744 lp->setLeftPrec(unspecified);
745 else
746 lp->setRightPrec(unspecified);
747 }
748 else if (nonnull(left)) {
749 lp->setLeftPrec(left);
750 }
751 else if (nonnull(right)) {
752 lp->setRightPrec(right);
753 }
754 else {
755 // Set a left, right or split preconditioner
756 TEUCHOS_TEST_FOR_EXCEPTION(
757 nonnull(left) && nonnull(right),std::logic_error
758 ,"Error, we can not currently handle split preconditioners!"
759 );
760 }
761 }
762 if(myPrec.get()) {
763 set_extra_data<RCP<PreconditionerBase<Scalar> > >(myPrec,"Belos::InternalPrec",
764 Teuchos::inOutArg(lp), Teuchos::POST_DESTROY, false);
765 }
766 else if(prec.get()) {
767 set_extra_data<RCP<const PreconditionerBase<Scalar> > >(prec,"Belos::ExternalPrec",
768 Teuchos::inOutArg(lp), Teuchos::POST_DESTROY, false);
769 }
770
771 //
772 // Generate the parameter list.
773 //
774
775 typedef Belos::SolverManager<Scalar,MV_t,LO_t> IterativeSolver_t;
776 RCP<IterativeSolver_t> iterativeSolver = Teuchos::null;
777 RCP<Teuchos::ParameterList> solverPL = Teuchos::rcp( new Teuchos::ParameterList() );
778
779 switch(solverType_) {
781 {
782 // Set the PL
783 if(paramList_.get()) {
784 Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
785 Teuchos::ParameterList &gmresPL = solverTypesPL.sublist(BlockGMRES_name);
786 solverPL = Teuchos::rcp( &gmresPL, false );
787 }
788 // Create the solver
789 if (oldIterSolver != Teuchos::null) {
790 iterativeSolver = oldIterSolver;
791 iterativeSolver->setProblem( lp );
792 iterativeSolver->setParameters( solverPL );
793 }
794 else {
795 iterativeSolver = rcp(new Belos::BlockGmresSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
796 }
797 break;
798 }
800 {
801 // Set the PL
802 if(paramList_.get()) {
803 Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
804 Teuchos::ParameterList &pbgmresPL = solverTypesPL.sublist(PseudoBlockGMRES_name);
805 solverPL = Teuchos::rcp( &pbgmresPL, false );
806 }
807 //
808 // Create the solver
809 //
810 if (oldIterSolver != Teuchos::null) {
811 iterativeSolver = oldIterSolver;
812 iterativeSolver->setProblem( lp );
813 iterativeSolver->setParameters( solverPL );
814 }
815 else {
816 iterativeSolver = rcp(new Belos::PseudoBlockGmresSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
817 }
818 break;
819 }
821 {
822 // Set the PL
823 if(paramList_.get()) {
824 Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
825 Teuchos::ParameterList &cgPL = solverTypesPL.sublist(BlockCG_name);
826 solverPL = Teuchos::rcp( &cgPL, false );
827 }
828 // Create the solver
829 if (oldIterSolver != Teuchos::null) {
830 iterativeSolver = oldIterSolver;
831 iterativeSolver->setProblem( lp );
832 iterativeSolver->setParameters( solverPL );
833 }
834 else {
835 iterativeSolver = rcp(new Belos::BlockCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
836 }
837 break;
838 }
840 {
841 // Set the PL
842 if(paramList_.get()) {
843 Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
844 Teuchos::ParameterList &pbcgPL = solverTypesPL.sublist(PseudoBlockCG_name);
845 solverPL = Teuchos::rcp( &pbcgPL, false );
846 }
847 //
848 // Create the solver
849 //
850 if (oldIterSolver != Teuchos::null) {
851 iterativeSolver = oldIterSolver;
852 iterativeSolver->setProblem( lp );
853 iterativeSolver->setParameters( solverPL );
854 }
855 else {
856 iterativeSolver = rcp(new Belos::PseudoBlockCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
857 }
858 break;
859 }
861 {
862 // Set the PL
863 if(paramList_.get()) {
864 Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
865 Teuchos::ParameterList &pbcgPL = solverTypesPL.sublist(PseudoBlockStochasticCG_name);
866 solverPL = Teuchos::rcp( &pbcgPL, false );
867 }
868 //
869 // Create the solver
870 //
871 if (oldIterSolver != Teuchos::null) {
872 iterativeSolver = oldIterSolver;
873 iterativeSolver->setProblem( lp );
874 iterativeSolver->setParameters( solverPL );
875 }
876 else {
877 iterativeSolver = rcp(new Belos::PseudoBlockStochasticCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
878 }
879 break;
880 }
882 {
883 // Set the PL
884 if(paramList_.get()) {
885 Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
886 Teuchos::ParameterList &gcrodrPL = solverTypesPL.sublist(GCRODR_name);
887 solverPL = Teuchos::rcp( &gcrodrPL, false );
888 }
889 // Create the solver
890 if (oldIterSolver != Teuchos::null) {
891 iterativeSolver = oldIterSolver;
892 iterativeSolver->setProblem( lp );
893 iterativeSolver->setParameters( solverPL );
894 }
895 else {
896 iterativeSolver = rcp(new Belos::GCRODRSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
897 }
898 break;
899 }
900 case SOLVER_TYPE_RCG:
901 {
902 // Set the PL
903 if(paramList_.get()) {
904 Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
905 Teuchos::ParameterList &rcgPL = solverTypesPL.sublist(RCG_name);
906 solverPL = Teuchos::rcp( &rcgPL, false );
907 }
908 // Create the solver
909 if (oldIterSolver != Teuchos::null) {
910 iterativeSolver = oldIterSolver;
911 iterativeSolver->setProblem( lp );
912 iterativeSolver->setParameters( solverPL );
913 }
914 else {
915 iterativeSolver = rcp(new Belos::RCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
916 }
917 break;
918 }
920 {
921 // Set the PL
922 if(paramList_.get()) {
923 Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
924 Teuchos::ParameterList &minresPL = solverTypesPL.sublist(MINRES_name);
925 solverPL = Teuchos::rcp( &minresPL, false );
926 }
927 // Create the solver
928 if (oldIterSolver != Teuchos::null) {
929 iterativeSolver = oldIterSolver;
930 iterativeSolver->setProblem( lp );
931 iterativeSolver->setParameters( solverPL );
932 }
933 else {
934 iterativeSolver = rcp(new Belos::MinresSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
935 }
936 break;
937 }
939 {
940 // Set the PL
941 if(paramList_.get()) {
942 Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
943 Teuchos::ParameterList &minresPL = solverTypesPL.sublist(TFQMR_name);
944 solverPL = Teuchos::rcp( &minresPL, false );
945 }
946 // Create the solver
947 if (oldIterSolver != Teuchos::null) {
948 iterativeSolver = oldIterSolver;
949 iterativeSolver->setProblem( lp );
950 iterativeSolver->setParameters( solverPL );
951 }
952 else {
953 iterativeSolver = rcp(new Belos::TFQMRSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
954 }
955 break;
956 }
958 {
959 // Set the PL
960 if(paramList_.get()) {
961 Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
962 Teuchos::ParameterList &bicgstabPL = solverTypesPL.sublist(BiCGStab_name);
963 solverPL = Teuchos::rcp( &bicgstabPL, false );
964 }
965 // Create the solver
966 if (oldIterSolver != Teuchos::null) {
967 iterativeSolver = oldIterSolver;
968 iterativeSolver->setProblem( lp );
969 iterativeSolver->setParameters( solverPL );
970 }
971 else {
972 iterativeSolver = rcp(new Belos::BiCGStabSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
973 }
974 break;
975 }
977 {
978 // Set the PL
979 if(paramList_.get()) {
980 Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
981 Teuchos::ParameterList &fixedPointPL = solverTypesPL.sublist(FixedPoint_name);
982 solverPL = Teuchos::rcp( &fixedPointPL, false );
983 }
984 // Create the solver
985 if (oldIterSolver != Teuchos::null) {
986 iterativeSolver = oldIterSolver;
987 iterativeSolver->setProblem( lp );
988 iterativeSolver->setParameters( solverPL );
989 }
990 else {
991 iterativeSolver = rcp(new Belos::FixedPointSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
992 }
993 break;
994 }
996 {
997 // Get the PL
998 if(paramList_.get()) {
999 Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
1000 Teuchos::ParameterList &tpetraGmresPL = solverTypesPL.sublist(TpetraGmres_name);
1001 solverPL = Teuchos::rcp( &tpetraGmresPL, false );
1002 }
1003 // Create the solver, always
1004 iterativeSolver = rcp(new Thyra::BelosTpetraGmres<Scalar,MV_t,LO_t>());
1005 // Set the Problem & PL
1006 iterativeSolver->setProblem( lp );
1007 iterativeSolver->setParameters( solverPL );
1008 break;
1009 }
1011 {
1012 // Get the PL
1013 if(paramList_.get()) {
1014 Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
1015 Teuchos::ParameterList &tpetraGmresPipelinePL = solverTypesPL.sublist(TpetraGmresPipeline_name);
1016 solverPL = Teuchos::rcp( &tpetraGmresPipelinePL, false );
1017 }
1018 // Create the solver, always
1019 iterativeSolver = rcp(new Thyra::BelosTpetraGmresPipeline<Scalar,MV_t,LO_t>());
1020 // Set the Problem & PL
1021 iterativeSolver->setProblem( lp );
1022 iterativeSolver->setParameters( solverPL );
1023 break;
1024 }
1026 {
1027 // Get the PL
1028 if(paramList_.get()) {
1029 Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
1030 Teuchos::ParameterList &tpetraGmresSingleReducePL = solverTypesPL.sublist(TpetraGmresSingleReduce_name);
1031 solverPL = Teuchos::rcp( &tpetraGmresSingleReducePL, false );
1032 }
1033 // Create the solver, always
1034 iterativeSolver = rcp(new Thyra::BelosTpetraGmresSingleReduce<Scalar,MV_t,LO_t>());
1035 // Set the Problem & PL
1036 iterativeSolver->setProblem( lp );
1037 iterativeSolver->setParameters( solverPL );
1038 break;
1039 }
1041 {
1042 // Get the PL
1043 if(paramList_.get()) {
1044 Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
1045 Teuchos::ParameterList &tpetraGmresSstepPL = solverTypesPL.sublist(TpetraGmresSstep_name);
1046 solverPL = Teuchos::rcp( &tpetraGmresSstepPL, false );
1047 }
1048 // Create the solver, always
1049 iterativeSolver = rcp(new Thyra::BelosTpetraGmresSstep<Scalar,MV_t,LO_t>());
1050 // Set the Problem & PL
1051 iterativeSolver->setProblem( lp );
1052 iterativeSolver->setParameters( solverPL );
1053 break;
1054 }
1055
1056 default:
1057 {
1058 TEUCHOS_TEST_FOR_EXCEPT(true);
1059 }
1060 }
1061
1062 //
1063 // Initialize the LOWS object
1064 //
1065
1066 belosOp->initialize(
1067 lp, solverPL, iterativeSolver,
1068 fwdOpSrc, prec, myPrec.get()==NULL, approxFwdOpSrc,
1069 supportSolveUse, convergenceTestFrequency_
1070 );
1071 belosOp->setOStream(out);
1072 belosOp->setVerbLevel(verbLevel);
1073#ifdef TEUCHOS_DEBUG
1074 if(paramList_.get()) {
1075 // Make sure we read the list correctly
1076 paramList_->validateParameters(*this->getValidParameters(),1); // Validate 0th and 1st level deep
1077 }
1078#endif
1079 if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW))
1080 *out << "\nLeaving Thyra::BelosLinearOpWithSolveFactory<"<<ST::name()<<">::initializeOpImpl(...) ...\n";
1081
1082}
1083
1084
1085} // namespace Thyra
1086
1087
1088#endif // THYRA_BELOS_LINEAR_OP_WITH_SOLVE_FACTORY_HPP
Thyra specializations of MultiVecTraits and OperatorTraits.
LinearOpWithSolveFactoryBase subclass implemented in terms of Belos.
bool isCompatible(const LinearOpSourceBase< Scalar > &fwdOpSrc) const
void setPreconditionerFactory(const Teuchos::RCP< PreconditionerFactoryBase< Scalar > > &precFactory, const std::string &precFactoryName)
BelosLinearOpWithSolveFactory()
Construct without preconditioner factory.