45 #include "Teuchos_Assert.hpp" 51 x_map = Teuchos::rcp(
new Epetra_Map(2, 0, *comm));
67 p_map = Teuchos::rcp(
new Epetra_LocalMap(1, 0, *comm));
74 p_names = Teuchos::rcp(
new Teuchos::Array<std::string>(1));
75 (*p_names)[0] =
"alpha";
78 graph = Teuchos::rcp(
new Epetra_CrsGraph(Copy, *
x_map, 2));
80 indices[0] = 0; indices[1] = 1;
82 graph->InsertGlobalIndices(0, 2, indices);
84 graph->InsertGlobalIndices(1, 2, indices);
85 graph->FillComplete();
86 graph->OptimizeStorage();
91 Teuchos::RCP<const Epetra_Map>
97 Teuchos::RCP<const Epetra_Map>
103 Teuchos::RCP<const Epetra_Map>
106 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
109 "Error! SimpleME::get_p_map(): " <<
110 "Invalid parameter index l = " << l << std::endl);
115 Teuchos::RCP<const Teuchos::Array<std::string> >
118 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
121 "Error! SimpleME::get_p_names(): " <<
122 "Invalid parameter index l = " << l << std::endl);
127 Teuchos::RCP<const Epetra_Vector>
133 Teuchos::RCP<const Epetra_Vector>
136 TEUCHOS_TEST_FOR_EXCEPTION(l != 0,
139 "Error! SimpleME::get_p_init(): " <<
140 "Invalid parameter index l = " << l << std::endl);
145 Teuchos::RCP<Epetra_Operator>
148 Teuchos::RCP<Epetra_CrsMatrix> A =
149 Teuchos::rcp(
new Epetra_CrsMatrix(Copy, *
graph));
151 A->OptimizeStorage();
155 EpetraExt::ModelEvaluator::InArgs
159 inArgs.setModelEvalDescription(
"Simple Model Evaluator");
162 inArgs.setSupports(IN_ARG_x,
true);
166 inArgs.setSupports(IN_ARG_x_sg,
true);
167 inArgs.setSupports(IN_ARG_p_sg, 0,
true);
168 inArgs.setSupports(IN_ARG_sg_basis,
true);
169 inArgs.setSupports(IN_ARG_sg_quadrature,
true);
170 inArgs.setSupports(IN_ARG_sg_expansion,
true);
175 EpetraExt::ModelEvaluator::OutArgs
178 OutArgsSetup outArgs;
179 outArgs.setModelEvalDescription(
"Simple Model Evaluator");
182 outArgs.set_Np_Ng(1, 0);
183 outArgs.setSupports(OUT_ARG_f,
true);
184 outArgs.setSupports(OUT_ARG_W,
true);
187 outArgs.setSupports(OUT_ARG_f_sg,
true);
188 outArgs.setSupports(OUT_ARG_W_sg,
true);
201 Teuchos::RCP<const Epetra_Vector>
x = inArgs.get_x();
202 if (
x != Teuchos::null)
204 double x0 = (*x_overlapped)[0];
205 double x1 = (*x_overlapped)[1];
208 Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(0);
209 if (p == Teuchos::null)
217 Teuchos::RCP<Epetra_Vector>
f = outArgs.get_f();
218 if (
f != Teuchos::null) {
221 if (
x_map->MyGID(0)) {
224 f->ReplaceGlobalValues(1, &v, &row);
226 if (
x_map->MyGID(1)) {
229 f->ReplaceGlobalValues(1, &v, &row);
236 Teuchos::RCP<Epetra_Operator> W = outArgs.get_W();
237 if (W != Teuchos::null) {
238 Teuchos::RCP<Epetra_CrsMatrix> jac =
239 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W,
true);
240 int indices[2] = { 0, 1 };
242 if (
x_map->MyGID(0)) {
243 values[0] = -1.0; values[1] = 0.0;
244 jac->ReplaceGlobalValues(0, 2, values, indices);
246 if (
x_map->MyGID(1)) {
247 values[0] = -1.0; values[1] = 2.0*x1;
248 jac->ReplaceGlobalValues(1, 2, values, indices);
257 Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> > basis =
258 inArgs.get_sg_basis();
259 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> > expn =
260 inArgs.get_sg_expansion();
263 InArgs::sg_const_vector_t x_sg = inArgs.get_x_sg();
265 if (x_sg != Teuchos::null) {
266 for (
int i=0; i<basis->size(); i++) {
268 x0_sg[i] = (*x_overlapped)[0];
269 x1_sg[i] = (*x_overlapped)[1];
274 InArgs::sg_const_vector_t p_sg = inArgs.get_p_sg(0);
276 if (p_sg != Teuchos::null) {
277 for (
int i=0; i<basis->size(); i++) {
278 a_sg[i] = (*p_sg)[i][0];
285 OutArgs::sg_vector_t f_sg = outArgs.get_f_sg();
287 if (f_sg != Teuchos::null) {
289 if (
x_map->MyGID(0)) {
291 expn->times(tmp0_sg, a_sg, a_sg);
292 expn->minus(tmp1_sg, tmp0_sg, x0_sg);
293 for (
int i=0; i<basis->size(); i++)
294 (*f_sg)[i].ReplaceGlobalValues(1, &tmp1_sg[i], &row);
296 if (
x_map->MyGID(1)) {
298 expn->times(tmp0_sg, x1_sg, x1_sg);
299 expn->minus(tmp1_sg, tmp0_sg, x0_sg);
300 for (
int i=0; i<basis->size(); i++)
301 (*f_sg)[i].ReplaceGlobalValues(1, &tmp1_sg[i], &row);
308 OutArgs::sg_operator_t W_sg = outArgs.get_W_sg();
309 if (W_sg != Teuchos::null) {
310 for (
int i=0; i<basis->size(); i++) {
311 Teuchos::RCP<Epetra_CrsMatrix> jac =
312 Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(W_sg->getCoeffPtr(i),
true);
313 int indices[2] = { 0, 1 };
315 if (
x_map->MyGID(0)) {
321 jac->ReplaceGlobalValues(0, 2, values, indices);
323 if (
x_map->MyGID(1)) {
328 values[1] = 2.0*x1_sg[i];
329 jac->ReplaceGlobalValues(1, 2, values, indices);
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
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.
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.
Teuchos::RCP< Epetra_Map > x_overlapped_map
Overlapped solution vector map.