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     // Number of arguments
00112     const int N = Expr<S>::num_args;
00113 
00114     // Compute partials
00115     LocalAccumOp< Expr<S> > op(x);
00116   
00117     // Compute each tangent direction
00118     for(int i=0; i<sz; ++i) {
00119       op.t = T(0.);
00120       op.i = i;
00121 
00122       // Automatically unrolled loop that computes
00123       // for (int j=0; j<N; j++)
00124       //   op.t += op.partials[j] * x.getTangent<j>(i);
00125       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00126 
00127       this->fastAccessDx(i) = op.t;
00128     }
00129 
00130   }
00131   
00132   // Compute value
00133   if (update_val_)
00134     this->val() = x.val();
00135 }
00136 
00137 template <typename T, typename Storage> 
00138 inline void 
00139 Sacado::ELRFad::GeneralFad<T,Storage>::diff(const int ith, const int n) 
00140 { 
00141   if (this->size() != n) 
00142     this->resize(n);
00143 
00144   this->zero();
00145   this->fastAccessDx(ith) = T(1.);
00146 
00147 }
00148 
00149 template <typename T, typename Storage> 
00150 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00151 Sacado::ELRFad::GeneralFad<T,Storage>::operator=(const T& v) 
00152 {
00153   this->val() = v;
00154 
00155   if (this->size()) {
00156     this->zero();    // We need to zero out the array for future resizes
00157     this->resize(0);
00158   }
00159 
00160   return *this;
00161 }
00162 
00163 template <typename T, typename Storage> 
00164 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00165 Sacado::ELRFad::GeneralFad<T,Storage>::operator=(
00166              const Sacado::ELRFad::GeneralFad<T,Storage>& x) 
00167 {
00168   // Copy value & dx_
00169   Storage::operator=(x);
00170   update_val_ = x.update_val_;
00171   
00172   return *this;
00173 }
00174 
00175 template <typename T, typename Storage> 
00176 template <typename S> 
00177 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00178 Sacado::ELRFad::GeneralFad<T,Storage>::operator=(const Expr<S>& x) 
00179 {
00180   int sz = x.size();
00181 
00182   if (sz != this->size()) 
00183     this->resize(sz);
00184 
00185   if (sz) {
00186 
00187     // Number of arguments
00188     const int N = Expr<S>::num_args;
00189 
00190     // Compute partials
00191     LocalAccumOp< Expr<S> > op(x);
00192   
00193     // Compute each tangent direction
00194     for(int i=0; i<sz; ++i) {
00195       op.t = T(0.);
00196       op.i = i;
00197 
00198       // Automatically unrolled loop that computes
00199       // for (int j=0; j<N; j++)
00200       //   op.t += op.partials[j] * x.getTangent<j>(i);
00201       Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00202 
00203       this->fastAccessDx(i) = op.t;
00204     }
00205 
00206   }
00207   
00208   // Compute value
00209   update_val_ = x.updateValue();
00210   if (update_val_)
00211     this->val() = x.val();
00212   
00213   return *this;
00214 }
00215 
00216 template <typename T, typename Storage> 
00217 inline  Sacado::ELRFad::GeneralFad<T,Storage>& 
00218 Sacado::ELRFad::GeneralFad<T,Storage>::operator += (const T& v)
00219 {
00220   if (update_val_)
00221     this->val() += v;
00222 
00223   return *this;
00224 }
00225 
00226 template <typename T, typename Storage> 
00227 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00228 Sacado::ELRFad::GeneralFad<T,Storage>::operator -= (const T& v)
00229 {
00230   if (update_val_)
00231     this->val() -= v;
00232 
00233   return *this;
00234 }
00235 
00236 template <typename T, typename Storage> 
00237 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00238 Sacado::ELRFad::GeneralFad<T,Storage>::operator *= (const T& v)
00239 {
00240   int sz = this->size();
00241 
00242   if (update_val_)
00243     this->val() *= v;
00244   for (int i=0; i<sz; ++i)
00245     this->fastAccessDx(i) *= v;
00246 
00247   return *this;
00248 }
00249 
00250 template <typename T, typename Storage> 
00251 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00252 Sacado::ELRFad::GeneralFad<T,Storage>::operator /= (const T& v)
00253 {
00254   int sz = this->size();
00255 
00256   if (update_val_)
00257     this->val() /= v;
00258   for (int i=0; i<sz; ++i)
00259     this->fastAccessDx(i) /= v;
00260 
00261   return *this;
00262 }
00263 
00264 template <typename T, typename Storage> 
00265 inline  Sacado::ELRFad::GeneralFad<T,Storage>& 
00266 Sacado::ELRFad::GeneralFad<T,Storage>::
00267 operator += (const typename Sacado::dummy<value_type,scalar_type>::type& v)
00268 {
00269   if (update_val_)
00270     this->val() += v;
00271 
00272   return *this;
00273 }
00274 
00275 template <typename T, typename Storage> 
00276 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00277 Sacado::ELRFad::GeneralFad<T,Storage>::
00278 operator -= (const typename Sacado::dummy<value_type,scalar_type>::type& v)
00279 {
00280   if (update_val_)
00281     this->val() -= v;
00282 
00283   return *this;
00284 }
00285 
00286 template <typename T, typename Storage> 
00287 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00288 Sacado::ELRFad::GeneralFad<T,Storage>::
00289 operator *= (const typename Sacado::dummy<value_type,scalar_type>::type& v)
00290 {
00291   int sz = this->size();
00292 
00293   if (update_val_)
00294     this->val() *= v;
00295   for (int i=0; i<sz; ++i)
00296     this->fastAccessDx(i) *= v;
00297 
00298   return *this;
00299 }
00300 
00301 template <typename T, typename Storage> 
00302 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00303 Sacado::ELRFad::GeneralFad<T,Storage>::
00304 operator /= (const typename Sacado::dummy<value_type,scalar_type>::type& v)
00305 {
00306   int sz = this->size();
00307 
00308   if (update_val_)
00309     this->val() /= v;
00310   for (int i=0; i<sz; ++i)
00311     this->fastAccessDx(i) /= v;
00312 
00313   return *this;
00314 }
00315 
00316 template <typename T, typename Storage> 
00317 template <typename S> 
00318 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00319 Sacado::ELRFad::GeneralFad<T,Storage>::operator += (
00320                const Sacado::ELRFad::Expr<S>& x)
00321 {
00322   int xsz = x.size(), sz = this->size();
00323 
00324 #ifdef SACADO_DEBUG
00325   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00326     throw "Fad Error:  Attempt to assign with incompatible sizes";
00327 #endif
00328 
00329   // Number of arguments
00330   const int N = Expr<S>::num_args;
00331 
00332   if (xsz) {
00333 
00334     // Compute partials
00335     LocalAccumOp< Expr<S> > op(x);
00336 
00337     if (sz) {
00338 
00339       // Compute each tangent direction
00340       for(int i=0; i<xsz; ++i) {
00341   op.t = T(0.);
00342   op.i = i;
00343   
00344   // Automatically unrolled loop that computes
00345   // for (int j=0; j<N; j++)
00346   //   op.t += op.partials[j] * x.getTangent<j>(i);
00347   Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00348   
00349   this->fastAccessDx(i) += op.t;
00350       }
00351 
00352     }
00353 
00354     else {
00355 
00356       this->resize(xsz);
00357 
00358       // Compute each tangent direction
00359       for(int i=0; i<xsz; ++i) {
00360   op.t = T(0.);
00361   op.i = i;
00362   
00363   // Automatically unrolled loop that computes
00364   // for (int j=0; j<N; j++)
00365   //   op.t += op.partials[j] * x.getTangent<j>(i);
00366   Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00367   
00368   this->fastAccessDx(i) = op.t;
00369       }
00370 
00371     }
00372 
00373   }
00374   
00375   // Compute value
00376   update_val_ = x.updateValue();
00377   if (update_val_)
00378     this->val() += x.val();
00379 
00380   return *this;
00381 }
00382 
00383 template <typename T, typename Storage> 
00384 template <typename S> 
00385 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00386 Sacado::ELRFad::GeneralFad<T,Storage>::operator -= (
00387                 const Sacado::ELRFad::Expr<S>& x)
00388 {
00389   int xsz = x.size(), sz = this->size();
00390 
00391 #ifdef SACADO_DEBUG
00392   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00393     throw "Fad Error:  Attempt to assign with incompatible sizes";
00394 #endif
00395 
00396   // Number of arguments
00397   const int N = Expr<S>::num_args;
00398 
00399   if (xsz) {
00400 
00401     // Compute partials
00402     LocalAccumOp< Expr<S> > op(x);
00403 
00404     if (sz) {
00405 
00406       // Compute each tangent direction
00407       for(int i=0; i<xsz; ++i) {
00408   op.t = T(0.);
00409   op.i = i;
00410   
00411   // Automatically unrolled loop that computes
00412   // for (int j=0; j<N; j++)
00413   //   op.t += op.partials[j] * x.getTangent<j>(i);
00414   Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00415   
00416   this->fastAccessDx(i) -= op.t;
00417       }
00418 
00419     }
00420 
00421     else {
00422 
00423       this->resize(xsz);
00424 
00425       // Compute each tangent direction
00426       for(int i=0; i<xsz; ++i) {
00427   op.t = T(0.);
00428   op.i = i;
00429   
00430   // Automatically unrolled loop that computes
00431   // for (int j=0; j<N; j++)
00432   //   op.t += op.partials[j] * x.getTangent<j>(i);
00433   Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00434   
00435   this->fastAccessDx(i) = -op.t;
00436       }
00437 
00438     }
00439 
00440   }
00441 
00442   update_val_ = x.updateValue();
00443   if (update_val_)
00444     this->val() -= x.val();
00445 
00446   return *this;
00447 }
00448 
00449 template <typename T, typename Storage> 
00450 template <typename S> 
00451 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00452 Sacado::ELRFad::GeneralFad<T,Storage>::operator *= (
00453                const Sacado::ELRFad::Expr<S>& x)
00454 {
00455   int xsz = x.size(), sz = this->size();
00456   T xval = x.val();
00457 
00458 #ifdef SACADO_DEBUG
00459   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00460     throw "Fad Error:  Attempt to assign with incompatible sizes";
00461 #endif
00462 
00463   // Number of arguments
00464   const int N = Expr<S>::num_args;
00465 
00466   if (xsz) {
00467 
00468     // Compute partials
00469     LocalAccumOp< Expr<S> > op(x);
00470 
00471     if (sz) {
00472 
00473       // Compute each tangent direction
00474       for(int i=0; i<xsz; ++i) {
00475   op.t = T(0.);
00476   op.i = i;
00477   
00478   // Automatically unrolled loop that computes
00479   // for (int j=0; j<N; j++)
00480   //   op.t += op.partials[j] * x.getTangent<j>(i);
00481   Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00482   
00483   this->fastAccessDx(i) = this->val() * op.t + this->fastAccessDx(i) * xval;
00484       }
00485 
00486     }
00487 
00488     else {
00489 
00490       this->resize(xsz);
00491 
00492       // Compute each tangent direction
00493       for(int i=0; i<xsz; ++i) {
00494   op.t = T(0.);
00495   op.i = i;
00496   
00497   // Automatically unrolled loop that computes
00498   // for (int j=0; j<N; j++)
00499   //   op.t += op.partials[j] * x.getTangent<j>(i);
00500   Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00501   
00502   this->fastAccessDx(i) = this->val() * op.t;
00503       }
00504 
00505     }
00506 
00507   }
00508 
00509   else {
00510 
00511     if (sz) {
00512       for (int i=0; i<sz; ++i)
00513   this->fastAccessDx(i) *= xval;
00514     }
00515 
00516   }
00517 
00518   update_val_ = x.updateValue();
00519   if (update_val_)
00520     this->val() *= xval;
00521 
00522   return *this;
00523 }
00524 
00525 template <typename T, typename Storage>
00526 template <typename S> 
00527 inline Sacado::ELRFad::GeneralFad<T,Storage>& 
00528 Sacado::ELRFad::GeneralFad<T,Storage>::operator /= (
00529                const Sacado::ELRFad::Expr<S>& x)
00530 {
00531   int xsz = x.size(), sz = this->size();
00532   T xval = x.val();
00533 
00534 #ifdef SACADO_DEBUG
00535   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00536     throw "Fad Error:  Attempt to assign with incompatible sizes";
00537 #endif
00538 
00539   // Number of arguments
00540   const int N = Expr<S>::num_args;
00541 
00542   if (xsz) {
00543 
00544     // Compute partials
00545     LocalAccumOp< Expr<S> > op(x);
00546 
00547     T xval2 = xval*xval;
00548 
00549     if (sz) {
00550 
00551       // Compute each tangent direction
00552       for(int i=0; i<xsz; ++i) {
00553   op.t = T(0.);
00554   op.i = i;
00555   
00556   // Automatically unrolled loop that computes
00557   // for (int j=0; j<N; j++)
00558   //   op.t += op.partials[j] * x.getTangent<j>(i);
00559   Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00560   
00561   this->fastAccessDx(i) = (this->fastAccessDx(i) * xval - this->val() * op.t) / xval2;
00562       }
00563 
00564     }
00565 
00566     else {
00567 
00568       this->resize(xsz);
00569 
00570       // Compute each tangent direction
00571       for(int i=0; i<xsz; ++i) {
00572   op.t = T(0.);
00573   op.i = i;
00574   
00575   // Automatically unrolled loop that computes
00576   // for (int j=0; j<N; j++)
00577   //   op.t += op.partials[j] * x.getTangent<j>(i);
00578   Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00579   
00580   this->fastAccessDx(i) = -this->val() * op.t / xval2;
00581       }
00582 
00583     }
00584 
00585   }
00586 
00587   else {
00588 
00589     if (sz) {
00590       for (int i=0; i<sz; ++i)
00591   this->fastAccessDx(i) /= xval;
00592     }
00593 
00594   }
00595 
00596   update_val_ = x.updateValue();
00597   if (update_val_)
00598     this->val() /= xval;
00599 
00600   return *this;
00601 }
00602 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines