Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_DirectionalFiniteDiffCalculator_def.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef THYRA_DIRECTIONAL_FINITE_DIFF_CALCULATOR_DEF_HPP
00030 #define THYRA_DIRECTIONAL_FINITE_DIFF_CALCULATOR_DEF_HPP
00031 
00032 
00033 #include "Thyra_DirectionalFiniteDiffCalculator_decl.hpp"
00034 #include "Thyra_ModelEvaluatorHelpers.hpp"
00035 #include "Thyra_DetachedVectorView.hpp"
00036 #include "Thyra_DetachedMultiVectorView.hpp"
00037 #include "Thyra_StateFuncModelEvaluatorBase.hpp"
00038 #include "Thyra_MultiVectorStdOps.hpp"
00039 #include "Thyra_VectorStdOps.hpp"
00040 #include "Teuchos_TimeMonitor.hpp"
00041 #include "Teuchos_VerboseObjectParameterListHelpers.hpp"
00042 
00043 
00044 namespace Thyra {
00045 
00046 
00047 namespace DirectionalFiniteDiffCalculatorTypes {
00048 
00049 
00050 //
00051 // Undocumented utility class for setting support for new derivative objects
00052 // on an OutArgs object!  Warning, users should not attempt to play these
00053 // tricks on their own!
00054 //
00055 // Note that because of the design of the OutArgs and OutArgsSetup classes,
00056 // you can only change the list of supported arguments in a subclass of
00057 // ModelEvaluatorBase since OutArgsSetup is a protected type.  The fact that
00058 // the only way to do this is convoluted is a feature!
00059 //
00060 template<class Scalar>
00061 class OutArgsCreator : public StateFuncModelEvaluatorBase<Scalar>
00062 {
00063 public:
00064   // Public functions overridden from ModelEvaulator.
00065   RCP<const VectorSpaceBase<Scalar> > get_x_space() const
00066     { TEST_FOR_EXCEPT(true); return Teuchos::null; }
00067   RCP<const VectorSpaceBase<Scalar> > get_f_space() const
00068     { TEST_FOR_EXCEPT(true); return Teuchos::null; }
00069   ModelEvaluatorBase::InArgs<Scalar> createInArgs() const
00070     { TEST_FOR_EXCEPT(true); return ModelEvaluatorBase::InArgs<Scalar>(); }
00071   ModelEvaluatorBase::OutArgs<Scalar> createOutArgs() const
00072     { TEST_FOR_EXCEPT(true); return ModelEvaluatorBase::OutArgs<Scalar>(); }
00073   void evalModel(
00074     const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00075     const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00076     ) const
00077     { TEST_FOR_EXCEPT(true); }
00078   // Static function that does the magic!
00079   static ModelEvaluatorBase::OutArgs<Scalar> createOutArgs(
00080     const ModelEvaluator<Scalar> &model,
00081     const SelectedDerivatives &fdDerivatives
00082     )
00083     {
00084       
00085       typedef ModelEvaluatorBase MEB;
00086 
00087       const MEB::OutArgs<Scalar> wrappedOutArgs = model.createOutArgs();
00088       const int Np = wrappedOutArgs.Np(), Ng = wrappedOutArgs.Ng();
00089       MEB::OutArgsSetup<Scalar> outArgs;
00090 
00091       outArgs.setModelEvalDescription(
00092         "DirectionalFiniteDiffCalculator: " + model.description()
00093         );
00094 
00095       // Start off by supporting everything that the underlying model supports
00096       // computing!
00097 
00098       outArgs.set_Np_Ng(Np,Ng);
00099       outArgs.setSupports(wrappedOutArgs);
00100 
00101       // Add support for finite difference DfDp(l) if asked
00102 
00103       const SelectedDerivatives::supports_DfDp_t
00104         &supports_DfDp = fdDerivatives.supports_DfDp_;
00105       for(
00106         SelectedDerivatives::supports_DfDp_t::const_iterator
00107           itr = supports_DfDp.begin();
00108         itr != supports_DfDp.end();
00109         ++itr
00110         )
00111       {
00112         const int l = *itr;
00113         assert_p_space(model,l);
00114         outArgs.setSupports(MEB::OUT_ARG_DfDp,l,MEB::DERIV_MV_BY_COL);
00115       }
00116 
00117       // Add support for finite difference DgDp(j,l) if asked
00118       
00119       const SelectedDerivatives::supports_DgDp_t
00120         &supports_DgDp = fdDerivatives.supports_DgDp_;
00121       for(
00122         SelectedDerivatives::supports_DgDp_t::const_iterator
00123           itr = supports_DgDp.begin();
00124         itr != supports_DgDp.end();
00125         ++itr
00126         )
00127       {
00128         const int j = itr->first;
00129         const int l = itr->second;
00130         assert_p_space(model,l);
00131         outArgs.setSupports(MEB::OUT_ARG_DgDp,j,l,MEB::DERIV_MV_BY_COL);
00132       }
00133 
00134       return outArgs;
00135 
00136     }
00137 
00138 private:
00139 
00140   static void assert_p_space( const ModelEvaluator<Scalar> &model, const int l )
00141     {
00142 #ifdef TEUCHOS_DEBUG
00143       const bool p_space_l_is_in_core = model.get_p_space(l)->hasInCoreView();
00144       TEST_FOR_EXCEPTION(
00145         !p_space_l_is_in_core, std::logic_error,
00146         "Error, for the model " << model.description()
00147         << ", the space p_space("<<l<<") must be in-core so that they can"
00148         " act as the domain space for the multi-vector derivative!"
00149         );
00150 #endif
00151     }
00152 
00153 };
00154 
00155 
00156 } // namespace DirectionalFiniteDiffCalculatorTypes
00157 
00158 
00159 // Private static data members
00160 
00161 
00162 template<class Scalar>
00163 const std::string DirectionalFiniteDiffCalculator<Scalar>::FDMethod_name = "FD Method";
00164 template<class Scalar>
00165 const RCP<
00166   Teuchos::StringToIntegralParameterEntryValidator<
00167     Thyra::DirectionalFiniteDiffCalculatorTypes::EFDMethodType
00168   >
00169 >
00170 DirectionalFiniteDiffCalculator<Scalar>::fdMethodValidator
00171 = Teuchos::rcp(
00172   new Teuchos::StringToIntegralParameterEntryValidator<Thyra::DirectionalFiniteDiffCalculatorTypes::EFDMethodType>(
00173     Teuchos::tuple<std::string>(
00174       "order-one"
00175       ,"order-two"
00176       ,"order-two-central"
00177       ,"order-two-auto"
00178       ,"order-four"
00179       ,"order-four-central"
00180       ,"order-four-auto"
00181       )
00182     ,Teuchos::tuple<std::string>(
00183       "Use O(eps) one sided finite differences (cramped bounds)"
00184       ,"Use O(eps^2) one sided finite differences (cramped bounds)"
00185       ,"Use O(eps^2) two sided central finite differences"
00186       ,"Use \"order-two-central\" when not cramped by bounds, otherwise use \"order-two\""
00187       ,"Use O(eps^4) one sided finite differences (cramped bounds)"
00188       ,"Use O(eps^4) two sided central finite differences"
00189       ,"Use \"order-four-central\" when not cramped by bounds, otherwise use \"order-four\""
00190       )
00191     ,Teuchos::tuple<Thyra::DirectionalFiniteDiffCalculatorTypes::EFDMethodType>(
00192       Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_ONE
00193       ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_TWO
00194       ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_TWO_CENTRAL
00195       ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_TWO_AUTO
00196       ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_FOUR
00197       ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_FOUR_CENTRAL
00198       ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_ORDER_FOUR_AUTO
00199       )
00200     ,""
00201     )
00202   );
00203 template<class Scalar>
00204 const std::string DirectionalFiniteDiffCalculator<Scalar>::FDMethod_default = "order-one";
00205 
00206 template<class Scalar>
00207 const std::string DirectionalFiniteDiffCalculator<Scalar>::FDStepLength_name = "FD Step Length";
00208 template<class Scalar>
00209 const double DirectionalFiniteDiffCalculator<Scalar>::FDStepLength_default = -1.0;
00210 
00211 template<class Scalar>
00212 const std::string DirectionalFiniteDiffCalculator<Scalar>::FDStepSelectType_name = "FD Step Select Type";
00213 template<class Scalar>
00214 const RCP<
00215   Teuchos::StringToIntegralParameterEntryValidator<
00216     Thyra::DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType
00217     >
00218   >
00219 DirectionalFiniteDiffCalculator<Scalar>::fdStepSelectTypeValidator
00220 = Teuchos::rcp(
00221   new Teuchos::StringToIntegralParameterEntryValidator<Thyra::DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType>(
00222     Teuchos::tuple<std::string>(
00223       "Absolute"
00224       ,"Relative"
00225       )
00226     ,Teuchos::tuple<std::string>(
00227       "Use absolute step size \""+FDStepLength_name+"\""
00228       ,"Use relative step size \""+FDStepLength_name+"\"*||xo||inf"
00229       )
00230     ,Teuchos::tuple<Thyra::DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType>(
00231       Thyra::DirectionalFiniteDiffCalculatorTypes::FD_STEP_ABSOLUTE
00232       ,Thyra::DirectionalFiniteDiffCalculatorTypes::FD_STEP_RELATIVE
00233       )
00234     ,""
00235     )
00236   );
00237 template<class Scalar>
00238 const std::string DirectionalFiniteDiffCalculator<Scalar>::FDStepSelectType_default = "Absolute";
00239 
00240 
00241 // Constructors/initializer
00242 
00243 
00244 template<class Scalar>
00245 DirectionalFiniteDiffCalculator<Scalar>::DirectionalFiniteDiffCalculator(
00246   EFDMethodType               fd_method_type_in
00247   ,EFDStepSelectType          fd_step_select_type_in
00248   ,ScalarMag                  fd_step_size_in
00249   ,ScalarMag                  fd_step_size_min_in
00250   )
00251   :fd_method_type_(fd_method_type_in)
00252   ,fd_step_select_type_(fd_step_select_type_in)
00253   ,fd_step_size_(fd_step_size_in)
00254   ,fd_step_size_min_(fd_step_size_min_in)
00255 {}
00256 
00257 
00258 // Overriden from ParameterListAcceptor
00259 
00260 
00261 template<class Scalar>
00262 void DirectionalFiniteDiffCalculator<Scalar>::setParameterList(
00263   RCP<ParameterList> const& paramList
00264   )
00265 {
00266   TEST_FOR_EXCEPT(paramList.get()==0);
00267   paramList->validateParameters(*getValidParameters());
00268   paramList_ = paramList;
00269   fd_method_type_ = fdMethodValidator->getIntegralValue(
00270     *paramList_,FDMethod_name,FDMethod_default);
00271   fd_step_select_type_ = fdStepSelectTypeValidator->getIntegralValue(
00272     *paramList_,FDStepSelectType_name,FDStepSelectType_default);
00273   fd_step_size_ = paramList_->get(
00274     FDStepLength_name,FDStepLength_default);
00275   Teuchos::readVerboseObjectSublist(&*paramList_,this);
00276 }
00277 
00278 
00279 template<class Scalar>
00280 RCP<ParameterList>
00281 DirectionalFiniteDiffCalculator<Scalar>::getNonconstParameterList()
00282 {
00283   return paramList_;
00284 }
00285 
00286 
00287 template<class Scalar>
00288 RCP<ParameterList>
00289 DirectionalFiniteDiffCalculator<Scalar>::unsetParameterList()
00290 {
00291   RCP<ParameterList> _paramList = paramList_;
00292   paramList_ = Teuchos::null;
00293   return _paramList;
00294 }
00295 
00296 
00297 template<class Scalar>
00298 RCP<const ParameterList>
00299 DirectionalFiniteDiffCalculator<Scalar>::getParameterList() const
00300 {
00301   return paramList_;
00302 }
00303 
00304 
00305 template<class Scalar>
00306 RCP<const ParameterList>
00307 DirectionalFiniteDiffCalculator<Scalar>::getValidParameters() const
00308 {
00309   using Teuchos::rcp_implicit_cast;
00310   typedef Teuchos::ParameterEntryValidator PEV;
00311   static RCP<ParameterList> pl;
00312   if(pl.get()==NULL) {
00313     pl = Teuchos::parameterList();
00314     pl->set(
00315       FDMethod_name, FDMethod_default,
00316       "The method used to compute the finite differences.",
00317       rcp_implicit_cast<const PEV>(fdMethodValidator)
00318       );
00319     pl->set(
00320       FDStepSelectType_name,FDStepSelectType_default,
00321       "Method used to select the finite difference step length.",
00322       rcp_implicit_cast<const PEV>(fdStepSelectTypeValidator)
00323       );
00324     pl->set(
00325       FDStepLength_name,FDStepLength_default
00326       ,"The length of the finite difference step to take.\n"
00327       "A value of < 0.0 means that the step length will be determined automatically."
00328       );
00329     Teuchos::setupVerboseObjectSublist(&*pl);
00330   }
00331   return pl;
00332 }
00333 
00334 
00335 template<class Scalar>
00336 ModelEvaluatorBase::OutArgs<Scalar>
00337 DirectionalFiniteDiffCalculator<Scalar>::createOutArgs(
00338   const ModelEvaluator<Scalar> &model,
00339   const SelectedDerivatives   &fdDerivatives
00340   )
00341 {
00342   return DirectionalFiniteDiffCalculatorTypes::OutArgsCreator<Scalar>::createOutArgs(
00343     model, fdDerivatives );
00344 }
00345 
00346 
00347 template<class Scalar>
00348 void DirectionalFiniteDiffCalculator<Scalar>::calcVariations(
00349   const ModelEvaluator<Scalar> &model,
00350   const ModelEvaluatorBase::InArgs<Scalar> &bp, // basePoint
00351   const ModelEvaluatorBase::InArgs<Scalar> &dir, // directions
00352   const ModelEvaluatorBase::OutArgs<Scalar> &bfunc, // baseFunctionValues
00353   const ModelEvaluatorBase::OutArgs<Scalar> &var // variations
00354   ) const
00355 {
00356 
00357   using std::string;
00358 
00359   TEUCHOS_FUNC_TIME_MONITOR(
00360     string("Thyra::DirectionalFiniteDiffCalculator<")+ST::name()+">::calcVariations(...)"
00361     );
00362   
00363   using std::setw;
00364   using std::endl;
00365   using std::right;
00366   using Teuchos::as;
00367   typedef ModelEvaluatorBase MEB;
00368   namespace DFDCT = DirectionalFiniteDiffCalculatorTypes;
00369   typedef VectorBase<Scalar> V;
00370   typedef RCP<VectorBase<Scalar> > VectorPtr;
00371   
00372   RCP<Teuchos::FancyOStream> out = this->getOStream();
00373   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00374   const bool trace = (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM));
00375   Teuchos::OSTab tab(out);
00376 
00377   if(out.get() && trace)
00378     *out << "\nEntering DirectionalFiniteDiffCalculator<Scalar>::calcVariations(...)\n";
00379 
00380   if(out.get() && trace)
00381     *out
00382       << "\nbasePoint=\n" << describe(bp,verbLevel)
00383       << "\ndirections=\n" << describe(dir,verbLevel)
00384       << "\nbaseFunctionValues=\n" << describe(bfunc,verbLevel)
00385       << "\nvariations=\n" << describe(var,Teuchos::VERB_LOW);
00386 
00387 #ifdef TEUCHOS_DEBUG
00388   TEST_FOR_EXCEPTION(
00389     var.isEmpty(), std::logic_error,
00390     "Error, all of the variations can not be null!"
00391     );
00392 #endif
00393   
00394   //
00395   // To illustrate the theory behind this implementation consider
00396   // the generic multi-variable function h(z) <: R^n -> R.  Now let's
00397   // consider we have the base point zo and the vector v to
00398   // perturb h(z) along.  First form the function h(zo+a*v) and then
00399   // let's compute d(h)/d(a) at a = 0:
00400   // 
00401   // (1) d(h(zo+a*v))/d(a) at a = 0
00402   //         = sum( d(h)/d(x(i)) * d(x(i))/d(a), i = 1...n)
00403   //         = sum( d(h)/d(x(i)) * v(i), i = 1...n)
00404   //         = d(h)/d(a) * v
00405   //
00406   // Now we can approximate (1) using, for example, central differences as:
00407   // 
00408   // (2) d(h(zo+a*v))/d(a) at a = 0
00409   //          \approx ( h(zo+h*v) - h(zo+h*v) ) / (2*h)
00410   //
00411   // If we equate (1) and (2) we have the approximation:
00412   // 
00413   // (3) d(h)/d(a) * v \approx ( h(zo+h*v) - g(zo+h*v) ) / (2*h)
00414   // 
00415   // 
00416 
00417   // /////////////////////////////////////////
00418   // Validate the input
00419 
00420   // ToDo: Validate input!
00421 
00422   switch(this->fd_method_type()) {
00423     case DFDCT::FD_ORDER_ONE:
00424       if(out.get()&&trace) *out<<"\nUsing one-sided, first-order finite differences ...\n";
00425       break;
00426     case DFDCT::FD_ORDER_TWO:
00427       if(out.get()&&trace) *out<<"\nUsing one-sided, second-order finite differences ...\n";
00428       break;
00429     case DFDCT::FD_ORDER_TWO_CENTRAL:
00430       if(out.get()&&trace) *out<<"\nUsing second-order central finite differences ...\n";
00431       break;
00432     case DFDCT::FD_ORDER_TWO_AUTO:
00433       if(out.get()&&trace) *out<<"\nUsing auto selection of some second-order finite difference method ...\n";
00434       break;
00435     case DFDCT::FD_ORDER_FOUR:
00436       if(out.get()&&trace) *out<<"\nUsing one-sided, fourth-order finite differences ...\n";
00437       break;
00438     case DFDCT::FD_ORDER_FOUR_CENTRAL:
00439       if(out.get()&&trace) *out<<"\nUsing fourth-order central finite differences ...\n";
00440       break;
00441     case DFDCT::FD_ORDER_FOUR_AUTO:
00442       if(out.get()&&trace) *out<<"\nUsing auto selection of some fourth-order finite difference method ...\n";
00443       break;
00444     default:
00445       TEST_FOR_EXCEPT(true); // Should not get here!
00446   }
00447 
00448   // ////////////////////////
00449   // Find the step size
00450 
00451   //
00452   // Get defaults for the optimal step sizes
00453   //
00454 
00455   const ScalarMag
00456     sqrt_epsilon = SMT::squareroot(SMT::eps()),
00457     u_optimal_1  = sqrt_epsilon,
00458     u_optimal_2  = SMT::squareroot(sqrt_epsilon),
00459     u_optimal_4  = SMT::squareroot(u_optimal_2);
00460 
00461   ScalarMag
00462     bp_norm = SMT::zero();
00463   // ToDo above: compute a reasonable norm somehow based on the base-point vector(s)!
00464   
00465   ScalarMag
00466     uh_opt = 0.0;
00467   switch(this->fd_method_type()) {
00468     case DFDCT::FD_ORDER_ONE:
00469       uh_opt = u_optimal_1 * ( fd_step_select_type() == DFDCT::FD_STEP_ABSOLUTE ? 1.0 : bp_norm + 1.0 );
00470       break;
00471     case DFDCT::FD_ORDER_TWO:
00472     case DFDCT::FD_ORDER_TWO_CENTRAL:
00473     case DFDCT::FD_ORDER_TWO_AUTO:
00474       uh_opt = u_optimal_2 * ( fd_step_select_type() == DFDCT::FD_STEP_ABSOLUTE ? 1.0 : bp_norm + 1.0 );
00475       break;
00476     case DFDCT::FD_ORDER_FOUR:
00477     case DFDCT::FD_ORDER_FOUR_CENTRAL:
00478     case DFDCT::FD_ORDER_FOUR_AUTO:
00479       uh_opt = u_optimal_4 * ( fd_step_select_type() == DFDCT::FD_STEP_ABSOLUTE ? 1.0 : bp_norm + 1.0 );
00480       break;
00481     default:
00482       TEST_FOR_EXCEPT(true); // Should not get here!
00483   }
00484 
00485   if(out.get()&&trace) *out<<"\nDefault optimal step length uh_opt = " << uh_opt << " ...\n";
00486 
00487   //
00488   // Set the step sizes used.
00489   //
00490 
00491   ScalarMag
00492     uh      = this->fd_step_size();
00493 
00494   if( uh < 0 )
00495     uh = uh_opt;
00496   else if( fd_step_select_type() == DFDCT::FD_STEP_RELATIVE )
00497     uh *= (bp_norm + 1.0);
00498 
00499   if(out.get()&&trace) *out<<"\nStep size to be used uh="<<uh<<"\n";
00500 
00501   //
00502   // Find step lengths that stay in bounds!
00503   //
00504 
00505   // ToDo: Add logic for variable bounds when needed!
00506 
00507   //
00508   // Set the actual method being used
00509   //
00510   // ToDo: Consider cramped bounds and method order.
00511   //
00512   
00513   DFDCT::EFDMethodType l_fd_method_type = this->fd_method_type();
00514   switch(l_fd_method_type) {
00515     case DFDCT::FD_ORDER_TWO_AUTO:
00516       l_fd_method_type = DFDCT::FD_ORDER_TWO_CENTRAL;
00517       break;
00518     case DFDCT::FD_ORDER_FOUR_AUTO:
00519       l_fd_method_type = DFDCT::FD_ORDER_FOUR_CENTRAL;
00520       break;
00521     default:
00522       break; // Okay
00523   }
00524 
00525   //if(out.get()&&trace) *out<<"\nStep size to fit in bounds: uh="<<uh"\n";
00526 
00527   int p_saved = -1;
00528   if(out.get())
00529     p_saved = out->precision();
00530 
00531   // ///////////////////////////////////////////////
00532   // Compute the directional derivatives
00533 
00534   const int Np = var.Np(), Ng = var.Ng();
00535   
00536   // Setup storage for perturbed variables
00537   VectorPtr per_x, per_x_dot;
00538   std::vector<VectorPtr> per_p(Np);
00539   MEB::InArgs<Scalar> pp = model.createInArgs();
00540   pp.setArgs(bp); // Set all args to start with
00541   if( bp.supports(MEB::IN_ARG_x) ) {
00542     if( dir.get_x().get() )
00543       pp.set_x(per_x=createMember(model.get_x_space()));
00544     else
00545       pp.set_x(bp.get_x());
00546   }
00547   if( bp.supports(MEB::IN_ARG_x_dot) ) {
00548     if( dir.get_x_dot().get() )
00549       pp.set_x_dot(per_x_dot=createMember(model.get_x_space()));
00550     else
00551       pp.set_x_dot(bp.get_x_dot());
00552   }
00553   for( int l = 0; l < Np; ++l ) {
00554     if( dir.get_p(l).get() )
00555       pp.set_p(l,per_p[l]=createMember(model.get_p_space(l)));
00556     else
00557       pp.set_p(l,bp.get_p(l));
00558   }
00559   if(out.get() && trace)
00560     *out
00561       << "\nperturbedPoint after initial setup (with some unintialized vectors) =\n"
00562       << describe(pp,verbLevel);
00563   
00564   // Setup storage for perturbed functions
00565   bool all_funcs_at_base_computed = true;
00566   MEB::OutArgs<Scalar> pfunc = model.createOutArgs();
00567   {
00568     VectorPtr f;
00569     if( var.supports(MEB::OUT_ARG_f) && (f=var.get_f()).get() ) {
00570       pfunc.set_f(createMember(model.get_f_space()));
00571       assign(f.ptr(),ST::zero());
00572       if(!bfunc.get_f().get()) all_funcs_at_base_computed = false;
00573     }
00574     for( int j = 0; j < Ng; ++j ) {
00575       VectorPtr g_j;
00576       if( (g_j=var.get_g(j)).get() ) {
00577         pfunc.set_g(j,createMember(model.get_g_space(j)));
00578         assign(g_j.ptr(),ST::zero());
00579         if(!bfunc.get_g(j).get()) all_funcs_at_base_computed = false;
00580       }
00581     }
00582   }
00583   if(out.get() && trace)
00584     *out
00585       << "\nperturbedFunctions after initial setup (with some unintialized vectors) =\n"
00586       << describe(pfunc,verbLevel);
00587   
00588   const int dbl_p = 15;
00589   if(out.get())
00590     *out << std::setprecision(dbl_p);
00591     
00592   //
00593   // Compute the weighted sum of the terms
00594   //
00595     
00596   int num_evals = 0;
00597   ScalarMag dwgt = SMT::zero();
00598   switch(l_fd_method_type) {
00599     case DFDCT::FD_ORDER_ONE: // may only need one eval if f(xo) etc is passed in
00600       num_evals = 2;
00601       dwgt      = ScalarMag(1.0);
00602       break;
00603     case DFDCT::FD_ORDER_TWO: // may only need two evals if c(xo) etc is passed in
00604       num_evals = 3;
00605       dwgt      = ScalarMag(2.0);
00606       break;
00607     case DFDCT::FD_ORDER_TWO_CENTRAL:
00608       num_evals = 2;
00609       dwgt      = ScalarMag(2.0);
00610       break;
00611     case DFDCT::FD_ORDER_FOUR:
00612       num_evals = 5;
00613       dwgt      = ScalarMag(12.0);
00614       break;
00615     case DFDCT::FD_ORDER_FOUR_CENTRAL:
00616       num_evals = 4;
00617       dwgt      = ScalarMag(12.0);
00618       break;
00619     default:
00620       TEST_FOR_EXCEPT(true); // Should not get here!
00621   }
00622   for( int eval_i = 1; eval_i <= num_evals; ++eval_i ) {
00623     // Set the step constant and the weighting constant
00624     ScalarMag
00625       uh_i   = 0.0,
00626       wgt_i  = 0.0;
00627     switch(l_fd_method_type) {
00628       case DFDCT::FD_ORDER_ONE: {
00629         switch(eval_i) {
00630           case 1:
00631             uh_i  = +0.0;
00632             wgt_i = -1.0;
00633             break;
00634           case 2:
00635             uh_i  = +1.0;
00636             wgt_i = +1.0;
00637             break;
00638         }
00639         break;
00640       }
00641       case DFDCT::FD_ORDER_TWO: {
00642         switch(eval_i) {
00643           case 1:
00644             uh_i  = +0.0;
00645             wgt_i = -3.0;
00646             break;
00647           case 2:
00648             uh_i  = +1.0;
00649             wgt_i = +4.0;
00650             break;
00651           case 3:
00652             uh_i  = +2.0;
00653             wgt_i = -1.0;
00654             break;
00655         }
00656         break;
00657       }
00658       case DFDCT::FD_ORDER_TWO_CENTRAL: {
00659         switch(eval_i) {
00660           case 1:
00661             uh_i  = -1.0;
00662             wgt_i = -1.0;
00663             break;
00664           case 2:
00665             uh_i  = +1.0;
00666             wgt_i = +1.0;
00667             break;
00668         }
00669         break;
00670       }
00671       case DFDCT::FD_ORDER_FOUR: {
00672         switch(eval_i) {
00673           case 1:
00674             uh_i  = +0.0;
00675             wgt_i = -25.0;
00676             break;
00677           case 2:
00678             uh_i  = +1.0;
00679             wgt_i = +48.0;
00680             break;
00681           case 3:
00682             uh_i  = +2.0;
00683             wgt_i = -36.0;
00684             break;
00685           case 4:
00686             uh_i  = +3.0;
00687             wgt_i = +16.0;
00688             break;
00689           case 5:
00690             uh_i  = +4.0;
00691             wgt_i = -3.0;
00692             break;
00693         }
00694         break;
00695       }
00696       case DFDCT::FD_ORDER_FOUR_CENTRAL: {
00697         switch(eval_i) {
00698           case 1:
00699             uh_i  = -2.0;
00700             wgt_i = +1.0;
00701             break;
00702           case 2:
00703             uh_i  = -1.0;
00704             wgt_i = -8.0;
00705             break;
00706           case 3:
00707             uh_i  = +1.0;
00708             wgt_i = +8.0;
00709             break;
00710           case 4:
00711             uh_i  = +2.0;
00712             wgt_i = -1.0;
00713             break;
00714         }
00715         break;
00716       }
00717       case DFDCT::FD_ORDER_TWO_AUTO:
00718       case DFDCT::FD_ORDER_FOUR_AUTO:
00719         break; // Okay
00720       default:
00721         TEST_FOR_EXCEPT(true);
00722     }
00723 
00724     if(out.get() && trace)
00725       *out << "\neval_i="<<eval_i<<", uh_i="<<uh_i<<", wgt_i="<<wgt_i<<"\n";
00726     Teuchos::OSTab tab2(out);
00727     
00728     // Compute the weighted term and add it to the sum
00729     if(uh_i == ST::zero()) {
00730       MEB::OutArgs<Scalar> bfuncall;
00731       if(!all_funcs_at_base_computed) {
00732         // Compute missing functions at the base point
00733         bfuncall = model.createOutArgs();
00734         if( pfunc.supports(MEB::OUT_ARG_f) && pfunc.get_f().get() && !bfunc.get_f().get() ) {
00735           bfuncall.set_f(createMember(model.get_f_space()));
00736         }
00737         for( int j = 0; j < Ng; ++j ) {
00738           if( pfunc.get_g(j).get() && !bfunc.get_g(j).get() ) {
00739             bfuncall.set_g(j,createMember(model.get_g_space(j)));
00740           }
00741         }
00742         model.evalModel(bp,bfuncall);
00743         bfuncall.setArgs(bfunc);
00744       }
00745       else {
00746         bfuncall = bfunc;
00747       }
00748       // Use functions at the base point
00749       if(out.get() && trace)
00750         *out << "\nSetting variations = wgt_i * basePoint ...\n";
00751       VectorPtr f;
00752       if( pfunc.supports(MEB::OUT_ARG_f) && (f=var.get_f()).get() ) {
00753         V_StV<Scalar>(f.ptr(), wgt_i, *bfuncall.get_f());
00754       }
00755       for( int j = 0; j < Ng; ++j ) {
00756         VectorPtr g_j;
00757         if( (g_j=var.get_g(j)).get() ) {
00758           V_StV<Scalar>(g_j.ptr(), wgt_i, *bfuncall.get_g(j));
00759         }
00760       }
00761     }
00762     else {
00763       if(out.get() && trace)
00764         *out << "\nSetting perturbedPoint = basePoint + uh_i*uh*direction ...\n";
00765       // z = zo + uh_i*uh*v
00766       {
00767         if ( dir.supports(MEB::IN_ARG_x) && dir.get_x().get() )
00768           V_StVpV(per_x.ptr(),as<Scalar>(uh_i*uh),*dir.get_x(),*bp.get_x());
00769         if ( dir.supports(MEB::IN_ARG_x_dot) && dir.get_x_dot().get() )
00770           V_StVpV(per_x_dot.ptr(), as<Scalar>(uh_i*uh), *dir.get_x_dot(), *bp.get_x_dot());
00771         for ( int l = 0; l < Np; ++l ) {
00772           if( dir.get_p(l).get() )
00773             V_StVpV(per_p[l].ptr(), as<Scalar>(uh_i*uh), *dir.get_p(l), *bp.get_p(l));
00774         }
00775       }
00776       if(out.get() && trace)
00777         *out << "\nperturbedPoint=\n" << describe(pp,verbLevel);
00778       // Compute perturbed function values h(zo+uh_i*uh)
00779       if(out.get() && trace)
00780         *out << "\nCompute perturbedFunctions at perturbedPoint...\n";
00781       model.evalModel(pp,pfunc);
00782       if(out.get() && trace)
00783         *out << "\nperturbedFunctions=\n" << describe(pfunc,verbLevel);
00784       // Sum perturbed function values into total variation
00785       {
00786         // var_h += wgt_i*perturbed_h
00787         if(out.get() && trace)
00788           *out << "\nComputing variations += wgt_i*perturbedfunctions ...\n";
00789         VectorPtr f;
00790         if( pfunc.supports(MEB::OUT_ARG_f) && (f=pfunc.get_f()).get() )
00791           Vp_StV<Scalar>(var.get_f().ptr(), wgt_i, *f);
00792         for( int j = 0; j < Ng; ++j ) {
00793           VectorPtr g_j;
00794           if( (g_j=pfunc.get_g(j)).get() )
00795             Vp_StV<Scalar>(var.get_g(j).ptr(), wgt_i, *g_j);
00796         }
00797       }
00798     }
00799     if(out.get() && trace)
00800       *out << "\nvariations=\n" << describe(var,verbLevel);
00801   }
00802 
00803   //
00804   // Multiply by the scaling factor!
00805   //
00806   
00807   {
00808     // var_h *= 1.0/(dwgt*uh)
00809     const Scalar alpha = ST::one()/(dwgt*uh);
00810     if(out.get() && trace)
00811       *out
00812         << "\nComputing variations *= (1.0)/(dwgt*uh),"
00813         << " where (1.0)/(dwgt*uh) = (1.0)/("<<dwgt<<"*"<<uh<<") = "<<alpha<<" ...\n";
00814     VectorPtr f;
00815     if( var.supports(MEB::OUT_ARG_f) && (f=var.get_f()).get() )
00816       Vt_S(f.ptr(),alpha);
00817     for( int j = 0; j < Ng; ++j ) {
00818       VectorPtr g_j;
00819       if( (g_j=var.get_g(j)).get() )
00820         Vt_S(g_j.ptr(),alpha);
00821     }
00822     if(out.get() && trace)
00823       *out << "\nFinal variations=\n" << describe(var,verbLevel);
00824   }
00825   
00826   if(out.get())
00827     *out << std::setprecision(p_saved);
00828 
00829   if(out.get() && trace)
00830     *out << "\nLeaving DirectionalFiniteDiffCalculator<Scalar>::calcVariations(...)\n";
00831   
00832 }
00833 
00834 
00835 template<class Scalar>
00836 void DirectionalFiniteDiffCalculator<Scalar>::calcDerivatives(
00837   const ModelEvaluator<Scalar> &model,
00838   const ModelEvaluatorBase::InArgs<Scalar> &bp, // basePoint
00839   const ModelEvaluatorBase::OutArgs<Scalar> &bfunc, // baseFunctionValues
00840   const ModelEvaluatorBase::OutArgs<Scalar> &deriv // derivatives
00841   ) const
00842 {
00843 
00844   using std::string;
00845   //typedef Teuchos::ScalarTraits<Scalar> ST;
00846 
00847   TEUCHOS_FUNC_TIME_MONITOR(
00848     string("Thyra::DirectionalFiniteDiffCalculator<")+ST::name()+">::calcDerivatives(...)"
00849     );
00850 
00851   typedef ModelEvaluatorBase MEB;
00852   typedef RCP<VectorBase<Scalar> > VectorPtr;
00853   typedef RCP<MultiVectorBase<Scalar> > MultiVectorPtr;
00854 
00855   RCP<Teuchos::FancyOStream> out = this->getOStream();
00856   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00857   const bool trace = (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM));
00858   Teuchos::OSTab tab(out);
00859 
00860   if(out.get() && trace)
00861     *out << "\nEntering DirectionalFiniteDiffCalculator<Scalar>::calcDerivatives(...)\n";
00862 
00863   if(out.get() && trace)
00864     *out
00865       << "\nbasePoint=\n" << describe(bp,verbLevel)
00866       << "\nbaseFunctionValues=\n" << describe(bfunc,verbLevel)
00867       << "\nderivatives=\n" << describe(deriv,Teuchos::VERB_LOW);
00868   
00869   //
00870   // We will only compute finite differences w.r.t. p(l) for now
00871   //
00872   const int Np = bp.Np(), Ng = bfunc.Ng();
00873   MEB::InArgs<Scalar> dir = model.createInArgs();
00874   MEB::OutArgs<Scalar> var = model.createOutArgs();
00875   MultiVectorPtr DfDp_l;
00876   std::vector<MEB::DerivativeMultiVector<Scalar> > DgDp_l(Ng);
00877   std::vector<VectorPtr> var_g(Ng);
00878   for( int l = 0; l < Np; ++l ) {
00879     if(out.get() && trace)
00880       *out << "\nComputing derivatives for parameter subvector p("<<l<<") ...\n";
00881     Teuchos::OSTab tab2(out);
00882     //
00883     // Load up OutArgs var object of derivative variations to compute
00884     //
00885     bool hasDerivObject = false;
00886     // DfDp(l)
00887     if (
00888       !deriv.supports(MEB::OUT_ARG_DfDp,l).none()
00889       && !deriv.get_DfDp(l).isEmpty()
00890       )
00891     {
00892       hasDerivObject = true;
00893       std::ostringstream name; name << "DfDp("<<l<<")";
00894       DfDp_l = get_mv(deriv.get_DfDp(l),name.str(),MEB::DERIV_MV_BY_COL);
00895     }
00896     else {
00897       DfDp_l = Teuchos::null;
00898     }
00899     // DgDp(j=1...Ng,l)
00900     for ( int j = 0; j < Ng; ++j ) {
00901       if (
00902         !deriv.supports(MEB::OUT_ARG_DgDp,j,l).none()
00903         &&
00904         !deriv.get_DgDp(j,l).isEmpty()
00905         )
00906       {
00907         hasDerivObject = true;
00908         std::ostringstream name; name << "DgDp("<<j<<","<<l<<")";
00909         DgDp_l[j] = get_dmv(deriv.get_DgDp(j,l),name.str());
00910         if( DgDp_l[j].getMultiVector().get() && !var_g[j].get() )
00911         {
00912           // Need a temporary vector for the variation for g(j)
00913           var_g[j] = createMember(model.get_g_space(j));
00914         }
00915       }
00916       else{
00917         DgDp_l[j] = MEB::DerivativeMultiVector<Scalar>();
00918         var_g[j] = Teuchos::null;
00919       }
00920     }
00921     //
00922     // Compute derivative variations by finite differences
00923     //
00924     if (hasDerivObject) {
00925       VectorPtr e_i = createMember(model.get_p_space(l));
00926       dir.set_p(l,e_i);
00927       assign(e_i.ptr(),ST::zero());
00928       const int np_l = e_i->space()->dim();
00929       for( int i = 0 ; i < np_l; ++ i ) {
00930         if(out.get() && trace)
00931           *out << "\nComputing derivatives for single variable p("<<l<<")("<<i<<") ...\n";
00932         Teuchos::OSTab tab3(out);
00933         if(DfDp_l.get()) var.set_f(DfDp_l->col(i)); // Compute DfDp(l)(i) in place!
00934         for(int j = 0; j < Ng; ++j) {
00935           MultiVectorPtr DgDp_j_l;
00936           if( (DgDp_j_l=DgDp_l[j].getMultiVector()).get() ) {
00937             var.set_g(j,var_g[j]); // Computes d(g(j))/d(p(l)(i))
00938           }
00939         }
00940         set_ele(i,ST::one(),e_i.ptr());
00941         this->calcVariations(
00942           model,bp,dir,bfunc,var
00943           );
00944         set_ele(i,ST::zero(),e_i.ptr());
00945         if (DfDp_l.get()) var.set_f(Teuchos::null);
00946         for (int j = 0; j < Ng; ++j) {
00947           MultiVectorPtr DgDp_j_l;
00948           if ( !is_null(DgDp_j_l=DgDp_l[j].getMultiVector()) ) {
00949             assign( DgDp_j_l->col(i).ptr(), *var_g[j] );
00950           }
00951         }
00952       }
00953       dir.set_p(l,Teuchos::null);
00954     }
00955   }
00956     
00957   if(out.get() && trace)
00958     *out
00959       << "\nderivatives=\n" << describe(deriv,verbLevel);
00960   
00961   if(out.get() && trace)
00962     *out << "\nLeaving DirectionalFiniteDiffCalculator<Scalar>::calcDerivatives(...)\n";
00963 
00964 }
00965 
00966 
00967 } // namespace Thyra
00968 
00969 
00970 #endif  // THYRA_DIRECTIONAL_FINITE_DIFF_CALCULATOR_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines