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;
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 =
627 Teuchos::rcp_dynamic_cast<Stokhos::MPPreconditioner>(WPrec_out,
true);
628 W_prec->setupPreconditioner(my_W, *my_x);
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 =
714 Teuchos::rcp_dynamic_cast<Stokhos::ProductEpetraOperator>(dfdp.getLinearOp(),
true);
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 =
726 Teuchos::rcp_dynamic_cast<Stokhos::BlockDiagonalOperator>(dfdp.getLinearOp(),
true);
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 =
734 Teuchos::rcp_dynamic_cast<Stokhos::ProductEpetraMultiVectorOperator>(dfdp_op_mp,
true);
735 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dfdp_mp =
736 mp_mv_op->productMultiVector();
737 if (me_outargs.supports(OUT_ARG_DfDp_mp,ii).supports(DERIV_MV_BY_COL))
738 me_outargs.set_DfDp_mp(
739 ii, MPDerivative(dfdp_mp, DERIV_MV_BY_COL));
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 =
770 Teuchos::rcp_dynamic_cast<Stokhos::BlockDiagonalOperator>(
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 =
778 Teuchos::rcp_dynamic_cast<Stokhos::ProductEpetraMultiVectorOperator>(mp_blocks,
true);
779 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdx_dot_mp =
780 mp_mv_op->productMultiVector();
781 if (me_outargs.supports(OUT_ARG_DgDx_dot_mp, ii).supports(DERIV_MV_BY_COL))
782 me_outargs.set_DgDx_dot_mp(ii, MPDerivative(dgdx_dot_mp,
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 =
801 Teuchos::rcp_dynamic_cast<Stokhos::BlockDiagonalOperator>(
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 =
809 Teuchos::rcp_dynamic_cast<Stokhos::ProductEpetraMultiVectorOperator>(mp_blocks,
true);
810 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdx_mp =
811 mp_mv_op->productMultiVector();
812 if (me_outargs.supports(OUT_ARG_DgDx_mp, ii).supports(DERIV_MV_BY_COL))
813 me_outargs.set_DgDx_mp(ii, MPDerivative(dgdx_mp,
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 =
848 Teuchos::rcp_dynamic_cast<Stokhos::ProductEpetraOperator>(dgdp.getLinearOp(),
true);
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 =
860 Teuchos::rcp_dynamic_cast<Stokhos::BlockDiagonalOperator>(dgdp.getLinearOp(),
true);
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 =
868 Teuchos::rcp_dynamic_cast<Stokhos::ProductEpetraMultiVectorOperator>(dgdp_op_mp,
true);
869 Teuchos::RCP<Stokhos::ProductEpetraMultiVector> dgdp_mp =
870 mp_mv_op->productMultiVector();
871 if (me_outargs.supports(OUT_ARG_DgDp_mp,ii,jj).supports(DERIV_MV_BY_COL))
872 me_outargs.set_DgDp_mp(
873 ii, jj, MPDerivative(dgdp_mp, DERIV_MV_BY_COL));
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 =
900 Teuchos::rcp_dynamic_cast<Stokhos::BlockDiagonalOperator>(W,
true);
901 W_mp->setupOperator(W_mp_blocks);
903 if (WPrec_out != Teuchos::null) {
904 Teuchos::RCP<Stokhos::MPPreconditioner> W_prec =
905 Teuchos::rcp_dynamic_cast<Stokhos::MPPreconditioner>(WPrec_out,
true);
906 W_prec->setupPreconditioner(W_mp, *my_x);