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