Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_SGModelEvaluator.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#include "Epetra_LocalMap.h"
53#include "Epetra_Export.h"
54#include "Epetra_Import.h"
55
57 const Teuchos::RCP<EpetraExt::ModelEvaluator>& me_,
58 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sg_basis_,
59 const Teuchos::RCP<const Stokhos::Quadrature<int,double> >& sg_quad_,
60 const Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> >& sg_exp_,
61 const Teuchos::RCP<const Stokhos::ParallelData>& sg_parallel_data_,
62 const Teuchos::RCP<Teuchos::ParameterList>& params_,
63 bool scaleOP_)
64 : me(me_),
65 sg_basis(sg_basis_),
66 sg_quad(sg_quad_),
67 sg_exp(sg_exp_),
68 params(params_),
69 num_sg_blocks(sg_basis->size()),
70 num_W_blocks(sg_basis->size()),
71 num_p_blocks(sg_basis->size()),
72 supports_x(false),
73 x_map(me->get_x_map()),
74 f_map(me->get_f_map()),
75 sg_parallel_data(sg_parallel_data_),
76 sg_comm(sg_parallel_data->getMultiComm()),
77 epetraCijk(sg_parallel_data->getEpetraCijk()),
78 serialCijk(),
79 sg_x_map(),
80 sg_overlapped_x_map(),
81 sg_f_map(),
82 sg_overlapped_f_map(),
83 sg_overlapped_x_importer(),
84 sg_overlapped_f_exporter(),
85 num_p(0),
86 num_p_sg(0),
87 sg_p_index_map(),
88 sg_p_map(),
89 sg_p_names(),
90 num_g(0),
91 num_g_sg(0),
92 sg_g_index_map(),
93 sg_g_map(),
94 x_dot_sg_blocks(),
95 x_sg_blocks(),
96 f_sg_blocks(),
97 W_sg_blocks(),
98 sg_x_init(),
99 sg_p_init(),
100 eval_W_with_f(false),
101 scaleOP(scaleOP_)
102{
103 if (x_map != Teuchos::null)
104 supports_x = true;
105
107 Teuchos::rcp(new Epetra_LocalMap(
108 static_cast<int>(num_sg_blocks), 0, *(sg_parallel_data->getStochasticComm())));
110 if (epetraCijk != Teuchos::null)
111 stoch_row_map = epetraCijk->getStochasticRowMap();
112
113 if (supports_x) {
114
115 // Create block SG x and f maps
116 sg_x_map =
117 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
120 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
122
123 sg_f_map =
124 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
127 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
129
130 // Create importer/exporter from/to overlapped distribution
132 Teuchos::rcp(new Epetra_Import(*sg_overlapped_x_map, *sg_x_map));
134 Teuchos::rcp(new Epetra_Export(*sg_overlapped_f_map, *sg_f_map));
135
136 // Create vector blocks
140
141 // Create default sg_x_init
143 if (sg_x_init->myGID(0))
144 (*sg_x_init)[0] = *(me->get_x_init());
145
146 // Preconditioner needs an x
147 my_x = Teuchos::rcp(new Epetra_Vector(*sg_x_map));
148
149 // Determine W expansion type
150 std::string W_expansion_type =
151 params->get("Jacobian Expansion Type", "Full");
152 if (W_expansion_type == "Linear")
153 num_W_blocks = sg_basis->dimension() + 1;
154 else
156 Teuchos::RCP<Epetra_BlockMap> W_overlap_map =
157 Teuchos::rcp(new Epetra_LocalMap(
158 static_cast<int>(num_W_blocks), 0,
159 *(sg_parallel_data->getStochasticComm())));
161 Teuchos::rcp(new Stokhos::EpetraOperatorOrthogPoly(
162 sg_basis, W_overlap_map, x_map, f_map, sg_f_map, sg_comm));
163 for (unsigned int i=0; i<num_W_blocks; i++)
164 W_sg_blocks->setCoeffPtr(i, me->create_W());
165
166 eval_W_with_f = params->get("Evaluate W with F", false);
167
168 Teuchos::RCP<Teuchos::ParameterList> sgPrecParams =
169 Teuchos::rcp(&(params->sublist("SG Preconditioner")), false);
171 Teuchos::rcp(new Stokhos::SGPreconditionerFactory (sgPrecParams));
172 }
173
174 // Parameters -- The idea here is to add new parameter vectors
175 // for the stochastic Galerkin components of the parameters
176
177 InArgs me_inargs = me->createInArgs();
178 OutArgs me_outargs = me->createOutArgs();
179 num_p = me_inargs.Np();
180
181 // Get the p_sg's supported and build index map
182 for (int i=0; i<num_p; i++) {
183 if (me_inargs.supports(IN_ARG_p_sg, i))
184 sg_p_index_map.push_back(i);
185 }
186 num_p_sg = sg_p_index_map.size();
187
188 sg_p_map.resize(num_p_sg);
189 sg_p_names.resize(num_p_sg);
190 sg_p_init.resize(num_p_sg);
191
192 // Determine parameter expansion type
193 std::string p_expansion_type =
194 params->get("Parameter Expansion Type", "Full");
195 if (p_expansion_type == "Linear")
196 num_p_blocks = sg_basis->dimension() + 1;
197 else
199
200 // Create parameter maps, names, and initial values
202 Teuchos::rcp(new Epetra_LocalMap(
203 static_cast<int>(num_p_blocks), 0,
204 *(sg_parallel_data->getStochasticComm())));
205 for (int i=0; i<num_p_sg; i++) {
206 Teuchos::RCP<const Epetra_Map> p_map = me->get_p_map(sg_p_index_map[i]);
207 sg_p_map[i] =
208 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
209 *p_map, *overlapped_stoch_p_map, *sg_comm));
210
211 Teuchos::RCP<const Teuchos::Array<std::string> > p_names =
212 me->get_p_names(sg_p_index_map[i]);
213 if (p_names != Teuchos::null) {
214 sg_p_names[i] =
215 Teuchos::rcp(new Teuchos::Array<std::string>(num_sg_blocks*(p_names->size())));
216 for (int j=0; j<p_names->size(); j++) {
217 std::stringstream ss;
218 ss << (*p_names)[j] << " -- SG Coefficient " << i;
219 (*sg_p_names[i])[j] = ss.str();
220 }
221 }
222
223 // Create default sg_p_init with an expansion equal to the mean parameter
225 sg_p_init[i]->init(0.0);
226 (*sg_p_init[i])[0] = *(me->get_p_init(sg_p_index_map[i]));
227 }
228
229 // Responses -- The idea here is to add new parameter vectors
230 // for the stochastic Galerkin components of the respones
231 num_g = me_outargs.Ng();
232
233 // Get the g_sg's supported and build index map
234 for (int i=0; i<num_g; i++) {
235 if (me_outargs.supports(OUT_ARG_g_sg, i))
236 sg_g_index_map.push_back(i);
237 }
238 num_g_sg = sg_g_index_map.size();
239
240 sg_g_map.resize(num_g_sg);
241
242 // Create response maps
243 for (int i=0; i<num_g_sg; i++) {
244 Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(sg_g_index_map[i]);
245 sg_g_map[i] =
246 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
248 }
249
250 // We don't support parallel for dgdx yet, so build a new EpetraCijk
251 if (supports_x) {
252 serialCijk =
253 Teuchos::rcp(new Stokhos::EpetraSparse3Tensor(sg_basis,
254 epetraCijk->getCijk(),
255 sg_comm,
257 }
258
259}
260
261// Overridden from EpetraExt::ModelEvaluator
262
263Teuchos::RCP<const Epetra_Map>
265{
266 return sg_x_map;
267}
268
269Teuchos::RCP<const Epetra_Map>
271{
272 return sg_f_map;
273}
274
275Teuchos::RCP<const Epetra_Map>
277{
278 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_sg, std::logic_error,
279 "Error! Invalid p map index " << l);
280 if (l < num_p)
281 return me->get_p_map(l);
282 else
283 return sg_p_map[l-num_p];
284
285 return Teuchos::null;
286}
287
288Teuchos::RCP<const Epetra_Map>
290{
291 TEUCHOS_TEST_FOR_EXCEPTION(j < 0 || j >= num_g_sg, std::logic_error,
292 "Error! Invalid g map index " << j);
293 return sg_g_map[j];
294}
295
296Teuchos::RCP<const Teuchos::Array<std::string> >
298{
299 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_sg, std::logic_error,
300 "Error! Invalid p map index " << l);
301 if (l < num_p)
302 return me->get_p_names(l);
303 else
304 return sg_p_names[l-num_p];
305
306 return Teuchos::null;
307}
308
309Teuchos::RCP<const Epetra_Vector>
311{
312 return sg_x_init->getBlockVector();
313}
314
315Teuchos::RCP<const Epetra_Vector>
317{
318 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_sg, std::logic_error,
319 "Error! Invalid p map index " << l);
320 if (l < num_p)
321 return me->get_p_init(l);
322 else
323 return sg_p_init[l-num_p]->getBlockVector();
324
325 return Teuchos::null;
326}
327
328Teuchos::RCP<Epetra_Operator>
330{
331 if (supports_x) {
332 Teuchos::RCP<Teuchos::ParameterList> sgOpParams =
333 Teuchos::rcp(&(params->sublist("SG Operator")), false);
334 Teuchos::RCP<Epetra_CrsMatrix> W_crs;
335 if (sgOpParams->get("Operator Method", "Matrix Free") ==
336 "Fully Assembled") {
337 W_crs = Teuchos::rcp_dynamic_cast<Epetra_CrsMatrix>(me->create_W(), true);
338 Teuchos::RCP<const Epetra_CrsGraph> W_graph =
339 Teuchos::rcp(&(W_crs->Graph()), false);
340 sgOpParams->set< Teuchos::RCP<const Epetra_CrsGraph> >("Base Graph",
341 W_graph);
342 }
343 Stokhos::SGOperatorFactory sg_op_factory(sgOpParams);
344 my_W = sg_op_factory.build(sg_comm, sg_basis, epetraCijk, x_map, f_map,
345 sg_x_map, sg_f_map);
346 my_W->setupOperator(W_sg_blocks);
347
348 return my_W;
349 }
350
351 return Teuchos::null;
352}
353
354Teuchos::RCP<EpetraExt::ModelEvaluator::Preconditioner>
356{
357 if (supports_x) {
358
359 Teuchos::RCP<Epetra_Operator> precOp =
360 sg_prec_factory->build(sg_comm, sg_basis, epetraCijk, x_map, sg_x_map);
361 return Teuchos::rcp(new EpetraExt::ModelEvaluator::Preconditioner(precOp,
362 true));
363 }
364
365 return Teuchos::null;
366}
367
368Teuchos::RCP<Epetra_Operator>
370{
371 TEUCHOS_TEST_FOR_EXCEPTION(j < 0 || j >= num_g_sg || !supports_x,
372 std::logic_error,
373 "Error: dg/dx index " << j << " is not supported!");
374
375 int jj = sg_g_index_map[j];
376 Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(jj);
377 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > sg_blocks;
378 OutArgs me_outargs = me->createOutArgs();
379 DerivativeSupport ds = me_outargs.supports(OUT_ARG_DgDx_sg, jj);
380 if (ds.supports(DERIV_LINEAR_OP)) {
381 sg_blocks =
382 Teuchos::rcp(new Stokhos::EpetraOperatorOrthogPoly(
383 sg_basis, overlapped_stoch_row_map, x_map,
384 g_map, sg_g_map[j], sg_comm));
385 for (unsigned int i=0; i<num_sg_blocks; i++)
386 sg_blocks->setCoeffPtr(i, me->create_DgDx_op(jj));
387 }
388 else if (ds.supports(DERIV_MV_BY_COL)) {
389 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_mv_blocks =
391 sg_basis, overlapped_stoch_row_map, g_map, sg_g_map[j],
392 sg_comm, x_map->NumMyElements()));
393 sg_blocks =
395 sg_mv_blocks, false));
396 }
397 else if (ds.supports(DERIV_TRANS_MV_BY_ROW)) {
398 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_mv_blocks =
400 sg_basis, overlapped_stoch_row_map, x_map, sg_x_map,
401 sg_comm, g_map->NumMyElements()));
402 sg_blocks =
404 sg_mv_blocks, true));
405 }
406 else
407 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
408 "Error! me_outargs.supports(OUT_ARG_DgDx_sg, " << jj
409 << ").none() is true!");
410
411 Teuchos::RCP<Teuchos::ParameterList> pl =
412 Teuchos::rcp(new Teuchos::ParameterList);
413 Teuchos::RCP<Stokhos::SGOperator> dgdx_sg =
414 Teuchos::rcp(new Stokhos::MatrixFreeOperator(
415 sg_comm, sg_basis, serialCijk, x_map, g_map,
416 sg_x_map, sg_g_map[j], pl));
417 dgdx_sg->setupOperator(sg_blocks);
418 return dgdx_sg;
419}
420
421Teuchos::RCP<Epetra_Operator>
423{
424 TEUCHOS_TEST_FOR_EXCEPTION(j < 0 || j >= num_g_sg || !supports_x,
425 std::logic_error,
426 "Error: dg/dx_dot index " << j << " is not supported!");
427
428 int jj = sg_g_index_map[j];
429 Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(jj);
430 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > sg_blocks;
431 OutArgs me_outargs = me->createOutArgs();
432 DerivativeSupport ds = me_outargs.supports(OUT_ARG_DgDx_dot_sg, jj);
433 if (ds.supports(DERIV_LINEAR_OP)) {
434 sg_blocks =
435 Teuchos::rcp(new Stokhos::EpetraOperatorOrthogPoly(
436 sg_basis, overlapped_stoch_row_map, x_map,
437 g_map, sg_g_map[j], sg_comm));
438 for (unsigned int i=0; i<num_sg_blocks; i++)
439 sg_blocks->setCoeffPtr(i, me->create_DgDx_dot_op(jj));
440 }
441 else if (ds.supports(DERIV_MV_BY_COL)) {
442 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_mv_blocks =
444 sg_basis, overlapped_stoch_row_map, g_map, sg_g_map[j],
445 sg_comm, x_map->NumMyElements()));
446 sg_blocks =
448 sg_mv_blocks, false));
449 }
450 else if (ds.supports(DERIV_TRANS_MV_BY_ROW)) {
451 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_mv_blocks =
453 sg_basis, overlapped_stoch_row_map, x_map, sg_x_map,
454 sg_comm, g_map->NumMyElements()));
455 sg_blocks =
457 sg_mv_blocks, true));
458 }
459 else
460 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
461 "Error! me_outargs.supports(OUT_ARG_DgDx_dot_sg, "
462 << jj << ").none() is true!");
463
464 Teuchos::RCP<Teuchos::ParameterList> pl =
465 Teuchos::rcp(new Teuchos::ParameterList);
466 Teuchos::RCP<Stokhos::SGOperator> dgdx_dot_sg =
467 Teuchos::rcp(new Stokhos::MatrixFreeOperator(
468 sg_comm, sg_basis, serialCijk, x_map, g_map,
469 sg_x_map, sg_g_map[j], pl));
470 dgdx_dot_sg->setupOperator(sg_blocks);
471 return dgdx_dot_sg;
472}
473
474Teuchos::RCP<Epetra_Operator>
476{
477 TEUCHOS_TEST_FOR_EXCEPTION(
478 j < 0 || j >= num_g_sg || i < 0 || i >= num_p+num_p_sg,
479 std::logic_error,
480 "Error: dg/dp index " << j << "," << i << " is not supported!");
481
482 int jj = sg_g_index_map[j];
483 Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(jj);
484 OutArgs me_outargs = me->createOutArgs();
485 if (i < num_p) {
486 if (me_outargs.supports(OUT_ARG_DgDp_sg,jj,i).supports(DERIV_LINEAR_OP)) {
487 Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly> sg_blocks =
488 Teuchos::rcp(new Stokhos::EpetraOperatorOrthogPoly(
489 sg_basis, overlapped_stoch_row_map,
490 me->get_p_map(i), me->get_g_map(jj),
491 sg_g_map[j], sg_comm));
492 for (unsigned int l=0; l<num_sg_blocks; l++)
493 sg_blocks->setCoeffPtr(l, me->create_DgDp_op(jj,i));
494 return sg_blocks;
495 }
496 else
497 TEUCHOS_TEST_FOR_EXCEPTION(
498 true, std::logic_error,
499 "Error: Underlying model evaluator must support DERIV_LINER_OP " <<
500 "to create operator for dg/dp index " << j << "," << i << "!");
501 }
502 else {
503 int ii = sg_p_index_map[i-num_p];
504 Teuchos::RCP<const Epetra_Map> p_map = me->get_p_map(ii);
505 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > sg_blocks;
506 DerivativeSupport ds = me_outargs.supports(OUT_ARG_DgDp_sg,jj,ii);
507 if (ds.supports(DERIV_LINEAR_OP)) {
508 sg_blocks =
509 Teuchos::rcp(new Stokhos::EpetraOperatorOrthogPoly(
510 sg_basis, overlapped_stoch_row_map,
511 p_map, g_map, sg_g_map[j], sg_comm));
512 for (unsigned int l=0; l<num_sg_blocks; l++)
513 sg_blocks->setCoeffPtr(l, me->create_DgDp_op(jj,ii));
514 }
515 else if (ds.supports(DERIV_MV_BY_COL)) {
516 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_mv_blocks =
518 sg_basis, overlapped_stoch_row_map, g_map, sg_g_map[j],
519 sg_comm, p_map->NumMyElements()));
520 sg_blocks =
522 sg_mv_blocks, false));
523 }
524 else if (ds.supports(DERIV_TRANS_MV_BY_ROW)) {
525 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_mv_blocks =
527 sg_basis, overlapped_stoch_row_map, p_map,
528 sg_p_map[i-num_p],
529 sg_comm, g_map->NumMyElements()));
530 sg_blocks =
532 sg_mv_blocks, true));
533 }
534 else
535 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
536 "Error! me_outargs.supports(OUT_ARG_DgDp_sg, " << jj
537 << "," << ii << ").none() is true!");
538
539 Teuchos::RCP<Teuchos::ParameterList> pl =
540 Teuchos::rcp(new Teuchos::ParameterList);
541 Teuchos::RCP<Stokhos::SGOperator> dgdp_sg =
542 Teuchos::rcp(new Stokhos::MatrixFreeOperator(
543 sg_comm, sg_basis, serialCijk,
544 p_map, g_map, sg_p_map[i-num_p], sg_g_map[j], pl));
545 dgdp_sg->setupOperator(sg_blocks);
546 return dgdp_sg;
547 }
548
549 return Teuchos::null;
550}
551
552Teuchos::RCP<Epetra_Operator>
554{
555 TEUCHOS_TEST_FOR_EXCEPTION(i < 0 || i >= num_p+num_p_sg,
556 std::logic_error,
557 "Error: df/dp index " << i << " is not supported!");
558
559 OutArgs me_outargs = me->createOutArgs();
560 if (i < num_p) {
561 if (me_outargs.supports(OUT_ARG_DfDp_sg,i).supports(DERIV_LINEAR_OP)) {
562 Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly> sg_blocks =
563 Teuchos::rcp(new Stokhos::EpetraOperatorOrthogPoly(
564 sg_basis, overlapped_stoch_row_map,
565 me->get_p_map(i), f_map,
566 sg_overlapped_f_map, sg_comm));
567 for (unsigned int l=0; l<num_sg_blocks; l++)
568 sg_blocks->setCoeffPtr(l, me->create_DfDp_op(i));
569 return sg_blocks;
570 }
571 else
572 TEUCHOS_TEST_FOR_EXCEPTION(
573 true, std::logic_error,
574 "Error: Underlying model evaluator must support DERIV_LINER_OP " <<
575 "to create operator for df/dp index " << i << "!");
576 }
577 else {
578 int ii = sg_p_index_map[i-num_p];
579 Teuchos::RCP<const Epetra_Map> p_map = me->get_p_map(ii);
580 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > sg_blocks;
581 DerivativeSupport ds = me_outargs.supports(OUT_ARG_DfDp_sg,ii);
582 if (ds.supports(DERIV_LINEAR_OP)) {
583 sg_blocks =
584 Teuchos::rcp(new Stokhos::EpetraOperatorOrthogPoly(
585 sg_basis, overlapped_stoch_row_map,
586 p_map, f_map, sg_overlapped_f_map, sg_comm));
587 for (unsigned int l=0; l<num_sg_blocks; l++)
588 sg_blocks->setCoeffPtr(l, me->create_DfDp_op(ii));
589 }
590 else if (ds.supports(DERIV_MV_BY_COL)) {
591 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_mv_blocks =
593 sg_basis, overlapped_stoch_row_map, f_map, sg_f_map,
594 sg_comm, p_map->NumMyElements()));
595 sg_blocks =
597 sg_mv_blocks, false));
598 }
599 else if (ds.supports(DERIV_TRANS_MV_BY_ROW)) {
600 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_mv_blocks =
602 sg_basis, overlapped_stoch_row_map, p_map,
603 sg_p_map[i-num_p],
604 sg_comm, f_map->NumMyElements()));
605 sg_blocks =
607 sg_mv_blocks, true));
608 }
609 else
610 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
611 "Error! me_outargs.supports(OUT_ARG_DfDp_sg, " << ii
612 << ").none() is true!");
613
614 Teuchos::RCP<Teuchos::ParameterList> pl =
615 Teuchos::rcp(new Teuchos::ParameterList);
616 Teuchos::RCP<Stokhos::SGOperator> dfdp_sg =
617 Teuchos::rcp(new Stokhos::MatrixFreeOperator(
618 sg_comm, sg_basis, epetraCijk,
619 p_map, f_map, sg_p_map[i-num_p], sg_f_map, pl));
620 dfdp_sg->setupOperator(sg_blocks);
621 return dfdp_sg;
622 }
623
624 return Teuchos::null;
625}
626
627EpetraExt::ModelEvaluator::InArgs
629{
630 InArgsSetup inArgs;
631 InArgs me_inargs = me->createInArgs();
632
633 inArgs.setModelEvalDescription(this->description());
634 inArgs.set_Np(num_p+num_p_sg);
635 inArgs.setSupports(IN_ARG_x_dot, me_inargs.supports(IN_ARG_x_dot_sg));
636 inArgs.setSupports(IN_ARG_x, me_inargs.supports(IN_ARG_x_sg));
637 inArgs.setSupports(IN_ARG_t, me_inargs.supports(IN_ARG_t));
638 inArgs.setSupports(IN_ARG_alpha, me_inargs.supports(IN_ARG_alpha));
639 inArgs.setSupports(IN_ARG_beta, me_inargs.supports(IN_ARG_beta));
640 inArgs.setSupports(IN_ARG_sg_basis, me_inargs.supports(IN_ARG_sg_basis));
641 inArgs.setSupports(IN_ARG_sg_quadrature,
642 me_inargs.supports(IN_ARG_sg_quadrature));
643 inArgs.setSupports(IN_ARG_sg_expansion,
644 me_inargs.supports(IN_ARG_sg_expansion));
645
646 return inArgs;
647}
648
649EpetraExt::ModelEvaluator::OutArgs
651{
652 OutArgsSetup outArgs;
653 OutArgs me_outargs = me->createOutArgs();
654
655 outArgs.setModelEvalDescription(this->description());
656 outArgs.set_Np_Ng(num_p+num_p_sg, num_g_sg);
657 outArgs.setSupports(OUT_ARG_f, me_outargs.supports(OUT_ARG_f_sg));
658 outArgs.setSupports(OUT_ARG_W, me_outargs.supports(OUT_ARG_W_sg));
659 outArgs.setSupports(
660 OUT_ARG_WPrec,
661 me_outargs.supports(OUT_ARG_W_sg) && sg_prec_factory->isPrecSupported());
662 for (int j=0; j<num_p; j++)
663 outArgs.setSupports(OUT_ARG_DfDp, j,
664 me_outargs.supports(OUT_ARG_DfDp_sg, j));
665 for (int j=0; j<num_p_sg; j++)
666 if (!me_outargs.supports(OUT_ARG_DfDp_sg, sg_p_index_map[j]).none())
667 outArgs.setSupports(OUT_ARG_DfDp, j+num_p, DERIV_LINEAR_OP);
668 for (int i=0; i<num_g_sg; i++) {
669 int ii = sg_g_index_map[i];
670 if (!me_outargs.supports(OUT_ARG_DgDx_dot_sg, ii).none())
671 outArgs.setSupports(OUT_ARG_DgDx_dot, i, DERIV_LINEAR_OP);
672 if (!me_outargs.supports(OUT_ARG_DgDx_sg, ii).none())
673 outArgs.setSupports(OUT_ARG_DgDx, i, DERIV_LINEAR_OP);
674 for (int j=0; j<num_p; j++)
675 outArgs.setSupports(OUT_ARG_DgDp, i, j,
676 me_outargs.supports(OUT_ARG_DgDp_sg, ii, j));
677 for (int j=0; j<num_p_sg; j++)
678 if (!me_outargs.supports(OUT_ARG_DgDp_sg, ii, sg_p_index_map[j]).none())
679 outArgs.setSupports(OUT_ARG_DgDp, i, j+num_p, DERIV_LINEAR_OP);
680 }
681
682 return outArgs;
683}
684
685void
687 const OutArgs& outArgs) const
688{
689 // Get the input arguments
690 Teuchos::RCP<const Epetra_Vector> x;
691 if (inArgs.supports(IN_ARG_x)) {
692 x = inArgs.get_x();
693 if (x != Teuchos::null)
694 *my_x = *x;
695 }
696 Teuchos::RCP<const Epetra_Vector> x_dot;
697 if (inArgs.supports(IN_ARG_x_dot))
698 x_dot = inArgs.get_x_dot();
699
700 // Get the output arguments
701 EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> f_out;
702 if (outArgs.supports(OUT_ARG_f))
703 f_out = outArgs.get_f();
704 Teuchos::RCP<Epetra_Operator> W_out;
705 if (outArgs.supports(OUT_ARG_W))
706 W_out = outArgs.get_W();
707 Teuchos::RCP<Epetra_Operator> WPrec_out;
708 if (outArgs.supports(OUT_ARG_WPrec))
709 WPrec_out = outArgs.get_WPrec();
710
711 // Check if we are using the "matrix-free" method for W and we are
712 // computing a preconditioner.
713 bool eval_prec = (W_out == Teuchos::null && WPrec_out != Teuchos::null);
714
715 // Here we are assuming a full W fill occurred previously which we can use
716 // for the preconditioner. Given the expense of computing the SG W blocks
717 // this saves significant computational cost
718 if (eval_prec) {
719 Teuchos::RCP<Stokhos::SGPreconditioner> W_prec =
720 Teuchos::rcp_dynamic_cast<Stokhos::SGPreconditioner>(WPrec_out, true);
721 W_prec->setupPreconditioner(my_W, *my_x);
722
723 // We can now quit unless a fill of f, g, or dg/dp was also requested
724 bool done = (f_out == Teuchos::null);
725 for (int i=0; i<outArgs.Ng(); i++) {
726 done = done && (outArgs.get_g(i) == Teuchos::null);
727 for (int j=0; j<outArgs.Np(); j++)
728 if (!outArgs.supports(OUT_ARG_DgDp, i, j).none())
729 done = done && (outArgs.get_DgDp(i,j).isEmpty());
730 }
731 if (done)
732 return;
733 }
734
735 // Create underlying inargs
736 InArgs me_inargs = me->createInArgs();
737 if (x != Teuchos::null) {
738 x_sg_blocks->getBlockVector()->Import(*x, *sg_overlapped_x_importer,
739 Insert);
740 me_inargs.set_x_sg(x_sg_blocks);
741 }
742 if (x_dot != Teuchos::null) {
743 x_dot_sg_blocks->getBlockVector()->Import(*x_dot, *sg_overlapped_x_importer,
744 Insert);
745 me_inargs.set_x_dot_sg(x_dot_sg_blocks);
746 }
747 if (me_inargs.supports(IN_ARG_alpha))
748 me_inargs.set_alpha(inArgs.get_alpha());
749 if (me_inargs.supports(IN_ARG_beta))
750 me_inargs.set_beta(inArgs.get_beta());
751 if (me_inargs.supports(IN_ARG_t))
752 me_inargs.set_t(inArgs.get_t());
753 if (me_inargs.supports(IN_ARG_sg_basis)) {
754 if (inArgs.get_sg_basis() != Teuchos::null)
755 me_inargs.set_sg_basis(inArgs.get_sg_basis());
756 else
757 me_inargs.set_sg_basis(sg_basis);
758 }
759 if (me_inargs.supports(IN_ARG_sg_quadrature)) {
760 if (inArgs.get_sg_quadrature() != Teuchos::null)
761 me_inargs.set_sg_quadrature(inArgs.get_sg_quadrature());
762 else
763 me_inargs.set_sg_quadrature(sg_quad);
764 }
765 if (me_inargs.supports(IN_ARG_sg_expansion)) {
766 if (inArgs.get_sg_expansion() != Teuchos::null)
767 me_inargs.set_sg_expansion(inArgs.get_sg_expansion());
768 else
769 me_inargs.set_sg_expansion(sg_exp);
770 }
771
772 // Pass parameters
773 for (int i=0; i<num_p; i++)
774 me_inargs.set_p(i, inArgs.get_p(i));
775 for (int i=0; i<num_p_sg; i++) {
776 Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(i+num_p);
777
778 // We always need to pass in the SG parameters, so just use
779 // the initial parameters if it is NULL
780 if (p == Teuchos::null)
781 p = sg_p_init[i]->getBlockVector();
782
783 // Convert block p to SG polynomial
784 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> p_sg =
785 create_p_sg(sg_p_index_map[i], View, p.get());
786 me_inargs.set_p_sg(sg_p_index_map[i], p_sg);
787 }
788
789 // Create underlying outargs
790 OutArgs me_outargs = me->createOutArgs();
791
792 // f
793 if (f_out != Teuchos::null) {
794 me_outargs.set_f_sg(f_sg_blocks);
795 if (eval_W_with_f)
796 me_outargs.set_W_sg(W_sg_blocks);
797 }
798
799 // W
800 if (W_out != Teuchos::null && !eval_W_with_f && !eval_prec)
801 me_outargs.set_W_sg(W_sg_blocks);
802
803 // df/dp -- deterministic p
804 for (int i=0; i<num_p; i++) {
805 if (!outArgs.supports(OUT_ARG_DfDp, i).none()) {
806 Derivative dfdp = outArgs.get_DfDp(i);
807 if (dfdp.getMultiVector() != Teuchos::null) {
808 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dfdp_sg;
809 if (dfdp.getMultiVectorOrientation() == DERIV_MV_BY_COL)
810 dfdp_sg =
812 sg_basis, overlapped_stoch_row_map,
813 me->get_f_map(), sg_overlapped_f_map, sg_comm,
814 me->get_p_map(i)->NumMyElements()));
815 else if (dfdp.getMultiVectorOrientation() == DERIV_TRANS_MV_BY_ROW)
816 dfdp_sg =
818 sg_basis, overlapped_stoch_row_map,
819 me->get_p_map(i), sg_comm,
820 me->get_f_map()->NumMyElements()));
821 me_outargs.set_DfDp_sg(
822 i, SGDerivative(dfdp_sg, dfdp.getMultiVectorOrientation()));
823 }
824 else if (dfdp.getLinearOp() != Teuchos::null) {
825 Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly> dfdp_sg =
826 Teuchos::rcp_dynamic_cast<Stokhos::EpetraOperatorOrthogPoly>(dfdp.getLinearOp(), true);
827 me_outargs.set_DfDp_sg(i, SGDerivative(dfdp_sg));
828 }
829 }
830 }
831
832 // dfdp -- stochastic p. Here we only support DERIV_LINEAR_OP
833 for (int i=0; i<num_p_sg; i++) {
834 if (!outArgs.supports(OUT_ARG_DfDp, i+num_p).none()) {
835 Derivative dfdp = outArgs.get_DfDp(i+num_p);
836 if (dfdp.getLinearOp() != Teuchos::null) {
837 Teuchos::RCP<Stokhos::MatrixFreeOperator> dfdp_op =
838 Teuchos::rcp_dynamic_cast<Stokhos::MatrixFreeOperator>(dfdp.getLinearOp(), true);
839 Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly > dfdp_op_sg =
840 dfdp_op->getSGPolynomial();
841 int ii = sg_p_index_map[i];
842 if (me_outargs.supports(OUT_ARG_DfDp_sg,ii).supports(DERIV_LINEAR_OP))
843 me_outargs.set_DfDp_sg(ii, SGDerivative(dfdp_op_sg));
844 else {
845 Teuchos::RCP<Stokhos::EpetraMultiVectorOperatorOrthogPoly> sg_mv_op =
846 Teuchos::rcp_dynamic_cast<Stokhos::EpetraMultiVectorOperatorOrthogPoly>(dfdp_op_sg, true);
847 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dfdp_sg =
848 sg_mv_op->multiVectorOrthogPoly();
849 if (me_outargs.supports(OUT_ARG_DfDp_sg,ii).supports(DERIV_MV_BY_COL))
850 me_outargs.set_DfDp_sg(
851 ii, SGDerivative(dfdp_sg, DERIV_MV_BY_COL));
852 else
853 me_outargs.set_DfDp_sg(
854 ii, SGDerivative(dfdp_sg, DERIV_TRANS_MV_BY_ROW));
855 }
856 }
857 TEUCHOS_TEST_FOR_EXCEPTION(
858 dfdp.getLinearOp() == Teuchos::null && dfdp.isEmpty() == false,
859 std::logic_error,
860 "Error! Stokhos::SGModelEvaluator::evalModel: " <<
861 "Operator form of df/dp(" << i+num_p << ") is required!");
862 }
863 }
864
865 // Responses (g, dg/dx, dg/dp, ...)
866 for (int i=0; i<num_g_sg; i++) {
867 int ii = sg_g_index_map[i];
868
869 // g
870 Teuchos::RCP<Epetra_Vector> g = outArgs.get_g(i);
871 if (g != Teuchos::null) {
872 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> g_sg =
873 create_g_sg(sg_g_index_map[i], View, g.get());
874 me_outargs.set_g_sg(ii, g_sg);
875 }
876
877 // dg/dxdot
878 if (outArgs.supports(OUT_ARG_DgDx_dot, i).supports(DERIV_LINEAR_OP)) {
879 Derivative dgdx_dot = outArgs.get_DgDx_dot(i);
880 if (dgdx_dot.getLinearOp() != Teuchos::null) {
881 Teuchos::RCP<Stokhos::SGOperator> op =
882 Teuchos::rcp_dynamic_cast<Stokhos::SGOperator>(
883 dgdx_dot.getLinearOp(), true);
884 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > sg_blocks =
885 op->getSGPolynomial();
886 if (me_outargs.supports(OUT_ARG_DgDx, ii).supports(DERIV_LINEAR_OP))
887 me_outargs.set_DgDx_dot_sg(ii, sg_blocks);
888 else {
889 Teuchos::RCP<Stokhos::EpetraMultiVectorOperatorOrthogPoly> sg_mv_op =
890 Teuchos::rcp_dynamic_cast<Stokhos::EpetraMultiVectorOperatorOrthogPoly>(sg_blocks, true);
891 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdx_dot_sg =
892 sg_mv_op->multiVectorOrthogPoly();
893 if (me_outargs.supports(OUT_ARG_DgDx_dot_sg, ii).supports(DERIV_MV_BY_COL))
894 me_outargs.set_DgDx_dot_sg(ii, SGDerivative(dgdx_dot_sg,
895 DERIV_MV_BY_COL));
896 else
897 me_outargs.set_DgDx_dot_sg(ii, SGDerivative(dgdx_dot_sg,
898 DERIV_TRANS_MV_BY_ROW));
899 }
900 }
901 TEUCHOS_TEST_FOR_EXCEPTION(
902 dgdx_dot.getLinearOp() == Teuchos::null && dgdx_dot.isEmpty() == false,
903 std::logic_error,
904 "Error! Stokhos::SGModelEvaluator::evalModel: " <<
905 "Operator form of dg/dxdot(" << i << ") is required!");
906 }
907
908 // dg/dx
909 if (outArgs.supports(OUT_ARG_DgDx, i).supports(DERIV_LINEAR_OP)) {
910 Derivative dgdx = outArgs.get_DgDx(i);
911 if (dgdx.getLinearOp() != Teuchos::null) {
912 Teuchos::RCP<Stokhos::SGOperator> op =
913 Teuchos::rcp_dynamic_cast<Stokhos::SGOperator>(
914 dgdx.getLinearOp(), true);
915 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > sg_blocks =
916 op->getSGPolynomial();
917 if (me_outargs.supports(OUT_ARG_DgDx, ii).supports(DERIV_LINEAR_OP))
918 me_outargs.set_DgDx_sg(ii, sg_blocks);
919 else {
920 Teuchos::RCP<Stokhos::EpetraMultiVectorOperatorOrthogPoly> sg_mv_op =
921 Teuchos::rcp_dynamic_cast<Stokhos::EpetraMultiVectorOperatorOrthogPoly>(sg_blocks, true);
922 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdx_sg =
923 sg_mv_op->multiVectorOrthogPoly();
924 if (me_outargs.supports(OUT_ARG_DgDx_sg, ii).supports(DERIV_MV_BY_COL))
925 me_outargs.set_DgDx_sg(ii, SGDerivative(dgdx_sg,
926 DERIV_MV_BY_COL));
927 else
928 me_outargs.set_DgDx_sg(ii, SGDerivative(dgdx_sg,
929 DERIV_TRANS_MV_BY_ROW));
930 }
931 }
932 TEUCHOS_TEST_FOR_EXCEPTION(
933 dgdx.getLinearOp() == Teuchos::null && dgdx.isEmpty() == false,
934 std::logic_error,
935 "Error! Stokhos::SGModelEvaluator::evalModel: " <<
936 "Operator form of dg/dx(" << i << ") is required!");
937 }
938
939 // dg/dp -- determinsitic p
940 for (int j=0; j<num_p; j++) {
941 if (!outArgs.supports(OUT_ARG_DgDp, i, j).none()) {
942 Derivative dgdp = outArgs.get_DgDp(i,j);
943 if (dgdp.getMultiVector() != Teuchos::null) {
944 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdp_sg;
945 if (dgdp.getMultiVectorOrientation() == DERIV_MV_BY_COL)
946 dgdp_sg =
948 sg_basis, overlapped_stoch_row_map,
949 me->get_g_map(ii), sg_g_map[i], sg_comm,
950 View, *(dgdp.getMultiVector())));
951 else if (dgdp.getMultiVectorOrientation() == DERIV_TRANS_MV_BY_ROW) {
952 Teuchos::RCP<const Epetra_BlockMap> product_map =
953 Teuchos::rcp(&(dgdp.getMultiVector()->Map()),false);
954 dgdp_sg =
956 sg_basis, overlapped_stoch_row_map,
957 me->get_p_map(j), product_map, sg_comm,
958 View, *(dgdp.getMultiVector())));
959 }
960 me_outargs.set_DgDp_sg(
961 ii, j, SGDerivative(dgdp_sg, dgdp.getMultiVectorOrientation()));
962 }
963 else if (dgdp.getLinearOp() != Teuchos::null) {
964 Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly> dgdp_sg =
965 Teuchos::rcp_dynamic_cast<Stokhos::EpetraOperatorOrthogPoly>(dgdp.getLinearOp(), true);
966 me_outargs.set_DgDp_sg(ii, j, SGDerivative(dgdp_sg));
967 }
968 }
969 }
970
971 // dgdp -- stochastic p. Here we only support DERIV_LINEAR_OP
972 for (int j=0; j<num_p_sg; j++) {
973 if (!outArgs.supports(OUT_ARG_DgDp, i, j+num_p).none()) {
974 Derivative dgdp = outArgs.get_DgDp(i,j+num_p);
975 if (dgdp.getLinearOp() != Teuchos::null) {
976 Teuchos::RCP<Stokhos::MatrixFreeOperator> dgdp_op =
977 Teuchos::rcp_dynamic_cast<Stokhos::MatrixFreeOperator>(dgdp.getLinearOp(), true);
978 Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly > dgdp_op_sg =
979 dgdp_op->getSGPolynomial();
980 int jj = sg_p_index_map[j];
981 if (me_outargs.supports(OUT_ARG_DgDp_sg,ii,jj).supports(DERIV_LINEAR_OP))
982 me_outargs.set_DgDp_sg(ii, jj, SGDerivative(dgdp_op_sg));
983 else {
984 Teuchos::RCP<Stokhos::EpetraMultiVectorOperatorOrthogPoly> sg_mv_op =
985 Teuchos::rcp_dynamic_cast<Stokhos::EpetraMultiVectorOperatorOrthogPoly>(dgdp_op_sg, true);
986 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdp_sg =
987 sg_mv_op->multiVectorOrthogPoly();
988 if (me_outargs.supports(OUT_ARG_DgDp_sg,ii,jj).supports(DERIV_MV_BY_COL))
989 me_outargs.set_DgDp_sg(
990 ii, jj, SGDerivative(dgdp_sg, DERIV_MV_BY_COL));
991 else
992 me_outargs.set_DgDp_sg(
993 ii, jj, SGDerivative(dgdp_sg, DERIV_TRANS_MV_BY_ROW));
994 }
995 }
996 TEUCHOS_TEST_FOR_EXCEPTION(
997 dgdp.getLinearOp() == Teuchos::null && dgdp.isEmpty() == false,
998 std::logic_error,
999 "Error! Stokhos::SGModelEvaluator::evalModel: " <<
1000 "Operator form of dg/dp(" << i << "," << j+num_p << ") is required!");
1001 }
1002 }
1003
1004 }
1005
1006 // Compute the functions
1007 me->evalModel(me_inargs, me_outargs);
1008
1009 // Copy block SG components for W
1010 if ((W_out != Teuchos::null || (eval_W_with_f && f_out != Teuchos::null))
1011 && !eval_prec) {
1012 Teuchos::RCP<Epetra_Operator> W;
1013 if (W_out != Teuchos::null)
1014 W = W_out;
1015 else
1016 W = my_W;
1017 Teuchos::RCP<Stokhos::SGOperator> W_sg =
1018 Teuchos::rcp_dynamic_cast<Stokhos::SGOperator>(W, true);
1019 W_sg->setupOperator(W_sg_blocks);
1020
1021 if (WPrec_out != Teuchos::null) {
1022 Teuchos::RCP<Stokhos::SGPreconditioner> W_prec =
1023 Teuchos::rcp_dynamic_cast<Stokhos::SGPreconditioner>(WPrec_out, true);
1024 W_prec->setupPreconditioner(W_sg, *my_x);
1025 }
1026 }
1027
1028 // f
1029 if (f_out!=Teuchos::null){
1030 if (!scaleOP)
1031 for (int i=0; i<sg_basis->size(); i++)
1032 (*f_sg_blocks)[i].Scale(sg_basis->norm_squared(i));
1033 f_out->Export(*(f_sg_blocks->getBlockVector()), *sg_overlapped_f_exporter,
1034 Insert);
1035 }
1036
1037 // df/dp -- deterministic p
1038 for (int i=0; i<num_p; i++) {
1039 if (!outArgs.supports(OUT_ARG_DfDp, i).none()) {
1040 Derivative dfdp = outArgs.get_DfDp(i);
1041 SGDerivative dfdp_sg = me_outargs.get_DfDp_sg(i);
1042 if (dfdp.getMultiVector() != Teuchos::null) {
1043 dfdp.getMultiVector()->Export(
1044 *(dfdp_sg.getMultiVector()->getBlockMultiVector()),
1045 *sg_overlapped_f_exporter, Insert);
1046 }
1047 // need to handle parallel df/dp operator case -- currently we have
1048 // no way to do communication of operators
1049 }
1050 }
1051
1052 // df/dp -- stochastic p
1053 // No need to handle this case as it is handled by the operator
1054
1055 // g, dg/dx, dg/dxdot, dg/dp -- these are all currently serial
1056}
1057
1058void
1060 const Stokhos::EpetraVectorOrthogPoly& x_sg_in)
1061{
1062 *sg_x_init = x_sg_in;
1063}
1064
1065Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly>
1067{
1068 return sg_x_init;
1069}
1070
1071void
1073 int i, const Stokhos::EpetraVectorOrthogPoly& p_sg_in)
1074{
1075 Teuchos::Array<int>::iterator it = std::find(sg_p_index_map.begin(),
1076 sg_p_index_map.end(),
1077 i);
1078 TEUCHOS_TEST_FOR_EXCEPTION(it == sg_p_index_map.end(), std::logic_error,
1079 "Error! Invalid p map index " << i);
1080 int ii = it - sg_p_index_map.begin();
1081 *sg_p_init[ii] = p_sg_in;
1082}
1083
1084Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly>
1086{
1087 Teuchos::Array<int>::const_iterator it = std::find(sg_p_index_map.begin(),
1088 sg_p_index_map.end(),
1089 l);
1090 TEUCHOS_TEST_FOR_EXCEPTION(it == sg_p_index_map.end(), std::logic_error,
1091 "Error! Invalid p map index " << l);
1092 int ll = it - sg_p_index_map.begin();
1093 return sg_p_init[ll];
1094}
1095
1096Teuchos::Array<int>
1098{
1099 return sg_p_index_map;
1100}
1101
1102Teuchos::Array<int>
1104{
1105 return sg_g_index_map;
1106}
1107
1108Teuchos::Array< Teuchos::RCP<const Epetra_Map> >
1110{
1111 Teuchos::Array< Teuchos::RCP<const Epetra_Map> > base_maps(num_g);
1112 for (int i=0; i<num_g; i++)
1113 base_maps[i] = me->get_g_map(i);
1114 return base_maps;
1115 }
1116
1117Teuchos::RCP<const Epetra_BlockMap>
1119{
1120 return overlapped_stoch_row_map;
1121}
1122
1123Teuchos::RCP<const Epetra_BlockMap>
1125{
1126 return sg_overlapped_x_map;
1127}
1128
1129Teuchos::RCP<const Epetra_Import>
1131{
1132 return sg_overlapped_x_importer;
1133}
1134
1135Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1137 const Epetra_Vector* v) const
1138{
1139 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x;
1140 if (v == NULL)
1141 sg_x = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1142 sg_basis, stoch_row_map, x_map, sg_x_map, sg_comm));
1143 else
1144 sg_x = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1145 sg_basis, stoch_row_map, x_map, sg_x_map, sg_comm,
1146 CV, *v));
1147 return sg_x;
1148}
1149
1150Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1152 const Epetra_Vector* v) const
1153{
1154 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x;
1155 if (v == NULL)
1156 sg_x = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1157 sg_basis, overlapped_stoch_row_map, x_map,
1158 sg_overlapped_x_map, sg_comm));
1159 else
1160 sg_x = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1161 sg_basis, overlapped_stoch_row_map, x_map,
1162 sg_overlapped_x_map, sg_comm, CV, *v));
1163 return sg_x;
1164}
1165
1166Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
1168 const Epetra_MultiVector* v) const
1169{
1170 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_x;
1171 if (v == NULL)
1172 sg_x = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
1173 sg_basis, stoch_row_map, x_map, sg_x_map, sg_comm,
1174 num_vecs));
1175 else
1176 sg_x = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
1177 sg_basis, stoch_row_map, x_map, sg_x_map, sg_comm,
1178 CV, *v));
1179 return sg_x;
1180}
1181
1182Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
1184 int num_vecs,
1186 const Epetra_MultiVector* v) const
1187{
1188 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_x;
1189 if (v == NULL)
1190 sg_x = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
1191 sg_basis, overlapped_stoch_row_map, x_map,
1192 sg_overlapped_x_map, sg_comm, num_vecs));
1193 else
1194 sg_x = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
1195 sg_basis, overlapped_stoch_row_map, x_map,
1196 sg_overlapped_x_map, sg_comm, CV, *v));
1197 return sg_x;
1198}
1199
1200Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1202 const Epetra_Vector* v) const
1203{
1204 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_p;
1205 Teuchos::Array<int>::const_iterator it = std::find(sg_p_index_map.begin(),
1206 sg_p_index_map.end(),
1207 l);
1208 TEUCHOS_TEST_FOR_EXCEPTION(it == sg_p_index_map.end(), std::logic_error,
1209 "Error! Invalid p map index " << l);
1210 int ll = it - sg_p_index_map.begin();
1211 if (v == NULL)
1212 sg_p = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1213 sg_basis, overlapped_stoch_p_map, me->get_p_map(l),
1214 sg_p_map[ll], sg_comm));
1215 else
1216 sg_p = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1217 sg_basis, overlapped_stoch_p_map, me->get_p_map(l),
1218 sg_p_map[ll], sg_comm, CV, *v));
1219 return sg_p;
1220}
1221
1222Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1224 const Epetra_Vector* v) const
1225{
1226 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_f;
1227 if (v == NULL)
1228 sg_f = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1229 sg_basis, stoch_row_map, f_map, sg_f_map, sg_comm));
1230 else
1231 sg_f = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1232 sg_basis, stoch_row_map, f_map, sg_f_map, sg_comm,
1233 CV, *v));
1234 return sg_f;
1235}
1236
1237Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1239 const Epetra_Vector* v) const
1240{
1241 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_f;
1242 if (v == NULL)
1243 sg_f = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1244 sg_basis, overlapped_stoch_row_map, f_map,
1245 sg_overlapped_f_map, sg_comm));
1246 else
1247 sg_f = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1248 sg_basis, overlapped_stoch_row_map, f_map,
1249 sg_overlapped_f_map, sg_comm, CV, *v));
1250 return sg_f;
1251}
1252
1253Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
1255 int num_vecs,
1257 const Epetra_MultiVector* v) const
1258{
1259 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_f;
1260 if (v == NULL)
1261 sg_f = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
1262 sg_basis, stoch_row_map, f_map, sg_f_map, sg_comm,
1263 num_vecs));
1264 else
1265 sg_f = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
1266 sg_basis, stoch_row_map, f_map, sg_f_map, sg_comm,
1267 CV, *v));
1268 return sg_f;
1269}
1270
1271Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
1273 int num_vecs,
1275 const Epetra_MultiVector* v) const
1276{
1277 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_f;
1278 if (v == NULL)
1279 sg_f = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
1280 sg_basis, overlapped_stoch_row_map, f_map,
1281 sg_overlapped_f_map, sg_comm, num_vecs));
1282 else
1283 sg_f = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
1284 sg_basis, overlapped_stoch_row_map, f_map,
1285 sg_overlapped_f_map, sg_comm, CV, *v));
1286 return sg_f;
1287}
1288
1289Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1291 const Epetra_Vector* v) const
1292{
1293 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_g;
1294 Teuchos::Array<int>::const_iterator it = std::find(sg_g_index_map.begin(),
1295 sg_g_index_map.end(),
1296 l);
1297 TEUCHOS_TEST_FOR_EXCEPTION(it == sg_g_index_map.end(), std::logic_error,
1298 "Error! Invalid g map index " << l);
1299 int ll = it - sg_g_index_map.begin();
1300 if (v == NULL)
1301 sg_g = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1302 sg_basis, overlapped_stoch_row_map,
1303 me->get_g_map(l),
1304 sg_g_map[ll], sg_comm));
1305 else
1306 sg_g = Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1307 sg_basis, overlapped_stoch_row_map,
1308 me->get_g_map(l),
1309 sg_g_map[ll], sg_comm, CV, *v));
1310 return sg_g;
1311}
1312
1313Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
1316 const Epetra_MultiVector* v) const
1317{
1318 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_g;
1319 Teuchos::Array<int>::const_iterator it = std::find(sg_g_index_map.begin(),
1320 sg_g_index_map.end(),
1321 l);
1322 TEUCHOS_TEST_FOR_EXCEPTION(it == sg_g_index_map.end(), std::logic_error,
1323 "Error! Invalid g map index " << l);
1324 int ll = it - sg_g_index_map.begin();
1325 if (v == NULL)
1326 sg_g = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
1327 sg_basis, overlapped_stoch_row_map,
1328 me->get_g_map(l),
1329 sg_g_map[ll], sg_comm, num_vecs));
1330 else
1331 sg_g = Teuchos::rcp(new Stokhos::EpetraMultiVectorOrthogPoly(
1332 sg_basis, overlapped_stoch_row_map,
1333 me->get_g_map(l),
1334 sg_g_map[ll], sg_comm, CV, *v));
1335 return sg_g;
1336}
1337
1338Teuchos::RCP<EpetraExt::BlockVector>
1340{
1341 Teuchos::RCP<EpetraExt::BlockVector> x_overlapped =
1342 Teuchos::rcp(new EpetraExt::BlockVector(*x_map, *sg_overlapped_x_map));
1343 x_overlapped->Import(x, *sg_overlapped_x_importer, Insert);
1344 return x_overlapped;
1345}
1346
1347Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1349{
1350 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x_overlapped =
1351 Teuchos::rcp(new Stokhos::EpetraVectorOrthogPoly(
1352 sg_basis, overlapped_stoch_row_map, x_map,
1353 sg_overlapped_x_map, sg_comm));
1354 sg_x_overlapped->getBlockVector()->Import(x, *sg_overlapped_x_importer,
1355 Insert);
1356 return sg_x_overlapped;
1357}
1358
1359Teuchos::RCP<EpetraExt::BlockVector>
1361{
1362 Teuchos::RCP<EpetraExt::BlockVector> x =
1363 Teuchos::rcp(new EpetraExt::BlockVector(*x_map, *sg_x_map));
1364 x->Export(x_overlapped, *sg_overlapped_x_importer, Insert);
1365 return x;
1366}
1367
1368Teuchos::RCP<EpetraExt::BlockVector>
1370{
1371 Teuchos::RCP<EpetraExt::BlockVector> f_overlapped =
1372 Teuchos::rcp(new EpetraExt::BlockVector(*f_map, *sg_overlapped_f_map));
1373 f_overlapped->Import(f, *sg_overlapped_f_exporter, Insert);
1374 return f_overlapped;
1375}
1376
1377Teuchos::RCP<EpetraExt::BlockVector>
1379{
1380 Teuchos::RCP<EpetraExt::BlockVector> f =
1381 Teuchos::rcp(new EpetraExt::BlockVector(*f_map, *sg_f_map));
1382 f->Export(f_overlapped, *sg_overlapped_f_exporter, Insert);
1383 return f;
1384}
Insert
Epetra_DataAccess
A container class storing an orthogonal polynomial whose coefficients are vectors,...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
An Epetra operator representing the block stochastic Galerkin operator.
Abstract base class for multivariate orthogonal polynomials.
Abstract base class for orthogonal polynomial-based expansions.
Abstract base class for quadrature methods.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > sg_basis
Stochastic Galerkin basis.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_x_sg_overlap(Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create vector orthog poly using x map and overlap sg map.
Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > create_f_mv_sg_overlap(int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create multi-vector orthog poly using f map and overlap sg map.
OutArgs createOutArgs() const
Create OutArgs.
Teuchos::RCP< EpetraExt::BlockVector > export_residual(const Epetra_Vector &f_overlapped) const
Export parallel residual vector.
Teuchos::RCP< const Stokhos::ParallelData > sg_parallel_data
Parallel SG data.
Teuchos::RCP< const Epetra_Map > sg_f_map
Block SG residual map.
int num_p_sg
Number of stochastic parameter vectors.
Teuchos::RCP< Epetra_Operator > create_DgDx_dot_op(int j) const
Create SG operator representing dg/dxdot.
Teuchos::RCP< EpetraExt::ModelEvaluator > me
Underlying model evaluator.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > sg_x_init
SG initial x.
int num_g_sg
Number of stochastic response vectors.
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
Teuchos::RCP< Stokhos::SGPreconditionerFactory > sg_prec_factory
Preconditioner factory.
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > serialCijk
Serial Epetra Cijk for dgdx*.
Teuchos::RCP< const Epetra_Map > get_g_map(int l) const
Return response map.
Teuchos::Array< int > get_g_sg_map_indices() const
Get indices of SG responses.
Teuchos::RCP< const Epetra_Map > f_map
Underlying residual map.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_f_sg_overlap(Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create vector orthog poly using f map and overlap sg map.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
unsigned int num_sg_blocks
Number of stochastic blocks.
Teuchos::RCP< const Epetra_Map > sg_overlapped_x_map
Block SG overlapped unknown map.
Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > create_f_mv_sg(int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create multi-vector orthog poly using f map and owned sg map.
Teuchos::RCP< const Epetra_BlockMap > overlapped_stoch_row_map
Overlapped map for stochastic blocks (local map)
InArgs createInArgs() const
Create InArgs.
Teuchos::RCP< const Epetra_Map > sg_x_map
Block SG unknown map.
Teuchos::RCP< EpetraExt::BlockVector > import_solution(const Epetra_Vector &x) const
Import parallel solution vector.
Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > create_x_mv_sg_overlap(int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create vector orthog poly using x map and overlap sg map.
Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > create_x_mv_sg(int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create vector orthog poly using x map and owned sg map.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
void set_x_sg_init(const Stokhos::EpetraVectorOrthogPoly &x_sg_in)
Set initial solution polynomial.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > sg_g_map
Block SG response map.
Teuchos::RCP< Epetra_Import > sg_overlapped_x_importer
Importer from SG to SG-overlapped maps.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_x_sg(Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create vector orthog poly using x map and owned sg map.
Teuchos::RCP< const Epetra_BlockMap > overlapped_stoch_p_map
Overlapped map for p stochastic blocks (local map)
bool eval_W_with_f
Whether to always evaluate W with f.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_g_sg(int l, Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create vector orthog poly using g map.
Teuchos::RCP< Teuchos::ParameterList > params
Algorithmic parameters.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > x_sg_blocks
x stochastic Galerkin components
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_p_sg(int l, Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create vector orthog poly using p map.
Teuchos::Array< int > sg_g_index_map
Index map between block-g and g_sg maps.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > sg_p_map
Block SG parameter map.
Teuchos::RCP< EpetraExt::BlockVector > export_solution(const Epetra_Vector &x_overlapped) const
Export parallel solution vector.
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
Teuchos::Array< Teuchos::RCP< Teuchos::Array< std::string > > > sg_p_names
SG coefficient parameter names.
Teuchos::RCP< Epetra_Vector > my_x
x pointer for evaluating preconditioner
Teuchos::RCP< Epetra_Export > sg_overlapped_f_exporter
Exporter from SG-overlapped to SG maps.
Teuchos::Array< int > sg_p_index_map
Index map between block-p and p_sg maps.
Teuchos::RCP< const Epetra_Import > get_x_sg_importer() const
Return x sg importer.
Teuchos::RCP< Epetra_Operator > create_DgDx_op(int j) const
Create SG operator representing dg/dx.
Teuchos::RCP< const EpetraExt::MultiComm > sg_comm
Parallel SG communicator.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > import_solution_poly(const Epetra_Vector &x) const
Import parallel solution vector.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
Teuchos::RCP< Epetra_Operator > create_DgDp_op(int j, int i) const
Create SG operator representing dg/dp.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > x_dot_sg_blocks
x_dot stochastic Galerkin components
Teuchos::RCP< EpetraExt::ModelEvaluator::Preconditioner > create_WPrec() const
Create preconditioner operator.
Teuchos::RCP< const Epetra_BlockMap > get_overlap_stochastic_map() const
Return overlap stochastic map.
Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > W_sg_blocks
W stochastic Galerkin components.
Teuchos::Array< int > get_p_sg_map_indices() const
Get indices of SG parameters.
Teuchos::RCP< const Epetra_BlockMap > stoch_row_map
Map for stochastic blocks.
SGModelEvaluator(const Teuchos::RCP< EpetraExt::ModelEvaluator > &me, const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &sg_basis, const Teuchos::RCP< const Stokhos::Quadrature< int, double > > &sg_quad, const Teuchos::RCP< Stokhos::OrthogPolyExpansion< int, double > > &sg_exp, const Teuchos::RCP< const Stokhos::ParallelData > &sg_parallel_data, const Teuchos::RCP< Teuchos::ParameterList > &params, bool scaleOP=true)
Teuchos::Array< Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > > sg_p_init
SG initial p.
Teuchos::RCP< const Stokhos::EpetraVectorOrthogPoly > get_x_sg_init() const
Return initial SG x.
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
void set_p_sg_init(int i, const Stokhos::EpetraVectorOrthogPoly &p_sg_in)
Set initial parameter polynomial.
Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > create_g_mv_sg(int l, int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create multi-vector orthog poly using g map.
Teuchos::RCP< const Epetra_BlockMap > get_x_sg_overlap_map() const
Return x sg overlap map.
int num_p
Number of parameter vectors of underlying model evaluator.
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > epetraCijk
Epetra Cijk.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > create_f_sg(Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create vector orthog poly using f map and owned sg map.
bool supports_x
Whether we support x (and thus f and W)
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > get_g_sg_base_maps() const
Get base maps of SG responses.
unsigned int num_W_blocks
Number of W stochastic blocks (may be smaller than num_sg_blocks)
Teuchos::RCP< EpetraExt::BlockVector > import_residual(const Epetra_Vector &f) const
Import parallel residual vector.
Teuchos::RCP< const Epetra_Map > x_map
Underlying unknown map.
Teuchos::RCP< Epetra_Operator > create_DfDp_op(int i) const
Create SG operator representing df/dp.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > f_sg_blocks
f stochastic Galerkin components
Teuchos::RCP< const Stokhos::EpetraVectorOrthogPoly > get_p_sg_init(int l) const
Return initial SG parameters.
unsigned int num_p_blocks
Number of p stochastic blocks (may be smaller than num_sg_blocks)
int num_g
Number of response vectors of underlying model evaluator.
Teuchos::RCP< const Epetra_Map > sg_overlapped_f_map
Block SG overlapped residual map.
Factory for generating stochastic Galerkin preconditioners.
virtual Teuchos::RCP< Stokhos::SGOperator > build(const Teuchos::RCP< const EpetraExt::MultiComm > &sg_comm, const Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > &sg_basis, const Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > &epetraCijk, const Teuchos::RCP< const Epetra_Map > &domain_base_map, const Teuchos::RCP< const Epetra_Map > &range_base_map, const Teuchos::RCP< const Epetra_Map > &domain_sg_map, const Teuchos::RCP< const Epetra_Map > &range_sg_map)
Build preconditioner operator.
Factory for generating stochastic Galerkin preconditioners.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)