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