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