C++ meb::OutArgs类代码示例

TEUCHOS_UNIT_TEST( Rythmos_ForwardSensitivityExplicitModelEvaluator, args ) {
  RCP<ForwardSensitivityExplicitModelEvaluator<double> > model =
    forwardSensitivityExplicitModelEvaluator<double>();
  RCP<SinCosModel> innerModel = sinCosModel();
  {
    RCP<ParameterList> pl = Teuchos::parameterList();
    pl->set("Accept model parameters",true);
    pl->set("Implicit model formulation",false);
    innerModel->setParameterList(pl);
  }
  model->initializeStructure(innerModel, 0 );
  typedef Thyra::ModelEvaluatorBase MEB;
  {
    MEB::InArgs<double> inArgs = model->createInArgs();
    TEST_EQUALITY_CONST( inArgs.supports(MEB::IN_ARG_t), true );
    TEST_EQUALITY_CONST( inArgs.supports(MEB::IN_ARG_x), true );
    TEST_EQUALITY_CONST( inArgs.supports(MEB::IN_ARG_x_dot), false );
    TEST_EQUALITY_CONST( inArgs.supports(MEB::IN_ARG_alpha), false );
    TEST_EQUALITY_CONST( inArgs.supports(MEB::IN_ARG_beta), true );
  }
  {
    MEB::OutArgs<double> outArgs = model->createOutArgs();
    TEST_EQUALITY_CONST( outArgs.supports(MEB::OUT_ARG_f), true );
    TEST_EQUALITY_CONST( outArgs.supports(MEB::OUT_ARG_W_op), false );
    TEST_EQUALITY_CONST( outArgs.supports(MEB::OUT_ARG_W), false );
  }
}

void eval_model_explicit(
    const Thyra::ModelEvaluator<Scalar> &model,
    Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
    const VectorBase<Scalar>& x_in,
    const typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag &t_in,
    const Ptr<VectorBase<Scalar> >& f_out
    )
{
  typedef Thyra::ModelEvaluatorBase MEB;
  MEB::InArgs<Scalar> inArgs = model.createInArgs();
  MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
  inArgs.setArgs(basePoint);
  inArgs.set_x(Teuchos::rcp(&x_in,false));
  if (inArgs.supports(MEB::IN_ARG_t)) {
    inArgs.set_t(t_in);
  }
  // For model evaluators whose state function f(x, x_dot, t) describes
  // an implicit ODE, and which accept an optional x_dot input argument,
  // make sure the latter is set to null in order to request the evaluation
  // of a state function corresponding to the explicit ODE formulation
  // x_dot = f(x, t)
  if (inArgs.supports(MEB::IN_ARG_x_dot)) {
    inArgs.set_x_dot(Teuchos::null);
  }
  outArgs.set_f(Teuchos::rcp(&*f_out,false));
  model.evalModel(inArgs,outArgs);
}

void assertValidModel(
  const StepperBase<Scalar>& stepper,
  const Thyra::ModelEvaluator<Scalar>& model
  )
{

  typedef Thyra::ModelEvaluatorBase MEB;

  TEUCHOS_ASSERT(stepper.acceptsModel());

  const MEB::InArgs<Scalar> inArgs = model.createInArgs();
  const MEB::OutArgs<Scalar> outArgs = model.createOutArgs();

  //TEUCHOS_ASSERT(inArgs.supports(MEB::IN_ARG_t));
  TEUCHOS_ASSERT(inArgs.supports(MEB::IN_ARG_x));
  TEUCHOS_ASSERT(outArgs.supports(MEB::OUT_ARG_f));
  
  if (stepper.isImplicit()) { // implicit stepper
    TEUCHOS_ASSERT( inArgs.supports(MEB::IN_ARG_x_dot) );
    TEUCHOS_ASSERT( inArgs.supports(MEB::IN_ARG_alpha) );
    TEUCHOS_ASSERT( inArgs.supports(MEB::IN_ARG_beta) );
    TEUCHOS_ASSERT( outArgs.supports(MEB::OUT_ARG_W) );
  } 
  //else { // explicit stepper
  //  TEUCHOS_ASSERT( !inArgs.supports(MEB::IN_ARG_x_dot) );
  //  TEUCHOS_ASSERT( !inArgs.supports(MEB::IN_ARG_alpha) );
  //  TEUCHOS_ASSERT( !inArgs.supports(MEB::IN_ARG_beta) );
  //  TEUCHOS_ASSERT( !outArgs.supports(MEB::OUT_ARG_W) );
  //}

}

int TriKota::ThyraDirectApplicInterface::derived_map_ac(const Dakota::String& ac_name)
{

  if (App != Teuchos::null) {

    // Test for consistency of problem definition between ModelEval and Dakota
    TEST_FOR_EXCEPTION(numVars > numParameters, std::logic_error,
                       "TriKota_Dakota Adapter Error: ");
    TEST_FOR_EXCEPTION(numFns > numResponses, std::logic_error,
                       "TriKota_Dakota Adapter Error: ");
    TEST_FOR_EXCEPTION(hessFlag, std::logic_error,
                       "TriKota_Dakota Adapter Error: ");

    MEB::InArgs<double> inArgs = App->createInArgs();
    MEB::OutArgs<double> outArgs = App->createOutArgs();

    TEST_FOR_EXCEPTION(gradFlag && !supportsSensitivities, std::logic_error,
                       "TriKota_Dakota Adapter Error: ");

    // Load parameters from Dakota to ModelEval data structure
    {
      Thyra::DetachedVectorView<double> my_p(model_p);
      for (unsigned int i=0; i<numVars; i++) my_p[i]=xC[i];
    }

    // Evaluate model
    inArgs.set_p(0,model_p);
    outArgs.set_g(0,model_g);
    if (gradFlag) outArgs.set_DgDp(0,0,
      MEB::DerivativeMultiVector<double>(model_dgdp,orientation));
    App->evalModel(inArgs, outArgs);

    Thyra::DetachedVectorView<double> my_g(model_g);
    for (unsigned int j=0; j<numFns; j++) fnVals[j]= my_g[j];

    if (gradFlag) {
      if (orientation == MEB::DERIV_MV_BY_COL) {
        for (unsigned int j=0; j<numVars; j++) {
          Thyra::DetachedVectorView<double>
             my_dgdp_j(model_dgdp->col(j));
          for (unsigned int i=0; i<numFns; i++)  fnGrads[i][j]= my_dgdp_j[i];
        }
      }
      else {
        for (unsigned int j=0; j<numFns; j++) {
          Thyra::DetachedVectorView<double>
             my_dgdp_j(model_dgdp->col(j));
          for (unsigned int i=0; i<numVars; i++) fnGrads[j][i]= my_dgdp_j[i]; 
        }
      }
    }
  }
  else {
    TEST_FOR_EXCEPTION(parallelLib.parallel_configuration().ea_parallel_level().server_intra_communicator()
               != MPI_COMM_NULL, std::logic_error,
              "\nTriKota Parallelism Error: ModelEvaluator=null, but analysis_comm != MPI_COMMM_NULL");
  }

  return 0;
}

void ExplicitModelEvaluator<Scalar>::
buildInverseMassMatrix() const
{
  typedef Thyra::ModelEvaluatorBase MEB;
  using Teuchos::RCP;
  using Thyra::createMember;
  
  RCP<const Thyra::ModelEvaluator<Scalar> > me = this->getUnderlyingModel();

  // first allocate space for the mass matrix
  RCP<Thyra::LinearOpBase<Scalar> > mass = me->create_W_op();

  // intialize a zero to get rid of the x-dot 
  if(zero_==Teuchos::null) {
    zero_ = Thyra::createMember(*me->get_x_space());
    Thyra::assign(zero_.ptr(),0.0);
  }
  
  // request only the mass matrix from the physics
  // Model evaluator builds: alpha*u_dot + beta*F(u) = 0
  MEB::InArgs<Scalar>  inArgs  = me->createInArgs();
  inArgs.set_x(createMember(me->get_x_space()));
  inArgs.set_x_dot(zero_);
  inArgs.set_alpha(-1.0);
  inArgs.set_beta(0.0);

  // set the one time beta to ensure dirichlet conditions
  // are correctly included in the mass matrix: do it for
  // both epetra and Tpetra. If a panzer model evaluator has
  // not been passed in...oh well you get what you asked for!
  if(panzerModel_!=Teuchos::null)
    panzerModel_->setOneTimeDirichletBeta(-1.0);
  else if(panzerEpetraModel_!=Teuchos::null)
    panzerEpetraModel_->setOneTimeDirichletBeta(-1.0);

  // set only the mass matrix
  MEB::OutArgs<Scalar> outArgs = me->createOutArgs();
  outArgs.set_W_op(mass);

  // this will fill the mass matrix operator 
  me->evalModel(inArgs,outArgs);

  if(!massLumping_) {
    invMassMatrix_ = Thyra::inverse<Scalar>(*me->get_W_factory(),mass);
  }
  else {
    // build lumped mass matrix (assumes all positive mass entries, does a simple sum)
    Teuchos::RCP<Thyra::VectorBase<Scalar> > ones = Thyra::createMember(*mass->domain());
    Thyra::assign(ones.ptr(),1.0);

    RCP<Thyra::VectorBase<Scalar> > invLumpMass = Thyra::createMember(*mass->range());
    Thyra::apply(*mass,Thyra::NOTRANS,*ones,invLumpMass.ptr());
    Thyra::reciprocal(*invLumpMass,invLumpMass.ptr());

    invMassMatrix_ = Thyra::diagonal(invLumpMass);
  }
}

Thyra::ModelEvaluatorBase::OutArgs<Scalar>
TimeDiscretizedBackwardEulerModelEvaluator<Scalar>::createOutArgsImpl() const
{
  typedef Thyra::ModelEvaluatorBase MEB;
  MEB::OutArgs<Scalar> daeOutArgs = daeModel_->createOutArgs();
  MEB::OutArgsSetup<Scalar> outArgs;
  outArgs.setModelEvalDescription(this->description());
  outArgs.setSupports(MEB::OUT_ARG_f);
  outArgs.setSupports(MEB::OUT_ARG_W_op);
  outArgs.set_W_properties(daeOutArgs.get_W_properties());
  return outArgs;
}

ModelEvaluatorBase::OutArgs<Scalar>
DefaultModelEvaluatorWithSolveFactory<Scalar>::createOutArgsImpl() const
{
  typedef ModelEvaluatorBase MEB;
  const RCP<const ModelEvaluator<Scalar> >
    thyraModel = this->getUnderlyingModel();
  const MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
  MEB::OutArgsSetup<Scalar> outArgs;
  outArgs.setModelEvalDescription(this->description());
  outArgs.set_Np_Ng(wrappedOutArgs.Np(),wrappedOutArgs.Ng());
  outArgs.setSupports(wrappedOutArgs);
  outArgs.setSupports(MEB::OUT_ARG_W,
    wrappedOutArgs.supports(MEB::OUT_ARG_W_op)&&W_factory_.get()!=NULL);
  return outArgs;
}

RCP<Thyra::VectorBase<Scalar> > eval_f_t(
    const Thyra::ModelEvaluator<Scalar>& me,
    Scalar t
    ) {
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef Thyra::ModelEvaluatorBase MEB;
  MEB::InArgs<Scalar> inArgs = me.createInArgs();
  inArgs.set_t(t);
  MEB::OutArgs<Scalar> outArgs = me.createOutArgs();
  RCP<Thyra::VectorBase<Scalar> > f_out = Thyra::createMember(me.get_f_space());
  V_S(outArg(*f_out),ST::zero());
  outArgs.set_f(f_out);
  me.evalModel(inArgs,outArgs);
  return f_out;
}

ModelEvaluatorBase::OutArgs<Scalar>
DefaultStateEliminationModelEvaluator<Scalar>::createOutArgsImpl() const
{
  typedef ModelEvaluatorBase MEB;
  const Teuchos::RCP<const ModelEvaluator<Scalar> >
    thyraModel = this->getUnderlyingModel();
  const MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
  const int Np = wrappedOutArgs.Np(), Ng = wrappedOutArgs.Ng();
  MEB::OutArgsSetup<Scalar> outArgs;
  outArgs.setModelEvalDescription(this->description());
  outArgs.set_Np_Ng(Np,Ng);
  outArgs.setSupports(wrappedOutArgs);
  outArgs.setUnsupportsAndRelated(MEB::IN_ARG_x); // wipe out DgDx ...
  outArgs.setUnsupportsAndRelated(MEB::OUT_ARG_f); // wipe out f, W, DfDp ...
  return outArgs;
}

RCP<LinearOpBase<Scalar> >
ModelEvaluatorDefaultBase<Scalar>::create_DgDp_op_impl(int j, int l) const
{
  typedef ModelEvaluatorBase MEB;
  MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
  TEUCHOS_TEST_FOR_EXCEPTION(
    outArgs.supports(MEB::OUT_ARG_DgDp,j,l).supports(MEB::DERIV_LINEAR_OP),
    std::logic_error,
    "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
    " supports the LinearOpBase form of DgDp("<<j<<","<<l<<")"
    " (as determined from its OutArgs object created by createOutArgsImpl())"
    " but this function create_DgDp_op_impl(...) has not been overriden"
    " to create such an object!"
    );
  return Teuchos::null;
}

void eval_model_explicit(
    const Thyra::ModelEvaluator<Scalar> &model,
    Thyra::ModelEvaluatorBase::InArgs<Scalar> &basePoint,
    const VectorBase<Scalar>& x_in,
    const typename Thyra::ModelEvaluatorBase::InArgs<Scalar>::ScalarMag &t_in,
    const Ptr<VectorBase<Scalar> >& f_out
    )
{
  typedef Thyra::ModelEvaluatorBase MEB;
  MEB::InArgs<Scalar> inArgs = model.createInArgs();
  MEB::OutArgs<Scalar> outArgs = model.createOutArgs();
  inArgs.setArgs(basePoint);
  inArgs.set_x(Teuchos::rcp(&x_in,false));
  if (inArgs.supports(MEB::IN_ARG_t)) {
    inArgs.set_t(t_in);
  }
  outArgs.set_f(Teuchos::rcp(&*f_out,false));
  model.evalModel(inArgs,outArgs);
}

void ForwardSensitivityExplicitModelEvaluator<Scalar>::computeDerivativeMatrices(
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
  ) const
{
  TEUCHOS_ASSERT( !is_null(stateModel_) );

  typedef Thyra::ModelEvaluatorBase MEB;
  typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;

  Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();

  MEB::InArgs<Scalar> inArgs = stateBasePoint_;
  MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
  
  if (is_null(DfDx_)) {
    DfDx_ = stateModel_->create_W_op();
  }
  if (inArgs.supports(MEB::IN_ARG_beta)) {
    inArgs.set_beta(1.0);
  }
  outArgs.set_W_op(DfDx_);

  if (is_null(DfDp_)) {
    DfDp_ = Thyra::create_DfDp_mv(
      *stateModel_,p_index_,
      MEB::DERIV_MV_BY_COL
      ).getMultiVector();
  }
  outArgs.set_DfDp(
    p_index_,
    MEB::Derivative<Scalar>(DfDp_,MEB::DERIV_MV_BY_COL)
    );
  
  VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
  stateModel_->evalModel(inArgs,outArgs);
  

}

void ModelEvaluatorDefaultBase<Scalar>::initializeDefaultBase()
{

  typedef ModelEvaluatorBase MEB;

  // In case we throw half way thorugh, set to uninitialized
  isInitialized_ = false;
  default_W_support_ = false;

  //
  // A) Get the InArgs and OutArgs from the subclass
  //
  
  const MEB::InArgs<Scalar> inArgs = this->createInArgs();
  const MEB::OutArgs<Scalar> outArgsImpl = this->createOutArgsImpl();

  //
  // B) Validate the subclasses InArgs and OutArgs
  //

#ifdef TEUCHOS_DEBUG
  assertInArgsOutArgsSetup( this->description(), inArgs, outArgsImpl );
#endif // TEUCHOS_DEBUG

  //
  // C) Set up support for default derivative objects and prototype OutArgs
  //

  const int l_Ng = outArgsImpl.Ng();
  const int l_Np = outArgsImpl.Np();

  // Set support for all outputs supported in the underly implementation
  MEB::OutArgsSetup<Scalar> outArgs;
  outArgs.setModelEvalDescription(this->description());
  outArgs.set_Np_Ng(l_Np,l_Ng);
  outArgs.setSupports(outArgsImpl);

  // DfDp
  DfDp_default_op_support_.clear();
  if (outArgs.supports(MEB::OUT_ARG_f)) {
    for ( int l = 0; l < l_Np; ++l ) {
      const MEB::DerivativeSupport DfDp_l_impl_support =
        outArgsImpl.supports(MEB::OUT_ARG_DfDp,l);
      const DefaultDerivLinearOpSupport DfDp_l_op_support =
        determineDefaultDerivLinearOpSupport(DfDp_l_impl_support);
      DfDp_default_op_support_.push_back(DfDp_l_op_support);
      outArgs.setSupports(
        MEB::OUT_ARG_DfDp, l,
        updateDefaultLinearOpSupport(
          DfDp_l_impl_support, DfDp_l_op_support
          )
        );
    }
  }

  // DgDx_dot
  DgDx_dot_default_op_support_.clear();
  for ( int j = 0; j < l_Ng; ++j ) {
    const MEB::DerivativeSupport DgDx_dot_j_impl_support =
      outArgsImpl.supports(MEB::OUT_ARG_DgDx_dot,j);
    const DefaultDerivLinearOpSupport DgDx_dot_j_op_support =
      determineDefaultDerivLinearOpSupport(DgDx_dot_j_impl_support);
    DgDx_dot_default_op_support_.push_back(DgDx_dot_j_op_support);
    outArgs.setSupports(
      MEB::OUT_ARG_DgDx_dot, j,
      updateDefaultLinearOpSupport(
        DgDx_dot_j_impl_support, DgDx_dot_j_op_support
        )
      );
  }

  // DgDx
  DgDx_default_op_support_.clear();
  for ( int j = 0; j < l_Ng; ++j ) {
    const MEB::DerivativeSupport DgDx_j_impl_support =
      outArgsImpl.supports(MEB::OUT_ARG_DgDx,j);
    const DefaultDerivLinearOpSupport DgDx_j_op_support =
      determineDefaultDerivLinearOpSupport(DgDx_j_impl_support);
    DgDx_default_op_support_.push_back(DgDx_j_op_support);
    outArgs.setSupports(
      MEB::OUT_ARG_DgDx, j,
      updateDefaultLinearOpSupport(
        DgDx_j_impl_support, DgDx_j_op_support
        )
      );
  }

  // DgDp
  DgDp_default_op_support_.clear();
  DgDp_default_mv_support_.clear();
  for ( int j = 0; j < l_Ng; ++j ) {
    DgDp_default_op_support_.push_back(Array<DefaultDerivLinearOpSupport>());
    DgDp_default_mv_support_.push_back(Array<DefaultDerivMvAdjointSupport>());
    for ( int l = 0; l < l_Np; ++l ) {
      const MEB::DerivativeSupport DgDp_j_l_impl_support =
        outArgsImpl.supports(MEB::OUT_ARG_DgDp,j,l);
      // LinearOpBase support
      const DefaultDerivLinearOpSupport DgDp_j_l_op_support =
        determineDefaultDerivLinearOpSupport(DgDp_j_l_impl_support);
      DgDp_default_op_support_[j].push_back(DgDp_j_l_op_support);
//.........这里部分代码省略.........

void DiagonalImplicitRKModelEvaluator<Scalar>::evalModelImpl(
  const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs_stage,
  const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs_stage
  ) const
{

  typedef ScalarTraits<Scalar> ST;
  typedef Thyra::ModelEvaluatorBase MEB;

  TEUCHOS_TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
      "Error!  initializeDIRKModel must be called before evalModel\n"
      );

  TEUCHOS_TEST_FOR_EXCEPTION( !setTimeStepPointCalled_, std::logic_error,
      "Error!  setTimeStepPoint must be called before evalModel"
      );

  TEUCHOS_TEST_FOR_EXCEPTION( currentStage_ == -1, std::logic_error,
      "Error!  setCurrentStage must be called before evalModel"
      );

  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
    "Rythmos::DiagonalImplicitRKModelEvaluator",inArgs_stage,outArgs_stage,daeModel_
    );

  //
  // A) Unwrap the inArgs and outArgs 
  //

  const RCP<const Thyra::VectorBase<Scalar> > x_in = inArgs_stage.get_x();
  const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs_stage.get_f();
  const RCP<Thyra::LinearOpBase<Scalar> > W_op_out = outArgs_stage.get_W_op();

  //
  // B) Assemble f_out and W_op_out for given stage
  //

  MEB::InArgs<Scalar> daeInArgs = daeModel_->createInArgs();
  MEB::OutArgs<Scalar> daeOutArgs = daeModel_->createOutArgs();
  const RCP<Thyra::VectorBase<Scalar> > x_i = createMember(daeModel_->get_x_space());
  daeInArgs.setArgs(basePoint_);
  
  // B.1) Setup the DAE's inArgs for stage f(currentStage_) ...
  V_V(stage_derivatives_->getNonconstVectorBlock(currentStage_).ptr(),*x_in);
  assembleIRKState( currentStage_, dirkButcherTableau_->A(), delta_t_, *x_old_, *stage_derivatives_, outArg(*x_i) );
  daeInArgs.set_x( x_i );
  daeInArgs.set_x_dot( x_in );
  daeInArgs.set_t( t_old_ + dirkButcherTableau_->c()(currentStage_) * delta_t_ );
  daeInArgs.set_alpha(ST::one());
  daeInArgs.set_beta( delta_t_ * dirkButcherTableau_->A()(currentStage_,currentStage_) );

  // B.2) Setup the DAE's outArgs for stage f(i) ...
  if (!is_null(f_out))
    daeOutArgs.set_f( f_out );
  if (!is_null(W_op_out))
    daeOutArgs.set_W_op(W_op_out);

  // B.3) Compute f_out(i) and/or W_op_out ...
  daeModel_->evalModel( daeInArgs, daeOutArgs );
  daeOutArgs.set_f(Teuchos::null);
  daeOutArgs.set_W_op(Teuchos::null);
  
  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
  
}

void ImplicitRKModelEvaluator<Scalar>::evalModelImpl(
  const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs_bar,
  const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs_bar
  ) const
{

  using Teuchos::rcp_dynamic_cast;
  typedef ScalarTraits<Scalar> ST;
  typedef Thyra::ModelEvaluatorBase MEB;
  typedef Thyra::VectorBase<Scalar> VB;
  typedef Thyra::ProductVectorBase<Scalar> PVB;
  typedef Thyra::BlockedLinearOpBase<Scalar> BLWB;

  TEST_FOR_EXCEPTION( !isInitialized_, std::logic_error,
      "Error!  initializeIRKModel must be called before evalModel\n"
      );

  TEST_FOR_EXCEPTION( !setTimeStepPointCalled_, std::logic_error,
      "Error!  setTimeStepPoint must be called before evalModel"
      );

  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
    "Rythmos::ImplicitRKModelEvaluator",inArgs_bar,outArgs_bar,daeModel_
    );

  //
  // A) Unwrap the inArgs and outArgs to get at product vectors and block op
  //

  const RCP<const PVB> x_bar = rcp_dynamic_cast<const PVB>(inArgs_bar.get_x(), true);
  const RCP<PVB> f_bar = rcp_dynamic_cast<PVB>(outArgs_bar.get_f(), true);
  const RCP<BLWB> W_op_bar = rcp_dynamic_cast<BLWB>(outArgs_bar.get_W_op(), true);

  //
  // B) Assemble f_bar and W_op_bar by looping over stages
  //

  MEB::InArgs<Scalar> daeInArgs = daeModel_->createInArgs();
  MEB::OutArgs<Scalar> daeOutArgs = daeModel_->createOutArgs();
  const RCP<VB> x_i = createMember(daeModel_->get_x_space());
  daeInArgs.setArgs(basePoint_);
  
  const int numStages = irkButcherTableau_->numStages();

  for ( int i = 0; i < numStages; ++i ) {

    // B.1) Setup the DAE's inArgs for stage f(i) ...
    assembleIRKState( i, irkButcherTableau_->A(), delta_t_, *x_old_, *x_bar, outArg(*x_i) );
    daeInArgs.set_x( x_i );
    daeInArgs.set_x_dot( x_bar->getVectorBlock(i) );
    daeInArgs.set_t( t_old_ + irkButcherTableau_->c()(i) * delta_t_ );
    Scalar alpha = ST::zero();
    if (i == 0) {
      alpha = ST::one();
    } else {
      alpha = ST::zero();
    }
    Scalar beta = delta_t_ * irkButcherTableau_->A()(i,0);
    daeInArgs.set_alpha( alpha );
    daeInArgs.set_beta( beta );

    // B.2) Setup the DAE's outArgs for stage f(i) ...
    if (!is_null(f_bar))
      daeOutArgs.set_f( f_bar->getNonconstVectorBlock(i) );
    if (!is_null(W_op_bar)) {
      daeOutArgs.set_W_op(W_op_bar->getNonconstBlock(i,0));
    }

    // B.3) Compute f_bar(i) and/or W_op_bar(i,0) ...
    daeModel_->evalModel( daeInArgs, daeOutArgs );
    daeOutArgs.set_f(Teuchos::null);
    daeOutArgs.set_W_op(Teuchos::null);
    
    // B.4) Evaluate the rest of the W_op_bar(i,j=1...numStages-1) ...
    if (!is_null(W_op_bar)) {
      for ( int j = 1; j < numStages; ++j ) {
        alpha = ST::zero();
        if (i == j) {
          alpha = ST::one();
        } else {
          alpha = ST::zero();
        }
        beta = delta_t_ * irkButcherTableau_->A()(i,j);
        daeInArgs.set_alpha( alpha );
        daeInArgs.set_beta( beta );
        daeOutArgs.set_W_op(W_op_bar->getNonconstBlock(i,j));
        daeModel_->evalModel( daeInArgs, daeOutArgs );
        daeOutArgs.set_W_op(Teuchos::null);
      }
    }

  }
  
  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
  
}