Stratimikos Version of the Day
Thyra_BelosLinearOpWithSolveFactory_def.hpp
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   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   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   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         "Performs block and single single-RHS GMRES as well as\n"
00342         "flexible GMRES by setting options in the \"Block GMRES\" sublist.",
00343 
00344         "GMRES solver that performs single-RHS GMRES on multiple RHSs taking\n"
00345         "advantage of operator multi-vector multiplication and the amortization\n"
00346         "of global communication.  Individual linear systems are deflated out as\n"
00347         "they are solved.",
00348 
00349         "CG solver that performs block and single-RHS CG.",
00350 
00351         "CG solver that performs single-RHS CG on multiple RHSs taking\n"
00352         "advantage of operator multi-vector multiplication and the amortization\n"
00353         "of global communication.  Individual linear systems are deflated out as\n"
00354         "they are solved.",
00355 
00356         "GMRES solver that performs subspace recycling between RHS and linear systems.",
00357   "CG solver that performs subspace recycling between RHS and linear systems.",
00358 
00359         "MINRES solver that performs single-RHS MINRES on multiple RHSs sequentially."
00360         ),
00361       tuple<EBelosSolverType>(
00362         SOLVER_TYPE_BLOCK_GMRES,
00363         SOLVER_TYPE_PSEUDO_BLOCK_GMRES,
00364         SOLVER_TYPE_BLOCK_CG,
00365         SOLVER_TYPE_PSEUDO_BLOCK_CG,
00366         SOLVER_TYPE_GCRODR,
00367   SOLVER_TYPE_RCG,
00368         SOLVER_TYPE_MINRES
00369         ),
00370       &*validParamList
00371       );
00372     validParamList->set(ConvergenceTestFrequency_name, as<int>(1),
00373       "Number of linear solver iterations to skip between applying"
00374       " user-defined convergence test.");
00375     Teuchos::ParameterList
00376       &solverTypesSL = validParamList->sublist(SolverTypes_name);
00377     {
00378       Belos::BlockGmresSolMgr<Scalar,MV_t,LO_t> mgr;
00379       solverTypesSL.sublist(BlockGMRES_name).setParameters(
00380         *mgr.getValidParameters()
00381         );
00382     }
00383     {
00384       Belos::PseudoBlockGmresSolMgr<Scalar,MV_t,LO_t> mgr;
00385       solverTypesSL.sublist(PseudoBlockGMRES_name).setParameters(
00386         *mgr.getValidParameters()
00387         );
00388     }
00389     {
00390       Belos::BlockCGSolMgr<Scalar,MV_t,LO_t> mgr;
00391       solverTypesSL.sublist(BlockCG_name).setParameters(
00392         *mgr.getValidParameters()
00393         );
00394     }
00395     {
00396       Belos::PseudoBlockCGSolMgr<Scalar,MV_t,LO_t> mgr;
00397       solverTypesSL.sublist(PseudoBlockCG_name).setParameters(
00398         *mgr.getValidParameters()
00399         );
00400     }
00401     {
00402       Belos::GCRODRSolMgr<Scalar,MV_t,LO_t> mgr;
00403       solverTypesSL.sublist(GCRODR_name).setParameters(
00404         *mgr.getValidParameters()
00405         );
00406     }
00407     {
00408       Belos::RCGSolMgr<Scalar,MV_t,LO_t> mgr;
00409       solverTypesSL.sublist(RCG_name).setParameters(
00410         *mgr.getValidParameters()
00411         );
00412     }
00413     {
00414       Belos::MinresSolMgr<Scalar,MV_t,LO_t> mgr;
00415       solverTypesSL.sublist(MINRES_name).setParameters(
00416         *mgr.getValidParameters()
00417         );
00418     }
00419   }
00420   return validParamList;
00421 }
00422 
00423 
00424 template<class Scalar>
00425 void BelosLinearOpWithSolveFactory<Scalar>::updateThisValidParamList()
00426 {
00427   thisValidParamList_ = Teuchos::rcp(
00428     new Teuchos::ParameterList(*generateAndGetValidParameters())
00429     );
00430   Teuchos::setupVerboseObjectSublist(&*thisValidParamList_);
00431 }
00432 
00433 
00434 template<class Scalar>
00435 void BelosLinearOpWithSolveFactory<Scalar>::initializeOpImpl(
00436   const RCP<const LinearOpSourceBase<Scalar> > &fwdOpSrc,
00437   const RCP<const LinearOpSourceBase<Scalar> > &approxFwdOpSrc,
00438   const RCP<const PreconditionerBase<Scalar> > &prec_in,
00439   const bool reusePrec,
00440   LinearOpWithSolveBase<Scalar> *Op,
00441   const ESupportSolveUse supportSolveUse
00442   ) const
00443 {
00444 
00445   using Teuchos::rcp;
00446   using Teuchos::set_extra_data;
00447   typedef Teuchos::ScalarTraits<Scalar> ST;
00448   typedef typename ST::magnitudeType ScalarMag;
00449   typedef MultiVectorBase<Scalar> MV_t;
00450   typedef LinearOpBase<Scalar> LO_t;
00451 
00452   const RCP<Teuchos::FancyOStream> out = this->getOStream();
00453   const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00454   Teuchos::OSTab tab(out);
00455   if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW))
00456     *out << "\nEntering Thyra::BelosLinearOpWithSolveFactory<"<<ST::name()<<">::initializeOpImpl(...) ...\n";
00457 
00458   // These lines are changing the verbosity of the preconditioner, which has its own verbose object list,
00459   // so I am commenting these out, as it is not the job of the linear solver to dictate preconditioner verbosity.
00460   //typedef Teuchos::VerboseObjectTempState<PreconditionerFactoryBase<Scalar> > VOTSPF;
00461   //VOTSPF precFactoryOutputTempState(precFactory_,out,verbLevel);
00462   
00463   TEST_FOR_EXCEPT(Op==NULL);
00464   TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
00465   TEST_FOR_EXCEPT(fwdOpSrc->getOp().get()==NULL);
00466   RCP<const LinearOpBase<Scalar> >
00467     fwdOp = fwdOpSrc->getOp(),
00468     approxFwdOp = ( approxFwdOpSrc.get() ? approxFwdOpSrc->getOp() : Teuchos::null );
00469 
00470   //
00471   // Get the BelosLinearOpWithSolve interface
00472   //
00473 
00474   BelosLinearOpWithSolve<Scalar>
00475     *belosOp = &Teuchos::dyn_cast<BelosLinearOpWithSolve<Scalar> >(*Op);
00476 
00477   //
00478   // Get/Create the preconditioner
00479   //
00480 
00481   RCP<PreconditionerBase<Scalar> > myPrec = Teuchos::null;
00482   RCP<const PreconditionerBase<Scalar> > prec = Teuchos::null;
00483   if(prec_in.get()) {
00484     // Use an externally defined preconditioner
00485     prec = prec_in;
00486   }
00487   else {
00488     // Try and generate a preconditioner on our own
00489     if(precFactory_.get()) {
00490       myPrec =
00491         ( !belosOp->isExternalPrec()
00492           ? Teuchos::rcp_const_cast<PreconditionerBase<Scalar> >(belosOp->extract_prec())
00493           : Teuchos::null
00494           );
00495       bool hasExistingPrec = false;
00496       if(myPrec.get()) {
00497         hasExistingPrec = true;
00498         // ToDo: Get the forward operator and validate that it is the same
00499         // operator that is used here!
00500       }
00501       else {
00502         hasExistingPrec = false;
00503         myPrec = precFactory_->createPrec();
00504       }
00505       if( hasExistingPrec && reusePrec ) {
00506         // Just reuse the existing preconditioner again!
00507       }
00508       else {
00509         // Update the preconditioner
00510         if(approxFwdOp.get())
00511           precFactory_->initializePrec(approxFwdOpSrc,&*myPrec);
00512         else
00513           precFactory_->initializePrec(fwdOpSrc,&*myPrec);
00514       }
00515       prec = myPrec;
00516     }
00517   }
00518 
00519   //
00520   // Uninitialize the current solver object
00521   //
00522 
00523   bool oldIsExternalPrec = false;
00524   RCP<Belos::LinearProblem<Scalar,MV_t,LO_t> > oldLP = Teuchos::null;
00525   RCP<Belos::SolverManager<Scalar,MV_t,LO_t> > oldIterSolver = Teuchos::null;
00526   RCP<const LinearOpSourceBase<Scalar> > oldFwdOpSrc = Teuchos::null;
00527   RCP<const LinearOpSourceBase<Scalar> > oldApproxFwdOpSrc = Teuchos::null;   
00528   ESupportSolveUse oldSupportSolveUse = SUPPORT_SOLVE_UNSPECIFIED;
00529 
00530   belosOp->uninitialize( &oldLP, NULL, &oldIterSolver, &oldFwdOpSrc,
00531     NULL, &oldIsExternalPrec, &oldApproxFwdOpSrc, &oldSupportSolveUse );
00532 
00533   //
00534   // Create the Belos linear problem
00535   // NOTE:  If one exists already, reuse it.
00536   //
00537 
00538   typedef Belos::LinearProblem<Scalar,MV_t,LO_t> LP_t;
00539   RCP<LP_t> lp;
00540   if (oldLP != Teuchos::null) {
00541     lp = oldLP;
00542   }
00543   else {
00544     lp = rcp(new LP_t());
00545   }
00546 
00547   //
00548   // Set the operator
00549   //
00550 
00551   lp->setOperator(fwdOp);
00552 
00553   //
00554   // Set the preconditioner
00555   //
00556 
00557   if(prec.get()) {
00558     RCP<const LinearOpBase<Scalar> > unspecified = prec->getUnspecifiedPrecOp();
00559     RCP<const LinearOpBase<Scalar> > left = prec->getLeftPrecOp();
00560     RCP<const LinearOpBase<Scalar> > right = prec->getRightPrecOp();
00561     TEST_FOR_EXCEPTION(
00562       !( left.get() || right.get() || unspecified.get() ), std::logic_error
00563       ,"Error, at least one preconditoner linear operator objects must be set!"
00564       );
00565     if(unspecified.get()) {
00566       lp->setRightPrec(unspecified);
00567       // ToDo: Allow user to determine whether this should be placed on the
00568       // left or on the right through a parameter in the parameter list!
00569     }
00570     else {
00571       // Set a left, right or split preconditioner
00572       TEST_FOR_EXCEPTION(
00573         left.get(),std::logic_error
00574         ,"Error, we can not currently handle a left preconditioner!"
00575         );
00576       lp->setRightPrec(right);
00577     }
00578   }
00579   if(myPrec.get()) {
00580     set_extra_data<RCP<PreconditionerBase<Scalar> > >(myPrec,"Belos::InternalPrec",
00581       Teuchos::inOutArg(lp), Teuchos::POST_DESTROY, false);
00582   }
00583   else if(prec.get()) {
00584     set_extra_data<RCP<const PreconditionerBase<Scalar> > >(prec,"Belos::ExternalPrec",
00585       Teuchos::inOutArg(lp), Teuchos::POST_DESTROY, false);
00586   }
00587 
00588   //
00589   // Generate the parameter list.
00590   //
00591 
00592   typedef Belos::SolverManager<Scalar,MV_t,LO_t> IterativeSolver_t;
00593   RCP<IterativeSolver_t> iterativeSolver = Teuchos::null;
00594   RCP<Teuchos::ParameterList> solverPL = Teuchos::rcp( new Teuchos::ParameterList() );
00595   
00596   switch(solverType_) {
00597     case SOLVER_TYPE_BLOCK_GMRES: 
00598     {
00599       // Set the PL
00600       if(paramList_.get()) {
00601         Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
00602         Teuchos::ParameterList &gmresPL = solverTypesPL.sublist(BlockGMRES_name);
00603         solverPL = Teuchos::rcp( &gmresPL, false );
00604       }
00605       // Create the solver
00606       if (oldIterSolver != Teuchos::null) {
00607         iterativeSolver = oldIterSolver;
00608         iterativeSolver->setProblem( lp );
00609         iterativeSolver->setParameters( solverPL );
00610       } 
00611       else {
00612         iterativeSolver = rcp(new Belos::BlockGmresSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
00613       }
00614       break;
00615     }
00616     case SOLVER_TYPE_PSEUDO_BLOCK_GMRES:
00617     {
00618       // Set the PL
00619       if(paramList_.get()) {
00620         Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
00621         Teuchos::ParameterList &pbgmresPL = solverTypesPL.sublist(PseudoBlockGMRES_name);
00622         solverPL = Teuchos::rcp( &pbgmresPL, false );
00623       }
00624       // 
00625       // Create the solver
00626       // 
00627       if (oldIterSolver != Teuchos::null) {
00628         iterativeSolver = oldIterSolver;
00629         iterativeSolver->setProblem( lp );
00630         iterativeSolver->setParameters( solverPL );
00631       }
00632       else {
00633         iterativeSolver = rcp(new Belos::PseudoBlockGmresSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
00634       }
00635       break;
00636     }
00637     case SOLVER_TYPE_BLOCK_CG:
00638     {
00639       // Set the PL
00640       if(paramList_.get()) {
00641         Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
00642         Teuchos::ParameterList &cgPL = solverTypesPL.sublist(BlockCG_name);
00643         solverPL = Teuchos::rcp( &cgPL, false );
00644       }
00645       // Create the solver
00646       if (oldIterSolver != Teuchos::null) {
00647         iterativeSolver = oldIterSolver;
00648         iterativeSolver->setProblem( lp );
00649         iterativeSolver->setParameters( solverPL );
00650       }
00651       else {
00652         iterativeSolver = rcp(new Belos::BlockCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
00653       }
00654       break;
00655     }
00656     case SOLVER_TYPE_PSEUDO_BLOCK_CG:
00657     {
00658       // Set the PL
00659       if(paramList_.get()) {
00660         Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
00661         Teuchos::ParameterList &pbcgPL = solverTypesPL.sublist(PseudoBlockCG_name);
00662         solverPL = Teuchos::rcp( &pbcgPL, false );
00663       }
00664       // 
00665       // Create the solver
00666       // 
00667       if (oldIterSolver != Teuchos::null) {
00668         iterativeSolver = oldIterSolver;
00669         iterativeSolver->setProblem( lp );
00670         iterativeSolver->setParameters( solverPL );
00671       }
00672       else {
00673         iterativeSolver = rcp(new Belos::PseudoBlockCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
00674       }
00675       break;
00676     }
00677     case SOLVER_TYPE_GCRODR:
00678     {
00679       // Set the PL
00680       if(paramList_.get()) {
00681         Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
00682         Teuchos::ParameterList &gcrodrPL = solverTypesPL.sublist(GCRODR_name);
00683         solverPL = Teuchos::rcp( &gcrodrPL, false );
00684       }
00685       // Create the solver
00686       if (oldIterSolver != Teuchos::null) {
00687         iterativeSolver = oldIterSolver;
00688         iterativeSolver->setProblem( lp );
00689         iterativeSolver->setParameters( solverPL );
00690       } 
00691       else {
00692         iterativeSolver = rcp(new Belos::GCRODRSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
00693       }
00694       break;
00695     }
00696     case SOLVER_TYPE_RCG:
00697     {
00698       // Set the PL
00699       if(paramList_.get()) {
00700         Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
00701         Teuchos::ParameterList &rcgPL = solverTypesPL.sublist(RCG_name);
00702         solverPL = Teuchos::rcp( &rcgPL, false );
00703       }
00704       // Create the solver
00705       if (oldIterSolver != Teuchos::null) {
00706         iterativeSolver = oldIterSolver;
00707         iterativeSolver->setProblem( lp );
00708         iterativeSolver->setParameters( solverPL );
00709       } 
00710       else {
00711         iterativeSolver = rcp(new Belos::RCGSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
00712       }
00713       break;
00714     }
00715     case SOLVER_TYPE_MINRES:
00716     {
00717       // Set the PL
00718       if(paramList_.get()) {
00719         Teuchos::ParameterList &solverTypesPL = paramList_->sublist(SolverTypes_name);
00720         Teuchos::ParameterList &minresPL = solverTypesPL.sublist(MINRES_name);
00721         solverPL = Teuchos::rcp( &minresPL, false );
00722       }
00723       // Create the solver
00724       if (oldIterSolver != Teuchos::null) {
00725         iterativeSolver = oldIterSolver;
00726         iterativeSolver->setProblem( lp );
00727         iterativeSolver->setParameters( solverPL );
00728       }
00729       else {
00730         iterativeSolver = rcp(new Belos::MinresSolMgr<Scalar,MV_t,LO_t>(lp,solverPL));
00731       }
00732       break;
00733     }
00734 
00735     default:
00736     {
00737       TEST_FOR_EXCEPT(true);
00738     }
00739   }
00740 
00741   //
00742   // Initialize the LOWS object
00743   //
00744 
00745   belosOp->initialize(
00746     lp, solverPL, iterativeSolver,
00747     fwdOpSrc, prec, myPrec.get()==NULL, approxFwdOpSrc,
00748     supportSolveUse, convergenceTestFrequency_
00749     );
00750   belosOp->setOStream(out);
00751   belosOp->setVerbLevel(verbLevel);
00752 #ifdef TEUCHOS_DEBUG
00753   if(paramList_.get()) {
00754     // Make sure we read the list correctly
00755     paramList_->validateParameters(*this->getValidParameters(),1); // Validate 0th and 1st level deep
00756   }
00757 #endif
00758   if(out.get() && static_cast<int>(verbLevel) > static_cast<int>(Teuchos::VERB_LOW))
00759     *out << "\nLeaving Thyra::BelosLinearOpWithSolveFactory<"<<ST::name()<<">::initializeOpImpl(...) ...\n";
00760   
00761 }
00762 
00763 
00764 } // namespace Thyra
00765 
00766 
00767 #endif // THYRA_BELOS_LINEAR_OP_WITH_SOLVE_FACTORY_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Friends