Sacado Package Browser (Single Doxygen Collection) Version of the Day
Sacado_ELRFad_GeneralFadImp.hpp
Go to the documentation of this file.
00001 // $Id$ 
00002 // $Source$ 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 //                           Sacado Package
00007 //                 Copyright (2006) Sandia Corporation
00008 // 
00009 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
00011 // 
00012 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //  
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
00027 // (etphipp@sandia.gov).
00028 // 
00029 // ***********************************************************************
00030 //
00031 // The forward-mode AD classes in Sacado are a derivative work of the
00032 // expression template classes in the Fad package by Nicolas Di Cesare.  
00033 // The following banner is included in the original Fad source code:
00034 //
00035 // ************ DO NOT REMOVE THIS BANNER ****************
00036 //
00037 //  Nicolas Di Cesare <Nicolas.Dicesare@ann.jussieu.fr>
00038 //  http://www.ann.jussieu.fr/~dicesare
00039 //
00040 //            CEMRACS 98 : C++ courses, 
00041 //         templates : new C++ techniques 
00042 //            for scientific computing 
00043 // 
00044 //********************************************************
00045 //
00046 //  A short implementation ( not all operators and 
00047 //  functions are overloaded ) of 1st order Automatic
00048 //  Differentiation in forward mode (FAD) using
00049 //  EXPRESSION TEMPLATES.
00050 //
00051 //********************************************************
00052 // @HEADER
00053 
00054 #include "Sacado_ConfigDefs.h"
00055 #include "Sacado_mpl_range_c.hpp"
00056 #include "Sacado_mpl_for_each.hpp"
00057 
00058 // template <typename T, typename Storage> 
00059 // template <typename S> 
00060 // inline Sacado::ELRFad::GeneralFad<T,Storage>::GeneralFad(const Expr<S>& x) :
00061 //   Storage(T(0.))
00062 // {
00063 //   int sz = x.size();
00064 
00065 //   if (sz != this->size()) 
00066 //     this->resize(sz);
00067 
00068 //   if (sz) {
00069 
00070 //     // Number of arguments
00071 //     const int N = Expr<S>::num_args;
00072 
00073 //     // Array to store partials
00074 //     T partials[N];
00075 
00076 //     // Array to store tangents
00077 //     T dots[N];
00078 
00079 //     // Compute partials
00080 //     x.computePartials(T(1.), partials);
00081   
00082 //     // Compute each tangent direction
00083 //     for(int i=0; i<sz; ++i) {
00084 //       T t = T(0.);
00085 //       x.getTangents(i,dots);
00086 //       for (int j=0; j<N; j++) {
00087 //  t += partials[j]*dots[j];
00088 //       }
00089 //       this->fastAccessDx(i) = t;
00090 //     }
00091 
00092 //   }
00093   
00094 //   // Compute value
00095 //   this->val() = x.val();
00096 // }
00097 
00098 template <typename T, typename Storage> 
00099 template <typename S> 
00100 inline Sacado::ELRFad::GeneralFad<T,Storage>::GeneralFad(const Expr<S>& x) :
00101   Storage(T(0.)),
00102   update_val_(x.updateValue())
00103 {
00104   int sz = x.size();
00105 
00106   if (sz != this->size()) 
00107     this->resize(sz);
00108 
00109   if (sz) {
00110 
00111     if (Expr<S>::is_linear) {
00112       if (x.hasFastAccess())
00113         for(int i=0; i<sz; ++i)
00114           this->fastAccessDx(i) = x.fastAccessDx(i);
00115       else
00116         for(int i=0; i<sz; ++i)
00117           this->fastAccessDx(i) = x.dx(i);
00118     }
00119     else {
00120 
00121       // Number of arguments
00122       const int N = Expr<S>::num_args;
00123 
00124       if (x.hasFastAccess()) {
00125   // Compute partials
00126   FastLocalAccumOp< Expr<S> > op(x);
00127   
00128   // Compute each tangent direction
00129   for(op.i=0; op.i<sz; ++op.i) {
00130     op.t = T(0.);
00131 
00132     // Automatically unrolled loop that computes
00133     // for (int j=0; j<N; j++)
00134     //   op.t += op.partials[j] * x.getTangent<j>(i);
00135     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00136     
00137     this->fastAccessDx(op.i) = op.t;
00138   }
00139       }
00140       else {
00141   // Compute partials
00142   SlowLocalAccumOp< Expr<S> > op(x);
00143   
00144   // Compute each tangent direction
00145   for(op.i=0; op.i<sz; ++op.i) {
00146     op.t = T(0.);
00147 
00148     // Automatically unrolled loop that computes
00149     // for (int j=0; j<N; j++)
00150     //   op.t += op.partials[j] * x.getTangent<j>(i);
00151     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00152     
00153     this->fastAccessDx(op.i) = op.t;
00154   }
00155       }
00156  
00157     }
00158 
00159   }
00160   
00161   // Compute value
00162   if (update_val_)
00163     this->val() = x.val();
00164 }
00165 
00166 template <typename T, typename Storage> 
00167 inline void 
00168 Sacado::ELRFad::GeneralFad<T,Storage>::diff(const int ith, const int n) 
00169 { 
00170   if (this->size() != n) 
00171     this->resize(n);
00172 
00173   this->zero();
00174   this->fastAccessDx(ith) = T(1.);
00175 
00176 }
00177 
00178 template <typename T, typename Storage> 
00179 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00180 Sacado::ELRFad::GeneralFad<T,Storage>::operator=(const T& v) 
00181 {
00182   this->val() = v;
00183 
00184   if (this->size()) {
00185     this->zero();    // We need to zero out the array for future resizes
00186     this->resize(0);
00187   }
00188 
00189   return *this;
00190 }
00191 
00192 template <typename T, typename Storage> 
00193 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00194 Sacado::ELRFad::GeneralFad<T,Storage>::operator=(
00195              const Sacado::ELRFad::GeneralFad<T,Storage>& x) 
00196 {
00197   // Copy value & dx_
00198   Storage::operator=(x);
00199   update_val_ = x.update_val_;
00200   
00201   return *this;
00202 }
00203 
00204 template <typename T, typename Storage> 
00205 template <typename S> 
00206 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00207 Sacado::ELRFad::GeneralFad<T,Storage>::operator=(const Expr<S>& x) 
00208 {
00209   int sz = x.size();
00210 
00211   if (sz != this->size()) 
00212     this->resize(sz);
00213 
00214   if (sz) {
00215  
00216     if (Expr<S>::is_linear) {
00217       if (x.hasFastAccess())
00218         for(int i=0; i<sz; ++i)
00219           this->fastAccessDx(i) = x.fastAccessDx(i);
00220       else
00221         for(int i=0; i<sz; ++i)
00222           this->fastAccessDx(i) = x.dx(i);
00223     }
00224     else {
00225 
00226       // Number of arguments
00227       const int N = Expr<S>::num_args;
00228 
00229       if (x.hasFastAccess()) {
00230   // Compute partials
00231   FastLocalAccumOp< Expr<S> > op(x);
00232   
00233   // Compute each tangent direction
00234   for(op.i=0; op.i<sz; ++op.i) {
00235     op.t = T(0.);
00236 
00237     // Automatically unrolled loop that computes
00238     // for (int j=0; j<N; j++)
00239     //   op.t += op.partials[j] * x.getTangent<j>(i);
00240     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00241     
00242     this->fastAccessDx(op.i) = op.t;
00243   }
00244       }
00245       else {
00246   // Compute partials
00247   SlowLocalAccumOp< Expr<S> > op(x);
00248   
00249   // Compute each tangent direction
00250   for(op.i=0; op.i<sz; ++op.i) {
00251     op.t = T(0.);
00252 
00253     // Automatically unrolled loop that computes
00254     // for (int j=0; j<N; j++)
00255     //   op.t += op.partials[j] * x.getTangent<j>(i);
00256     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00257     
00258     this->fastAccessDx(op.i) = op.t;
00259   }
00260       }
00261     }
00262   }
00263   
00264   // Compute value
00265   update_val_ = x.updateValue();
00266   if (update_val_)
00267     this->val() = x.val();
00268   
00269   return *this;
00270 }
00271 
00272 template <typename T, typename Storage> 
00273 inline  Sacado::ELRFad::GeneralFad<T,Storage>& 
00274 Sacado::ELRFad::GeneralFad<T,Storage>::operator += (const T& v)
00275 {
00276   if (update_val_)
00277     this->val() += v;
00278 
00279   return *this;
00280 }
00281 
00282 template <typename T, typename Storage> 
00283 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00284 Sacado::ELRFad::GeneralFad<T,Storage>::operator -= (const T& v)
00285 {
00286   if (update_val_)
00287     this->val() -= v;
00288 
00289   return *this;
00290 }
00291 
00292 template <typename T, typename Storage> 
00293 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00294 Sacado::ELRFad::GeneralFad<T,Storage>::operator *= (const T& v)
00295 {
00296   int sz = this->size();
00297 
00298   if (update_val_)
00299     this->val() *= v;
00300   for (int i=0; i<sz; ++i)
00301     this->fastAccessDx(i) *= v;
00302 
00303   return *this;
00304 }
00305 
00306 template <typename T, typename Storage> 
00307 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00308 Sacado::ELRFad::GeneralFad<T,Storage>::operator /= (const T& v)
00309 {
00310   int sz = this->size();
00311 
00312   if (update_val_)
00313     this->val() /= v;
00314   for (int i=0; i<sz; ++i)
00315     this->fastAccessDx(i) /= v;
00316 
00317   return *this;
00318 }
00319 
00320 template <typename T, typename Storage> 
00321 inline  Sacado::ELRFad::GeneralFad<T,Storage>& 
00322 Sacado::ELRFad::GeneralFad<T,Storage>::
00323 operator += (const typename Sacado::dummy<value_type,scalar_type>::type& v)
00324 {
00325   if (update_val_)
00326     this->val() += v;
00327 
00328   return *this;
00329 }
00330 
00331 template <typename T, typename Storage> 
00332 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00333 Sacado::ELRFad::GeneralFad<T,Storage>::
00334 operator -= (const typename Sacado::dummy<value_type,scalar_type>::type& v)
00335 {
00336   if (update_val_)
00337     this->val() -= v;
00338 
00339   return *this;
00340 }
00341 
00342 template <typename T, typename Storage> 
00343 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00344 Sacado::ELRFad::GeneralFad<T,Storage>::
00345 operator *= (const typename Sacado::dummy<value_type,scalar_type>::type& v)
00346 {
00347   int sz = this->size();
00348 
00349   if (update_val_)
00350     this->val() *= v;
00351   for (int i=0; i<sz; ++i)
00352     this->fastAccessDx(i) *= v;
00353 
00354   return *this;
00355 }
00356 
00357 template <typename T, typename Storage> 
00358 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00359 Sacado::ELRFad::GeneralFad<T,Storage>::
00360 operator /= (const typename Sacado::dummy<value_type,scalar_type>::type& v)
00361 {
00362   int sz = this->size();
00363 
00364   if (update_val_)
00365     this->val() /= v;
00366   for (int i=0; i<sz; ++i)
00367     this->fastAccessDx(i) /= v;
00368 
00369   return *this;
00370 }
00371 
00372 template <typename T, typename Storage> 
00373 template <typename S> 
00374 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00375 Sacado::ELRFad::GeneralFad<T,Storage>::operator += (
00376                const Sacado::ELRFad::Expr<S>& x)
00377 {
00378   int xsz = x.size(), sz = this->size();
00379 
00380 #ifdef SACADO_DEBUG
00381   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00382     throw "Fad Error:  Attempt to assign with incompatible sizes";
00383 #endif
00384 
00385   if (Expr<S>::is_linear) {
00386     if (xsz) {
00387       if (sz) {
00388   if (x.hasFastAccess())
00389     for (int i=0; i<sz; ++i)
00390       this->fastAccessDx(i) += x.fastAccessDx(i);
00391   else
00392     for (int i=0; i<sz; ++i)
00393       this->fastAccessDx(i) += x.dx(i);
00394       }
00395       else {
00396   this->resize(xsz);
00397   if (x.hasFastAccess())
00398     for (int i=0; i<xsz; ++i)
00399       this->fastAccessDx(i) = x.fastAccessDx(i);
00400   else
00401     for (int i=0; i<xsz; ++i)
00402       this->fastAccessDx(i) = x.dx(i);
00403       }
00404     }
00405   }
00406   else {
00407 
00408     // Number of arguments
00409     const int N = Expr<S>::num_args;
00410     
00411     if (xsz) {
00412 
00413       if (sz != xsz) {
00414   this->resize(xsz);
00415   this->zero();
00416       }
00417 
00418       if (x.hasFastAccess()) {
00419   // Compute partials
00420   FastLocalAccumOp< Expr<S> > op(x);
00421   
00422   // Compute each tangent direction
00423   for(op.i=0; op.i<xsz; ++op.i) {
00424     op.t = T(0.);
00425     
00426     // Automatically unrolled loop that computes
00427     // for (int j=0; j<N; j++)
00428     //   op.t += op.partials[j] * x.getTangent<j>(i);
00429     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00430     
00431     this->fastAccessDx(op.i) += op.t;
00432   }
00433       }
00434       else {
00435   // Compute partials
00436   SlowLocalAccumOp< Expr<S> > op(x);
00437   
00438   // Compute each tangent direction
00439   for(op.i=0; op.i<xsz; ++op.i) {
00440     op.t = T(0.);
00441 
00442     // Automatically unrolled loop that computes
00443     // for (int j=0; j<N; j++)
00444     //   op.t += op.partials[j] * x.getTangent<j>(i);
00445     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00446     
00447     this->fastAccessDx(op.i) += op.t;
00448   }
00449       }
00450       
00451     }
00452 
00453   }
00454   
00455   // Compute value
00456   update_val_ = x.updateValue();
00457   if (update_val_)
00458     this->val() += x.val();
00459 
00460   return *this;
00461 }
00462 
00463 template <typename T, typename Storage> 
00464 template <typename S> 
00465 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00466 Sacado::ELRFad::GeneralFad<T,Storage>::operator -= (
00467                 const Sacado::ELRFad::Expr<S>& x)
00468 {
00469   int xsz = x.size(), sz = this->size();
00470 
00471 #ifdef SACADO_DEBUG
00472   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00473     throw "Fad Error:  Attempt to assign with incompatible sizes";
00474 #endif
00475 
00476   if (Expr<S>::is_linear) {
00477     if (xsz) {
00478       if (sz) {
00479   if (x.hasFastAccess())
00480     for(int i=0; i<sz; ++i)
00481       this->fastAccessDx(i) -= x.fastAccessDx(i);
00482   else
00483     for (int i=0; i<sz; ++i)
00484       this->fastAccessDx(i) -= x.dx(i);
00485       }
00486       else {
00487   this->resize(xsz);
00488   if (x.hasFastAccess())
00489     for(int i=0; i<xsz; ++i)
00490       this->fastAccessDx(i) = -x.fastAccessDx(i);
00491   else
00492     for (int i=0; i<xsz; ++i)
00493       this->fastAccessDx(i) = -x.dx(i);
00494       }
00495     }
00496   }
00497   else {
00498     // Number of arguments
00499     const int N = Expr<S>::num_args;
00500     
00501     if (xsz) {
00502       
00503       if (sz != xsz) {
00504   this->resize(xsz);
00505   this->zero();
00506       }
00507 
00508       if (x.hasFastAccess()) {
00509   // Compute partials
00510   FastLocalAccumOp< Expr<S> > op(x);
00511   
00512   // Compute each tangent direction
00513   for(op.i=0; op.i<xsz; ++op.i) {
00514     op.t = T(0.);
00515     
00516     // Automatically unrolled loop that computes
00517     // for (int j=0; j<N; j++)
00518     //   op.t += op.partials[j] * x.getTangent<j>(i);
00519     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00520     
00521     this->fastAccessDx(op.i) -= op.t;
00522   }
00523       }
00524       else {
00525   // Compute partials
00526   SlowLocalAccumOp< Expr<S> > op(x);
00527   
00528   // Compute each tangent direction
00529   for(op.i=0; op.i<xsz; ++op.i) {
00530     op.t = T(0.);
00531 
00532     // Automatically unrolled loop that computes
00533     // for (int j=0; j<N; j++)
00534     //   op.t += op.partials[j] * x.getTangent<j>(i);
00535     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00536     
00537     this->fastAccessDx(op.i) -= op.t;
00538   } 
00539       }
00540     }
00541 
00542   }
00543 
00544   update_val_ = x.updateValue();
00545   if (update_val_)
00546     this->val() -= x.val();
00547 
00548   return *this;
00549 }
00550 
00551 template <typename T, typename Storage> 
00552 template <typename S> 
00553 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00554 Sacado::ELRFad::GeneralFad<T,Storage>::operator *= (
00555                const Sacado::ELRFad::Expr<S>& x)
00556 {
00557   int xsz = x.size(), sz = this->size();
00558   T xval = x.val();
00559   T v = this->val();
00560 
00561 #ifdef SACADO_DEBUG
00562   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00563     throw "Fad Error:  Attempt to assign with incompatible sizes";
00564 #endif
00565 
00566   if (Expr<S>::is_linear) {
00567     if (xsz) {
00568       if (sz) {
00569   if (x.hasFastAccess())
00570     for(int i=0; i<sz; ++i)
00571       this->fastAccessDx(i) = v*x.fastAccessDx(i) + this->fastAccessDx(i)*xval;
00572   else
00573     for (int i=0; i<sz; ++i)
00574       this->fastAccessDx(i) = v*x.dx(i) + this->fastAccessDx(i)*xval;
00575       }
00576       else {
00577   this->resize(xsz);
00578   if (x.hasFastAccess())
00579     for(int i=0; i<xsz; ++i)
00580       this->fastAccessDx(i) = v*x.fastAccessDx(i);
00581   else
00582     for (int i=0; i<xsz; ++i)
00583       this->fastAccessDx(i) = v*x.dx(i);
00584       }
00585     }
00586     else {
00587       if (sz) {
00588   for (int i=0; i<sz; ++i)
00589     this->fastAccessDx(i) *= xval;
00590       }
00591     }
00592   }
00593   else {
00594 
00595     // Number of arguments
00596     const int N = Expr<S>::num_args;
00597     
00598     if (xsz) {
00599       
00600       if (sz) {
00601 
00602   if (x.hasFastAccess()) {
00603     // Compute partials
00604     FastLocalAccumOp< Expr<S> > op(x);
00605     
00606     // Compute each tangent direction
00607     for(op.i=0; op.i<xsz; ++op.i) {
00608       op.t = T(0.);
00609       
00610       // Automatically unrolled loop that computes
00611       // for (int j=0; j<N; j++)
00612       //   op.t += op.partials[j] * x.getTangent<j>(i);
00613       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00614     
00615       this->fastAccessDx(op.i) = 
00616         v * op.t + this->fastAccessDx(op.i) * xval;
00617     }
00618   }
00619   else {
00620     // Compute partials
00621     SlowLocalAccumOp< Expr<S> > op(x);
00622     
00623     // Compute each tangent direction
00624     for(op.i=0; op.i<xsz; ++op.i) {
00625       op.t = T(0.);
00626       
00627       // Automatically unrolled loop that computes
00628       // for (int j=0; j<N; j++)
00629       //   op.t += op.partials[j] * x.getTangent<j>(i);
00630       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00631     
00632       this->fastAccessDx(op.i) = 
00633         v * op.t + this->fastAccessDx(op.i) * xval;
00634     }
00635   }
00636 
00637       }
00638       
00639       else {
00640   
00641   this->resize(xsz);
00642   
00643   if (x.hasFastAccess()) {
00644     // Compute partials
00645     FastLocalAccumOp< Expr<S> > op(x);
00646     
00647     // Compute each tangent direction
00648     for(op.i=0; op.i<xsz; ++op.i) {
00649       op.t = T(0.);
00650       
00651       // Automatically unrolled loop that computes
00652       // for (int j=0; j<N; j++)
00653       //   op.t += op.partials[j] * x.getTangent<j>(i);
00654       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00655     
00656       this->fastAccessDx(op.i) = v * op.t;
00657     }
00658   }
00659   else {
00660     // Compute partials
00661     SlowLocalAccumOp< Expr<S> > op(x);
00662     
00663     // Compute each tangent direction
00664     for(op.i=0; op.i<xsz; ++op.i) {
00665       op.t = T(0.);
00666       
00667       // Automatically unrolled loop that computes
00668       // for (int j=0; j<N; j++)
00669       //   op.t += op.partials[j] * x.getTangent<j>(i);
00670       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00671     
00672       this->fastAccessDx(op.i) = v * op.t;
00673     }
00674   }
00675   
00676       }
00677       
00678     }
00679     
00680     else {
00681       
00682       if (sz) {
00683   for (int i=0; i<sz; ++i)
00684     this->fastAccessDx(i) *= xval;
00685       }
00686       
00687     }
00688 
00689   }
00690 
00691   update_val_ = x.updateValue();
00692   if (update_val_)
00693     this->val() *= xval;
00694 
00695   return *this;
00696 }
00697 
00698 template <typename T, typename Storage>
00699 template <typename S> 
00700 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00701 Sacado::ELRFad::GeneralFad<T,Storage>::operator /= (
00702                const Sacado::ELRFad::Expr<S>& x)
00703 {
00704   int xsz = x.size(), sz = this->size();
00705   T xval = x.val();
00706   T v = this->val();
00707 
00708 #ifdef SACADO_DEBUG
00709   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00710     throw "Fad Error:  Attempt to assign with incompatible sizes";
00711 #endif
00712 
00713   if (Expr<S>::is_linear) {
00714     if (xsz) {
00715       if (sz) {
00716   if (x.hasFastAccess())
00717     for(int i=0; i<sz; ++i)
00718       this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.fastAccessDx(i) )/ (xval*xval);
00719   else
00720     for (int i=0; i<sz; ++i)
00721       this->fastAccessDx(i) = ( this->fastAccessDx(i)*xval - v*x.dx(i) )/ (xval*xval);
00722       }
00723       else {
00724   this->resize(xsz);
00725   if (x.hasFastAccess())
00726     for(int i=0; i<xsz; ++i)
00727       this->fastAccessDx(i) = - v*x.fastAccessDx(i) / (xval*xval);
00728   else
00729     for (int i=0; i<xsz; ++i)
00730       this->fastAccessDx(i) = -v*x.dx(i) / (xval*xval);
00731       }
00732     }
00733     else {
00734       if (sz) {
00735   for (int i=0; i<sz; ++i)
00736     this->fastAccessDx(i) /= xval;
00737       }
00738     }
00739   }
00740   else {
00741     // Number of arguments
00742     const int N = Expr<S>::num_args;
00743     
00744     if (xsz) {
00745       
00746       T xval2 = xval*xval;
00747       
00748       if (sz) {
00749 
00750   if (x.hasFastAccess()) {
00751     // Compute partials
00752     FastLocalAccumOp< Expr<S> > op(x);
00753     
00754     // Compute each tangent direction
00755     for(op.i=0; op.i<xsz; ++op.i) {
00756       op.t = T(0.);
00757       
00758       // Automatically unrolled loop that computes
00759       // for (int j=0; j<N; j++)
00760       //   op.t += op.partials[j] * x.getTangent<j>(i);
00761       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00762     
00763       this->fastAccessDx(op.i) = 
00764         (this->fastAccessDx(op.i) * xval - v * op.t) / xval2;
00765     }
00766   }
00767   else {
00768     // Compute partials
00769     SlowLocalAccumOp< Expr<S> > op(x);
00770     
00771     // Compute each tangent direction
00772     for(op.i=0; op.i<xsz; ++op.i) {
00773       op.t = T(0.);
00774       
00775       // Automatically unrolled loop that computes
00776       // for (int j=0; j<N; j++)
00777       //   op.t += op.partials[j] * x.getTangent<j>(i);
00778       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00779     
00780       this->fastAccessDx(op.i) = 
00781         (this->fastAccessDx(op.i) * xval - v * op.t) / xval2;
00782     }
00783   }
00784   
00785       }
00786       
00787       else {
00788   
00789   this->resize(xsz);
00790   
00791   if (x.hasFastAccess()) {
00792     // Compute partials
00793     FastLocalAccumOp< Expr<S> > op(x);
00794     
00795     // Compute each tangent direction
00796     for(op.i=0; op.i<xsz; ++op.i) {
00797       op.t = T(0.);
00798       
00799       // Automatically unrolled loop that computes
00800       // for (int j=0; j<N; j++)
00801       //   op.t += op.partials[j] * x.getTangent<j>(i);
00802       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00803     
00804       this->fastAccessDx(op.i) = (-v * op.t) / xval2;
00805     }
00806   }
00807   else {
00808     // Compute partials
00809     SlowLocalAccumOp< Expr<S> > op(x);
00810     
00811     // Compute each tangent direction
00812     for(op.i=0; op.i<xsz; ++op.i) {
00813       op.t = T(0.);
00814       
00815       // Automatically unrolled loop that computes
00816       // for (int j=0; j<N; j++)
00817       //   op.t += op.partials[j] * x.getTangent<j>(i);
00818       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00819     
00820       this->fastAccessDx(op.i) = (-v * op.t) / xval2;
00821     }
00822   }
00823   
00824       }
00825       
00826     }
00827     
00828     else {
00829       
00830       if (sz) {
00831   for (int i=0; i<sz; ++i)
00832     this->fastAccessDx(i) /= xval;
00833       }
00834       
00835     }
00836 
00837   }
00838 
00839   update_val_ = x.updateValue();
00840   if (update_val_)
00841     this->val() /= xval;
00842 
00843   return *this;
00844 }
00845 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines