45 #include "Teuchos_Assert.hpp" 46 #include "EpetraExt_BlockUtility.h" 47 #include "EpetraExt_BlockMultiVector.h" 54 const Teuchos::RCP<EpetraExt::ModelEvaluator>& me_,
55 const Teuchos::RCP<const EpetraExt::MultiComm>& mp_comm_,
56 const Teuchos::RCP<const Epetra_Map>& mp_block_map_,
57 const Teuchos::RCP<Teuchos::ParameterList>& params_)
59 num_mp_blocks(mp_block_map_->NumMyElements()),
61 mp_block_map(mp_block_map_),
64 x_map(me->get_x_map()),
65 f_map(me->get_f_map()),
82 if (
x_map != Teuchos::null)
89 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
93 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
104 if (
me->get_x_dot_init() != Teuchos::null) {
121 InArgs me_inargs =
me->createInArgs();
122 num_p = me_inargs.Np();
125 for (
int i=0; i<
num_p; i++) {
126 if (me_inargs.supports(IN_ARG_p_mp, i))
139 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
142 Teuchos::RCP<const Teuchos::Array<std::string> > p_names =
144 if (p_names != Teuchos::null) {
146 Teuchos::rcp(
new Teuchos::Array<std::string>(
num_mp_blocks*(p_names->size())));
147 for (
int j=0;
j<p_names->size();
j++) {
148 std::stringstream ss;
149 ss << (*p_names)[
j] <<
" -- MP Coefficient " << i;
168 OutArgs me_outargs =
me->createOutArgs();
169 num_g = me_outargs.Ng();
172 for (
int i=0; i<
num_g; i++) {
173 if (me_outargs.supports(OUT_ARG_g_mp, i))
184 Teuchos::rcp(EpetraExt::BlockUtility::GenerateBlockMap(
191 Teuchos::RCP<const Epetra_Map>
197 Teuchos::RCP<const Epetra_Map>
203 Teuchos::RCP<const Epetra_Map>
206 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_mp, std::logic_error,
207 "Error! Invalid p map index " << l);
209 return me->get_p_map(l);
211 return mp_p_map[l-num_p];
213 return Teuchos::null;
216 Teuchos::RCP<const Epetra_Map>
219 TEUCHOS_TEST_FOR_EXCEPTION(j < 0 || j >= num_g_mp, std::logic_error,
220 "Error! Invalid g map index " <<
j);
224 Teuchos::RCP<const Teuchos::Array<std::string> >
227 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_mp, std::logic_error,
228 "Error! Invalid p map index " << l);
230 return me->get_p_names(l);
232 return mp_p_names[l-num_p];
234 return Teuchos::null;
237 Teuchos::RCP<const Epetra_Vector>
240 return mp_x_init->getBlockVector();
243 Teuchos::RCP<const Epetra_Vector>
246 return mp_x_dot_init->getBlockVector();
249 Teuchos::RCP<const Epetra_Vector>
252 TEUCHOS_TEST_FOR_EXCEPTION(l < 0 || l >= num_p + num_p_mp, std::logic_error,
253 "Error! Invalid p map index " << l);
255 return me->get_p_init(l);
256 else if (l < num_p + num_p_mp)
257 return mp_p_init[l-num_p]->getBlockVector();
259 return Teuchos::null;
262 Teuchos::RCP<Epetra_Operator>
269 mp_x_map, mp_f_map));
270 my_W->setupOperator(W_mp_blocks);
275 return Teuchos::null;
278 Teuchos::RCP<EpetraExt::ModelEvaluator::Preconditioner>
282 Teuchos::RCP<Teuchos::ParameterList> mpPrecParams =
283 Teuchos::rcp(&(params->sublist(
"MP Preconditioner")),
false);
285 Teuchos::RCP<Epetra_Operator> precOp =
286 mp_prec_factory.
build(mp_comm, num_mp_blocks, x_map, mp_x_map);
287 return Teuchos::rcp(
new EpetraExt::ModelEvaluator::Preconditioner(precOp,
291 return Teuchos::null;
294 Teuchos::RCP<Epetra_Operator>
297 TEUCHOS_TEST_FOR_EXCEPTION(j < 0 || j >= num_g_mp || !supports_x,
299 "Error: dg/dx index " <<
j <<
" is not supported!");
301 int jj = mp_g_index_map[
j];
302 Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(jj);
303 Teuchos::RCP< Stokhos::ProductEpetraOperator > mp_blocks;
304 OutArgs me_outargs = me->createOutArgs();
305 DerivativeSupport ds = me_outargs.supports(OUT_ARG_DgDx_mp, jj);
306 if (ds.supports(DERIV_LINEAR_OP)) {
309 mp_block_map, x_map, g_map, mp_g_map[
j], mp_comm));
310 for (
unsigned int i=0; i<num_mp_blocks; i++)
311 mp_blocks->setCoeffPtr(i, me->create_DgDx_op(jj));
313 else if (ds.supports(DERIV_MV_BY_COL)) {
314 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
316 mp_block_map, g_map, mp_g_map[
j],
317 mp_comm, x_map->NumMyElements()));
320 mp_mv_blocks,
false));
322 else if (ds.supports(DERIV_TRANS_MV_BY_ROW)) {
323 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
325 mp_block_map, x_map, mp_x_map,
326 mp_comm, g_map->NumMyElements()));
329 mp_mv_blocks,
true));
332 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
333 "Error! me_outargs.supports(OUT_ARG_DgDx_mp, " <<
j 334 <<
").none() is true!");
336 Teuchos::RCP<Stokhos::BlockDiagonalOperator> dgdx_mp =
338 mp_comm, num_mp_blocks, x_map, g_map,
339 mp_x_map, mp_g_map[
j]));
340 dgdx_mp->setupOperator(mp_blocks);
344 Teuchos::RCP<Epetra_Operator>
347 TEUCHOS_TEST_FOR_EXCEPTION(j < 0 || j >= num_g_mp || !supports_x,
349 "Error: dg/dx_dot index " <<
j <<
" is not supported!");
351 int jj = mp_g_index_map[
j];
352 Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(jj);
353 Teuchos::RCP< Stokhos::ProductEpetraOperator > mp_blocks;
354 OutArgs me_outargs = me->createOutArgs();
355 DerivativeSupport ds = me_outargs.supports(OUT_ARG_DgDx_dot_mp, jj);
356 if (ds.supports(DERIV_LINEAR_OP)) {
359 mp_block_map, x_map, g_map, mp_g_map[
j],
361 for (
unsigned int i=0; i<num_mp_blocks; i++)
362 mp_blocks->setCoeffPtr(i, me->create_DgDx_dot_op(jj));
364 else if (ds.supports(DERIV_MV_BY_COL)) {
365 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
367 mp_block_map, g_map, mp_g_map[
j],
368 mp_comm, x_map->NumMyElements()));
371 mp_mv_blocks,
false));
373 else if (ds.supports(DERIV_TRANS_MV_BY_ROW)) {
374 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
376 mp_block_map, x_map, mp_x_map,
377 mp_comm, g_map->NumMyElements()));
380 mp_mv_blocks,
true));
383 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
384 "Error! me_outargs.supports(OUT_ARG_DgDx_dot_mp, " 385 <<
j <<
").none() is true!");
387 Teuchos::RCP<Stokhos::BlockDiagonalOperator> dgdx_dot_mp =
389 mp_comm, num_mp_blocks, x_map, g_map,
390 mp_x_map, mp_g_map[
j]));
391 dgdx_dot_mp->setupOperator(mp_blocks);
395 Teuchos::RCP<Epetra_Operator>
398 TEUCHOS_TEST_FOR_EXCEPTION(
399 j < 0 || j >= num_g_mp || i < 0 || i >= num_p+num_p_mp,
401 "Error: dg/dp index " <<
j <<
"," << i <<
" is not supported!");
403 OutArgs me_outargs = me->createOutArgs();
404 int jj = mp_g_index_map[
j];
405 Teuchos::RCP<const Epetra_Map> g_map = me->get_g_map(jj);
407 if (me_outargs.supports(OUT_ARG_DgDp_mp,jj,i).supports(DERIV_LINEAR_OP)) {
408 Teuchos::RCP<Stokhos::ProductEpetraOperator> mp_blocks =
410 mp_block_map, me->get_p_map(i), g_map,
411 mp_g_map[
j], mp_comm));
412 for (
unsigned int l=0; l<num_mp_blocks; l++)
413 mp_blocks->setCoeffPtr(l, me->create_DgDp_op(i,
j));
417 TEUCHOS_TEST_FOR_EXCEPTION(
418 true, std::logic_error,
419 "Error: Underlying model evaluator must support DERIV_LINER_OP " <<
420 "to create operator for dg/dp index " <<
j <<
"," << i <<
"!");
423 int ii = mp_p_index_map[i-num_p];
424 Teuchos::RCP<const Epetra_Map> p_map = me->get_p_map(ii);
425 Teuchos::RCP< Stokhos::ProductEpetraOperator> mp_blocks;
426 DerivativeSupport ds = me_outargs.supports(OUT_ARG_DgDp_mp,jj,ii);
427 if (ds.supports(DERIV_LINEAR_OP)) {
430 mp_block_map, p_map, g_map,
431 mp_g_map[
j], mp_comm));
432 for (
unsigned int l=0; l<num_mp_blocks; l++)
433 mp_blocks->setCoeffPtr(l, me->create_DgDp_op(jj,ii));
435 else if (ds.supports(DERIV_MV_BY_COL)) {
436 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
438 mp_block_map, g_map, mp_g_map[
j],
439 mp_comm, p_map->NumMyElements()));
442 mp_mv_blocks,
false));
444 else if (ds.supports(DERIV_TRANS_MV_BY_ROW)) {
445 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
447 mp_block_map, p_map, mp_p_map[i-num_p],
448 mp_comm, g_map->NumMyElements()));
451 mp_mv_blocks,
true));
454 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
455 "Error! me_outargs.supports(OUT_ARG_DgDp_mp, " << jj
456 <<
"," << ii <<
").none() is true!");
458 Teuchos::RCP<Stokhos::BlockDiagonalOperator> dgdp_mp =
460 mp_comm, num_mp_blocks, p_map, g_map,
461 mp_p_map[i-num_p], mp_g_map[
j]));
462 dgdp_mp->setupOperator(mp_blocks);
466 return Teuchos::null;
469 Teuchos::RCP<Epetra_Operator>
472 TEUCHOS_TEST_FOR_EXCEPTION(i < 0 || i >= num_p+num_p_mp,
474 "Error: df/dp index " << i <<
" is not supported!");
476 OutArgs me_outargs = me->createOutArgs();
478 if (me_outargs.supports(OUT_ARG_DfDp_mp,i).supports(DERIV_LINEAR_OP)) {
479 Teuchos::RCP<Stokhos::ProductEpetraOperator> mp_blocks =
481 mp_block_map, me->get_p_map(i), me->get_f_map(),
483 for (
unsigned int l=0; l<num_mp_blocks; l++)
484 mp_blocks->setCoeffPtr(l, me->create_DfDp_op(i));
488 TEUCHOS_TEST_FOR_EXCEPTION(
489 true, std::logic_error,
490 "Error: Underlying model evaluator must support DERIV_LINER_OP " <<
491 "to create operator for df/dp index " << i <<
"!");
494 int ii = mp_p_index_map[i-num_p];
495 Teuchos::RCP<const Epetra_Map> p_map = me->get_p_map(ii);
496 Teuchos::RCP<Stokhos::ProductEpetraOperator> mp_blocks;
497 DerivativeSupport ds = me_outargs.supports(OUT_ARG_DfDp_mp,ii);
498 if (ds.supports(DERIV_LINEAR_OP)) {
501 mp_block_map, me->get_p_map(ii), me->get_f_map(),
503 for (
unsigned int l=0; l<num_mp_blocks; l++)
504 mp_blocks->setCoeffPtr(l, me->create_DfDp_op(ii));
506 else if (ds.supports(DERIV_MV_BY_COL)) {
507 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
509 mp_block_map, f_map, mp_f_map,
510 mp_comm, p_map->NumMyElements()));
513 mp_mv_blocks,
false));
515 else if (ds.supports(DERIV_TRANS_MV_BY_ROW)) {
516 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_mv_blocks =
518 mp_block_map, p_map, mp_p_map[i-num_p],
519 mp_comm, f_map->NumMyElements()));
522 mp_mv_blocks,
true));
525 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
526 "Error! me_outargs.supports(OUT_ARG_DfDp_mp, " << ii
527 <<
").none() is true!");
530 Teuchos::RCP<Stokhos::BlockDiagonalOperator> dfdp_mp =
532 mp_comm, num_mp_blocks,
533 p_map, f_map, mp_p_map[i-num_p], mp_f_map));
534 dfdp_mp->setupOperator(mp_blocks);
538 return Teuchos::null;
541 EpetraExt::ModelEvaluator::InArgs
545 InArgs me_inargs = me->createInArgs();
547 inArgs.setModelEvalDescription(this->description());
548 inArgs.set_Np(num_p + num_p_mp);
549 inArgs.setSupports(IN_ARG_x_dot, me_inargs.supports(IN_ARG_x_dot));
550 inArgs.setSupports(IN_ARG_x, me_inargs.supports(IN_ARG_x));
551 inArgs.setSupports(IN_ARG_t, me_inargs.supports(IN_ARG_t));
552 inArgs.setSupports(IN_ARG_alpha, me_inargs.supports(IN_ARG_alpha));
553 inArgs.setSupports(IN_ARG_beta, me_inargs.supports(IN_ARG_beta));
558 EpetraExt::ModelEvaluator::OutArgs
561 OutArgsSetup outArgs;
562 OutArgs me_outargs = me->createOutArgs();
564 outArgs.setModelEvalDescription(this->description());
565 outArgs.set_Np_Ng(num_p + num_p_mp, num_g_mp);
566 outArgs.setSupports(OUT_ARG_f, me_outargs.supports(OUT_ARG_f));
567 outArgs.setSupports(OUT_ARG_W, me_outargs.supports(OUT_ARG_W));
568 outArgs.setSupports(OUT_ARG_WPrec, me_outargs.supports(OUT_ARG_W));
569 for (
int j=0;
j<num_p;
j++)
570 outArgs.setSupports(OUT_ARG_DfDp,
j,
571 me_outargs.supports(OUT_ARG_DfDp,
j));
572 for (
int j=0;
j<num_p_mp;
j++)
573 if (!me_outargs.supports(OUT_ARG_DfDp_mp, mp_p_index_map[
j]).none())
574 outArgs.setSupports(OUT_ARG_DfDp,
j+num_p, DERIV_LINEAR_OP);
575 for (
int i=0; i<num_g_mp; i++) {
576 int ii = mp_g_index_map[i];
577 if (!me_outargs.supports(OUT_ARG_DgDx_dot_mp, ii).none())
578 outArgs.setSupports(OUT_ARG_DgDx_dot, i, DERIV_LINEAR_OP);
579 if (!me_outargs.supports(OUT_ARG_DgDx_mp, ii).none())
580 outArgs.setSupports(OUT_ARG_DgDx, i, DERIV_LINEAR_OP);
581 for (
int j=0;
j<num_p;
j++)
582 outArgs.setSupports(OUT_ARG_DgDp, i,
j,
583 me_outargs.supports(OUT_ARG_DgDp_mp, ii,
j));
584 for (
int j=0;
j<num_p_mp;
j++)
585 if (!me_outargs.supports(OUT_ARG_DgDp_mp, ii, mp_p_index_map[
j]).none())
586 outArgs.setSupports(OUT_ARG_DgDp, i,
j+num_p, DERIV_LINEAR_OP);
594 const OutArgs& outArgs)
const 597 Teuchos::RCP<const Epetra_Vector>
x;
598 if (inArgs.supports(IN_ARG_x)) {
600 if (
x != Teuchos::null)
603 Teuchos::RCP<const Epetra_Vector> x_dot;
604 if (inArgs.supports(IN_ARG_x_dot))
605 x_dot = inArgs.get_x_dot();
608 EpetraExt::ModelEvaluator::Evaluation<Epetra_Vector> f_out;
609 if (outArgs.supports(OUT_ARG_f))
610 f_out = outArgs.get_f();
611 Teuchos::RCP<Epetra_Operator> W_out;
612 if (outArgs.supports(OUT_ARG_W))
613 W_out = outArgs.get_W();
614 Teuchos::RCP<Epetra_Operator> WPrec_out;
615 if (outArgs.supports(OUT_ARG_WPrec))
616 WPrec_out = outArgs.get_WPrec();
620 bool eval_prec = (W_out == Teuchos::null && WPrec_out != Teuchos::null);
626 Teuchos::RCP<Stokhos::MPPreconditioner> W_prec =
631 bool done = (f_out == Teuchos::null);
632 for (
int i=0; i<outArgs.Ng(); i++) {
633 done = done && (outArgs.get_g(i) == Teuchos::null);
634 for (
int j=0;
j<outArgs.Np();
j++)
635 if (!outArgs.supports(OUT_ARG_DgDp, i,
j).none())
636 done = done && (outArgs.get_DgDp(i,
j).isEmpty());
643 InArgs me_inargs = me->createInArgs();
644 if (
x != Teuchos::null) {
645 Teuchos::RCP<Stokhos::ProductEpetraVector> x_mp =
646 create_x_mp(
View,
x.get());
647 me_inargs.set_x_mp(x_mp);
649 if (x_dot != Teuchos::null) {
650 Teuchos::RCP<Stokhos::ProductEpetraVector> x_dot_mp =
651 create_x_mp(
View, x_dot.get());
652 me_inargs.set_x_dot_mp(x_dot_mp);
654 if (me_inargs.supports(IN_ARG_alpha))
655 me_inargs.set_alpha(inArgs.get_alpha());
656 if (me_inargs.supports(IN_ARG_beta))
657 me_inargs.set_beta(inArgs.get_beta());
658 if (me_inargs.supports(IN_ARG_t))
659 me_inargs.set_t(inArgs.get_t());
662 for (
int i=0; i<num_p; i++)
663 me_inargs.set_p(i, inArgs.get_p(i));
664 for (
int i=0; i<num_p_mp; i++) {
665 Teuchos::RCP<const Epetra_Vector> p = inArgs.get_p(i+num_p);
669 if (p == Teuchos::null)
670 p = mp_p_init[i]->getBlockVector();
673 Teuchos::RCP<Stokhos::ProductEpetraVector> p_mp =
674 create_p_mp(mp_p_index_map[i],
View, p.get());
675 me_inargs.set_p_mp(mp_p_index_map[i], p_mp);
679 OutArgs me_outargs = me->createOutArgs();
682 if (f_out != Teuchos::null) {
683 Teuchos::RCP<Stokhos::ProductEpetraVector> f_mp =
684 create_f_mp(
View, f_out.get());
685 me_outargs.set_f_mp(f_mp);
689 if (W_out != Teuchos::null && !eval_prec)
690 me_outargs.set_W_mp(W_mp_blocks);
693 for (
int i=0; i<num_p; i++) {
694 if (!outArgs.supports(OUT_ARG_DfDp, i).none()) {
695 Derivative dfdp = outArgs.get_DfDp(i);
696 if (dfdp.getMultiVector() != Teuchos::null) {
697 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dfdp_mp;
698 if (dfdp.getMultiVectorOrientation() == DERIV_MV_BY_COL)
701 mp_block_map, me->get_f_map(), mp_f_map, mp_comm,
702 View, *(dfdp.getMultiVector())));
703 else if (dfdp.getMultiVectorOrientation() == DERIV_TRANS_MV_BY_ROW)
706 mp_block_map, me->get_p_map(i), mp_p_map[i], mp_comm,
707 View, *(dfdp.getMultiVector())));
708 me_outargs.set_DfDp_mp(i,
709 MPDerivative(dfdp_mp,
710 dfdp.getMultiVectorOrientation()));
712 else if (dfdp.getLinearOp() != Teuchos::null) {
713 Teuchos::RCP<Stokhos::ProductEpetraOperator> dfdp_mp =
715 me_outargs.set_DfDp_mp(i, MPDerivative(dfdp_mp));
721 for (
int i=0; i<num_p_mp; i++) {
722 if (!outArgs.supports(OUT_ARG_DfDp, i+num_p).none()) {
723 Derivative dfdp = outArgs.get_DfDp(i+num_p);
724 if (dfdp.getLinearOp() != Teuchos::null) {
725 Teuchos::RCP<Stokhos::BlockDiagonalOperator> dfdp_op =
727 Teuchos::RCP<Stokhos::ProductEpetraOperator> dfdp_op_mp =
729 int ii = mp_p_index_map[i];
730 if (me_outargs.supports(OUT_ARG_DfDp_mp,ii).supports(DERIV_LINEAR_OP))
731 me_outargs.set_DfDp_mp(ii, MPDerivative(dfdp_op_mp));
733 Teuchos::RCP<Stokhos::ProductEpetraMultiVectorOperator> mp_mv_op =
735 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dfdp_mp =
737 if (me_outargs.supports(OUT_ARG_DfDp_mp,ii).supports(DERIV_MV_BY_COL))
738 me_outargs.set_DfDp_mp(
739 ii, MPDerivative(dfdp_mp, DERIV_MV_BY_COL));
741 me_outargs.set_DfDp_mp(
742 ii, MPDerivative(dfdp_mp, DERIV_TRANS_MV_BY_ROW));
745 TEUCHOS_TEST_FOR_EXCEPTION(
746 dfdp.getLinearOp() == Teuchos::null && dfdp.isEmpty() ==
false,
748 "Error! Stokhos::MPModelEvaluator::evalModel: " <<
749 "Operator form of df/dp(" << i+num_p <<
") is required!");
754 for (
int i=0; i<num_g_mp; i++) {
755 int ii = mp_g_index_map[i];
758 Teuchos::RCP<Epetra_Vector>
g = outArgs.get_g(i);
759 if (
g != Teuchos::null) {
760 Teuchos::RCP<Stokhos::ProductEpetraVector> g_mp =
761 create_g_mp(ii,
View,
g.get());
762 me_outargs.set_g_mp(ii, g_mp);
766 if (outArgs.supports(OUT_ARG_DgDx_dot, i).supports(DERIV_LINEAR_OP)) {
767 Derivative dgdx_dot = outArgs.get_DgDx_dot(i);
768 if (dgdx_dot.getLinearOp() != Teuchos::null) {
769 Teuchos::RCP<Stokhos::BlockDiagonalOperator> op =
771 dgdx_dot.getLinearOp(),
true);
772 Teuchos::RCP<Stokhos::ProductEpetraOperator> mp_blocks =
774 if (me_outargs.supports(OUT_ARG_DgDx, ii).supports(DERIV_LINEAR_OP))
775 me_outargs.set_DgDx_dot_mp(ii, mp_blocks);
777 Teuchos::RCP<Stokhos::ProductEpetraMultiVectorOperator> mp_mv_op =
779 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdx_dot_mp =
781 if (me_outargs.supports(OUT_ARG_DgDx_dot_mp, ii).supports(DERIV_MV_BY_COL))
782 me_outargs.set_DgDx_dot_mp(ii, MPDerivative(dgdx_dot_mp,
785 me_outargs.set_DgDx_dot_mp(ii, MPDerivative(dgdx_dot_mp,
786 DERIV_TRANS_MV_BY_ROW));
789 TEUCHOS_TEST_FOR_EXCEPTION(dgdx_dot.getLinearOp() == Teuchos::null &&
790 dgdx_dot.isEmpty() ==
false,
792 "Error! Stokhos::MPModelEvaluator::evalModel: " <<
793 "Operator form of dg/dxdot is required!");
797 if (outArgs.supports(OUT_ARG_DgDx, i).supports(DERIV_LINEAR_OP)) {
798 Derivative dgdx = outArgs.get_DgDx(i);
799 if (dgdx.getLinearOp() != Teuchos::null) {
800 Teuchos::RCP<Stokhos::BlockDiagonalOperator> op =
802 dgdx.getLinearOp(),
true);
803 Teuchos::RCP<Stokhos::ProductEpetraOperator> mp_blocks =
805 if (me_outargs.supports(OUT_ARG_DgDx, ii).supports(DERIV_LINEAR_OP))
806 me_outargs.set_DgDx_mp(ii, mp_blocks);
808 Teuchos::RCP<Stokhos::ProductEpetraMultiVectorOperator> mp_mv_op =
810 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdx_mp =
812 if (me_outargs.supports(OUT_ARG_DgDx_mp, ii).supports(DERIV_MV_BY_COL))
813 me_outargs.set_DgDx_mp(ii, MPDerivative(dgdx_mp,
816 me_outargs.set_DgDx_mp(ii, MPDerivative(dgdx_mp,
817 DERIV_TRANS_MV_BY_ROW));
820 TEUCHOS_TEST_FOR_EXCEPTION(dgdx.getLinearOp() == Teuchos::null &&
821 dgdx.isEmpty() ==
false,
823 "Error! Stokhos::MPModelEvaluator::evalModel: " <<
824 "Operator form of dg/dxdot is required!");
828 for (
int j=0;
j<num_p;
j++) {
829 if (!outArgs.supports(OUT_ARG_DgDp, i,
j).none()) {
830 Derivative dgdp = outArgs.get_DgDp(i,
j);
831 if (dgdp.getMultiVector() != Teuchos::null) {
832 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdp_mp;
833 if (dgdp.getMultiVectorOrientation() == DERIV_MV_BY_COL)
836 mp_block_map, me->get_g_map(ii), mp_g_map[i],
837 mp_comm,
View, *(dgdp.getMultiVector())));
838 else if (dgdp.getMultiVectorOrientation() == DERIV_TRANS_MV_BY_ROW)
841 mp_block_map, me->get_p_map(
j), mp_p_map[
j],
842 mp_comm,
View, *(dgdp.getMultiVector())));
843 me_outargs.set_DgDp_mp(
844 ii,
j, MPDerivative(dgdp_mp, dgdp.getMultiVectorOrientation()));
846 else if (dgdp.getLinearOp() != Teuchos::null) {
847 Teuchos::RCP<Stokhos::ProductEpetraOperator> dgdp_mp =
849 me_outargs.set_DgDp_mp(ii,
j, MPDerivative(dgdp_mp));
855 for (
int j=0;
j<num_p_mp;
j++) {
856 if (!outArgs.supports(OUT_ARG_DgDp, i,
j+num_p).none()) {
857 Derivative dgdp = outArgs.get_DgDp(i,
j+num_p);
858 if (dgdp.getLinearOp() != Teuchos::null) {
859 Teuchos::RCP<Stokhos::BlockDiagonalOperator> dgdp_op =
861 Teuchos::RCP<Stokhos::ProductEpetraOperator> dgdp_op_mp =
863 int jj = mp_p_index_map[
j];
864 if (me_outargs.supports(OUT_ARG_DgDp_mp,ii,jj).supports(DERIV_LINEAR_OP))
865 me_outargs.set_DgDp_mp(ii, jj, MPDerivative(dgdp_op_mp));
867 Teuchos::RCP<Stokhos::ProductEpetraMultiVectorOperator> mp_mv_op =
869 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdp_mp =
871 if (me_outargs.supports(OUT_ARG_DgDp_mp,ii,jj).supports(DERIV_MV_BY_COL))
872 me_outargs.set_DgDp_mp(
873 ii, jj, MPDerivative(dgdp_mp, DERIV_MV_BY_COL));
875 me_outargs.set_DgDp_mp(
876 ii, jj, MPDerivative(dgdp_mp, DERIV_TRANS_MV_BY_ROW));
879 TEUCHOS_TEST_FOR_EXCEPTION(
880 dgdp.getLinearOp() == Teuchos::null && dgdp.isEmpty() ==
false,
882 "Error! Stokhos::MPModelEvaluator::evalModel: " <<
883 "Operator form of dg/dp(" << i <<
"," <<
j+num_p <<
") is required!");
890 me->evalModel(me_inargs, me_outargs);
893 if (W_out != Teuchos::null && !eval_prec) {
894 Teuchos::RCP<Epetra_Operator> W;
895 if (W_out != Teuchos::null)
899 Teuchos::RCP<Stokhos::BlockDiagonalOperator> W_mp =
903 if (WPrec_out != Teuchos::null) {
904 Teuchos::RCP<Stokhos::MPPreconditioner> W_prec =
915 *mp_x_init = x_mp_in;
922 *mp_x_dot_init = x_dot_mp_in;
925 Teuchos::RCP<const Stokhos::ProductEpetraVector>
931 Teuchos::RCP<const Stokhos::ProductEpetraVector>
934 return mp_x_dot_init;
941 Teuchos::Array<int>::iterator it = std::find(mp_p_index_map.begin(),
942 mp_p_index_map.end(),
944 TEUCHOS_TEST_FOR_EXCEPTION(it == mp_p_index_map.end(), std::logic_error,
945 "Error! Invalid p map index " << i);
946 int ii = it - mp_p_index_map.begin();
947 *mp_p_init[ii] = p_mp_in;
950 Teuchos::RCP<const Stokhos::ProductEpetraVector>
953 Teuchos::Array<int>::const_iterator it = std::find(mp_p_index_map.begin(),
954 mp_p_index_map.end(),
956 TEUCHOS_TEST_FOR_EXCEPTION(it == mp_p_index_map.end(), std::logic_error,
957 "Error! Invalid p map index " << l);
958 int ll = it - mp_p_index_map.begin();
959 return mp_p_init[ll];
965 return mp_p_index_map;
971 return mp_g_index_map;
974 Teuchos::Array< Teuchos::RCP<const Epetra_Map> >
977 Teuchos::Array< Teuchos::RCP<const Epetra_Map> > base_maps(num_g);
978 for (
int i=0; i<num_g; i++)
979 base_maps[i] = me->get_g_map(i);
983 Teuchos::RCP<Stokhos::ProductEpetraVector>
987 Teuchos::RCP<Stokhos::ProductEpetraVector> mp_x;
990 mp_block_map, x_map, mp_x_map, mp_comm));
993 mp_block_map, x_map, mp_x_map, mp_comm,
998 Teuchos::RCP<Stokhos::ProductEpetraMultiVector>
1002 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_x;
1005 mp_block_map, x_map, mp_x_map, mp_comm,
1009 mp_block_map, x_map, mp_x_map, mp_comm,
1014 Teuchos::RCP<Stokhos::ProductEpetraVector>
1018 Teuchos::RCP<Stokhos::ProductEpetraVector> mp_p;
1019 Teuchos::Array<int>::const_iterator it = std::find(mp_p_index_map.begin(),
1020 mp_p_index_map.end(),
1022 TEUCHOS_TEST_FOR_EXCEPTION(it == mp_p_index_map.end(), std::logic_error,
1023 "Error! Invalid p map index " << l);
1024 int ll = it - mp_p_index_map.begin();
1027 mp_block_map, me->get_p_map(l),
1028 mp_p_map[ll], mp_comm));
1031 mp_block_map, me->get_p_map(l),
1032 mp_p_map[ll], mp_comm, CV, *v));
1036 Teuchos::RCP<Stokhos::ProductEpetraMultiVector>
1041 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_p;
1042 Teuchos::Array<int>::const_iterator it = std::find(mp_p_index_map.begin(),
1043 mp_p_index_map.end(),
1045 TEUCHOS_TEST_FOR_EXCEPTION(it == mp_p_index_map.end(), std::logic_error,
1046 "Error! Invalid p map index " << l);
1047 int ll = it - mp_p_index_map.begin();
1050 mp_block_map, me->get_p_map(l),
1051 mp_p_map[ll], mp_comm, num_vecs));
1054 mp_block_map, me->get_p_map(l),
1055 mp_p_map[ll], mp_comm, CV, *v));
1059 Teuchos::RCP<Stokhos::ProductEpetraVector>
1063 Teuchos::RCP<Stokhos::ProductEpetraVector> mp_f;
1066 mp_block_map, f_map, mp_f_map, mp_comm));
1069 mp_block_map, f_map, mp_f_map, mp_comm,
1074 Teuchos::RCP<Stokhos::ProductEpetraMultiVector>
1080 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_f;
1083 mp_block_map, f_map, mp_f_map, mp_comm,
1087 mp_block_map, f_map, mp_f_map, mp_comm,
1092 Teuchos::RCP<Stokhos::ProductEpetraVector>
1096 Teuchos::RCP<Stokhos::ProductEpetraVector> mp_g;
1097 Teuchos::Array<int>::const_iterator it = std::find(mp_g_index_map.begin(),
1098 mp_g_index_map.end(),
1100 TEUCHOS_TEST_FOR_EXCEPTION(it == mp_g_index_map.end(), std::logic_error,
1101 "Error! Invalid g map index " << l);
1102 int ll = it - mp_g_index_map.begin();
1107 mp_g_map[ll], mp_comm));
1112 mp_g_map[ll], mp_comm, CV, *v));
1116 Teuchos::RCP<Stokhos::ProductEpetraMultiVector>
1121 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> mp_g;
1122 Teuchos::Array<int>::const_iterator it = std::find(mp_g_index_map.begin(),
1123 mp_g_index_map.end(),
1125 TEUCHOS_TEST_FOR_EXCEPTION(it == mp_g_index_map.end(), std::logic_error,
1126 "Error! Invalid g map index " << l);
1127 int ll = it - mp_g_index_map.begin();
1132 mp_g_map[ll], mp_comm, num_vecs));
1137 mp_g_map[ll], mp_comm, CV, *v));
Teuchos::RCP< const Epetra_Map > get_p_map(int l) const
Return parameter vector map.
Factory for generating stochastic Galerkin preconditioners.
Teuchos::RCP< const Epetra_Map > mp_x_map
Block MP unknown map.
Teuchos::Array< int > get_p_mp_map_indices() const
Get indices of MP parameters.
A container class for products of Epetra_Vector's.
bool supports_x
Whether we support x (and thus f and W)
virtual Teuchos::RCP< Stokhos::MPPreconditioner > build(const Teuchos::RCP< const EpetraExt::MultiComm > &mp_comm, int num_mp_blocks, const Teuchos::RCP< const Epetra_Map > &base_map, const Teuchos::RCP< const Epetra_Map > &mp_map)
Build preconditioner operator.
void set_x_mp_init(const Stokhos::ProductEpetraVector &x_mp_in)
Set initial multi-point solution.
int num_p_mp
Number of multi-point parameter vectors.
int num_p
Number of parameter vectors of underlying model evaluator.
Teuchos::RCP< const Epetra_Map > get_x_map() const
Return solution vector map.
OutArgs createOutArgs() const
Create OutArgs.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > mp_g_map
Block MP response map.
void set_x_dot_mp_init(const Stokhos::ProductEpetraVector &x_dot_mp_in)
Teuchos::Array< int > mp_p_index_map
Index map between block-p and p_mp maps.
Teuchos::RCP< Epetra_Operator > create_DgDx_op(int j) const
Create MP operator representing dg/dx.
virtual void setupOperator(const Teuchos::RCP< Stokhos::ProductEpetraOperator > &ops)
Setup operator.
Teuchos::RCP< const Epetra_Map > get_g_map(int l) const
Return response map.
InArgs createInArgs() const
Create InArgs.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Teuchos::RCP< Epetra_Operator > create_DgDx_dot_op(int j) const
Create MP operator representing dg/dxdot.
Teuchos::RCP< Stokhos::ProductEpetraMultiVector > create_g_mv_mp(int l, int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create multi-point multi-vector using g map.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > get_g_mp_base_maps() const
Get base maps of MP responses.
Teuchos::RCP< Stokhos::ProductEpetraVector > create_x_mp(Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create multi-point vector using x map and owned mp map.
A container class for products of Epetra_Vector's.
Teuchos::RCP< EpetraExt::ModelEvaluator > me
Underlying model evaluator.
Teuchos::RCP< Stokhos::ProductEpetraOperator > W_mp_blocks
W multi-point components.
Teuchos::RCP< Stokhos::ProductEpetraVector > create_f_mp(Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create multi-point vector using f map.
Teuchos::RCP< const Epetra_Map > get_f_map() const
Return residual vector map.
An Epetra operator representing the block stochastic Galerkin operator.
Teuchos::RCP< Stokhos::ProductEpetraVector > mp_x_dot_init
Teuchos::RCP< const Epetra_Map > f_map
Underlying residual map.
Teuchos::RCP< const Stokhos::ProductEpetraVector > get_x_dot_mp_init() const
A container class for products of Epetra_Vector's.
A container class storing products of Epetra_MultiVector's.
Teuchos::RCP< const Stokhos::ProductEpetraVector > get_x_mp_init() const
Return initial SG x.
Teuchos::RCP< const Epetra_Map > mp_block_map
Map for layout of parallel MP blocks.
MPModelEvaluator(const Teuchos::RCP< EpetraExt::ModelEvaluator > &me, const Teuchos::RCP< const EpetraExt::MultiComm > &mp_comm, const Teuchos::RCP< const Epetra_Map > &mp_block_map, const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Teuchos::Array< Teuchos::RCP< ProductEpetraVector > > mp_p_init
MP initial p.
Teuchos::RCP< Epetra_Vector > my_x
x pointer for evaluating preconditioner
Teuchos::Array< int > mp_g_index_map
Index map between block-g and g_mp maps.
Teuchos::RCP< const Epetra_Map > mp_f_map
Block MP residual map.
Teuchos::RCP< ProductEpetraMultiVector > productMultiVector() const
Get product multi vector.
Teuchos::RCP< const Epetra_Map > x_map
Underlying unknown map.
Teuchos::RCP< const Epetra_Vector > get_p_init(int l) const
Return initial parameters.
Teuchos::RCP< Stokhos::ProductEpetraMultiVector > create_x_mv_mp(int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create multi-point vector using x map.
virtual Teuchos::RCP< Stokhos::ProductEpetraOperator > getMPOps()
Get multi-point ops.
int num_g
Number of response vectors of underlying model evaluator.
Teuchos::RCP< Stokhos::ProductEpetraMultiVector > create_p_mv_mp(int l, int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create multi-point vector using p map.
int num_g_mp
Number of multi-point response vectors.
Teuchos::RCP< Epetra_Operator > create_W() const
Create W = alpha*M + beta*J matrix.
Teuchos::Array< int > get_g_mp_map_indices() const
Get indices of MP responses.
virtual void setupPreconditioner(const Teuchos::RCP< Stokhos::BlockDiagonalOperator > &mp_op, const Epetra_Vector &x)=0
Setup preconditioner.
An abstract class to represent a generic stochastic Galerkin preconditioner as an Epetra_Operator...
Teuchos::RCP< const Epetra_Vector > get_x_init() const
Return initial solution.
Teuchos::RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
Return array of parameter names.
Teuchos::RCP< Stokhos::ProductEpetraVector > create_g_mp(int l, Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create multi-point vector using g map.
Teuchos::RCP< Stokhos::ProductEpetraMultiVector > create_f_mv_mp(int num_vecs, Epetra_DataAccess CV=Copy, const Epetra_MultiVector *v=NULL) const
Create multi-point multi-vector using f map.
Teuchos::RCP< const EpetraExt::MultiComm > mp_comm
Parallel MP communicator.
Teuchos::RCP< Stokhos::ProductEpetraVector > create_p_mp(int l, Epetra_DataAccess CV=Copy, const Epetra_Vector *v=NULL) const
Create multi-point vector using p map.
Teuchos::Array< Teuchos::RCP< const Epetra_Map > > mp_p_map
Block MP parameter map.
unsigned int num_mp_blocks
Number of blocks.
void set_p_mp_init(int i, const Stokhos::ProductEpetraVector &p_mp_in)
Set initial multi-point parameter.
Teuchos::RCP< Epetra_Operator > create_DfDp_op(int i) const
Create MP operator representing df/dp.
void evalModel(const InArgs &inArgs, const OutArgs &outArgs) const
Evaluate model on InArgs.
Teuchos::RCP< Stokhos::ProductEpetraVector > mp_x_init
MP initial x.
Teuchos::RCP< Epetra_Operator > create_DgDp_op(int j, int i) const
Create MP operator representing dg/dp.
Teuchos::RCP< const Stokhos::ProductEpetraVector > get_p_mp_init(int l) const
Return initial SG parameters.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
Teuchos::RCP< const Epetra_Vector > get_x_dot_init() const
Teuchos::RCP< EpetraExt::ModelEvaluator::Preconditioner > create_WPrec() const
Create preconditioner operator.
Teuchos::Array< Teuchos::RCP< Teuchos::Array< std::string > > > mp_p_names
MP coefficient parameter names.