45 #include "Teuchos_Assert.hpp" 46 #include "EpetraExt_BlockUtility.h" 47 #include "EpetraExt_BlockMultiVector.h" 57 const Teuchos::RCP<EpetraExt::ModelEvaluator>& me_,
61 const Teuchos::RCP<const Stokhos::ParallelData>& sg_parallel_data_,
62 const Teuchos::RCP<Teuchos::ParameterList>& params_,
69 num_sg_blocks(sg_basis->size()),
70 num_W_blocks(sg_basis->size()),
71 num_p_blocks(sg_basis->size()),
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()),
80 sg_overlapped_x_map(),
82 sg_overlapped_f_map(),
83 sg_overlapped_x_importer(),
84 sg_overlapped_f_exporter(),
100 eval_W_with_f(false),
103 if (
x_map != Teuchos::null)
117 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
120 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
124 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
127 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
150 std::string W_expansion_type =
151 params->get(
"Jacobian Expansion Type",
"Full");
152 if (W_expansion_type ==
"Linear")
156 Teuchos::RCP<Epetra_BlockMap> W_overlap_map =
168 Teuchos::RCP<Teuchos::ParameterList> sgPrecParams =
169 Teuchos::rcp(&(
params->sublist(
"SG Preconditioner")),
false);
177 InArgs me_inargs =
me->createInArgs();
178 OutArgs me_outargs =
me->createOutArgs();
179 num_p = me_inargs.Np();
182 for (
int i=0; i<
num_p; i++) {
183 if (me_inargs.supports(IN_ARG_p_sg, i))
193 std::string p_expansion_type =
194 params->get(
"Parameter Expansion Type",
"Full");
195 if (p_expansion_type ==
"Linear")
208 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
211 Teuchos::RCP<const Teuchos::Array<std::string> > p_names =
213 if (p_names != Teuchos::null) {
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;
231 num_g = me_outargs.Ng();
234 for (
int i=0; i<
num_g; i++) {
235 if (me_outargs.supports(OUT_ARG_g_sg, i))
246 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
263 Teuchos::RCP<const Epetra_Map>
269 Teuchos::RCP<const Epetra_Map>
275 Teuchos::RCP<const Epetra_Map>
278 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_sg, std::logic_error,
279 "Error! Invalid p map index " << l);
281 return me->get_p_map(l);
283 return sg_p_map[l-num_p];
285 return Teuchos::null;
288 Teuchos::RCP<const Epetra_Map>
291 TEUCHOS_TEST_FOR_EXCEPTION(j < 0 || j >= num_g_sg, std::logic_error,
292 "Error! Invalid g map index " <<
j);
296 Teuchos::RCP<const Teuchos::Array<std::string> >
299 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_sg, std::logic_error,
300 "Error! Invalid p map index " << l);
302 return me->get_p_names(l);
304 return sg_p_names[l-num_p];
306 return Teuchos::null;
309 Teuchos::RCP<const Epetra_Vector>
312 return sg_x_init->getBlockVector();
315 Teuchos::RCP<const Epetra_Vector>
318 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_sg, std::logic_error,
319 "Error! Invalid p map index " << l);
321 return me->get_p_init(l);
323 return sg_p_init[l-num_p]->getBlockVector();
325 return Teuchos::null;
328 Teuchos::RCP<Epetra_Operator>
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") ==
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",
344 my_W = sg_op_factory.
build(sg_comm, sg_basis, epetraCijk, x_map, f_map,
346 my_W->setupOperator(W_sg_blocks);
351 return Teuchos::null;
354 Teuchos::RCP<EpetraExt::ModelEvaluator::Preconditioner>
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,
365 return Teuchos::null;
368 Teuchos::RCP<Epetra_Operator>
371 TEUCHOS_TEST_FOR_EXCEPTION(j < 0 || j >= num_g_sg || !supports_x,
373 "Error: dg/dx index " <<
j <<
" is not supported!");
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)) {
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));
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()));
395 sg_mv_blocks,
false));
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()));
404 sg_mv_blocks,
true));
407 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
408 "Error! me_outargs.supports(OUT_ARG_DgDx_sg, " << jj
409 <<
").none() is true!");
411 Teuchos::RCP<Teuchos::ParameterList> pl =
412 Teuchos::rcp(
new Teuchos::ParameterList);
413 Teuchos::RCP<Stokhos::SGOperator> dgdx_sg =
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);
421 Teuchos::RCP<Epetra_Operator>
424 TEUCHOS_TEST_FOR_EXCEPTION(j < 0 || j >= num_g_sg || !supports_x,
426 "Error: dg/dx_dot index " <<
j <<
" is not supported!");
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)) {
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));
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()));
448 sg_mv_blocks,
false));
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()));
457 sg_mv_blocks,
true));
460 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
461 "Error! me_outargs.supports(OUT_ARG_DgDx_dot_sg, " 462 << jj <<
").none() is true!");
464 Teuchos::RCP<Teuchos::ParameterList> pl =
465 Teuchos::rcp(
new Teuchos::ParameterList);
466 Teuchos::RCP<Stokhos::SGOperator> dgdx_dot_sg =
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);
474 Teuchos::RCP<Epetra_Operator>
477 TEUCHOS_TEST_FOR_EXCEPTION(
478 j < 0 || j >= num_g_sg || i < 0 || i >= num_p+num_p_sg,
480 "Error: dg/dp index " <<
j <<
"," << i <<
" is not supported!");
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();
486 if (me_outargs.supports(OUT_ARG_DgDp_sg,jj,i).supports(DERIV_LINEAR_OP)) {
487 Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly> sg_blocks =
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));
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 <<
"!");
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)) {
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));
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()));
522 sg_mv_blocks,
false));
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,
529 sg_comm, g_map->NumMyElements()));
532 sg_mv_blocks,
true));
535 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
536 "Error! me_outargs.supports(OUT_ARG_DgDp_sg, " << jj
537 <<
"," << ii <<
").none() is true!");
539 Teuchos::RCP<Teuchos::ParameterList> pl =
540 Teuchos::rcp(
new Teuchos::ParameterList);
541 Teuchos::RCP<Stokhos::SGOperator> dgdp_sg =
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);
549 return Teuchos::null;
552 Teuchos::RCP<Epetra_Operator>
555 TEUCHOS_TEST_FOR_EXCEPTION(i < 0 || i >= num_p+num_p_sg,
557 "Error: df/dp index " << i <<
" is not supported!");
559 OutArgs me_outargs = me->createOutArgs();
561 if (me_outargs.supports(OUT_ARG_DfDp_sg,i).supports(DERIV_LINEAR_OP)) {
562 Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly> sg_blocks =
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));
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 <<
"!");
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)) {
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));
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()));
597 sg_mv_blocks,
false));
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,
604 sg_comm, f_map->NumMyElements()));
607 sg_mv_blocks,
true));
610 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
611 "Error! me_outargs.supports(OUT_ARG_DfDp_sg, " << ii
612 <<
").none() is true!");
614 Teuchos::RCP<Teuchos::ParameterList> pl =
615 Teuchos::rcp(
new Teuchos::ParameterList);
616 Teuchos::RCP<Stokhos::SGOperator> dfdp_sg =
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);
624 return Teuchos::null;
627 EpetraExt::ModelEvaluator::InArgs
631 InArgs me_inargs = me->createInArgs();
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));
649 EpetraExt::ModelEvaluator::OutArgs
652 OutArgsSetup outArgs;
653 OutArgs me_outargs = me->createOutArgs();
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));
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);
687 const OutArgs& outArgs)
const 690 Teuchos::RCP<const Epetra_Vector>
x;
691 if (inArgs.supports(IN_ARG_x)) {
693 if (
x != Teuchos::null)
696 Teuchos::RCP<const Epetra_Vector> x_dot;
697 if (inArgs.supports(IN_ARG_x_dot))
698 x_dot = inArgs.get_x_dot();
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();
713 bool eval_prec = (W_out == Teuchos::null && WPrec_out != Teuchos::null);
719 Teuchos::RCP<Stokhos::SGPreconditioner> W_prec =
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());
736 InArgs me_inargs = me->createInArgs();
737 if (
x != Teuchos::null) {
738 x_sg_blocks->getBlockVector()->Import(*
x, *sg_overlapped_x_importer,
740 me_inargs.set_x_sg(x_sg_blocks);
742 if (x_dot != Teuchos::null) {
743 x_dot_sg_blocks->getBlockVector()->Import(*x_dot, *sg_overlapped_x_importer,
745 me_inargs.set_x_dot_sg(x_dot_sg_blocks);
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());
757 me_inargs.set_sg_basis(sg_basis);
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());
763 me_inargs.set_sg_quadrature(sg_quad);
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());
769 me_inargs.set_sg_expansion(sg_exp);
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);
780 if (p == Teuchos::null)
781 p = sg_p_init[i]->getBlockVector();
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);
790 OutArgs me_outargs = me->createOutArgs();
793 if (f_out != Teuchos::null) {
794 me_outargs.set_f_sg(f_sg_blocks);
796 me_outargs.set_W_sg(W_sg_blocks);
800 if (W_out != Teuchos::null && !eval_W_with_f && !eval_prec)
801 me_outargs.set_W_sg(W_sg_blocks);
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)
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)
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()));
824 else if (dfdp.getLinearOp() != Teuchos::null) {
825 Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly> dfdp_sg =
827 me_outargs.set_DfDp_sg(i, SGDerivative(dfdp_sg));
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 =
839 Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly > dfdp_op_sg =
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));
845 Teuchos::RCP<Stokhos::EpetraMultiVectorOperatorOrthogPoly> sg_mv_op =
847 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dfdp_sg =
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));
853 me_outargs.set_DfDp_sg(
854 ii, SGDerivative(dfdp_sg, DERIV_TRANS_MV_BY_ROW));
857 TEUCHOS_TEST_FOR_EXCEPTION(
858 dfdp.getLinearOp() == Teuchos::null && dfdp.isEmpty() ==
false,
860 "Error! Stokhos::SGModelEvaluator::evalModel: " <<
861 "Operator form of df/dp(" << i+num_p <<
") is required!");
866 for (
int i=0; i<num_g_sg; i++) {
867 int ii = sg_g_index_map[i];
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);
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 =
883 dgdx_dot.getLinearOp(),
true);
884 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > sg_blocks =
886 if (me_outargs.supports(OUT_ARG_DgDx, ii).supports(DERIV_LINEAR_OP))
887 me_outargs.set_DgDx_dot_sg(ii, sg_blocks);
889 Teuchos::RCP<Stokhos::EpetraMultiVectorOperatorOrthogPoly> sg_mv_op =
891 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdx_dot_sg =
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,
897 me_outargs.set_DgDx_dot_sg(ii, SGDerivative(dgdx_dot_sg,
898 DERIV_TRANS_MV_BY_ROW));
901 TEUCHOS_TEST_FOR_EXCEPTION(
902 dgdx_dot.getLinearOp() == Teuchos::null && dgdx_dot.isEmpty() ==
false,
904 "Error! Stokhos::SGModelEvaluator::evalModel: " <<
905 "Operator form of dg/dxdot(" << i <<
") is required!");
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 =
914 dgdx.getLinearOp(),
true);
915 Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > sg_blocks =
917 if (me_outargs.supports(OUT_ARG_DgDx, ii).supports(DERIV_LINEAR_OP))
918 me_outargs.set_DgDx_sg(ii, sg_blocks);
920 Teuchos::RCP<Stokhos::EpetraMultiVectorOperatorOrthogPoly> sg_mv_op =
922 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdx_sg =
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,
928 me_outargs.set_DgDx_sg(ii, SGDerivative(dgdx_sg,
929 DERIV_TRANS_MV_BY_ROW));
932 TEUCHOS_TEST_FOR_EXCEPTION(
933 dgdx.getLinearOp() == Teuchos::null && dgdx.isEmpty() ==
false,
935 "Error! Stokhos::SGModelEvaluator::evalModel: " <<
936 "Operator form of dg/dx(" << i <<
") is required!");
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)
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);
956 sg_basis, overlapped_stoch_row_map,
957 me->get_p_map(
j), product_map, sg_comm,
958 View, *(dgdp.getMultiVector())));
960 me_outargs.set_DgDp_sg(
961 ii,
j, SGDerivative(dgdp_sg, dgdp.getMultiVectorOrientation()));
963 else if (dgdp.getLinearOp() != Teuchos::null) {
964 Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly> dgdp_sg =
966 me_outargs.set_DgDp_sg(ii,
j, SGDerivative(dgdp_sg));
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 =
978 Teuchos::RCP<Stokhos::EpetraOperatorOrthogPoly > dgdp_op_sg =
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));
984 Teuchos::RCP<Stokhos::EpetraMultiVectorOperatorOrthogPoly> sg_mv_op =
986 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> dgdp_sg =
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));
992 me_outargs.set_DgDp_sg(
993 ii, jj, SGDerivative(dgdp_sg, DERIV_TRANS_MV_BY_ROW));
996 TEUCHOS_TEST_FOR_EXCEPTION(
997 dgdp.getLinearOp() == Teuchos::null && dgdp.isEmpty() ==
false,
999 "Error! Stokhos::SGModelEvaluator::evalModel: " <<
1000 "Operator form of dg/dp(" << i <<
"," <<
j+num_p <<
") is required!");
1007 me->evalModel(me_inargs, me_outargs);
1010 if ((W_out != Teuchos::null || (eval_W_with_f && f_out != Teuchos::null))
1012 Teuchos::RCP<Epetra_Operator> W;
1013 if (W_out != Teuchos::null)
1017 Teuchos::RCP<Stokhos::SGOperator> W_sg =
1021 if (WPrec_out != Teuchos::null) {
1022 Teuchos::RCP<Stokhos::SGPreconditioner> W_prec =
1029 if (f_out!=Teuchos::null){
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,
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);
1062 *sg_x_init = x_sg_in;
1065 Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly>
1075 Teuchos::Array<int>::iterator it = std::find(sg_p_index_map.begin(),
1076 sg_p_index_map.end(),
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;
1084 Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly>
1087 Teuchos::Array<int>::const_iterator it = std::find(sg_p_index_map.begin(),
1088 sg_p_index_map.end(),
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];
1099 return sg_p_index_map;
1105 return sg_g_index_map;
1108 Teuchos::Array< Teuchos::RCP<const Epetra_Map> >
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);
1117 Teuchos::RCP<const Epetra_BlockMap>
1120 return overlapped_stoch_row_map;
1123 Teuchos::RCP<const Epetra_BlockMap>
1126 return sg_overlapped_x_map;
1129 Teuchos::RCP<const Epetra_Import>
1132 return sg_overlapped_x_importer;
1135 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1139 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x;
1142 sg_basis, stoch_row_map, x_map, sg_x_map, sg_comm));
1145 sg_basis, stoch_row_map, x_map, sg_x_map, sg_comm,
1150 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1154 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x;
1157 sg_basis, overlapped_stoch_row_map, x_map,
1158 sg_overlapped_x_map, sg_comm));
1161 sg_basis, overlapped_stoch_row_map, x_map,
1162 sg_overlapped_x_map, sg_comm, CV, *v));
1166 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
1170 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_x;
1173 sg_basis, stoch_row_map, x_map, sg_x_map, sg_comm,
1177 sg_basis, stoch_row_map, x_map, sg_x_map, sg_comm,
1182 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
1188 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_x;
1191 sg_basis, overlapped_stoch_row_map, x_map,
1192 sg_overlapped_x_map, sg_comm, num_vecs));
1195 sg_basis, overlapped_stoch_row_map, x_map,
1196 sg_overlapped_x_map, sg_comm, CV, *v));
1200 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
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(),
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();
1213 sg_basis, overlapped_stoch_p_map, me->get_p_map(l),
1214 sg_p_map[ll], sg_comm));
1217 sg_basis, overlapped_stoch_p_map, me->get_p_map(l),
1218 sg_p_map[ll], sg_comm, CV, *v));
1222 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1226 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_f;
1229 sg_basis, stoch_row_map, f_map, sg_f_map, sg_comm));
1232 sg_basis, stoch_row_map, f_map, sg_f_map, sg_comm,
1237 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1241 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_f;
1244 sg_basis, overlapped_stoch_row_map, f_map,
1245 sg_overlapped_f_map, sg_comm));
1248 sg_basis, overlapped_stoch_row_map, f_map,
1249 sg_overlapped_f_map, sg_comm, CV, *v));
1253 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
1259 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_f;
1262 sg_basis, stoch_row_map, f_map, sg_f_map, sg_comm,
1266 sg_basis, stoch_row_map, f_map, sg_f_map, sg_comm,
1271 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
1277 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> sg_f;
1280 sg_basis, overlapped_stoch_row_map, f_map,
1281 sg_overlapped_f_map, sg_comm, num_vecs));
1284 sg_basis, overlapped_stoch_row_map, f_map,
1285 sg_overlapped_f_map, sg_comm, CV, *v));
1289 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
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(),
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();
1302 sg_basis, overlapped_stoch_row_map,
1304 sg_g_map[ll], sg_comm));
1307 sg_basis, overlapped_stoch_row_map,
1309 sg_g_map[ll], sg_comm, CV, *v));
1313 Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly>
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(),
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();
1327 sg_basis, overlapped_stoch_row_map,
1329 sg_g_map[ll], sg_comm, num_vecs));
1332 sg_basis, overlapped_stoch_row_map,
1334 sg_g_map[ll], sg_comm, CV, *v));
1338 Teuchos::RCP<EpetraExt::BlockVector>
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;
1347 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly>
1350 Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> sg_x_overlapped =
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,
1356 return sg_x_overlapped;
1359 Teuchos::RCP<EpetraExt::BlockVector>
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);
1368 Teuchos::RCP<EpetraExt::BlockVector>
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;
1377 Teuchos::RCP<EpetraExt::BlockVector>
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);
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > epetraCijk
Epetra Cijk.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > import_solution_poly(const Epetra_Vector &x) const
Import parallel solution vector.
Teuchos::RCP< EpetraExt::BlockVector > import_residual(const Epetra_Vector &f) const
Import parallel residual vector.
Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > W_sg_blocks
W stochastic Galerkin components.
An Epetra operator representing the block stochastic Galerkin operator.
Teuchos::RCP< const Epetra_Map > sg_overlapped_f_map
Block SG overlapped residual map.
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
virtual void setupOperator(const Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > &poly)=0
Setup operator.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > sg_x_init
SG initial x.
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
Teuchos::Array< int > sg_g_index_map
Index map between block-g and g_sg maps.
unsigned int num_sg_blocks
Number of stochastic blocks.
Teuchos::RCP< const Epetra_Map > x_map
Underlying unknown map.
Teuchos::RCP< const Stokhos::EpetraVectorOrthogPoly > get_p_sg_init(int l) const
Return initial SG parameters.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > x_sg_blocks
x stochastic Galerkin components
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > f_sg_blocks
f stochastic Galerkin components
int num_p_sg
Number of stochastic parameter vectors.
bool eval_W_with_f
Whether to always evaluate W with f.
Teuchos::RCP< const Epetra_Map > sg_overlapped_x_map
Block SG overlapped unknown map.
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 > ¶ms, bool scaleOP=true)
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
Teuchos::RCP< const Epetra_Map > get_g_map(int l) const
Return response map.
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > sg_g_map
Block SG response 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 > get_overlap_stochastic_map() const
Return overlap stochastic map.
Teuchos::RCP< EpetraExt::ModelEvaluator::Preconditioner > create_WPrec() const
Create preconditioner operator.
unsigned int num_p_blocks
Number of p stochastic blocks (may be smaller than num_sg_blocks)
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
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< 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< const Stokhos::ParallelData > sg_parallel_data
Parallel SG data.
Teuchos::RCP< const Epetra_BlockMap > overlapped_stoch_p_map
Overlapped map for p stochastic blocks (local map)
Teuchos::RCP< const Stokhos::EpetraVectorOrthogPoly > get_x_sg_init() const
Return initial SG x.
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< Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > > sg_p_init
SG initial p.
Teuchos::RCP< const Stokhos::EpetraSparse3Tensor > serialCijk
Serial Epetra Cijk for dgdx*.
Teuchos::RCP< Epetra_Operator > create_DfDp_op(int i) const
Create SG operator representing df/dp.
Teuchos::RCP< const Epetra_Map > sg_f_map
Block SG residual map.
Abstract base class for orthogonal polynomial-based expansions.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
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.
Teuchos::RCP< const Epetra_BlockMap > get_x_sg_overlap_map() const
Return x sg overlap 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.
int num_g_sg
Number of stochastic response vectors.
Teuchos::RCP< const Epetra_BlockMap > stoch_row_map
Map for stochastic blocks.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
An abstract class to represent a generic stochastic Galerkin operator as an Epetra_Operator.
Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > x_dot_sg_blocks
x_dot stochastic Galerkin components
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > get_g_sg_base_maps() const
Get base maps of SG responses.
Abstract base class for quadrature methods.
Teuchos::RCP< Stokhos::SGPreconditionerFactory > sg_prec_factory
Preconditioner factory.
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< Epetra_Operator > create_DgDx_dot_op(int j) const
Create SG operator representing dg/dxdot.
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
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.
virtual Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > getSGPolynomial()=0
Get SG polynomial.
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< int, double > > sg_basis
Stochastic Galerkin basis.
Teuchos::RCP< EpetraExt::BlockVector > import_solution(const Epetra_Vector &x) const
Import parallel solution vector.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > sg_p_map
Block SG parameter map.
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
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.
bool supports_x
Whether we support x (and thus f and W)
void set_p_sg_init(int i, const Stokhos::EpetraVectorOrthogPoly &p_sg_in)
Set initial parameter polynomial.
Factory for generating stochastic Galerkin preconditioners.
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
Teuchos::RCP< Epetra_Export > sg_overlapped_f_exporter
Exporter from SG-overlapped to SG maps.
virtual Teuchos::RCP< Stokhos::EpetraOperatorOrthogPoly > getSGPolynomial()
Get SG polynomial.
Teuchos::RCP< EpetraExt::BlockVector > export_residual(const Epetra_Vector &f_overlapped) const
Export parallel residual vector.
Teuchos::RCP< EpetraExt::BlockVector > export_solution(const Epetra_Vector &x_overlapped) const
Export parallel solution vector.
Teuchos::Array< Teuchos::RCP< Teuchos::Array< std::string > > > sg_p_names
SG coefficient parameter names.
Teuchos::RCP< EpetraExt::ModelEvaluator > me
Underlying model evaluator.
virtual void setupPreconditioner(const Teuchos::RCP< Stokhos::SGOperator > &sg_op, const Epetra_Vector &x)=0
Setup preconditioner.
unsigned int num_W_blocks
Number of W stochastic blocks (may be smaller than num_sg_blocks)
OutArgs createOutArgs() const
Create OutArgs.
void set_x_sg_init(const Stokhos::EpetraVectorOrthogPoly &x_sg_in)
Set initial solution polynomial.
InArgs createInArgs() const
Create InArgs.
Teuchos::Array< int > get_p_sg_map_indices() const
Get indices of SG parameters.
Teuchos::RCP< Epetra_Operator > create_DgDp_op(int j, int i) const
Create SG operator representing dg/dp.
Teuchos::RCP< const Epetra_BlockMap > overlapped_stoch_row_map
Overlapped map for stochastic blocks (local map)
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
int num_g
Number of response vectors of underlying model evaluator.
int num_p
Number of parameter vectors of underlying model evaluator.
Teuchos::RCP< Epetra_Import > sg_overlapped_x_importer
Importer from SG to SG-overlapped maps.
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< 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.
Factory for generating stochastic Galerkin preconditioners.
Teuchos::RCP< const Epetra_Import > get_x_sg_importer() const
Return x sg importer.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
Teuchos::RCP< Epetra_Vector > my_x
x pointer for evaluating preconditioner
Teuchos::RCP< Teuchos::ParameterList > params
Algorithmic parameters.
An abstract class to represent a generic stochastic Galerkin preconditioner as an Epetra_Operator...
Teuchos::Array< int > sg_p_index_map
Index map between block-p and p_sg maps.
Teuchos::RCP< const Epetra_Map > f_map
Underlying residual map.
Teuchos::RCP< EpetraMultiVectorOrthogPoly > multiVectorOrthogPoly() const
Get multi vector orthog poly.
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< Epetra_Operator > create_DgDx_op(int j) const
Create SG operator representing dg/dx.
Teuchos::Array< int > get_g_sg_map_indices() const
Get indices of SG responses.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
Teuchos::RCP< const Epetra_Map > sg_x_map
Block SG unknown map.
Teuchos::RCP< const EpetraExt::MultiComm > sg_comm
Parallel SG communicator.
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.