42#include "Thyra_EpetraModelEvaluator.hpp"
43#include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
45#include "Thyra_EpetraLinearOp.hpp"
46#include "Thyra_DetachedMultiVectorView.hpp"
47#include "Thyra_ModelEvaluatorDelegatorBase.hpp"
48#include "EpetraExt_ModelEvaluatorScalingTools.h"
49#include "Epetra_RowMatrix.h"
50#include "Teuchos_Time.hpp"
51#include "Teuchos_implicit_cast.hpp"
52#include "Teuchos_Assert.hpp"
53#include "Teuchos_StandardParameterEntryValidators.hpp"
54#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
60const std::string StateFunctionScaling_name =
"State Function Scaling";
63 Thyra::EpetraModelEvaluator::EStateFunctionScaling
66stateFunctionScalingValidator;
67const std::string StateFunctionScaling_default =
"None";
73 const EpetraExt::ModelEvaluator::OutArgs &epetraOutArgs
77 const RCP<Epetra_Operator>
78 eW = epetraOutArgs.get_W();
79 const RCP<Epetra_RowMatrix>
80 ermW = Teuchos::rcp_dynamic_cast<Epetra_RowMatrix>(eW,
false);
82 is_null(ermW), std::logic_error,
83 "Thyra::EpetraModelEvaluator::evalModel(...): Error, if\n"
84 "scaling is turned on, the underlying Epetra_Operator created\n"
85 "an initialized by the underlying epetra model evaluator\n"
86 "\"" << epetraOutArgs.modelEvalDescription() <<
"\"\n"
87 "must support the Epetra_RowMatrix interface through a dynamic cast.\n"
97 const EpetraExt::ModelEvaluator &epetraModel
102 eW = epetraModel.create_W();
105 "Error, the call to create_W() returned null on the "
106 "EpetraExt::ModelEvaluator object "
107 "\"" << epetraModel.description() <<
"\". This may mean that "
108 "the underlying model does not support more than one copy of "
124 :nominalValuesAndBoundsAreUpdated_(false), stateFunctionScaling_(STATE_FUNC_SCALING_NONE),
125 currentInArgsOutArgs_(false), finalPointWasSolved_(false)
133 :nominalValuesAndBoundsAreUpdated_(false), stateFunctionScaling_(STATE_FUNC_SCALING_NONE),
134 currentInArgsOutArgs_(false), finalPointWasSolved_(false)
148 epetraModel_ = epetraModel;
150 W_factory_ = W_factory;
152 x_map_ = epetraModel_->get_x_map();
153 f_map_ = epetraModel_->get_f_map();
154 if (!is_null(x_map_)) {
159 EpetraExt::ModelEvaluator::InArgs inArgs = epetraModel_->createInArgs();
160 p_map_.
resize(inArgs.Np()); p_space_.
resize(inArgs.Np());
161 p_map_is_local_.resize(inArgs.Np(),
false);
162 for(
int l = 0; l < implicit_cast<int>(p_space_.
size()); ++l ) {
164 p_map_l = ( p_map_[l] = epetraModel_->get_p_map(l) );
167 is_null(p_map_l), std::logic_error,
168 "Error, the the map p["<<l<<
"] for the model \""
169 <<epetraModel->description()<<
"\" can not be null!");
172 p_map_is_local_[l] = !p_map_l->DistributedGlobal();
176 EpetraExt::ModelEvaluator::OutArgs outArgs = epetraModel_->createOutArgs();
177 g_map_.
resize(outArgs.Ng()); g_space_.
resize(outArgs.Ng());
178 g_map_is_local_.resize(outArgs.Ng(),
false);
179 for(
int j = 0; j < implicit_cast<int>(g_space_.
size()); ++j ) {
181 g_map_j = ( g_map_[j] = epetraModel_->get_g_map(j) );
182 g_map_is_local_[j] = !g_map_j->DistributedGlobal();
186 epetraInArgsScaling_ = epetraModel_->createInArgs();
187 epetraOutArgsScaling_ = epetraModel_->createOutArgs();
188 nominalValuesAndBoundsAreUpdated_ =
false;
189 finalPointWasSolved_ =
false;
190 stateFunctionScalingVec_ = Teuchos::null;
192 currentInArgsOutArgs_ =
false;
207 nominalValues_.
setArgs(nominalValues);
221 invStateVariableScalingVec_ = Teuchos::null;
222 nominalValuesAndBoundsAreUpdated_ =
false;
229 return stateVariableScalingVec_;
236 updateNominalValuesAndBounds();
237 return invStateVariableScalingVec_;
245 stateFunctionScalingVec_ = stateFunctionScalingVec;
252 return stateFunctionScalingVec_;
261 if(epetraModel) *epetraModel = epetraModel_;
262 if(W_factory) *W_factory = W_factory_;
263 epetraModel_ = Teuchos::null;
264 W_factory_ = Teuchos::null;
265 stateFunctionScalingVec_ = Teuchos::null;
266 stateVariableScalingVec_ = Teuchos::null;
267 invStateVariableScalingVec_ = Teuchos::null;
268 currentInArgsOutArgs_ =
false;
281 return finalPointWasSolved_;
290 std::ostringstream oss;
291 oss <<
"Thyra::EpetraModelEvaluator{";
292 oss <<
"epetraModel=";
293 if(epetraModel_.
get())
294 oss <<
"\'"<<epetraModel_->description()<<
"\'";
297 oss <<
",W_factory=";
299 oss <<
"\'"<<W_factory_->description()<<
"\'";
316 paramList_ = paramList;
317 const EStateFunctionScaling stateFunctionScaling_old = stateFunctionScaling_;
318 stateFunctionScaling_ = stateFunctionScalingValidator->getIntegralValue(
319 *paramList_, StateFunctionScaling_name, StateFunctionScaling_default
321 if( stateFunctionScaling_ != stateFunctionScaling_old )
322 stateFunctionScalingVec_ = Teuchos::null;
323 Teuchos::readVerboseObjectSublist(&*paramList_,
this);
341 paramList_ = Teuchos::null;
358 using Teuchos::tuple;
359 using Teuchos::rcp_implicit_cast;
362 if(is_null(validPL)) {
365 stateFunctionScalingValidator =
rcp(
366 new StringToIntegralParameterEntryValidator<EStateFunctionScaling>(
372 "Do not scale the state function f(...) in this class.",
374 "Scale the state function f(...) and all its derivatives\n"
375 "using the row sum scaling from the initial Jacobian\n"
376 "W=d(f)/d(x). Note, this only works with Epetra_CrsMatrix\n"
379 tuple<EStateFunctionScaling>(
380 STATE_FUNC_SCALING_NONE,
381 STATE_FUNC_SCALING_ROW_SUM
383 StateFunctionScaling_name
386 pl->set(StateFunctionScaling_name,StateFunctionScaling_default,
387 "Determines if and how the state function f(...) and all of its\n"
388 "derivatives are scaled. The scaling is done explicitly so there should\n"
389 "be no impact on the meaning of inner products or tolerances for\n"
391 rcp_implicit_cast<const PEV>(stateFunctionScalingValidator)
393 Teuchos::setupVerboseObjectSublist(&*pl);
405 return p_space_.
size();
411 return g_space_.
size();
445 return epetraModel_->get_p_names(l);
463 return epetraModel_->get_g_names(j);
470 updateNominalValuesAndBounds();
471 return nominalValues_;
478 updateNominalValuesAndBounds();
486 updateNominalValuesAndBounds();
494 return this->create_epetra_W_op();
501 return Teuchos::null;
514 if (!currentInArgsOutArgs_)
515 updateInArgsOutArgs();
516 return prototypeInArgs_;
526 finalPoint_.
setArgs(finalPoint);
527 finalPointWasSolved_ = wasSolved;
535EpetraModelEvaluator::create_DfDp_op_impl(
int )
const
542RCP<LinearOpBase<double> >
543EpetraModelEvaluator::create_DgDx_dot_op_impl(
int )
const
550RCP<LinearOpBase<double> >
551EpetraModelEvaluator::create_DgDx_op_impl(
int )
const
558RCP<LinearOpBase<double> >
559EpetraModelEvaluator::create_DgDp_op_impl(
int ,
int )
const
566ModelEvaluatorBase::OutArgs<double>
567EpetraModelEvaluator::createOutArgsImpl()
const
569 if (!currentInArgsOutArgs_) {
570 updateInArgsOutArgs();
572 return prototypeOutArgs_;
576void EpetraModelEvaluator::evalModelImpl(
577 const ModelEvaluatorBase::InArgs<double>& inArgs_in,
578 const ModelEvaluatorBase::OutArgs<double>& outArgs
583 using Teuchos::rcp_const_cast;
584 using Teuchos::rcp_dynamic_cast;
587 typedef EpetraExt::ModelEvaluator EME;
594 this->updateNominalValuesAndBounds();
605 inArgs.setArgs(inArgs_in);
612 if (
is_null(inArgs_in.get_x_dot()))
613 inArgs.set_x_dot(Teuchos::null);
617 typedef double Scalar;
618 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_GEN_BEGIN(
619 "Thyra::EpetraModelEvaluator",inArgs,outArgs,Teuchos::null
623 const bool firstTimeStateFuncScaling
625 stateFunctionScaling_ != STATE_FUNC_SCALING_NONE
626 &&
is_null(stateFunctionScalingVec_)
630 VOTSLOWSF W_factory_outputTempState(W_factory_,out,verbLevel);
644 *out <<
"\nSetting-up/creating input arguments ...\n";
648 EME::InArgs epetraScaledInArgs = epetraModel_->createInArgs();
649 convertInArgsFromThyraToEpetra( inArgs, &epetraScaledInArgs );
653 EME::InArgs epetraInArgs = epetraModel_->createInArgs();
654 EpetraExt::unscaleModelVars(
655 epetraScaledInArgs, epetraInArgsScaling_, &epetraInArgs,
661 OSTab(out).o() <<
"\nTime to setup InArgs = "<<timer.totalElapsedTime()<<
" sec\n";
668 *out <<
"\nSetting-up/creating output arguments ...\n";
673 EME::OutArgs epetraUnscaledOutArgs = epetraModel_->createOutArgs();
677 RCP<LinearOpBase<double> > W_op;
678 RCP<EpetraLinearOp> efwdW;
679 RCP<Epetra_Operator> eW;
683 convertOutArgsFromThyraToEpetra(
685 &epetraUnscaledOutArgs,
693 if (firstTimeStateFuncScaling) {
694 preEvalScalingSetup(&epetraInArgs,&epetraUnscaledOutArgs,out,verbLevel);
700 <<
"\nTime to setup OutArgs = "
701 << timer.totalElapsedTime() <<
" sec\n";
708 *out <<
"\nEvaluating the Epetra output functions ...\n";
711 epetraModel_->evalModel(epetraInArgs,epetraUnscaledOutArgs);
716 <<
"\nTime to evaluate Epetra output functions = "
717 << timer.totalElapsedTime() <<
" sec\n";
728 *out <<
"\nCompute scale factors if needed ...\n";
731 if (firstTimeStateFuncScaling) {
732 postEvalScalingSetup(epetraUnscaledOutArgs,out,verbLevel);
738 <<
"\nTime to compute scale factors = "
739 << timer.totalElapsedTime() <<
" sec\n";
746 *out <<
"\nScale the output objects ...\n";
749 EME::OutArgs epetraOutArgs = epetraModel_->createOutArgs();
750 bool allFuncsWhereScaled =
false;
751 EpetraExt::scaleModelFuncs(
752 epetraUnscaledOutArgs, epetraInArgsScaling_, epetraOutArgsScaling_,
753 &epetraOutArgs, &allFuncsWhereScaled,
757 if (stateFunctionScaling_ != STATE_FUNC_SCALING_NONE)
759 !allFuncsWhereScaled, std::logic_error,
760 "Error, we can not currently handle epetra output objects that could not be"
761 " scaled. Special code will have to be added to handle this (i.e. using"
762 " implicit diagonal and multiplied linear operators to implicitly do"
769 <<
"\nTime to scale the output objects = "
770 << timer.totalElapsedTime() <<
" sec\n";
778 *out <<
"\nFinish processing and wrapping the output objects ...\n";
781 finishConvertingOutArgsFromEpetraToThyra(
782 epetraOutArgs, W_op, efwdW, eW,
789 <<
"\nTime to finish processing and wrapping the output objects = "
790 << timer.totalElapsedTime() <<
" sec\n";
796 THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
804void EpetraModelEvaluator::convertInArgsFromEpetraToThyra(
805 const EpetraExt::ModelEvaluator::InArgs &epetraInArgs,
806 ModelEvaluatorBase::InArgs<double> *inArgs
815 if(inArgs->supports(MEB::IN_ARG_x)) {
816 inArgs->set_x(
create_Vector( epetraInArgs.get_x(), x_space_ ) );
819 if(inArgs->supports(MEB::IN_ARG_x_dot)) {
820 inArgs->set_x_dot(
create_Vector( epetraInArgs.get_x_dot(), x_space_ ) );
823 if(inArgs->supports(MEB::IN_ARG_x_mp)) {
824 inArgs->set_x_mp( epetraInArgs.get_x_mp() );
827 if(inArgs->supports(MEB::IN_ARG_x_dot_mp)) {
828 inArgs->set_x_dot_mp( epetraInArgs.get_x_dot_mp() );
831 const int l_Np = inArgs->Np();
832 for(
int l = 0; l < l_Np; ++l ) {
833 inArgs->set_p( l,
create_Vector( epetraInArgs.get_p(l), p_space_[l] ) );
835 for(
int l = 0; l < l_Np; ++l ) {
836 if(inArgs->supports(MEB::IN_ARG_p_mp, l))
837 inArgs->set_p_mp( l, epetraInArgs.get_p_mp(l) );
840 if(inArgs->supports(MEB::IN_ARG_t)) {
841 inArgs->set_t(epetraInArgs.get_t());
847void EpetraModelEvaluator::convertInArgsFromThyraToEpetra(
848 const ModelEvaluatorBase::InArgs<double> &inArgs,
849 EpetraExt::ModelEvaluator::InArgs *epetraInArgs
854 using Teuchos::rcp_const_cast;
855#ifdef HAVE_THYRA_ME_POLYNOMIAL
862 RCP<const VectorBase<double> > x_dot;
863 if( inArgs.supports(
IN_ARG_x_dot) && (x_dot = inArgs.get_x_dot()).get() ) {
865 epetraInArgs->set_x_dot(e_x_dot);
867 RCP<const Stokhos::ProductEpetraVector > x_dot_mp;
868 if( inArgs.supports(
IN_ARG_x_dot_mp) && (x_dot_mp = inArgs.get_x_dot_mp()).get() ) {
869 epetraInArgs->set_x_dot_mp(x_dot_mp);
872 RCP<const VectorBase<double> > x;
873 if( inArgs.supports(
IN_ARG_x) && (x = inArgs.get_x()).get() ) {
875 epetraInArgs->set_x(e_x);
877 RCP<const Stokhos::ProductEpetraVector > x_mp;
878 if( inArgs.supports(
IN_ARG_x_mp) && (x_mp = inArgs.get_x_mp()).get() ) {
879 epetraInArgs->set_x_mp(x_mp);
882 RCP<const VectorBase<double> > p_l;
883 for(
int l = 0; l < inArgs.Np(); ++l ) {
884 p_l = inArgs.get_p(l);
887 RCP<const Stokhos::ProductEpetraVector > p_mp_l;
888 for(
int l = 0; l < inArgs.Np(); ++l ) {
890 p_mp_l = inArgs.get_p_mp(l);
891 if(p_mp_l.get()) epetraInArgs->set_p_mp(l,p_mp_l);
895#ifdef HAVE_THYRA_ME_POLYNOMIAL
897 RCP<const Polynomial< VectorBase<double> > > x_dot_poly;
898 RCP<Epetra_Vector> epetra_ptr;
901 && (x_dot_poly = inArgs.get_x_dot_poly()).get()
904 RCP<Polynomial<Epetra_Vector> > epetra_x_dot_poly =
905 rcp(
new Polynomial<Epetra_Vector>(x_dot_poly->degree()));
906 for (
unsigned int i=0; i<=x_dot_poly->degree(); i++) {
907 epetra_ptr = rcp_const_cast<Epetra_Vector>(
909 epetra_x_dot_poly->setCoefficientPtr(i,epetra_ptr);
911 epetraInArgs->set_x_dot_poly(epetra_x_dot_poly);
914 RCP<const Polynomial< VectorBase<double> > > x_poly;
917 && (x_poly = inArgs.get_x_poly()).get()
920 RCP<Polynomial<Epetra_Vector> > epetra_x_poly =
921 rcp(
new Polynomial<Epetra_Vector>(x_poly->degree()));
922 for (
unsigned int i=0; i<=x_poly->degree(); i++) {
923 epetra_ptr = rcp_const_cast<Epetra_Vector>(
925 epetra_x_poly->setCoefficientPtr(i,epetra_ptr);
927 epetraInArgs->set_x_poly(epetra_x_poly);
933 epetraInArgs->set_t(inArgs.get_t());
936 epetraInArgs->set_alpha(inArgs.get_alpha());
939 epetraInArgs->set_beta(inArgs.get_beta());
942 epetraInArgs->set_step_size(inArgs.get_step_size());
945 epetraInArgs->set_stage_number(inArgs.get_stage_number());
949void EpetraModelEvaluator::convertOutArgsFromThyraToEpetra(
950 const ModelEvaluatorBase::OutArgs<double> &outArgs,
951 EpetraExt::ModelEvaluator::OutArgs *epetraUnscaledOutArgs_inout,
952 RCP<LinearOpBase<double> > *W_op_out,
953 RCP<EpetraLinearOp> *efwdW_out,
954 RCP<Epetra_Operator> *eW_out
959 using Teuchos::rcp_const_cast;
960 using Teuchos::rcp_dynamic_cast;
975 EpetraExt::ModelEvaluator::OutArgs &epetraUnscaledOutArgs = *epetraUnscaledOutArgs_inout;
976 RCP<LinearOpBase<double> > &W_op = *W_op_out;
977 RCP<EpetraLinearOp> &efwdW = *efwdW_out;
978 RCP<Epetra_Operator> &eW = *eW_out;
982 RCP<VectorBase<double> > f;
983 if( outArgs.supports(
OUT_ARG_f) && (f = outArgs.get_f()).get() )
985 RCP<Stokhos::ProductEpetraVector > f_mp;
986 if( outArgs.supports(
OUT_ARG_f_mp) && (f_mp = outArgs.get_f_mp()).get() )
987 epetraUnscaledOutArgs.set_f_mp(f_mp);
992 RCP<VectorBase<double> > g_j;
993 for(
int j = 0; j < outArgs.Ng(); ++j ) {
994 g_j = outArgs.get_g(j);
997 RCP<Stokhos::ProductEpetraVector > g_mp_j;
998 for(
int j = 0; j < outArgs.Ng(); ++j ) {
1000 g_mp_j = outArgs.get_g_mp(j);
1001 if(g_mp_j.get()) epetraUnscaledOutArgs.set_g_mp(j,g_mp_j);
1008 RCP<Stokhos::ProductEpetraOperator > W_mp;
1009 if( outArgs.supports(
OUT_ARG_W_mp) && (W_mp = outArgs.get_W_mp()).get() )
1010 epetraUnscaledOutArgs.set_W_mp(W_mp);
1018 efwdW = rcp_const_cast<EpetraLinearOp>(
1019 rcp_dynamic_cast<const EpetraLinearOp>(W_op,
true));
1028 eW = efwdW->epetra_op();
1029 epetraUnscaledOutArgs.set_W(eW);
1038 Derivative<double> DfDp_l;
1039 for(
int l = 0; l < outArgs.Np(); ++l ) {
1041 && !(DfDp_l = outArgs.get_DfDp(l)).isEmpty() )
1043 epetraUnscaledOutArgs.set_DfDp(l,
convert(DfDp_l,f_map_,p_map_[l]));
1046 MPDerivative DfDp_mp_l;
1047 for(
int l = 0; l < outArgs.Np(); ++l ) {
1049 && !(DfDp_mp_l = outArgs.get_DfDp_mp(l)).isEmpty() )
1051 epetraUnscaledOutArgs.set_DfDp_mp(l,
convert(DfDp_mp_l,f_map_,p_map_[l]));
1058 Derivative<double> DgDx_dot_j;
1059 for(
int j = 0; j < outArgs.Ng(); ++j ) {
1061 && !(DgDx_dot_j = outArgs.get_DgDx_dot(j)).isEmpty() )
1063 epetraUnscaledOutArgs.set_DgDx_dot(j,
convert(DgDx_dot_j,g_map_[j],x_map_));
1066 MPDerivative DgDx_dot_mp_j;
1067 for(
int j = 0; j < outArgs.Ng(); ++j ) {
1069 && !(DgDx_dot_mp_j = outArgs.get_DgDx_dot_mp(j)).isEmpty() )
1071 epetraUnscaledOutArgs.set_DgDx_dot_mp(j,
convert(DgDx_dot_mp_j,g_map_[j],x_map_));
1078 Derivative<double> DgDx_j;
1079 for(
int j = 0; j < outArgs.Ng(); ++j ) {
1081 && !(DgDx_j = outArgs.get_DgDx(j)).isEmpty() )
1083 epetraUnscaledOutArgs.set_DgDx(j,
convert(DgDx_j,g_map_[j],x_map_));
1086 MPDerivative DgDx_mp_j;
1087 for(
int j = 0; j < outArgs.Ng(); ++j ) {
1089 && !(DgDx_mp_j = outArgs.get_DgDx_mp(j)).isEmpty() )
1091 epetraUnscaledOutArgs.set_DgDx_mp(j,
convert(DgDx_mp_j,g_map_[j],x_map_));
1098 DerivativeSupport DgDp_j_l_support;
1099 Derivative<double> DgDp_j_l;
1100 for (
int j = 0; j < outArgs.Ng(); ++j ) {
1101 for (
int l = 0; l < outArgs.Np(); ++l ) {
1102 if (!(DgDp_j_l_support = outArgs.supports(
OUT_ARG_DgDp,j,l)).none()
1103 && !(DgDp_j_l = outArgs.get_DgDp(j,l)).isEmpty() )
1105 epetraUnscaledOutArgs.set_DgDp(j,l,
convert(DgDp_j_l,g_map_[j],p_map_[l]));
1109 DerivativeSupport DgDp_mp_j_l_support;
1110 MPDerivative DgDp_mp_j_l;
1111 for (
int j = 0; j < outArgs.Ng(); ++j ) {
1112 for (
int l = 0; l < outArgs.Np(); ++l ) {
1113 if (!(DgDp_mp_j_l_support = outArgs.supports(
OUT_ARG_DgDp_mp,j,l)).none()
1114 && !(DgDp_mp_j_l = outArgs.get_DgDp_mp(j,l)).isEmpty() )
1116 epetraUnscaledOutArgs.set_DgDp_mp(j,l,
convert(DgDp_mp_j_l,g_map_[j],p_map_[l]));
1122#ifdef HAVE_THYRA_ME_POLYNOMIAL
1125 RCP<const Teuchos::Polynomial< VectorBase<double> > > f_poly;
1126 if( outArgs.supports(
OUT_ARG_f_poly) && (f_poly = outArgs.get_f_poly()).get() )
1128 RCP<Teuchos::Polynomial<Epetra_Vector> > epetra_f_poly =
1130 for (
unsigned int i=0; i<=f_poly->degree(); i++) {
1131 RCP<Epetra_Vector> epetra_ptr
1133 f_poly->getCoefficient(i)));
1134 epetra_f_poly->setCoefficientPtr(i,epetra_ptr);
1136 epetraUnscaledOutArgs.set_f_poly(epetra_f_poly);
1144void EpetraModelEvaluator::preEvalScalingSetup(
1145 EpetraExt::ModelEvaluator::InArgs *epetraInArgs_inout,
1146 EpetraExt::ModelEvaluator::OutArgs *epetraUnscaledOutArgs_inout,
1147 const RCP<Teuchos::FancyOStream> &out,
1152 typedef EpetraExt::ModelEvaluator EME;
1159 EpetraExt::ModelEvaluator::InArgs
1160 &epetraInArgs = *epetraInArgs_inout;
1161 EpetraExt::ModelEvaluator::OutArgs
1162 &epetraUnscaledOutArgs = *epetraUnscaledOutArgs_inout;
1165 ( stateFunctionScaling_ == STATE_FUNC_SCALING_ROW_SUM )
1168 epetraUnscaledOutArgs.supports(EME::OUT_ARG_f)
1170 epetraUnscaledOutArgs.funcOrDerivesAreSet(EME::OUT_ARG_f)
1174 epetraUnscaledOutArgs.supports(EME::OUT_ARG_W)
1176 is_null(epetraUnscaledOutArgs.get_W())
1191 <<
"\nCreating a temporary Epetra W to compute scale factors"
1192 <<
" for f(...) ...\n";
1193 epetraUnscaledOutArgs.set_W(create_and_assert_W(*epetraModel_));
1194 if( epetraInArgs.supports(EME::IN_ARG_beta) )
1195 epetraInArgs.set_beta(1.0);
1196 if( epetraInArgs.supports(EME::IN_ARG_alpha) )
1197 epetraInArgs.set_alpha(0.0);
1198 if( epetraInArgs.supports(EME::IN_ARG_step_size) )
1199 epetraInArgs.set_step_size(0.0);
1200 if( epetraInArgs.supports(EME::IN_ARG_stage_number) )
1201 epetraInArgs.set_stage_number(1.0);
1207void EpetraModelEvaluator::postEvalScalingSetup(
1208 const EpetraExt::ModelEvaluator::OutArgs &epetraUnscaledOutArgs,
1209 const RCP<Teuchos::FancyOStream> &out,
1216 using Teuchos::rcp_const_cast;
1220 switch(stateFunctionScaling_) {
1222 case STATE_FUNC_SCALING_ROW_SUM: {
1226 const RCP<Epetra_RowMatrix>
1227 ermW = get_Epetra_RowMatrix(epetraUnscaledOutArgs);
1238 ermW->InvRowSums(*invRowSums);
1242 <<
"\nComputed inverse row sum scaling from W that"
1243 " will be used to scale f(...) and its derivatives:\n";
1244 double minVal = 0, maxVal = 0, avgVal = 0;
1245 invRowSums->MinValue(&minVal);
1246 invRowSums->MaxValue(&maxVal);
1247 invRowSums->MeanValue(&avgVal);
1250 <<
"min(invRowSums) = " << minVal <<
"\n"
1251 <<
"max(invRowSums) = " << maxVal <<
"\n"
1252 <<
"avg(invRowSums) = " << avgVal <<
"\n";
1255 stateFunctionScalingVec_ = invRowSums;
1266 epetraOutArgsScaling_ = epetraModel_->createOutArgs();
1268 epetraOutArgsScaling_.set_f(
1269 rcp_const_cast<Epetra_Vector>(stateFunctionScalingVec_) );
1274void EpetraModelEvaluator::finishConvertingOutArgsFromEpetraToThyra(
1275 const EpetraExt::ModelEvaluator::OutArgs &,
1276 RCP<LinearOpBase<double> > &W_op,
1277 RCP<EpetraLinearOp> &efwdW,
1278 RCP<Epetra_Operator> &,
1279 const ModelEvaluatorBase::OutArgs<double> &
1283 using Teuchos::rcp_dynamic_cast;
1287 efwdW->setFullyInitialized(
true);
1292 if (W_op.shares_resource(efwdW)) {
1296 rcp_dynamic_cast<EpetraLinearOp>(W_op,
true)->setFullyInitialized(
true);
1303void EpetraModelEvaluator::updateNominalValuesAndBounds()
const
1309 typedef EpetraExt::ModelEvaluator EME;
1311 if( !nominalValuesAndBoundsAreUpdated_ ) {
1315 EME::InArgs epetraOrigNominalValues;
1316 EpetraExt::gatherModelNominalValues(
1317 *epetraModel_, &epetraOrigNominalValues );
1319 EME::InArgs epetraOrigLowerBounds;
1320 EME::InArgs epetraOrigUpperBounds;
1321 EpetraExt::gatherModelBounds(
1322 *epetraModel_, &epetraOrigLowerBounds, &epetraOrigUpperBounds );
1326 epetraInArgsScaling_ = epetraModel_->createInArgs();
1328 if( !
is_null(stateVariableScalingVec_) ) {
1329 invStateVariableScalingVec_
1330 = EpetraExt::createInverseModelScalingVector(stateVariableScalingVec_);
1331 if( epetraOrigNominalValues.supports(EME::IN_ARG_x_dot) ) {
1332 epetraInArgsScaling_.set_x_dot(invStateVariableScalingVec_);
1334 if( epetraOrigNominalValues.supports(EME::IN_ARG_x) ) {
1335 epetraInArgsScaling_.set_x(invStateVariableScalingVec_);
1341 EME::InArgs epetraScaledNominalValues = epetraModel_->createInArgs();
1342 EpetraExt::scaleModelVars(
1343 epetraOrigNominalValues, epetraInArgsScaling_, &epetraScaledNominalValues
1346 EME::InArgs epetraScaledLowerBounds = epetraModel_->createInArgs();
1347 EME::InArgs epetraScaledUpperBounds = epetraModel_->createInArgs();
1348 EpetraExt::scaleModelBounds(
1349 epetraOrigLowerBounds, epetraOrigUpperBounds, epetraModel_->getInfBound(),
1350 epetraInArgsScaling_,
1351 &epetraScaledLowerBounds, &epetraScaledUpperBounds
1359 convertInArgsFromEpetraToThyra(epetraScaledNominalValues, &nominalValues_);
1360 convertInArgsFromEpetraToThyra(epetraScaledLowerBounds, &lowerBounds_);
1361 convertInArgsFromEpetraToThyra(epetraScaledUpperBounds, &upperBounds_);
1363 nominalValuesAndBoundsAreUpdated_ =
true;
1376void EpetraModelEvaluator::updateInArgsOutArgs()
const
1379 typedef EpetraExt::ModelEvaluator EME;
1381 const EpetraExt::ModelEvaluator &epetraModel = *epetraModel_;
1382 EME::InArgs epetraInArgs = epetraModel.createInArgs();
1383 EME::OutArgs epetraOutArgs = epetraModel.createOutArgs();
1384 const int l_Np = epetraOutArgs.Np();
1385 const int l_Ng = epetraOutArgs.Ng();
1391 InArgsSetup<double> inArgs;
1392 inArgs.setModelEvalDescription(this->
description());
1393 inArgs.set_Np(epetraInArgs.Np());
1394 inArgs.setSupports(
IN_ARG_x_dot, epetraInArgs.supports(EME::IN_ARG_x_dot));
1395 inArgs.setSupports(
IN_ARG_x, epetraInArgs.supports(EME::IN_ARG_x));
1396 inArgs.setSupports(
IN_ARG_x_dot_mp, epetraInArgs.supports(EME::IN_ARG_x_dot_mp));
1397 inArgs.setSupports(
IN_ARG_x_mp, epetraInArgs.supports(EME::IN_ARG_x_mp));
1398#ifdef HAVE_THYRA_ME_POLYNOMIAL
1400 epetraInArgs.supports(EME::IN_ARG_x_dot_poly));
1401 inArgs.setSupports(
IN_ARG_x_poly, epetraInArgs.supports(EME::IN_ARG_x_poly));
1403 inArgs.setSupports(
IN_ARG_t, epetraInArgs.supports(EME::IN_ARG_t));
1404 inArgs.setSupports(
IN_ARG_alpha, epetraInArgs.supports(EME::IN_ARG_alpha));
1405 inArgs.setSupports(
IN_ARG_beta, epetraInArgs.supports(EME::IN_ARG_beta));
1406 inArgs.setSupports(
IN_ARG_step_size, epetraInArgs.supports(EME::IN_ARG_step_size));
1408 for(
int l=0; l<l_Np; ++l) {
1409 inArgs.setSupports(
IN_ARG_p_mp, l, epetraInArgs.supports(EME::IN_ARG_p_mp, l));
1411 prototypeInArgs_ = inArgs;
1417 OutArgsSetup<double> outArgs;
1418 outArgs.setModelEvalDescription(this->
description());
1419 outArgs.set_Np_Ng(l_Np, l_Ng);
1421 outArgs.setSupports(
OUT_ARG_f, epetraOutArgs.supports(EME::OUT_ARG_f));
1422 outArgs.setSupports(
OUT_ARG_f_mp, epetraOutArgs.supports(EME::OUT_ARG_f_mp));
1425 outArgs.setSupports(
OUT_ARG_W_op, epetraOutArgs.supports(EME::OUT_ARG_W));
1426 outArgs.set_W_properties(
convert(epetraOutArgs.get_W_properties()));
1428 for(
int l=0; l<l_Np; ++l) {
1430 convert(epetraOutArgs.supports(EME::OUT_ARG_DfDp, l)));
1432 outArgs.set_DfDp_properties(l,
1433 convert(epetraOutArgs.get_DfDp_properties(l)));
1437 convert(epetraOutArgs.supports(EME::OUT_ARG_DfDp_mp, l)));
1439 outArgs.set_DfDp_mp_properties(l,
1440 convert(epetraOutArgs.get_DfDp_mp_properties(l)));
1445 for(
int j=0; j<l_Ng; ++j) {
1448 convert(epetraOutArgs.supports(EME::OUT_ARG_DgDx_dot, j)));
1450 outArgs.set_DgDx_dot_properties(j,
1451 convert(epetraOutArgs.get_DgDx_dot_properties(j)));
1454 convert(epetraOutArgs.supports(EME::OUT_ARG_DgDx, j)));
1456 outArgs.set_DgDx_properties(j,
1457 convert(epetraOutArgs.get_DgDx_properties(j)));
1460 convert(epetraOutArgs.supports(EME::OUT_ARG_DgDx_dot_mp, j)));
1462 outArgs.set_DgDx_dot_mp_properties(j,
1463 convert(epetraOutArgs.get_DgDx_dot_mp_properties(j)));
1466 convert(epetraOutArgs.supports(EME::OUT_ARG_DgDx_mp, j)));
1468 outArgs.set_DgDx_mp_properties(j,
1469 convert(epetraOutArgs.get_DgDx_mp_properties(j)));
1472 for(
int j=0; j < l_Ng; ++j)
for(
int l=0; l < l_Np; ++l) {
1473 const EME::DerivativeSupport epetra_DgDp_j_l_support =
1474 epetraOutArgs.
supports(EME::OUT_ARG_DgDp, j, l);
1476 convert(epetra_DgDp_j_l_support));
1478 outArgs.set_DgDp_properties(j, l,
1479 convert(epetraOutArgs.get_DgDp_properties(j, l)));
1480 const EME::DerivativeSupport epetra_DgDp_mp_j_l_support =
1481 epetraOutArgs.supports(EME::OUT_ARG_DgDp_mp, j, l);
1483 convert(epetra_DgDp_mp_j_l_support));
1485 outArgs.set_DgDp_mp_properties(j, l,
1486 convert(epetraOutArgs.get_DgDp_mp_properties(j, l)));
1488 for(
int j=0; j<l_Ng; ++j) {
1489 outArgs.setSupports(
OUT_ARG_g_mp, j, epetraOutArgs.supports(EME::OUT_ARG_g_mp, j));
1491#ifdef HAVE_THYRA_ME_POLYNOMIAL
1493 epetraOutArgs.supports(EME::OUT_ARG_f_poly));
1495 prototypeOutArgs_ = outArgs;
1498 currentInArgsOutArgs_ =
true;
1504EpetraModelEvaluator::create_epetra_W_op()
const
1506 return Thyra::partialNonconstEpetraLinearOp(
1508 create_and_assert_W(*epetraModel_)
1522Thyra::epetraModelEvaluator(
1523 const RCP<const EpetraExt::ModelEvaluator> &epetraModel,
1524 const RCP<LinearOpWithSolveFactoryBase<double> > &W_factory
1527 return Teuchos::rcp(
new EpetraModelEvaluator(epetraModel,W_factory));
1533 const EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation &mvOrientation
1536 switch(mvOrientation) {
1537 case EpetraExt::ModelEvaluator::DERIV_MV_BY_COL :
1538 return ModelEvaluatorBase::DERIV_MV_BY_COL;
1539 case EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW :
1540 return ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW;
1548EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation
1550 const ModelEvaluatorBase::EDerivativeMultiVectorOrientation &mvOrientation
1553 switch(mvOrientation) {
1554 case ModelEvaluatorBase::DERIV_MV_BY_COL :
1555 return EpetraExt::ModelEvaluator::DERIV_MV_BY_COL;
1556 case ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW :
1557 return EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW;
1567 const EpetraExt::ModelEvaluator::DerivativeProperties &derivativeProperties
1570 ModelEvaluatorBase::EDerivativeLinearity linearity;
1571 switch(derivativeProperties.linearity) {
1572 case EpetraExt::ModelEvaluator::DERIV_LINEARITY_UNKNOWN:
1573 linearity = ModelEvaluatorBase::DERIV_LINEARITY_UNKNOWN;
1575 case EpetraExt::ModelEvaluator::DERIV_LINEARITY_CONST:
1576 linearity = ModelEvaluatorBase::DERIV_LINEARITY_CONST;
1578 case EpetraExt::ModelEvaluator::DERIV_LINEARITY_NONCONST:
1579 linearity = ModelEvaluatorBase::DERIV_LINEARITY_NONCONST;
1584 ModelEvaluatorBase::ERankStatus rank;
1585 switch(derivativeProperties.rank) {
1586 case EpetraExt::ModelEvaluator::DERIV_RANK_UNKNOWN:
1587 rank = ModelEvaluatorBase::DERIV_RANK_UNKNOWN;
1589 case EpetraExt::ModelEvaluator::DERIV_RANK_FULL:
1590 rank = ModelEvaluatorBase::DERIV_RANK_FULL;
1592 case EpetraExt::ModelEvaluator::DERIV_RANK_DEFICIENT:
1593 rank = ModelEvaluatorBase::DERIV_RANK_DEFICIENT;
1598 return ModelEvaluatorBase::DerivativeProperties(
1599 linearity,rank,derivativeProperties.supportsAdjoint);
1605 const EpetraExt::ModelEvaluator::DerivativeSupport &derivativeSupport
1608 ModelEvaluatorBase::DerivativeSupport ds;
1609 if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_LINEAR_OP))
1610 ds.plus(ModelEvaluatorBase::DERIV_LINEAR_OP);
1611 if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_MV_BY_COL))
1612 ds.plus(ModelEvaluatorBase::DERIV_MV_BY_COL);
1613 if(derivativeSupport.supports(EpetraExt::ModelEvaluator::DERIV_TRANS_MV_BY_ROW))
1614 ds.plus(ModelEvaluatorBase::DERIV_TRANS_MV_BY_ROW);
1619EpetraExt::ModelEvaluator::Derivative
1621 const ModelEvaluatorBase::Derivative<double> &derivative,
1622 const RCP<const Epetra_Map> &fnc_map,
1623 const RCP<const Epetra_Map> &var_map
1626 typedef ModelEvaluatorBase MEB;
1627 if(derivative.getLinearOp().get()) {
1628 return EpetraExt::ModelEvaluator::Derivative(
1629 Teuchos::rcp_const_cast<Epetra_Operator>(
1630 Teuchos::dyn_cast<const EpetraLinearOp>(*derivative.getLinearOp()).epetra_op()
1634 else if(derivative.getDerivativeMultiVector().getMultiVector().get()) {
1635 return EpetraExt::ModelEvaluator::Derivative(
1636 EpetraExt::ModelEvaluator::DerivativeMultiVector(
1638 ( derivative.getDerivativeMultiVector().getOrientation() == MEB::DERIV_MV_BY_COL
1642 ,derivative.getDerivativeMultiVector().getMultiVector()
1644 ,convert(derivative.getDerivativeMultiVector().getOrientation())
1648 return EpetraExt::ModelEvaluator::Derivative();
1650EpetraExt::ModelEvaluator::MPDerivative
1652 const ModelEvaluatorBase::MPDerivative &derivative,
1653 const RCP<const Epetra_Map> &,
1654 const RCP<const Epetra_Map> &
1658 if(derivative.getLinearOp().get()) {
1659 return EpetraExt::ModelEvaluator::MPDerivative(
1660 derivative.getLinearOp()
1663 else if(derivative.getDerivativeMultiVector().getMultiVector().get()) {
1664 return EpetraExt::ModelEvaluator::MPDerivative(
1665 EpetraExt::ModelEvaluator::MPDerivativeMultiVector(
1666 derivative.getDerivativeMultiVector().getMultiVector()
1667 ,convert(derivative.getDerivativeMultiVector().getOrientation())
1671 return EpetraExt::ModelEvaluator::MPDerivative();
void resize(size_type new_size, const value_type &x=value_type())
const RCP< T > & assert_not_null() const
RCP< Teuchos::ParameterList > unsetParameterList()
RCP< Teuchos::ParameterList > getNonconstParameterList()
bool finalPointWasSolved() const
void setStateFunctionScalingVec(const RCP< const Epetra_Vector > &stateFunctionScalingVec)
Set the state function scaling vector s_f (see above).
RCP< const Epetra_Vector > getStateFunctionScalingVec() const
Get the state function scaling vector s_f (see above).
std::string description() const
void initialize(const RCP< const EpetraExt::ModelEvaluator > &epetraModel, const RCP< LinearOpWithSolveFactoryBase< double > > &W_factory)
RCP< const VectorSpaceBase< double > > get_f_space() const
void setStateVariableScalingVec(const RCP< const Epetra_Vector > &stateVariableScalingVec)
Set the state variable scaling vector s_x (see above).
RCP< LinearOpBase< double > > create_W_op() const
ModelEvaluatorBase::InArgs< double > getUpperBounds() const
void setNominalValues(const ModelEvaluatorBase::InArgs< double > &nominalValues)
Set the nominal values.
RCP< PreconditionerBase< double > > create_W_prec() const
Returns null currently.
RCP< const LinearOpWithSolveFactoryBase< double > > get_W_factory() const
ModelEvaluatorBase::InArgs< double > getNominalValues() const
RCP< const Teuchos::ParameterList > getParameterList() const
void reportFinalPoint(const ModelEvaluatorBase::InArgs< double > &finalPoint, const bool wasSolved)
ModelEvaluatorBase::InArgs< double > createInArgs() const
void setParameterList(RCP< Teuchos::ParameterList > const ¶mList)
RCP< const VectorSpaceBase< double > > get_g_space(int j) const
RCP< const EpetraExt::ModelEvaluator > getEpetraModel() const
void uninitialize(RCP< const EpetraExt::ModelEvaluator > *epetraModel=NULL, RCP< LinearOpWithSolveFactoryBase< double > > *W_factory=NULL)
RCP< const VectorSpaceBase< double > > get_x_space() const
RCP< const Epetra_Vector > getStateVariableInvScalingVec() const
Get the state variable scaling vector s_x (see above).
ModelEvaluatorBase::InArgs< double > getLowerBounds() const
RCP< const Epetra_Vector > getStateVariableScalingVec() const
Get the inverse state variable scaling vector inv_s_x (see above).
RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::ArrayView< const std::string > get_g_names(int j) const
RCP< const Teuchos::Array< std::string > > get_p_names(int l) const
RCP< const VectorSpaceBase< double > > get_p_space(int l) const
const ModelEvaluatorBase::InArgs< double > & getFinalPoint() const
ModelEvaluatorBase::EDerivativeMultiVectorOrientation convert(const EpetraExt::ModelEvaluator::EDerivativeMultiVectorOrientation &mvOrientation)
Factory interface for creating LinearOpWithSolveBase objects from compatible LinearOpBase objects.
Determines the forms of a general derivative that are supported.
Concrete aggregate class for all input arguments computable by a ModelEvaluator subclass object.
void setArgs(const InArgs< Scalar > &inArgs, bool ignoreUnsupported=false, bool cloneObjects=false)
Set non-null arguments (does not overwrite non-NULLs with NULLs) .
bool supports(EInArgsMembers arg) const
Determines if an input argument is supported or not.
Base subclass for ModelEvaluator that defines some basic types.
ModelEvaluatorBase()
constructor
EDerivativeMultiVectorOrientation
RCP< VectorBase< double > > create_Vector(const RCP< Epetra_Vector > &epetra_v, const RCP< const VectorSpaceBase< double > > &space=Teuchos::null)
Create a non-const VectorBase object from a non-const Epetra_Vector object.
RCP< const VectorSpaceBase< double > > create_VectorSpace(const RCP< const Epetra_Map > &epetra_map)
Create an VectorSpaceBase object given an Epetra_Map object.
RCP< Epetra_Vector > get_Epetra_Vector(const Epetra_Map &map, const RCP< VectorBase< double > > &v)
Get a non-const Epetra_Vector view from a non-const VectorBase object if possible.
RCP< Epetra_MultiVector > get_Epetra_MultiVector(const Epetra_Map &map, const RCP< MultiVectorBase< double > > &mv)
Get a non-const Epetra_MultiVector view from a non-const MultiVectorBase object if possible.
#define TEUCHOS_ASSERT(assertion_test)
#define TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(index, lower_inclusive, upper_exclusive)
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
bool is_null(const std::shared_ptr< T > &p)
bool nonnull(const std::shared_ptr< T > &p)
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
TypeTo implicit_cast(const TypeFrom &t)
std::string typeName(const T &t)
basic_OSTab< char > OSTab
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
TEUCHOSCORE_LIB_DLL_EXPORT bool includesVerbLevel(const EVerbosityLevel verbLevel, const EVerbosityLevel requestedVerbLevel, const bool isDefaultLevel=false)
Simple public strict containing properties of a derivative object.