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