Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_MPModelEvaluator.cpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
43
44#include <algorithm>
45#include "Teuchos_Assert.hpp"
46#include "EpetraExt_BlockUtility.h"
47#include "EpetraExt_BlockMultiVector.h"
52
54 const Teuchos::RCP<EpetraExt::ModelEvaluator>& me_,
55 const Teuchos::RCP<const EpetraExt::MultiComm>& mp_comm_,
56 const Teuchos::RCP<const Epetra_Map>& mp_block_map_,
57 const Teuchos::RCP<Teuchos::ParameterList>& params_)
58 : me(me_),
59 num_mp_blocks(mp_block_map_->NumMyElements()),
60 mp_comm(mp_comm_),
61 mp_block_map(mp_block_map_),
62 params(params_),
63 supports_x(false),
64 x_map(me->get_x_map()),
65 f_map(me->get_f_map()),
66 mp_x_map(),
67 mp_f_map(),
68 num_p(0),
69 num_p_mp(0),
70 mp_p_index_map(),
71 mp_p_map(),
72 mp_p_names(),
73 num_g(0),
74 num_g_mp(0),
75 mp_g_index_map(),
76 mp_g_map(),
77 W_mp_blocks(),
78 mp_p_init(),
79 my_W(),
80 my_x()
81{
82 if (x_map != Teuchos::null)
83 supports_x = true;
84
85 if (supports_x) {
86
87 // Create block MP x and f maps
88 mp_x_map =
89 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
91
92 mp_f_map =
93 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
95
96 // Create default mp_x_init
97 mp_x_init = Teuchos::rcp(new ProductEpetraVector(
99 for (unsigned int i=0; i<num_mp_blocks; i++)
100 (*mp_x_init)[i] = *(me->get_x_init());
101
102 mp_x_dot_init = Teuchos::rcp(new ProductEpetraVector(
104 if (me->get_x_dot_init() != Teuchos::null) {
105 for (unsigned int i=0; i<num_mp_blocks; i++)
106 (*mp_x_dot_init)[i] = *(me->get_x_dot_init());
107 }
108 // Preconditioner needs an x
109 my_x = Teuchos::rcp(new Epetra_Vector(*mp_x_map));
110
111 W_mp_blocks =
112 Teuchos::rcp(new Stokhos::ProductEpetraOperator(
114 for (unsigned int i=0; i<num_mp_blocks; i++)
115 W_mp_blocks->setCoeffPtr(i, me->create_W());
116 }
117
118 // Parameters -- The idea here is to add new parameter vectors
119 // for the stochastic Galerkin components of the parameters
120
121 InArgs me_inargs = me->createInArgs();
122 num_p = me_inargs.Np();
123
124 // Get the p_mp's supported and build index map
125 for (int i=0; i<num_p; i++) {
126 if (me_inargs.supports(IN_ARG_p_mp, i))
127 mp_p_index_map.push_back(i);
128 }
129 num_p_mp = mp_p_index_map.size();
130
131 mp_p_map.resize(num_p_mp);
132 mp_p_names.resize(num_p_mp);
133 mp_p_init.resize(num_p_mp);
134
135 // Create parameter maps, names, and initial values
136 for (int i=0; i<num_p_mp; i++) {
137 Teuchos::RCP<const Epetra_Map> p_map = me->get_p_map(mp_p_index_map[i]);
138 mp_p_map[i] =
139 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
140 *p_map, *mp_block_map, *mp_comm));
141
142 Teuchos::RCP<const Teuchos::Array<std::string> > p_names =
143 me->get_p_names(mp_p_index_map[i]);
144 if (p_names != Teuchos::null) {
145 mp_p_names[i] =
146 Teuchos::rcp(new Teuchos::Array<std::string>(num_mp_blocks*(p_names->size())));
147 for (int j=0; j<p_names->size(); j++) {
148 std::stringstream ss;
149 ss << (*p_names)[j] << " -- MP Coefficient " << i;
150 (*mp_p_names[i])[j] = ss.str();
151 }
152 }
153
154 // Create default mp_p_init
155 mp_p_init[i] =
156 Teuchos::rcp(new ProductEpetraVector(
157 mp_block_map, p_map, mp_p_map[i], mp_comm));
158 mp_p_init[i]->init(0.0);
159 // for (unsigned int j=0; j<num_mp_blocks; j++)
160 // (*mp_p_init[i])[j] = *(me->get_p_init(i));
161
162 }
163
164 // Responses -- The idea here is to add new parameter vectors
165 // for the stochastic Galerkin components of the respones
166
167 // Get number of MP parameters model supports derivatives w.r.t.
168 OutArgs me_outargs = me->createOutArgs();
169 num_g = me_outargs.Ng();
170
171 // Get the g_mp's supported and build index map
172 for (int i=0; i<num_g; i++) {
173 if (me_outargs.supports(OUT_ARG_g_mp, i))
174 mp_g_index_map.push_back(i);
175 }
176 num_g_mp = mp_g_index_map.size();
177
178 mp_g_map.resize(num_g_mp);
179
180 // Create response maps
181 for (int i=0; i<num_g_mp; i++) {
182 Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(mp_g_index_map[i]);
183 mp_g_map[i] =
184 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
185 *g_map, *mp_block_map, *mp_comm));
186 }
187}
188
189// Overridden from EpetraExt::ModelEvaluator
190
191Teuchos::RCP<const Epetra_Map>
193{
194 return mp_x_map;
195}
196
197Teuchos::RCP<const Epetra_Map>
199{
200 return mp_f_map;
201}
202
203Teuchos::RCP<const Epetra_Map>
205{
206 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_mp, std::logic_error,
207 "Error! Invalid p map index " << l);
208 if (l < num_p)
209 return me->get_p_map(l);
210 else
211 return mp_p_map[l-num_p];
212
213 return Teuchos::null;
214}
215
216Teuchos::RCP<const Epetra_Map>
218{
219 TEUCHOS_TEST_FOR_EXCEPTION(j < 0 || j >= num_g_mp, std::logic_error,
220 "Error! Invalid g map index " << j);
221 return mp_g_map[j];
222}
223
224Teuchos::RCP<const Teuchos::Array<std::string> >
226{
227 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_mp, std::logic_error,
228 "Error! Invalid p map index " << l);
229 if (l < num_p)
230 return me->get_p_names(l);
231 else
232 return mp_p_names[l-num_p];
233
234 return Teuchos::null;
235}
236
237Teuchos::RCP<const Epetra_Vector>
239{
240 return mp_x_init->getBlockVector();
241}
242
243Teuchos::RCP<const Epetra_Vector>
245{
246 return mp_x_dot_init->getBlockVector();
247}
248
249Teuchos::RCP<const Epetra_Vector>
251{
252 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_mp, std::logic_error,
253 "Error! Invalid p map index " << l);
254 if (l < num_p)
255 return me->get_p_init(l);
256 else if (l < num_p + num_p_mp)
257 return mp_p_init[l-num_p]->getBlockVector();
258
259 return Teuchos::null;
260}
261
262Teuchos::RCP<Epetra_Operator>
264{
265 if (supports_x) {
266 my_W =
267 Teuchos::rcp(new Stokhos::BlockDiagonalOperator(mp_comm, num_mp_blocks,
268 x_map, f_map,
269 mp_x_map, mp_f_map));
270 my_W->setupOperator(W_mp_blocks);
271
272 return my_W;
273 }
274
275 return Teuchos::null;
276}
277
278Teuchos::RCP<EpetraExt::ModelEvaluator::Preconditioner>
280{
281 if (supports_x) {
282 Teuchos::RCP<Teuchos::ParameterList> mpPrecParams =
283 Teuchos::rcp(&(params->sublist("MP Preconditioner")), false);
284 Stokhos::MPPreconditionerFactory mp_prec_factory(mpPrecParams);
285 Teuchos::RCP<Epetra_Operator> precOp =
286 mp_prec_factory.build(mp_comm, num_mp_blocks, x_map, mp_x_map);
287 return Teuchos::rcp(new EpetraExt::ModelEvaluator::Preconditioner(precOp,
288 true));
289 }
290
291 return Teuchos::null;
292}
293
294Teuchos::RCP<Epetra_Operator>
296{
297 TEUCHOS_TEST_FOR_EXCEPTION(j < 0 || j >= num_g_mp || !supports_x,
298 std::logic_error,
299 "Error: dg/dx index " << j << " is not supported!");
300
301 int jj = mp_g_index_map[j];
302 Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(jj);
303 Teuchos::RCP< Stokhos::ProductEpetraOperator > mp_blocks;
304 OutArgs me_outargs = me->createOutArgs();
305 DerivativeSupport ds = me_outargs.supports(OUT_ARG_DgDx_mp, jj);
306 if (ds.supports(DERIV_LINEAR_OP)) {
307 mp_blocks =
308 Teuchos::rcp(new Stokhos::ProductEpetraOperator(
309 mp_block_map, x_map, g_map, mp_g_map[j], mp_comm));
310 for (unsigned int i=0; i<num_mp_blocks; i++)
311 mp_blocks->setCoeffPtr(i, me->create_DgDx_op(jj));
312 }
313 else if (ds.supports(DERIV_MV_BY_COL)) {
314 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
315 Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
316 mp_block_map, g_map, mp_g_map[j],
317 mp_comm, x_map->NumMyElements()));
318 mp_blocks =
320 mp_mv_blocks, false));
321 }
322 else if (ds.supports(DERIV_TRANS_MV_BY_ROW)) {
323 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
324 Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
325 mp_block_map, x_map, mp_x_map,
326 mp_comm, g_map->NumMyElements()));
327 mp_blocks =
329 mp_mv_blocks, true));
330 }
331 else
332 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
333 "Error! me_outargs.supports(OUT_ARG_DgDx_mp, " << j
334 << ").none() is true!");
335
336 Teuchos::RCP<Stokhos::BlockDiagonalOperator> dgdx_mp =
337 Teuchos::rcp(new Stokhos::BlockDiagonalOperator(
338 mp_comm, num_mp_blocks, x_map, g_map,
339 mp_x_map, mp_g_map[j]));
340 dgdx_mp->setupOperator(mp_blocks);
341 return dgdx_mp;
342}
343
344Teuchos::RCP<Epetra_Operator>
346{
347 TEUCHOS_TEST_FOR_EXCEPTION(j < 0 || j >= num_g_mp || !supports_x,
348 std::logic_error,
349 "Error: dg/dx_dot index " << j << " is not supported!");
350
351 int jj = mp_g_index_map[j];
352 Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(jj);
353 Teuchos::RCP< Stokhos::ProductEpetraOperator > mp_blocks;
354 OutArgs me_outargs = me->createOutArgs();
355 DerivativeSupport ds = me_outargs.supports(OUT_ARG_DgDx_dot_mp, jj);
356 if (ds.supports(DERIV_LINEAR_OP)) {
357 mp_blocks =
358 Teuchos::rcp(new Stokhos::ProductEpetraOperator(
359 mp_block_map, x_map, g_map, mp_g_map[j],
360 mp_comm));
361 for (unsigned int i=0; i<num_mp_blocks; i++)
362 mp_blocks->setCoeffPtr(i, me->create_DgDx_dot_op(jj));
363 }
364 else if (ds.supports(DERIV_MV_BY_COL)) {
365 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
366 Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
367 mp_block_map, g_map, mp_g_map[j],
368 mp_comm, x_map->NumMyElements()));
369 mp_blocks =
371 mp_mv_blocks, false));
372 }
373 else if (ds.supports(DERIV_TRANS_MV_BY_ROW)) {
374 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
375 Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
376 mp_block_map, x_map, mp_x_map,
377 mp_comm, g_map->NumMyElements()));
378 mp_blocks =
380 mp_mv_blocks, true));
381 }
382 else
383 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
384 "Error! me_outargs.supports(OUT_ARG_DgDx_dot_mp, "
385 << j << ").none() is true!");
386
387 Teuchos::RCP<Stokhos::BlockDiagonalOperator> dgdx_dot_mp =
388 Teuchos::rcp(new Stokhos::BlockDiagonalOperator(
389 mp_comm, num_mp_blocks, x_map, g_map,
390 mp_x_map, mp_g_map[j]));
391 dgdx_dot_mp->setupOperator(mp_blocks);
392 return dgdx_dot_mp;
393}
394
395Teuchos::RCP<Epetra_Operator>
397{
398 TEUCHOS_TEST_FOR_EXCEPTION(
399 j < 0 || j >= num_g_mp || i < 0 || i >= num_p+num_p_mp,
400 std::logic_error,
401 "Error: dg/dp index " << j << "," << i << " is not supported!");
402
403 OutArgs me_outargs = me->createOutArgs();
404 int jj = mp_g_index_map[j];
405 Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(jj);
406 if (i < num_p) {
407 if (me_outargs.supports(OUT_ARG_DgDp_mp,jj,i).supports(DERIV_LINEAR_OP)) {
408 Teuchos::RCP<Stokhos::ProductEpetraOperator> mp_blocks =
409 Teuchos::rcp(new Stokhos::ProductEpetraOperator(
410 mp_block_map, me->get_p_map(i), g_map,
411 mp_g_map[j], mp_comm));
412 for (unsigned int l=0; l<num_mp_blocks; l++)
413 mp_blocks->setCoeffPtr(l, me->create_DgDp_op(i,j));
414 return mp_blocks;
415 }
416 else
417 TEUCHOS_TEST_FOR_EXCEPTION(
418 true, std::logic_error,
419 "Error: Underlying model evaluator must support DERIV_LINER_OP " <<
420 "to create operator for dg/dp index " << j << "," << i << "!");
421 }
422 else {
423 int ii = mp_p_index_map[i-num_p];
424 Teuchos::RCP<const Epetra_Map> p_map = me->get_p_map(ii);
425 Teuchos::RCP< Stokhos::ProductEpetraOperator> mp_blocks;
426 DerivativeSupport ds = me_outargs.supports(OUT_ARG_DgDp_mp,jj,ii);
427 if (ds.supports(DERIV_LINEAR_OP)) {
428 mp_blocks =
429 Teuchos::rcp(new Stokhos::ProductEpetraOperator(
430 mp_block_map, p_map, g_map,
431 mp_g_map[j], mp_comm));
432 for (unsigned int l=0; l<num_mp_blocks; l++)
433 mp_blocks->setCoeffPtr(l, me->create_DgDp_op(jj,ii));
434 }
435 else if (ds.supports(DERIV_MV_BY_COL)) {
436 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
437 Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
438 mp_block_map, g_map, mp_g_map[j],
439 mp_comm, p_map->NumMyElements()));
440 mp_blocks =
442 mp_mv_blocks, false));
443 }
444 else if (ds.supports(DERIV_TRANS_MV_BY_ROW)) {
445 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
446 Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
447 mp_block_map, p_map, mp_p_map[i-num_p],
448 mp_comm, g_map->NumMyElements()));
449 mp_blocks =
451 mp_mv_blocks, true));
452 }
453 else
454 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
455 "Error! me_outargs.supports(OUT_ARG_DgDp_mp, " << jj
456 << "," << ii << ").none() is true!");
457
458 Teuchos::RCP<Stokhos::BlockDiagonalOperator> dgdp_mp =
459 Teuchos::rcp(new Stokhos::BlockDiagonalOperator(
460 mp_comm, num_mp_blocks, p_map, g_map,
461 mp_p_map[i-num_p], mp_g_map[j]));
462 dgdp_mp->setupOperator(mp_blocks);
463 return dgdp_mp;
464 }
465
466 return Teuchos::null;
467}
468
469Teuchos::RCP<Epetra_Operator>
471{
472 TEUCHOS_TEST_FOR_EXCEPTION(i < 0 || i >= num_p+num_p_mp,
473 std::logic_error,
474 "Error: df/dp index " << i << " is not supported!");
475
476 OutArgs me_outargs = me->createOutArgs();
477 if (i < num_p) {
478 if (me_outargs.supports(OUT_ARG_DfDp_mp,i).supports(DERIV_LINEAR_OP)) {
479 Teuchos::RCP<Stokhos::ProductEpetraOperator> mp_blocks =
480 Teuchos::rcp(new Stokhos::ProductEpetraOperator(
481 mp_block_map, me->get_p_map(i), me->get_f_map(),
482 mp_f_map, mp_comm));
483 for (unsigned int l=0; l<num_mp_blocks; l++)
484 mp_blocks->setCoeffPtr(l, me->create_DfDp_op(i));
485 return mp_blocks;
486 }
487 else
488 TEUCHOS_TEST_FOR_EXCEPTION(
489 true, std::logic_error,
490 "Error: Underlying model evaluator must support DERIV_LINER_OP " <<
491 "to create operator for df/dp index " << i << "!");
492 }
493 else {
494 int ii = mp_p_index_map[i-num_p];
495 Teuchos::RCP<const Epetra_Map> p_map = me->get_p_map(ii);
496 Teuchos::RCP<Stokhos::ProductEpetraOperator> mp_blocks;
497 DerivativeSupport ds = me_outargs.supports(OUT_ARG_DfDp_mp,ii);
498 if (ds.supports(DERIV_LINEAR_OP)) {
499 mp_blocks =
500 Teuchos::rcp(new Stokhos::ProductEpetraOperator(
501 mp_block_map, me->get_p_map(ii), me->get_f_map(),
502 mp_f_map, mp_comm));
503 for (unsigned int l=0; l<num_mp_blocks; l++)
504 mp_blocks->setCoeffPtr(l, me->create_DfDp_op(ii));
505 }
506 else if (ds.supports(DERIV_MV_BY_COL)) {
507 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
508 Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
509 mp_block_map, f_map, mp_f_map,
510 mp_comm, p_map->NumMyElements()));
511 mp_blocks =
513 mp_mv_blocks, false));
514 }
515 else if (ds.supports(DERIV_TRANS_MV_BY_ROW)) {
516 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
517 Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
518 mp_block_map, p_map, mp_p_map[i-num_p],
519 mp_comm, f_map->NumMyElements()));
520 mp_blocks =
522 mp_mv_blocks, true));
523 }
524 else
525 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
526 "Error! me_outargs.supports(OUT_ARG_DfDp_mp, " << ii
527 << ").none() is true!");
528
529
530 Teuchos::RCP<Stokhos::BlockDiagonalOperator> dfdp_mp =
531 Teuchos::rcp(new Stokhos::BlockDiagonalOperator(
532 mp_comm, num_mp_blocks,
533 p_map, f_map, mp_p_map[i-num_p], mp_f_map));
534 dfdp_mp->setupOperator(mp_blocks);
535 return dfdp_mp;
536 }
537
538 return Teuchos::null;
539}
540
541EpetraExt::ModelEvaluator::InArgs
543{
544 InArgsSetup inArgs;
545 InArgs me_inargs = me->createInArgs();
546
547 inArgs.setModelEvalDescription(this->description());
548 inArgs.set_Np(num_p + num_p_mp);
549 inArgs.setSupports(IN_ARG_x_dot, me_inargs.supports(IN_ARG_x_dot));
550 inArgs.setSupports(IN_ARG_x, me_inargs.supports(IN_ARG_x));
551 inArgs.setSupports(IN_ARG_t, me_inargs.supports(IN_ARG_t));
552 inArgs.setSupports(IN_ARG_alpha, me_inargs.supports(IN_ARG_alpha));
553 inArgs.setSupports(IN_ARG_beta, me_inargs.supports(IN_ARG_beta));
554
555 return inArgs;
556}
557
558EpetraExt::ModelEvaluator::OutArgs
560{
561 OutArgsSetup outArgs;
562 OutArgs me_outargs = me->createOutArgs();
563
564 outArgs.setModelEvalDescription(this->description());
565 outArgs.set_Np_Ng(num_p + num_p_mp, num_g_mp);
566 outArgs.setSupports(OUT_ARG_f, me_outargs.supports(OUT_ARG_f));
567 outArgs.setSupports(OUT_ARG_W, me_outargs.supports(OUT_ARG_W));
568 outArgs.setSupports(OUT_ARG_WPrec, me_outargs.supports(OUT_ARG_W));
569 for (int j=0; j<num_p; j++)
570 outArgs.setSupports(OUT_ARG_DfDp, j,
571 me_outargs.supports(OUT_ARG_DfDp, j));
572 for (int j=0; j<num_p_mp; j++)
573 if (!me_outargs.supports(OUT_ARG_DfDp_mp, mp_p_index_map[j]).none())
574 outArgs.setSupports(OUT_ARG_DfDp, j+num_p, DERIV_LINEAR_OP);
575 for (int i=0; i<num_g_mp; i++) {
576 int ii = mp_g_index_map[i];
577 if (!me_outargs.supports(OUT_ARG_DgDx_dot_mp, ii).none())
578 outArgs.setSupports(OUT_ARG_DgDx_dot, i, DERIV_LINEAR_OP);
579 if (!me_outargs.supports(OUT_ARG_DgDx_mp, ii).none())
580 outArgs.setSupports(OUT_ARG_DgDx, i, DERIV_LINEAR_OP);
581 for (int j=0; j<num_p; j++)
582 outArgs.setSupports(OUT_ARG_DgDp, i, j,
583 me_outargs.supports(OUT_ARG_DgDp_mp, ii, j));
584 for (int j=0; j<num_p_mp; j++)
585 if (!me_outargs.supports(OUT_ARG_DgDp_mp, ii, mp_p_index_map[j]).none())
586 outArgs.setSupports(OUT_ARG_DgDp, i, j+num_p, DERIV_LINEAR_OP);
587 }
588
589 return outArgs;
590}
591
592void
594 const OutArgs& outArgs) const
595{
596 // Get the input arguments
597 Teuchos::RCP<const Epetra_Vector> x;
598 if (inArgs.supports(IN_ARG_x)) {
599 x = inArgs.get_x();
600 if (x != Teuchos::null)
601 *my_x = *x;
602 }
603 Teuchos::RCP<const Epetra_Vector> x_dot;
604 if (inArgs.supports(IN_ARG_x_dot))
605 x_dot = inArgs.get_x_dot();
606
607 // Get the output arguments
608 EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> f_out;
609 if (outArgs.supports(OUT_ARG_f))
610 f_out = outArgs.get_f();
611 Teuchos::RCP<Epetra_Operator> W_out;
612 if (outArgs.supports(OUT_ARG_W))
613 W_out = outArgs.get_W();
614 Teuchos::RCP<Epetra_Operator> WPrec_out;
615 if (outArgs.supports(OUT_ARG_WPrec))
616 WPrec_out = outArgs.get_WPrec();
617
618 // Check if we are using the "matrix-free" method for W and we are
619 // computing a preconditioner.
620 bool eval_prec = (W_out == Teuchos::null && WPrec_out != Teuchos::null);
621
622 // Here we are assuming a full W fill occurred previously which we can use
623 // for the preconditioner. Given the expense of computing the MP W blocks
624 // this saves significant computational cost
625 if (eval_prec) {
626 Teuchos::RCP<Stokhos::MPPreconditioner> W_prec =
627 Teuchos::rcp_dynamic_cast<Stokhos::MPPreconditioner>(WPrec_out, true);
628 W_prec->setupPreconditioner(my_W, *my_x);
629
630 // We can now quit unless a fill of f, g, or dg/dp was also requested
631 bool done = (f_out == Teuchos::null);
632 for (int i=0; i<outArgs.Ng(); i++) {
633 done = done && (outArgs.get_g(i) == Teuchos::null);
634 for (int j=0; j<outArgs.Np(); j++)
635 if (!outArgs.supports(OUT_ARG_DgDp, i, j).none())
636 done = done && (outArgs.get_DgDp(i,j).isEmpty());
637 }
638 if (done)
639 return;
640 }
641
642 // Create underlying inargs
643 InArgs me_inargs = me->createInArgs();
644 if (x != Teuchos::null) {
645 Teuchos::RCP<Stokhos::ProductEpetraVector> x_mp =
646 create_x_mp(View, x.get());
647 me_inargs.set_x_mp(x_mp);
648 }
649 if (x_dot != Teuchos::null) {
650 Teuchos::RCP<Stokhos::ProductEpetraVector> x_dot_mp =
651 create_x_mp(View, x_dot.get());
652 me_inargs.set_x_dot_mp(x_dot_mp);
653 }
654 if (me_inargs.supports(IN_ARG_alpha))
655 me_inargs.set_alpha(inArgs.get_alpha());
656 if (me_inargs.supports(IN_ARG_beta))
657 me_inargs.set_beta(inArgs.get_beta());
658 if (me_inargs.supports(IN_ARG_t))
659 me_inargs.set_t(inArgs.get_t());
660
661 // Pass parameters
662 for (int i=0; i<num_p; i++)
663 me_inargs.set_p(i, inArgs.get_p(i));
664 for (int i=0; i<num_p_mp; i++) {
665 Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(i+num_p);
666
667 // We always need to pass in the MP parameters, so just use
668 // the initial parameters if it is NULL
669 if (p == Teuchos::null)
670 p = mp_p_init[i]->getBlockVector();
671
672 // Convert block p to MP polynomial
673 Teuchos::RCP<Stokhos::ProductEpetraVector> p_mp =
674 create_p_mp(mp_p_index_map[i], View, p.get());
675 me_inargs.set_p_mp(mp_p_index_map[i], p_mp);
676 }
677
678 // Create underlying outargs
679 OutArgs me_outargs = me->createOutArgs();
680
681 // f
682 if (f_out != Teuchos::null) {
683 Teuchos::RCP<Stokhos::ProductEpetraVector> f_mp =
684 create_f_mp(View, f_out.get());
685 me_outargs.set_f_mp(f_mp);
686 }
687
688 // W
689 if (W_out != Teuchos::null && !eval_prec)
690 me_outargs.set_W_mp(W_mp_blocks);
691
692 // df/dp -- scalar p
693 for (int i=0; i<num_p; i++) {
694 if (!outArgs.supports(OUT_ARG_DfDp, i).none()) {
695 Derivative dfdp = outArgs.get_DfDp(i);
696 if (dfdp.getMultiVector() != Teuchos::null) {
697 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dfdp_mp;
698 if (dfdp.getMultiVectorOrientation() == DERIV_MV_BY_COL)
699 dfdp_mp =
700 Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
701 mp_block_map, me->get_f_map(), mp_f_map, mp_comm,
702 View, *(dfdp.getMultiVector())));
703 else if (dfdp.getMultiVectorOrientation() == DERIV_TRANS_MV_BY_ROW)
704 dfdp_mp =
705 Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
706 mp_block_map, me->get_p_map(i), mp_p_map[i], mp_comm,
707 View, *(dfdp.getMultiVector())));
708 me_outargs.set_DfDp_mp(i,
709 MPDerivative(dfdp_mp,
710 dfdp.getMultiVectorOrientation()));
711 }
712 else if (dfdp.getLinearOp() != Teuchos::null) {
713 Teuchos::RCP<Stokhos::ProductEpetraOperator> dfdp_mp =
714 Teuchos::rcp_dynamic_cast<Stokhos::ProductEpetraOperator>(dfdp.getLinearOp(), true);
715 me_outargs.set_DfDp_mp(i, MPDerivative(dfdp_mp));
716 }
717 }
718 }
719
720 // dfdp -- block p. Here we only support DERIV_LINEAR_OP
721 for (int i=0; i<num_p_mp; i++) {
722 if (!outArgs.supports(OUT_ARG_DfDp, i+num_p).none()) {
723 Derivative dfdp = outArgs.get_DfDp(i+num_p);
724 if (dfdp.getLinearOp() != Teuchos::null) {
725 Teuchos::RCP<Stokhos::BlockDiagonalOperator> dfdp_op =
726 Teuchos::rcp_dynamic_cast<Stokhos::BlockDiagonalOperator>(dfdp.getLinearOp(), true);
727 Teuchos::RCP<Stokhos::ProductEpetraOperator> dfdp_op_mp =
728 dfdp_op->getMPOps();
729 int ii = mp_p_index_map[i];
730 if (me_outargs.supports(OUT_ARG_DfDp_mp,ii).supports(DERIV_LINEAR_OP))
731 me_outargs.set_DfDp_mp(ii, MPDerivative(dfdp_op_mp));
732 else {
733 Teuchos::RCP<Stokhos::ProductEpetraMultiVectorOperator> mp_mv_op =
734 Teuchos::rcp_dynamic_cast<Stokhos::ProductEpetraMultiVectorOperator>(dfdp_op_mp, true);
735 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dfdp_mp =
736 mp_mv_op->productMultiVector();
737 if (me_outargs.supports(OUT_ARG_DfDp_mp,ii).supports(DERIV_MV_BY_COL))
738 me_outargs.set_DfDp_mp(
739 ii, MPDerivative(dfdp_mp, DERIV_MV_BY_COL));
740 else
741 me_outargs.set_DfDp_mp(
742 ii, MPDerivative(dfdp_mp, DERIV_TRANS_MV_BY_ROW));
743 }
744 }
745 TEUCHOS_TEST_FOR_EXCEPTION(
746 dfdp.getLinearOp() == Teuchos::null && dfdp.isEmpty() == false,
747 std::logic_error,
748 "Error! Stokhos::MPModelEvaluator::evalModel: " <<
749 "Operator form of df/dp(" << i+num_p << ") is required!");
750 }
751 }
752
753 // Responses (g, dg/dx, dg/dp, ...)
754 for (int i=0; i<num_g_mp; i++) {
755 int ii = mp_g_index_map[i];
756
757 // g
758 Teuchos::RCP<Epetra_Vector> g = outArgs.get_g(i);
759 if (g != Teuchos::null) {
760 Teuchos::RCP<Stokhos::ProductEpetraVector> g_mp =
761 create_g_mp(ii, View, g.get());
762 me_outargs.set_g_mp(ii, g_mp);
763 }
764
765 // dg/dxdot
766 if (outArgs.supports(OUT_ARG_DgDx_dot, i).supports(DERIV_LINEAR_OP)) {
767 Derivative dgdx_dot = outArgs.get_DgDx_dot(i);
768 if (dgdx_dot.getLinearOp() != Teuchos::null) {
769 Teuchos::RCP<Stokhos::BlockDiagonalOperator> op =
770 Teuchos::rcp_dynamic_cast<Stokhos::BlockDiagonalOperator>(
771 dgdx_dot.getLinearOp(), true);
772 Teuchos::RCP<Stokhos::ProductEpetraOperator> mp_blocks =
773 op->getMPOps();
774 if (me_outargs.supports(OUT_ARG_DgDx, ii).supports(DERIV_LINEAR_OP))
775 me_outargs.set_DgDx_dot_mp(ii, mp_blocks);
776 else {
777 Teuchos::RCP<Stokhos::ProductEpetraMultiVectorOperator> mp_mv_op =
778 Teuchos::rcp_dynamic_cast<Stokhos::ProductEpetraMultiVectorOperator>(mp_blocks, true);
779 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdx_dot_mp =
780 mp_mv_op->productMultiVector();
781 if (me_outargs.supports(OUT_ARG_DgDx_dot_mp, ii).supports(DERIV_MV_BY_COL))
782 me_outargs.set_DgDx_dot_mp(ii, MPDerivative(dgdx_dot_mp,
783 DERIV_MV_BY_COL));
784 else
785 me_outargs.set_DgDx_dot_mp(ii, MPDerivative(dgdx_dot_mp,
786 DERIV_TRANS_MV_BY_ROW));
787 }
788 }
789 TEUCHOS_TEST_FOR_EXCEPTION(dgdx_dot.getLinearOp() == Teuchos::null &&
790 dgdx_dot.isEmpty() == false,
791 std::logic_error,
792 "Error! Stokhos::MPModelEvaluator::evalModel: " <<
793 "Operator form of dg/dxdot is required!");
794 }
795
796 // dg/dx
797 if (outArgs.supports(OUT_ARG_DgDx, i).supports(DERIV_LINEAR_OP)) {
798 Derivative dgdx = outArgs.get_DgDx(i);
799 if (dgdx.getLinearOp() != Teuchos::null) {
800 Teuchos::RCP<Stokhos::BlockDiagonalOperator> op =
801 Teuchos::rcp_dynamic_cast<Stokhos::BlockDiagonalOperator>(
802 dgdx.getLinearOp(), true);
803 Teuchos::RCP<Stokhos::ProductEpetraOperator> mp_blocks =
804 op->getMPOps();
805 if (me_outargs.supports(OUT_ARG_DgDx, ii).supports(DERIV_LINEAR_OP))
806 me_outargs.set_DgDx_mp(ii, mp_blocks);
807 else {
808 Teuchos::RCP<Stokhos::ProductEpetraMultiVectorOperator> mp_mv_op =
809 Teuchos::rcp_dynamic_cast<Stokhos::ProductEpetraMultiVectorOperator>(mp_blocks, true);
810 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdx_mp =
811 mp_mv_op->productMultiVector();
812 if (me_outargs.supports(OUT_ARG_DgDx_mp, ii).supports(DERIV_MV_BY_COL))
813 me_outargs.set_DgDx_mp(ii, MPDerivative(dgdx_mp,
814 DERIV_MV_BY_COL));
815 else
816 me_outargs.set_DgDx_mp(ii, MPDerivative(dgdx_mp,
817 DERIV_TRANS_MV_BY_ROW));
818 }
819 }
820 TEUCHOS_TEST_FOR_EXCEPTION(dgdx.getLinearOp() == Teuchos::null &&
821 dgdx.isEmpty() == false,
822 std::logic_error,
823 "Error! Stokhos::MPModelEvaluator::evalModel: " <<
824 "Operator form of dg/dxdot is required!");
825 }
826
827 // dg/dp -- scalar p
828 for (int j=0; j<num_p; j++) {
829 if (!outArgs.supports(OUT_ARG_DgDp, i, j).none()) {
830 Derivative dgdp = outArgs.get_DgDp(i,j);
831 if (dgdp.getMultiVector() != Teuchos::null) {
832 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdp_mp;
833 if (dgdp.getMultiVectorOrientation() == DERIV_MV_BY_COL)
834 dgdp_mp =
835 Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
836 mp_block_map, me->get_g_map(ii), mp_g_map[i],
837 mp_comm, View, *(dgdp.getMultiVector())));
838 else if (dgdp.getMultiVectorOrientation() == DERIV_TRANS_MV_BY_ROW)
839 dgdp_mp =
840 Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
841 mp_block_map, me->get_p_map(j), mp_p_map[j],
842 mp_comm, View, *(dgdp.getMultiVector())));
843 me_outargs.set_DgDp_mp(
844 ii, j, MPDerivative(dgdp_mp, dgdp.getMultiVectorOrientation()));
845 }
846 else if (dgdp.getLinearOp() != Teuchos::null) {
847 Teuchos::RCP<Stokhos::ProductEpetraOperator> dgdp_mp =
848 Teuchos::rcp_dynamic_cast<Stokhos::ProductEpetraOperator>(dgdp.getLinearOp(), true);
849 me_outargs.set_DgDp_mp(ii, j, MPDerivative(dgdp_mp));
850 }
851 }
852 }
853
854 // dgdp -- block p. Here we only support DERIV_LINEAR_OP
855 for (int j=0; j<num_p_mp; j++) {
856 if (!outArgs.supports(OUT_ARG_DgDp, i, j+num_p).none()) {
857 Derivative dgdp = outArgs.get_DgDp(i,j+num_p);
858 if (dgdp.getLinearOp() != Teuchos::null) {
859 Teuchos::RCP<Stokhos::BlockDiagonalOperator> dgdp_op =
860 Teuchos::rcp_dynamic_cast<Stokhos::BlockDiagonalOperator>(dgdp.getLinearOp(), true);
861 Teuchos::RCP<Stokhos::ProductEpetraOperator> dgdp_op_mp =
862 dgdp_op->getMPOps();
863 int jj = mp_p_index_map[j];
864 if (me_outargs.supports(OUT_ARG_DgDp_mp,ii,jj).supports(DERIV_LINEAR_OP))
865 me_outargs.set_DgDp_mp(ii, jj, MPDerivative(dgdp_op_mp));
866 else {
867 Teuchos::RCP<Stokhos::ProductEpetraMultiVectorOperator> mp_mv_op =
868 Teuchos::rcp_dynamic_cast<Stokhos::ProductEpetraMultiVectorOperator>(dgdp_op_mp, true);
869 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdp_mp =
870 mp_mv_op->productMultiVector();
871 if (me_outargs.supports(OUT_ARG_DgDp_mp,ii,jj).supports(DERIV_MV_BY_COL))
872 me_outargs.set_DgDp_mp(
873 ii, jj, MPDerivative(dgdp_mp, DERIV_MV_BY_COL));
874 else
875 me_outargs.set_DgDp_mp(
876 ii, jj, MPDerivative(dgdp_mp, DERIV_TRANS_MV_BY_ROW));
877 }
878 }
879 TEUCHOS_TEST_FOR_EXCEPTION(
880 dgdp.getLinearOp() == Teuchos::null && dgdp.isEmpty() == false,
881 std::logic_error,
882 "Error! Stokhos::MPModelEvaluator::evalModel: " <<
883 "Operator form of dg/dp(" << i << "," << j+num_p << ") is required!");
884 }
885 }
886
887 }
888
889 // Compute the functions
890 me->evalModel(me_inargs, me_outargs);
891
892 // Copy block MP components for W
893 if (W_out != Teuchos::null && !eval_prec) {
894 Teuchos::RCP<Epetra_Operator> W;
895 if (W_out != Teuchos::null)
896 W = W_out;
897 else
898 W = my_W;
899 Teuchos::RCP<Stokhos::BlockDiagonalOperator> W_mp =
900 Teuchos::rcp_dynamic_cast<Stokhos::BlockDiagonalOperator>(W, true);
901 W_mp->setupOperator(W_mp_blocks);
902
903 if (WPrec_out != Teuchos::null) {
904 Teuchos::RCP<Stokhos::MPPreconditioner> W_prec =
905 Teuchos::rcp_dynamic_cast<Stokhos::MPPreconditioner>(WPrec_out, true);
906 W_prec->setupPreconditioner(W_mp, *my_x);
907 }
908 }
909}
910
911void
913 const Stokhos::ProductEpetraVector& x_mp_in)
914{
915 *mp_x_init = x_mp_in;
916}
917
918void
920 const Stokhos::ProductEpetraVector& x_dot_mp_in)
921{
922 *mp_x_dot_init = x_dot_mp_in;
923}
924
925Teuchos::RCP<const Stokhos::ProductEpetraVector>
927{
928 return mp_x_init;
929}
930
931Teuchos::RCP<const Stokhos::ProductEpetraVector>
933{
934 return mp_x_dot_init;
935}
936
937void
939 int i, const Stokhos::ProductEpetraVector& p_mp_in)
940{
941 Teuchos::Array<int>::iterator it = std::find(mp_p_index_map.begin(),
942 mp_p_index_map.end(),
943 i);
944 TEUCHOS_TEST_FOR_EXCEPTION(it == mp_p_index_map.end(), std::logic_error,
945 "Error! Invalid p map index " << i);
946 int ii = it - mp_p_index_map.begin();
947 *mp_p_init[ii] = p_mp_in;
948}
949
950Teuchos::RCP<const Stokhos::ProductEpetraVector>
952{
953 Teuchos::Array<int>::const_iterator it = std::find(mp_p_index_map.begin(),
954 mp_p_index_map.end(),
955 l);
956 TEUCHOS_TEST_FOR_EXCEPTION(it == mp_p_index_map.end(), std::logic_error,
957 "Error! Invalid p map index " << l);
958 int ll = it - mp_p_index_map.begin();
959 return mp_p_init[ll];
960}
961
962Teuchos::Array<int>
964{
965 return mp_p_index_map;
966}
967
968Teuchos::Array<int>
970{
971 return mp_g_index_map;
972}
973
974Teuchos::Array< Teuchos::RCP<const Epetra_Map> >
976{
977 Teuchos::Array< Teuchos::RCP<const Epetra_Map> > base_maps(num_g);
978 for (int i=0; i<num_g; i++)
979 base_maps[i] = me->get_g_map(i);
980 return base_maps;
981 }
982
983Teuchos::RCP<Stokhos::ProductEpetraVector>
985 const Epetra_Vector* v) const
986{
987 Teuchos::RCP<Stokhos::ProductEpetraVector> mp_x;
988 if (v == NULL)
989 mp_x = Teuchos::rcp(new Stokhos::ProductEpetraVector(
990 mp_block_map, x_map, mp_x_map, mp_comm));
991 else
992 mp_x = Teuchos::rcp(new Stokhos::ProductEpetraVector(
993 mp_block_map, x_map, mp_x_map, mp_comm,
994 CV, *v));
995 return mp_x;
996}
997
998Teuchos::RCP<Stokhos::ProductEpetraMultiVector>
1000 const Epetra_MultiVector* v) const
1001{
1002 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_x;
1003 if (v == NULL)
1004 mp_x = Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
1005 mp_block_map, x_map, mp_x_map, mp_comm,
1006 num_vecs));
1007 else
1008 mp_x = Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
1009 mp_block_map, x_map, mp_x_map, mp_comm,
1010 CV, *v));
1011 return mp_x;
1012}
1013
1014Teuchos::RCP<Stokhos::ProductEpetraVector>
1016 const Epetra_Vector* v) const
1017{
1018 Teuchos::RCP<Stokhos::ProductEpetraVector> mp_p;
1019 Teuchos::Array<int>::const_iterator it = std::find(mp_p_index_map.begin(),
1020 mp_p_index_map.end(),
1021 l);
1022 TEUCHOS_TEST_FOR_EXCEPTION(it == mp_p_index_map.end(), std::logic_error,
1023 "Error! Invalid p map index " << l);
1024 int ll = it - mp_p_index_map.begin();
1025 if (v == NULL)
1026 mp_p = Teuchos::rcp(new Stokhos::ProductEpetraVector(
1027 mp_block_map, me->get_p_map(l),
1028 mp_p_map[ll], mp_comm));
1029 else
1030 mp_p = Teuchos::rcp(new Stokhos::ProductEpetraVector(
1031 mp_block_map, me->get_p_map(l),
1032 mp_p_map[ll], mp_comm, CV, *v));
1033 return mp_p;
1034}
1035
1036Teuchos::RCP<Stokhos::ProductEpetraMultiVector>
1039 const Epetra_MultiVector* v) const
1040{
1041 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_p;
1042 Teuchos::Array<int>::const_iterator it = std::find(mp_p_index_map.begin(),
1043 mp_p_index_map.end(),
1044 l);
1045 TEUCHOS_TEST_FOR_EXCEPTION(it == mp_p_index_map.end(), std::logic_error,
1046 "Error! Invalid p map index " << l);
1047 int ll = it - mp_p_index_map.begin();
1048 if (v == NULL)
1049 mp_p = Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
1050 mp_block_map, me->get_p_map(l),
1051 mp_p_map[ll], mp_comm, num_vecs));
1052 else
1053 mp_p = Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
1054 mp_block_map, me->get_p_map(l),
1055 mp_p_map[ll], mp_comm, CV, *v));
1056 return mp_p;
1057}
1058
1059Teuchos::RCP<Stokhos::ProductEpetraVector>
1061 const Epetra_Vector* v) const
1062{
1063 Teuchos::RCP<Stokhos::ProductEpetraVector> mp_f;
1064 if (v == NULL)
1065 mp_f = Teuchos::rcp(new Stokhos::ProductEpetraVector(
1066 mp_block_map, f_map, mp_f_map, mp_comm));
1067 else
1068 mp_f = Teuchos::rcp(new Stokhos::ProductEpetraVector(
1069 mp_block_map, f_map, mp_f_map, mp_comm,
1070 CV, *v));
1071 return mp_f;
1072}
1073
1074Teuchos::RCP<Stokhos::ProductEpetraMultiVector>
1076 int num_vecs,
1078 const Epetra_MultiVector* v) const
1079{
1080 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_f;
1081 if (v == NULL)
1082 mp_f = Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
1083 mp_block_map, f_map, mp_f_map, mp_comm,
1084 num_vecs));
1085 else
1086 mp_f = Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
1087 mp_block_map, f_map, mp_f_map, mp_comm,
1088 CV, *v));
1089 return mp_f;
1090}
1091
1092Teuchos::RCP<Stokhos::ProductEpetraVector>
1094 const Epetra_Vector* v) const
1095{
1096 Teuchos::RCP<Stokhos::ProductEpetraVector> mp_g;
1097 Teuchos::Array<int>::const_iterator it = std::find(mp_g_index_map.begin(),
1098 mp_g_index_map.end(),
1099 l);
1100 TEUCHOS_TEST_FOR_EXCEPTION(it == mp_g_index_map.end(), std::logic_error,
1101 "Error! Invalid g map index " << l);
1102 int ll = it - mp_g_index_map.begin();
1103 if (v == NULL)
1104 mp_g = Teuchos::rcp(new Stokhos::ProductEpetraVector(
1105 mp_block_map,
1106 me->get_g_map(l),
1107 mp_g_map[ll], mp_comm));
1108 else
1109 mp_g = Teuchos::rcp(new Stokhos::ProductEpetraVector(
1110 mp_block_map,
1111 me->get_g_map(l),
1112 mp_g_map[ll], mp_comm, CV, *v));
1113 return mp_g;
1114}
1115
1116Teuchos::RCP<Stokhos::ProductEpetraMultiVector>
1119 const Epetra_MultiVector* v) const
1120{
1121 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_g;
1122 Teuchos::Array<int>::const_iterator it = std::find(mp_g_index_map.begin(),
1123 mp_g_index_map.end(),
1124 l);
1125 TEUCHOS_TEST_FOR_EXCEPTION(it == mp_g_index_map.end(), std::logic_error,
1126 "Error! Invalid g map index " << l);
1127 int ll = it - mp_g_index_map.begin();
1128 if (v == NULL)
1129 mp_g = Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
1130 mp_block_map,
1131 me->get_g_map(l),
1132 mp_g_map[ll], mp_comm, num_vecs));
1133 else
1134 mp_g = Teuchos::rcp(new Stokhos::ProductEpetraMultiVector(
1135 mp_block_map,
1136 me->get_g_map(l),
1137 mp_g_map[ll], mp_comm, CV, *v));
1138 return mp_g;
1139}
Epetra_DataAccess
An Epetra operator representing the block stochastic Galerkin operator.
Teuchos::RCP< EpetraExt::ModelEvaluator::Preconditioner > create_WPrec() const
Create preconditioner operator.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > get_g_mp_base_maps() const
Get base maps of MP responses.
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
OutArgs createOutArgs() const
Create OutArgs.
Teuchos::RCP< Stokhos::ProductEpetraMultiVector > create_g_mv_mp(int l, int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create multi-point multi-vector using g map.
void set_p_mp_init(int i, const Stokhos::ProductEpetraVector &p_mp_in)
Set initial multi-point parameter.
Teuchos::RCP< const Epetra_Map > x_map
Underlying unknown map.
Teuchos::RCP< const Epetra_Vector > get_x_dot_init() const
Teuchos::RCP< Stokhos::ProductEpetraVector > create_x_mp(Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create multi-point vector using x map and owned mp map.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
Teuchos::RCP< Stokhos::ProductEpetraVector > create_f_mp(Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create multi-point vector using f map.
bool supports_x
Whether we support x (and thus f and W)
Teuchos::Array< int > mp_g_index_map
Index map between block-g and g_mp maps.
Teuchos::RCP< Stokhos::ProductEpetraMultiVector > create_x_mv_mp(int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create multi-point vector using x map.
Teuchos::RCP< const Epetra_Map > f_map
Underlying residual map.
void set_x_mp_init(const Stokhos::ProductEpetraVector &x_mp_in)
Set initial multi-point solution.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
Teuchos::Array< Teuchos::RCP< Teuchos::Array< std::string > > > mp_p_names
MP coefficient parameter names.
int num_g_mp
Number of multi-point response vectors.
Teuchos::RCP< Stokhos::ProductEpetraOperator > W_mp_blocks
W multi-point components.
MPModelEvaluator(const Teuchos::RCP< EpetraExt::ModelEvaluator > &me, const Teuchos::RCP< const EpetraExt::MultiComm > &mp_comm, const Teuchos::RCP< const Epetra_Map > &mp_block_map, const Teuchos::RCP< Teuchos::ParameterList > &params)
Teuchos::RCP< const Epetra_Map > mp_block_map
Map for layout of parallel MP blocks.
Teuchos::RCP< const Epetra_Map > get_g_map(int l) const
Return response map.
Teuchos::RCP< Stokhos::ProductEpetraVector > mp_x_dot_init
Teuchos::RCP< Epetra_Operator > create_DgDp_op(int j, int i) const
Create MP operator representing dg/dp.
Teuchos::RCP< Epetra_Operator > create_DfDp_op(int i) const
Create MP operator representing df/dp.
Teuchos::Array< int > get_p_mp_map_indices() const
Get indices of MP parameters.
Teuchos::RCP< Epetra_Operator > create_DgDx_op(int j) const
Create MP operator representing dg/dx.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > mp_g_map
Block MP response map.
Teuchos::RCP< EpetraExt::ModelEvaluator > me
Underlying model evaluator.
Teuchos::Array< Teuchos::RCP< ProductEpetraVector > > mp_p_init
MP initial p.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
Teuchos::RCP< Stokhos::ProductEpetraMultiVector > create_f_mv_mp(int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create multi-point multi-vector using f map.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
void set_x_dot_mp_init(const Stokhos::ProductEpetraVector &x_dot_mp_in)
Teuchos::RCP< Epetra_Operator > create_DgDx_dot_op(int j) const
Create MP operator representing dg/dxdot.
Teuchos::RCP< Stokhos::ProductEpetraMultiVector > create_p_mv_mp(int l, int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create multi-point vector using p map.
Teuchos::RCP< const Stokhos::ProductEpetraVector > get_x_dot_mp_init() const
Teuchos::RCP< const Stokhos::ProductEpetraVector > get_x_mp_init() const
Return initial SG x.
unsigned int num_mp_blocks
Number of blocks.
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
Teuchos::RCP< const EpetraExt::MultiComm > mp_comm
Parallel MP communicator.
Teuchos::RCP< const Epetra_Map > mp_f_map
Block MP residual map.
Teuchos::RCP< Stokhos::ProductEpetraVector > create_p_mp(int l, Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create multi-point vector using p map.
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
Teuchos::Array< int > get_g_mp_map_indices() const
Get indices of MP responses.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > mp_p_map
Block MP parameter map.
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
Teuchos::RCP< const Stokhos::ProductEpetraVector > get_p_mp_init(int l) const
Return initial SG parameters.
int num_p_mp
Number of multi-point parameter vectors.
InArgs createInArgs() const
Create InArgs.
Teuchos::RCP< Epetra_Vector > my_x
x pointer for evaluating preconditioner
int num_g
Number of response vectors of underlying model evaluator.
Teuchos::RCP< Stokhos::ProductEpetraVector > mp_x_init
MP initial x.
Teuchos::RCP< const Epetra_Map > mp_x_map
Block MP unknown map.
Teuchos::Array< int > mp_p_index_map
Index map between block-p and p_mp maps.
int num_p
Number of parameter vectors of underlying model evaluator.
Teuchos::RCP< Stokhos::ProductEpetraVector > create_g_mp(int l, Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create multi-point vector using g map.
Factory for generating stochastic Galerkin preconditioners.
A container class for products of Epetra_Vector's.
A container class storing products of Epetra_MultiVector's.
A container class for products of Epetra_Vector's.
A container class for products of Epetra_Vector's.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)