C++ modelevaluatorbase::InArgs类代码示例

void
Piro::LOCASolver<Scalar>::evalModelImpl(
    const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
    const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs) const
{
  const int l = 0; // TODO: Allow user to select parameter index
  const Teuchos::RCP<const Thyra::VectorBase<Scalar> > p_inargs = inArgs.get_p(l);

  // Forward parameter values to the LOCA stepper
  {
    const Teuchos::RCP<const Thyra::VectorBase<Scalar> > p_inargs_or_nominal =
      Teuchos::nonnull(p_inargs) ? p_inargs : this->getNominalValues().get_p(l);
    const Thyra::ConstDetachedVectorView<Scalar> p_init_values(p_inargs_or_nominal);
    const Teuchos_Ordinal p_entry_count = p_init_values.subDim();
    TEUCHOS_ASSERT(p_entry_count == Teuchos::as<Teuchos_Ordinal>(paramVector_.length()));

    for (Teuchos_Ordinal k = 0; k < p_entry_count; ++k) {
      paramVector_[k] = p_init_values[k];
    }

    group_->setParams(paramVector_);
  }

  stepper_->reset(globalData_, group_, locaStatusTests_, noxStatusTests_, piroParams_);
  const LOCA::Abstract::Iterator::IteratorStatus status = stepper_->run();

  if (status == LOCA::Abstract::Iterator::Finished) {
    std::cerr << "Continuation Stepper Finished.\n";
  } else if (status == LOCA::Abstract::Iterator::NotFinished) {
    std::cerr << "Continuation Stepper did not reach final value.\n";
  } else {
    std::cerr << "Nonlinear solver failed to converge.\n";
    outArgs.setFailed();
  }

  const Teuchos::RCP<Thyra::VectorBase<Scalar> > x_outargs = outArgs.get_g(this->num_g());
  const Teuchos::RCP<Thyra::VectorBase<Scalar> > x_final =
    Teuchos::nonnull(x_outargs) ? x_outargs : Thyra::createMember(this->get_g_space(this->num_g()));

  {
    // Deep copy final solution from LOCA group
    NOX::Thyra::Vector finalSolution(x_final);
    finalSolution = group_->getX();
  }

  // Compute responses for the final solution
  {
    Thyra::ModelEvaluatorBase::InArgs<Scalar> modelInArgs =
      this->getModel().createInArgs();
    {
      modelInArgs.set_x(x_final);
      modelInArgs.set_p(l, p_inargs);
    }

    this->evalConvergedModel(modelInArgs, outArgs);
  }
}

Thyra::ModelEvaluatorBase::InArgs<Scalar>
Piro::VelocityVerletSolver<Scalar>::getNominalValues() const
{
  Thyra::ModelEvaluatorBase::InArgs<Scalar> result = this->createInArgs();
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> modelNominalValues = model->getNominalValues();
  for (int l = 0; l < num_p; ++l) {
    result.set_p(l, modelNominalValues.get_p(l));
  }
  return result;
}

Thyra::ModelEvaluatorBase::InArgs<Scalar>
Piro::SteadyStateSolver<Scalar>::getNominalValues() const
{
  Thyra::ModelEvaluatorBase::InArgs<Scalar> result = this->createInArgsImpl();
  result.setArgs(
      model_->getNominalValues(),
      /* ignoreUnsupported = */ true,
      /* cloneObjects = */ false);
  return result;
}

Thyra::ModelEvaluatorBase::InArgs<Scalar>
DiagonalROME<Scalar>::getNominalValues() const
{
    Thyra::ModelEvaluatorBase::InArgs<Scalar> initialGuess =
        this->createInArgs();
    RCP<Thyra::VectorBase<Scalar> > p_init =
        Thyra::createMember<Scalar>(p_space_);
    Thyra::V_S( p_init.ptr(), 1.5 );
    initialGuess.set_p(0, p_init);
    return initialGuess;
}

void ImplicitRKStepper<Scalar>::setInitialCondition(
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initialCondition
  )
{

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

  basePoint_ = initialCondition;

  // x

  RCP<const Thyra::VectorBase<Scalar> >
    x_init = initialCondition.get_x();

#ifdef HAVE_RYTHMOS_DEBUG
  TEUCHOS_TEST_FOR_EXCEPTION(
    is_null(x_init), std::logic_error,
    "Error, if the client passes in an intial condition to setInitialCondition(...),\n"
    "then x can not be null!" );
#endif

  x_ = x_init->clone_v();

  // x_dot

  x_dot_ = createMember(x_->space());

  RCP<const Thyra::VectorBase<Scalar> >
    x_dot_init = initialCondition.get_x_dot();

  if (!is_null(x_dot_init))
    assign(x_dot_.ptr(),*x_dot_init);
  else
    assign(x_dot_.ptr(),ST::zero());
  
  // t

  const Scalar t =
    (
      initialCondition.supports(MEB::IN_ARG_t)
      ? initialCondition.get_t()
      : ST::zero()
      );

  timeRange_ = timeRange(t,t);

  // x_old
  x_old_ = x_->clone_v();

  haveInitialCondition_ = true;

}

TEUCHOS_UNIT_TEST( Rythmos_ExplicitRKStepper, basePoint ) {
  RCP<SinCosModel> model = sinCosModel(false);
  {
    RCP<ParameterList> pl = Teuchos::parameterList();
    pl->set("Accept model parameters",true);
    model->setParameterList(pl);
  }
  Thyra::ModelEvaluatorBase::InArgs<double> ic = model->getNominalValues();
  // t_ic
  double t_ic = 1.0; // not used
  // x_ic
  RCP<VectorBase<double> > x_ic = Thyra::createMember(*model->get_x_space());
  {
    Thyra::DetachedVectorView<double> x_ic_view( *x_ic );
    x_ic_view[0] = 5.0;
    x_ic_view[1] = 6.0;
  }
  // parameter 0 ic
  RCP<VectorBase<double> > p_ic = Thyra::createMember(*model->get_p_space(0));
  {
    Thyra::DetachedVectorView<double> p_ic_view( *p_ic );
    p_ic_view[0] = 2.0; // a
    p_ic_view[1] = 3.0; // f 
    p_ic_view[2] = 4.0; // L
  }
  ic.set_p(0,p_ic); 
  ic.set_x(x_ic);
  ic.set_t(t_ic);
  RCP<ExplicitRKStepper<double> > stepper = explicitRKStepper<double>();
  stepper->setModel(model);
  stepper->setInitialCondition(ic);
  stepper->setRKButcherTableau(createRKBT<double>("Forward Euler"));
  double dt = 0.2;
  double dt_taken;
  dt_taken = stepper->takeStep(dt,STEP_TYPE_FIXED);
  TEST_EQUALITY_CONST( dt_taken, 0.2 );
  const StepStatus<double> status = stepper->getStepStatus();
  TEST_ASSERT( !is_null(status.solution) );
  double tol = 1.0e-10;
  {
    Thyra::ConstDetachedVectorView<double> x_new_view( *(status.solution) );
    TEST_FLOATING_EQUALITY( x_new_view[0], 5.0 + 0.2*(6.0), tol );
    TEST_FLOATING_EQUALITY( x_new_view[1], 6.0 + 0.2*( (3.0/4.0)*(3.0/4.0)*(2.0-5.0) ), tol );
  }
}

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

    using Teuchos::as;
    using Teuchos::outArg;
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Thyra::get_mv;
    using Thyra::ConstDetachedSpmdVectorView;
    using Thyra::DetachedSpmdVectorView;
    typedef Thyra::Ordinal Ordinal;
    typedef Thyra::ModelEvaluatorBase MEB;
    typedef MEB::DerivativeMultiVector<Scalar> DMV;

    const ConstDetachedSpmdVectorView<Scalar> p(inArgs.get_p(0));
    const ConstDetachedSpmdVectorView<Scalar> ps(ps_);
    const ConstDetachedSpmdVectorView<Scalar> diag(diag_);
    const ConstDetachedSpmdVectorView<Scalar> s_bar(s_bar_);

    // g(p)
    if (!is_null(outArgs.get_g(0))) {
        Scalar g_val = ST::zero();
        for (Ordinal i = 0; i < p.subDim(); ++i) {
            const Scalar p_ps = p[i] - ps[i];
            g_val += diag[i] * p_ps*p_ps;
            if (nonlinearTermFactor_ != ST::zero()) {
                g_val += nonlinearTermFactor_ * p_ps * p_ps * p_ps;
            }
        }
        Scalar global_g_val;
        Teuchos::reduceAll<Ordinal, Scalar>(*comm_, Teuchos::REDUCE_SUM,
                                            g_val, outArg(global_g_val) );
        DetachedSpmdVectorView<Scalar>(outArgs.get_g(0))[0] =
            as<Scalar>(0.5) * global_g_val + g_offset_;
    }

    // DgDp[i]
    if (!outArgs.get_DgDp(0,0).isEmpty()) {
        const RCP<Thyra::MultiVectorBase<Scalar> > DgDp_trans_mv =
            get_mv<Scalar>(outArgs.get_DgDp(0,0), "DgDp^T", MEB::DERIV_TRANS_MV_BY_ROW);
        const DetachedSpmdVectorView<Scalar> DgDp_grad(DgDp_trans_mv->col(0));
        for (Thyra::Ordinal i = 0; i < p.subDim(); ++i) {
            const Scalar p_ps = p[i] - ps[i];
            Scalar DgDp_grad_i = diag[i] * p_ps;
            if (nonlinearTermFactor_ != ST::zero()) {
                DgDp_grad_i += as<Scalar>(1.5) * nonlinearTermFactor_ * p_ps * p_ps;
            }
            DgDp_grad[i] = DgDp_grad_i / s_bar[i];

        }
    }

}

void TimeDiscretizedBackwardEulerModelEvaluator<Scalar>::initialize(
  const RCP<const Thyra::ModelEvaluator<Scalar> > &daeModel,
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &initCond,
  const Scalar finalTime,
  const int numTimeSteps,
  const RCP<Thyra::LinearOpWithSolveFactoryBase<Scalar> > &W_bar_factory
  )
{

  TEST_FOR_EXCEPT(is_null(daeModel));
  TEST_FOR_EXCEPT(is_null(initCond.get_x()));
  TEST_FOR_EXCEPT(is_null(initCond.get_x_dot()));
  TEST_FOR_EXCEPT(finalTime <= initCond.get_t());
  TEST_FOR_EXCEPT(numTimeSteps <= 0);
  // ToDo: Validate that daeModel is of the right form!

  daeModel_ = daeModel;
  initCond_ = initCond;
  finalTime_ = finalTime;
  numTimeSteps_ = numTimeSteps;

  initTime_ = initCond.get_t();
  delta_t_ = (finalTime_ - initTime_) / numTimeSteps_;

  x_bar_space_ = productVectorSpace(daeModel_->get_x_space(),numTimeSteps_);
  f_bar_space_ = productVectorSpace(daeModel_->get_f_space(),numTimeSteps_);

  if (!is_null(W_bar_factory)) {
    W_bar_factory_ = W_bar_factory;
  }
  else {
    W_bar_factory_ =
      Thyra::defaultBlockedTriangularLinearOpWithSolveFactory<Scalar>(
        daeModel_->get_W_factory()
        );
  }
  
}

void Simple2DModelEvaluator<Scalar>::evalModelImpl(
  const Thyra::ModelEvaluatorBase::InArgs<Scalar> &inArgs,
  const Thyra::ModelEvaluatorBase::OutArgs<Scalar> &outArgs
  ) const
{
  using Teuchos::rcp_dynamic_cast;
  const Scalar one = 1.0, two = 2.0, zero = 0.0;

  const ConstDetachedVectorView<Scalar> x(inArgs.get_x());

  const RCP<Thyra::VectorBase<Scalar> > f_out = outArgs.get_f();
  const RCP<Thyra::LinearOpBase< Scalar > > W_op_out = outArgs.get_W_op();
  const RCP<Thyra::PreconditionerBase< Scalar > > W_prec_out = outArgs.get_W_prec();

  if (nonnull(f_out)) {
    const DetachedVectorView<Scalar> f(f_out);
    f[0] = x[0] + x[1] * x[1] - p_[0];
    f[1] = d_ * (x[0] * x[0] - x[1] - p_[1]);
  }

  if (nonnull(W_op_out)) {
    const RCP<SimpleDenseLinearOp<Scalar> > W =
      rcp_dynamic_cast<SimpleDenseLinearOp<Scalar> >(W_op_out, true);
    const RCP<MultiVectorBase<Scalar> > W_mv = W->getNonconstMultiVector();
    Thyra::DetachedMultiVectorView<Scalar> W_dmvv(W_mv);
    W_dmvv(0, 0) = one;
    W_dmvv(0, 1) = two * x[1];
    W_dmvv(1, 0) = d_ * two * x[0];
    W_dmvv(1, 1) = -d_;
  }

  if (nonnull(W_prec_out)) {
    const RCP<SimpleDenseLinearOp<Scalar> > W_prec_op =
      rcp_dynamic_cast<SimpleDenseLinearOp<Scalar> >(
        W_prec_out->getNonconstUnspecifiedPrecOp(), true);
    const RCP<MultiVectorBase<Scalar> > W_prec_mv = W_prec_op->getNonconstMultiVector();
    Thyra::DetachedMultiVectorView<Scalar> W_prec_dmvv(W_prec_mv);
    // Diagonal inverse of W (see W above)
    W_prec_dmvv(0, 0) = one;
    W_prec_dmvv(0, 1) = zero;
    W_prec_dmvv(1, 0) = zero;
    W_prec_dmvv(1, 1) = -one/d_;
  }
  
}

void
Albany::ModelEvaluatorT::evalModelImpl(
    const Thyra::ModelEvaluatorBase::InArgs<ST>& inArgsT,
    const Thyra::ModelEvaluatorBase::OutArgs<ST>& outArgsT) const
{

  Teuchos::TimeMonitor Timer(*timer); //start timer
  //
  // Get the input arguments
  //
  const Teuchos::RCP<const Tpetra_Vector> xT =
    ConverterT::getConstTpetraVector(inArgsT.get_x());

  const Teuchos::RCP<const Tpetra_Vector> x_dotT =
    Teuchos::nonnull(inArgsT.get_x_dot()) ?
    ConverterT::getConstTpetraVector(inArgsT.get_x_dot()) :
    Teuchos::null;

  // AGS: x_dotdot time integrators not imlemented in Thyra ME yet
  //const Teuchos::RCP<const Tpetra_Vector> x_dotdotT =
  //  Teuchos::nonnull(inArgsT.get_x_dotdot()) ?
  //  ConverterT::getConstTpetraVector(inArgsT.get_x_dotdot()) :
  //  Teuchos::null;
  const Teuchos::RCP<const Tpetra_Vector> x_dotdotT = Teuchos::null;


  const double alpha = (Teuchos::nonnull(x_dotT) || Teuchos::nonnull(x_dotdotT)) ? inArgsT.get_alpha() : 0.0;
  // AGS: x_dotdot time integrators not imlemented in Thyra ME yet
  // const double omega = (Teuchos::nonnull(x_dotT) || Teuchos::nonnull(x_dotdotT)) ? inArgsT.get_omega() : 0.0;
  const double omega = 0.0;
  const double beta = (Teuchos::nonnull(x_dotT) || Teuchos::nonnull(x_dotdotT)) ? inArgsT.get_beta() : 1.0;
  const double curr_time = (Teuchos::nonnull(x_dotT) || Teuchos::nonnull(x_dotdotT)) ? inArgsT.get_t() : 0.0;

  for (int l = 0; l < inArgsT.Np(); ++l) {
    const Teuchos::RCP<const Thyra::VectorBase<ST> > p = inArgsT.get_p(l);
    if (Teuchos::nonnull(p)) {
      const Teuchos::RCP<const Tpetra_Vector> pT = ConverterT::getConstTpetraVector(p);
      const Teuchos::ArrayRCP<const ST> pT_constView = pT->get1dView();

      ParamVec &sacado_param_vector = sacado_param_vec[l];
      for (unsigned int k = 0; k < sacado_param_vector.size(); ++k) {
        sacado_param_vector[k].baseValue = pT_constView[k];
      }
    }
  }

  //
  // Get the output arguments
  //
  const Teuchos::RCP<Tpetra_Vector> fT_out =
    Teuchos::nonnull(outArgsT.get_f()) ?
    ConverterT::getTpetraVector(outArgsT.get_f()) :
    Teuchos::null;

  const Teuchos::RCP<Tpetra_Operator> W_op_outT =
    Teuchos::nonnull(outArgsT.get_W_op()) ?
    ConverterT::getTpetraOperator(outArgsT.get_W_op()) :
    Teuchos::null;

  // Cast W to a CrsMatrix, throw an exception if this fails
  const Teuchos::RCP<Tpetra_CrsMatrix> W_op_out_crsT =
    Teuchos::nonnull(W_op_outT) ?
    Teuchos::rcp_dynamic_cast<Tpetra_CrsMatrix>(W_op_outT, true) :
    Teuchos::null;

  //
  // Compute the functions
  //
  bool f_already_computed = false;

  // W matrix
  if (Teuchos::nonnull(W_op_out_crsT)) {
    app->computeGlobalJacobianT(
        alpha, beta, omega, curr_time, x_dotT.get(), x_dotdotT.get(),  *xT,
        sacado_param_vec, fT_out.get(), *W_op_out_crsT);
    f_already_computed = true;
  }

  // df/dp
  for (int l = 0; l < outArgsT.Np(); ++l) {
    const Teuchos::RCP<Thyra::MultiVectorBase<ST> > dfdp_out =
      outArgsT.get_DfDp(l).getMultiVector();

    const Teuchos::RCP<Tpetra_MultiVector> dfdp_outT =
      Teuchos::nonnull(dfdp_out) ?
      ConverterT::getTpetraMultiVector(dfdp_out) :
      Teuchos::null;

    if (Teuchos::nonnull(dfdp_outT)) {
      const Teuchos::RCP<ParamVec> p_vec = Teuchos::rcpFromRef(sacado_param_vec[l]);

      app->computeGlobalTangentT(
          0.0, 0.0, 0.0, curr_time, false, x_dotT.get(), x_dotdotT.get(), *xT,
          sacado_param_vec, p_vec.get(),
          NULL, NULL, NULL, NULL, fT_out.get(), NULL,
          dfdp_outT.get());

      f_already_computed = true;
    }
  }
//.........这里部分代码省略.........

void
Albany::ModelEvaluatorT::evalModelImpl(
    const Thyra::ModelEvaluatorBase::InArgs<ST>& inArgsT,
    const Thyra::ModelEvaluatorBase::OutArgs<ST>& outArgsT) const
{

  #ifdef OUTPUT_TO_SCREEN
    std::cout << "DEBUG: " << __PRETTY_FUNCTION__ << "\n";
  #endif

  Teuchos::TimeMonitor Timer(*timer); //start timer
  //
  // Get the input arguments
  //
  const Teuchos::RCP<const Tpetra_Vector> xT =
    ConverterT::getConstTpetraVector(inArgsT.get_x());

  const Teuchos::RCP<const Tpetra_Vector> x_dotT =
    (supports_xdot && Teuchos::nonnull(inArgsT.get_x_dot())) ?
    ConverterT::getConstTpetraVector(inArgsT.get_x_dot()) :
    Teuchos::null;


  const Teuchos::RCP<const Tpetra_Vector> x_dotdotT =
    (supports_xdotdot && Teuchos::nonnull(this->get_x_dotdot())) ?
    ConverterT::getConstTpetraVector(this->get_x_dotdot()) :
    Teuchos::null;

  const double alpha = (Teuchos::nonnull(x_dotT) || Teuchos::nonnull(x_dotdotT)) ? inArgsT.get_alpha() : 0.0;
  const double omega = Teuchos::nonnull(x_dotdotT) ? this->get_omega() : 0.0;
  const double beta = (Teuchos::nonnull(x_dotT) || Teuchos::nonnull(x_dotdotT)) ? inArgsT.get_beta() : 1.0;
  const double curr_time = (Teuchos::nonnull(x_dotT) || Teuchos::nonnull(x_dotdotT)) ? inArgsT.get_t() : 0.0;

  for (int l = 0; l < inArgsT.Np(); ++l) {
    const Teuchos::RCP<const Thyra::VectorBase<ST> > p = inArgsT.get_p(l);
    if (Teuchos::nonnull(p)) {
      const Teuchos::RCP<const Tpetra_Vector> pT = ConverterT::getConstTpetraVector(p);
      const Teuchos::ArrayRCP<const ST> pT_constView = pT->get1dView();

      ParamVec &sacado_param_vector = sacado_param_vec[l];
      for (unsigned int k = 0; k < sacado_param_vector.size(); ++k) {
        sacado_param_vector[k].baseValue = pT_constView[k];
      }
    }
  }

  //
  // Get the output arguments
  //
  const Teuchos::RCP<Tpetra_Vector> fT_out =
    Teuchos::nonnull(outArgsT.get_f()) ?
    ConverterT::getTpetraVector(outArgsT.get_f()) :
    Teuchos::null;

  const Teuchos::RCP<Tpetra_Operator> W_op_outT =
    Teuchos::nonnull(outArgsT.get_W_op()) ?
    ConverterT::getTpetraOperator(outArgsT.get_W_op()) :
    Teuchos::null;

#ifdef WRITE_MASS_MATRIX_TO_MM_FILE
  //IK, 4/24/15: adding object to hold mass matrix to be written to matrix market file
  const Teuchos::RCP<Tpetra_Operator> Mass =
    Teuchos::nonnull(outArgsT.get_W_op()) ?
    ConverterT::getTpetraOperator(outArgsT.get_W_op()) :
    Teuchos::null;
  //IK, 4/24/15: needed for writing mass matrix out to matrix market file
  const Teuchos::RCP<Tpetra_Vector> ftmp =
    Teuchos::nonnull(outArgsT.get_f()) ?
    ConverterT::getTpetraVector(outArgsT.get_f()) :
    Teuchos::null;
#endif

  // Cast W to a CrsMatrix, throw an exception if this fails
  const Teuchos::RCP<Tpetra_CrsMatrix> W_op_out_crsT =
    Teuchos::nonnull(W_op_outT) ?
    Teuchos::rcp_dynamic_cast<Tpetra_CrsMatrix>(W_op_outT, true) :
    Teuchos::null;

#ifdef WRITE_MASS_MATRIX_TO_MM_FILE
  //IK, 4/24/15: adding object to hold mass matrix to be written to matrix market file
  const Teuchos::RCP<Tpetra_CrsMatrix> Mass_crs =
    Teuchos::nonnull(Mass) ?
    Teuchos::rcp_dynamic_cast<Tpetra_CrsMatrix>(Mass, true) :
    Teuchos::null;
#endif

  //
  // Compute the functions
  //
  bool f_already_computed = false;

  // W matrix
  if (Teuchos::nonnull(W_op_out_crsT)) {
    app->computeGlobalJacobianT(
        alpha, beta, omega, curr_time, x_dotT.get(), x_dotdotT.get(),  *xT,
        sacado_param_vec, fT_out.get(), *W_op_out_crsT);
    f_already_computed = true;
#ifdef WRITE_MASS_MATRIX_TO_MM_FILE
    //IK, 4/24/15: write mass matrix to matrix market file
    //Warning: to read this in to MATLAB correctly, code must be run in serial.
//.........这里部分代码省略.........

//.........这里部分代码省略.........
#endif
#endif
   
   //set previousSolution (used as initial guess for next time step) to final Albany solution. 
   previousSolution = Teuchos::rcp(new Tpetra_Vector(*albanyApp->getDiscretization()->getSolutionFieldT())); 
   nElementsActivePrevious = nElementsActive;   
 
   //std::cout << "Final solution: " << *albanyApp->getDiscretization()->getSolutionField() << std::endl;  
    // ---------------------------------------------------------------------------------------------------
    // Compute sensitivies / responses and perform regression tests
    // IK, 12/9/13: how come this is turned off in mpas branch? 
    // ---------------------------------------------------------------------------------------------------
 
    if (debug_output_verbosity != 0 & mpiCommT->getRank() == 0) 
      std::cout << "Computing responses and sensitivities..." << std::endl;
    int status=0; // 0 = pass, failures are incremented
#ifdef CISM_USE_EPETRA
    Teuchos::Array<Teuchos::RCP<const Epetra_Vector> > responses;
    Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Epetra_MultiVector> > > sensitivities;
    epetraFromThyra(mpiComm, thyraResponses, thyraSensitivities, responses, sensitivities);
#else
    Teuchos::Array<Teuchos::RCP<const Tpetra_Vector> > responses;
    Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Tpetra_MultiVector> > > sensitivities;
    tpetraFromThyra(thyraResponses, thyraSensitivities, responses, sensitivities);
#endif

    const int num_p = solver->Np(); // Number of *vectors* of parameters
    const int num_g = solver->Ng(); // Number of *vectors* of responses

   if (debug_output_verbosity != 0) {
    *out << "Finished eval of first model: Params, Responses "
      << std::setprecision(12) << std::endl;
   }
   const Thyra::ModelEvaluatorBase::InArgs<double> nominal = solver->getNominalValues();

   if (debug_output_verbosity != 0) {
    for (int i=0; i<num_p; i++) {
#ifdef CISM_USE_EPETRA
      const Teuchos::RCP<const Epetra_Vector> p_init = epetraVectorFromThyra(mpiComm, nominal.get_p(i));
      p_init->Print(*out << "\nParameter vector " << i << ":\n");
#else
      Albany::printTpetraVector(*out << "\nParameter vector " << i << ":\n",
           ConverterT::getConstTpetraVector(nominal.get_p(i)));
#endif
    }
   }

    for (int i=0; i<num_g-1; i++) {
#ifdef CISM_USE_EPETRA
      const Teuchos::RCP<const Epetra_Vector> g = responses[i];
#else
      const Teuchos::RCP<const Tpetra_Vector> g = responses[i];
#endif
      bool is_scalar = true;

      if (albanyApp != Teuchos::null)
        is_scalar = albanyApp->getResponse(i)->isScalarResponse();

      if (is_scalar) {
        if (debug_output_verbosity != 0) {
#ifdef CISM_USE_EPETRA
         g->Print(*out << "\nResponse vector " << i << ":\n");
#else
         Albany::printTpetraVector(*out << "\nResponse vector " << i << ":\n", g);
#endif
        }

// hide the original parental method AMET->evalModelImpl():
void
Aeras::HVDecorator::evalModelImpl(
    const Thyra::ModelEvaluatorBase::InArgs<ST>& inArgsT,
    const Thyra::ModelEvaluatorBase::OutArgs<ST>& outArgsT) const
{
#ifdef OUTPUT_TO_SCREEN
  std::cout << "DEBUG WHICH HVDecorator: " << __PRETTY_FUNCTION__ << "\n";
#endif
	
  Teuchos::TimeMonitor Timer(*timer); //start timer

  //
  // Get the input arguments
  //
  // Thyra vectors
  const Teuchos::RCP<const Thyra_Vector> x = inArgsT.get_x();
  const Teuchos::RCP<const Thyra_Vector> x_dot =
      (supports_xdot ? inArgsT.get_x_dot() : Teuchos::null);
  const Teuchos::RCP<const Thyra_Vector> x_dotdot =
      (supports_xdotdot ? inArgsT.get_x_dot_dot() : Teuchos::null);

  const double alpha = (Teuchos::nonnull(x_dot) || Teuchos::nonnull(x_dotdot)) ? inArgsT.get_alpha() : 0.0;
  // AGS: x_dotdot time integrators not imlemented in Thyra ME yet
  // const double omega = (Teuchos::nonnull(x_dot) || Teuchos::nonnull(x_dotdot)) ? inArgsT.get_omega() : 0.0;
  const double omega = 0.0;
  const double beta = (Teuchos::nonnull(x_dot) || Teuchos::nonnull(x_dotdot)) ? inArgsT.get_beta() : 1.0;
  const double curr_time = (Teuchos::nonnull(x_dot) || Teuchos::nonnull(x_dotdot)) ? inArgsT.get_t() : 0.0;

  for (int l = 0; l < inArgsT.Np(); ++l) {
    const Teuchos::RCP<const Thyra_Vector> p = inArgsT.get_p(l);
    if (Teuchos::nonnull(p)) {
      const Teuchos::RCP<const Tpetra_Vector> pT = Albany::getConstTpetraVector(p);
      const Teuchos::ArrayRCP<const ST> pT_constView = pT->get1dView();

      ParamVec &sacado_param_vector = sacado_param_vec[l];
      for (unsigned int k = 0; k < sacado_param_vector.size(); ++k) {
        sacado_param_vector[k].baseValue = pT_constView[k];
      }
    }
  }

  //
  // Get the output arguments
  //
  auto f    = outArgsT.get_f();
  auto W_op = outArgsT.get_W_op();

  //
  // Compute the functions
  //
  bool f_already_computed = false;

  // W matrix
  if (Teuchos::nonnull(W_op)) {
    app->computeGlobalJacobian(
        alpha, beta, omega, curr_time, x, x_dot, x_dotdot,
        sacado_param_vec, f, W_op);
    f_already_computed = true;
  }

  // df/dp
  for (int l = 0; l < outArgsT.Np(); ++l) {
    const Teuchos::RCP<Thyra_MultiVector> df_dp = outArgsT.get_DfDp(l).getMultiVector();

    if (Teuchos::nonnull(df_dp)) {
      const Teuchos::RCP<ParamVec> p_vec = Teuchos::rcpFromRef(sacado_param_vec[l]);

      app->computeGlobalTangent(
          0.0, 0.0, 0.0, curr_time, false, x, x_dot, x_dotdot,
          sacado_param_vec, p_vec.get(),
          Teuchos::null, Teuchos::null, Teuchos::null, Teuchos::null,
          f, Teuchos::null, df_dp);

      f_already_computed = true;
    }
  }

  // f
  if (app->is_adjoint) {
    const Thyra_Derivative f_deriv(f, Thyra::ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW);
    const Thyra_Derivative dummy_deriv;

    const int response_index = 0; // need to add capability for sending this in
    app->evaluateResponseDerivative(
        response_index, curr_time, x, x_dot, x_dotdot,
        sacado_param_vec, NULL,
        Teuchos::null, f_deriv, dummy_deriv, dummy_deriv, dummy_deriv);
  } else {
    if (Teuchos::nonnull(f) && !f_already_computed) {
      app->computeGlobalResidual(
          curr_time, x, x_dot, x_dotdot,
          sacado_param_vec, f);
    }
  }

  //compute xtilde 
  applyLinvML(x, xtilde); 

#ifdef WRITE_TO_MATRIX_MARKET_TO_MM_FILE
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {

  int status=0; // 0 = pass, failures are incremented
  bool success = true;

#ifdef ALBANY_DEBUG
  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
#else // bypass printing process startup info
  Teuchos::GlobalMPISession mpiSession(&argc, &argv, NULL);
#endif

  Kokkos::initialize(argc, argv);

#ifdef ALBANY_FLUSH_DENORMALS
  _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
  _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);
#endif

#ifdef ALBANY_CHECK_FPE
   // Catch FPEs. Follow Main_SolveT.cpp's approach to checking for floating
   // point exceptions.
   //_mm_setcsr(_MM_MASK_MASK &~ (_MM_MASK_OVERFLOW | _MM_MASK_INVALID | _MM_MASK_DIV_ZERO) );
   _MM_SET_EXCEPTION_MASK(_MM_GET_EXCEPTION_MASK() & ~_MM_MASK_INVALID);
#endif

  using Teuchos::RCP;
  using Teuchos::rcp;

  RCP<Teuchos::FancyOStream> out(Teuchos::VerboseObjectBase::getDefaultOStream());

  // Command-line argument for input file
  Albany::CmdLineArgs cmd;
  cmd.parse_cmdline(argc, argv, *out);

  try {

    RCP<Teuchos::Time> totalTime =
      Teuchos::TimeMonitor::getNewTimer("Albany: ***Total Time***");

    RCP<Teuchos::Time> setupTime =
      Teuchos::TimeMonitor::getNewTimer("Albany: Setup Time");
    Teuchos::TimeMonitor totalTimer(*totalTime); //start timer
    Teuchos::TimeMonitor setupTimer(*setupTime); //start timer

    RCP<const Teuchos_Comm> comm =
      Tpetra::DefaultPlatform::getDefaultPlatform().getComm();

    // Connect vtune for performance profiling
    if (cmd.vtune) {
      Albany::connect_vtune(comm->getRank());
    }

    Albany::SolverFactory slvrfctry(cmd.xml_filename, comm);
    RCP<Epetra_Comm> appComm = Albany::createEpetraCommFromTeuchosComm(comm);
    RCP<Albany::Application> app;
    const RCP<Thyra::ModelEvaluator<double> > solver =
      slvrfctry.createThyraSolverAndGetAlbanyApp(app, appComm, appComm);

    setupTimer.~TimeMonitor();

//    PHX::InitializeKokkosDevice();
   
    Teuchos::ParameterList &solveParams =
      slvrfctry.getAnalysisParameters().sublist("Solve", /*mustAlreadyExist =*/ false);
    // By default, request the sensitivities if not explicitly disabled
    solveParams.get("Compute Sensitivities", true);

    Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<double> > > thyraResponses;
    Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<double> > > > thyraSensitivities;

       // The PoissonSchrodinger_SchroPo and PoissonSchroMosCap1D tests seg fault as albanyApp is null -
       // For now, do not resize the response vectors. FIXME sort out this issue.
    if(Teuchos::nonnull(app))
      Piro::PerformSolveBase(*solver, solveParams, thyraResponses, thyraSensitivities, app->getAdaptSolMgr()->getSolObserver());
    else
      Piro::PerformSolveBase(*solver, solveParams, thyraResponses, thyraSensitivities);

    Teuchos::Array<Teuchos::RCP<const Epetra_Vector> > responses;
    Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Epetra_MultiVector> > > sensitivities;
    epetraFromThyra(appComm, thyraResponses, thyraSensitivities, responses, sensitivities);

    const int num_p = solver->Np(); // Number of *vectors* of parameters
    const int num_g = solver->Ng(); // Number of *vectors* of responses

    *out << "Finished eval of first model: Params, Responses "
      << std::setprecision(12) << std::endl;

    Teuchos::ParameterList& parameterParams = slvrfctry.getParameters().sublist("Problem").sublist("Parameters");
    int num_param_vecs = (parameterParams.isType<int>("Number")) ?
        int(parameterParams.get("Number", 0) > 0) :
        parameterParams.get("Number of Parameter Vectors", 0);

    const Thyra::ModelEvaluatorBase::InArgs<double> nominal = solver->getNominalValues();
    double norm2;
    for (int i=0; i<num_p; i++) {
      const Teuchos::RCP<const Epetra_Vector> p_init = epetraVectorFromThyra(appComm, nominal.get_p(i));
      if(i < num_param_vecs)
        p_init->Print(*out << "\nParameter vector " << i << ":\n");
      else { //distributed parameters, we print only 2-norm
        p_init->Norm2(&norm2);
//.........这里部分代码省略.........

TEUCHOS_UNIT_TEST( Rythmos_GlobalErrorEstimator, SinCos ) {
  typedef Teuchos::ScalarTraits<double> ST;
  // Forward Solve, storing data in linear interpolation buffer
  int storageLimit = 100;
  double finalTime = 0.1;
  double dt = 0.1;
  RCP<IntegratorBuilder<double> > ib = integratorBuilder<double>();
  {
    RCP<ParameterList> ibPL = Teuchos::parameterList();
    ibPL->sublist("Integrator Settings").sublist("Integrator Selection").set("Integrator Type","Default Integrator");
    ibPL->sublist("Integrator Settings").set("Final Time",finalTime);
    ibPL->sublist("Integration Control Strategy Selection").set("Integration Control Strategy Type","Simple Integration Control Strategy");
    ibPL->sublist("Integration Control Strategy Selection").sublist("Simple Integration Control Strategy").set("Take Variable Steps",false);
    ibPL->sublist("Integration Control Strategy Selection").sublist("Simple Integration Control Strategy").set("Fixed dt",dt);

    ibPL->sublist("Stepper Settings").sublist("Stepper Selection").set("Stepper Type","Backward Euler");
    //ibPL->sublist("Stepper Settings").sublist("Stepper Selection").set("Stepper Type","Implicit RK");
    //ibPL->sublist("Stepper Settings").sublist("Runge Kutta Butcher Tableau Selection").set("Runge Kutta Butcher Tableau Type","Backward Euler");
    ibPL->sublist("Interpolation Buffer Settings").sublist("Trailing Interpolation Buffer Selection").set("Interpolation Buffer Type","Interpolation Buffer");
    ibPL->sublist("Interpolation Buffer Settings").sublist("Trailing Interpolation Buffer Selection").sublist("Interpolation Buffer").set("StorageLimit",storageLimit);
    ibPL->sublist("Interpolation Buffer Settings").sublist("Interpolator Selection").set("Interpolator Type","Linear Interpolator");
    ib->setParameterList(ibPL);
  }
  RCP<SinCosModel> fwdModel = sinCosModel(true); // implicit formulation
  Thyra::ModelEvaluatorBase::InArgs<double> fwdIC = fwdModel->getNominalValues();
  RCP<Thyra::NonlinearSolverBase<double> > fwdNLSolver = timeStepNonlinearSolver<double>();
  RCP<IntegratorBase<double> > fwdIntegrator = ib->create(fwdModel,fwdIC,fwdNLSolver);
  RCP<const VectorBase<double> > x_final;
  {
    Array<double> time_vec;
    time_vec.push_back(finalTime);
    Array<RCP<const Thyra::VectorBase<double> > > x_final_array;
    fwdIntegrator->getFwdPoints(time_vec,&x_final_array,NULL,NULL);
    x_final = x_final_array[0];
  }
  // Verify x_final is correct
  {
    // Defaults from SinCos Model:
    double f = 1.0;
    double L = 1.0;
    double a = 0.0;
    double x_ic_0 = 0.0;
    double x_ic_1 = 1.0;
    double x_0 = dt/(1.0+std::pow(dt*f/L,2))*(x_ic_0/dt+x_ic_1+dt*std::pow(f/L,2)*a);
    double x_1 = dt/(1.0+std::pow(dt*f/L,2))*(-std::pow(f/L,2)*x_ic_0+x_ic_1/dt+std::pow(f/L,2)*a);
    double tol = 1.0e-10;
    Thyra::ConstDetachedVectorView<double> x_final_view( *x_final );
    TEST_FLOATING_EQUALITY( x_final_view[0], x_0, tol );
    TEST_FLOATING_EQUALITY( x_final_view[1], x_1, tol );
  }
  // Copy InterpolationBuffer data into Cubic Spline interpolation buffer for use in Adjoint Solve
  TimeRange<double> fwdTimeRange; 
  RCP<InterpolationBufferBase<double> > fwdCubicSplineInterpBuffer;
  {
    RCP<PointwiseInterpolationBufferAppender<double> > piba = pointwiseInterpolationBufferAppender<double>();
    RCP<InterpolationBuffer<double> > sinkInterpBuffer = interpolationBuffer<double>();
    sinkInterpBuffer->setStorage(storageLimit);
    RCP<CubicSplineInterpolator<double> > csi = cubicSplineInterpolator<double>();
    sinkInterpBuffer->setInterpolator(csi);
    RCP<const InterpolationBufferBase<double> > sourceInterpBuffer;
    {
      RCP<TrailingInterpolationBufferAcceptingIntegratorBase<double> > tibaib = 
        Teuchos::rcp_dynamic_cast<TrailingInterpolationBufferAcceptingIntegratorBase<double> >(fwdIntegrator,true);
      sourceInterpBuffer = tibaib->getTrailingInterpolationBuffer();
    }
    fwdTimeRange = sourceInterpBuffer->getTimeRange();
    piba->append(*sourceInterpBuffer, fwdTimeRange, Teuchos::outArg(*sinkInterpBuffer));
    fwdCubicSplineInterpBuffer = sinkInterpBuffer;

    TimeRange<double> sourceRange = sourceInterpBuffer->getTimeRange();
    TimeRange<double> sinkRange = sinkInterpBuffer->getTimeRange();
    TEST_EQUALITY( sourceRange.lower(), sinkRange.lower() );
    TEST_EQUALITY( sourceRange.upper(), sinkRange.upper() );
  }
  // Adjoint Solve, reading forward solve data from Cubic Spline interpolation buffer
  {
    RCP<ParameterList> ibPL = Teuchos::parameterList();
    ibPL->sublist("Integrator Settings").sublist("Integrator Selection").set("Integrator Type","Default Integrator");
    ibPL->sublist("Integrator Settings").set("Final Time",finalTime);
    ibPL->sublist("Integration Control Strategy Selection").set("Integration Control Strategy Type","Simple Integration Control Strategy");
    ibPL->sublist("Integration Control Strategy Selection").sublist("Simple Integration Control Strategy").set("Take Variable Steps",false);
    ibPL->sublist("Integration Control Strategy Selection").sublist("Simple Integration Control Strategy").set("Fixed dt",dt);

    ibPL->sublist("Stepper Settings").sublist("Stepper Selection").set("Stepper Type","Backward Euler");
    //ibPL->sublist("Stepper Settings").sublist("Stepper Selection").set("Stepper Type","Implicit RK");
    //ibPL->sublist("Stepper Settings").sublist("Runge Kutta Butcher Tableau Selection").set("Runge Kutta Butcher Tableau Type","Implicit 1 Stage 2nd order Gauss");
    ibPL->sublist("Interpolation Buffer Settings").sublist("Trailing Interpolation Buffer Selection").set("Interpolation Buffer Type","Interpolation Buffer");
    ibPL->sublist("Interpolation Buffer Settings").sublist("Trailing Interpolation Buffer Selection").sublist("Interpolation Buffer").set("StorageLimit",storageLimit);
    ibPL->sublist("Interpolation Buffer Settings").sublist("Interpolator Selection").set("Interpolator Type","Linear Interpolator");
    ib->setParameterList(ibPL);
  }
  RCP<Thyra::ModelEvaluator<double> > adjModel;
  {
    RCP<Rythmos::AdjointModelEvaluator<double> > model = 
      Rythmos::adjointModelEvaluator<double>(
          fwdModel, fwdTimeRange
          );
    //model->setFwdStateSolutionBuffer(fwdCubicSplineInterpBuffer);
    adjModel = model;
  }
//.........这里部分代码省略.........