42 #ifndef THYRA_DIRECTIONAL_FINITE_DIFF_CALCULATOR_DEF_HPP 43 #define THYRA_DIRECTIONAL_FINITE_DIFF_CALCULATOR_DEF_HPP 46 #include "Thyra_DirectionalFiniteDiffCalculator_decl.hpp" 47 #include "Thyra_ModelEvaluatorHelpers.hpp" 48 #include "Thyra_DetachedVectorView.hpp" 49 #include "Thyra_DetachedMultiVectorView.hpp" 50 #include "Thyra_StateFuncModelEvaluatorBase.hpp" 51 #include "Thyra_MultiVectorStdOps.hpp" 52 #include "Thyra_VectorStdOps.hpp" 53 #include "Teuchos_TimeMonitor.hpp" 54 #include "Teuchos_VerboseObjectParameterListHelpers.hpp" 60 namespace DirectionalFiniteDiffCalculatorTypes {
73 template<
class Scalar>
74 class OutArgsCreator :
public StateFuncModelEvaluatorBase<Scalar>
78 RCP<const VectorSpaceBase<Scalar> > get_x_space()
const 79 { TEUCHOS_TEST_FOR_EXCEPT(
true);
return Teuchos::null; }
80 RCP<const VectorSpaceBase<Scalar> > get_f_space()
const 81 { TEUCHOS_TEST_FOR_EXCEPT(
true);
return Teuchos::null; }
82 ModelEvaluatorBase::InArgs<Scalar> createInArgs()
const 83 { TEUCHOS_TEST_FOR_EXCEPT(
true);
return ModelEvaluatorBase::InArgs<Scalar>(); }
84 ModelEvaluatorBase::OutArgs<Scalar> createOutArgs()
const 85 { TEUCHOS_TEST_FOR_EXCEPT(
true);
return ModelEvaluatorBase::OutArgs<Scalar>(); }
87 const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
88 const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
90 { TEUCHOS_TEST_FOR_EXCEPT(
true); }
92 static ModelEvaluatorBase::OutArgs<Scalar> createOutArgs(
93 const ModelEvaluator<Scalar> &model,
94 const SelectedDerivatives &fdDerivatives
100 const MEB::OutArgs<Scalar> wrappedOutArgs = model.createOutArgs();
101 const int Np = wrappedOutArgs.Np(),
Ng = wrappedOutArgs.Ng();
102 MEB::OutArgsSetup<Scalar> outArgs;
104 outArgs.setModelEvalDescription(
105 "DirectionalFiniteDiffCalculator: " + model.description()
111 outArgs.set_Np_Ng(Np,Ng);
112 outArgs.setSupports(wrappedOutArgs);
116 const SelectedDerivatives::supports_DfDp_t
117 &supports_DfDp = fdDerivatives.supports_DfDp_;
119 SelectedDerivatives::supports_DfDp_t::const_iterator
120 itr = supports_DfDp.begin();
121 itr != supports_DfDp.end();
126 assert_p_space(model,l);
127 outArgs.setSupports(MEB::OUT_ARG_DfDp,l,MEB::DERIV_MV_BY_COL);
132 const SelectedDerivatives::supports_DgDp_t
133 &supports_DgDp = fdDerivatives.supports_DgDp_;
135 SelectedDerivatives::supports_DgDp_t::const_iterator
136 itr = supports_DgDp.begin();
137 itr != supports_DgDp.end();
141 const int j = itr->first;
142 const int l = itr->second;
143 assert_p_space(model,l);
144 outArgs.setSupports(MEB::OUT_ARG_DgDp,j,l,MEB::DERIV_MV_BY_COL);
153 static void assert_p_space(
const ModelEvaluator<Scalar> &model,
const int l )
156 const bool p_space_l_is_in_core = model.get_p_space(l)->hasInCoreView();
157 TEUCHOS_TEST_FOR_EXCEPTION(
158 !p_space_l_is_in_core, std::logic_error,
159 "Error, for the model " << model.description()
160 <<
", the space p_space("<<l<<
") must be in-core so that they can" 161 " act as the domain space for the multi-vector derivative!" 175 template<
class Scalar>
176 const std::string& DirectionalFiniteDiffCalculator<Scalar>::FDMethod_name()
178 static std::string loc_FDMethod_name =
"FD Method";
179 return loc_FDMethod_name;
183 template<
class Scalar>
185 Teuchos::StringToIntegralParameterEntryValidator<
186 Thyra::DirectionalFiniteDiffCalculatorTypes::EFDMethodType
189 DirectionalFiniteDiffCalculator<Scalar>::fdMethodValidator()
191 static RCP<Teuchos::StringToIntegralParameterEntryValidator<EFDMethodType> >
192 loc_fdMethodValidator
194 new Teuchos::StringToIntegralParameterEntryValidator<Thyra::DirectionalFiniteDiffCalculatorTypes::EFDMethodType>(
195 Teuchos::tuple<std::string>(
201 ,
"order-four-central" 204 ,Teuchos::tuple<std::string>(
205 "Use O(eps) one sided finite differences (cramped bounds)" 206 ,
"Use O(eps^2) one sided finite differences (cramped bounds)" 207 ,
"Use O(eps^2) two sided central finite differences" 208 ,
"Use \"order-two-central\" when not cramped by bounds, otherwise use \"order-two\"" 209 ,
"Use O(eps^4) one sided finite differences (cramped bounds)" 210 ,
"Use O(eps^4) two sided central finite differences" 211 ,
"Use \"order-four-central\" when not cramped by bounds, otherwise use \"order-four\"" 213 ,Teuchos::tuple<Thyra::DirectionalFiniteDiffCalculatorTypes::EFDMethodType>(
214 Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_ONE
215 ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_TWO
216 ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_TWO_CENTRAL
217 ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_TWO_AUTO
218 ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_FOUR
219 ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_FOUR_CENTRAL
220 ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_FOUR_AUTO
225 return loc_fdMethodValidator;
229 template<
class Scalar>
231 DirectionalFiniteDiffCalculator<Scalar>::FDMethod_default()
233 static std::string loc_FDMethod_default =
"order-one";
234 return loc_FDMethod_default;
238 template<
class Scalar>
240 DirectionalFiniteDiffCalculator<Scalar>::FDStepSelectType_name()
242 static std::string loc_FDStepSelectType_name =
"FD Step Select Type";
243 return loc_FDStepSelectType_name;
247 template<
class Scalar>
249 Teuchos::StringToIntegralParameterEntryValidator<
250 Thyra::DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType
253 DirectionalFiniteDiffCalculator<Scalar>::fdStepSelectTypeValidator()
255 static const RCP<Teuchos::StringToIntegralParameterEntryValidator<EFDStepSelectType> >
256 loc_fdStepSelectTypeValidator
258 new Teuchos::StringToIntegralParameterEntryValidator<Thyra::DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType>(
259 Teuchos::tuple<std::string>(
263 ,Teuchos::tuple<std::string>(
264 "Use absolute step size \""+FDStepLength_name()+
"\"" 265 ,
"Use relative step size \""+FDStepLength_name()+
"\"*||xo||inf" 267 ,Teuchos::tuple<Thyra::DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType>(
268 Thyra::DirectionalFiniteDiffCalculatorTypes::FD_STEP_ABSOLUTE
269 ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_STEP_RELATIVE
274 return loc_fdStepSelectTypeValidator;
278 template<
class Scalar>
280 DirectionalFiniteDiffCalculator<Scalar>::FDStepSelectType_default()
282 static std::string loc_FDStepSelectType_default =
"Absolute";
283 return loc_FDStepSelectType_default;
287 template<
class Scalar>
289 DirectionalFiniteDiffCalculator<Scalar>::FDStepLength_name()
291 static std::string loc_FDStepLength_name =
"FD Step Length";
292 return loc_FDStepLength_name;
296 template<
class Scalar>
298 DirectionalFiniteDiffCalculator<Scalar>::FDStepLength_default()
300 static double loc_FDStepLength_default = -1.0;
301 return loc_FDStepLength_default;
308 template<
class Scalar>
315 :fd_method_type_(fd_method_type_in)
316 ,fd_step_select_type_(fd_step_select_type_in)
317 ,fd_step_size_(fd_step_size_in)
318 ,fd_step_size_min_(fd_step_size_min_in)
325 template<
class Scalar>
327 RCP<ParameterList>
const& paramList
330 TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==0);
331 paramList->validateParameters(*getValidParameters());
332 paramList_ = paramList;
333 fd_method_type_ = fdMethodValidator()->getIntegralValue(
334 *paramList_, FDMethod_name(), FDMethod_default());
335 fd_step_select_type_ = fdStepSelectTypeValidator()->getIntegralValue(
336 *paramList_, FDStepSelectType_name(), FDStepSelectType_default());
337 fd_step_size_ = paramList_->get(
338 FDStepLength_name(),FDStepLength_default());
339 Teuchos::readVerboseObjectSublist(&*paramList_,
this);
343 template<
class Scalar>
351 template<
class Scalar>
355 RCP<ParameterList> _paramList = paramList_;
356 paramList_ = Teuchos::null;
361 template<
class Scalar>
362 RCP<const ParameterList>
369 template<
class Scalar>
370 RCP<const ParameterList>
373 using Teuchos::rcp_implicit_cast;
374 typedef Teuchos::ParameterEntryValidator PEV;
375 static RCP<ParameterList> pl;
377 pl = Teuchos::parameterList();
379 FDMethod_name(), FDMethod_default(),
380 "The method used to compute the finite differences.",
381 rcp_implicit_cast<const PEV>(fdMethodValidator())
384 FDStepSelectType_name(), FDStepSelectType_default(),
385 "Method used to select the finite difference step length.",
386 rcp_implicit_cast<const PEV>(fdStepSelectTypeValidator())
389 FDStepLength_name(), FDStepLength_default()
390 ,
"The length of the finite difference step to take.\n" 391 "A value of < 0.0 means that the step length will be determined automatically." 393 Teuchos::setupVerboseObjectSublist(&*pl);
399 template<
class Scalar>
406 return DirectionalFiniteDiffCalculatorTypes::OutArgsCreator<Scalar>::createOutArgs(
407 model, fdDerivatives );
411 template<
class Scalar>
423 THYRA_FUNC_TIME_MONITOR(
424 string(
"Thyra::DirectionalFiniteDiffCalculator<")+ST::name()+
">::calcVariations(...)" 432 namespace DFDCT = DirectionalFiniteDiffCalculatorTypes;
433 typedef RCP<VectorBase<Scalar> > VectorPtr;
435 RCP<Teuchos::FancyOStream> out = this->getOStream();
436 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
437 const bool trace = (
static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM));
438 Teuchos::OSTab tab(out);
440 if(out.get() && trace)
441 *out <<
"\nEntering DirectionalFiniteDiffCalculator<Scalar>::calcVariations(...)\n";
443 if(out.get() && trace)
445 <<
"\nbasePoint=\n" << describe(bp,verbLevel)
446 <<
"\ndirections=\n" << describe(dir,verbLevel)
447 <<
"\nbaseFunctionValues=\n" << describe(bfunc,verbLevel)
448 <<
"\nvariations=\n" << describe(var,Teuchos::VERB_LOW);
451 TEUCHOS_TEST_FOR_EXCEPTION(
452 var.
isEmpty(), std::logic_error,
453 "Error, all of the variations can not be null!" 485 switch(this->fd_method_type()) {
486 case DFDCT::FD_ORDER_ONE:
487 if(out.get()&&trace) *out<<
"\nUsing one-sided, first-order finite differences ...\n";
489 case DFDCT::FD_ORDER_TWO:
490 if(out.get()&&trace) *out<<
"\nUsing one-sided, second-order finite differences ...\n";
492 case DFDCT::FD_ORDER_TWO_CENTRAL:
493 if(out.get()&&trace) *out<<
"\nUsing second-order central finite differences ...\n";
495 case DFDCT::FD_ORDER_TWO_AUTO:
496 if(out.get()&&trace) *out<<
"\nUsing auto selection of some second-order finite difference method ...\n";
498 case DFDCT::FD_ORDER_FOUR:
499 if(out.get()&&trace) *out<<
"\nUsing one-sided, fourth-order finite differences ...\n";
501 case DFDCT::FD_ORDER_FOUR_CENTRAL:
502 if(out.get()&&trace) *out<<
"\nUsing fourth-order central finite differences ...\n";
504 case DFDCT::FD_ORDER_FOUR_AUTO:
505 if(out.get()&&trace) *out<<
"\nUsing auto selection of some fourth-order finite difference method ...\n";
508 TEUCHOS_TEST_FOR_EXCEPT(
true);
519 sqrt_epsilon = SMT::squareroot(SMT::eps()),
520 u_optimal_1 = sqrt_epsilon,
521 u_optimal_2 = SMT::squareroot(sqrt_epsilon),
522 u_optimal_4 = SMT::squareroot(u_optimal_2);
525 bp_norm = SMT::zero();
530 switch(this->fd_method_type()) {
531 case DFDCT::FD_ORDER_ONE:
532 uh_opt = u_optimal_1 * ( fd_step_select_type() == DFDCT::FD_STEP_ABSOLUTE ? 1.0 : bp_norm + 1.0 );
534 case DFDCT::FD_ORDER_TWO:
535 case DFDCT::FD_ORDER_TWO_CENTRAL:
536 case DFDCT::FD_ORDER_TWO_AUTO:
537 uh_opt = u_optimal_2 * ( fd_step_select_type() == DFDCT::FD_STEP_ABSOLUTE ? 1.0 : bp_norm + 1.0 );
539 case DFDCT::FD_ORDER_FOUR:
540 case DFDCT::FD_ORDER_FOUR_CENTRAL:
541 case DFDCT::FD_ORDER_FOUR_AUTO:
542 uh_opt = u_optimal_4 * ( fd_step_select_type() == DFDCT::FD_STEP_ABSOLUTE ? 1.0 : bp_norm + 1.0 );
545 TEUCHOS_TEST_FOR_EXCEPT(
true);
548 if(out.get()&&trace) *out<<
"\nDefault optimal step length uh_opt = " << uh_opt <<
" ...\n";
555 uh = this->fd_step_size();
559 else if( fd_step_select_type() == DFDCT::FD_STEP_RELATIVE )
560 uh *= (bp_norm + 1.0);
562 if(out.get()&&trace) *out<<
"\nStep size to be used uh="<<uh<<
"\n";
576 DFDCT::EFDMethodType l_fd_method_type = this->fd_method_type();
577 switch(l_fd_method_type) {
578 case DFDCT::FD_ORDER_TWO_AUTO:
579 l_fd_method_type = DFDCT::FD_ORDER_TWO_CENTRAL;
581 case DFDCT::FD_ORDER_FOUR_AUTO:
582 l_fd_method_type = DFDCT::FD_ORDER_FOUR_CENTRAL;
592 p_saved = out->precision();
597 const int Np = var.
Np(), Ng = var.
Ng();
600 VectorPtr per_x, per_x_dot;
601 std::vector<VectorPtr> per_p(Np);
605 if( dir.
get_x().get() )
608 pp.set_x(bp.
get_x());
610 if( bp.
supports(MEB::IN_ARG_x_dot) ) {
612 pp.set_x_dot(per_x_dot=createMember(model.
get_x_space()));
616 for(
int l = 0; l < Np; ++l ) {
617 if( dir.
get_p(l).get() )
618 pp.set_p(l,per_p[l]=createMember(model.
get_p_space(l)));
620 pp.set_p(l,bp.
get_p(l));
622 if(out.get() && trace)
624 <<
"\nperturbedPoint after initial setup (with some unintialized vectors) =\n" 625 << describe(pp,verbLevel);
628 bool all_funcs_at_base_computed =
true;
632 if( var.
supports(MEB::OUT_ARG_f) && (f=var.
get_f()).
get() ) {
634 assign(f.ptr(),ST::zero());
635 if(!bfunc.
get_f().get()) all_funcs_at_base_computed =
false;
637 for(
int j = 0; j < Ng; ++j ) {
639 if( (g_j=var.
get_g(j)).
get() ) {
641 assign(g_j.ptr(),ST::zero());
642 if(!bfunc.
get_g(j).get()) all_funcs_at_base_computed =
false;
646 if(out.get() && trace)
648 <<
"\nperturbedFunctions after initial setup (with some unintialized vectors) =\n" 649 << describe(pfunc,verbLevel);
651 const int dbl_p = 15;
653 *out << std::setprecision(dbl_p);
661 switch(l_fd_method_type) {
662 case DFDCT::FD_ORDER_ONE:
666 case DFDCT::FD_ORDER_TWO:
670 case DFDCT::FD_ORDER_TWO_CENTRAL:
674 case DFDCT::FD_ORDER_FOUR:
678 case DFDCT::FD_ORDER_FOUR_CENTRAL:
683 TEUCHOS_TEST_FOR_EXCEPT(
true);
685 for(
int eval_i = 1; eval_i <= num_evals; ++eval_i ) {
690 switch(l_fd_method_type) {
691 case DFDCT::FD_ORDER_ONE: {
704 case DFDCT::FD_ORDER_TWO: {
721 case DFDCT::FD_ORDER_TWO_CENTRAL: {
734 case DFDCT::FD_ORDER_FOUR: {
759 case DFDCT::FD_ORDER_FOUR_CENTRAL: {
780 case DFDCT::FD_ORDER_TWO_AUTO:
781 case DFDCT::FD_ORDER_FOUR_AUTO:
784 TEUCHOS_TEST_FOR_EXCEPT(
true);
787 if(out.get() && trace)
788 *out <<
"\neval_i="<<eval_i<<
", uh_i="<<uh_i<<
", wgt_i="<<wgt_i<<
"\n";
789 Teuchos::OSTab tab2(out);
792 if(uh_i == ST::zero()) {
793 MEB::OutArgs<Scalar> bfuncall;
794 if(!all_funcs_at_base_computed) {
797 if( pfunc.supports(MEB::OUT_ARG_f) && pfunc.get_f().get() && !bfunc.
get_f().get() ) {
800 for(
int j = 0; j < Ng; ++j ) {
801 if( pfunc.get_g(j).get() && !bfunc.
get_g(j).get() ) {
802 bfuncall.set_g(j,createMember(model.
get_g_space(j)));
806 bfuncall.setArgs(bfunc);
812 if(out.get() && trace)
813 *out <<
"\nSetting variations = wgt_i * basePoint ...\n";
815 if( pfunc.supports(MEB::OUT_ARG_f) && (f=var.
get_f()).
get() ) {
816 V_StV<Scalar>(f.ptr(), wgt_i, *bfuncall.get_f());
818 for(
int j = 0; j < Ng; ++j ) {
820 if( (g_j=var.
get_g(j)).
get() ) {
821 V_StV<Scalar>(g_j.ptr(), wgt_i, *bfuncall.get_g(j));
826 if(out.get() && trace)
827 *out <<
"\nSetting perturbedPoint = basePoint + uh_i*uh*direction ...\n";
831 V_StVpV(per_x.ptr(),as<Scalar>(uh_i*uh),*dir.
get_x(),*bp.
get_x());
834 for (
int l = 0; l < Np; ++l ) {
835 if( dir.
get_p(l).get() )
836 V_StVpV(per_p[l].ptr(), as<Scalar>(uh_i*uh), *dir.
get_p(l), *bp.
get_p(l));
839 if(out.get() && trace)
840 *out <<
"\nperturbedPoint=\n" << describe(pp,verbLevel);
842 if(out.get() && trace)
843 *out <<
"\nCompute perturbedFunctions at perturbedPoint...\n";
845 if(out.get() && trace)
846 *out <<
"\nperturbedFunctions=\n" << describe(pfunc,verbLevel);
850 if(out.get() && trace)
851 *out <<
"\nComputing variations += wgt_i*perturbedfunctions ...\n";
853 if( pfunc.supports(MEB::OUT_ARG_f) && (f=pfunc.get_f()).
get() )
854 Vp_StV<Scalar>(var.
get_f().ptr(), wgt_i, *f);
855 for(
int j = 0; j < Ng; ++j ) {
857 if( (g_j=pfunc.get_g(j)).
get() )
858 Vp_StV<Scalar>(var.
get_g(j).ptr(), wgt_i, *g_j);
862 if(out.get() && trace)
863 *out <<
"\nvariations=\n" << describe(var,verbLevel);
872 const Scalar alpha = ST::one()/(dwgt*uh);
873 if(out.get() && trace)
875 <<
"\nComputing variations *= (1.0)/(dwgt*uh)," 876 <<
" where (1.0)/(dwgt*uh) = (1.0)/("<<dwgt<<
"*"<<uh<<
") = "<<alpha<<
" ...\n";
880 for(
int j = 0; j < Ng; ++j ) {
882 if( (g_j=var.
get_g(j)).
get() )
883 Vt_S(g_j.ptr(),alpha);
885 if(out.get() && trace)
886 *out <<
"\nFinal variations=\n" << describe(var,verbLevel);
890 *out << std::setprecision(p_saved);
892 if(out.get() && trace)
893 *out <<
"\nLeaving DirectionalFiniteDiffCalculator<Scalar>::calcVariations(...)\n";
898 template<
class Scalar>
910 THYRA_FUNC_TIME_MONITOR(
911 string(
"Thyra::DirectionalFiniteDiffCalculator<")+ST::name()+
">::calcDerivatives(...)" 915 typedef RCP<VectorBase<Scalar> > VectorPtr;
916 typedef RCP<MultiVectorBase<Scalar> > MultiVectorPtr;
918 RCP<Teuchos::FancyOStream> out = this->getOStream();
919 Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
920 const bool trace = (
static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM));
921 Teuchos::OSTab tab(out);
923 if(out.get() && trace)
924 *out <<
"\nEntering DirectionalFiniteDiffCalculator<Scalar>::calcDerivatives(...)\n";
926 if(out.get() && trace)
928 <<
"\nbasePoint=\n" << describe(bp,verbLevel)
929 <<
"\nbaseFunctionValues=\n" << describe(bfunc,verbLevel)
930 <<
"\nderivatives=\n" << describe(deriv,Teuchos::VERB_LOW);
935 const int Np = bp.
Np(), Ng = bfunc.
Ng();
938 MultiVectorPtr DfDp_l;
939 std::vector<MEB::DerivativeMultiVector<Scalar> > DgDp_l(Ng);
940 std::vector<VectorPtr> var_g(Ng);
941 for(
int l = 0; l < Np; ++l ) {
942 if(out.get() && trace)
943 *out <<
"\nComputing derivatives for parameter subvector p("<<l<<
") ...\n";
944 Teuchos::OSTab tab2(out);
948 bool hasDerivObject =
false;
951 !deriv.
supports(MEB::OUT_ARG_DfDp,l).none()
955 hasDerivObject =
true;
956 std::ostringstream name; name <<
"DfDp("<<l<<
")";
957 DfDp_l = get_mv(deriv.
get_DfDp(l),name.str(),MEB::DERIV_MV_BY_COL);
960 DfDp_l = Teuchos::null;
963 for (
int j = 0; j < Ng; ++j ) {
965 !deriv.
supports(MEB::OUT_ARG_DgDp,j,l).none()
970 hasDerivObject =
true;
971 std::ostringstream name; name <<
"DgDp("<<j<<
","<<l<<
")";
972 DgDp_l[j] = get_dmv(deriv.
get_DgDp(j,l),name.str());
973 if( DgDp_l[j].getMultiVector().get() && !var_g[j].get() )
980 DgDp_l[j] = MEB::DerivativeMultiVector<Scalar>();
981 var_g[j] = Teuchos::null;
987 if (hasDerivObject) {
988 VectorPtr e_i = createMember(model.
get_p_space(l));
990 assign(e_i.ptr(),ST::zero());
991 const int np_l = e_i->space()->dim();
992 for(
int i = 0 ; i < np_l; ++ i ) {
993 if(out.get() && trace)
994 *out <<
"\nComputing derivatives for single variable p("<<l<<
")("<<i<<
") ...\n";
995 Teuchos::OSTab tab3(out);
996 if(DfDp_l.get()) var.set_f(DfDp_l->col(i));
997 for(
int j = 0; j < Ng; ++j) {
998 MultiVectorPtr DgDp_j_l;
999 if( (DgDp_j_l=DgDp_l[j].getMultiVector()).get() ) {
1000 var.set_g(j,var_g[j]);
1003 set_ele(i,ST::one(),e_i.ptr());
1004 this->calcVariations(
1005 model,bp,dir,bfunc,var
1007 set_ele(i,ST::zero(),e_i.ptr());
1008 if (DfDp_l.get()) var.set_f(Teuchos::null);
1009 for (
int j = 0; j < Ng; ++j) {
1010 MultiVectorPtr DgDp_j_l;
1011 if ( !is_null(DgDp_j_l=DgDp_l[j].getMultiVector()) ) {
1012 assign( DgDp_j_l->col(i).ptr(), *var_g[j] );
1016 dir.set_p(l,Teuchos::null);
1020 if(out.get() && trace)
1022 <<
"\nderivatives=\n" << describe(deriv,verbLevel);
1024 if(out.get() && trace)
1025 *out <<
"\nLeaving DirectionalFiniteDiffCalculator<Scalar>::calcDerivatives(...)\n";
1033 #endif // THYRA_DIRECTIONAL_FINITE_DIFF_CALCULATOR_DEF_HPP virtual ModelEvaluatorBase::OutArgs< Scalar > createOutArgs() const =0
Create an empty output functions/derivatives object that can be set up and passed to evalModel()...
virtual RCP< const VectorSpaceBase< Scalar > > get_x_space() const =0
Return the vector space for the state variables x <: RE^n_x.
ModelEvaluatorBase::OutArgs< Scalar > createOutArgs(const ModelEvaluator< Scalar > &model, const SelectedDerivatives &fdDerivatives)
Create an augmented out args object for holding finite difference objects.
Pure abstract base interface for evaluating a stateless "model" that can be mapped into a number of d...
Simple utility class used to select finite difference derivatives for OutArgs object.
void setParameterList(RCP< ParameterList > const ¶mList)
void calcDerivatives(const ModelEvaluator< Scalar > &model, const ModelEvaluatorBase::InArgs< Scalar > &basePoint, const ModelEvaluatorBase::OutArgs< Scalar > &baseFunctionValues, const ModelEvaluatorBase::OutArgs< Scalar > &derivatives) const
Compute entire derivative objects using finite differences.
DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType EFDStepSelectType
bool supports(EInArgsMembers arg) const
Determines if an input argument is supported or not.
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object...
Derivative< Scalar > get_DfDp(int l) const
Precondition: supports(OUT_ARG_DfDp,l)==true.
Derivative< Scalar > get_DgDp(int j, int l) const
Precondition: supports(OUT_ARG_DgDp,j,l)==true.
virtual void evalModel(const ModelEvaluatorBase::InArgs< Scalar > &inArgs, const ModelEvaluatorBase::OutArgs< Scalar > &outArgs) const =0
Compute all of the requested functions/derivatives at the given point.
RCP< const VectorBase< Scalar > > get_x() const
Precondition: supports(IN_ARG_x)==true.
virtual ModelEvaluatorBase::InArgs< Scalar > createInArgs() const =0
Create an empty input arguments object that can be set up and passed to evalModel().
RCP< const VectorBase< Scalar > > get_p(int l) const
Get p(l) where 0 <= l && l < this->Np().
DirectionalFiniteDiffCalculatorTypes::EFDMethodType EFDMethodType
RCP< const VectorBase< Scalar > > get_x_dot() const
Precondition: supports(IN_ARG_x_dot)==true.
bool supports(EOutArgsMembers arg) const
Determine if an input argument is supported or not.
virtual RCP< const VectorSpaceBase< Scalar > > get_p_space(int l) const =0
Return the vector space for the auxiliary parameters p(l) <: RE^n_p_l.
Evaluation< VectorBase< Scalar > > get_g(int j) const
Precondition: supports(OUT_ARG_g)==true..
int Ng() const
Return the number of axillary response functions g(j)(...) supported (Ng >= 0).
Base subclass for ModelEvaluator that defines some basic types.
virtual RCP< const VectorSpaceBase< Scalar > > get_f_space() const =0
Return the vector space for the state function f(...) <: RE^n_x.
DirectionalFiniteDiffCalculator(EFDMethodType fd_method_type=DirectionalFiniteDiffCalculatorTypes::FD_ORDER_FOUR_AUTO, EFDStepSelectType fd_step_select_type=DirectionalFiniteDiffCalculatorTypes::FD_STEP_ABSOLUTE, ScalarMag fd_step_size=-1.0, ScalarMag fd_step_size_min=-1.0)
void calcVariations(const ModelEvaluator< Scalar > &model, const ModelEvaluatorBase::InArgs< Scalar > &basePoint, const ModelEvaluatorBase::InArgs< Scalar > &directions, const ModelEvaluatorBase::OutArgs< Scalar > &baseFunctionValues, const ModelEvaluatorBase::OutArgs< Scalar > &variations) const
Compute variations using directional finite differences..
int Np() const
Return the number of parameter subvectors p(l) supported (Np >= 0).
RCP< const ParameterList > getParameterList() const
virtual RCP< const VectorSpaceBase< Scalar > > get_g_space(int j) const =0
Return the vector space for the auxiliary response functions g(j) <: RE^n_g_j.
ST::magnitudeType ScalarMag
Evaluation< VectorBase< Scalar > > get_f() const
Precondition: supports(OUT_ARG_f)==true.
RCP< const ParameterList > getValidParameters() const
RCP< ParameterList > unsetParameterList()
ModelEvaluatorBase()
constructor
RCP< ParameterList > getNonconstParameterList()
int Np() const
Return the number of parameter subvectors p(l) supported (Np >= 0).
Concrete aggregate class for all output arguments computable by a ModelEvaluator subclass object...