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