Thyra_DirectionalFiniteDiffCalculator.hpp

00001 
00002 #ifndef THYRA_DIRECTIONAL_FINITE_DIFF_CALCULATOR_HPP
00003 #define THYRA_DIRECTIONAL_FINITE_DIFF_CALCULATOR_HPP
00004 
00005 #include "Thyra_ModelEvaluator.hpp"
00006 #include "Thyra_ModelEvaluatorHelpers.hpp"
00007 #include "Thyra_DetachedVectorView.hpp"
00008 #include "Thyra_DetachedMultiVectorView.hpp"
00009 #include "Teuchos_VerboseObject.hpp"
00010 #include "Teuchos_ParameterListAcceptor.hpp"
00011 #include "Teuchos_StandardMemberCompositionMacros.hpp"
00012 
00013 namespace Thyra {
00014 
00015 namespace DirectionalFiniteDiffCalculatorTypes {
00016 
00018 enum EFDMethodType {
00019   FD_ORDER_ONE           
00020   ,FD_ORDER_TWO          
00021   ,FD_ORDER_TWO_CENTRAL  
00022   ,FD_ORDER_TWO_AUTO     
00023   ,FD_ORDER_FOUR         
00024   ,FD_ORDER_FOUR_CENTRAL 
00025   ,FD_ORDER_FOUR_AUTO    
00026 };
00027 
00029 enum EFDStepSelectType {
00030   FD_STEP_ABSOLUTE      
00031   ,FD_STEP_RELATIVE     
00032 };
00033 
00034 } // namespace DirectionalFiniteDiffCalculatorTypes
00035 
00050 template<class Scalar>
00051 class DirectionalFiniteDiffCalculator
00052   : public Teuchos::VerboseObject<DirectionalFiniteDiffCalculator<Scalar> >
00053 //, public Teuchos::ParameterListAccetor
00054 {
00055 public:
00056 
00058   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00059 
00061   typedef DirectionalFiniteDiffCalculatorTypes::EFDMethodType EFDMethodType;
00062 
00064   typedef DirectionalFiniteDiffCalculatorTypes::EFDStepSelectType EFDStepSelectType;
00065 
00067   STANDARD_MEMBER_COMPOSITION_MEMBERS( EFDMethodType, fd_method_type )
00068 
00069   
00070   STANDARD_MEMBER_COMPOSITION_MEMBERS( EFDStepSelectType, fd_step_select_type )
00071 
00078   STANDARD_MEMBER_COMPOSITION_MEMBERS( ScalarMag, fd_step_size )
00079 
00088   STANDARD_MEMBER_COMPOSITION_MEMBERS( ScalarMag, fd_step_size_min )
00089 
00091   DirectionalFiniteDiffCalculator(
00092     EFDMethodType               fd_method_type        = DirectionalFiniteDiffCalculatorTypes::FD_ORDER_FOUR_AUTO
00093     ,EFDStepSelectType          fd_step_select_type   = DirectionalFiniteDiffCalculatorTypes::FD_STEP_ABSOLUTE
00094     ,ScalarMag                  fd_step_size          = -1.0
00095     ,ScalarMag                  fd_step_size_min      = -1.0
00096     );
00097 
00110   void calcVariations(
00111     const ModelEvaluator<Scalar>                 &model
00112     ,const ModelEvaluatorBase::InArgs<Scalar>    &basePoint
00113     ,const ModelEvaluatorBase::InArgs<Scalar>    &directions
00114     ,const ModelEvaluatorBase::OutArgs<Scalar>   &baseFunctionValues
00115     ,const ModelEvaluatorBase::OutArgs<Scalar>   &variations
00116     ) const;
00117 
00120   void calcDerivatives(
00121     const ModelEvaluator<Scalar>                 &model
00122     ,const ModelEvaluatorBase::InArgs<Scalar>    &basePoint
00123     ,const ModelEvaluatorBase::OutArgs<Scalar>   &baseFunctionValues
00124     ,const ModelEvaluatorBase::OutArgs<Scalar>   &derivatives
00125     ) const;
00126 
00127 };
00128 
00129 // //////////////////////////////
00130 // Implementations
00131 
00132 template<class Scalar>
00133 DirectionalFiniteDiffCalculator<Scalar>::DirectionalFiniteDiffCalculator(
00134   EFDMethodType               fd_method_type
00135   ,EFDStepSelectType          fd_step_select_type
00136   ,ScalarMag                  fd_step_size
00137   ,ScalarMag                  fd_step_size_min
00138   )
00139   :fd_method_type_(fd_method_type)
00140   ,fd_step_select_type_(fd_step_select_type)
00141   ,fd_step_size_(fd_step_size)
00142   ,fd_step_size_min_(fd_step_size_min)
00143 {}
00144 
00145 template<class Scalar>
00146 void DirectionalFiniteDiffCalculator<Scalar>::calcVariations(
00147   const ModelEvaluator<Scalar>                 &model
00148   ,const ModelEvaluatorBase::InArgs<Scalar>    &bp         // basePoint
00149   ,const ModelEvaluatorBase::InArgs<Scalar>    &dir        // directions
00150   ,const ModelEvaluatorBase::OutArgs<Scalar>   &bfunc      // baseFunctionValues
00151   ,const ModelEvaluatorBase::OutArgs<Scalar>   &var        // variations
00152   ) const
00153 {
00154   
00155   using std::setw;
00156   using std::endl;
00157   using std::right;
00158   typedef Teuchos::ScalarTraits<Scalar> ST;
00159   typedef typename ST::magnitudeType ScalarMag;
00160   typedef ModelEvaluatorBase MEB;
00161   namespace DFDCT = DirectionalFiniteDiffCalculatorTypes;
00162   typedef VectorBase<Scalar> V;
00163   typedef Teuchos::RefCountPtr<VectorBase<Scalar> > VectorPtr;
00164   
00165   Teuchos::RefCountPtr<Teuchos::FancyOStream> out = this->getOStream();
00166   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00167   const bool trace = (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM));
00168   Teuchos::OSTab tab(out);
00169 
00170   if(out.get() && trace)
00171     *out << "\nEntering DirectionalFiniteDiffCalculator<Scalar>::calcVariations(...)\n";
00172 
00173   if(out.get() && trace)
00174     *out
00175       << "\nbasePoint=\n" << describe(bp,verbLevel)
00176       << "\ndirections=\n" << describe(dir,verbLevel)
00177       << "\nbaseFunctionValues=\n" << describe(bfunc,verbLevel)
00178       << "\nvariations=\n" << describe(var,Teuchos::VERB_LOW);
00179   
00180   //
00181   // To illustrate the theory behind this implementation consider
00182   // the generic multi-variable function h(z) <: R^n -> R.  Now let's
00183   // consider we have the base point zo and the vector v to
00184   // perturb h(z) along.  First form the function h(zo+a*v) and then
00185   // let's compute d(h)/d(a) at a = 0:
00186   // 
00187   // (1) d(h(zo+a*v))/d(a) at a = 0
00188   //         = sum( d(h)/d(x(i)) * d(x(i))/d(a), i = 1...n)
00189   //         = sum( d(h)/d(x(i)) * v(i), i = 1...n)
00190   //         = d(h)/d(a) * v
00191   //
00192   // Now we can approximate (1) using, for example, central differences as:
00193   // 
00194   // (2) d(h(zo+a*v))/d(a) at a = 0
00195   //          \approx ( h(zo+h*v) - h(zo+h*v) ) / (2*h)
00196   //
00197   // If we equate (1) and (2) we have the approximation:
00198   // 
00199   // (3) d(h)/d(a) * v \approx ( h(zo+h*v) - g(zo+h*v) ) / (2*h)
00200   // 
00201   // 
00202 
00203   // /////////////////////////////////////////
00204   // Validate the input
00205 
00206   // ToDo: Validate input!
00207 
00208   switch(this->fd_method_type()) {
00209     case DFDCT::FD_ORDER_ONE:
00210       if(out.get()&&trace) *out<<"\nUsing one-sided, first-order finite differences ...\n";
00211       break;
00212     case DFDCT::FD_ORDER_TWO:
00213       if(out.get()&&trace) *out<<"\nUsing one-sided, second-order finite differences ...\n";
00214       break;
00215     case DFDCT::FD_ORDER_TWO_CENTRAL:
00216       if(out.get()&&trace) *out<<"\nUsing second-order central finite differences ...\n";
00217       break;
00218     case DFDCT::FD_ORDER_TWO_AUTO:
00219       if(out.get()&&trace) *out<<"\nUsing auto selection of some second-order finite difference method ...\n";
00220       break;
00221     case DFDCT::FD_ORDER_FOUR:
00222       if(out.get()&&trace) *out<<"\nUsing one-sided, fourth-order finite differences ...\n";
00223       break;
00224     case DFDCT::FD_ORDER_FOUR_CENTRAL:
00225       if(out.get()&&trace) *out<<"\nUsing fourth-order central finite differences ...\n";
00226       break;
00227     case DFDCT::FD_ORDER_FOUR_AUTO:
00228       if(out.get()&&trace) *out<<"\nUsing auto selection of some fourth-order finite difference method ...\n";
00229       break;
00230     default:
00231       TEST_FOR_EXCEPT(true); // Should not get here!
00232   }
00233 
00234   // ////////////////////////
00235   // Find the step size
00236 
00237   //
00238   // Get defaults for the optimal step sizes
00239   //
00240 
00241   const ScalarMag
00242     sqrt_epsilon = ST::squareroot(ST::eps()),
00243     u_optimal_1  = sqrt_epsilon,
00244     u_optimal_2  = ST::squareroot(sqrt_epsilon),
00245     u_optimal_4  = ST::squareroot(u_optimal_2);
00246 
00247   ScalarMag
00248     bp_norm = ST::zero();
00249   // ToDo above: compute a reasonable norm somehow based on the base-point vector(s)!
00250   
00251   ScalarMag
00252     uh_opt = 0.0;
00253   switch(this->fd_method_type()) {
00254     case DFDCT::FD_ORDER_ONE:
00255       uh_opt = u_optimal_1 * ( fd_step_select_type() == DFDCT::FD_STEP_ABSOLUTE ? 1.0 : bp_norm + 1.0 );
00256       break;
00257     case DFDCT::FD_ORDER_TWO:
00258     case DFDCT::FD_ORDER_TWO_CENTRAL:
00259     case DFDCT::FD_ORDER_TWO_AUTO:
00260       uh_opt = u_optimal_2 * ( fd_step_select_type() == DFDCT::FD_STEP_ABSOLUTE ? 1.0 : bp_norm + 1.0 );
00261       break;
00262     case DFDCT::FD_ORDER_FOUR:
00263     case DFDCT::FD_ORDER_FOUR_CENTRAL:
00264     case DFDCT::FD_ORDER_FOUR_AUTO:
00265       uh_opt = u_optimal_4 * ( fd_step_select_type() == DFDCT::FD_STEP_ABSOLUTE ? 1.0 : bp_norm + 1.0 );
00266       break;
00267     default:
00268       TEST_FOR_EXCEPT(true); // Should not get here!
00269   }
00270 
00271   if(out.get()&&trace) *out<<"\nDefault optimal step length uh_opt = " << uh_opt << " ...\n";
00272 
00273   //
00274   // Set the step sizes used.
00275   //
00276 
00277   ScalarMag
00278     uh      = this->fd_step_size();
00279 
00280   if( uh < 0 )
00281     uh = uh_opt;
00282   else if( fd_step_select_type() == DFDCT::FD_STEP_RELATIVE )
00283     uh *= (bp_norm + 1.0);
00284 
00285   if(out.get()&&trace) *out<<"\nStep size to be used uh="<<uh<<"\n";
00286 
00287   //
00288   // Find step lengths that stay in bounds!
00289   //
00290 
00291   // ToDo: Add logic for variable bounds when needed!
00292 
00293   //
00294   // Set the actual method being used
00295   //
00296   // ToDo: Consider cramped bounds and method order.
00297   //
00298   
00299   DFDCT::EFDMethodType fd_method_type = this->fd_method_type();
00300   switch(fd_method_type) {
00301     case DFDCT::FD_ORDER_TWO_AUTO:
00302       fd_method_type = DFDCT::FD_ORDER_TWO_CENTRAL;
00303       break;
00304     case DFDCT::FD_ORDER_FOUR_AUTO:
00305       fd_method_type = DFDCT::FD_ORDER_FOUR_CENTRAL;
00306       break;
00307     default:
00308       break; // Okay
00309   }
00310 
00311   //if(out.get()&&trace) *out<<"\nStep size to fit in bounds: uh="<<uh"\n";
00312 
00313   int p_saved;
00314   if(out.get())
00315     p_saved = out->precision();
00316 
00317   // ///////////////////////////////////////////////
00318   // Compute the directional derivatives
00319 
00320   const int Np = var.Np(), Ng = var.Ng();
00321   
00322   // Setup storage for perturbed variables
00323   VectorPtr per_x, per_x_dot;
00324   std::vector<VectorPtr> per_p(Np);
00325   MEB::InArgs<Scalar> pp = model.createInArgs();
00326   if( bp.supports(MEB::IN_ARG_x) ) {
00327     if( dir.get_x().get() )
00328       pp.set_x(per_x=createMember(model.get_x_space()));
00329     else
00330       pp.set_x(bp.get_x());
00331   }
00332   if( bp.supports(MEB::IN_ARG_x_dot) ) {
00333     if( dir.get_x_dot().get() )
00334       pp.set_x_dot(per_x_dot=createMember(model.get_x_space()));
00335     else
00336       pp.set_x_dot(bp.get_x_dot());
00337   }
00338   for( int l = 0; l < Np; ++l ) {
00339     if( dir.get_p(l).get() )
00340       pp.set_p(l,per_p[l]=createMember(model.get_p_space(l)));
00341     else
00342       pp.set_p(l,bp.get_p(l));
00343   }
00344   if(out.get() && trace)
00345     *out << "\nperturbedPoint after initial setup =\n" << describe(pp,verbLevel);
00346   
00347   // Setup storage for perturbed functions
00348   bool all_funcs_at_base_computed = true;
00349   MEB::OutArgs<Scalar> pfunc = model.createOutArgs();
00350   if(1) {
00351     VectorPtr f;
00352     if( var.supports(MEB::OUT_ARG_f) && (f=var.get_f()).get() ) {
00353       pfunc.set_f(createMember(model.get_f_space()));
00354       assign(&*f,ST::zero());
00355       if(!bfunc.get_f().get()) all_funcs_at_base_computed = false;
00356     }
00357     for( int j = 0; j < Ng; ++j ) {
00358       VectorPtr g_j;
00359       if( (g_j=var.get_g(j)).get() ) {
00360         pfunc.set_g(j,createMember(model.get_g_space(j)));
00361         assign(&*g_j,ST::zero());
00362         if(!bfunc.get_g(j).get()) all_funcs_at_base_computed = false;
00363       }
00364     }
00365   }
00366   if(out.get() && trace)
00367     *out << "\nperturbedFunctions after initial setup =\n" << describe(pfunc,verbLevel);
00368   
00369   const int dbl_p = 15;
00370   if(out.get())
00371     *out << std::setprecision(dbl_p);
00372     
00373   //
00374   // Compute the weighted sum of the terms
00375   //
00376     
00377   int         num_evals  = 0;
00378   ScalarMag   dwgt       = ST::zero();
00379   switch(fd_method_type) {
00380     case DFDCT::FD_ORDER_ONE: // may only need one eval if f(xo) etc is passed in
00381       num_evals = 2;
00382       dwgt      = ScalarMag(1.0);
00383       break;
00384     case DFDCT::FD_ORDER_TWO: // may only need two evals if c(xo) etc is passed in
00385       num_evals = 3;
00386       dwgt      = ScalarMag(2.0);
00387       break;
00388     case DFDCT::FD_ORDER_TWO_CENTRAL:
00389       num_evals = 2;
00390       dwgt      = ScalarMag(2.0);
00391       break;
00392     case DFDCT::FD_ORDER_FOUR:
00393       num_evals = 5;
00394       dwgt      = ScalarMag(12.0);
00395       break;
00396     case DFDCT::FD_ORDER_FOUR_CENTRAL:
00397       num_evals = 4;
00398       dwgt      = ScalarMag(12.0);
00399       break;
00400     default:
00401       TEST_FOR_EXCEPT(true); // Should not get here!
00402   }
00403   for( int eval_i = 1; eval_i <= num_evals; ++eval_i ) {
00404     // Set the step constant and the weighting constant
00405     ScalarMag
00406       uh_i   = 0.0,
00407       wgt_i  = 0.0;
00408     switch(fd_method_type) {
00409       case DFDCT::FD_ORDER_ONE: {
00410         switch(eval_i) {
00411           case 1:
00412             uh_i  = +0.0;
00413             wgt_i = -1.0;
00414             break;
00415           case 2:
00416             uh_i  = +1.0;
00417             wgt_i = +1.0;
00418             break;
00419         }
00420         break;
00421       }
00422       case DFDCT::FD_ORDER_TWO: {
00423         switch(eval_i) {
00424           case 1:
00425             uh_i  = +0.0;
00426             wgt_i = -3.0;
00427             break;
00428           case 2:
00429             uh_i  = +1.0;
00430             wgt_i = +4.0;
00431             break;
00432           case 3:
00433             uh_i  = +2.0;
00434             wgt_i = -1.0;
00435             break;
00436         }
00437         break;
00438       }
00439       case DFDCT::FD_ORDER_TWO_CENTRAL: {
00440         switch(eval_i) {
00441           case 1:
00442             uh_i  = -1.0;
00443             wgt_i = -1.0;
00444             break;
00445           case 2:
00446             uh_i  = +1.0;
00447             wgt_i = +1.0;
00448             break;
00449         }
00450         break;
00451       }
00452       case DFDCT::FD_ORDER_FOUR: {
00453         switch(eval_i) {
00454           case 1:
00455             uh_i  = +0.0;
00456             wgt_i = -25.0;
00457             break;
00458           case 2:
00459             uh_i  = +1.0;
00460             wgt_i = +48.0;
00461             break;
00462           case 3:
00463             uh_i  = +2.0;
00464             wgt_i = -36.0;
00465             break;
00466           case 4:
00467             uh_i  = +3.0;
00468             wgt_i = +16.0;
00469             break;
00470           case 5:
00471             uh_i  = +4.0;
00472             wgt_i = -3.0;
00473             break;
00474         }
00475         break;
00476       }
00477       case DFDCT::FD_ORDER_FOUR_CENTRAL: {
00478         switch(eval_i) {
00479           case 1:
00480             uh_i  = -2.0;
00481             wgt_i = +1.0;
00482             break;
00483           case 2:
00484             uh_i  = -1.0;
00485             wgt_i = -8.0;
00486             break;
00487           case 3:
00488             uh_i  = +1.0;
00489             wgt_i = +8.0;
00490             break;
00491           case 4:
00492             uh_i  = +2.0;
00493             wgt_i = -1.0;
00494             break;
00495         }
00496         break;
00497       }
00498       case DFDCT::FD_ORDER_TWO_AUTO:
00499       case DFDCT::FD_ORDER_FOUR_AUTO:
00500         break; // Okay
00501       default:
00502         TEST_FOR_EXCEPT(true);
00503     }
00504 
00505     if(out.get() && trace)
00506       *out << "\neval_i="<<eval_i<<", uh_i="<<uh_i<<", wgt_i="<<wgt_i<<"\n";
00507     Teuchos::OSTab tab(out);
00508     
00509     // Compute the weighted term and add it to the sum
00510     if(uh_i == ST::zero()) {
00511       MEB::OutArgs<Scalar> bfuncall;
00512       if(!all_funcs_at_base_computed) {
00513         // Compute missing functions at the base point
00514         bfuncall = model.createOutArgs();
00515         if( pfunc.supports(MEB::OUT_ARG_f) && pfunc.get_f().get() && !bfunc.get_f().get() ) {
00516           bfuncall.set_f(createMember(model.get_f_space()));
00517         }
00518         for( int j = 0; j < Ng; ++j ) {
00519           if( pfunc.get_g(j).get() && !bfunc.get_g(j).get() ) {
00520             bfuncall.set_g(j,createMember(model.get_g_space(j)));
00521           }
00522         }
00523         model.evalModel(bp,bfuncall);
00524         bfuncall.setArgs(bfunc);
00525       }
00526       else {
00527         bfuncall = bfunc;
00528       }
00529       // Use functions at the base point
00530       if(out.get() && trace)
00531         *out << "\nSetting variations = wgt_i * basePoint ...\n";
00532       VectorPtr f;
00533       if( pfunc.supports(MEB::OUT_ARG_f) && (f=var.get_f()).get() ) {
00534         V_StV(&*f,wgt_i,*bfuncall.get_f());
00535       }
00536       for( int j = 0; j < Ng; ++j ) {
00537         VectorPtr g_j;
00538         if( (g_j=var.get_g(j)).get() ) {
00539           V_StV(&*g_j,wgt_i,*bfuncall.get_g(j));
00540         }
00541       }
00542     }
00543     else {
00544       if(out.get() && trace)
00545         *out << "\nSetting perturbedPoint = basePoint + uh_i*uh*direction ...\n";
00546       // z = zo + uh_i*uh*v
00547       if(1) {
00548         if( dir.supports(MEB::IN_ARG_x) && dir.get_x().get() )
00549           V_StVpV(&*per_x,Scalar(uh_i*uh),*dir.get_x(),*bp.get_x());
00550         if( dir.supports(MEB::IN_ARG_x_dot) && dir.get_x_dot().get() )
00551           V_StVpV(&*per_x_dot,Scalar(uh_i*uh),*dir.get_x_dot(),*bp.get_x_dot());
00552         for( int l = 0; l < Np; ++l ) {
00553           if( dir.get_p(l).get() )
00554             V_StVpV(&*per_p[l],Scalar(uh_i*uh),*dir.get_p(l),*bp.get_p(l));
00555         }
00556       }
00557       if(out.get() && trace)
00558         *out << "\nperturbedPoint=\n" << describe(pp,verbLevel);
00559       // Compute perturbed function values h(zo+uh_i*uh)
00560       if(out.get() && trace)
00561         *out << "\nCompute perturnedFunctions at perturbedPoint...\n";
00562       model.evalModel(pp,pfunc);
00563       if(out.get() && trace)
00564         *out << "\nperturnedFunctions=\n" << describe(pfunc,verbLevel);
00565       // Sum perturbed function values into total variation
00566       if(1) {
00567         // var_h += wgt_i*perturbed_h
00568         if(out.get() && trace)
00569           *out << "\nComputing variations += wgt_i*perturbedfunctions ...\n";
00570         VectorPtr f;
00571         if( pfunc.supports(MEB::OUT_ARG_f) && (f=pfunc.get_f()).get() )
00572           Vp_StV(&*var.get_f(),wgt_i,*f);
00573         for( int j = 0; j < Ng; ++j ) {
00574           VectorPtr g_j;
00575           if( (g_j=pfunc.get_g(j)).get() )
00576             Vp_StV(&*var.get_g(j),wgt_i,*g_j);
00577         }
00578       }
00579     }
00580     if(out.get() && trace)
00581       *out << "\nvariations=\n" << describe(var,verbLevel);
00582   }
00583   //
00584   // Multiply by the scaling factor!
00585   //
00586   
00587   if(1) {
00588     // var_h *= 1.0/(dwgt*uh)
00589     const Scalar alpha = ST::one()/(dwgt*uh);
00590     if(out.get() && trace)
00591       *out
00592         << "\nComputing variations *= (1.0)/(dwgt*uh),"
00593         << " where (1.0)/(dwgt*uh) = (1.0)/("<<dwgt<<"*"<<uh<<") = "<<alpha<<" ...\n";
00594     VectorPtr f;
00595     if( var.supports(MEB::OUT_ARG_f) && (f=var.get_f()).get() )
00596       Vt_S(&*f,alpha);
00597     for( int j = 0; j < Ng; ++j ) {
00598       VectorPtr g_j;
00599       if( (g_j=var.get_g(j)).get() )
00600         Vt_S(&*g_j,alpha);
00601     }
00602     if(out.get() && trace)
00603       *out << "\nFinal variations=\n" << describe(var,verbLevel);
00604   }
00605   
00606   if(out.get())
00607     *out << std::setprecision(p_saved);
00608 
00609   if(out.get() && trace)
00610     *out << "\nLeaving DirectionalFiniteDiffCalculator<Scalar>::calcVariations(...)\n";
00611   
00612 }
00613 
00614 template<class Scalar>
00615 void DirectionalFiniteDiffCalculator<Scalar>::calcDerivatives(
00616   const ModelEvaluator<Scalar>                 &model
00617   ,const ModelEvaluatorBase::InArgs<Scalar>    &bp      // basePoint
00618   ,const ModelEvaluatorBase::OutArgs<Scalar>   &bfunc   // baseFunctionValues
00619   ,const ModelEvaluatorBase::OutArgs<Scalar>   &deriv   // derivatives
00620   ) const
00621 {
00622   typedef Teuchos::ScalarTraits<Scalar> ST;
00623   typedef ModelEvaluatorBase MEB;
00624   typedef Teuchos::RefCountPtr<VectorBase<Scalar> > VectorPtr;
00625   typedef Teuchos::RefCountPtr<MultiVectorBase<Scalar> > MultiVectorPtr;
00626 
00627   Teuchos::RefCountPtr<Teuchos::FancyOStream> out = this->getOStream();
00628   Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
00629   const bool trace = (static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_MEDIUM));
00630   Teuchos::OSTab tab(out);
00631 
00632   if(out.get() && trace)
00633     *out << "\nEntering DirectionalFiniteDiffCalculator<Scalar>::calcDerivatives(...)\n";
00634 
00635   if(out.get() && trace)
00636     *out
00637       << "\nbasePoint=\n" << describe(bp,verbLevel)
00638       << "\nbaseFunctionValues=\n" << describe(bfunc,verbLevel)
00639       << "\nderivatives=\n" << describe(deriv,Teuchos::VERB_LOW);
00640   
00641   //
00642   // We will only compute finite differences w.r.t. p(l) for now
00643   //
00644   const int Np = bp.Np(), Ng = bfunc.Ng();
00645   MEB::InArgs<Scalar> dir = model.createInArgs();
00646   MEB::OutArgs<Scalar> var = model.createOutArgs();
00647   MultiVectorPtr DfDp_l;
00648   std::vector<MEB::DerivativeMultiVector<Scalar> > DgDp_l(Ng);
00649   std::vector<VectorPtr> var_g(Ng);
00650   for( int l = 0; l < Np; ++l ) {
00651     if(out.get() && trace)
00652       *out << "\nComputing deriatives for parameter subvector p("<<l<<") ...\n";
00653     Teuchos::OSTab tab(out);
00654     //
00655     // Load up OutArgs var object of derivative variations to compute
00656     //
00657     bool hasDerivObject = false;
00658     // DfDp(l)
00659     if( deriv.supports(MEB::OUT_ARG_DfDp,l).none()==false
00660         && deriv.get_DfDp(l).isEmpty()==false )
00661     {
00662       hasDerivObject = true;
00663       std::ostringstream name; name << "DfDp("<<l<<")";
00664       DfDp_l = get_mv(deriv.get_DfDp(l),name.str(),MEB::DERIV_MV_BY_COL);
00665     }
00666     else {
00667       DfDp_l = Teuchos::null;
00668     }
00669     // DgDp(j=1...Ng,l)
00670     for( int j = 0; j < Ng; ++j ) {
00671       if(
00672         deriv.supports(MEB::OUT_ARG_DgDp,j,l).none()==false
00673         && deriv.get_DgDp(j,l).isEmpty()==false
00674         )
00675       {
00676         hasDerivObject = true;
00677         std::ostringstream name; name << "DgDp("<<j<<","<<l<<")";
00678         DgDp_l[j] = get_dmv(deriv.get_DgDp(j,l),name.str());
00679         if( DgDp_l[j].getMultiVector().get()
00680             && DgDp_l[j].getOrientation()==MEB::DERIV_TRANS_MV_BY_ROW
00681             && !var_g[j].get() )
00682         {
00683           // Need a temporary vector for the variation for g(j)
00684           var_g[j] = createMember(model.get_g_space(j));
00685         }
00686       }
00687       else{
00688         DgDp_l[j] = MEB::DerivativeMultiVector<Scalar>();
00689         var_g[j] = Teuchos::null;
00690       }
00691     }
00692     //
00693     // Compute derivative variations by finite differences
00694     //
00695     if(hasDerivObject) {
00696       VectorPtr e_i = createMember(model.get_p_space(l));
00697       dir.set_p(l,e_i);
00698       assign(&*e_i,ST::zero());
00699       const int np_l = e_i->space()->dim();
00700       for( int i = 0 ; i < np_l; ++ i ) {
00701         if(out.get() && trace)
00702           *out << "\nComputing derivatives for single variable p("<<l<<")("<<i<<") ...\n";
00703         Teuchos::OSTab tab(out);
00704         if(DfDp_l.get()) var.set_f(DfDp_l->col(i)); // Compute in place!
00705         for(int j = 0; j < Ng; ++j) {
00706           MultiVectorPtr DgDp_j_l;
00707           if( (DgDp_j_l=DgDp_l[j].getMultiVector()).get() ) {
00708             if(DgDp_l[j].getOrientation()==MEB::DERIV_TRANS_MV_BY_ROW) {
00709               var.set_g(j,var_g[j]); // Computes d(g(j))/d(p(l)(i))
00710             }
00711             else {
00712               TEST_FOR_EXCEPTION(
00713                 true, std::logic_error
00714                 ,"Error, do not support the nontransposed version of DgDp yet!"
00715                 );
00716             }
00717           }
00718         }
00719         set_ele(i,ST::one(),&*e_i);
00720         this->calcVariations(
00721           model,bp,dir,bfunc,var
00722           );
00723         set_ele(i,ST::zero(),&*e_i);
00724         if(DfDp_l.get()) var.set_f(Teuchos::null);
00725         for(int j = 0; j < Ng; ++j) {
00726           MultiVectorPtr DgDp_j_l;
00727           if( (DgDp_j_l=DgDp_l[j].getMultiVector()).get() ) {
00728             if(DgDp_l[j].getOrientation()==MEB::DERIV_TRANS_MV_BY_ROW) {
00729               ConstDetachedVectorView<Scalar> _var_g_j(var_g[j]); //d(g(j))/d(p(l)(i))
00730               DetachedMultiVectorView<Scalar> _DgDp_j_l_i(*DgDp_j_l,Range1D(i,i),Range1D());
00731               const int ng = _var_g_j.subDim();
00732               for( int k = 0; k < ng; ++k )
00733                 _DgDp_j_l_i(0,k) = _var_g_j(k);
00734             }
00735             else {
00736               TEST_FOR_EXCEPTION(
00737                 true, std::logic_error
00738                 ,"Error, do not support the nontransposed version of DgDp yet!"
00739                 );
00740             }
00741           }
00742         }
00743       }
00744       dir.set_p(l,Teuchos::null);
00745     }
00746   }
00747     
00748   if(out.get() && trace)
00749     *out
00750       << "\nderivatives=\n" << describe(deriv,verbLevel);
00751   
00752   if(out.get() && trace)
00753     *out << "\nEntering DirectionalFiniteDiffCalculator<Scalar>::calcDerivatives(...)\n";
00754   
00755 }
00756 
00757 } // end namespace Thyra
00758 
00759 #endif  // THYRA_DIRECTIONAL_FINITE_DIFF_CALCULATOR_HPP

Generated on Thu Sep 18 12:32:48 2008 for Thyra Nonlinear Model Evaluator Support by doxygen 1.3.9.1