43 #include "Teuchos_Assert.hpp" 50 template <
typename ScalarA,
typename ScalarX>
51 void func(
const ScalarA& a,
const ScalarX
x[2], ScalarX
y[2]) {
53 y[1] =
x[1]*
x[1] -
x[0];
84 p_names = Teuchos::rcp(
new Teuchos::Array<std::string>(1));
85 (*p_names)[0] =
"alpha";
90 indices[0] = 0; indices[1] = 1;
92 graph->InsertGlobalIndices(0, 2, indices);
94 graph->InsertGlobalIndices(1, 2, indices);
95 graph->FillComplete();
96 graph->OptimizeStorage();
101 Teuchos::RCP<const Epetra_Map>
107 Teuchos::RCP<const Epetra_Map>
113 Teuchos::RCP<const Epetra_Map>
116 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
119 "Error! SimpleME::get_p_map(): " <<
120 "Invalid parameter index l = " << l << std::endl);
125 Teuchos::RCP<const Teuchos::Array<std::string> >
128 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
131 "Error! SimpleME::get_p_names(): " <<
132 "Invalid parameter index l = " << l << std::endl);
137 Teuchos::RCP<const Epetra_Vector>
143 Teuchos::RCP<const Epetra_Vector>
146 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
149 "Error! SimpleME::get_p_init(): " <<
150 "Invalid parameter index l = " << l << std::endl);
155 Teuchos::RCP<Epetra_Operator>
158 Teuchos::RCP<Epetra_CrsMatrix> A =
161 A->OptimizeStorage();
165 EpetraExt::ModelEvaluator::InArgs
169 inArgs.setModelEvalDescription(
"Simple Model Evaluator");
172 inArgs.setSupports(IN_ARG_x,
true);
176 inArgs.setSupports(IN_ARG_x_sg,
true);
177 inArgs.setSupports(IN_ARG_p_sg, 0,
true);
178 inArgs.setSupports(IN_ARG_sg_basis,
true);
179 inArgs.setSupports(IN_ARG_sg_quadrature,
true);
180 inArgs.setSupports(IN_ARG_sg_expansion,
true);
185 EpetraExt::ModelEvaluator::OutArgs
188 OutArgsSetup outArgs;
189 outArgs.setModelEvalDescription(
"Simple Model Evaluator");
192 outArgs.set_Np_Ng(1, 0);
193 outArgs.setSupports(OUT_ARG_f,
true);
194 outArgs.setSupports(OUT_ARG_W,
true);
197 outArgs.setSupports(OUT_ARG_f_sg,
true);
198 outArgs.setSupports(OUT_ARG_W_sg,
true);
211 Teuchos::RCP<const Epetra_Vector>
x = inArgs.get_x();
212 if (
x != Teuchos::null) {
214 double x0 = (*x_overlapped)[0];
215 double x1 = (*x_overlapped)[1];
218 Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(0);
219 if (p == Teuchos::null)
227 Teuchos::RCP<Epetra_Vector>
f = outArgs.get_f();
228 if (
f != Teuchos::null) {
229 double x[2] = { x0, x1 };
233 if (
x_map->MyGID(0)) {
235 f->ReplaceGlobalValues(1, &
y[0], &row);
237 if (
x_map->MyGID(1)) {
239 f->ReplaceGlobalValues(1, &
y[1], &row);
246 Teuchos::RCP<Epetra_Operator> W = outArgs.get_W();
247 if (W != Teuchos::null) {
248 typedef Sacado::Fad::SFad<double,2>
fad_type;
254 Teuchos::RCP<Epetra_CrsMatrix> jac =
256 int indices[2] = { 0, 1 };
257 if (
x_map->MyGID(0)) {
261 if (
x_map->MyGID(1)) {
263 jac->ReplaceGlobalValues(row, 2,
y[1].
dx(), indices);
273 InArgs::sg_const_vector_t x_sg = inArgs.get_x_sg();
274 if (x_sg != Teuchos::null) {
277 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> > basis =
278 inArgs.get_sg_basis();
279 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expn =
280 inArgs.get_sg_expansion();
285 for (
int i=0; i<basis->size(); i++) {
287 x0.fastAccessCoeff(i) = (*x_overlapped)[0];
288 x1.fastAccessCoeff(i) = (*x_overlapped)[1];
292 InArgs::sg_const_vector_t p_sg = inArgs.get_p_sg(0);
294 if (p_sg != Teuchos::null) {
295 for (
int i=0; i<basis->size(); i++) {
296 a.fastAccessCoeff(i) = (*p_sg)[i][0];
303 OutArgs::sg_vector_t f_sg = outArgs.get_f_sg();
304 if (f_sg != Teuchos::null) {
309 if (
x_map->MyGID(0)) {
311 for (
int i=0; i<basis->size(); i++) {
312 double c =
y[0].coeff(i);
313 (*f_sg)[i].ReplaceGlobalValues(1, &c, &row);
316 if (
x_map->MyGID(1)) {
318 for (
int i=0; i<basis->size(); i++) {
319 double c =
y[1].coeff(i);
320 (*f_sg)[i].ReplaceGlobalValues(1, &c, &row);
328 OutArgs::sg_operator_t W_sg = outArgs.get_W_sg();
329 if (W_sg != Teuchos::null) {
330 typedef Sacado::Fad::SFad<pce_type,2>
fad_type;
336 for (
int i=0; i<basis->size(); i++) {
337 Teuchos::RCP<Epetra_CrsMatrix> jac =
340 int indices[2] = { 0, 1 };
341 if (
x_map->MyGID(0)) {
343 double values[2] = {
y[0].dx(0).coeff(i),
y[0].dx(1).coeff(i) };
344 jac->ReplaceGlobalValues(row, 2, values, indices);
346 if (
x_map->MyGID(1)) {
348 double values[2] = {
y[1].dx(0).coeff(i),
y[1].dx(1).coeff(i) };
349 jac->ReplaceGlobalValues(row, 2, values, indices);
Stokhos::StandardStorage< int, double > storage_type
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
Sacado::ETPCE::OrthogPoly< double, Stokhos::StandardStorage< int, double > > pce_type
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
Teuchos::RCP< Epetra_Vector > x_init
Initial guess.
SimpleME(const Teuchos::RCP< const Epetra_Comm > &comm)
Constructor.
Teuchos::RCP< Epetra_Vector > p_init
Initial parameters.
Teuchos::RCP< Epetra_Vector > x_overlapped
Overlapped solution vector.
Teuchos::RCP< Epetra_Map > x_map
Solution vector map.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
InArgs createInArgs() const
Create InArgs.
Sacado::Fad::DFad< double > fad_type
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Teuchos::RCP< Epetra_Map > p_map
Parameter vector map.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
Teuchos::RCP< Teuchos::Array< std::string > > p_names
Parameter names.
Teuchos::RCP< Epetra_CrsGraph > graph
Jacobian graph.
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
Teuchos::RCP< Epetra_Import > importer
Importer to overlapped distribution.
OutArgs createOutArgs() const
Create OutArgs.
virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
Sacado::Fad::DFad< double > fad_type
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Teuchos::RCP< Epetra_Map > x_overlapped_map
Overlapped solution vector map.