Rythmos_ForwardSensitivityImplicitModelEvaluator.hpp

//@HEADER
// ***********************************************************************
//
//                           Rythmos Package
//                 Copyright (2006) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Todd S. Coffey (tscoffe@sandia.gov)
//
// ***********************************************************************
//@HEADER

#ifndef RYTHMOS_FORWARD_SENSITIVITY_IMPLICIT_MODEL_EVALUATOR_HPP
#define RYTHMOS_FORWARD_SENSITIVITY_IMPLICIT_MODEL_EVALUATOR_HPP


#include "Rythmos_IntegratorBase.hpp"
#include "Rythmos_ForwardSensitivityModelEvaluatorBase.hpp"
#include "Thyra_ModelEvaluator.hpp" // Interface
#include "Thyra_StateFuncModelEvaluatorBase.hpp" // Implementation
#include "Thyra_DefaultProductVectorSpace.hpp"
#include "Thyra_PhysicallyBlockedLinearOpWithSolveBase.hpp" // Interface
#include "Thyra_DefaultBlockedTriangularLinearOpWithSolve.hpp" // Implementation
#include "Thyra_ModelEvaluatorDelegatorBase.hpp"
#include "Thyra_ModelEvaluatorHelpers.hpp"
#include "Thyra_DefaultMultiVectorProductVectorSpace.hpp"
#include "Thyra_DefaultMultiVectorProductVector.hpp"
#include "Thyra_DefaultMultiVectorLinearOpWithSolve.hpp"
#include "Teuchos_implicit_cast.hpp"
#include "Teuchos_Assert.hpp"


namespace Rythmos {


template<class Scalar>
class ForwardSensitivityImplicitModelEvaluator
  : virtual public ForwardSensitivityModelEvaluatorBase<Scalar>
{
public:


  ForwardSensitivityImplicitModelEvaluator();

  void initializeStructure(
    const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
    const int p_index
    );
  
  RCP<const Thyra::ModelEvaluator<Scalar> >
  getStateModel() const;
  
  int get_p_index() const;

  void setStateIntegrator(
    const RCP<IntegratorBase<Scalar> > &stateIntegrator,
    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint
    );

  void initializePointState(
    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
    const RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
    const Scalar &coeff_x_dot,
    const Scalar &coeff_x,
    const RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot = Teuchos::null,
    const RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp = Teuchos::null
    );

  // 2007/05/22: rabartl: ToDo: Add an InterpolationBufferBase
  // stateInterpBuffer object to the initailizeState(...) function that can be
  // used to get x and x_dot at different points in time t.  Then, modify the
  // logic to recompute all of the needed matrices if t != t_base (as passed
  // in through stateBasePoint).  The values of x(t) and xdot(t) can then be
  // gotten from the stateInterpBuffer object!
  


  RCP<const Thyra::VectorSpaceBase<Scalar> > get_x_space() const;
  RCP<const Thyra::VectorSpaceBase<Scalar> > get_f_space() const;
  Thyra::ModelEvaluatorBase::InArgs<Scalar> getNominalValues() const;
  RCP<Thyra::LinearOpWithSolveBase<Scalar> > create_W() const;
  Thyra::ModelEvaluatorBase::InArgs<Scalar> createInArgs() const;



  void initializeState(
    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
    const RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
    const Scalar &coeff_x_dot,
    const Scalar &coeff_x,
    const RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot = Teuchos::null,
    const RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp = Teuchos::null
    )
    {
      initializePointState( stateBasePoint, W_tilde, coeff_x_dot, coeff_x,
        DfDx_dot, DfDp );
    }


private:


  Thyra::ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const;
  void evalModelImpl(
    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
    const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
    ) const;


private:

  // /////////////////////////
  // Private data members

  RCP<const Thyra::ModelEvaluator<Scalar> > stateModel_;
  int p_index_;
  int np_;

  RCP<IntegratorBase<Scalar> > stateIntegrator_;

  RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > s_bar_space_;
  RCP<const Thyra::DefaultMultiVectorProductVectorSpace<Scalar> > f_sens_space_;

  Thyra::ModelEvaluatorBase::InArgs<Scalar> nominalValues_;

  mutable Thyra::ModelEvaluatorBase::InArgs<Scalar> stateBasePoint_;

  mutable RCP<const Thyra::LinearOpWithSolveBase<Scalar> > W_tilde_;
  mutable Scalar coeff_x_dot_;
  mutable Scalar coeff_x_;
  mutable RCP<const Thyra::LinearOpBase<Scalar> > DfDx_dot_;
  mutable RCP<const Thyra::MultiVectorBase<Scalar> > DfDp_;

  mutable RCP<Thyra::LinearOpWithSolveBase<Scalar> > W_tilde_compute_;
  mutable RCP<Thyra::LinearOpBase<Scalar> > DfDx_dot_compute_;
  mutable RCP<Thyra::MultiVectorBase<Scalar> > DfDp_compute_;

  // /////////////////////////
  // Private member functions

  void wrapNominalValuesAndBounds();

  void computeDerivativeMatrices(
    const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
    ) const;
  
};


// /////////////////////////////////
// Implementations


// Constructors/Intializers/Accessors

template<class Scalar>
RCP<ForwardSensitivityImplicitModelEvaluator<Scalar> > forwardSensitivityImplicitModelEvaluator()
{
  return rcp(new ForwardSensitivityImplicitModelEvaluator<Scalar>());
}

template<class Scalar>
ForwardSensitivityImplicitModelEvaluator<Scalar>::ForwardSensitivityImplicitModelEvaluator()
  : p_index_(0), np_(-1)
{}


template<class Scalar>
void ForwardSensitivityImplicitModelEvaluator<Scalar>::initializeStructure(
  const RCP<const Thyra::ModelEvaluator<Scalar> > &stateModel,
  const int p_index
  )
{

  typedef Thyra::ModelEvaluatorBase MEB;

  //
  // Validate input
  //

  TEST_FOR_EXCEPT( is_null(stateModel) );
  TEST_FOR_EXCEPTION(
    !( 0 <= p_index && p_index < stateModel->Np() ), std::logic_error,
    "Error, p_index does not fall in the range [0,"<<(stateModel->Np()-1)<<"]!" );
  // ToDo: Validate support for DfDp!

  //
  // Set the input objects
  //

  stateModel_ = stateModel;
  p_index_ = p_index;
  np_ = stateModel_->get_p_space(p_index)->dim();

  //
  // Create the structure of the model
  //

  s_bar_space_ = Thyra::multiVectorProductVectorSpace(
    stateModel_->get_x_space(), np_
    );

  f_sens_space_ = Thyra::multiVectorProductVectorSpace(
    stateModel_->get_f_space(), np_
    );

  nominalValues_ = this->createInArgs();

  //
  // Wipe out matrix storage
  //

  stateBasePoint_ = MEB::InArgs<Scalar>();
  W_tilde_ = Teuchos::null;
  coeff_x_dot_ = 0.0;
  coeff_x_ = 0.0;
  DfDx_dot_ = Teuchos::null;
  DfDp_ = Teuchos::null;
  W_tilde_compute_ = Teuchos::null;
  DfDx_dot_compute_ = Teuchos::null;
  DfDp_compute_ = Teuchos::null;

}


template<class Scalar>
RCP<const Thyra::ModelEvaluator<Scalar> >
ForwardSensitivityImplicitModelEvaluator<Scalar>::getStateModel() const
{
  return stateModel_;
}


template<class Scalar>
int ForwardSensitivityImplicitModelEvaluator<Scalar>::get_p_index() const
{
  return p_index_;
}


template<class Scalar>
void ForwardSensitivityImplicitModelEvaluator<Scalar>::setStateIntegrator(
  const RCP<IntegratorBase<Scalar> > &stateIntegrator,
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint
  )
{
  stateIntegrator_ = stateIntegrator.assert_not_null();
  stateBasePoint_ = stateBasePoint;
}


template<class Scalar>
void ForwardSensitivityImplicitModelEvaluator<Scalar>::initializePointState(
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &stateBasePoint,
  const RCP<const Thyra::LinearOpWithSolveBase<Scalar> > &W_tilde,
  const Scalar &coeff_x_dot,
  const Scalar &coeff_x,
  const RCP<const Thyra::LinearOpBase<Scalar> > &DfDx_dot,
  const RCP<const Thyra::MultiVectorBase<Scalar> > &DfDp
  )
{

  typedef Thyra::ModelEvaluatorBase MEB;

#ifdef RYTHMOS_DEBUG
  TEST_FOR_EXCEPTION(
    is_null(stateModel_), std::logic_error,
    "Error, you must call intializeStructure(...) before you call initializeState(...)"
    );
  TEST_FOR_EXCEPT( is_null(stateBasePoint.get_x()) );
  TEST_FOR_EXCEPT( is_null(stateBasePoint.get_x_dot()) );
  TEST_FOR_EXCEPT( is_null(stateBasePoint.get_p(p_index_)) );
  // What about the other parameter values?  We really can't say anything
  // about these and we can't check them.  They can be null just fine.
  if (!is_null(W_tilde)) {
    THYRA_ASSERT_VEC_SPACES("",*W_tilde->domain(),*stateModel_->get_x_space());
    THYRA_ASSERT_VEC_SPACES("",*W_tilde->range(),*stateModel_->get_f_space());
  }
  if (!is_null(DfDx_dot)) {
    THYRA_ASSERT_VEC_SPACES("",*DfDx_dot->domain(),*stateModel_->get_x_space());
    THYRA_ASSERT_VEC_SPACES("",*DfDx_dot->range(),*stateModel_->get_f_space());
  }
  if (!is_null(DfDp)) {
    THYRA_ASSERT_VEC_SPACES("",*DfDp->domain(),*stateModel_->get_p_space(p_index_));
    THYRA_ASSERT_VEC_SPACES("",*DfDp->range(),*stateModel_->get_f_space());
  }
#endif

  stateBasePoint_ = stateBasePoint;

  // Set whatever derivatives where passed in.  If an input in null, then the
  // member will be null and the null linear operators will be computed later
  // just in time.

  W_tilde_ = W_tilde;
  coeff_x_dot_ = coeff_x_dot;
  coeff_x_ = coeff_x;
  DfDx_dot_ = DfDx_dot;
  DfDp_ = DfDp;

  wrapNominalValuesAndBounds();

}


// Public functions overridden from ModelEvaulator


template<class Scalar>
RCP<const Thyra::VectorSpaceBase<Scalar> >
ForwardSensitivityImplicitModelEvaluator<Scalar>::get_x_space() const
{
  return s_bar_space_;
}


template<class Scalar>
RCP<const Thyra::VectorSpaceBase<Scalar> >
ForwardSensitivityImplicitModelEvaluator<Scalar>::get_f_space() const
{
  return f_sens_space_;
}


template<class Scalar>
Thyra::ModelEvaluatorBase::InArgs<Scalar>
ForwardSensitivityImplicitModelEvaluator<Scalar>::getNominalValues() const
{
  return nominalValues_;
}


template<class Scalar>
RCP<Thyra::LinearOpWithSolveBase<Scalar> >
ForwardSensitivityImplicitModelEvaluator<Scalar>::create_W() const
{
  return Thyra::multiVectorLinearOpWithSolve<Scalar>();
}


template<class Scalar>
Thyra::ModelEvaluatorBase::InArgs<Scalar>
ForwardSensitivityImplicitModelEvaluator<Scalar>::createInArgs() const
{
  typedef Thyra::ModelEvaluatorBase MEB;
  MEB::InArgs<Scalar> stateModelInArgs = stateModel_->createInArgs();
  MEB::InArgsSetup<Scalar> inArgs;
  inArgs.setModelEvalDescription(this->description());
  inArgs.setSupports( MEB::IN_ARG_x_dot,
    stateModelInArgs.supports(MEB::IN_ARG_x_dot) );
  inArgs.setSupports( MEB::IN_ARG_x );
  inArgs.setSupports( MEB::IN_ARG_t );
  inArgs.setSupports( MEB::IN_ARG_alpha,
    stateModelInArgs.supports(MEB::IN_ARG_alpha) );
  inArgs.setSupports( MEB::IN_ARG_beta,
    stateModelInArgs.supports(MEB::IN_ARG_beta) );
  return inArgs;
}


// Private functions overridden from ModelEvaulatorDefaultBase


template<class Scalar>
Thyra::ModelEvaluatorBase::OutArgs<Scalar>
ForwardSensitivityImplicitModelEvaluator<Scalar>::createOutArgsImpl() const
{

  typedef Thyra::ModelEvaluatorBase MEB;

  MEB::OutArgs<Scalar> stateModelOutArgs = stateModel_->createOutArgs();
  MEB::OutArgsSetup<Scalar> outArgs;

  outArgs.setModelEvalDescription(this->description());

  outArgs.setSupports(MEB::OUT_ARG_f);

  if (stateModelOutArgs.supports(MEB::OUT_ARG_W) ) {
    outArgs.setSupports(MEB::OUT_ARG_W);
    outArgs.set_W_properties(stateModelOutArgs.get_W_properties());
  }

  return outArgs;

}


template<class Scalar>
void ForwardSensitivityImplicitModelEvaluator<Scalar>::evalModelImpl(
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
  ) const
{

  using Teuchos::rcp_dynamic_cast;
  typedef Teuchos::ScalarTraits<Scalar> ST;
  typedef Thyra::ModelEvaluatorBase MEB;
  typedef Teuchos::VerboseObjectTempState<Thyra::ModelEvaluatorBase> VOTSME;

  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
    "ForwardSensitivityImplicitModelEvaluator", inArgs, outArgs, Teuchos::null );

  //
  // Update the derivative matrices if they are not already updated for the
  // given time!.
  //
  
  {
#ifdef ENABLE_RYTHMOS_TIMERS
    TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityImplicitModelEvaluator::evalModel: computeMatrices");
#endif
    computeDerivativeMatrices(inArgs);
  }

  //
  // InArgs
  //

  RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
    s_bar = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
      inArgs.get_x().assert_not_null(), true
      );
  RCP<const Thyra::DefaultMultiVectorProductVector<Scalar> >
    s_bar_dot = rcp_dynamic_cast<const Thyra::DefaultMultiVectorProductVector<Scalar> >(
      inArgs.get_x_dot().assert_not_null(), true
      );
  const Scalar
    alpha = inArgs.get_alpha();
  const Scalar
    beta = inArgs.get_beta();

  RCP<const Thyra::MultiVectorBase<Scalar> >
    S = s_bar->getMultiVector();
  RCP<const Thyra::MultiVectorBase<Scalar> >
    S_dot = s_bar_dot->getMultiVector();
  
  //
  // OutArgs
  //

  RCP<Thyra::DefaultMultiVectorProductVector<Scalar> >
    f_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorProductVector<Scalar> >(
      outArgs.get_f(), true
      );
  RCP<Thyra::DefaultMultiVectorLinearOpWithSolve<Scalar> >
    W_sens = rcp_dynamic_cast<Thyra::DefaultMultiVectorLinearOpWithSolve<Scalar> >(
      outArgs.get_W(), true
      );

  RCP<Thyra::MultiVectorBase<Scalar> >
    F_sens = f_sens->getNonconstMultiVector().assert_not_null();

  //
  // Compute the requested functions
  //

  if(!is_null(F_sens)) {

#ifdef ENABLE_RYTHMOS_TIMERS
    TEUCHOS_FUNC_TIME_MONITOR("Rythmos:ForwardSensitivityImplicitModelEvaluator::evalModel: computeSens");
#endif

    // S_diff =  -(coeff_x_dot/coeff_x)*S + S_dot
    RCP<Thyra::MultiVectorBase<Scalar> >
      S_diff = createMembers( stateModel_->get_x_space(), np_ );
    V_StVpV( &*S_diff, Scalar(-coeff_x_dot_/coeff_x_), *S, *S_dot );
    // F_sens = (1/coeff_x) * W_tilde * S
    Thyra::apply(
      *W_tilde_, Thyra::NOTRANS,
      *S, &*F_sens,
      Scalar(1.0/coeff_x_), ST::zero()
      );
    // F_sens += d(f)/d(x_dot) * S_diff
    Thyra::apply(
      *DfDx_dot_, Thyra::NOTRANS,
      *S_diff, &*F_sens,
      ST::one(), ST::one()
      );
    // F_sens += d(f)/d(p)
    Vp_V( &*F_sens, *DfDp_ );
  }
  
  if(!is_null(W_sens)) {
    TEST_FOR_EXCEPTION(
      alpha != coeff_x_dot_, std::logic_error,
      "Error, alpha="<<alpha<<" != coeff_x_dot="<<coeff_x_dot_
      <<" with difference = "<<(alpha-coeff_x_dot_)<<"!" );
    TEST_FOR_EXCEPTION(
      beta != coeff_x_, std::logic_error,
      "Error, beta="<<beta<<" != coeff_x="<<coeff_x_
      <<" with difference = "<<(beta-coeff_x_)<<"!" );
    W_sens->initialize( W_tilde_, s_bar_space_, f_sens_space_ );
    
  }

  THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();

}


// private


template<class Scalar>
void ForwardSensitivityImplicitModelEvaluator<Scalar>::wrapNominalValuesAndBounds()
{

  using Teuchos::rcp_dynamic_cast;
  typedef Thyra::ModelEvaluatorBase MEB;

  // nominalValues_.clear(); // ToDo: Implement this!

  nominalValues_.set_t(stateModel_->getNominalValues().get_t());
  
  // 2007/05/22: rabartl: Note: Currently there is not much of a reason to set
  // an initial condition here since the initial condition for the
  // sensitivities is really being set in the ForwardSensitivityStepper
  // object!  In the future, a more general use of this class might benefit
  // from setting the initial condition here.

}


template<class Scalar>
void ForwardSensitivityImplicitModelEvaluator<Scalar>::computeDerivativeMatrices(
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &point
  ) const
{

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

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

  const Scalar
    t_base = stateBasePoint_.get_t(),
    t = point.get_t();

  TEUCHOS_ASSERT_EQUALITY( t , t_base );

  if (is_null(W_tilde_)) {
    TEST_FOR_EXCEPT_MSG(true, "ToDo: compute W_tilde from scratch!");
  }
  
  if ( is_null(DfDx_dot_) || is_null(DfDp_) ) {

    MEB::InArgs<Scalar> inArgs = stateBasePoint_;
    MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
    
    Teuchos::RCP<Thyra::LinearOpBase<Scalar> > DfDx_dot_compute;
    if (is_null(DfDx_dot_)) {
      DfDx_dot_compute = stateModel_->create_W_op();
      inArgs.set_alpha(1.0);
      inArgs.set_beta(0.0);
      outArgs.set_W_op(DfDx_dot_compute);
    }

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

  TEST_FOR_EXCEPT_MSG( !is_null(stateIntegrator_),
    "ToDo: Update for using the stateIntegrator!" );

/* 2007/12/11: rabartl: ToDo: Update the code below to work for the general
 * case.  It does not work in its current form!
 */

/*

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

  RCP<Teuchos::FancyOStream> out = this->getOStream();
  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
  
  const Scalar t = point.get_t();

  //
  // A) Set up the base point at time t and determine what needs to be
  // computed here.
  //

  bool update_W_tilde = false;
  bool update_DfDx_dot = false;
  bool update_DfDp = false;

  if (!is_null(stateIntegrator_)) {
    // Get x and x_dot at t and flag to recompute all matrices!
    RCP<const Thyra::VectorBase<Scalar> > x, x_dot;
    get_fwd_x_and_x_dot( *stateIntegrator_, t, &x, &x_dot );
    stateBasePoint_.set_t(t);
    stateBasePoint_.set_x(x);
    stateBasePoint_.set_x_dot(x_dot);
    update_W_tilde = true;
    update_DfDx_dot = true;
    update_DfDp = true;
  }
  else {
    // The time point should be the same so only compute matrices that have
    // not already been computed.
    TEUCHOS_ASSERT_EQUALITY( t , stateBasePoint_.get_t() );
    if (is_null(W_tilde_))
      update_W_tilde = true;
    if (is_null(DfDx_dot_))
      update_DfDx_dot = true;
    if (is_null(DfDp_))
      update_DfDx_dot = true;
  }

  //
  // B) Compute DfDx_dot and DfDp at the same time if needed
  //

  if ( update_DfDx_dot || update_DfDp ) {

    // B.1) Allocate storage for matrices.  Note: All of these matrices are
    // needed so any of these that is null must be coputed!

    if ( is_null(DfDx_dot_) || is_null(DfDx_dot_compute_) )
      DfDx_dot_compute_ = stateModel_->create_W_op();

    if ( is_null(DfDp_) || is_null(DfDp_compute_) )
      DfDp_compute_ = Thyra::create_DfDp_mv(
        *stateModel_,p_index_,
        MEB::DERIV_MV_BY_COL
        ).getMultiVector();

    // B.2) Setup the InArgs and OutArgs

    MEB::InArgs<Scalar> inArgs = stateBasePoint_;
    MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
    
    if (update_DfDx_dot) {
      inArgs.set_alpha(1.0);
      inArgs.set_beta(0.0);
      outArgs.set_W_op(DfDx_dot_compute_);
    }
    // 2007/12/07: rabartl: ToDo: I need to change the structure of the
    // derivative properties in OutArgs to keep track of whether DfDx_dot is
    // constant or not separate from W.  Right now I can't do this.  This is
    // one reason to add a new DfDx_dot output object to OutArgs.  A better
    // solution is a more general data structure giving the dependence of f()
    // and g[j]() on x, x_dot, p[l], t, etc ...

    if (update_DfDp) {
      outArgs.set_DfDp(
        p_index_,
        MEB::Derivative<Scalar>(DfDp_compute_, MEB::DERIV_MV_BY_COL)
        );
    }

    // B.3) Evaluate the outputs
    
    VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
    stateModel_->evalModel(inArgs,outArgs);

    // B.4) Set the outputs

    if (!is_null(DfDx_dot_compute_))
      DfDx_dot_ = DfDx_dot_compute_;

    if (!is_null(DfDp_compute_))
      DfDp_ = DfDp_compute_;
  
  }

  //
  // C) Compute W_tilde separately if needed
  //

  if ( update_W_tilde ) {

    // C.1) Allocate storage for matrices.  Note: All of these matrices are
    // needed so any of these that is null must be coputed!

    if ( is_null(W_tilde_) || is_null(W_tilde_compute_) )
      W_tilde_compute_ = stateModel_->create_W();

    // C.2) Setup the InArgs and OutArgs

    MEB::InArgs<Scalar> inArgs = stateBasePoint_;
    MEB::OutArgs<Scalar> outArgs = stateModel_->createOutArgs();
    
    if (is_null(W_tilde_)) {
      coeff_x_dot_ = point.get_alpha();
      coeff_x_ = point.get_beta();
      inArgs.set_alpha(coeff_x_dot_);
      inArgs.set_beta(coeff_x_);
      outArgs.set_W(W_tilde_compute_);
    }

    // C.3) Evaluate the outputs
    
    VOTSME stateModel_outputTempState(stateModel_,out,verbLevel);
    stateModel_->evalModel(inArgs,outArgs);

    // B.4) Set the outputs

    if (!is_null(W_tilde_compute_))
      W_tilde_ = W_tilde_compute_;

  }

  // 2007/12/07: rabartl: Note: Above, we see an example of where it would be
  // good to have a separate OutArg for DfDx_dot (as a LOB object) so that we
  // can compute W and DfDx_dot (and any other output quantitiy) at the same
  // time.

*/


}


} // namespace Rythmos


#endif // RYTHMOS_FORWARD_SENSITIVITY_IMPLICIT_MODEL_EVALUATOR_HPP

---