Thyra Version of the Day
Loading...
Searching...
No Matches
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"
54#include "Teuchos_VerboseObjectParameterListHelpers.hpp"
55
56
57namespace Thyra {
58
59
60namespace 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//
73template<class Scalar>
74class OutArgsCreator : public StateFuncModelEvaluatorBase<Scalar>
75{
76public:
77 // Public functions overridden from ModelEvaulator.
78 RCP<const VectorSpaceBase<Scalar> > get_x_space() const
80 RCP<const VectorSpaceBase<Scalar> > get_f_space() const
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
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
151private:
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
178template<class Scalar>
179const std::string& DirectionalFiniteDiffCalculator<Scalar>::FDMethod_name()
180{
181 static std::string loc_FDMethod_name = "FD Method";
182 return loc_FDMethod_name;
183}
184
185
186template<class Scalar>
187const RCP<
189 Thyra::DirectionalFiniteDiffCalculatorTypes::EFDMethodType
190 >
191>&
192DirectionalFiniteDiffCalculator<Scalar>::fdMethodValidator()
193{
194 static RCP<Teuchos::StringToIntegralParameterEntryValidator<EFDMethodType> >
195 loc_fdMethodValidator
196 = Teuchos::rcp(
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 )
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 )
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
232template<class Scalar>
233const std::string&
234DirectionalFiniteDiffCalculator<Scalar>::FDMethod_default()
235{
236 static std::string loc_FDMethod_default = "order-one";
237 return loc_FDMethod_default;
238}
239
240
241template<class Scalar>
242const std::string&
243DirectionalFiniteDiffCalculator<Scalar>::FDStepSelectType_name()
244{
245 static std::string loc_FDStepSelectType_name = "FD Step Select Type";
246 return loc_FDStepSelectType_name;
247}
248
249
250template<class Scalar>
251const RCP<
253 Thyra::DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType
254 >
255 >&
256DirectionalFiniteDiffCalculator<Scalar>::fdStepSelectTypeValidator()
257{
258 static const RCP<Teuchos::StringToIntegralParameterEntryValidator<EFDStepSelectType> >
259 loc_fdStepSelectTypeValidator
260 = Teuchos::rcp(
263 "Absolute"
264 ,"Relative"
265 )
267 "Use absolute step size \""+FDStepLength_name()+"\""
268 ,"Use relative step size \""+FDStepLength_name()+"\"*||xo||inf"
269 )
271 Thyra::DirectionalFiniteDiffCalculatorTypes::FD_STEP_ABSOLUTE
272 ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_STEP_RELATIVE
273 )
274 ,""
275 )
276 );
277 return loc_fdStepSelectTypeValidator;
278}
279
280
281template<class Scalar>
282const std::string&
283DirectionalFiniteDiffCalculator<Scalar>::FDStepSelectType_default()
284{
285 static std::string loc_FDStepSelectType_default = "Absolute";
286 return loc_FDStepSelectType_default;
287}
288
289
290template<class Scalar>
291const std::string&
292DirectionalFiniteDiffCalculator<Scalar>::FDStepLength_name()
293{
294 static std::string loc_FDStepLength_name = "FD Step Length";
295 return loc_FDStepLength_name;
296}
297
298
299template<class Scalar>
300const double&
301DirectionalFiniteDiffCalculator<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
311template<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
328template<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
346template<class Scalar>
352
353
354template<class Scalar>
357{
358 RCP<ParameterList> _paramList = paramList_;
359 paramList_ = Teuchos::null;
360 return _paramList;
361}
362
363
364template<class Scalar>
367{
368 return paramList_;
369}
370
371
372template<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
402template<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
414template<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
528 bp_norm = SMT::zero();
529 // ToDo above: compute a reasonable norm somehow based on the base-point vector(s)!
530
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
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
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
909template<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
T * get() const
Simple utility class used to select finite difference derivatives for OutArgs object.
void setParameterList(RCP< ParameterList > const &paramList)
ModelEvaluatorBase::OutArgs< Scalar > createOutArgs(const ModelEvaluator< Scalar > &model, const SelectedDerivatives &fdDerivatives)
Create an augmented out args object for holding finite difference objects.
DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType EFDStepSelectType
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.
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..
DirectionalFiniteDiffCalculatorTypes::EFDMethodType EFDMethodType
Concrete aggregate class for all input arguments computable by a ModelEvaluator subclass object.
RCP< const VectorBase< Scalar > > get_p(int l) const
Get p(l) where 0 <= l && l < this->Np().
RCP< const VectorBase< Scalar > > get_x() const
Precondition: supports(IN_ARG_x)==true.
RCP< const VectorBase< Scalar > > get_x_dot() const
Precondition: supports(IN_ARG_x_dot)==true.
bool supports(EInArgsMembers arg) const
Determines if an input argument is supported or not.
int Np() const
Return the number of parameter subvectors p(l) supported (Np >= 0).
RCP< const VectorBase< Scalar > > get_x_dot_dot() const
Precondition: supports(IN_ARG_x_dot_dot)==true.
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.
int Ng() const
Return the number of axillary response functions g(j)(...) supported (Ng >= 0).
Derivative< Scalar > get_DgDp(int j, int l) const
Precondition: supports(OUT_ARG_DgDp,j,l)==true.
Evaluation< VectorBase< Scalar > > get_f() const
Precondition: supports(OUT_ARG_f)==true.
bool supports(EOutArgsMembers arg) const
Determine if an input argument is supported or not.
int Np() const
Return the number of parameter subvectors p(l) supported (Np >= 0).
Evaluation< VectorBase< Scalar > > get_g(int j) const
Precondition: supports(OUT_ARG_g)==true..
Base subclass for ModelEvaluator that defines some basic types.
Pure abstract base interface for evaluating a stateless "model" that can be mapped into a number of d...
virtual RCP< const VectorSpaceBase< Scalar > > get_f_space() const =0
Return the vector space for the state function f(...) <: RE^n_x.
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_g_space(int j) const =0
Return the vector space for the auxiliary response functions g(j) <: RE^n_g_j.
virtual ModelEvaluatorBase::InArgs< Scalar > createInArgs() const =0
Create an empty input arguments object that can be set up and passed to evalModel().
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.
virtual RCP< const VectorSpaceBase< Scalar > > get_x_space() const =0
Return the vector space for the state variables x <: RE^n_x.
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.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TypeTo as(const TypeFrom &t)
#define TEUCHOS_UNREACHABLE_RETURN(dummyReturnVal)
T_To & dyn_cast(T_From &from)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)