Sacado Package Browser (Single Doxygen Collection): Sacado_CacheFad

00001 // $Id: Sacado_CacheFad_GeneralFadImp.hpp,v 1.8 2007/07/24 00:08:04 etphipp Exp $ 
00002 // $Source: /space/CVS/Trilinos/packages/sacado/src/Sacado_CacheFad_GeneralFadImp.hpp,v $ 
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   s_(T(0.))
00060 {
00061   int sz = x.size();
00062 
00063   T xval = x.val(); 
00064 
00065   if (sz != s_.size()) 
00066     s_.resize(sz);
00067 
00068   if (sz) {
00069     if (x.hasFastAccess())
00070       for(int i=0; i<sz; ++i) 
00071   s_.dx_[i] = x.fastAccessDx(i);
00072     else
00073       for(int i=0; i<sz; ++i) 
00074   s_.dx_[i] = x.dx(i);
00075   }
00076 
00077   s_.val_ = xval;
00078 }
00079 
00080 
00081 template <typename T, typename Storage> 
00082 inline void 
00083 Sacado::CacheFad::GeneralFad<T,Storage>::diff(const int ith, const int n) 
00084 { 
00085   if (s_.size() == 0) 
00086     s_.resize(n);
00087 
00088   s_.zero();
00089   s_.dx_[ith] = T(1.);
00090 
00091 }
00092 
00093 template <typename T, typename Storage> 
00094 inline Sacado::CacheFad::GeneralFad<T,Storage>& 
00095 Sacado::CacheFad::GeneralFad<T,Storage>::operator=(const T& val) 
00096 {
00097   s_.val_ = val;
00098 
00099   if (s_.size()) {
00100     s_.zero();    // We need to zero out the array for future resizes
00101     s_.resize(0);
00102   }
00103 
00104   return *this;
00105 }
00106 
00107 template <typename T, typename Storage> 
00108 inline Sacado::CacheFad::GeneralFad<T,Storage>& 
00109 Sacado::CacheFad::GeneralFad<T,Storage>::operator=(
00110            const Sacado::CacheFad::GeneralFad<T,Storage>& x) 
00111 {
00112   // Copy val_ and dx_
00113   s_.operator=(x.s_);
00114   
00115   return *this;
00116 }
00117 
00118 template <typename T, typename Storage> 
00119 template <typename S> 
00120 inline Sacado::CacheFad::GeneralFad<T,Storage>& 
00121 Sacado::CacheFad::GeneralFad<T,Storage>::operator=(const Expr<S>& x) 
00122 {
00123   int sz = x.size();
00124 
00125   T xval = x.val();
00126 
00127   if (sz != s_.size()) 
00128     s_.resize(sz);
00129 
00130   if (sz) {
00131     if (x.hasFastAccess())
00132       for(int i=0; i<sz; ++i)
00133   s_.dx_[i] = x.fastAccessDx(i);
00134     else
00135       for(int i=0; i<sz; ++i)
00136   s_.dx_[i] = x.dx(i);
00137   }
00138   
00139   s_.val_ = xval;
00140   
00141   return *this;
00142 }
00143 
00144 template <typename T, typename Storage> 
00145 inline  Sacado::CacheFad::GeneralFad<T,Storage>& 
00146 Sacado::CacheFad::GeneralFad<T,Storage>::operator += (const T& val)
00147 {
00148   s_.val_ += val;
00149 
00150   return *this;
00151 }
00152 
00153 template <typename T, typename Storage> 
00154 inline Sacado::CacheFad::GeneralFad<T,Storage>& 
00155 Sacado::CacheFad::GeneralFad<T,Storage>::operator -= (const T& val)
00156 {
00157   s_.val_ -= val;
00158 
00159   return *this;
00160 }
00161 
00162 template <typename T, typename Storage> 
00163 inline Sacado::CacheFad::GeneralFad<T,Storage>& 
00164 Sacado::CacheFad::GeneralFad<T,Storage>::operator *= (const T& val)
00165 {
00166   int sz = s_.size();
00167 
00168   s_.val_ *= val;
00169   for (int i=0; i<sz; ++i)
00170     s_.dx_[i] *= val;
00171 
00172   return *this;
00173 }
00174 
00175 template <typename T, typename Storage> 
00176 inline Sacado::CacheFad::GeneralFad<T,Storage>& 
00177 Sacado::CacheFad::GeneralFad<T,Storage>::operator /= (const T& val)
00178 {
00179   int sz = s_.size();
00180 
00181   s_.val_ /= val;
00182   for (int i=0; i<sz; ++i)
00183     s_.dx_[i] /= val;
00184 
00185   return *this;
00186 }
00187 
00188 template <typename T, typename Storage> 
00189 template <typename S> 
00190 inline Sacado::CacheFad::GeneralFad<T,Storage>& 
00191 Sacado::CacheFad::GeneralFad<T,Storage>::operator += (const Sacado::CacheFad::Expr<S>& x)
00192 {
00193   int xsz = x.size(), sz = s_.size();
00194 
00195 #ifdef SACADO_DEBUG
00196   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00197     throw "Fad Error:  Attempt to assign with incompatible sizes";
00198 #endif
00199 
00200   T xval = x.val();
00201 
00202   if (xsz) {
00203     if (sz) {
00204       if (x.hasFastAccess())
00205   for (int i=0; i<sz; ++i)
00206     s_.dx_[i] += x.fastAccessDx(i);
00207       else
00208   for (int i=0; i<sz; ++i)
00209     s_.dx_[i] += x.dx(i);
00210     }
00211     else {
00212       s_.resize(xsz);
00213       if (x.hasFastAccess())
00214   for (int i=0; i<xsz; ++i)
00215     s_.dx_[i] = x.fastAccessDx(i);
00216       else
00217   for (int i=0; i<xsz; ++i)
00218     s_.dx_[i] = x.dx(i);
00219     }
00220   }
00221 
00222   s_.val_ += xval;
00223 
00224   return *this;
00225 }
00226 
00227 template <typename T, typename Storage> 
00228 template <typename S> 
00229 inline Sacado::CacheFad::GeneralFad<T,Storage>& 
00230 Sacado::CacheFad::GeneralFad<T,Storage>::operator -= (const Sacado::CacheFad::Expr<S>& x)
00231 {
00232   int xsz = x.size(), sz = s_.size();
00233 
00234 #ifdef SACADO_DEBUG
00235   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00236     throw "Fad Error:  Attempt to assign with incompatible sizes";
00237 #endif
00238 
00239   T xval = x.val();
00240 
00241   if (xsz) {
00242     if (sz) {
00243       if (x.hasFastAccess())
00244   for(int i=0; i<sz; ++i)
00245     s_.dx_[i] -= x.fastAccessDx(i);
00246       else
00247   for (int i=0; i<sz; ++i)
00248     s_.dx_[i] -= x.dx(i);
00249     }
00250     else {
00251       s_.resize(xsz);
00252       if (x.hasFastAccess())
00253   for(int i=0; i<xsz; ++i)
00254     s_.dx_[i] = -x.fastAccessDx(i);
00255       else
00256   for (int i=0; i<xsz; ++i)
00257     s_.dx_[i] = -x.dx(i);
00258     }
00259   }
00260 
00261   s_.val_ -= xval;
00262 
00263 
00264   return *this;
00265 }
00266 
00267 template <typename T, typename Storage> 
00268 template <typename S> 
00269 inline Sacado::CacheFad::GeneralFad<T,Storage>& 
00270 Sacado::CacheFad::GeneralFad<T,Storage>::operator *= (const Sacado::CacheFad::Expr<S>& x)
00271 {
00272   int xsz = x.size(), sz = s_.size();
00273   T xval = x.val();
00274 
00275 #ifdef SACADO_DEBUG
00276   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00277     throw "Fad Error:  Attempt to assign with incompatible sizes";
00278 #endif
00279 
00280   if (xsz) {
00281     if (sz) {
00282       if (x.hasFastAccess())
00283   for(int i=0; i<sz; ++i)
00284     s_.dx_[i] = s_.val_ * x.fastAccessDx(i) + s_.dx_[i] * xval;
00285       else
00286   for (int i=0; i<sz; ++i)
00287     s_.dx_[i] = s_.val_ * x.dx(i) + s_.dx_[i] * xval;
00288     }
00289     else {
00290       s_.resize(xsz);
00291       if (x.hasFastAccess())
00292   for(int i=0; i<xsz; ++i)
00293     s_.dx_[i] = s_.val_ * x.fastAccessDx(i);
00294       else
00295   for (int i=0; i<xsz; ++i)
00296     s_.dx_[i] = s_.val_ * x.dx(i);
00297     }
00298   }
00299   else {
00300     if (sz) {
00301       for (int i=0; i<sz; ++i)
00302   s_.dx_[i] *= xval;
00303     }
00304   }
00305 
00306   s_.val_ *= xval;
00307 
00308   return *this;
00309 }
00310 
00311 template <typename T, typename Storage>
00312 template <typename S> 
00313 inline Sacado::CacheFad::GeneralFad<T,Storage>& 
00314 Sacado::CacheFad::GeneralFad<T,Storage>::operator /= (const Sacado::CacheFad::Expr<S>& x)
00315 {
00316   int xsz = x.size(), sz = s_.size();
00317   T xval = x.val();
00318 
00319 #ifdef SACADO_DEBUG
00320   if ((xsz != sz) && (xsz != 0) && (sz != 0))
00321     throw "Fad Error:  Attempt to assign with incompatible sizes";
00322 #endif
00323 
00324   if (xsz) {
00325     if (sz) {
00326       if (x.hasFastAccess())
00327   for(int i=0; i<sz; ++i)
00328     s_.dx_[i] = ( s_.dx_[i]*xval - s_.val_*x.fastAccessDx(i) )/ (xval*xval);
00329       else
00330   for (int i=0; i<sz; ++i)
00331     s_.dx_[i] = ( s_.dx_[i]*xval - s_.val_*x.dx(i) )/ (xval*xval);
00332     }
00333     else {
00334       s_.resize(xsz);
00335       if (x.hasFastAccess())
00336   for(int i=0; i<xsz; ++i)
00337     s_.dx_[i] = - s_.val_*x.fastAccessDx(i) / (xval*xval);
00338       else
00339   for (int i=0; i<xsz; ++i)
00340     s_.dx_[i] = -s_.val_ * x.dx(i) / (xval*xval);
00341     }
00342   }
00343   else {
00344     if (sz) {
00345       for (int i=0; i<sz; ++i)
00346   s_.dx_[i] /= xval;
00347     }
00348   }
00349 
00350   s_.val_ /= xval;
00351 
00352   return *this;
00353 }
00354