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