|
Sacado Package Browser (Single Doxygen Collection) Version of the Day
|
|
Sacado Package Browser (Single Doxygen Collection) Version of the Day
|
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
1.7.4
