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