Sacado Package Browser (Single Doxygen Collection) Version of the Day
Sacado_ELRCacheFad_SFadImp.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, int Num> 
00059 inline Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00060 Expr(const int sz, const T & x) : val_(x), update_val_(true)
00061 {
00062 #ifdef SACADO_DEBUG
00063   if (sz != Num)
00064     throw "ELRCacheFad::SFad() Error:  Supplied derivative dimension does not match compile time length.";
00065 #endif
00066 
00067   ss_array<T>::zero(dx_, Num); 
00068 }
00069 
00070 template <typename T, int Num> 
00071 inline Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00072 Expr(const int sz, const int i, const T & x) : val_(x), update_val_(true)
00073 {
00074 #ifdef SACADO_DEBUG
00075   if (sz != Num)
00076     throw "ELRCacheFad::SFad() Error:  Supplied derivative dimension does not match compile time length.";
00077   if (i >= Num)
00078     throw "ELRCacheFad::SFad() Error:  Invalid derivative index.";
00079 #endif
00080 
00081   ss_array<T>::zero(dx_, Num);
00082   dx_[i]=1.; 
00083 }
00084 
00085 template <typename T, int Num> 
00086 template <typename S> 
00087 inline Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00088 Expr(const Expr<S>& x) : update_val_(x.updateValue())
00089 {
00090 #ifdef SACADO_DEBUG
00091   if (x.size() != Num)
00092     throw "ELRCacheFad::SFad() Error:  Attempt to assign with incompatible sizes";
00093 #endif
00094 
00095   x.cache();
00096 
00097   if (update_val_)
00098     this->val() = x.val();
00099 
00100   // Number of arguments
00101   const int N = Expr<S>::num_args;
00102 
00103   // Compute partials
00104   LocalAccumOp< Expr<S> > op(x);
00105   
00106   // Compute each tangent direction
00107   for(int i=0; i<Num; ++i) {
00108     op.t = T(0.);
00109     op.i = i;
00110 
00111     // Automatically unrolled loop that computes
00112     // for (int j=0; j<N; j++)
00113     //   op.t += op.partials[j] * x.getTangent<j>(i);
00114     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00115 
00116     dx_[i] = op.t;
00117 
00118   }
00119 }
00120 
00121 
00122 template <typename T, int Num> 
00123 inline void 
00124 Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00125 diff(const int ith, const int n) 
00126 { 
00127 #ifdef SACADO_DEBUG
00128   if (n != Num)
00129     throw "ELRCacheFad::diff() Error:  Supplied derivative dimension does not match compile time length.";
00130 #endif
00131 
00132   ss_array<T>::zero(dx_, Num);
00133   dx_[ith] = T(1.);
00134 }
00135 
00136 template <typename T, int Num> 
00137 inline void 
00138 Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00139 resize(int sz)
00140 {
00141 #ifdef SACADO_DEBUG
00142   if (sz != Num)
00143     throw "ELRCacheFad::resize() Error:  Cannot resize fixed derivative array dimension";
00144 #endif
00145 }
00146 
00147 template <typename T, int Num> 
00148 inline Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >& 
00149 Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00150 operator=(const T& v) 
00151 {
00152   val_ = v;
00153   ss_array<T>::zero(dx_, Num);
00154 
00155   return *this;
00156 }
00157 
00158 template <typename T, int Num> 
00159 inline Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >& 
00160 Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00161 operator=(const Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >& x) 
00162 {
00163   // Copy value
00164   val_ = x.val_;
00165 
00166   // Copy dx_
00167   //ss_array<T>::copy(x.dx_, dx_, Num);
00168   for (int i=0; i<Num; i++)
00169     dx_[i] = x.dx_[i];
00170 
00171   update_val_ = x.update_val_;
00172   
00173   return *this;
00174 }
00175 
00176 template <typename T, int Num> 
00177 template <typename S> 
00178 inline Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >& 
00179 Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00180 operator=(const Expr<S>& x) 
00181 {
00182 #ifdef SACADO_DEBUG
00183   if (x.size() != Num)
00184     throw "ELRCacheFad::operator=() Error:  Attempt to assign with incompatible sizes";
00185 #endif
00186 
00187   x.cache();
00188 
00189   // Number of arguments
00190   const int N = Expr<S>::num_args;
00191 
00192   // Compute partials
00193   LocalAccumOp< Expr<S> > op(x);
00194   
00195   // Compute each tangent direction
00196   for(int i=0; i<Num; ++i) {
00197     op.t = T(0.);
00198     op.i = i;
00199 
00200     // Automatically unrolled loop that computes
00201     // for (int j=0; j<N; j++)
00202     //   op.t += op.partials[j] * x.getTangent<j>(i);
00203     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00204     
00205     dx_[i] = op.t;
00206   }
00207   
00208   update_val_ = x.updateValue();
00209   if (update_val_)
00210     this->val() = x.val();
00211   
00212   return *this;
00213 }
00214 
00215 template <typename T, int Num> 
00216 inline  Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >& 
00217 Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00218 operator += (const T& v)
00219 {
00220   if (update_val_)
00221     val_ += v;
00222 
00223   return *this;
00224 }
00225 
00226 template <typename T, int Num> 
00227 inline Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >& 
00228 Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00229 operator -= (const T& v)
00230 {
00231   if (update_val_)
00232     val_ -= v;
00233 
00234   return *this;
00235 }
00236 
00237 template <typename T, int Num> 
00238 inline Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >& 
00239 Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00240 operator *= (const T& v)
00241 {
00242   if (update_val_)
00243     val_ *= v;
00244 
00245   for (int i=0; i<Num; ++i)
00246     dx_[i] *= v;
00247 
00248   return *this;
00249 }
00250 
00251 template <typename T, int Num> 
00252 inline Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >& 
00253 Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00254 operator /= (const T& v)
00255 {
00256   if (update_val_)
00257     val_ /= v;
00258 
00259   for (int i=0; i<Num; ++i)
00260     dx_[i] /= v;
00261 
00262   return *this;
00263 }
00264 
00265 template <typename T, int Num> 
00266 template <typename S> 
00267 inline Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >& 
00268 Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00269 operator += (const Sacado::ELRCacheFad::Expr<S>& x)
00270 {
00271 #ifdef SACADO_DEBUG
00272   if (x.size() != Num)
00273     throw "ELRCacheFad::operator+=() Error:  Attempt to assign with incompatible sizes";
00274 #endif
00275 
00276   x.cache();
00277 
00278   // Number of arguments
00279   const int N = Expr<S>::num_args;
00280 
00281   // Compute partials
00282   LocalAccumOp< Expr<S> > op(x);
00283 
00284   // Compute each tangent direction
00285   for(int i=0; i<Num; ++i) {
00286     op.t = T(0.);
00287     op.i = i;
00288   
00289     // Automatically unrolled loop that computes
00290     // for (int j=0; j<N; j++)
00291     //   op.t += op.partials[j] * x.getTangent<j>(i);
00292     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00293     
00294     dx_[i] += op.t;
00295   }
00296  
00297   update_val_ = x.updateValue();
00298   if (update_val_)
00299     this->val() += x.val();
00300 
00301   return *this;
00302 }
00303 
00304 template <typename T, int Num> 
00305 template <typename S> 
00306 inline Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >& 
00307 Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00308 operator -= (const Sacado::ELRCacheFad::Expr<S>& x)
00309 {
00310 #ifdef SACADO_DEBUG
00311   if (x.size() != Num)
00312     throw "ELRCacheFad::operator-=() Error:  Attempt to assign with incompatible sizes";
00313 #endif
00314 
00315   x.cache();
00316 
00317   // Number of arguments
00318   const int N = Expr<S>::num_args;
00319 
00320   // Compute partials
00321   LocalAccumOp< Expr<S> > op(x);
00322 
00323   // Compute each tangent direction
00324   for(int i=0; i<Num; ++i) {
00325     op.t = T(0.);
00326     op.i = i;
00327   
00328     // Automatically unrolled loop that computes
00329     // for (int j=0; j<N; j++)
00330     //   op.t += op.partials[j] * x.getTangent<j>(i);
00331     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00332   
00333     dx_[i] -= op.t;
00334   }
00335 
00336   update_val_ = x.updateValue();
00337   if (update_val_)
00338     this->val() -= x.val();
00339 
00340   return *this;
00341 }
00342 
00343 template <typename T, int Num> 
00344 template <typename S> 
00345 inline Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >& 
00346 Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00347 operator *= (const Sacado::ELRCacheFad::Expr<S>& x)
00348 {
00349   x.cache();
00350 
00351   T xval = x.val();
00352 
00353 #ifdef SACADO_DEBUG
00354   if (x.size() != Num)
00355     throw "ELRCacheFad::operator*=() Error:  Attempt to assign with incompatible sizes";
00356 #endif
00357 
00358   // Number of arguments
00359   const int N = Expr<S>::num_args;
00360 
00361   // Compute partials
00362   LocalAccumOp< Expr<S> > op(x);
00363 
00364   // Compute each tangent direction
00365   for(int i=0; i<Num; ++i) {
00366     op.t = T(0.);
00367     op.i = i;
00368   
00369     // Automatically unrolled loop that computes
00370     // for (int j=0; j<N; j++)
00371     //   op.t += op.partials[j] * x.getTangent<j>(i);
00372     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00373     
00374     dx_[i] = val_ * op.t + dx_[i] * xval;
00375   }
00376  
00377   // Compute value
00378   update_val_ = x.updateValue();
00379   if (update_val_)
00380     val_ *= xval;
00381 
00382   return *this;
00383 }
00384 
00385 template <typename T, int Num>
00386 template <typename S> 
00387 inline Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >& 
00388 Sacado::ELRCacheFad::Expr< Sacado::ELRCacheFad::SFadExprTag<T,Num> >::
00389 operator /= (const Sacado::ELRCacheFad::Expr<S>& x)
00390 {
00391   x.cache();
00392 
00393   T xval = x.val();
00394 
00395 #ifdef SACADO_DEBUG
00396   if (x.size() != Num)
00397     throw "ELRCacheFad::operator/=() Error:  Attempt to assign with incompatible sizes";
00398 #endif
00399 
00400   // Number of arguments
00401   const int N = Expr<S>::num_args;
00402 
00403   // Compute partials
00404   LocalAccumOp< Expr<S> > op(x);
00405 
00406   T xval2 = xval*xval;
00407     
00408   // Compute each tangent direction
00409   for(int i=0; i<Num; ++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     dx_[i] = (dx_[i] * xval - val_ * op.t) / xval2;
00419   }
00420 
00421   // Compute value
00422   update_val_ = x.updateValue();
00423   if (update_val_)
00424     val_ /= xval;
00425 
00426   return *this;
00427 }
00428 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines