Stratimikos Package Browser (Single Doxygen Collection) Version of the Day
Thyra_BelosLinearOpWithSolveFactory_def.hpp
Go to the documentation of this file.
00001 
00002 #ifndef THYRA_BELOS_LINEAR_OP_WITH_SOLVE_FACTORY_HPP
00003 #define THYRA_BELOS_LINEAR_OP_WITH_SOLVE_FACTORY_HPP
00004 
00005 
00006 #include "Thyra_BelosLinearOpWithSolveFactory_decl.hpp"
00007 #include "Thyra_BelosLinearOpWithSolve.hpp"
00008 #include "Thyra_ScaledAdjointLinearOpBase.hpp"
00009 #include "BelosBlockGmresSolMgr.hpp"
00010 #include "BelosPseudoBlockGmresSolMgr.hpp"
00011 #include "BelosBlockCGSolMgr.hpp"
00012 #include "BelosPseudoBlockCGSolMgr.hpp"
00013 #include "BelosGCRODRSolMgr.hpp"
00014 #include "BelosRCGSolMgr.hpp"
00015 #include "BelosMinresSolMgr.hpp"
00016 #include "BelosThyraAdapter.hpp"
00017 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00018 #include "Teuchos_StandardParameterEntryValidators.hpp"
00019 #include "Teuchos_ParameterList.hpp"
00020 #include "Teuchos_dyn_cast.hpp"
00021 #include "Teuchos_ValidatorXMLConverterDB.hpp"
00022 #include "Teuchos_StandardValidatorXMLConverters.hpp"
00023 
00024 
00025 namespace Thyra {
00026 
00027 
00028 // Parameter names for Paramter List
00029 
00030 template<class Scalar>
00031 const std::string BelosLinearOpWithSolveFactory<Scalar>::SolverType_name = "Solver Type";
00032 template<class Scalar>
00033 const std::string BelosLinearOpWithSolveFactory<Scalar>::SolverType_default = "Pseudo Block GMRES";
00034 template<class Scalar>
00035 const std::string BelosLinearOpWithSolveFactory<Scalar>::SolverTypes_name = "Solver Types";
00036 template<class Scalar>
00037 const std::string BelosLinearOpWithSolveFactory<Scalar>::BlockGMRES_name = "Block GMRES";
00038 template<class Scalar>
00039 const std::string BelosLinearOpWithSolveFactory<Scalar>::PseudoBlockGMRES_name = "Pseudo Block GMRES";
00040 template<class Scalar>
00041 const std::string BelosLinearOpWithSolveFactory<Scalar>::BlockCG_name = "Block CG";
00042 template<class Scalar>
00043 const std::string BelosLinearOpWithSolveFactory<Scalar>::PseudoBlockCG_name = "Pseudo Block CG";
00044 template<class Scalar>
00045 const std::string BelosLinearOpWithSolveFactory<Scalar>::GCRODR_name = "GCRODR";
00046 template<class Scalar>
00047 const std::string BelosLinearOpWithSolveFactory<Scalar>::RCG_name = "RCG";
00048 template<class Scalar>
00049 const std::string BelosLinearOpWithSolveFactory<Scalar>::MINRES_name = "MINRES";
00050 template<class Scalar>
00051 const std::string BelosLinearOpWithSolveFactory<Scalar>::ConvergenceTestFrequency_name = "Convergence Test Frequency";
00052 
00053 // Constructors/initializers/accessors
00054 
00055 
00056 template<class Scalar>
00057 BelosLinearOpWithSolveFactory<Scalar>::BelosLinearOpWithSolveFactory()
00058   :solverType_(SOLVER_TYPE_PSEUDO_BLOCK_GMRES),
00059    convergenceTestFrequency_(1)
00060 {
00061   updateThisValidParamList();
00062 }
00063 
00064 
00065 template<class Scalar>
00066 BelosLinearOpWithSolveFactory<Scalar>::BelosLinearOpWithSolveFactory(
00067   const RCP<PreconditionerFactoryBase<Scalar> > &precFactory
00068   )
00069   :solverType_(SOLVER_TYPE_PSEUDO_BLOCK_GMRES)
00070 {
00071   this->setPreconditionerFactory(precFactory, "");
00072 }
00073 
00074 
00075 // Overridden from LinearOpWithSolveFactoryBase
00076 
00077 
00078 template<class Scalar>
00079 bool BelosLinearOpWithSolveFactory<Scalar>::acceptsPreconditionerFactory() const
00080 {
00081   return true;
00082 }
00083 
00084 
00085 template<class Scalar>
00086 void BelosLinearOpWithSolveFactory<Scalar>::setPreconditionerFactory(
00087   const RCP<PreconditionerFactoryBase<Scalar> > &precFactory,
00088   const std::string &precFactoryName
00089   )
00090 {
00091   TEUCHOS_TEST_FOR_EXCEPT(!precFactory.get());
00092   RCP<const Teuchos::ParameterList>
00093     precFactoryValidPL = precFactory->getValidParameters();
00094   const std::string _precFactoryName =
00095     ( precFactoryName != ""
00096       ? precFactoryName
00097       : ( precFactoryValidPL.get() ? precFactoryValidPL->name() : "GENERIC PRECONDITIONER FACTORY" )
00098       );
00099   precFactory_ = precFactory;
00100   precFactoryName_ = _precFactoryName;
00101   updateThisValidParamList();
00102 }
00103 
00104 
00105 template<class Scalar>
00106 RCP<PreconditionerFactoryBase<Scalar> >
00107 BelosLinearOpWithSolveFactory<Scalar>::getPreconditionerFactory() const
00108 {
00109   return precFactory_;
00110 }
00111 
00112 
00113 template<class Scalar>
00114 void BelosLinearOpWithSolveFactory<Scalar>::unsetPreconditionerFactory(
00115   RCP<PreconditionerFactoryBase<Scalar> > *precFactory,
00116   std::string *precFactoryName
00117   )
00118 {
00119   if(precFactory) *precFactory = precFactory_;
00120   if(precFactoryName) *precFactoryName = precFactoryName_;
00121   precFactory_ = Teuchos::null;
00122   precFactoryName_ = "";
00123   updateThisValidParamList();
00124 }
00125 
00126 
00127 template<class Scalar>
00128 bool BelosLinearOpWithSolveFactory<Scalar>::isCompatible(
00129   const LinearOpSourceBase<Scalar> &fwdOpSrc
00130   ) const
00131 {
00132   if(precFactory_.get())
00133     return precFactory_->isCompatible(fwdOpSrc);
00134   return true; // Without a preconditioner, we are compatible with all linear operators!
00135 }
00136 
00137 
00138 template<class Scalar>
00139 RCP<LinearOpWithSolveBase<Scalar> >
00140 BelosLinearOpWithSolveFactory<Scalar>::createOp() const
00141 {
00142   return Teuchos::rcp(new BelosLinearOpWithSolve<Scalar>());
00143 }
00144 
00145 
00146 template<class Scalar>
00147 void BelosLinearOpWithSolveFactory<Scalar>::initializeOp(
00148   const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
00149   LinearOpWithSolveBase<Scalar> *Op,
00150   const ESupportSolveUse supportSolveUse
00151   ) const
00152 {
00153   using Teuchos::null;
00154   initializeOpImpl(fwdOpSrc,null,null,false,Op,supportSolveUse);
00155 }
00156 
00157 
00158 template<class Scalar>
00159 void BelosLinearOpWithSolveFactory<Scalar>::initializeAndReuseOp(
00160   const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
00161   LinearOpWithSolveBase<Scalar> *Op
00162   ) const
00163 {
00164   using Teuchos::null;
00165   initializeOpImpl(fwdOpSrc,null,null,true,Op,SUPPORT_SOLVE_UNSPECIFIED);
00166 }
00167 
00168 
00169 template<class Scalar>
00170 bool BelosLinearOpWithSolveFactory<Scalar>::supportsPreconditionerInputType(
00171   const EPreconditionerInputType precOpType
00172   ) const
00173 {
00174   if(precFactory_.get())
00175     return true;
00176   return (precOpType==PRECONDITIONER_INPUT_TYPE_AS_OPERATOR);
00177 }
00178 
00179 
00180 template<class Scalar>
00181 void BelosLinearOpWithSolveFactory<Scalar>::initializePreconditionedOp(
00182   const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
00183   const RCP<const PreconditionerBase<Scalar> > &prec,
00184   LinearOpWithSolveBase<Scalar> *Op,
00185   const ESupportSolveUse supportSolveUse
00186   ) const
00187 {
00188   using Teuchos::null;
00189   initializeOpImpl(fwdOpSrc,null,prec,false,Op,supportSolveUse);
00190 }
00191 
00192 
00193 template<class Scalar>
00194 void BelosLinearOpWithSolveFactory<Scalar>::initializeApproxPreconditionedOp(
00195   const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
00196   const RCP<const LinearOpSourceBase<Scalar> > &approxFwdOpSrc,
00197   LinearOpWithSolveBase<Scalar> *Op,
00198   const ESupportSolveUse supportSolveUse
00199   ) const
00200 {
00201   using Teuchos::null;
00202   initializeOpImpl(fwdOpSrc,approxFwdOpSrc,null,false,Op,supportSolveUse);
00203 }
00204 
00205 
00206 template<class Scalar>
00207 void BelosLinearOpWithSolveFactory<Scalar>::uninitializeOp(
00208   LinearOpWithSolveBase<Scalar> *Op,
00209   RCP<const LinearOpSourceBase<Scalar> > *fwdOpSrc,
00210   RCP<const PreconditionerBase<Scalar> > *prec,
00211   RCP<const LinearOpSourceBase<Scalar> > *approxFwdOpSrc,
00212   ESupportSolveUse *supportSolveUse
00213   ) const
00214 {
00215 #ifdef TEUCHOS_DEBUG
00216   TEUCHOS_TEST_FOR_EXCEPT(Op==NULL);
00217 #endif
00218   BelosLinearOpWithSolve<Scalar>
00219     &belosOp = Teuchos::dyn_cast<BelosLinearOpWithSolve<Scalar> >(*Op);
00220   RCP<const LinearOpSourceBase<Scalar> > 
00221     _fwdOpSrc = belosOp.extract_fwdOpSrc();
00222   RCP<const PreconditionerBase<Scalar> >
00223     _prec = ( belosOp.isExternalPrec() ? belosOp.extract_prec() : Teuchos::null );
00224   // Note: above we only extract the preconditioner if it was passed in
00225   // externally.  Otherwise, we need to hold on to it so that we can reuse it
00226   // in the next initialization.
00227   RCP<const LinearOpSourceBase<Scalar> >
00228     _approxFwdOpSrc = belosOp.extract_approxFwdOpSrc();
00229   ESupportSolveUse
00230     _supportSolveUse = belosOp.supportSolveUse();
00231   if(fwdOpSrc) *fwdOpSrc = _fwdOpSrc;
00232   if(prec) *prec = _prec;
00233   if(approxFwdOpSrc) *approxFwdOpSrc = _approxFwdOpSrc;
00234   if(supportSolveUse) *supportSolveUse = _supportSolveUse;
00235 }
00236 
00237 
00238 // Overridden from ParameterListAcceptor
00239 
00240 
00241 template<class Scalar>
00242 void BelosLinearOpWithSolveFactory<Scalar>::setParameterList(
00243   RCP<Teuchos::ParameterList> const& paramList
00244   )
00245 {
00246   TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==NULL);
00247   paramList->validateParametersAndSetDefaults(*this->getValidParameters(), 1);
00248   paramList_ = paramList;
00249   solverType_ =
00250     Teuchos::getIntegralValue<EBelosSolverType>(*paramList_, SolverType_name);
00251   convergenceTestFrequency_ =
00252     Teuchos::getParameter<int>(*paramList_, ConvergenceTestFrequency_name);
00253   Teuchos::readVerboseObjectSublist(&*paramList_,this);
00254 }
00255 
00256 
00257 template<class Scalar>
00258 RCP<Teuchos::ParameterList>
00259 BelosLinearOpWithSolveFactory<Scalar>::getNonconstParameterList()
00260 {
00261   return paramList_;
00262 }
00263 
00264 
00265 template<class Scalar>
00266 RCP<Teuchos::ParameterList>
00267 BelosLinearOpWithSolveFactory<Scalar>::unsetParameterList()
00268 {
00269   RCP<Teuchos::ParameterList> _paramList = paramList_;
00270   paramList_ = Teuchos::null;
00271   return _paramList;
00272 }
00273 
00274 
00275 template<class Scalar>
00276 RCP<const Teuchos::ParameterList>
00277 BelosLinearOpWithSolveFactory<Scalar>::getParameterList() const
00278 {
00279   return paramList_;
00280 }
00281 
00282 
00283 template<class Scalar>
00284 RCP<const Teuchos::ParameterList>
00285 BelosLinearOpWithSolveFactory<Scalar>::getValidParameters() const
00286 {
00287   return thisValidParamList_;
00288 }
00289 
00290 
00291 // Public functions overridden from Teuchos::Describable
00292 
00293 
00294 template<class Scalar>
00295 std::string BelosLinearOpWithSolveFactory<Scalar>::description() const
00296 {
00297   std::ostringstream oss;
00298   oss << "Thyra::BelosLinearOpWithSolveFactory";
00299   //oss << "{";
00300   // ToDo: Fill this in some!
00301   //oss << "}";
00302   return oss.str();
00303 }
00304 
00305 
00306 // private
00307 
00308 
00309 template<class Scalar>
00310 RCP<const Teuchos::ParameterList>
00311 BelosLinearOpWithSolveFactory<Scalar>::generateAndGetValidParameters()
00312 {
00313   using Teuchos::as;
00314   using Teuchos::tuple;
00315   using Teuchos::setStringToIntegralParameter;
00316 Teuchos::ValidatorXMLConverterDB::addConverter(
00317   Teuchos::DummyObjectGetter<
00318     Teuchos::StringToIntegralParameterEntryValidator<EBelosSolverType> 
00319   >::getDummyObject(),
00320   Teuchos::DummyObjectGetter<Teuchos::StringToIntegralValidatorXMLConverter<
00321     EBelosSolverType> >::getDummyObject());
00322 
00323   typedef MultiVectorBase<Scalar> MV_t;
00324   typedef LinearOpBase<Scalar> LO_t;
00325   static RCP<Teuchos::ParameterList> validParamList;
00326   if(validParamList.get()==NULL) {
00327     validParamList = Teuchos::rcp(new Teuchos::ParameterList("BelosLinearOpWithSolveFactory"));
00328     setStringToIntegralParameter<EBelosSolverType>(
00329       SolverType_name, SolverType_default,
00330       "Type of linear solver algorithm to use.",
00331       tuple<std::string>(
00332         "Block GMRES",
00333         "Pseudo Block GMRES",
00334         "Block CG",
00335         "Pseudo Block CG",
00336         "GCRODR",
00337   "RCG",
00338         "MINRES"
00339         ),
00340       tuple<std::string>(
00341         "Block GMRES solver for nonsymmetric linear systems.  It can also solve "
00342   "single right-hand side systems, and can also perform Flexible GMRES "
00343   "(where the preconditioner may change at every iteration, for example "
00344   "for inner-outer iterations), by setting options in the \"Block GMRES\" "
00345   "sublist.",
00346 
00347         "GMRES solver for nonsymmetric linear systems, that performs single "
00348   "right-hand side solves on multiple right-hand sides at once.  It "
00349   "exploits operator multivector multiplication in order to amortize "
00350         "global communication costs.  Individual linear systems are deflated "
00351   "out as they are solved.",
00352 
00353         "Block CG solver for symmetric (Hermitian in complex arithmetic) "
00354   "positive definite linear systems.  It can also solve single "
00355   "right-hand-side systems.",
00356 
00357         "CG solver that performs single right-hand side CG on multiple right-hand "
00358   "sides at once.  It exploits operator multivector multiplication in order "
00359   "to amortize global communication costs.  Individual linear systems are "
00360   "deflated out as they are solved.",
00361 
00362         "GMRES solver for nonsymmetric linear systems, that performs subspace "
00363   "recycling to accelerate convergence for sequences of related linear "
00364   "systems.",
00365 
00366   "CG solver for symmetric (Hermitian in complex arithmetic) positive "
00367   "definite linear systems, that performs subspace recycling to "
00368   "accelerate convergence for sequences of related linear systems.",
00369 
00370         "MINRES solver for symmetric indefinite linear systems.  It performs "
00371   "single-right-hand-side solves on multiple right-hand sides sequentially."
00372         ),
00373       tuple<EBelosSolverType>(
00374         SOLVER_TYPE_BLOCK_GMRES,
00375         SOLVER_TYPE_PSEUDO_BLOCK_GMRES,
00376         SOLVER_TYPE_BLOCK_CG,
00377         SOLVER_TYPE_PSEUDO_BLOCK_CG,
00378         SOLVER_TYPE_GCRODR,
00379   SOLVER_TYPE_RCG,
00380         SOLVER_TYPE_MINRES
00381         ),
00382       &*validParamList
00383       );
00384     validParamList->set(ConvergenceTestFrequency_name, as<int>(1),
00385       "Number of linear solver iterations to skip between applying"
00386       " user-defined convergence test.");
00387     Teuchos::ParameterList
00388       &solverTypesSL = validParamList->sublist(SolverTypes_name);
00389     {
00390       Belos::BlockGmresSolMgr<Scalar,MV_t,LO_t> mgr;
00391       solverTypesSL.sublist(BlockGMRES_name).setParameters(
00392         *mgr.getValidParameters()
00393         );
00394     }
00395     {
00396       Belos::PseudoBlockGmresSolMgr<Scalar,MV_t,LO_t> mgr;
00397       solverTypesSL.sublist(PseudoBlockGMRES_name).setParameters(
00398         *mgr.getValidParameters()
00399         );
00400     }
00401     {
00402       Belos::BlockCGSolMgr<Scalar,MV_t,LO_t> mgr;
00403       solverTypesSL.sublist(BlockCG_name).setParameters(
00404         *mgr.getValidParameters()
00405         );
00406     }
00407     {
00408       Belos::PseudoBlockCGSolMgr<Scalar,MV_t,LO_t> mgr;
00409       solverTypesSL.sublist(PseudoBlockCG_name).setParameters(
00410         *mgr.getValidParameters()
00411         );
00412     }
00413     {
00414       Belos::GCRODRSolMgr<Scalar,MV_t,LO_t> mgr;
00415       solverTypesSL.sublist(GCRODR_name).setParameters(
00416         *mgr.getValidParameters()
00417         );
00418     }
00419     {
00420       Belos::RCGSolMgr<Scalar,MV_t,LO_t> mgr;
00421       solverTypesSL.sublist(RCG_name).setParameters(
00422         *mgr.getValidParameters()
00423         );
00424     }
00425     {
00426       Belos::MinresSolMgr<Scalar,MV_t,LO_t> mgr;
00427       solverTypesSL.sublist(MINRES_name).setParameters(
00428         *mgr.getValidParameters()
00429         );
00430     }
00431   }
00432   return validParamList;
00433 }
00434 
00435 
00436 template<class Scalar>
00437 void BelosLinearOpWithSolveFactory<Scalar>::updateThisValidParamList()
00438 {
00439   thisValidParamList_ = Teuchos::rcp(
00440     new Teuchos::ParameterList(*generateAndGetValidParameters())
00441     );
00442   Teuchos::setupVerboseObjectSublist(&*thisValidParamList_);
00443 }
00444 
00445 
00446 template<class Scalar>
00447 void BelosLinearOpWithSolveFactory<Scalar>::initializeOpImpl(
00448   const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
00449   const RCP<const LinearOpSourceBase<Scalar> > &approxFwdOpSrc,
00450   const RCP<const PreconditionerBase<Scalar> > &prec_in,
00451   const bool reusePrec,
00452   LinearOpWithSolveBase<Scalar> *Op,
00453   const ESupportSolveUse supportSolveUse
00454   ) const
00455 {
00456 
00457   using Teuchos::rcp;
00458   using Teuchos::set_extra_data;
00459   typedef Teuchos::ScalarTraits<Scalar> ST;
00460   typedef typename ST::magnitudeType ScalarMag;
00461   typedef MultiVectorBase<Scalar> MV_t;
00462   typedef LinearOpBase<Scalar> LO_t;
00463 
00464   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00465   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00466   Teuchos::OSTab tab(out);
00467   if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW))
00468     *out << "\nEntering Thyra::BelosLinearOpWithSolveFactory<"<<ST::name()<<">::initializeOpImpl(...) ...\n";
00469 
00470   // These lines are changing the verbosity of the preconditioner, which has its own verbose object list,
00471   // so I am commenting these out, as it is not the job of the linear solver to dictate preconditioner verbosity.
00472   //typedef Teuchos::VerboseObjectTempState<PreconditionerFactoryBase<Scalar> > VOTSPF;
00473   //VOTSPF precFactoryOutputTempState(precFactory_,out,verbLevel);
00474   
00475   TEUCHOS_TEST_FOR_EXCEPT(Op==NULL);
00476   TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
00477   TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc->getOp().get()==NULL);
00478   RCP<const LinearOpBase<Scalar> >
00479     fwdOp = fwdOpSrc->getOp(),
00480     approxFwdOp = ( approxFwdOpSrc.get() ? approxFwdOpSrc->getOp() : Teuchos::null );
00481 
00482   //
00483   // Get the BelosLinearOpWithSolve interface
00484   //
00485 
00486   BelosLinearOpWithSolve<Scalar>
00487     *belosOp = &Teuchos::dyn_cast<BelosLinearOpWithSolve<Scalar> >(*Op);
00488 
00489   //
00490   // Get/Create the preconditioner
00491   //
00492 
00493   RCP<PreconditionerBase<Scalar> > myPrec = Teuchos::null;
00494   RCP<const PreconditionerBase<Scalar> > prec = Teuchos::null;
00495   if(prec_in.get()) {
00496     // Use an externally defined preconditioner
00497     prec = prec_in;
00498   }
00499   else {
00500     // Try and generate a preconditioner on our own
00501     if(precFactory_.get()) {
00502       myPrec =
00503         ( !belosOp->isExternalPrec()
00504           ? Teuchos::rcp_const_cast<PreconditionerBase<Scalar> >(belosOp->extract_prec())
00505           : Teuchos::null
00506           );
00507       bool hasExistingPrec = false;
00508       if(myPrec.get()) {
00509         hasExistingPrec = true;
00510         // ToDo: Get the forward operator and validate that it is the same
00511         // operator that is used here!
00512       }
00513       else {
00514         hasExistingPrec = false;
00515         myPrec = precFactory_->createPrec();
00516       }
00517       if( hasExistingPrec && reusePrec ) {
00518         // Just reuse the existing preconditioner again!
00519       }
00520       else {
00521         // Update the preconditioner
00522         if(approxFwdOp.get())
00523           precFactory_->initializePrec(approxFwdOpSrc,&*myPrec);
00524         else
00525           precFactory_->initializePrec(fwdOpSrc,&*myPrec);
00526       }
00527       prec = myPrec;
00528     }
00529   }
00530 
00531   //
00532   // Uninitialize the current solver object
00533   //
00534 
00535   bool oldIsExternalPrec = false;
00536   RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > oldLP = Teuchos::null;
00537   RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > oldIterSolver = Teuchos::null;
00538   RCP<const LinearOpSourceBase<Scalar> > oldFwdOpSrc = Teuchos::null;
00539   RCP<const LinearOpSourceBase<Scalar> > oldApproxFwdOpSrc = Teuchos::null;   
00540   ESupportSolveUse oldSupportSolveUse = SUPPORT_SOLVE_UNSPECIFIED;
00541 
00542   belosOp->uninitialize( &oldLP, NULL, &oldIterSolver, &oldFwdOpSrc,
00543     NULL, &oldIsExternalPrec, &oldApproxFwdOpSrc, &oldSupportSolveUse );
00544 
00545   //
00546   // Create the Belos linear problem
00547   // NOTE:  If one exists already, reuse it.
00548   //
00549 
00550   typedef Belos::LinearProblem<Scalar,MV_t,LO_t> LP_t;
00551   RCP<LP_t> lp;
00552   if (oldLP != Teuchos::null) {
00553     lp = oldLP;
00554   }
00555   else {
00556     lp = rcp(new LP_t());
00557   }
00558 
00559   //
00560   // Set the operator
00561   //
00562 
00563   lp->setOperator(fwdOp);
00564 
00565   //
00566   // Set the preconditioner
00567   //
00568 
00569   if(prec.get()) {
00570     RCP<const LinearOpBase<Scalar> > unspecified = prec->getUnspecifiedPrecOp();
00571     RCP<const LinearOpBase<Scalar> > left = prec->getLeftPrecOp();
00572     RCP<const LinearOpBase<Scalar> > right = prec->getRightPrecOp();
00573     TEUCHOS_TEST_FOR_EXCEPTION(
00574       !( left.get() || right.get() || unspecified.get() ), std::logic_error
00575       ,"Error, at least one preconditoner linear operator objects must be set!"
00576       );
00577     if(unspecified.get()) {
00578       lp->setRightPrec(unspecified);
00579       // ToDo: Allow user to determine whether this should be placed on the
00580       // left or on the right through a parameter in the parameter list!
00581     }
00582     else {
00583       // Set a left, right or split preconditioner
00584       TEUCHOS_TEST_FOR_EXCEPTION(
00585         left.get(),std::logic_error
00586         ,"Error, we can not currently handle a left preconditioner!"
00587         );
00588       lp->setRightPrec(right);
00589     }
00590   }
00591   if(myPrec.get()) {
00592     set_extra_data<RCP<PreconditionerBase<Scalar> > >(myPrec,"Belos::InternalPrec",
00593       Teuchos::inOutArg(lp), Teuchos::POST_DESTROY, false);
00594   }
00595   else if(prec.get()) {
00596     set_extra_data<RCP<const PreconditionerBase<Scalar> > >(prec,"Belos::ExternalPrec",
00597       Teuchos::inOutArg(lp), Teuchos::POST_DESTROY, false);
00598   }
00599 
00600   //
00601   // Generate the parameter list.
00602   //
00603 
00604   typedef Belos::SolverManager<Scalar,MV_t,LO_t> IterativeSolver_t;
00605   RCP<IterativeSolver_t> iterativeSolver = Teuchos::null;
00606   RCP<Teuchos::ParameterList> solverPL = Teuchos::rcp( new Teuchos::ParameterList() );
00607   
00608   switch(solverType_) {
00609     case SOLVER_TYPE_BLOCK_GMRES: 
00610     {
00611       // Set the PL
00612       if(paramList_.get()) {
00613         Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
00614         Teuchos::ParameterList &gmresPL = solverTypesPL.sublist(BlockGMRES_name);
00615         solverPL = Teuchos::rcp( &gmresPL, false );
00616       }
00617       // Create the solver
00618       if (oldIterSolver != Teuchos::null) {
00619         iterativeSolver = oldIterSolver;
00620         iterativeSolver->setProblem( lp );
00621         iterativeSolver->setParameters( solverPL );
00622       } 
00623       else {
00624         iterativeSolver = rcp(new Belos::BlockGmresSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
00625       }
00626       break;
00627     }
00628     case SOLVER_TYPE_PSEUDO_BLOCK_GMRES:
00629     {
00630       // Set the PL
00631       if(paramList_.get()) {
00632         Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
00633         Teuchos::ParameterList &pbgmresPL = solverTypesPL.sublist(PseudoBlockGMRES_name);
00634         solverPL = Teuchos::rcp( &pbgmresPL, false );
00635       }
00636       // 
00637       // Create the solver
00638       // 
00639       if (oldIterSolver != Teuchos::null) {
00640         iterativeSolver = oldIterSolver;
00641         iterativeSolver->setProblem( lp );
00642         iterativeSolver->setParameters( solverPL );
00643       }
00644       else {
00645         iterativeSolver = rcp(new Belos::PseudoBlockGmresSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
00646       }
00647       break;
00648     }
00649     case SOLVER_TYPE_BLOCK_CG:
00650     {
00651       // Set the PL
00652       if(paramList_.get()) {
00653         Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
00654         Teuchos::ParameterList &cgPL = solverTypesPL.sublist(BlockCG_name);
00655         solverPL = Teuchos::rcp( &cgPL, false );
00656       }
00657       // Create the solver
00658       if (oldIterSolver != Teuchos::null) {
00659         iterativeSolver = oldIterSolver;
00660         iterativeSolver->setProblem( lp );
00661         iterativeSolver->setParameters( solverPL );
00662       }
00663       else {
00664         iterativeSolver = rcp(new Belos::BlockCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
00665       }
00666       break;
00667     }
00668     case SOLVER_TYPE_PSEUDO_BLOCK_CG:
00669     {
00670       // Set the PL
00671       if(paramList_.get()) {
00672         Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
00673         Teuchos::ParameterList &pbcgPL = solverTypesPL.sublist(PseudoBlockCG_name);
00674         solverPL = Teuchos::rcp( &pbcgPL, false );
00675       }
00676       // 
00677       // Create the solver
00678       // 
00679       if (oldIterSolver != Teuchos::null) {
00680         iterativeSolver = oldIterSolver;
00681         iterativeSolver->setProblem( lp );
00682         iterativeSolver->setParameters( solverPL );
00683       }
00684       else {
00685         iterativeSolver = rcp(new Belos::PseudoBlockCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
00686       }
00687       break;
00688     }
00689     case SOLVER_TYPE_GCRODR:
00690     {
00691       // Set the PL
00692       if(paramList_.get()) {
00693         Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
00694         Teuchos::ParameterList &gcrodrPL = solverTypesPL.sublist(GCRODR_name);
00695         solverPL = Teuchos::rcp( &gcrodrPL, false );
00696       }
00697       // Create the solver
00698       if (oldIterSolver != Teuchos::null) {
00699         iterativeSolver = oldIterSolver;
00700         iterativeSolver->setProblem( lp );
00701         iterativeSolver->setParameters( solverPL );
00702       } 
00703       else {
00704         iterativeSolver = rcp(new Belos::GCRODRSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
00705       }
00706       break;
00707     }
00708     case SOLVER_TYPE_RCG:
00709     {
00710       // Set the PL
00711       if(paramList_.get()) {
00712         Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
00713         Teuchos::ParameterList &rcgPL = solverTypesPL.sublist(RCG_name);
00714         solverPL = Teuchos::rcp( &rcgPL, false );
00715       }
00716       // Create the solver
00717       if (oldIterSolver != Teuchos::null) {
00718         iterativeSolver = oldIterSolver;
00719         iterativeSolver->setProblem( lp );
00720         iterativeSolver->setParameters( solverPL );
00721       } 
00722       else {
00723         iterativeSolver = rcp(new Belos::RCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
00724       }
00725       break;
00726     }
00727     case SOLVER_TYPE_MINRES:
00728     {
00729       // Set the PL
00730       if(paramList_.get()) {
00731         Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
00732         Teuchos::ParameterList &minresPL = solverTypesPL.sublist(MINRES_name);
00733         solverPL = Teuchos::rcp( &minresPL, false );
00734       }
00735       // Create the solver
00736       if (oldIterSolver != Teuchos::null) {
00737         iterativeSolver = oldIterSolver;
00738         iterativeSolver->setProblem( lp );
00739         iterativeSolver->setParameters( solverPL );
00740       }
00741       else {
00742         iterativeSolver = rcp(new Belos::MinresSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
00743       }
00744       break;
00745     }
00746 
00747     default:
00748     {
00749       TEUCHOS_TEST_FOR_EXCEPT(true);
00750     }
00751   }
00752 
00753   //
00754   // Initialize the LOWS object
00755   //
00756 
00757   belosOp->initialize(
00758     lp, solverPL, iterativeSolver,
00759     fwdOpSrc, prec, myPrec.get()==NULL, approxFwdOpSrc,
00760     supportSolveUse, convergenceTestFrequency_
00761     );
00762   belosOp->setOStream(out);
00763   belosOp->setVerbLevel(verbLevel);
00764 #ifdef TEUCHOS_DEBUG
00765   if(paramList_.get()) {
00766     // Make sure we read the list correctly
00767     paramList_->validateParameters(*this->getValidParameters(),1); // Validate 0th and 1st level deep
00768   }
00769 #endif
00770   if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW))
00771     *out << "\nLeaving Thyra::BelosLinearOpWithSolveFactory<"<<ST::name()<<">::initializeOpImpl(...) ...\n";
00772   
00773 }
00774 
00775 
00776 } // namespace Thyra
00777 
00778 
00779 #endif // THYRA_BELOS_LINEAR_OP_WITH_SOLVE_FACTORY_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines