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