C++ InArgs类代码示例

void 
Stokhos::MPInverseModelEvaluator::evalModel(const InArgs& inArgs,
					    const OutArgs& outArgs) const
{
  // Create underlying inargs
  InArgs me_inargs = me->createInArgs();

  // Pass parameters
  for (int i=0; i<num_p; i++)
    me_inargs.set_p(i, inArgs.get_p(i));

  // Pass MP parameters
  for (int i=0; i<num_p_mp; i++) {
    mp_const_vector_t p_mp = inArgs.get_p_mp(mp_p_index_map[i]);
    if (p_mp != Teuchos::null) {
      me_inargs.set_p(i+num_p, p_mp->getBlockVector());
    }
  }

  // Create underlying outargs
  OutArgs me_outargs = me->createOutArgs();


  // MP Responses
  for (int i=0; i<num_g_mp; i++) {
    int ii = mp_g_index_map[i];

    // g
    mp_vector_t g_mp = outArgs.get_g_mp(ii);
    if (g_mp != Teuchos::null) {
      me_outargs.set_g(i, Teuchos::rcp_dynamic_cast<Epetra_Vector>(g_mp->getBlockVector()));
    }

    // dg/dp
    for (int j=0; j<num_p; j++) {
      if (!outArgs.supports(OUT_ARG_DgDp_mp, ii, j).none()) {
	MPDerivative dgdp_mp = outArgs.get_DgDp_mp(ii,j);
	Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdp_mp_mv = 
	  dgdp_mp.getMultiVector();
	Teuchos::RCP<Epetra_Operator> dgdp_mp_op =
	  dgdp_mp.getLinearOp();
	if (dgdp_mp_mv != Teuchos::null) {
	  me_outargs.set_DgDp(
	    i, j, Derivative(dgdp_mp_mv->getBlockMultiVector(),
			     dgdp_mp.getMultiVectorOrientation()));
	}
	else if (dgdp_mp_op != Teuchos::null) {
	  me_outargs.set_DgDp(i, j, Derivative(dgdp_mp_op));
	}
      }
    }

  }

  // Compute the functions
  me->evalModel(me_inargs, me_outargs);

}

Path::Path(const std::string &path, const PathArgument &a1, const PathArgument &a2, const PathArgument &a3, const PathArgument &a4, const PathArgument &a5) {
    InArgs in;
    in.push_back( &a1 );
    in.push_back( &a2 );
    in.push_back( &a3 );
    in.push_back( &a4 );
    in.push_back( &a5 );
    makePath( path, in );
}

void Piro::Epetra::MatrixFreeDecorator::evalModel( const InArgs& inArgs,
                                     const OutArgs& outArgs ) const
{
  using Teuchos::RCP;
  using Teuchos::rcp;

  RCP<Epetra_Operator> W_out = outArgs.get_W();

  if (W_out == Teuchos::null) {
    // Just pass through as is: nothing to Decorate
    model->evalModel(inArgs, outArgs);
  }
  else {

    RCP<Piro::Epetra::MatrixFreeOperator> W_mfo =
      Teuchos::rcp_dynamic_cast<Piro::Epetra::MatrixFreeOperator>(W_out);

    TEUCHOS_TEST_FOR_EXCEPTION(W_mfo==Teuchos::null, std::logic_error, 
      "Epetra_Operator sent as W to Piro::Epetra::MatrixFreeDecorator\n" 
      "be of type Piro::Epetra::MatrixFreeOperator");
   
    // Do base case for MatrixFree: set f instead of W
    OutArgs modelOutArgs(outArgs);
    InArgs modelInArgs(inArgs);

    // Store f_out in case it was also requested
    RCP<Epetra_Vector> f_out = outArgs.get_f();

    modelOutArgs.set_f(fBase);
    modelOutArgs.set_W(Teuchos::null);

    //Evaluate the underlying model
    model->evalModel(modelInArgs, modelOutArgs);

    // If f_out was requested, return it.
    if (f_out != Teuchos::null) *f_out = *fBase;

    // Save unperturbed solution (deep copy inArgs, shallow f)
    InArgs clonedInArgs = inArgs;
    for (int l = 0; l < inArgs.Np(); ++l) {
      const RCP<const Epetra_Vector> p_l = inArgs.get_p(l);
      if (nonnull(p_l))
        clonedInArgs.set_p(l, Teuchos::rcp(new Epetra_Vector(*p_l)));
    }
    clonedInArgs.set_x(Teuchos::rcp(new Epetra_Vector(*inArgs.get_x())));

    bool haveXdot = false;
    if (inArgs.supports(IN_ARG_x_dot)) {
      RCP<const Epetra_Vector> xdot = inArgs.get_x_dot();
      if (nonnull(xdot)) {
        clonedInArgs.set_x_dot(Teuchos::rcp(new Epetra_Vector(*inArgs.get_x_dot())));
        haveXdot = true;
      }
    }
    W_mfo->setBase(clonedInArgs, fBase, haveXdot);
  }
}

  TEUCHOS_UNIT_TEST(thyra_model_evaluator, basic)
  {
    using Teuchos::RCP;

    bool parameter_on = true;
    Teuchos::RCP<panzer::FieldManagerBuilder> fmb;  
    Teuchos::RCP<panzer::ResponseLibrary<panzer::Traits> > rLibrary; 
    Teuchos::RCP<panzer::LinearObjFactory<panzer::Traits> > lof;
    Teuchos::RCP<panzer::GlobalData> gd;
  
    buildAssemblyPieces(parameter_on,fmb,rLibrary,gd,lof);

    // Test a transient me
    {
      typedef Thyra::ModelEvaluatorBase MEB;
      typedef Thyra::ModelEvaluatorBase::InArgs<double> InArgs;
      typedef Thyra::ModelEvaluatorBase::OutArgs<double> OutArgs;
      typedef Thyra::VectorBase<double> VectorType;
      typedef Thyra::LinearOpBase<double> OperatorType;
      typedef panzer::ModelEvaluator<double> PME;

      std::vector<Teuchos::RCP<Teuchos::Array<std::string> > > p_names;
      bool build_transient_support = true;
    
      RCP<PME> me = Teuchos::rcp(new PME(fmb,rLibrary,lof,p_names,Teuchos::null,gd,build_transient_support,0.0));

      InArgs in_args = me->createInArgs();
      OutArgs out_args = me->createOutArgs();
      
      TEST_ASSERT(in_args.supports(MEB::IN_ARG_x));
      TEST_ASSERT(in_args.supports(MEB::IN_ARG_x_dot));
      TEST_ASSERT(in_args.supports(MEB::IN_ARG_alpha));
      TEST_ASSERT(in_args.supports(MEB::IN_ARG_beta));
      TEST_ASSERT(out_args.supports(MEB::OUT_ARG_f));
      TEST_ASSERT(out_args.supports(MEB::OUT_ARG_W_op));

      InArgs nomValues = me->getNominalValues();
      RCP<const VectorType> x = nomValues.get_x();
      RCP<VectorType> x_dot = Thyra::createMember(*me->get_x_space());
      Thyra::assign(x_dot.ptr(),0.0);
      in_args.set_x(x);
      in_args.set_x_dot(x_dot);
      in_args.set_alpha(0.0);
      in_args.set_beta(1.0);
      
      RCP<VectorType> f = Thyra::createMember(*me->get_f_space());
      RCP<OperatorType> J_tmp = me->create_W_op();
      out_args.set_f(f);
      out_args.set_W_op(J_tmp);

      me->evalModel(in_args, out_args);
    }
  }

Path::Path(const JSONCPP_STRING& path,
           const PathArgument& a1,
           const PathArgument& a2,
           const PathArgument& a3,
           const PathArgument& a4,
           const PathArgument& a5) {
  InArgs in;
  in.push_back(&a1);
  in.push_back(&a2);
  in.push_back(&a3);
  in.push_back(&a4);
  in.push_back(&a5);
  makePath(path, in);
}

void Simple_ModelEval::evalModel( const InArgs& inArgs,
                              const OutArgs& outArgs ) const
{

  // Parse InArgs
  Teuchos::RCP<const Epetra_Vector> p_in = inArgs.get_p(0);
  if (!p_in.get()) cout << "ERROR: Simple_ModelEval requires p as inargs" << endl;
  int numParameters = p_in->GlobalLength();

  // Parse OutArgs

  Teuchos::RCP<Epetra_Vector> g_out = outArgs.get_g(0); 

  // Parse out-args for sensitivity calculation
  Teuchos::RCP<Epetra_MultiVector> dgdp_out;

  dgdp_out = outArgs.get_DgDp(0,0).getMultiVector();

  if (!is_null(g_out)) {
    (*g_out)[0] = 1.0 - (*p_in)[0];
    (*g_out)[1] = 1.2 - (*p_in)[1];
    (*g_out)[2] = 4.0 - (*p_in)[2] - 0.5* (1.0 - (*p_in)[0]);
  }

  if (dgdp_out != Teuchos::null) {
     // Must initialize since Thyra will init with NaN
     dgdp_out->PutScalar(0.0);
     // Set gradient of above g equations (derived by hand)
     for (int i=0; i<numParameters; i++) {
       (*dgdp_out)[i][i] = -1.0;
     }
     (*dgdp_out)[0][2] = 0.5; //DERIV_BY_COL: [p][g]
   }
} 

 void Path::makePath(const std::string &path, const InArgs &in) {
   const char *current = path.c_str();
   const char *end = current + path.length();
   InArgs::const_iterator itInArg = in.begin();
   while (current != end) {
     if (*current == '[') {
       ++current;
       if (*current == '%')
         addPathInArg(path, in, itInArg, PathArgument::kindIndex);
       else {
         Value::UInt index = 0;
         for (; current != end && *current >= '0' && *current <= '9'; ++current)
           index = index * 10 + Value::UInt(*current - '0');
         args_.push_back(index);
       }
       if (current == end || *current++ != ']')
         invalidPath(path, int(current - path.c_str()));
     } else if (*current == '%') {
       addPathInArg(path, in, itInArg, PathArgument::kindKey);
       ++current;
     } else if (*current == '.') {
       ++current;
     } else {
       const char *beginName = current;
       while (current != end && !strchr("[.", *current))
         ++current;
       args_.push_back(std::string(beginName, current));
     }
   }
 }

// Evaluate model on InArgs
void trilinosModelEvaluator::evalModel(const InArgs& inArgs, const OutArgs& outArgs) const{

  //double nrm;
  // Get the solution vector x from inArgs and residual vector from outArgs
  RCP<const Epetra_Vector> x = inArgs.get_x();
  RCP<Epetra_Vector> f = outArgs.get_f();
  
  if (x == Teuchos::null) throw "trilinosModelEvaluator::evalModel: x was NOT specified!";

  // Save the current solution, which makes it initial guess for next nonlienar solve
  *xVec = *x;

  if (f != Teuchos::null) {
    f->PutScalar(0.0);
    calc_F(x->Values(), f->Values(), N, blackbox_res);

    //f->Norm2(&nrm);
    //cout << "AGS  Resid norm in eval_model total " << nrm << endl;
  }

  RCP<Epetra_Operator> WPrec = outArgs.get_WPrec();
  if (WPrec != Teuchos::null) {
     //cout << "evalModel called for WPrec -- doing nothing " <<  endl;
  }
}

 void Path::addPathInArg(const std::string &path,
                         const InArgs &in,
                         InArgs::const_iterator &itInArg,
                         PathArgument::Kind kind) {
   if (itInArg == in.end()) {
     // Error: missing argument %d
   } else if ((*itInArg)->kind_ != kind) {
     // Error: bad argument type
   } else {
     args_.push_back(**itInArg);
   }
 }

void Path::addPathInArg(const JSONCPP_STRING& /*path*/,
                        const InArgs& in,
                        InArgs::const_iterator& itInArg,
                        PathArgument::Kind kind) {
  if (itInArg == in.end()) {
    // Error: missing argument %d
  } else if ((*itInArg)->kind_ != kind) {
    // Error: bad argument type
  } else {
    args_.push_back(**itInArg++);
  }
}

void ExampleApplication1Dfem::evalModel( const InArgs& inArgs, const OutArgs& outArgs ) const
{
  Teuchos::RCP<const Epetra_Vector> x = inArgs.get_x();
  Teuchos::RCP<const Epetra_Vector> xdot = inArgs.get_x_dot();
#ifdef EXAMPLEAPPLICATION_DEBUG
  std::cout << "ExampleApplication1Dfem::evalModel ---------------------------{" << std::endl;
  std::cout << "x = " << std::endl;
  x->Print(std::cout);
  std::cout << "xdot = " << std::endl;
  xdot->Print(std::cout);
#endif // EXAMPLEAPPLICATION_DEBUG
  Teuchos::RCP<Epetra_Vector> f;
  if( (f = outArgs.get_f()).get() ) 
  {
    NOX::Epetra::Interface::Required::FillType flag = NOX::Epetra::Interface::Required::Residual;
    problemInterfacePtr_->evaluate(flag,&*x,&*xdot,0.0,0.0,&*f,NULL);
#ifdef EXAMPLEAPPLICATION_DEBUG
    std::cout << "f = " << std::endl;
    f->Print(std::cout);
#endif // EXAMPLEAPPLICATION_DEBUG
  }
  Teuchos::RCP<Epetra_Operator> W;
  if( (W = outArgs.get_W()).get() ) 
  {
    const double alpha = inArgs.get_alpha();
    const double beta = inArgs.get_beta();
    NOX::Epetra::Interface::Required::FillType flag = NOX::Epetra::Interface::Required::Jac;
//    Epetra_CrsMatrix& jacobian = problemInterfacePtr_->getJacobian();
    Epetra_CrsMatrix& jac = Teuchos::dyn_cast<Epetra_CrsMatrix>(*W);
    problemInterfacePtr_->evaluate(flag,&*x,&*xdot,alpha,beta,NULL,&jac);
#ifdef EXAMPLEAPPLICATION_DEBUG
    std::cout << "jac = " << std::endl;
    jac.Print(std::cout);
#endif // EXAMPLEAPPLICATION_DEBUG
  }
#ifdef EXAMPLEAPPLICATION_DEBUG
  std::cout << "ExampleApplication1Dfem::evalModel ---------------------------}" << std::endl;
#endif // EXAMPLEAPPLICATION_DEBUG
}

void LorenzModel::evalModel( const InArgs& inArgs, const OutArgs& outArgs ) const
{
  const Epetra_Vector &yin = *(inArgs.get_x());
  const double t = inArgs.get_t(); // ignored
#ifdef LORENZMODEL_DEBUG
      std::cout << "----------------------------------------------------------------------" << std::endl;
      std::cout << "LorenzModel::evalModel yin = " << std::endl;
      yin.Print(std::cout);
#endif
  Epetra_Vector &yout = *outArgs.get_f();

  yout[0] = -param0 * yin[0] + param0 * yin[1];
  yout[1] = param1 * yin[0] - yin[1] - yin[0]*yin[2];
  yout[2] = -param2*yin[2] + yin[0]*yin[1];

#ifdef LORENZMODEL_DEBUG
  std::cout << "LorenzModel::evalModel (explicit) f = " << std::endl;
  yout.Print(std::cout);
#endif
#ifdef LORENZMODEL_DEBUG
  std::cout << "----------------------------------------------------------------------" << std::endl;
#endif
}

void Piro::Epetra::TrapezoidDecorator::evalModel( const InArgs& inArgs,
                                     const OutArgs& outArgs ) const
{
  using Teuchos::RCP;
  using Teuchos::rcp;

  // Copy outArgs; add time term 
  OutArgs modelOutArgs(outArgs);
  InArgs modelInArgs(inArgs);

  xDotDot->Update(fdt2, *inArgs.get_x(), -fdt2, *x_pred, 0.0);
  // WARNING:: Model must use  x_dot  as xDotDot!!!
  modelInArgs.set_x_dot(xDotDot); 

  modelInArgs.set_alpha(fdt2);  // fdt2 = 4/(dt)^2
  modelInArgs.set_beta(1.0);
  modelInArgs.set_t(time);

  //Evaluate the underlying model
  model->evalModel(modelInArgs, modelOutArgs);
}

void MockModelEval_A::evalModel( const InArgs& inArgs,
                              const OutArgs& outArgs ) const
{

  // Parse InArgs
  RCP<const Epetra_Vector> p_in = inArgs.get_p(0);
  if (!p_in.get()) cout << "ERROR: MockModelEval_A requires p as inargs" << endl;
  //int numParameters = p_in->GlobalLength();

  RCP<const Epetra_Vector> x_in = inArgs.get_x();
  if (!x_in.get()) cout << "ERROR: MockModelEval_A requires x as inargs" << endl;
  int vecLength = x_in->GlobalLength();
  int myVecLength = x_in->MyLength();

  // Parse OutArgs

  RCP<Epetra_Vector> f_out = outArgs.get_f(); 
  RCP<Epetra_Vector> g_out = outArgs.get_g(0); 
  Teuchos::RCP<Epetra_Operator> W_out = outArgs.get_W();
  Teuchos::RCP<Epetra_MultiVector> dfdp_out;
  if (outArgs.Np() > 0)
    dfdp_out = outArgs.get_DfDp(0).getMultiVector();
  RCP<Epetra_MultiVector> dgdp_out;
  dgdp_out = outArgs.get_DgDp(0,0).getMultiVector();
  RCP<Epetra_MultiVector> dgdx_out;
  dgdx_out = outArgs.get_DgDx(0).getMultiVector();

  if (f_out != Teuchos::null) {
    for (int i=0; i<myVecLength; i++) {
      int gid = x_in->Map().GID(i);
      if (gid==0) // x_0^2 = p_0
       (*f_out)[i] = (*x_in)[i] * (*x_in)[i] -  (*p_in)[i];
      else // x^2 = (i+p_1)^2
       (*f_out)[i] = (*x_in)[i] * (*x_in)[i] - (gid + (*p_in)[1])*(gid + (*p_in)[1]);
    }
  }
  if (W_out != Teuchos::null) {
    Teuchos::RCP<Epetra_CrsMatrix> W_out_crs =
      Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_out, true);
    W_out_crs->PutScalar(0.0);

    double diag=0.0;
    for (int i=0; i<myVecLength; i++) {
      diag = 2.0 * (*x_in)[i];
      W_out_crs->ReplaceMyValues(i, 1, &diag, &i);
    }
  }

  if (dfdp_out != Teuchos::null) {
    dfdp_out->PutScalar(0.0);
    for (int i=0; i<myVecLength; i++) {
      int gid = x_in->Map().GID(i);
      if   (gid==0) (*dfdp_out)[0][i] = -1.0;
      else          (*dfdp_out)[1][i] =  -2.0* (gid + (*p_in)[1]);
    }
  }

  // ObjFn = 0.5*(Sum(x)-Sum(p)-12)^2 + 0.5*(p0-1)^2:  min at 1,3

  double term1, term2;
  x_in->MeanValue(&term1); 
  term1 =  vecLength * term1 - ((*p_in)[0] + (*p_in)[1]) - 12.0;
  term2 = (*p_in)[0] - 1.0;
  
  if (!is_null(g_out)) {
    (*g_out)[0] = 0.5*term1*term1 + 0.5*term2*term2;
  }

  if (dgdx_out != Teuchos::null) {
     dgdx_out->PutScalar(term1);
   }
  if (dgdp_out != Teuchos::null) {
     dgdp_out->PutScalar(0.0);
     (*dgdp_out)[0][0] = -term1 + term2;
     (*dgdp_out)[0][1] = -term1;
   }

  // Modify for time dependent (implicit timeintegration or eigensolves
  // Check if time dependent
  RCP<const Epetra_Vector> x_dot = inArgs.get_x_dot();

  if (x_dot.get()) {
    double alpha =  inArgs.get_alpha();
    double beta  =  inArgs.get_beta();
    if (alpha==0.0 && beta==0.0) {
      cout << "MockModelEval Warning: alpha=beta=0 -- setting beta=1" << endl;
      beta = 1.0;
    }

    if (f_out != Teuchos::null) {
      for (int i=0; i<myVecLength; i++) {
         (*f_out)[i] = -alpha*(*x_dot)[i] + beta * (*f_out)[i];
      }
    }
    if (dfdp_out != Teuchos::null) {
      dfdp_out->Scale(beta);
    }
    if (W_out != Teuchos::null) {
      Teuchos::RCP<Epetra_CrsMatrix> W_out_crs =
        Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_out, true);
//.........这里部分代码省略.........

void VanDerPolOscillator::evalModel( const InArgs& inArgs, 
				     const OutArgs& outArgs ) const
{
  // compute f(x)
  Teuchos::RCP<const Epetra_Vector> x = inArgs.get_x();
  Teuchos::RCP<Epetra_Vector> f = outArgs.get_f();
  if ( (x != Teuchos::null) && (f != Teuchos::null) ) {
    evalVField((*x)[0],(*x)[1],(*f)[0],(*f)[1]);
  }

  // compute f([x])
  Teuchos::RCP<const Teuchos::Polynomial<Epetra_Vector> > x_poly = 
    inArgs.get_x_poly();
  Teuchos::RCP<Teuchos::Polynomial<Epetra_Vector> > f_poly = 
    outArgs.get_f_poly();
  if ( (x_poly != Teuchos::null) && (f_poly != Teuchos::null) ) {
    unsigned int d = x_poly->degree();
    Sacado::Tay::Taylor<double> x1(d,0.0);
    Sacado::Tay::Taylor<double> x2(d,0.0);
    Sacado::Tay::Taylor<double> f1(d,0.0);
    Sacado::Tay::Taylor<double> f2(d,0.0);

    for (unsigned int i=0; i<=d; i++) {
      x1.fastAccessCoeff(i) = (*(x_poly->getCoefficient(i)))[0];
      x2.fastAccessCoeff(i) = (*(x_poly->getCoefficient(i)))[1];
    }

    evalVField(x1,x2,f1,f2);

    for (unsigned int i=0; i<=d; i++) {
      (*(f_poly->getCoefficient(i)))[0] = f1.coeff(i);
      (*(f_poly->getCoefficient(i)))[1] = f2.coeff(i);
    }
  }

  // compute W
  Teuchos::RCP<Epetra_Operator> W = outArgs.get_W();
  if (W != Teuchos::null) {
    const double alpha = inArgs.get_alpha();
    const double beta = inArgs.get_beta();
    Epetra_CrsMatrix &crsW = Teuchos::dyn_cast<Epetra_CrsMatrix>(*W);
    const int dim = 2;
    double values_1[2];
    double values_2[2];
    int indices[] = {0,1};

    Sacado::Fad::DFad<double> x1(dim,0,(*x)[0]);
    Sacado::Fad::DFad<double> x2(dim,1,(*x)[1]);
    Sacado::Fad::DFad<double> f1;
    Sacado::Fad::DFad<double> f2;

    evalVField(x1,x2,f1,f2);

    values_1[0] = alpha * f1.fastAccessDx(0) - beta;
    values_1[1] = alpha * f1.fastAccessDx(1);
    values_2[0] = alpha * f2.fastAccessDx(0);
    values_2[0] = alpha * f2.fastAccessDx(1) - beta;
   
    crsW.ReplaceGlobalValues(0,dim,values_1,indices);
    crsW.ReplaceGlobalValues(1,dim,values_2,indices);
  }
}