Thyra  Version of the Day
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); 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>(); }
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();
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!"
162  );
163 #endif
164  }
165 
166 };
167 
168 
169 } // namespace DirectionalFiniteDiffCalculatorTypes
170 
171 
172 // Private static data members
173 
174 
175 template<class Scalar>
176 const std::string& DirectionalFiniteDiffCalculator<Scalar>::FDMethod_name()
177 {
178  static std::string loc_FDMethod_name = "FD Method";
179  return loc_FDMethod_name;
180 }
181 
182 
183 template<class Scalar>
184 const RCP<
185  Teuchos::StringToIntegralParameterEntryValidator<
186  Thyra::DirectionalFiniteDiffCalculatorTypes::EFDMethodType
187  >
188 >&
189 DirectionalFiniteDiffCalculator<Scalar>::fdMethodValidator()
190 {
191  static RCP<Teuchos::StringToIntegralParameterEntryValidator<EFDMethodType> >
192  loc_fdMethodValidator
193  = Teuchos::rcp(
194  new Teuchos::StringToIntegralParameterEntryValidator<Thyra::DirectionalFiniteDiffCalculatorTypes::EFDMethodType>(
195  Teuchos::tuple<std::string>(
196  "order-one"
197  ,"order-two"
198  ,"order-two-central"
199  ,"order-two-auto"
200  ,"order-four"
201  ,"order-four-central"
202  ,"order-four-auto"
203  )
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\""
212  )
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
221  )
222  ,""
223  )
224  );
225  return loc_fdMethodValidator;
226 }
227 
228 
229 template<class Scalar>
230 const std::string&
231 DirectionalFiniteDiffCalculator<Scalar>::FDMethod_default()
232 {
233  static std::string loc_FDMethod_default = "order-one";
234  return loc_FDMethod_default;
235 }
236 
237 
238 template<class Scalar>
239 const std::string&
240 DirectionalFiniteDiffCalculator<Scalar>::FDStepSelectType_name()
241 {
242  static std::string loc_FDStepSelectType_name = "FD Step Select Type";
243  return loc_FDStepSelectType_name;
244 }
245 
246 
247 template<class Scalar>
248 const RCP<
249  Teuchos::StringToIntegralParameterEntryValidator<
250  Thyra::DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType
251  >
252  >&
253 DirectionalFiniteDiffCalculator<Scalar>::fdStepSelectTypeValidator()
254 {
255  static const RCP<Teuchos::StringToIntegralParameterEntryValidator<EFDStepSelectType> >
256  loc_fdStepSelectTypeValidator
257  = Teuchos::rcp(
258  new Teuchos::StringToIntegralParameterEntryValidator<Thyra::DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType>(
259  Teuchos::tuple<std::string>(
260  "Absolute"
261  ,"Relative"
262  )
263  ,Teuchos::tuple<std::string>(
264  "Use absolute step size \""+FDStepLength_name()+"\""
265  ,"Use relative step size \""+FDStepLength_name()+"\"*||xo||inf"
266  )
267  ,Teuchos::tuple<Thyra::DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType>(
268  Thyra::DirectionalFiniteDiffCalculatorTypes::FD_STEP_ABSOLUTE
269  ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_STEP_RELATIVE
270  )
271  ,""
272  )
273  );
274  return loc_fdStepSelectTypeValidator;
275 }
276 
277 
278 template<class Scalar>
279 const std::string&
280 DirectionalFiniteDiffCalculator<Scalar>::FDStepSelectType_default()
281 {
282  static std::string loc_FDStepSelectType_default = "Absolute";
283  return loc_FDStepSelectType_default;
284 }
285 
286 
287 template<class Scalar>
288 const std::string&
289 DirectionalFiniteDiffCalculator<Scalar>::FDStepLength_name()
290 {
291  static std::string loc_FDStepLength_name = "FD Step Length";
292  return loc_FDStepLength_name;
293 }
294 
295 
296 template<class Scalar>
297 const double&
298 DirectionalFiniteDiffCalculator<Scalar>::FDStepLength_default()
299 {
300  static double loc_FDStepLength_default = -1.0;
301  return loc_FDStepLength_default;
302 }
303 
304 
305 // Constructors/initializer
306 
307 
308 template<class Scalar>
310  EFDMethodType fd_method_type_in
311  ,EFDStepSelectType fd_step_select_type_in
312  ,ScalarMag fd_step_size_in
313  ,ScalarMag fd_step_size_min_in
314  )
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)
319 {}
320 
321 
322 // Overriden from ParameterListAcceptor
323 
324 
325 template<class Scalar>
327  RCP<ParameterList> const& paramList
328  )
329 {
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);
340 }
341 
342 
343 template<class Scalar>
344 RCP<ParameterList>
346 {
347  return paramList_;
348 }
349 
350 
351 template<class Scalar>
352 RCP<ParameterList>
354 {
355  RCP<ParameterList> _paramList = paramList_;
356  paramList_ = Teuchos::null;
357  return _paramList;
358 }
359 
360 
361 template<class Scalar>
362 RCP<const ParameterList>
364 {
365  return paramList_;
366 }
367 
368 
369 template<class Scalar>
370 RCP<const ParameterList>
372 {
373  using Teuchos::rcp_implicit_cast;
374  typedef Teuchos::ParameterEntryValidator PEV;
375  static RCP<ParameterList> pl;
376  if(pl.get()==NULL) {
377  pl = Teuchos::parameterList();
378  pl->set(
379  FDMethod_name(), FDMethod_default(),
380  "The method used to compute the finite differences.",
381  rcp_implicit_cast<const PEV>(fdMethodValidator())
382  );
383  pl->set(
384  FDStepSelectType_name(), FDStepSelectType_default(),
385  "Method used to select the finite difference step length.",
386  rcp_implicit_cast<const PEV>(fdStepSelectTypeValidator())
387  );
388  pl->set(
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."
392  );
393  Teuchos::setupVerboseObjectSublist(&*pl);
394  }
395  return pl;
396 }
397 
398 
399 template<class Scalar>
402  const ModelEvaluator<Scalar> &model,
403  const SelectedDerivatives &fdDerivatives
404  )
405 {
406  return DirectionalFiniteDiffCalculatorTypes::OutArgsCreator<Scalar>::createOutArgs(
407  model, fdDerivatives );
408 }
409 
410 
411 template<class Scalar>
413  const ModelEvaluator<Scalar> &model,
414  const ModelEvaluatorBase::InArgs<Scalar> &bp, // basePoint
415  const ModelEvaluatorBase::InArgs<Scalar> &dir, // directions
416  const ModelEvaluatorBase::OutArgs<Scalar> &bfunc, // baseFunctionValues
417  const ModelEvaluatorBase::OutArgs<Scalar> &var // variations
418  ) const
419 {
420 
421  using std::string;
422 
423  THYRA_FUNC_TIME_MONITOR(
424  string("Thyra::DirectionalFiniteDiffCalculator<")+ST::name()+">::calcVariations(...)"
425  );
426 
427  using std::setw;
428  using std::endl;
429  using std::right;
430  using Teuchos::as;
431  typedef ModelEvaluatorBase MEB;
432  namespace DFDCT = DirectionalFiniteDiffCalculatorTypes;
433  typedef RCP<VectorBase<Scalar> > VectorPtr;
434 
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);
439 
440  if(out.get() && trace)
441  *out << "\nEntering DirectionalFiniteDiffCalculator<Scalar>::calcVariations(...)\n";
442 
443  if(out.get() && trace)
444  *out
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);
449 
450 #ifdef TEUCHOS_DEBUG
451  TEUCHOS_TEST_FOR_EXCEPTION(
452  var.isEmpty(), std::logic_error,
453  "Error, all of the variations can not be null!"
454  );
455 #endif
456 
457  //
458  // To illustrate the theory behind this implementation consider
459  // the generic multi-variable function h(z) <: R^n -> R. Now let's
460  // consider we have the base point zo and the vector v to
461  // perturb h(z) along. First form the function h(zo+a*v) and then
462  // let's compute d(h)/d(a) at a = 0:
463  //
464  // (1) d(h(zo+a*v))/d(a) at a = 0
465  // = sum( d(h)/d(x(i)) * d(x(i))/d(a), i = 1...n)
466  // = sum( d(h)/d(x(i)) * v(i), i = 1...n)
467  // = d(h)/d(a) * v
468  //
469  // Now we can approximate (1) using, for example, central differences as:
470  //
471  // (2) d(h(zo+a*v))/d(a) at a = 0
472  // \approx ( h(zo+h*v) - h(zo+h*v) ) / (2*h)
473  //
474  // If we equate (1) and (2) we have the approximation:
475  //
476  // (3) d(h)/d(a) * v \approx ( h(zo+h*v) - g(zo+h*v) ) / (2*h)
477  //
478  //
479 
480  // /////////////////////////////////////////
481  // Validate the input
482 
483  // ToDo: Validate input!
484 
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";
488  break;
489  case DFDCT::FD_ORDER_TWO:
490  if(out.get()&&trace) *out<<"\nUsing one-sided, second-order finite differences ...\n";
491  break;
492  case DFDCT::FD_ORDER_TWO_CENTRAL:
493  if(out.get()&&trace) *out<<"\nUsing second-order central finite differences ...\n";
494  break;
495  case DFDCT::FD_ORDER_TWO_AUTO:
496  if(out.get()&&trace) *out<<"\nUsing auto selection of some second-order finite difference method ...\n";
497  break;
498  case DFDCT::FD_ORDER_FOUR:
499  if(out.get()&&trace) *out<<"\nUsing one-sided, fourth-order finite differences ...\n";
500  break;
501  case DFDCT::FD_ORDER_FOUR_CENTRAL:
502  if(out.get()&&trace) *out<<"\nUsing fourth-order central finite differences ...\n";
503  break;
504  case DFDCT::FD_ORDER_FOUR_AUTO:
505  if(out.get()&&trace) *out<<"\nUsing auto selection of some fourth-order finite difference method ...\n";
506  break;
507  default:
508  TEUCHOS_TEST_FOR_EXCEPT(true); // Should not get here!
509  }
510 
511  // ////////////////////////
512  // Find the step size
513 
514  //
515  // Get defaults for the optimal step sizes
516  //
517 
518  const ScalarMag
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);
523 
524  ScalarMag
525  bp_norm = SMT::zero();
526  // ToDo above: compute a reasonable norm somehow based on the base-point vector(s)!
527 
528  ScalarMag
529  uh_opt = 0.0;
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 );
533  break;
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 );
538  break;
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 );
543  break;
544  default:
545  TEUCHOS_TEST_FOR_EXCEPT(true); // Should not get here!
546  }
547 
548  if(out.get()&&trace) *out<<"\nDefault optimal step length uh_opt = " << uh_opt << " ...\n";
549 
550  //
551  // Set the step sizes used.
552  //
553 
554  ScalarMag
555  uh = this->fd_step_size();
556 
557  if( uh < 0 )
558  uh = uh_opt;
559  else if( fd_step_select_type() == DFDCT::FD_STEP_RELATIVE )
560  uh *= (bp_norm + 1.0);
561 
562  if(out.get()&&trace) *out<<"\nStep size to be used uh="<<uh<<"\n";
563 
564  //
565  // Find step lengths that stay in bounds!
566  //
567 
568  // ToDo: Add logic for variable bounds when needed!
569 
570  //
571  // Set the actual method being used
572  //
573  // ToDo: Consider cramped bounds and method order.
574  //
575 
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;
580  break;
581  case DFDCT::FD_ORDER_FOUR_AUTO:
582  l_fd_method_type = DFDCT::FD_ORDER_FOUR_CENTRAL;
583  break;
584  default:
585  break; // Okay
586  }
587 
588  //if(out.get()&&trace) *out<<"\nStep size to fit in bounds: uh="<<uh"\n";
589 
590  int p_saved = -1;
591  if(out.get())
592  p_saved = out->precision();
593 
594  // ///////////////////////////////////////////////
595  // Compute the directional derivatives
596 
597  const int Np = var.Np(), Ng = var.Ng();
598 
599  // Setup storage for perturbed variables
600  VectorPtr per_x, per_x_dot;
601  std::vector<VectorPtr> per_p(Np);
602  MEB::InArgs<Scalar> pp = model.createInArgs();
603  pp.setArgs(bp); // Set all args to start with
604  if( bp.supports(MEB::IN_ARG_x) ) {
605  if( dir.get_x().get() )
606  pp.set_x(per_x=createMember(model.get_x_space()));
607  else
608  pp.set_x(bp.get_x());
609  }
610  if( bp.supports(MEB::IN_ARG_x_dot) ) {
611  if( dir.get_x_dot().get() )
612  pp.set_x_dot(per_x_dot=createMember(model.get_x_space()));
613  else
614  pp.set_x_dot(bp.get_x_dot());
615  }
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)));
619  else
620  pp.set_p(l,bp.get_p(l));
621  }
622  if(out.get() && trace)
623  *out
624  << "\nperturbedPoint after initial setup (with some unintialized vectors) =\n"
625  << describe(pp,verbLevel);
626 
627  // Setup storage for perturbed functions
628  bool all_funcs_at_base_computed = true;
629  MEB::OutArgs<Scalar> pfunc = model.createOutArgs();
630  {
631  VectorPtr f;
632  if( var.supports(MEB::OUT_ARG_f) && (f=var.get_f()).get() ) {
633  pfunc.set_f(createMember(model.get_f_space()));
634  assign(f.ptr(),ST::zero());
635  if(!bfunc.get_f().get()) all_funcs_at_base_computed = false;
636  }
637  for( int j = 0; j < Ng; ++j ) {
638  VectorPtr g_j;
639  if( (g_j=var.get_g(j)).get() ) {
640  pfunc.set_g(j,createMember(model.get_g_space(j)));
641  assign(g_j.ptr(),ST::zero());
642  if(!bfunc.get_g(j).get()) all_funcs_at_base_computed = false;
643  }
644  }
645  }
646  if(out.get() && trace)
647  *out
648  << "\nperturbedFunctions after initial setup (with some unintialized vectors) =\n"
649  << describe(pfunc,verbLevel);
650 
651  const int dbl_p = 15;
652  if(out.get())
653  *out << std::setprecision(dbl_p);
654 
655  //
656  // Compute the weighted sum of the terms
657  //
658 
659  int num_evals = 0;
660  ScalarMag dwgt = SMT::zero();
661  switch(l_fd_method_type) {
662  case DFDCT::FD_ORDER_ONE: // may only need one eval if f(xo) etc is passed in
663  num_evals = 2;
664  dwgt = ScalarMag(1.0);
665  break;
666  case DFDCT::FD_ORDER_TWO: // may only need two evals if c(xo) etc is passed in
667  num_evals = 3;
668  dwgt = ScalarMag(2.0);
669  break;
670  case DFDCT::FD_ORDER_TWO_CENTRAL:
671  num_evals = 2;
672  dwgt = ScalarMag(2.0);
673  break;
674  case DFDCT::FD_ORDER_FOUR:
675  num_evals = 5;
676  dwgt = ScalarMag(12.0);
677  break;
678  case DFDCT::FD_ORDER_FOUR_CENTRAL:
679  num_evals = 4;
680  dwgt = ScalarMag(12.0);
681  break;
682  default:
683  TEUCHOS_TEST_FOR_EXCEPT(true); // Should not get here!
684  }
685  for( int eval_i = 1; eval_i <= num_evals; ++eval_i ) {
686  // Set the step constant and the weighting constant
687  ScalarMag
688  uh_i = 0.0,
689  wgt_i = 0.0;
690  switch(l_fd_method_type) {
691  case DFDCT::FD_ORDER_ONE: {
692  switch(eval_i) {
693  case 1:
694  uh_i = +0.0;
695  wgt_i = -1.0;
696  break;
697  case 2:
698  uh_i = +1.0;
699  wgt_i = +1.0;
700  break;
701  }
702  break;
703  }
704  case DFDCT::FD_ORDER_TWO: {
705  switch(eval_i) {
706  case 1:
707  uh_i = +0.0;
708  wgt_i = -3.0;
709  break;
710  case 2:
711  uh_i = +1.0;
712  wgt_i = +4.0;
713  break;
714  case 3:
715  uh_i = +2.0;
716  wgt_i = -1.0;
717  break;
718  }
719  break;
720  }
721  case DFDCT::FD_ORDER_TWO_CENTRAL: {
722  switch(eval_i) {
723  case 1:
724  uh_i = -1.0;
725  wgt_i = -1.0;
726  break;
727  case 2:
728  uh_i = +1.0;
729  wgt_i = +1.0;
730  break;
731  }
732  break;
733  }
734  case DFDCT::FD_ORDER_FOUR: {
735  switch(eval_i) {
736  case 1:
737  uh_i = +0.0;
738  wgt_i = -25.0;
739  break;
740  case 2:
741  uh_i = +1.0;
742  wgt_i = +48.0;
743  break;
744  case 3:
745  uh_i = +2.0;
746  wgt_i = -36.0;
747  break;
748  case 4:
749  uh_i = +3.0;
750  wgt_i = +16.0;
751  break;
752  case 5:
753  uh_i = +4.0;
754  wgt_i = -3.0;
755  break;
756  }
757  break;
758  }
759  case DFDCT::FD_ORDER_FOUR_CENTRAL: {
760  switch(eval_i) {
761  case 1:
762  uh_i = -2.0;
763  wgt_i = +1.0;
764  break;
765  case 2:
766  uh_i = -1.0;
767  wgt_i = -8.0;
768  break;
769  case 3:
770  uh_i = +1.0;
771  wgt_i = +8.0;
772  break;
773  case 4:
774  uh_i = +2.0;
775  wgt_i = -1.0;
776  break;
777  }
778  break;
779  }
780  case DFDCT::FD_ORDER_TWO_AUTO:
781  case DFDCT::FD_ORDER_FOUR_AUTO:
782  break; // Okay
783  default:
784  TEUCHOS_TEST_FOR_EXCEPT(true);
785  }
786 
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);
790 
791  // Compute the weighted term and add it to the sum
792  if(uh_i == ST::zero()) {
793  MEB::OutArgs<Scalar> bfuncall;
794  if(!all_funcs_at_base_computed) {
795  // Compute missing functions at the base point
796  bfuncall = model.createOutArgs();
797  if( pfunc.supports(MEB::OUT_ARG_f) && pfunc.get_f().get() && !bfunc.get_f().get() ) {
798  bfuncall.set_f(createMember(model.get_f_space()));
799  }
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)));
803  }
804  }
805  model.evalModel(bp,bfuncall);
806  bfuncall.setArgs(bfunc);
807  }
808  else {
809  bfuncall = bfunc;
810  }
811  // Use functions at the base point
812  if(out.get() && trace)
813  *out << "\nSetting variations = wgt_i * basePoint ...\n";
814  VectorPtr f;
815  if( pfunc.supports(MEB::OUT_ARG_f) && (f=var.get_f()).get() ) {
816  V_StV<Scalar>(f.ptr(), wgt_i, *bfuncall.get_f());
817  }
818  for( int j = 0; j < Ng; ++j ) {
819  VectorPtr g_j;
820  if( (g_j=var.get_g(j)).get() ) {
821  V_StV<Scalar>(g_j.ptr(), wgt_i, *bfuncall.get_g(j));
822  }
823  }
824  }
825  else {
826  if(out.get() && trace)
827  *out << "\nSetting perturbedPoint = basePoint + uh_i*uh*direction ...\n";
828  // z = zo + uh_i*uh*v
829  {
830  if ( dir.supports(MEB::IN_ARG_x) && dir.get_x().get() )
831  V_StVpV(per_x.ptr(),as<Scalar>(uh_i*uh),*dir.get_x(),*bp.get_x());
832  if ( dir.supports(MEB::IN_ARG_x_dot) && dir.get_x_dot().get() )
833  V_StVpV(per_x_dot.ptr(), as<Scalar>(uh_i*uh), *dir.get_x_dot(), *bp.get_x_dot());
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));
837  }
838  }
839  if(out.get() && trace)
840  *out << "\nperturbedPoint=\n" << describe(pp,verbLevel);
841  // Compute perturbed function values h(zo+uh_i*uh)
842  if(out.get() && trace)
843  *out << "\nCompute perturbedFunctions at perturbedPoint...\n";
844  model.evalModel(pp,pfunc);
845  if(out.get() && trace)
846  *out << "\nperturbedFunctions=\n" << describe(pfunc,verbLevel);
847  // Sum perturbed function values into total variation
848  {
849  // var_h += wgt_i*perturbed_h
850  if(out.get() && trace)
851  *out << "\nComputing variations += wgt_i*perturbedfunctions ...\n";
852  VectorPtr f;
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 ) {
856  VectorPtr g_j;
857  if( (g_j=pfunc.get_g(j)).get() )
858  Vp_StV<Scalar>(var.get_g(j).ptr(), wgt_i, *g_j);
859  }
860  }
861  }
862  if(out.get() && trace)
863  *out << "\nvariations=\n" << describe(var,verbLevel);
864  }
865 
866  //
867  // Multiply by the scaling factor!
868  //
869 
870  {
871  // var_h *= 1.0/(dwgt*uh)
872  const Scalar alpha = ST::one()/(dwgt*uh);
873  if(out.get() && trace)
874  *out
875  << "\nComputing variations *= (1.0)/(dwgt*uh),"
876  << " where (1.0)/(dwgt*uh) = (1.0)/("<<dwgt<<"*"<<uh<<") = "<<alpha<<" ...\n";
877  VectorPtr f;
878  if( var.supports(MEB::OUT_ARG_f) && (f=var.get_f()).get() )
879  Vt_S(f.ptr(),alpha);
880  for( int j = 0; j < Ng; ++j ) {
881  VectorPtr g_j;
882  if( (g_j=var.get_g(j)).get() )
883  Vt_S(g_j.ptr(),alpha);
884  }
885  if(out.get() && trace)
886  *out << "\nFinal variations=\n" << describe(var,verbLevel);
887  }
888 
889  if(out.get())
890  *out << std::setprecision(p_saved);
891 
892  if(out.get() && trace)
893  *out << "\nLeaving DirectionalFiniteDiffCalculator<Scalar>::calcVariations(...)\n";
894 
895 }
896 
897 
898 template<class Scalar>
900  const ModelEvaluator<Scalar> &model,
901  const ModelEvaluatorBase::InArgs<Scalar> &bp, // basePoint
902  const ModelEvaluatorBase::OutArgs<Scalar> &bfunc, // baseFunctionValues
903  const ModelEvaluatorBase::OutArgs<Scalar> &deriv // derivatives
904  ) const
905 {
906 
907  using std::string;
908  //typedef Teuchos::ScalarTraits<Scalar> ST;
909 
910  THYRA_FUNC_TIME_MONITOR(
911  string("Thyra::DirectionalFiniteDiffCalculator<")+ST::name()+">::calcDerivatives(...)"
912  );
913 
914  typedef ModelEvaluatorBase MEB;
915  typedef RCP<VectorBase<Scalar> > VectorPtr;
916  typedef RCP<MultiVectorBase<Scalar> > MultiVectorPtr;
917 
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);
922 
923  if(out.get() && trace)
924  *out << "\nEntering DirectionalFiniteDiffCalculator<Scalar>::calcDerivatives(...)\n";
925 
926  if(out.get() && trace)
927  *out
928  << "\nbasePoint=\n" << describe(bp,verbLevel)
929  << "\nbaseFunctionValues=\n" << describe(bfunc,verbLevel)
930  << "\nderivatives=\n" << describe(deriv,Teuchos::VERB_LOW);
931 
932  //
933  // We will only compute finite differences w.r.t. p(l) for now
934  //
935  const int Np = bp.Np(), Ng = bfunc.Ng();
936  MEB::InArgs<Scalar> dir = model.createInArgs();
937  MEB::OutArgs<Scalar> var = model.createOutArgs();
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);
945  //
946  // Load up OutArgs var object of derivative variations to compute
947  //
948  bool hasDerivObject = false;
949  // DfDp(l)
950  if (
951  !deriv.supports(MEB::OUT_ARG_DfDp,l).none()
952  && !deriv.get_DfDp(l).isEmpty()
953  )
954  {
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);
958  }
959  else {
960  DfDp_l = Teuchos::null;
961  }
962  // DgDp(j=1...Ng,l)
963  for ( int j = 0; j < Ng; ++j ) {
964  if (
965  !deriv.supports(MEB::OUT_ARG_DgDp,j,l).none()
966  &&
967  !deriv.get_DgDp(j,l).isEmpty()
968  )
969  {
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() )
974  {
975  // Need a temporary vector for the variation for g(j)
976  var_g[j] = createMember(model.get_g_space(j));
977  }
978  }
979  else{
980  DgDp_l[j] = MEB::DerivativeMultiVector<Scalar>();
981  var_g[j] = Teuchos::null;
982  }
983  }
984  //
985  // Compute derivative variations by finite differences
986  //
987  if (hasDerivObject) {
988  VectorPtr e_i = createMember(model.get_p_space(l));
989  dir.set_p(l,e_i);
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)); // Compute DfDp(l)(i) in place!
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]); // Computes d(g(j))/d(p(l)(i))
1001  }
1002  }
1003  set_ele(i,ST::one(),e_i.ptr());
1004  this->calcVariations(
1005  model,bp,dir,bfunc,var
1006  );
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] );
1013  }
1014  }
1015  }
1016  dir.set_p(l,Teuchos::null);
1017  }
1018  }
1019 
1020  if(out.get() && trace)
1021  *out
1022  << "\nderivatives=\n" << describe(deriv,verbLevel);
1023 
1024  if(out.get() && trace)
1025  *out << "\nLeaving DirectionalFiniteDiffCalculator<Scalar>::calcDerivatives(...)\n";
1026 
1027 }
1028 
1029 
1030 } // namespace Thyra
1031 
1032 
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 &paramList)
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).
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.
Evaluation< VectorBase< Scalar > > get_f() const
Precondition: supports(OUT_ARG_f)==true.
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...