C++ InArgs::get_p方法代码示例

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;
}

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];

        }
    }

}

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

  int status=0; // 0 = pass, failures are incremented
  bool success = true;
  Teuchos::GlobalMPISession mpiSession(&argc,&argv);

#ifdef ENABLE_CHECK_FPE
   // Catch FPEs
   _mm_setcsr(_MM_MASK_MASK &~
		(_MM_MASK_OVERFLOW | _MM_MASK_INVALID | _MM_MASK_DIV_ZERO) );
#endif

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

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

  // Command-line argument for input file
  std::string xmlfilename;
  if(argc > 1){

    if(!strcmp(argv[1],"--help")){
      printf("albany [inputfile.xml]\n");
      exit(1);
    }
    else
      xmlfilename = argv[1];

  }
  else
    xmlfilename = "input.xml";

  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

    Albany::SolverFactory slvrfctry(xmlfilename, Albany_MPI_COMM_WORLD);
    RCP<Epetra_Comm> appComm = Albany::createEpetraCommFromMpiComm(Albany_MPI_COMM_WORLD);
    RCP<Albany::Application> app;
    const RCP<Thyra::ModelEvaluator<double> > solver =
      slvrfctry.createThyraSolverAndGetAlbanyApp(app, appComm, appComm);

    setupTimer.~TimeMonitor();

    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);

    *out << "After main solve" << std::endl;

    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;

    const Thyra::ModelEvaluatorBase::InArgs<double> nominal = solver->getNominalValues();
    for (int i=0; i<num_p; i++) {
      const Teuchos::RCP<const Epetra_Vector> p_init = epetraVectorFromThyra(appComm, nominal.get_p(i));
      p_init->Print(*out << "\nParameter vector " << i << ":\n");
    }

    for (int i=0; i<num_g-1; i++) {
      const RCP<const Epetra_Vector> g = responses[i];
      bool is_scalar = true;

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

      if (is_scalar) {
        g->Print(*out << "\nResponse vector " << i << ":\n");

        if (num_p == 0) {
          // Just calculate regression data
          status += slvrfctry.checkSolveTestResults(i, 0, g.get(), NULL);
        } else {
          for (int j=0; j<num_p; j++) {
            const RCP<const Epetra_MultiVector> dgdp = sensitivities[i][j];
            if (Teuchos::nonnull(dgdp)) {
              dgdp->Print(*out << "\nSensitivities (" << i << "," << j << "):!\n");
//.........这里部分代码省略.........

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.
//.........这里部分代码省略.........

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;
    }
  }
//.........这里部分代码省略.........

// 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
//.........这里部分代码省略.........

//.........这里部分代码省略.........
   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
        }

        if (num_p == 0 && cur_time_yr == final_time) {
          // Just calculate regression data -- only if in final time step
#ifdef CISM_USE_EPETRA
          status += slvrfctry->checkSolveTestResults(i, 0, g.get(), NULL);

// hide the original parental method AMET->evalModelImpl():
void
Aeras::HVDecorator::evalModelImpl(
    const Thyra::ModelEvaluatorBase::InArgs<ST>& inArgsT,
    const Thyra::ModelEvaluatorBase::OutArgs<ST>& outArgsT) const
{

  std::cout << "DEBUG WHICH HVDecorator: " << __PRETTY_FUNCTION__ << "\n";
	
  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;

#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;
//.........这里部分代码省略.........

  virtual
  void evalModelImpl(
      const Thyra::ModelEvaluatorBase::InArgs<double> &in_args,
      const Thyra::ModelEvaluatorBase::OutArgs<double> &out_args
      ) const
  {
    const auto & x_in = in_args.get_x();
#ifndef NDEBUG
    TEUCHOS_ASSERT(!x_in.is_null());
#endif
    // create corresponding tpetra vector
    auto x_in_tpetra =
      Thyra::TpetraOperatorVectorExtraction<double,int,int>::getConstTpetraVector(
          x_in
          );

    // compute F
    const auto & f_out = out_args.get_f();
    if (!f_out.is_null()) {
      // Dissect in_args.get_p(0) into parameter sublists.
      const auto & p_in = in_args.get_p(0);
#ifndef NDEBUG
      TEUCHOS_ASSERT(!p_in.is_null());
      // Make sure the parameters aren't NaNs.
      for (int k = 0; k < p_in->space()->dim(); k++) {
        TEUCHOS_ASSERT(!std::isnan(Thyra::get_ele(*p_in, k)));
      }
#endif
      // Fill the parameters into a std::map.
      const auto param_names = this->get_p_names(0);
      const double alph = Thyra::get_ele(*p_in, 0);

      auto f_out_tpetra =
        Thyra::TpetraOperatorVectorExtraction<double,int,int>::getTpetraVector(
            f_out
            );
      const auto x_data = x_in_tpetra->getData();
      auto f_data = f_out_tpetra->getDataNonConst();
      for (size_t i = 0; i < f_data.size(); i++) {
        f_data[i] = x_data[i] * x_data[i] - alph;
      }
    }

    // Compute df/dp.
    const auto & derivMv = out_args.get_DfDp(0).getDerivativeMultiVector();
    const auto & dfdp_out = derivMv.getMultiVector();
    if (!dfdp_out.is_null()) {
      auto dfdp_out_tpetra =
        Thyra::TpetraOperatorVectorExtraction<double,int,int>::getTpetraMultiVector(
            dfdp_out
            );

      TEUCHOS_ASSERT_EQUALITY(dfdp_out_tpetra->getNumVectors(), 1);
      auto out = dfdp_out_tpetra->getVectorNonConst(0);
      auto out_data = out->getDataNonConst();
      for (size_t k = 0; k < out_data.size(); k++) {
        out_data[k] = -1.0;
      }
    }

    // Fill Jacobian.
    const auto & W_out = out_args.get_W_op();
    if(!W_out.is_null()) {
      auto W_outT =
        Thyra::TpetraOperatorVectorExtraction<double,int,int>::getTpetraOperator(
            W_out
            );
      const auto & jac = Teuchos::rcp_dynamic_cast<jac_sqrt_alpha>(W_outT, true);
      jac->set_x0(*x_in_tpetra);
    }

    return;
  }

void Piro::RythmosSolver<Scalar>::evalModelImpl(
#endif
    const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
    const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs) const
{
    using Teuchos::RCP;
    using Teuchos::rcp;

    // TODO: Support more than 1 parameter and 1 response
    const int j = 0;
    const int l = 0;

    // Parse InArgs
    RCP<const Thyra::VectorBase<Scalar> > p_in;
    if (num_p > 0) {
        p_in = inArgs.get_p(l);
    }
    RCP<const Thyra::VectorBase<Scalar> > p_in2;  //JF add for multipoint
    if (num_p > 1) {
        p_in2 = inArgs.get_p(l+1);
    }

    // Parse OutArgs
    RCP<Thyra::VectorBase<Scalar> > g_out;
    if (num_g > 0) {
        g_out = outArgs.get_g(j);
    }
    const RCP<Thyra::VectorBase<Scalar> > gx_out = outArgs.get_g(num_g);

    Thyra::ModelEvaluatorBase::InArgs<Scalar> state_ic = model->getNominalValues();

    // Set initial time in ME if needed

    if(t_initial > 0.0 && state_ic.supports(Thyra::ModelEvaluatorBase::IN_ARG_t))

        state_ic.set_t(t_initial);

    if (Teuchos::nonnull(initialConditionModel)) {
        // The initial condition depends on the parameter
        // It is found by querying the auxiliary model evaluator as the last response
        const RCP<Thyra::VectorBase<Scalar> > initialState =
            Thyra::createMember(model->get_x_space());

        {
            Thyra::ModelEvaluatorBase::InArgs<Scalar> initCondInArgs = initialConditionModel->createInArgs();
            if (num_p > 0) {
                initCondInArgs.set_p(l, inArgs.get_p(l));
            }

            Thyra::ModelEvaluatorBase::OutArgs<Scalar> initCondOutArgs = initialConditionModel->createOutArgs();
            initCondOutArgs.set_g(initCondOutArgs.Ng() - 1, initialState);

            initialConditionModel->evalModel(initCondInArgs, initCondOutArgs);
        }

        state_ic.set_x(initialState);
    }

    // Set paramters p_in as part of initial conditions
    if (num_p > 0) {
        if (Teuchos::nonnull(p_in)) {
            state_ic.set_p(l, p_in);
        }
    }
    if (num_p > 1) { //JF added for multipoint
        if (Teuchos::nonnull(p_in2)) {
            state_ic.set_p(l+1, p_in2);
        }
    }

    *out << "\nstate_ic:\n" << Teuchos::describe(state_ic, solnVerbLevel);

    //JF  may need a version of the following for multipoint, i.e. num_p>1, l+1, if we want sensitivities
    RCP<Thyra::MultiVectorBase<Scalar> > dgxdp_out;
    Thyra::ModelEvaluatorBase::Derivative<Scalar> dgdp_deriv_out;
    if (num_p > 0) {
        const Thyra::ModelEvaluatorBase::DerivativeSupport dgxdp_support =
            outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, num_g, l);
        if (dgxdp_support.supports(Thyra::ModelEvaluatorBase::DERIV_MV_JACOBIAN_FORM)) {
            const Thyra::ModelEvaluatorBase::Derivative<Scalar> dgxdp_deriv =
                outArgs.get_DgDp(num_g, l);
            dgxdp_out = dgxdp_deriv.getMultiVector();
        }

        if (num_g > 0) {
            const Thyra::ModelEvaluatorBase::DerivativeSupport dgdp_support =
                outArgs.supports(Thyra::ModelEvaluatorBase::OUT_ARG_DgDp, j, l);
            if (!dgdp_support.none()) {
                dgdp_deriv_out = outArgs.get_DgDp(j, l);
            }
        }
    }

    const bool requestedSensitivities =
        Teuchos::nonnull(dgxdp_out) || !dgdp_deriv_out.isEmpty();

    RCP<const Thyra::VectorBase<Scalar> > finalSolution;
    if (!requestedSensitivities) {
        //
        *out << "\nE) Solve the forward problem ...\n";
//.........这里部分代码省略.........

void
Piro::VelocityVerletSolver<Scalar, LocalOrdinal, GlobalOrdinal, Node>::evalModelImpl(
    const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
    const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs) const
{
  using Teuchos::RCP;
  using Teuchos::rcp;

  // TODO: Support more than 1 parameter and 1 response
  const int j = 0;
  const int l = 0;

  // Parse InArgs
  RCP<const Thyra::VectorBase<Scalar> > p_in;
  if (num_p > 0) {
    p_in = inArgs.get_p(l);
  }

  // Parse OutArgs
  RCP<Thyra::VectorBase<Scalar> > g_out;
  if (num_g > 0) {
    g_out = outArgs.get_g(j);
  }
  const RCP<Thyra::VectorBase<Scalar> > gx_out = outArgs.get_g(num_g);

// create a new vector and fill it with the contents of model->get_x()
  // Build a multivector holding x (0th vector), v (1st vector), and a (2nd vector)
//  Teuchos::RCP<Thyra::MultiVectorBase<Scalar> > soln = createMembers(model->get_x_space(), 3);

// create a new vector and fill it with the contents of model->get_x()
/*
  Teuchos::RCP<Thyra::VectorBase<Scalar> > x = soln->col(0);
  assign(x.ptr(), *model->getNominalValues().get_x());
  Teuchos::RCP<Thyra::VectorBase<Scalar> > v = soln->col(1);
  assign(v.ptr(), *model->getNominalValues().get_x_dot());
  Teuchos::RCP<Thyra::VectorBase<Scalar> > a = soln->col(2);
  assign(a.ptr(), *model->get_x_dotdot());
*/

  Teuchos::RCP<Thyra::VectorBase<Scalar> > x = model->getNominalValues().get_x()->clone_v();
  Teuchos::RCP<Thyra::VectorBase<Scalar> > v = model->getNominalValues().get_x_dot()->clone_v();

// Note that Thyra doesn't have x_dotdot - go get it from the transient decorator around the Albany model
//  Teuchos::RCP<Thyra::VectorBase<Scalar> > a = model->get_x_dotdot()->clone_v();

  Teuchos::RCP<Thyra::DefaultModelEvaluatorWithSolveFactory<Scalar> >
      DMEWSF(Teuchos::rcp_dynamic_cast<Thyra::DefaultModelEvaluatorWithSolveFactory<Scalar> >(model));

  Teuchos::RCP<const Piro::TransientDecorator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > dec =
       Teuchos::rcp_dynamic_cast<const Piro::TransientDecorator<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
           (DMEWSF->getUnderlyingModel());

  TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(dec), std::logic_error,
      "Underlying model in VelovityVerletSolver does not cast to a Piro::TransientDecorator<Scalar, LocalOrdinal, GlobalOrdinal, Node>"
      << std::endl);

  Teuchos::RCP<Thyra::VectorBase<Scalar> > a = dec->get_x_dotdot()->clone_v();

  RCP<Thyra::VectorBase<Scalar> > finalSolution;

  // Zero out the acceleration vector
  put_scalar(0.0, a.ptr());

  TEUCHOS_TEST_FOR_EXCEPTION(v == Teuchos::null || x == Teuchos::null,
                     Teuchos::Exceptions::InvalidParameter,
                     std::endl << "Error in Piro::VelocityVerletSolver " <<
                     "Requires initial x and x_dot: " << std::endl);

  Scalar t = t_init;

  // Observe initial condition
//  if (observer != Teuchos::null) observer->observeSolution(*soln, t);

  Scalar vo = norm_2(*v);
  *out << "Initial Velocity = " << vo << std::endl;

   if (Teuchos::VERB_MEDIUM <= solnVerbLevel) *out << std::endl;

   Thyra::ModelEvaluatorBase::InArgs<Scalar> model_inargs = model->createInArgs();
   Thyra::ModelEvaluatorBase::OutArgs<Scalar> model_outargs = model->createOutArgs();
   model_inargs.set_x(x);
   if (num_p > 0)  model_inargs.set_p(0, p_in);

   model_outargs.set_f(a);
   if (g_out != Teuchos::null) model_outargs.set_g(0, g_out);

   Scalar ddt = 0.5 * delta_t * delta_t;

   // Calculate acceleration at time 0
   model->evalModel(model_inargs, model_outargs);

   for (int timeStep = 1; timeStep <= numTimeSteps; timeStep++) {

//     x->Update(delta_t, *v, ddt, *a, 1.0);
     Vp_StV(x.ptr(), delta_t, *v);
     Vp_StV(x.ptr(), ddt, *a);

     t += delta_t;
     model_inargs.set_t(t);

//.........这里部分代码省略.........

void
Piro::LOCAAdaptiveSolver<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 LOCAAdaptive 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];
    }

//    solMgr_->getSolutionGroup()->setParams(paramVector_);
    Teuchos::rcp_dynamic_cast< ::Thyra::LOCAAdaptiveState >(solMgr_->getState())
                 ->getSolutionGroup()->setParams(paramVector_);
  }

  LOCA::Abstract::Iterator::IteratorStatus status;

  status = stepper_->run();

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

  // The time spent
  globalData_->locaUtils->out() << std::endl <<
    "#### Statistics ########" << std::endl;

  // Check number of steps
  int numSteps = stepper_->getStepNumber();
  globalData_->locaUtils->out() << std::endl <<
    " Number of continuation Steps = " << numSteps << std::endl;

  // Check number of failed steps
  int numFailedSteps = stepper_->getNumFailedSteps();
  globalData_->locaUtils->out() << std::endl <<
    " Number of failed continuation Steps = " << numFailedSteps << std::endl;

  globalData_->locaUtils->out() << std::endl;


  // Note: the last g is used to store the final solution. It can be null - if it is just
  // skip the store. If adaptation has occurred, g is not the correct size.

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

  int x_args_dim = 0;
  int f_sol_dim = 0;

  // Pardon the nasty cast to resize the last g in outArgs - need to fit the solution
  Thyra::ModelEvaluatorBase::OutArgs<Scalar>* mutable_outArgsPtr =
    const_cast<Thyra::ModelEvaluatorBase::OutArgs<Scalar>* >(&outArgs);

  if(Teuchos::nonnull(x_outargs)){ // g has been allocated, calculate the sizes of g and the solution

    x_args_dim = x_outargs->space()->dim();
//    f_sol_dim = solMgr_->getSolutionGroup()->getX().length();
    f_sol_dim = Teuchos::rcp_dynamic_cast< ::Thyra::LOCAAdaptiveState >(solMgr_->getState())
          ->getSolutionGroup()->getX().length();


  }

  if(Teuchos::is_null(x_outargs) || (x_args_dim != f_sol_dim)){ // g is not the right size

      x_final = Thyra::createMember(this->get_g_space(this->num_g()));

      mutable_outArgsPtr->set_g(this->num_g(), x_final);

  }
  else { // g is OK, use it
    x_final = x_outargs;
  }

  {
    // Deep copy final solution from LOCA group
    NOX::Thyra::Vector finalSolution(x_final);
//    finalSolution = solMgr_->getSolutionGroup()->getX();
    finalSolution = Teuchos::rcp_dynamic_cast< ::Thyra::LOCAAdaptiveState >(solMgr_->getState())
                      ->getSolutionGroup()->getX();

  }

  // If the arrays need resizing
  if(x_args_dim != f_sol_dim){
//.........这里部分代码省略.........

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);
//.........这里部分代码省略.........

void Piro::NOXSolver<Scalar>::evalModelImpl(
    const Thyra::ModelEvaluatorBase::InArgs<Scalar>& inArgs,
    const Thyra::ModelEvaluatorBase::OutArgs<Scalar>& outArgs ) const
{
    using Teuchos::RCP;
    using Teuchos::rcp;

    //*out << "In eval Modle " << endl;

    // Parse InArgs

    RCP<const Thyra::VectorBase<Scalar> > p_in;
    if (num_p > 0) p_in = inArgs.get_p(0);

    // Parse OutArgs: always 1 extra
    RCP< Thyra::VectorBase<Scalar> > g_out;
    if (num_g > 0) g_out = outArgs.get_g(0);
    RCP< Thyra::VectorBase<Scalar> > gx_out = outArgs.get_g(num_g);

    // Parse out-args for sensitivity calculation


    ::Thyra::SolveCriteria<double> solve_criteria;
    ::Thyra::SolveStatus<double> solve_status;

    RCP< ::Thyra::VectorBase<double> >
    initial_guess = model->getNominalValues().get_x()->clone_v();

    solve_status = solver->solve(initial_guess.get(), &solve_criteria, NULL);

//   if (solve_status.solveStatus == ::Thyra::SOLVE_STATUS_CONVERGED)
//     std::cout << "Test passed!" << std::endl;

    // return the final solution as an additional g-vector, if requested
    RCP<const Thyra::VectorBase<Scalar> > finalSolution = solver->get_current_x();

    if (gx_out != Teuchos::null)  Thyra::copy(*finalSolution, gx_out.ptr());

    if (g_out != Teuchos::null) {
        // As post-processing step, calc responses at final solution
        Thyra::ModelEvaluatorBase::InArgs<Scalar>  model_inargs = model->createInArgs();
        Thyra::ModelEvaluatorBase::OutArgs<Scalar> model_outargs = model->createOutArgs();
        model_inargs.set_x(finalSolution);
        if (num_p > 0)  model_inargs.set_p(0, p_in);
        if (g_out != Teuchos::null) {
            Thyra::put_scalar(0.0,g_out.ptr());
            model_outargs.set_g(0, g_out);
        }

        model->evalModel(model_inargs, model_outargs);
    }

    /*********************  NEED TO CONVERT TO THYRA *******************
      RCP< Thyra::MultiVectorBase<Scalar> > dgdp_out;
      if (num_p>0 && num_g>0)
        dgdp_out = outArgs.get_DgDp(0,0).getMultiVector();

      if (dgdp_out == Teuchos::null) {

         Teuchos::RCP<Epetra_MultiVector> dgdx
              = Teuchos::rcp(new Epetra_MultiVector(finalSolution->Map(),
                                                       dgdp_out->GlobalLength()));
         Teuchos::Array<int> p_indexes =
           outArgs.get_DgDp(0,0).getDerivativeMultiVector().getParamIndexes();

         EpetraExt::ModelEvaluator::DerivativeMultiVector dmv_dgdp(dgdp_out,
                                                                   DERIV_MV_BY_COL,
                                                                   p_indexes);

         EpetraExt::ModelEvaluator::InArgs model_inargs = model->createInArgs();
         EpetraExt::ModelEvaluator::OutArgs model_outargs = model->createOutArgs();
         model_inargs.set_x(finalSolution);
         model_inargs.set_p(0, p_in);

         if (g_out != Teuchos::null) {
           g_out->PutScalar(0.0);
           model_outargs.set_g(0, g_out);
         }
         model_outargs.set_DgDp(0,0,dmv_dgdp);
         model_outargs.set_DgDx(0,dgdx);

         model->evalModel(model_inargs, model_outargs);


         // (3) Calculate dg/dp = dg/dx*dx/dp + dg/dp
         // This may be the transpose of what we want since we specified
         // we want dg/dp by column in createOutArgs().
         // In this case just interchange the order of dgdx and dxdp
         // We should really probably check what the underlying ME does

         if (Teuchos::VERB_MEDIUM <= solnVerbLevel) cout << " dgdx \n" << *dgdx << endl;
         if (Teuchos::VERB_MEDIUM <= solnVerbLevel) cout << " dxdp \n" << *dxdp << endl;

         dgdp_out->Multiply('T', 'N', 1.0, *dgdx, *dxdp, 1.0);
       }
    *********************/
}