Rythmos - Transient Integration for Differential Equations Version of the Day
Rythmos_ImplicitBDFStepper_def.hpp
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //                           Rythmos Package
00005 //                 Copyright (2006) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 //@HEADER
00028 
00029 #ifndef Rythmos_IMPLICITBDF_STEPPER_DEF_H
00030 #define Rythmos_IMPLICITBDF_STEPPER_DEF_H
00031 
00032 #include "Rythmos_ImplicitBDFStepper_decl.hpp"
00033 #include "Rythmos_StepperHelpers.hpp"
00034 #include "Rythmos_ImplicitBDFStepperStepControl.hpp"
00035 
00036 namespace Rythmos {
00037 
00038 // ////////////////////////////
00039 // Defintions
00040 
00041 // Nonmember constructor
00042 template<class Scalar>
00043 RCP<ImplicitBDFStepper<Scalar> > implicitBDFStepper() {
00044   RCP<ImplicitBDFStepper<Scalar> > stepper = rcp(new ImplicitBDFStepper<Scalar>() );
00045   return stepper;
00046 }
00047 
00048 template<class Scalar>
00049 RCP<ImplicitBDFStepper<Scalar> > implicitBDFStepper(
00050   const RCP<Thyra::ModelEvaluator<Scalar> >& model,
00051   const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver,
00052   const RCP<Teuchos::ParameterList>& parameterList
00053   )
00054 {
00055   RCP<ImplicitBDFStepper<Scalar> > stepper = Teuchos::rcp(new ImplicitBDFStepper<Scalar>(model,solver,parameterList));
00056   return stepper;
00057 }
00058 
00059 // Constructors, intializers, Misc.
00060 
00061 
00062 template<class Scalar>
00063 ImplicitBDFStepper<Scalar>::ImplicitBDFStepper()
00064 {
00065   this->defaultInitializeAll_();
00066   haveInitialCondition_ = false;
00067   isInitialized_=false;
00068 }
00069 
00070 
00071 template<class Scalar>
00072 ImplicitBDFStepper<Scalar>::ImplicitBDFStepper(
00073   const RCP<Thyra::ModelEvaluator<Scalar> >& model
00074   ,const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
00075   ,const RCP<Teuchos::ParameterList>& parameterList
00076   )
00077 {
00078   this->defaultInitializeAll_();
00079   this->setParameterList(parameterList);
00080   // Now we instantiate the model and the solver
00081   setModel(model);
00082   setSolver(solver);
00083   haveInitialCondition_ = false;
00084   isInitialized_=false;
00085 }
00086 
00087 
00088 template<class Scalar>
00089 ImplicitBDFStepper<Scalar>::ImplicitBDFStepper(
00090   const RCP<Thyra::ModelEvaluator<Scalar> >& model
00091   ,const RCP<Thyra::NonlinearSolverBase<Scalar> >& solver
00092   )
00093 {
00094   this->defaultInitializeAll_();
00095   // Now we instantiate the model and the solver
00096   setModel(model);
00097   setSolver(solver);
00098   haveInitialCondition_ = false;
00099   isInitialized_=false;
00100 }
00101 
00102 template<class Scalar>
00103 const Thyra::VectorBase<Scalar>& 
00104   ImplicitBDFStepper<Scalar>::getxHistory(int index) const
00105 {
00106   TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,
00107       "Error, attempting to call getxHistory before initialization!\n");
00108   TEUCHOS_TEST_FOR_EXCEPT( !(( 0 <= index ) && ( index <= maxOrder_ )) );
00109   TEUCHOS_TEST_FOR_EXCEPT( !( index <= usedOrder_+1 ) );
00110   return(*(xHistory_[index]));
00111 }
00112 
00113 template<class Scalar>
00114 void ImplicitBDFStepper<Scalar>::setStepControlStrategy(const RCP<StepControlStrategyBase<Scalar> >& stepControl)
00115 {
00116   TEUCHOS_TEST_FOR_EXCEPTION(stepControl == Teuchos::null,std::logic_error,"Error, stepControl == Teuchos::null!\n");
00117   stepControl_ = stepControl;    
00118 }
00119 
00120 template<class Scalar>
00121 RCP<StepControlStrategyBase<Scalar> > ImplicitBDFStepper<Scalar>::getNonconstStepControlStrategy() 
00122 {
00123   return(stepControl_);
00124 }
00125 
00126 template<class Scalar>
00127 RCP<const StepControlStrategyBase<Scalar> > ImplicitBDFStepper<Scalar>::getStepControlStrategy() const
00128 {
00129   return(stepControl_);
00130 }
00131 
00132 
00133 // Overridden from SolverAcceptingStepperBase
00134 
00135 
00136 template<class Scalar>
00137 void ImplicitBDFStepper<Scalar>::setSolver(const RCP<Thyra::NonlinearSolverBase<Scalar> > &solver)
00138 {
00139   TEUCHOS_TEST_FOR_EXCEPT(solver == Teuchos::null)
00140     solver_ = solver;
00141 }
00142 
00143 
00144 template<class Scalar>
00145 RCP<Thyra::NonlinearSolverBase<Scalar> >
00146 ImplicitBDFStepper<Scalar>::getNonconstSolver()
00147 {
00148   return (solver_);
00149 }
00150 
00151 
00152 template<class Scalar>
00153 RCP<const Thyra::NonlinearSolverBase<Scalar> >
00154 ImplicitBDFStepper<Scalar>::getSolver() const
00155 {
00156   return (solver_);
00157 }
00158 
00159 
00160 // Overridden from StepperBase
00161 
00162 
00163 template<class Scalar>
00164 bool ImplicitBDFStepper<Scalar>::isImplicit() const
00165 {
00166   return true;
00167 }
00168 
00169 template<class Scalar>
00170 bool ImplicitBDFStepper<Scalar>::supportsCloning() const
00171 {
00172   return true;
00173 }
00174 
00175 
00176 template<class Scalar>
00177 RCP<StepperBase<Scalar> >
00178 ImplicitBDFStepper<Scalar>::cloneStepperAlgorithm() const
00179 {
00180 
00181   // Just use the interface to clone the algorithm in an basically
00182   // uninitialized state
00183 
00184   RCP<ImplicitBDFStepper<Scalar> >
00185     stepper = Teuchos::rcp(new ImplicitBDFStepper<Scalar>());
00186 
00187   if (!is_null(model_))
00188     stepper->setModel(model_); // Shallow copy is okay!
00189 
00190   if (!is_null(solver_))
00191     stepper->setSolver(solver_->cloneNonlinearSolver().assert_not_null());
00192   
00193   if (!is_null(parameterList_))
00194     stepper->setParameterList(Teuchos::parameterList(*parameterList_));
00195 
00196   if (!is_null(stepControl_)) {
00197     if (stepControl_->supportsCloning())
00198       stepper->setStepControlStrategy(
00199         stepControl_->cloneStepControlStrategyAlgorithm().assert_not_null());
00200   }
00201   
00202   // At this point, we should have a valid algorithm state.  What might be
00203   // missing is the initial condition if it was not given in *model_ but was
00204   // set explicitly.  However, the specification for this function does not
00205   // guarantee that the full state will be copied in any case!
00206 
00207   return stepper;
00208 
00209 }
00210 
00211 
00212 template<class Scalar>
00213 void ImplicitBDFStepper<Scalar>::setModel(
00214   const RCP<const Thyra::ModelEvaluator<Scalar> >& model
00215   )
00216 {
00217   typedef Teuchos::ScalarTraits<Scalar> ST;
00218   TEUCHOS_TEST_FOR_EXCEPT( is_null(model) );
00219   assertValidModel( *this, *model );
00220   model_ = model;
00221 }
00222 
00223 template<class Scalar>
00224 void ImplicitBDFStepper<Scalar>::setNonconstModel(
00225   const RCP<Thyra::ModelEvaluator<Scalar> >& model
00226   )
00227 {
00228   this->setModel(model); // TODO 09/09/09 tscoffe:  use ConstNonconstObjectContainer!
00229 }
00230 
00231 
00232 template<class Scalar>
00233 RCP<const Thyra::ModelEvaluator<Scalar> >
00234 ImplicitBDFStepper<Scalar>::getModel() const
00235 {
00236   return model_;
00237 }
00238 
00239 
00240 template<class Scalar>
00241 RCP<Thyra::ModelEvaluator<Scalar> >
00242 ImplicitBDFStepper<Scalar>::getNonconstModel()
00243 {
00244   return Teuchos::null;
00245 }
00246 
00247 
00248 template<class Scalar>
00249 void ImplicitBDFStepper<Scalar>::setInitialCondition(
00250   const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
00251   )
00252 {
00253   typedef Teuchos::ScalarTraits<Scalar> ST;
00254   typedef Thyra::ModelEvaluatorBase MEB;
00255   TEUCHOS_TEST_FOR_EXCEPT(is_null(initialCondition.get_x()));
00256   TEUCHOS_TEST_FOR_EXCEPT(is_null(initialCondition.get_x_dot()));
00257   basePoint_ = initialCondition;
00258   xn0_ = initialCondition.get_x()->clone_v();
00259   xpn0_ = initialCondition.get_x_dot()->clone_v(); 
00260   time_ = initialCondition.get_t();
00261   // Generate vectors for use in calculations
00262   x_dot_base_ = createMember(xpn0_->space());
00263   V_S(x_dot_base_.ptr(),ST::zero());
00264   ee_ = createMember(xn0_->space());
00265   V_S(ee_.ptr(),ST::zero());
00266   // x history
00267   xHistory_.clear();
00268   xHistory_.push_back(xn0_->clone_v());
00269   xHistory_.push_back(xpn0_->clone_v());
00270 
00271   haveInitialCondition_ = true;
00272   if (isInitialized_) {
00273     initialize_();
00274   }
00275 }
00276 
00277 
00278 template<class Scalar>
00279 Thyra::ModelEvaluatorBase::InArgs<Scalar> 
00280 ImplicitBDFStepper<Scalar>::getInitialCondition() const
00281 {
00282   return basePoint_;
00283 }
00284 
00285 
00286 template<class Scalar>
00287 Scalar ImplicitBDFStepper<Scalar>::takeStep(Scalar dt, StepSizeType stepType)
00288 {
00289   
00290   RYTHMOS_FUNC_TIME_MONITOR("Rythmos::ImplicitBDFStepper::takeStep");
00291   
00292   using Teuchos::as;
00293   using Teuchos::incrVerbLevel;
00294   typedef Teuchos::ScalarTraits<Scalar> ST;
00295   typedef typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag TScalarMag;
00296   typedef Thyra::NonlinearSolverBase<Scalar> NSB;
00297   typedef Teuchos::VerboseObjectTempState<NSB> VOTSNSB;
00298 
00299   RCP<Teuchos::FancyOStream> out = this->getOStream();
00300   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00301   Teuchos::OSTab ostab(out,1,"takeStep");
00302   VOTSNSB solver_outputTempState(solver_,out,incrVerbLevel(verbLevel,-1));
00303 
00304   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00305     *out
00306       << "\nEntering " << this->Teuchos::Describable::description()
00307       << "::takeStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
00308   }
00309 
00310   if (!isInitialized_) {
00311     initialize_(); 
00312   }
00313 
00314   stepControl_->setOStream(out);
00315   stepControl_->setVerbLevel(verbLevel);
00316   stepControl_->setRequestedStepSize(*this,dt,stepType);
00317 
00318   AttemptedStepStatusFlag status;
00319   while (1) {
00320     // Set up problem coefficients (and handle first step)
00321     Scalar hh_old = hh_;
00322     int desiredOrder;
00323     stepControl_->nextStepSize(*this,&hh_,&stepType,&desiredOrder);
00324     TEUCHOS_TEST_FOR_EXCEPT(!((1 <= desiredOrder) && (desiredOrder <= maxOrder_)));
00325     TEUCHOS_TEST_FOR_EXCEPT(!(desiredOrder <= usedOrder_+1));
00326     currentOrder_ = desiredOrder;
00327     if (numberOfSteps_ == 0) {
00328       psi_[0] = hh_;
00329       if (nef_ == 0) {
00330         Vt_S(xHistory_[1].ptr(),hh_);
00331       } else {
00332         Vt_S(xHistory_[1].ptr(),hh_/hh_old);
00333       }
00334     }
00335     this->updateCoeffs_();
00336     // compute predictor
00337     obtainPredictor_();
00338     // solve nonlinear problem (as follows)
00339     
00340     //
00341     // Setup the nonlinear equations:
00342     //
00343     //   f_bar( x_dot_coeff * x_bar + x_dot_base, x_coeff * x_bar + x_base, t_base ) = 0
00344     //   x_dot_coeff = -alpha_s/dt
00345     //   x_dot_base = x_prime_pred + (alpha_s/dt) * x_pred
00346     //   x_coeff = 1
00347     //   x_base = 0
00348     //   t_base = tn+dt
00349     //
00350     Scalar coeff_x_dot = Scalar(-ST::one())*alpha_s_/hh_;
00351     V_StVpStV( x_dot_base_.ptr(), ST::one(), *xpn0_, alpha_s_/hh_, *xn0_ );
00352     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00353       *out << "model_ = " << std::endl;
00354       model_->describe(*out,verbLevel);
00355       *out << "basePoint_ = " << std::endl;
00356       basePoint_.describe(*out,verbLevel);
00357       *out << "coeff_x_dot = " << coeff_x_dot << std::endl;
00358       *out << "x_dot_base_ = " << std::endl;
00359       x_dot_base_->describe(*out,verbLevel);
00360       *out << "time_+hh_ = " << time_+hh_ << std::endl;
00361       *out << "xn0_ = " << std::endl;
00362       xn0_->describe(*out,verbLevel);
00363     }
00364     neModel_.initializeSingleResidualModel(
00365       model_, basePoint_,
00366       coeff_x_dot, x_dot_base_,
00367       ST::one(), Teuchos::null,
00368       time_+hh_,
00369       xn0_
00370       );
00371     //
00372     // Solve the implicit nonlinear system to a tolerance of ???
00373     // 
00374     if (solver_->getModel().get() != &neModel_) {
00375       solver_->setModel( Teuchos::rcpFromRef(neModel_) );
00376     }
00377     /* // Rythmos::TimeStepNonlinearSolver uses a built in solveCriteria, so you can't pass one in.
00378     // I believe this is the correct solveCriteria for IDA though.
00379     Thyra::SolveMeasureType nonlinear_solve_measure_type(Thyra::SOLVE_MEASURE_NORM_RESIDUAL,Thyra::SOLVE_MEASURE_ONE); 
00380     TScalarMag tolerance = relErrTol_/TScalarMag(10.0); // This should be changed to match the condition in IDA
00381     Thyra::SolveCriteria<Scalar> nonlinearSolveCriteria(nonlinear_solve_measure_type, tolerance);
00382     Thyra::SolveStatus<Scalar> nonlinearSolveStatus = solver_->solve( &*xn0_, &nonlinearSolveCriteria, &*delta_ ); 
00383     */
00384     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00385       *out << "xn0 = " << std::endl;
00386       xn0_->describe(*out,verbLevel);
00387       *out << "ee = " << std::endl;
00388       ee_->describe(*out,verbLevel);
00389     }
00390     Thyra::SolveStatus<Scalar> nonlinearSolveStatus = solver_->solve( &*xn0_, NULL, &*ee_ ); 
00391     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_EXTREME) ) {
00392       *out << "xn0 = " << std::endl;
00393       xn0_->describe(*out,verbLevel);
00394       *out << "ee = " << std::endl;
00395       ee_->describe(*out,verbLevel);
00396     }
00397     // In the above solve, on input *xn0_ is the initial guess that comes from
00398     // the predictor.  On output, *xn0_ is the solved for solution and *ee_ is
00399     // the difference computed from the intial guess in *xn0_ to the final
00400     // solved value of *xn0_.  This is needed for basic numerical stability.
00401     if (nonlinearSolveStatus.solveStatus == Thyra::SOLVE_STATUS_CONVERGED)  {
00402       newtonConvergenceStatus_ = 0;
00403     }
00404     else {
00405       newtonConvergenceStatus_ = -1;
00406     }
00407 
00408     // check error and evaluate LTE
00409     stepControl_->setCorrection(*this,xn0_,ee_,newtonConvergenceStatus_);
00410     bool stepPass = stepControl_->acceptStep(*this,&LETvalue_);
00411     
00412     if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00413       *out << "xn0_ = " << std::endl;
00414       xn0_->describe(*out,verbLevel);
00415       *out << "xpn0_ = " << std::endl;
00416       xpn0_->describe(*out,verbLevel);
00417       *out << "ee_ = " << std::endl;
00418       ee_->describe(*out,verbLevel);
00419       for (int i=0; i<std::max(2,maxOrder_); ++i) {
00420         *out << "xHistory_[" << i << "] = " << std::endl;
00421         xHistory_[i]->describe(*out,verbLevel);
00422       }
00423     }
00424 
00425     // Check LTE here:
00426     if (!stepPass) { // stepPass = false
00427       stepLETStatus_ = STEP_LET_STATUS_FAILED;
00428       status = stepControl_->rejectStep(*this);
00429       nef_++;
00430       if (status == CONTINUE_ANYWAY) {
00431         break;
00432       } else {
00433         restoreHistory_();
00434       }
00435     } else { // stepPass = true
00436       stepLETStatus_ = STEP_LET_STATUS_PASSED;
00437       break;
00438     }
00439   }
00440 
00441   // 08/22/07 the history array must be updated before stepControl_->completeStep.
00442   completeStep_();   
00443   stepControl_->completeStep(*this);
00444 
00445   if ( !is_null(out) && as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW) ) {
00446     *out
00447       << "\nLeaving " << this->Teuchos::Describable::description()
00448       << "::takeStep("<<dt<<","<<toString(stepType)<<") ...\n"; 
00449   }
00450 
00451   return(usedStep_);
00452 
00453 }
00454 
00455 
00456 template<class Scalar>
00457 const StepStatus<Scalar> ImplicitBDFStepper<Scalar>::getStepStatus() const
00458 {
00459 
00460   // 2007/08/24: rabartl: We agreed that getStepStatus() would be free
00461   // so I have commented out removed all code that is not free
00462 
00463   typedef Teuchos::ScalarTraits<Scalar> ST;
00464   StepStatus<Scalar> stepStatus;
00465   if (!isInitialized_) {
00466     stepStatus.message = "This stepper is uninitialized.";
00467     stepStatus.stepStatus = STEP_STATUS_UNINITIALIZED;
00468     stepStatus.stepSize = Scalar(-ST::one());
00469     stepStatus.order = -1;
00470     stepStatus.time = Scalar(-ST::one());
00471     return(stepStatus);
00472   }
00473 
00474   if (numberOfSteps_ > 0) {
00475     stepStatus.stepStatus = STEP_STATUS_CONVERGED; 
00476   } else {
00477     stepStatus.stepStatus = STEP_STATUS_UNKNOWN;
00478   }
00479   stepStatus.stepLETStatus = stepLETStatus_;
00480   stepStatus.stepSize = usedStep_; 
00481   stepStatus.order = usedOrder_;
00482   stepStatus.time = time_;
00483   stepStatus.stepLETValue = LETvalue_; 
00484   stepStatus.solution = xHistory_[0];
00485   stepStatus.solutionDot = Teuchos::null; // This is *not* free!
00486   stepStatus.residual = Teuchos::null; // This is *not* free!
00487 
00488   return(stepStatus);
00489 
00490 }
00491 
00492 
00493 // Overridden from InterpolationBufferBase
00494 
00495 
00496 template<class Scalar>
00497 RCP<const Thyra::VectorSpaceBase<Scalar> >
00498 ImplicitBDFStepper<Scalar>::get_x_space() const
00499 {
00500   //TEUCHOS_TEST_FOR_EXCEPTION(!isInitialized_,std::logic_error,"Error, attempting to call get_x_space before initialization!\n");
00501   return ( !is_null(model_) ? model_->get_x_space() : Teuchos::null );
00502 }
00503 
00504 
00505 template<class Scalar>
00506 void ImplicitBDFStepper<Scalar>::addPoints(
00507   const Array<Scalar>& time_vec,
00508   const Array<RCP<const Thyra::VectorBase<Scalar> > >& x_vec,
00509   const Array<RCP<const Thyra::VectorBase<Scalar> > >& xdot_vec
00510   )
00511 {
00512   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
00513     "Error, addPoints is not implemented for ImplicitBDFStepper.\n");
00514 }
00515 
00516 
00517 template<class Scalar>
00518 TimeRange<Scalar> ImplicitBDFStepper<Scalar>::getTimeRange() const
00519 {
00520   if ( !isInitialized_ && haveInitialCondition_ )
00521     return timeRange<Scalar>(time_,time_);
00522   if ( !isInitialized_ && !haveInitialCondition_ )
00523     return invalidTimeRange<Scalar>();
00524   return timeRange<Scalar>(time_-usedStep_,time_);
00525 }
00526 
00527 
00528 template<class Scalar>
00529 void ImplicitBDFStepper<Scalar>::getPoints(
00530   const Array<Scalar>& time_vec
00531   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* x_vec
00532   ,Array<RCP<const Thyra::VectorBase<Scalar> > >* xdot_vec
00533   ,Array<ScalarMag>* accuracy_vec) const
00534 {
00535   using Teuchos::as;
00536   using Teuchos::constOptInArg;
00537   using Teuchos::null;
00538   using Teuchos::ptr;
00539   typedef Teuchos::ScalarTraits<Scalar> ST;
00540   typedef typename ST::magnitudeType mScalarMag;
00541 
00542   TEUCHOS_ASSERT(haveInitialCondition_);
00543   // Only do this if we're being called pre-initialization to get the IC.
00544   if ( (numberOfSteps_ == -1) && 
00545        (time_vec.length() == 1) &&
00546        (compareTimeValues<Scalar>(time_vec[0],time_)==0) ) {
00547     defaultGetPoints<Scalar>(
00548         time_, constOptInArg(*xn0_), constOptInArg(*xpn0_),
00549         time_, constOptInArg(*xn0_), constOptInArg(*xpn0_),
00550         time_vec, ptr(x_vec), ptr(xdot_vec), ptr(accuracy_vec),
00551         Ptr<InterpolatorBase<Scalar> >(null)
00552         );
00553     return;
00554   }
00555   TEUCHOS_ASSERT(isInitialized_);
00556   RYTHMOS_FUNC_TIME_MONITOR("Rythmos::ImplicitBDFStepper::getPoints");
00557   if (x_vec)
00558     x_vec->clear();
00559   if (xdot_vec)
00560     xdot_vec->clear();
00561   for (Teuchos::Ordinal i=0 ; i<time_vec.size() ; ++i) {
00562     RCP<Thyra::VectorBase<Scalar> >
00563       x_temp = createMember(xn0_->space());
00564     RCP<Thyra::VectorBase<Scalar> >
00565       xdot_temp = createMember(xn0_->space());
00566     mScalarMag accuracy = -ST::zero();
00567     interpolateSolution_(
00568       time_vec[i], &*x_temp, &*xdot_temp,
00569       accuracy_vec ? &accuracy : 0
00570       );
00571     if (x_vec)
00572       x_vec->push_back(x_temp);
00573     if (xdot_vec)
00574       xdot_vec->push_back(xdot_temp);
00575     if (accuracy_vec)
00576       accuracy_vec->push_back(accuracy);
00577   }
00578   if ( as<int>(this->getVerbLevel()) >= as<int>(Teuchos::VERB_HIGH) ) {
00579     RCP<Teuchos::FancyOStream> out = this->getOStream();
00580     Teuchos::OSTab ostab(out,1,"getPoints");
00581     *out << "Passing out the interpolated values:" << std::endl;
00582     for (Teuchos::Ordinal i=0; i<time_vec.size() ; ++i) {
00583       *out << "time_[" << i << "] = " << time_vec[i] << std::endl;
00584       if (x_vec) {
00585         *out << "x_vec[" << i << "] = " << std::endl;
00586         (*x_vec)[i]->describe(*out,this->getVerbLevel());
00587       }
00588       if (xdot_vec) {
00589         *out << "xdot_vec[" << i << "] = ";
00590         if ( (*xdot_vec)[i] == Teuchos::null) {
00591           *out << "Teuchos::null" << std::endl;
00592         }
00593         else {
00594           *out << std::endl << Teuchos::describe(*(*xdot_vec)[i],this->getVerbLevel());
00595         }
00596       }
00597       if (accuracy_vec)
00598         *out << "accuracy[" << i << "] = " << (*accuracy_vec)[i] << std::endl;
00599     }
00600   }
00601 }
00602 
00603 
00604 template<class Scalar>
00605 void ImplicitBDFStepper<Scalar>::getNodes(Array<Scalar>* time_vec) const
00606 {
00607   TEUCHOS_ASSERT( time_vec != NULL );
00608   time_vec->clear();
00609   if (!haveInitialCondition_) {
00610     return;
00611   }
00612   if (numberOfSteps_ > 0) {
00613     time_vec->push_back(time_-usedStep_);
00614   }
00615   time_vec->push_back(time_);
00616 }
00617 
00618 
00619 template<class Scalar>
00620 void ImplicitBDFStepper<Scalar>::removeNodes(Array<Scalar>& time_vec) 
00621 {
00622   TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
00623     "Error, removeNodes is not implemented for ImplicitBDFStepper.\n");
00624 }
00625 
00626 
00627 template<class Scalar>
00628 int ImplicitBDFStepper<Scalar>::getOrder() const
00629 {
00630   if (!isInitialized_) {
00631     return(-1);
00632   }
00633   return(usedOrder_);
00634 }
00635 
00636 
00637 // Overridden from Teuchos::ParameterListAcceptor
00638 
00639 
00640 template<class Scalar>
00641 void ImplicitBDFStepper<Scalar>::setParameterList(
00642   RCP<Teuchos::ParameterList> const& paramList
00643   )
00644 {
00645   TEUCHOS_TEST_FOR_EXCEPT(paramList == Teuchos::null);
00646   paramList->validateParameters(*this->getValidParameters(),0);
00647   parameterList_ = paramList;
00648   Teuchos::readVerboseObjectSublist(&*parameterList_,this);
00649 }
00650 
00651 
00652 template<class Scalar>
00653 RCP<Teuchos::ParameterList> ImplicitBDFStepper<Scalar>::getNonconstParameterList()
00654 {
00655   return(parameterList_);
00656 }
00657 
00658 
00659 template<class Scalar>
00660 RCP<Teuchos::ParameterList>
00661 ImplicitBDFStepper<Scalar>::unsetParameterList()
00662 {
00663   RCP<Teuchos::ParameterList> temp_param_list = parameterList_;
00664   parameterList_ = Teuchos::null;
00665   return(temp_param_list);
00666 }
00667 
00668 
00669 template<class Scalar>
00670 RCP<const Teuchos::ParameterList>
00671 ImplicitBDFStepper<Scalar>::getValidParameters() const
00672 {
00673 
00674   static RCP<Teuchos::ParameterList> validPL;
00675 
00676   if (is_null(validPL)) {
00677 
00678     RCP<Teuchos::ParameterList>
00679       pl = Teuchos::parameterList();
00680 
00681     pl->sublist(RythmosStepControlSettings_name);
00682 
00683     Teuchos::setupVerboseObjectSublist(&*pl);
00684 
00685     validPL = pl;
00686 
00687   }
00688 
00689   RCP<Teuchos::FancyOStream> out = this->getOStream();
00690   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00691   Teuchos::OSTab ostab(out,1,"getValidParameters");
00692   if (Teuchos::as<int>(verbLevel) == Teuchos::VERB_HIGH) {
00693     *out << "Setting up valid parameterlist." << std::endl;
00694     validPL->print(*out);
00695   }
00696 
00697   return (validPL);
00698   
00699 }
00700 
00701 
00702 // Overridden from Teuchos::Describable
00703 
00704 
00705 template<class Scalar>
00706 std::string ImplicitBDFStepper<Scalar>::description() const
00707 {
00708   std::ostringstream out;
00709   out << this->Teuchos::Describable::description();
00710   const TimeRange<Scalar> timeRange = this->getTimeRange();
00711   if (timeRange.isValid())
00712     out << " (timeRange="<<timeRange<<")";
00713   else
00714     out << " (This stepper is not initialized yet)";
00715   out << std::endl;
00716   return out.str();
00717 }
00718 
00719 
00720 template<class Scalar>
00721 void ImplicitBDFStepper<Scalar>::describe(
00722   Teuchos::FancyOStream &out,
00723   const Teuchos::EVerbosityLevel verbLevel
00724   ) const
00725 {
00726 
00727   using Teuchos::as;
00728 
00729   if (!isInitialized_) {
00730     out << this->description();
00731     return;
00732   }
00733   
00734   if ( (as<int>(verbLevel) == as<int>(Teuchos::VERB_DEFAULT) ) ||
00735     (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)     )
00736     )
00737   {
00738     out << this->description() << std::endl;
00739     out << "model_ = " << Teuchos::describe(*model_,verbLevel);
00740     out << "solver_ = " << Teuchos::describe(*solver_,verbLevel);
00741     out << "neModel_ = " << Teuchos::describe(neModel_,verbLevel);
00742   }
00743   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_LOW)) {
00744     out << "time_ = " << time_ << std::endl;
00745     out << "hh_ = " << hh_ << std::endl;
00746     out << "currentOrder_ = " << currentOrder_ << std::endl;
00747   }
00748   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH)) {
00749     out << "xn0_ = " << Teuchos::describe(*xn0_,verbLevel);
00750     out << "xpn0_ = " << Teuchos::describe(*xpn0_,verbLevel);
00751     out << "x_dot_base_ = " << Teuchos::describe(*x_dot_base_,verbLevel);
00752     for (int i=0 ; i < std::max(2,maxOrder_) ; ++i) {
00753       out << "xHistory_[" << i << "] = "
00754           << Teuchos::describe(*xHistory_[i],verbLevel);
00755     }
00756     out << "ee_ = " << Teuchos::describe(*ee_,verbLevel);
00757   }
00758 }
00759 
00760 
00761 // private
00762 
00763 
00764 // 2007/08/24: rabartl: Belos: We really should initialize all of this data in
00765 // a member initialization list but since there are like three constructors
00766 // this would mean that we would have to duplicate code (which is error prone)
00767 // or use a macro (which is not easy to debug).  We really should remove all
00768 // but the default constructor which then would set this data once in the
00769 // initialization list.
00770 
00771 template<class Scalar>
00772 void ImplicitBDFStepper<Scalar>::defaultInitializeAll_()
00773 {
00774   typedef Teuchos::ScalarTraits<Scalar> ST;
00775   const Scalar nan = ST::nan(), one = ST::one(), zero = ST::zero();
00776   // Initialize some data members to their rightful default values
00777   haveInitialCondition_ = false;
00778   isInitialized_ = false;
00779   currentOrder_ = 1;
00780   usedOrder_ = 1;
00781   usedStep_ = zero;
00782   // Initialize the rest of the private data members to invalid values to
00783   // avoid uninitialed memory
00784   time_ = nan;
00785   hh_ = nan;
00786   maxOrder_ = -1;
00787   LETvalue_ = -one;
00788   stepLETStatus_ = STEP_LET_STATUS_UNKNOWN;
00789   alpha_s_ = -one;
00790   numberOfSteps_ = -1;
00791   nef_ = -1;
00792   nscsco_ = -1;
00793   newtonConvergenceStatus_ = -1;
00794 }
00795 
00796 
00797 template<class Scalar>
00798 void ImplicitBDFStepper<Scalar>::obtainPredictor_()
00799 {
00800 
00801   using Teuchos::as;
00802   typedef Teuchos::ScalarTraits<Scalar> ST;
00803 
00804   if (!isInitialized_) {
00805     return;
00806   }
00807 
00808   RCP<Teuchos::FancyOStream> out = this->getOStream();
00809   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00810   Teuchos::OSTab ostab(out,1,"obtainPredictor_");
00811   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00812     *out << "currentOrder_ = " << currentOrder_ << std::endl;
00813   }
00814   
00815   // prepare history array for prediction
00816   for (int i=nscsco_;i<=currentOrder_;++i) {
00817     Vt_S(xHistory_[i].ptr(),beta_[i]);
00818   }
00819   
00820   // evaluate predictor
00821   V_V(xn0_.ptr(),*xHistory_[0]);
00822   V_S(xpn0_.ptr(),ST::zero());
00823   for (int i=1;i<=currentOrder_;++i) {
00824     Vp_V(xn0_.ptr(),*xHistory_[i]);
00825     Vp_StV(xpn0_.ptr(),gamma_[i],*xHistory_[i]);
00826   }
00827   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00828     *out << "xn0_ = " << std::endl;
00829     xn0_->describe(*out,verbLevel);
00830     *out << "xpn0_ = " << std::endl;
00831     xpn0_->describe(*out,verbLevel);
00832   }
00833 }
00834 
00835 
00836 template<class Scalar>
00837 void ImplicitBDFStepper<Scalar>::interpolateSolution_(
00838   const Scalar& timepoint,
00839   Thyra::VectorBase<Scalar>* x_ptr,
00840   Thyra::VectorBase<Scalar>* xdot_ptr,
00841   ScalarMag* accuracy_ptr
00842   ) const
00843 {
00844 
00845   typedef std::numeric_limits<Scalar> NL;
00846   typedef Teuchos::ScalarTraits<Scalar> ST;
00847 
00848 #ifdef HAVE_RYTHMOS_DEBUG
00849   TEUCHOS_TEST_FOR_EXCEPTION(
00850     !isInitialized_,std::logic_error,
00851     "Error, attempting to call interpolateSolution before initialization!\n");
00852   const TimeRange<Scalar> currTimeRange = this->getTimeRange();
00853   TEUCHOS_TEST_FOR_EXCEPTION(
00854     !currTimeRange.isInRange(timepoint), std::logic_error,
00855     "Error, timepoint = " << timepoint << " is not in the time range "
00856     << currTimeRange << "!" );
00857 #endif
00858 
00859   const Scalar tn = time_;
00860   const int kused = usedOrder_;
00861 
00862   // order of interpolation
00863   int kord = kused;
00864   if ( (kused == 0) || (timepoint == tn) )  {
00865     kord = 1;
00866   }
00867 
00868   // Initialize interploation
00869   Thyra::V_V(ptr(x_ptr),*xHistory_[0]);
00870   Thyra::V_S(ptr(xdot_ptr),ST::zero());
00871   
00872   // Add history array contributions
00873   const Scalar delt = timepoint - tn;
00874   Scalar c = ST::one(); // coefficient for interpolation of x
00875   Scalar d = ST::zero(); // coefficient for interpolation of xdot
00876   Scalar gam = delt/psi_[0]; // coefficient for interpolation
00877   for (int j=1 ; j <= kord ; ++j) {
00878     d = d*gam + c/psi_[j-1];
00879     c = c*gam;
00880     gam = (delt + psi_[j-1])/psi_[j];
00881     Thyra::Vp_StV(ptr(x_ptr),c,*xHistory_[j]);
00882     Thyra::Vp_StV(ptr(xdot_ptr),d,*xHistory_[j]);
00883   }
00884 
00885   // Set approximate accuracy
00886   if (accuracy_ptr) {
00887     *accuracy_ptr = pow(usedStep_,kord);
00888   }
00889   
00890 }
00891 
00892 
00893 template<class Scalar>
00894 void ImplicitBDFStepper<Scalar>::updateHistory_()
00895 {
00896 
00897   using Teuchos::as;
00898 
00899   // Save Newton correction for potential order increase on next step.
00900   if (usedOrder_ < maxOrder_)  {
00901     assign( xHistory_[usedOrder_+1].ptr(), *ee_ );
00902   }
00903   // Update history arrays
00904   Vp_V( xHistory_[usedOrder_].ptr(), *ee_ );
00905   for (int j=usedOrder_-1;j>=0;j--) {
00906     Vp_V( xHistory_[j].ptr(), *xHistory_[j+1] );
00907   }
00908   RCP<Teuchos::FancyOStream> out = this->getOStream();
00909   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00910   Teuchos::OSTab ostab(out,1,"updateHistory_");
00911   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00912     for (int i=0;i<std::max(2,maxOrder_);++i) {
00913       *out << "xHistory_[" << i << "] = " << std::endl;
00914       xHistory_[i]->describe(*out,verbLevel);
00915     }
00916   }
00917 
00918 }
00919 
00920 
00921 template<class Scalar>
00922 void ImplicitBDFStepper<Scalar>::restoreHistory_()
00923 {
00924 
00925   using Teuchos::as;
00926   typedef Teuchos::ScalarTraits<Scalar> ST;
00927 
00928   // undo preparation of history array for prediction
00929   for (int i=nscsco_;i<=currentOrder_;++i) {
00930     Vt_S( xHistory_[i].ptr(), ST::one()/beta_[i] );
00931   }
00932   for (int i=1;i<=currentOrder_;++i) {
00933     psi_[i-1] = psi_[i] - hh_;
00934   }
00935   RCP<Teuchos::FancyOStream> out = this->getOStream();
00936   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00937   Teuchos::OSTab ostab(out,1,"restoreHistory_");
00938   if (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00939     for (int i=0;i<maxOrder_;++i) {
00940       *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
00941     }
00942     for (int i=0;i<maxOrder_;++i) {
00943       *out << "xHistory_[" << i << "] = " << std::endl;
00944       xHistory_[i]->describe(*out,verbLevel);
00945     }
00946   }
00947 
00948 } 
00949 
00950 
00951 template<class Scalar>
00952 void ImplicitBDFStepper<Scalar>::updateCoeffs_()
00953 {
00954 
00955   using Teuchos::as;
00956   typedef Teuchos::ScalarTraits<Scalar> ST;
00957 
00958   // If the number of steps taken with constant order and constant stepsize is
00959   // more than the current order + 1 then we don't bother to update the
00960   // coefficients because we've reached a constant step-size formula.  When
00961   // this is is not true, then we update the coefficients for the variable
00962   // step-sizes. 
00963   if ((hh_ != usedStep_) || (currentOrder_ != usedOrder_)) {
00964     nscsco_ = 0;
00965   }
00966   nscsco_ = std::min(nscsco_+1,usedOrder_+2);
00967   if (currentOrder_+1 >= nscsco_) {
00968     beta_[0] = ST::one();
00969     alpha_[0] = ST::one();
00970     Scalar temp1 = hh_;
00971     gamma_[0] = ST::zero();
00972     for (int i=1;i<=currentOrder_;++i) {
00973       Scalar temp2 = psi_[i-1];
00974       psi_[i-1] = temp1;
00975       beta_[i] = beta_[i-1]*psi_[i-1]/temp2;
00976       temp1 = temp2 + hh_;
00977       alpha_[i] = hh_/temp1;
00978       gamma_[i] = gamma_[i-1]+alpha_[i-1]/hh_;
00979     }
00980     psi_[currentOrder_] = temp1;
00981   }
00982   alpha_s_ = ST::zero();
00983   for (int i=0;i<currentOrder_;++i) {
00984     alpha_s_ = alpha_s_ - Scalar(ST::one()/(i+ST::one()));
00985   }
00986   RCP<Teuchos::FancyOStream> out = this->getOStream();
00987   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00988   Teuchos::OSTab ostab(out,1,"updateCoeffs_");
00989   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
00990     for (int i=0;i<=maxOrder_;++i) {
00991       *out << "alpha_[" << i << "] = " << alpha_[i] << std::endl;
00992       *out << "beta_[" << i << "] = " << beta_[i] << std::endl;
00993       *out << "gamma_[" << i << "] = " << gamma_[i] << std::endl;
00994       *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
00995       *out << "alpha_s_ = " << alpha_s_ << std::endl;
00996     }
00997   }
00998 }
00999 
01000 
01001 template<class Scalar>
01002 void ImplicitBDFStepper<Scalar>::initialize_()
01003 {
01004 
01005   using Teuchos::as;
01006   typedef Teuchos::ScalarTraits<Scalar> ST;
01007   using Thyra::createMember;
01008 
01009   RCP<Teuchos::FancyOStream> out = this->getOStream();
01010   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01011   const bool doTrace = (as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH));
01012   Teuchos::OSTab ostab(out,1,"initialize_");
01013 
01014   if (doTrace) {
01015     *out
01016       << "\nEntering " << this->Teuchos::Describable::description()
01017       << "::initialize_()...\n";
01018   }
01019 
01020   TEUCHOS_TEST_FOR_EXCEPT(model_ == Teuchos::null);
01021   TEUCHOS_TEST_FOR_EXCEPT(solver_ == Teuchos::null);
01022   TEUCHOS_ASSERT(haveInitialCondition_);
01023 
01024   // Initialize Parameter List if none provided.
01025   if (parameterList_ == Teuchos::null) {
01026     RCP<Teuchos::ParameterList> emptyParameterList = Teuchos::rcp(new Teuchos::ParameterList);
01027     this->setParameterList(emptyParameterList);
01028   }
01029 
01030   // Initialize StepControl
01031   if (stepControl_ == Teuchos::null) {
01032     RCP<ImplicitBDFStepperStepControl<Scalar> > implicitBDFStepperStepControl =
01033       Teuchos::rcp(new ImplicitBDFStepperStepControl<Scalar>());
01034     RCP<Teuchos::ParameterList> stepControlPL = 
01035       Teuchos::sublist(parameterList_, RythmosStepControlSettings_name);
01036     implicitBDFStepperStepControl->setParameterList(stepControlPL);
01037     this->setStepControlStrategy(implicitBDFStepperStepControl);
01038   }
01039   stepControl_->initialize(*this);
01040 
01041   maxOrder_ = stepControl_->getMaxOrder(); // maximum order
01042   TEUCHOS_TEST_FOR_EXCEPTION(
01043       !((1 <= maxOrder_) && (maxOrder_ <= 5)), std::logic_error,
01044       "Error, maxOrder returned from stepControl_->getMaxOrder() = " << maxOrder_ << " is outside range of [1,5]!\n"
01045       );
01046 
01047   Scalar zero = ST::zero();
01048 
01049   currentOrder_ = 1; // Current order of integration
01050   usedOrder_ = 1;  // order used in current step (used after currentOrder_ is updated)
01051   usedStep_ = zero;
01052   nscsco_  =  0;
01053   LETvalue_ = zero;
01054 
01055   alpha_.clear();  // $\alpha_j(n)=h_n/\psi_j(n)$ coefficient used in local error test
01056   // note:   $h_n$ = current step size, n = current time step
01057   gamma_.clear();  // calculate time derivative of history array for predictor 
01058   beta_.clear();   // coefficients used to evaluate predictor from history array
01059   psi_.clear();    // $\psi_j(n) = t_n-t_{n-j}$ intermediary variable used to 
01060   // compute $\beta_j(n):$
01061   for (int i=0 ; i<=maxOrder_ ; ++i) {
01062     alpha_.push_back(zero);
01063     beta_.push_back(zero);
01064     gamma_.push_back(zero);
01065     psi_.push_back(zero);
01066   }
01067   alpha_s_=Scalar(-ST::one());  // $\alpha_s$ fixed-leading coefficient of this BDF method
01068   hh_=zero;
01069   numberOfSteps_=0;   // number of total time integration steps taken
01070   nef_ = 0;
01071 
01072   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
01073     *out << "alpha_s_ = " << alpha_s_ << std::endl;
01074     for (int i=0 ; i<=maxOrder_ ; ++i) {
01075       *out << "alpha_[" << i << "] = " << alpha_[i] << std::endl;
01076       *out << "beta_[" << i << "] = " << beta_[i] << std::endl;
01077       *out << "gamma_[" << i << "] = " << gamma_[i] << std::endl;
01078       *out << "psi_[" << i << "] = " << psi_[i] << std::endl;
01079     }
01080     *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
01081   }
01082 
01083   // setInitialCondition initialized xHistory with xn0, xpn0.  Now we add the rest of the vectors.
01084   // Store maxOrder_+1 vectors
01085   for (int i=2 ; i<=maxOrder_ ; ++i) {
01086     xHistory_.push_back(createMember(xn0_->space())); 
01087     V_S(xHistory_[i].ptr(),zero);
01088   }
01089 
01090   isInitialized_ = true;
01091 
01092   if (doTrace) {
01093     *out
01094       << "\nLeaving " << this->Teuchos::Describable::description()
01095       << "::initialize_()...\n";
01096   }
01097 
01098 }
01099 
01100 
01101 template<class Scalar>
01102 void ImplicitBDFStepper<Scalar>::completeStep_()
01103 {
01104 
01105   using Teuchos::as;
01106   typedef Teuchos::ScalarTraits<Scalar> ST;
01107 
01108 #ifdef HAVE_RYTHMOS_DEBUG
01109   TEUCHOS_TEST_FOR_EXCEPT(ST::isnaninf(hh_));
01110 #endif  
01111 
01112   numberOfSteps_ ++;
01113   nef_ = 0;
01114   usedStep_ = hh_;
01115   usedOrder_ = currentOrder_;
01116   time_ += hh_;
01117   RCP<Teuchos::FancyOStream> out = this->getOStream();
01118   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
01119   Teuchos::OSTab ostab(out,1,"completeStep_");
01120 
01121   if ( as<int>(verbLevel) >= as<int>(Teuchos::VERB_HIGH) ) {
01122     *out << "numberOfSteps_ = " << numberOfSteps_ << std::endl;
01123     *out << "time_ = " << time_ << std::endl;
01124   }
01125   
01126   // 03/22/04 tscoffe:  Note that updating the history has nothing to do with
01127   // the step-size and everything to do with the newton correction vectors.
01128   updateHistory_();
01129 
01130 }
01131 
01132 template<class Scalar>
01133 void ImplicitBDFStepper<Scalar>::setStepControlData(const StepperBase<Scalar> & stepper)
01134 {
01135   if (!isInitialized_) {
01136     initialize_();
01137   }
01138   stepControl_->setStepControlData(stepper);
01139 }
01140 
01141 // 
01142 // Explicit Instantiation macro
01143 //
01144 // Must be expanded from within the Rythmos namespace!
01145 //
01146 
01147 #define RYTHMOS_IMPLICITBDF_STEPPER_INSTANT(SCALAR) \
01148   \
01149   template class ImplicitBDFStepper< SCALAR >; \
01150   \
01151   template RCP< ImplicitBDFStepper< SCALAR > > \
01152   implicitBDFStepper();  \
01153   \
01154   template RCP< ImplicitBDFStepper< SCALAR > > \
01155   implicitBDFStepper( \
01156     const RCP<Thyra::ModelEvaluator< SCALAR > >& model, \
01157     const RCP<Thyra::NonlinearSolverBase< SCALAR > >& solver, \
01158     const RCP<Teuchos::ParameterList>& parameterList \
01159     ); \
01160 
01161 
01162 } // namespace Rythmos
01163 
01164 
01165 #endif //Rythmos_IMPLICITBDF_STEPPER_DEF_H
 All Classes Functions Variables Typedefs Friends