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 =
720 Teuchos::rcp_dynamic_cast<Stokhos::SGPreconditioner>(WPrec_out,
true);
721 W_prec->setupPreconditioner(my_W, *my_x);
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 =
826 Teuchos::rcp_dynamic_cast<Stokhos::EpetraOperatorOrthogPoly>(dfdp.getLinearOp(),
true);
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 =
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));
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));
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 =
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);
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,
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 =
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);
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,
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 =
965 Teuchos::rcp_dynamic_cast<Stokhos::EpetraOperatorOrthogPoly>(dgdp.getLinearOp(),
true);
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 =
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));
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));
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 =
1018 Teuchos::rcp_dynamic_cast<Stokhos::SGOperator>(W,
true);
1019 W_sg->setupOperator(W_sg_blocks);
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);
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);