Thyra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Thyra_DirectionalFiniteDiffCalculator_def.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Thyra: Interfaces and Support for Abstract Numerical Algorithms
5 // Copyright (2004) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef THYRA_DIRECTIONAL_FINITE_DIFF_CALCULATOR_DEF_HPP
43 #define THYRA_DIRECTIONAL_FINITE_DIFF_CALCULATOR_DEF_HPP
44 
45 
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"
55 
56 
57 namespace Thyra {
58 
59 
60 namespace DirectionalFiniteDiffCalculatorTypes {
61 
62 
63 //
64 // Undocumented utility class for setting support for new derivative objects
65 // on an OutArgs object! Warning, users should not attempt to play these
66 // tricks on their own!
67 //
68 // Note that because of the design of the OutArgs and OutArgsSetup classes,
69 // you can only change the list of supported arguments in a subclass of
70 // ModelEvaluatorBase since OutArgsSetup is a protected type. The fact that
71 // the only way to do this is convoluted is a feature!
72 //
73 template<class Scalar>
74 class OutArgsCreator : public StateFuncModelEvaluatorBase<Scalar>
75 {
76 public:
77  // Public functions overridden from ModelEvaulator.
78  RCP<const VectorSpaceBase<Scalar> > get_x_space() const
79  { TEUCHOS_TEST_FOR_EXCEPT(true); TEUCHOS_UNREACHABLE_RETURN(Teuchos::null); }
80  RCP<const VectorSpaceBase<Scalar> > get_f_space() const
81  { TEUCHOS_TEST_FOR_EXCEPT(true); TEUCHOS_UNREACHABLE_RETURN(Teuchos::null); }
82  ModelEvaluatorBase::InArgs<Scalar> createInArgs() const
83  { TEUCHOS_TEST_FOR_EXCEPT(true); TEUCHOS_UNREACHABLE_RETURN(ModelEvaluatorBase::InArgs<Scalar>()); }
84  ModelEvaluatorBase::OutArgs<Scalar> createOutArgs() const
85  { TEUCHOS_TEST_FOR_EXCEPT(true); TEUCHOS_UNREACHABLE_RETURN(ModelEvaluatorBase::OutArgs<Scalar>()); }
86  void evalModel(
87  const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
88  const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
89  ) const
90  { TEUCHOS_TEST_FOR_EXCEPT(true); }
91  // Static function that does the magic!
92  static ModelEvaluatorBase::OutArgs<Scalar> createOutArgs(
93  const ModelEvaluator<Scalar> &model,
94  const SelectedDerivatives &fdDerivatives
95  )
96  {
97 
98  typedef ModelEvaluatorBase MEB;
99 
100  const MEB::OutArgs<Scalar> wrappedOutArgs = model.createOutArgs();
101  const int Np = wrappedOutArgs.Np(), Ng = wrappedOutArgs.Ng();
102  MEB::OutArgsSetup<Scalar> outArgs;
103 
104  outArgs.setModelEvalDescription(
105  "DirectionalFiniteDiffCalculator: " + model.description()
106  );
107 
108  // Start off by supporting everything that the underlying model supports
109  // computing!
110 
111  outArgs.set_Np_Ng(Np,Ng);
112  outArgs.setSupports(wrappedOutArgs);
113 
114  // Add support for finite difference DfDp(l) if asked
115 
116  const SelectedDerivatives::supports_DfDp_t
117  &supports_DfDp = fdDerivatives.supports_DfDp_;
118  for(
119  SelectedDerivatives::supports_DfDp_t::const_iterator
120  itr = supports_DfDp.begin();
121  itr != supports_DfDp.end();
122  ++itr
123  )
124  {
125  const int l = *itr;
126  assert_p_space(model,l);
127  outArgs.setSupports(MEB::OUT_ARG_DfDp,l,MEB::DERIV_MV_BY_COL);
128  }
129 
130  // Add support for finite difference DgDp(j,l) if asked
131 
132  const SelectedDerivatives::supports_DgDp_t
133  &supports_DgDp = fdDerivatives.supports_DgDp_;
134  for(
135  SelectedDerivatives::supports_DgDp_t::const_iterator
136  itr = supports_DgDp.begin();
137  itr != supports_DgDp.end();
138  ++itr
139  )
140  {
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);
145  }
146 
147  return outArgs;
148 
149  }
150 
151 private:
152 
153  static void assert_p_space( const ModelEvaluator<Scalar> &model, const int l )
154  {
155 #ifdef TEUCHOS_DEBUG
156  const bool p_space_l_is_in_core = model.get_p_space(l)->hasInCoreView();
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!"
162  );
163 #else
164  (void)model;
165  (void)l;
166 #endif
167  }
168 
169 };
170 
171 
172 } // namespace DirectionalFiniteDiffCalculatorTypes
173 
174 
175 // Private static data members
176 
177 
178 template<class Scalar>
179 const std::string& DirectionalFiniteDiffCalculator<Scalar>::FDMethod_name()
180 {
181  static std::string loc_FDMethod_name = "FD Method";
182  return loc_FDMethod_name;
183 }
184 
185 
186 template<class Scalar>
187 const RCP<
189  Thyra::DirectionalFiniteDiffCalculatorTypes::EFDMethodType
190  >
191 >&
192 DirectionalFiniteDiffCalculator<Scalar>::fdMethodValidator()
193 {
194  static RCP<Teuchos::StringToIntegralParameterEntryValidator<EFDMethodType> >
195  loc_fdMethodValidator
196  = Teuchos::rcp(
198  Teuchos::tuple<std::string>(
199  "order-one"
200  ,"order-two"
201  ,"order-two-central"
202  ,"order-two-auto"
203  ,"order-four"
204  ,"order-four-central"
205  ,"order-four-auto"
206  )
207  ,Teuchos::tuple<std::string>(
208  "Use O(eps) one sided finite differences (cramped bounds)"
209  ,"Use O(eps^2) one sided finite differences (cramped bounds)"
210  ,"Use O(eps^2) two sided central finite differences"
211  ,"Use \"order-two-central\" when not cramped by bounds, otherwise use \"order-two\""
212  ,"Use O(eps^4) one sided finite differences (cramped bounds)"
213  ,"Use O(eps^4) two sided central finite differences"
214  ,"Use \"order-four-central\" when not cramped by bounds, otherwise use \"order-four\""
215  )
216  ,Teuchos::tuple<Thyra::DirectionalFiniteDiffCalculatorTypes::EFDMethodType>(
217  Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_ONE
218  ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_TWO
219  ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_TWO_CENTRAL
220  ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_TWO_AUTO
221  ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_FOUR
222  ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_FOUR_CENTRAL
223  ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_FOUR_AUTO
224  )
225  ,""
226  )
227  );
228  return loc_fdMethodValidator;
229 }
230 
231 
232 template<class Scalar>
233 const std::string&
234 DirectionalFiniteDiffCalculator<Scalar>::FDMethod_default()
235 {
236  static std::string loc_FDMethod_default = "order-one";
237  return loc_FDMethod_default;
238 }
239 
240 
241 template<class Scalar>
242 const std::string&
243 DirectionalFiniteDiffCalculator<Scalar>::FDStepSelectType_name()
244 {
245  static std::string loc_FDStepSelectType_name = "FD Step Select Type";
246  return loc_FDStepSelectType_name;
247 }
248 
249 
250 template<class Scalar>
251 const RCP<
253  Thyra::DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType
254  >
255  >&
256 DirectionalFiniteDiffCalculator<Scalar>::fdStepSelectTypeValidator()
257 {
258  static const RCP<Teuchos::StringToIntegralParameterEntryValidator<EFDStepSelectType> >
259  loc_fdStepSelectTypeValidator
260  = Teuchos::rcp(
262  Teuchos::tuple<std::string>(
263  "Absolute"
264  ,"Relative"
265  )
266  ,Teuchos::tuple<std::string>(
267  "Use absolute step size \""+FDStepLength_name()+"\""
268  ,"Use relative step size \""+FDStepLength_name()+"\"*||xo||inf"
269  )
270  ,Teuchos::tuple<Thyra::DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType>(
271  Thyra::DirectionalFiniteDiffCalculatorTypes::FD_STEP_ABSOLUTE
272  ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_STEP_RELATIVE
273  )
274  ,""
275  )
276  );
277  return loc_fdStepSelectTypeValidator;
278 }
279 
280 
281 template<class Scalar>
282 const std::string&
283 DirectionalFiniteDiffCalculator<Scalar>::FDStepSelectType_default()
284 {
285  static std::string loc_FDStepSelectType_default = "Absolute";
286  return loc_FDStepSelectType_default;
287 }
288 
289 
290 template<class Scalar>
291 const std::string&
292 DirectionalFiniteDiffCalculator<Scalar>::FDStepLength_name()
293 {
294  static std::string loc_FDStepLength_name = "FD Step Length";
295  return loc_FDStepLength_name;
296 }
297 
298 
299 template<class Scalar>
300 const double&
301 DirectionalFiniteDiffCalculator<Scalar>::FDStepLength_default()
302 {
303  static double loc_FDStepLength_default = -1.0;
304  return loc_FDStepLength_default;
305 }
306 
307 
308 // Constructors/initializer
309 
310 
311 template<class Scalar>
313  EFDMethodType fd_method_type_in
314  ,EFDStepSelectType fd_step_select_type_in
315  ,ScalarMag fd_step_size_in
316  ,ScalarMag fd_step_size_min_in
317  )
318  :fd_method_type_(fd_method_type_in)
319  ,fd_step_select_type_(fd_step_select_type_in)
320  ,fd_step_size_(fd_step_size_in)
321  ,fd_step_size_min_(fd_step_size_min_in)
322 {}
323 
324 
325 // Overriden from ParameterListAcceptor
326 
327 
328 template<class Scalar>
330  RCP<ParameterList> const& paramList
331  )
332 {
333  TEUCHOS_TEST_FOR_EXCEPT(paramList.get()==0);
334  paramList->validateParameters(*getValidParameters());
335  paramList_ = paramList;
336  fd_method_type_ = fdMethodValidator()->getIntegralValue(
337  *paramList_, FDMethod_name(), FDMethod_default());
338  fd_step_select_type_ = fdStepSelectTypeValidator()->getIntegralValue(
339  *paramList_, FDStepSelectType_name(), FDStepSelectType_default());
340  fd_step_size_ = paramList_->get(
341  FDStepLength_name(),FDStepLength_default());
342  Teuchos::readVerboseObjectSublist(&*paramList_,this);
343 }
344 
345 
346 template<class Scalar>
349 {
350  return paramList_;
351 }
352 
353 
354 template<class Scalar>
357 {
358  RCP<ParameterList> _paramList = paramList_;
359  paramList_ = Teuchos::null;
360  return _paramList;
361 }
362 
363 
364 template<class Scalar>
367 {
368  return paramList_;
369 }
370 
371 
372 template<class Scalar>
375 {
376  using Teuchos::rcp_implicit_cast;
378  static RCP<ParameterList> pl;
379  if(pl.get()==NULL) {
380  pl = Teuchos::parameterList();
381  pl->set(
382  FDMethod_name(), FDMethod_default(),
383  "The method used to compute the finite differences.",
384  rcp_implicit_cast<const PEV>(fdMethodValidator())
385  );
386  pl->set(
387  FDStepSelectType_name(), FDStepSelectType_default(),
388  "Method used to select the finite difference step length.",
389  rcp_implicit_cast<const PEV>(fdStepSelectTypeValidator())
390  );
391  pl->set(
392  FDStepLength_name(), FDStepLength_default()
393  ,"The length of the finite difference step to take.\n"
394  "A value of < 0.0 means that the step length will be determined automatically."
395  );
396  Teuchos::setupVerboseObjectSublist(&*pl);
397  }
398  return pl;
399 }
400 
401 
402 template<class Scalar>
405  const ModelEvaluator<Scalar> &model,
406  const SelectedDerivatives &fdDerivatives
407  )
408 {
409  return DirectionalFiniteDiffCalculatorTypes::OutArgsCreator<Scalar>::createOutArgs(
410  model, fdDerivatives );
411 }
412 
413 
414 template<class Scalar>
416  const ModelEvaluator<Scalar> &model,
417  const ModelEvaluatorBase::InArgs<Scalar> &bp, // basePoint
418  const ModelEvaluatorBase::InArgs<Scalar> &dir, // directions
419  const ModelEvaluatorBase::OutArgs<Scalar> &bfunc, // baseFunctionValues
420  const ModelEvaluatorBase::OutArgs<Scalar> &var // variations
421  ) const
422 {
423 
424  using std::string;
425 
426  THYRA_FUNC_TIME_MONITOR(
427  string("Thyra::DirectionalFiniteDiffCalculator<")+ST::name()+">::calcVariations(...)"
428  );
429 
430  using std::setw;
431  using std::endl;
432  using std::right;
433  using Teuchos::as;
434  typedef ModelEvaluatorBase MEB;
435  namespace DFDCT = DirectionalFiniteDiffCalculatorTypes;
436  typedef RCP<VectorBase<Scalar> > VectorPtr;
437 
438  RCP<Teuchos::FancyOStream> out = this->getOStream();
439  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
440  const bool trace = (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM));
441  Teuchos::OSTab tab(out);
442 
443  if(out.get() && trace)
444  *out << "\nEntering DirectionalFiniteDiffCalculator<Scalar>::calcVariations(...)\n";
445 
446  if(out.get() && trace)
447  *out
448  << "\nbasePoint=\n" << describe(bp,verbLevel)
449  << "\ndirections=\n" << describe(dir,verbLevel)
450  << "\nbaseFunctionValues=\n" << describe(bfunc,verbLevel)
451  << "\nvariations=\n" << describe(var,Teuchos::VERB_LOW);
452 
453 #ifdef TEUCHOS_DEBUG
455  var.isEmpty(), std::logic_error,
456  "Error, all of the variations can not be null!"
457  );
458 #endif
459 
460  //
461  // To illustrate the theory behind this implementation consider
462  // the generic multi-variable function h(z) <: R^n -> R. Now let's
463  // consider we have the base point zo and the vector v to
464  // perturb h(z) along. First form the function h(zo+a*v) and then
465  // let's compute d(h)/d(a) at a = 0:
466  //
467  // (1) d(h(zo+a*v))/d(a) at a = 0
468  // = sum( d(h)/d(x(i)) * d(x(i))/d(a), i = 1...n)
469  // = sum( d(h)/d(x(i)) * v(i), i = 1...n)
470  // = d(h)/d(a) * v
471  //
472  // Now we can approximate (1) using, for example, central differences as:
473  //
474  // (2) d(h(zo+a*v))/d(a) at a = 0
475  // \approx ( h(zo+h*v) - h(zo+h*v) ) / (2*h)
476  //
477  // If we equate (1) and (2) we have the approximation:
478  //
479  // (3) d(h)/d(a) * v \approx ( h(zo+h*v) - g(zo+h*v) ) / (2*h)
480  //
481  //
482 
483  // /////////////////////////////////////////
484  // Validate the input
485 
486  // ToDo: Validate input!
487 
488  switch(this->fd_method_type()) {
489  case DFDCT::FD_ORDER_ONE:
490  if(out.get()&&trace) *out<<"\nUsing one-sided, first-order finite differences ...\n";
491  break;
492  case DFDCT::FD_ORDER_TWO:
493  if(out.get()&&trace) *out<<"\nUsing one-sided, second-order finite differences ...\n";
494  break;
495  case DFDCT::FD_ORDER_TWO_CENTRAL:
496  if(out.get()&&trace) *out<<"\nUsing second-order central finite differences ...\n";
497  break;
498  case DFDCT::FD_ORDER_TWO_AUTO:
499  if(out.get()&&trace) *out<<"\nUsing auto selection of some second-order finite difference method ...\n";
500  break;
501  case DFDCT::FD_ORDER_FOUR:
502  if(out.get()&&trace) *out<<"\nUsing one-sided, fourth-order finite differences ...\n";
503  break;
504  case DFDCT::FD_ORDER_FOUR_CENTRAL:
505  if(out.get()&&trace) *out<<"\nUsing fourth-order central finite differences ...\n";
506  break;
507  case DFDCT::FD_ORDER_FOUR_AUTO:
508  if(out.get()&&trace) *out<<"\nUsing auto selection of some fourth-order finite difference method ...\n";
509  break;
510  default:
511  TEUCHOS_TEST_FOR_EXCEPT(true); // Should not get here!
512  }
513 
514  // ////////////////////////
515  // Find the step size
516 
517  //
518  // Get defaults for the optimal step sizes
519  //
520 
521  const ScalarMag
522  sqrt_epsilon = SMT::squareroot(SMT::eps()),
523  u_optimal_1 = sqrt_epsilon,
524  u_optimal_2 = SMT::squareroot(sqrt_epsilon),
525  u_optimal_4 = SMT::squareroot(u_optimal_2);
526 
527  ScalarMag
528  bp_norm = SMT::zero();
529  // ToDo above: compute a reasonable norm somehow based on the base-point vector(s)!
530 
531  ScalarMag
532  uh_opt = 0.0;
533  switch(this->fd_method_type()) {
534  case DFDCT::FD_ORDER_ONE:
535  uh_opt = u_optimal_1 * ( fd_step_select_type() == DFDCT::FD_STEP_ABSOLUTE ? 1.0 : bp_norm + 1.0 );
536  break;
537  case DFDCT::FD_ORDER_TWO:
538  case DFDCT::FD_ORDER_TWO_CENTRAL:
539  case DFDCT::FD_ORDER_TWO_AUTO:
540  uh_opt = u_optimal_2 * ( fd_step_select_type() == DFDCT::FD_STEP_ABSOLUTE ? 1.0 : bp_norm + 1.0 );
541  break;
542  case DFDCT::FD_ORDER_FOUR:
543  case DFDCT::FD_ORDER_FOUR_CENTRAL:
544  case DFDCT::FD_ORDER_FOUR_AUTO:
545  uh_opt = u_optimal_4 * ( fd_step_select_type() == DFDCT::FD_STEP_ABSOLUTE ? 1.0 : bp_norm + 1.0 );
546  break;
547  default:
548  TEUCHOS_TEST_FOR_EXCEPT(true); // Should not get here!
549  }
550 
551  if(out.get()&&trace) *out<<"\nDefault optimal step length uh_opt = " << uh_opt << " ...\n";
552 
553  //
554  // Set the step sizes used.
555  //
556 
557  ScalarMag
558  uh = this->fd_step_size();
559 
560  if( uh < 0 )
561  uh = uh_opt;
562  else if( fd_step_select_type() == DFDCT::FD_STEP_RELATIVE )
563  uh *= (bp_norm + 1.0);
564 
565  if(out.get()&&trace) *out<<"\nStep size to be used uh="<<uh<<"\n";
566 
567  //
568  // Find step lengths that stay in bounds!
569  //
570 
571  // ToDo: Add logic for variable bounds when needed!
572 
573  //
574  // Set the actual method being used
575  //
576  // ToDo: Consider cramped bounds and method order.
577  //
578 
579  DFDCT::EFDMethodType l_fd_method_type = this->fd_method_type();
580  switch(l_fd_method_type) {
581  case DFDCT::FD_ORDER_TWO_AUTO:
582  l_fd_method_type = DFDCT::FD_ORDER_TWO_CENTRAL;
583  break;
584  case DFDCT::FD_ORDER_FOUR_AUTO:
585  l_fd_method_type = DFDCT::FD_ORDER_FOUR_CENTRAL;
586  break;
587  default:
588  break; // Okay
589  }
590 
591  //if(out.get()&&trace) *out<<"\nStep size to fit in bounds: uh="<<uh"\n";
592 
593  int p_saved = -1;
594  if(out.get())
595  p_saved = out->precision();
596 
597  // ///////////////////////////////////////////////
598  // Compute the directional derivatives
599 
600  const int Np = var.Np(), Ng = var.Ng();
601 
602  // Setup storage for perturbed variables
603  VectorPtr per_x, per_x_dot, per_x_dot_dot;
604  std::vector<VectorPtr> per_p(Np);
605  MEB::InArgs<Scalar> pp = model.createInArgs();
606  pp.setArgs(bp); // Set all args to start with
607  if( bp.supports(MEB::IN_ARG_x) ) {
608  if( dir.get_x().get() )
609  pp.set_x(per_x=createMember(model.get_x_space()));
610  else
611  pp.set_x(bp.get_x());
612  }
613  if( bp.supports(MEB::IN_ARG_x_dot) ) {
614  if( dir.get_x_dot().get() )
615  pp.set_x_dot(per_x_dot=createMember(model.get_x_space()));
616  else
617  pp.set_x_dot(bp.get_x_dot());
618  }
619  if( bp.supports(MEB::IN_ARG_x_dot_dot) ) {
620  if( dir.get_x_dot_dot().get() )
621  pp.set_x_dot_dot(per_x_dot_dot=createMember(model.get_x_space()));
622  else
623  pp.set_x_dot_dot(bp.get_x_dot_dot());
624  }
625  for( int l = 0; l < Np; ++l ) {
626  if( dir.get_p(l).get() )
627  pp.set_p(l,per_p[l]=createMember(model.get_p_space(l)));
628  else
629  pp.set_p(l,bp.get_p(l));
630  }
631  if(out.get() && trace)
632  *out
633  << "\nperturbedPoint after initial setup (with some uninitialized vectors) =\n"
634  << describe(pp,verbLevel);
635 
636  // Setup storage for perturbed functions
637  bool all_funcs_at_base_computed = true;
638  MEB::OutArgs<Scalar> pfunc = model.createOutArgs();
639  {
640  VectorPtr f;
641  if( var.supports(MEB::OUT_ARG_f) && (f=var.get_f()).get() ) {
642  pfunc.set_f(createMember(model.get_f_space()));
643  assign(f.ptr(),ST::zero());
644  if(!bfunc.get_f().get()) all_funcs_at_base_computed = false;
645  }
646  for( int j = 0; j < Ng; ++j ) {
647  VectorPtr g_j;
648  if( (g_j=var.get_g(j)).get() ) {
649  pfunc.set_g(j,createMember(model.get_g_space(j)));
650  assign(g_j.ptr(),ST::zero());
651  if(!bfunc.get_g(j).get()) all_funcs_at_base_computed = false;
652  }
653  }
654  }
655  if(out.get() && trace)
656  *out
657  << "\nperturbedFunctions after initial setup (with some uninitialized vectors) =\n"
658  << describe(pfunc,verbLevel);
659 
660  const int dbl_p = 15;
661  if(out.get())
662  *out << std::setprecision(dbl_p);
663 
664  //
665  // Compute the weighted sum of the terms
666  //
667 
668  int num_evals = 0;
669  ScalarMag dwgt = SMT::zero();
670  switch(l_fd_method_type) {
671  case DFDCT::FD_ORDER_ONE: // may only need one eval if f(xo) etc is passed in
672  num_evals = 2;
673  dwgt = ScalarMag(1.0);
674  break;
675  case DFDCT::FD_ORDER_TWO: // may only need two evals if c(xo) etc is passed in
676  num_evals = 3;
677  dwgt = ScalarMag(2.0);
678  break;
679  case DFDCT::FD_ORDER_TWO_CENTRAL:
680  num_evals = 2;
681  dwgt = ScalarMag(2.0);
682  break;
683  case DFDCT::FD_ORDER_FOUR:
684  num_evals = 5;
685  dwgt = ScalarMag(12.0);
686  break;
687  case DFDCT::FD_ORDER_FOUR_CENTRAL:
688  num_evals = 4;
689  dwgt = ScalarMag(12.0);
690  break;
691  default:
692  TEUCHOS_TEST_FOR_EXCEPT(true); // Should not get here!
693  }
694  for( int eval_i = 1; eval_i <= num_evals; ++eval_i ) {
695  // Set the step constant and the weighting constant
696  ScalarMag
697  uh_i = 0.0,
698  wgt_i = 0.0;
699  switch(l_fd_method_type) {
700  case DFDCT::FD_ORDER_ONE: {
701  switch(eval_i) {
702  case 1:
703  uh_i = +0.0;
704  wgt_i = -1.0;
705  break;
706  case 2:
707  uh_i = +1.0;
708  wgt_i = +1.0;
709  break;
710  }
711  break;
712  }
713  case DFDCT::FD_ORDER_TWO: {
714  switch(eval_i) {
715  case 1:
716  uh_i = +0.0;
717  wgt_i = -3.0;
718  break;
719  case 2:
720  uh_i = +1.0;
721  wgt_i = +4.0;
722  break;
723  case 3:
724  uh_i = +2.0;
725  wgt_i = -1.0;
726  break;
727  }
728  break;
729  }
730  case DFDCT::FD_ORDER_TWO_CENTRAL: {
731  switch(eval_i) {
732  case 1:
733  uh_i = -1.0;
734  wgt_i = -1.0;
735  break;
736  case 2:
737  uh_i = +1.0;
738  wgt_i = +1.0;
739  break;
740  }
741  break;
742  }
743  case DFDCT::FD_ORDER_FOUR: {
744  switch(eval_i) {
745  case 1:
746  uh_i = +0.0;
747  wgt_i = -25.0;
748  break;
749  case 2:
750  uh_i = +1.0;
751  wgt_i = +48.0;
752  break;
753  case 3:
754  uh_i = +2.0;
755  wgt_i = -36.0;
756  break;
757  case 4:
758  uh_i = +3.0;
759  wgt_i = +16.0;
760  break;
761  case 5:
762  uh_i = +4.0;
763  wgt_i = -3.0;
764  break;
765  }
766  break;
767  }
768  case DFDCT::FD_ORDER_FOUR_CENTRAL: {
769  switch(eval_i) {
770  case 1:
771  uh_i = -2.0;
772  wgt_i = +1.0;
773  break;
774  case 2:
775  uh_i = -1.0;
776  wgt_i = -8.0;
777  break;
778  case 3:
779  uh_i = +1.0;
780  wgt_i = +8.0;
781  break;
782  case 4:
783  uh_i = +2.0;
784  wgt_i = -1.0;
785  break;
786  }
787  break;
788  }
789  case DFDCT::FD_ORDER_TWO_AUTO:
790  case DFDCT::FD_ORDER_FOUR_AUTO:
791  break; // Okay
792  default:
794  }
795 
796  if(out.get() && trace)
797  *out << "\neval_i="<<eval_i<<", uh_i="<<uh_i<<", wgt_i="<<wgt_i<<"\n";
798  Teuchos::OSTab tab2(out);
799 
800  // Compute the weighted term and add it to the sum
801  if(uh_i == ST::zero()) {
802  MEB::OutArgs<Scalar> bfuncall;
803  if(!all_funcs_at_base_computed) {
804  // Compute missing functions at the base point
805  bfuncall = model.createOutArgs();
806  if( pfunc.supports(MEB::OUT_ARG_f) && pfunc.get_f().get() && !bfunc.get_f().get() ) {
807  bfuncall.set_f(createMember(model.get_f_space()));
808  }
809  for( int j = 0; j < Ng; ++j ) {
810  if( pfunc.get_g(j).get() && !bfunc.get_g(j).get() ) {
811  bfuncall.set_g(j,createMember(model.get_g_space(j)));
812  }
813  }
814  model.evalModel(bp,bfuncall);
815  bfuncall.setArgs(bfunc);
816  }
817  else {
818  bfuncall = bfunc;
819  }
820  // Use functions at the base point
821  if(out.get() && trace)
822  *out << "\nSetting variations = wgt_i * basePoint ...\n";
823  VectorPtr f;
824  if( pfunc.supports(MEB::OUT_ARG_f) && (f=var.get_f()).get() ) {
825  V_StV<Scalar>(f.ptr(), wgt_i, *bfuncall.get_f());
826  }
827  for( int j = 0; j < Ng; ++j ) {
828  VectorPtr g_j;
829  if( (g_j=var.get_g(j)).get() ) {
830  V_StV<Scalar>(g_j.ptr(), wgt_i, *bfuncall.get_g(j));
831  }
832  }
833  }
834  else {
835  if(out.get() && trace)
836  *out << "\nSetting perturbedPoint = basePoint + uh_i*uh*direction ...\n";
837  // z = zo + uh_i*uh*v
838  {
839  if ( dir.supports(MEB::IN_ARG_x) && dir.get_x().get() )
840  V_StVpV(per_x.ptr(),as<Scalar>(uh_i*uh),*dir.get_x(),*bp.get_x());
841  if ( dir.supports(MEB::IN_ARG_x_dot) && dir.get_x_dot().get() )
842  V_StVpV(per_x_dot.ptr(), as<Scalar>(uh_i*uh), *dir.get_x_dot(), *bp.get_x_dot());
843  if ( dir.supports(MEB::IN_ARG_x_dot_dot) && dir.get_x_dot_dot().get() )
844  V_StVpV(per_x_dot_dot.ptr(), as<Scalar>(uh_i*uh), *dir.get_x_dot_dot(), *bp.get_x_dot_dot());
845  for ( int l = 0; l < Np; ++l ) {
846  if( dir.get_p(l).get() )
847  V_StVpV(per_p[l].ptr(), as<Scalar>(uh_i*uh), *dir.get_p(l), *bp.get_p(l));
848  }
849  }
850  if(out.get() && trace)
851  *out << "\nperturbedPoint=\n" << describe(pp,verbLevel);
852  // Compute perturbed function values h(zo+uh_i*uh)
853  if(out.get() && trace)
854  *out << "\nCompute perturbedFunctions at perturbedPoint...\n";
855  model.evalModel(pp,pfunc);
856  if(out.get() && trace)
857  *out << "\nperturbedFunctions=\n" << describe(pfunc,verbLevel);
858  // Sum perturbed function values into total variation
859  {
860  // var_h += wgt_i*perturbed_h
861  if(out.get() && trace)
862  *out << "\nComputing variations += wgt_i*perturbedfunctions ...\n";
863  VectorPtr f;
864  if( pfunc.supports(MEB::OUT_ARG_f) && (f=pfunc.get_f()).get() )
865  Vp_StV<Scalar>(var.get_f().ptr(), wgt_i, *f);
866  for( int j = 0; j < Ng; ++j ) {
867  VectorPtr g_j;
868  if( (g_j=pfunc.get_g(j)).get() )
869  Vp_StV<Scalar>(var.get_g(j).ptr(), wgt_i, *g_j);
870  }
871  }
872  }
873  if(out.get() && trace)
874  *out << "\nvariations=\n" << describe(var,verbLevel);
875  }
876 
877  //
878  // Multiply by the scaling factor!
879  //
880 
881  {
882  // var_h *= 1.0/(dwgt*uh)
883  const Scalar alpha = ST::one()/(dwgt*uh);
884  if(out.get() && trace)
885  *out
886  << "\nComputing variations *= (1.0)/(dwgt*uh),"
887  << " where (1.0)/(dwgt*uh) = (1.0)/("<<dwgt<<"*"<<uh<<") = "<<alpha<<" ...\n";
888  VectorPtr f;
889  if( var.supports(MEB::OUT_ARG_f) && (f=var.get_f()).get() )
890  Vt_S(f.ptr(),alpha);
891  for( int j = 0; j < Ng; ++j ) {
892  VectorPtr g_j;
893  if( (g_j=var.get_g(j)).get() )
894  Vt_S(g_j.ptr(),alpha);
895  }
896  if(out.get() && trace)
897  *out << "\nFinal variations=\n" << describe(var,verbLevel);
898  }
899 
900  if(out.get())
901  *out << std::setprecision(p_saved);
902 
903  if(out.get() && trace)
904  *out << "\nLeaving DirectionalFiniteDiffCalculator<Scalar>::calcVariations(...)\n";
905 
906 }
907 
908 
909 template<class Scalar>
911  const ModelEvaluator<Scalar> &model,
912  const ModelEvaluatorBase::InArgs<Scalar> &bp, // basePoint
913  const ModelEvaluatorBase::OutArgs<Scalar> &bfunc, // baseFunctionValues
914  const ModelEvaluatorBase::OutArgs<Scalar> &deriv // derivatives
915  ) const
916 {
917 
918  using std::string;
919  //typedef Teuchos::ScalarTraits<Scalar> ST;
920 
921  THYRA_FUNC_TIME_MONITOR(
922  string("Thyra::DirectionalFiniteDiffCalculator<")+ST::name()+">::calcDerivatives(...)"
923  );
924 
925  typedef ModelEvaluatorBase MEB;
926  typedef RCP<VectorBase<Scalar> > VectorPtr;
927  typedef RCP<MultiVectorBase<Scalar> > MultiVectorPtr;
928 
929  RCP<Teuchos::FancyOStream> out = this->getOStream();
930  Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
931  const bool trace = (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM));
932  Teuchos::OSTab tab(out);
933 
934  if(out.get() && trace)
935  *out << "\nEntering DirectionalFiniteDiffCalculator<Scalar>::calcDerivatives(...)\n";
936 
937  if(out.get() && trace)
938  *out
939  << "\nbasePoint=\n" << describe(bp,verbLevel)
940  << "\nbaseFunctionValues=\n" << describe(bfunc,verbLevel)
941  << "\nderivatives=\n" << describe(deriv,Teuchos::VERB_LOW);
942 
943  //
944  // We will only compute finite differences w.r.t. p(l) for now
945  //
946  const int Np = bp.Np(), Ng = bfunc.Ng();
947  MEB::InArgs<Scalar> dir = model.createInArgs();
948  MEB::OutArgs<Scalar> var = model.createOutArgs();
949  MultiVectorPtr DfDp_l;
950  std::vector<MEB::DerivativeMultiVector<Scalar> > DgDp_l(Ng);
951  std::vector<VectorPtr> var_g(Ng);
952  for( int l = 0; l < Np; ++l ) {
953  if(out.get() && trace)
954  *out << "\nComputing derivatives for parameter subvector p("<<l<<") ...\n";
955  Teuchos::OSTab tab2(out);
956  //
957  // Load up OutArgs var object of derivative variations to compute
958  //
959  bool hasDerivObject = false;
960  // DfDp(l)
961  if (
962  !deriv.supports(MEB::OUT_ARG_DfDp,l).none()
963  && !deriv.get_DfDp(l).isEmpty()
964  )
965  {
966  hasDerivObject = true;
967  std::ostringstream name; name << "DfDp("<<l<<")";
968  DfDp_l = get_mv(deriv.get_DfDp(l),name.str(),MEB::DERIV_MV_BY_COL);
969  }
970  else {
971  DfDp_l = Teuchos::null;
972  }
973  // DgDp(j=1...Ng,l)
974  for ( int j = 0; j < Ng; ++j ) {
975  if (
976  !deriv.supports(MEB::OUT_ARG_DgDp,j,l).none()
977  &&
978  !deriv.get_DgDp(j,l).isEmpty()
979  )
980  {
981  hasDerivObject = true;
982  std::ostringstream name; name << "DgDp("<<j<<","<<l<<")";
983  DgDp_l[j] = get_dmv(deriv.get_DgDp(j,l),name.str());
984  if( DgDp_l[j].getMultiVector().get() && !var_g[j].get() )
985  {
986  // Need a temporary vector for the variation for g(j)
987  var_g[j] = createMember(model.get_g_space(j));
988  }
989  }
990  else{
991  DgDp_l[j] = MEB::DerivativeMultiVector<Scalar>();
992  var_g[j] = Teuchos::null;
993  }
994  }
995  //
996  // Compute derivative variations by finite differences
997  //
998  if (hasDerivObject) {
999  VectorPtr e_i = createMember(model.get_p_space(l));
1000  dir.set_p(l,e_i);
1001  assign(e_i.ptr(),ST::zero());
1002  const int np_l = e_i->space()->dim();
1003  for( int i = 0 ; i < np_l; ++ i ) {
1004  if(out.get() && trace)
1005  *out << "\nComputing derivatives for single variable p("<<l<<")("<<i<<") ...\n";
1006  Teuchos::OSTab tab3(out);
1007  if(DfDp_l.get()) var.set_f(DfDp_l->col(i)); // Compute DfDp(l)(i) in place!
1008  for(int j = 0; j < Ng; ++j) {
1009  MultiVectorPtr DgDp_j_l;
1010  if( (DgDp_j_l=DgDp_l[j].getMultiVector()).get() ) {
1011  var.set_g(j,var_g[j]); // Computes d(g(j))/d(p(l)(i))
1012  }
1013  }
1014  set_ele(i,ST::one(),e_i.ptr());
1015  this->calcVariations(
1016  model,bp,dir,bfunc,var
1017  );
1018  set_ele(i,ST::zero(),e_i.ptr());
1019  if (DfDp_l.get()) var.set_f(Teuchos::null);
1020  for (int j = 0; j < Ng; ++j) {
1021  MultiVectorPtr DgDp_j_l;
1022  if ( !is_null(DgDp_j_l=DgDp_l[j].getMultiVector()) ) {
1023  assign( DgDp_j_l->col(i).ptr(), *var_g[j] );
1024  }
1025  }
1026  }
1027  dir.set_p(l,Teuchos::null);
1028  }
1029  }
1030 
1031  if(out.get() && trace)
1032  *out
1033  << "\nderivatives=\n" << describe(deriv,verbLevel);
1034 
1035  if(out.get() && trace)
1036  *out << "\nLeaving DirectionalFiniteDiffCalculator<Scalar>::calcDerivatives(...)\n";
1037 
1038 }
1039 
1040 
1041 } // namespace Thyra
1042 
1043 
1044 #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()...
int Ng() const
Return the number of axillary response functions g(j)(...) supported (Ng &gt;= 0).
virtual RCP< const VectorSpaceBase< Scalar > > get_x_space() const =0
Return the vector space for the state variables x &lt;: 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 &quot;model&quot; 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 &paramList)
bool is_null(const boost::shared_ptr< T > &p)
DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType EFDStepSelectType
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.
RCP< const VectorBase< Scalar > > get_x_dot() const
Precondition: supports(IN_ARG_x_dot)==true.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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.
Evaluation< VectorBase< Scalar > > get_g(int j) const
Precondition: supports(OUT_ARG_g)==true..
Evaluation< VectorBase< Scalar > > get_f() const
Precondition: supports(OUT_ARG_f)==true.
T * get() const
virtual ModelEvaluatorBase::InArgs< Scalar > createInArgs() const =0
Create an empty input arguments object that can be set up and passed to evalModel().
bool supports(EOutArgsMembers arg) const
Determine if an input argument is supported or not.
Derivative< Scalar > get_DgDp(int j, int l) const
Precondition: supports(OUT_ARG_DgDp,j,l)==true.
IntegralType getIntegralValue(ParameterList const &paramList, std::string const &paramName)
DirectionalFiniteDiffCalculatorTypes::EFDMethodType EFDMethodType
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 &gt;= 0).
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual RCP< const VectorSpaceBase< Scalar > > get_p_space(int l) const =0
Return the vector space for the auxiliary parameters p(l) &lt;: RE^n_p_l.
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.
void validateParameters(ParameterList const &validParamList, int const depth=1000, EValidateUsed const validateUsed=VALIDATE_USED_ENABLED, EValidateDefaults const validateDefaults=VALIDATE_DEFAULTS_ENABLED) const
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(...) &lt;: 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)
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
TypeTo as(const TypeFrom &t)
virtual RCP< const VectorSpaceBase< Scalar > > get_g_space(int j) const =0
Return the vector space for the auxiliary response functions g(j) &lt;: RE^n_g_j.
int Np() const
Return the number of parameter subvectors p(l) supported (Np &gt;= 0).
RCP< const VectorBase< Scalar > > get_x_dot_dot() const
Precondition: supports(IN_ARG_x_dot_dot)==true.
bool supports(EInArgsMembers arg) const
Determines if an input argument is supported or not.
RCP< const VectorBase< Scalar > > get_x() const
Precondition: supports(IN_ARG_x)==true.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
RCP< const VectorBase< Scalar > > get_p(int l) const
Get p(l) where 0 &lt;= l &amp;&amp; l &lt; this-&gt;Np().
Concrete aggregate class for all input arguments computable by a ModelEvaluator subclass object...