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)
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) 
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) : val_(x.val())
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 
00116 
00117 template <typename T, int Num> 
00118 inline void 
00119 Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >::
00120 diff(const int ith, const int n) 
00121 { 
00122 #ifdef SACADO_DEBUG
00123   if (n != Num)
00124     throw "SELRFad::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::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >::
00134 resize(int sz)
00135 {
00136 #ifdef SACADO_DEBUG
00137   if (sz != Num)
00138     throw "SELRFad::resize() Error:  Cannot resize fixed derivative array dimension";
00139 #endif
00140 }
00141 
00142 template <typename T, int Num> 
00143 inline Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >& 
00144 Sacado::ELRFad::Expr< Sacado::ELRFad::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::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >& 
00155 Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >::
00156 operator=(const Sacado::ELRFad::Expr< Sacado::ELRFad::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::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >& 
00172 Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >::
00173 operator=(const Expr<S>& x) 
00174 {
00175 #ifdef SACADO_DEBUG
00176   if (x.size() != Num)
00177     throw "SELRFad::operator=() Error:  Attempt to assign with incompatible sizes";
00178 #endif
00179 
00180   // Number of arguments
00181   const int N = Expr<S>::num_args;
00182 
00183   // Compute partials
00184   LocalAccumOp< Expr<S> > op(x);
00185   
00186   // Compute each tangent direction
00187   for(int i=0; i<Num; ++i) {
00188     op.t = T(0.);
00189     op.i = i;
00190 
00191     // Automatically unrolled loop that computes
00192     // for (int j=0; j<N; j++)
00193     //   op.t += op.partials[j] * x.getTangent<j>(i);
00194     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00195     
00196     dx_[i] = op.t;
00197   }
00198   
00199   // Compute value
00200   val_ = x.val();
00201   
00202   return *this;
00203 }
00204 
00205 template <typename T, int Num> 
00206 inline  Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >& 
00207 Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >::
00208 operator += (const T& v)
00209 {
00210   val_ += v;
00211 
00212   return *this;
00213 }
00214 
00215 template <typename T, int Num> 
00216 inline Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >& 
00217 Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >::
00218 operator -= (const T& v)
00219 {
00220   val_ -= v;
00221 
00222   return *this;
00223 }
00224 
00225 template <typename T, int Num> 
00226 inline Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >& 
00227 Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >::
00228 operator *= (const T& v)
00229 {
00230   val_ *= v;
00231 
00232   for (int i=0; i<Num; ++i)
00233     dx_[i] *= v;
00234 
00235   return *this;
00236 }
00237 
00238 template <typename T, int Num> 
00239 inline Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >& 
00240 Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >::
00241 operator /= (const T& v)
00242 {
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 template <typename S> 
00253 inline Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >& 
00254 Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >::
00255 operator += (const Sacado::ELRFad::Expr<S>& x)
00256 {
00257 #ifdef SACADO_DEBUG
00258   if (x.size() != Num)
00259     throw "SELRFad::operator+=() Error:  Attempt to assign with incompatible sizes";
00260 #endif
00261 
00262   // Number of arguments
00263   const int N = Expr<S>::num_args;
00264 
00265   // Compute partials
00266   LocalAccumOp< Expr<S> > op(x);
00267 
00268   // Compute each tangent direction
00269   for(int i=0; i<Num; ++i) {
00270     op.t = T(0.);
00271     op.i = i;
00272   
00273     // Automatically unrolled loop that computes
00274     // for (int j=0; j<N; j++)
00275     //   op.t += op.partials[j] * x.getTangent<j>(i);
00276     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00277     
00278     dx_[i] += op.t;
00279   }
00280  
00281   // Compute value
00282   val_ += x.val();
00283 
00284   return *this;
00285 }
00286 
00287 template <typename T, int Num> 
00288 template <typename S> 
00289 inline Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >& 
00290 Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >::
00291 operator -= (const Sacado::ELRFad::Expr<S>& x)
00292 {
00293 #ifdef SACADO_DEBUG
00294   if (x.size() != Num)
00295     throw "SELRFad::operator-=() Error:  Attempt to assign with incompatible sizes";
00296 #endif
00297 
00298   // Number of arguments
00299   const int N = Expr<S>::num_args;
00300 
00301   // Compute partials
00302   LocalAccumOp< Expr<S> > op(x);
00303 
00304   // Compute each tangent direction
00305   for(int i=0; i<Num; ++i) {
00306     op.t = T(0.);
00307     op.i = i;
00308   
00309     // Automatically unrolled loop that computes
00310     // for (int j=0; j<N; j++)
00311     //   op.t += op.partials[j] * x.getTangent<j>(i);
00312     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00313   
00314     dx_[i] -= op.t;
00315   }
00316 
00317   // Compute value
00318   val_ -= x.val();
00319 
00320   return *this;
00321 }
00322 
00323 template <typename T, int Num> 
00324 template <typename S> 
00325 inline Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >& 
00326 Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >::
00327 operator *= (const Sacado::ELRFad::Expr<S>& x)
00328 {
00329   T xval = x.val();
00330 
00331 #ifdef SACADO_DEBUG
00332   if (x.size() != Num)
00333     throw "SELRFad::operator*=() Error:  Attempt to assign with incompatible sizes";
00334 #endif
00335 
00336   // Number of arguments
00337   const int N = Expr<S>::num_args;
00338 
00339   // Compute partials
00340   LocalAccumOp< Expr<S> > op(x);
00341 
00342   // Compute each tangent direction
00343   for(int i=0; i<Num; ++i) {
00344     op.t = T(0.);
00345     op.i = i;
00346   
00347     // Automatically unrolled loop that computes
00348     // for (int j=0; j<N; j++)
00349     //   op.t += op.partials[j] * x.getTangent<j>(i);
00350     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00351     
00352     dx_[i] = val_ * op.t + dx_[i] * xval;
00353   }
00354  
00355   // Compute value
00356   val_ *= xval;
00357 
00358   return *this;
00359 }
00360 
00361 template <typename T, int Num>
00362 template <typename S> 
00363 inline Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >& 
00364 Sacado::ELRFad::Expr< Sacado::ELRFad::SFadExprTag<T,Num> >::
00365 operator /= (const Sacado::ELRFad::Expr<S>& x)
00366 {
00367   T xval = x.val();
00368 
00369 #ifdef SACADO_DEBUG
00370   if (x.size() != Num)
00371     throw "SELRFad::operator/=() Error:  Attempt to assign with incompatible sizes";
00372 #endif
00373 
00374   // Number of arguments
00375   const int N = Expr<S>::num_args;
00376 
00377   // Compute partials
00378   LocalAccumOp< Expr<S> > op(x);
00379 
00380   T xval2 = xval*xval;
00381     
00382   // Compute each tangent direction
00383   for(int i=0; i<Num; ++i) {
00384     op.t = T(0.);
00385     op.i = i;
00386     
00387     // Automatically unrolled loop that computes
00388     // for (int j=0; j<N; j++)
00389     //   op.t += op.partials[j] * x.getTangent<j>(i);
00390     Sacado::mpl::for_each< mpl::range_c< int, 0, N > > f(op);
00391     
00392     dx_[i] = (dx_[i] * xval - val_ * op.t) / xval2;
00393   }
00394 
00395   // Compute value
00396   val_ /= xval;
00397 
00398   return *this;
00399 }
00400 

Generated on Wed May 12 21:39:33 2010 for Sacado Package Browser (Single Doxygen Collection) by  doxygen 1.4.7