|
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 ValT, typename LogT> 00057 template <typename S> 00058 inline Sacado::LFad::LogicalSparseImp<ValT,LogT>:: 00059 LogicalSparseImp(const Expr<S>& x) : 00060 Storage(value_type(0)) 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 this->val() = x.val(); 00077 } 00078 00079 00080 template <typename ValT, typename LogT> 00081 inline void 00082 Sacado::LFad::LogicalSparseImp<ValT,LogT>:: 00083 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) = logical_type(1); 00090 00091 } 00092 00093 template <typename ValT, typename LogT> 00094 inline Sacado::LFad::LogicalSparseImp<ValT,LogT>& 00095 Sacado::LFad::LogicalSparseImp<ValT,LogT>:: 00096 operator=(const ValT& v) 00097 { 00098 this->val() = v; 00099 00100 if (this->size()) { 00101 this->zero(); // We need to zero out the array for future resizes 00102 this->resize(0); 00103 } 00104 00105 return *this; 00106 } 00107 00108 template <typename ValT, typename LogT> 00109 inline Sacado::LFad::LogicalSparseImp<ValT,LogT>& 00110 Sacado::LFad::LogicalSparseImp<ValT,LogT>:: 00111 operator=(const Sacado::LFad::LogicalSparseImp<ValT,LogT>& x) 00112 { 00113 // Copy value & dx_ 00114 Storage::operator=(x); 00115 00116 return *this; 00117 } 00118 00119 template <typename ValT, typename LogT> 00120 template <typename S> 00121 inline Sacado::LFad::LogicalSparseImp<ValT,LogT>& 00122 Sacado::LFad::LogicalSparseImp<ValT,LogT>:: 00123 operator=(const Expr<S>& x) 00124 { 00125 int sz = x.size(); 00126 00127 if (sz != this->size()) 00128 this->resize(sz); 00129 00130 if (sz) { 00131 if (x.hasFastAccess()) 00132 for(int i=0; i<sz; ++i) 00133 this->fastAccessDx(i) = x.fastAccessDx(i); 00134 else 00135 for(int i=0; i<sz; ++i) 00136 this->fastAccessDx(i) = x.dx(i); 00137 } 00138 00139 this->val() = x.val(); 00140 00141 return *this; 00142 } 00143 00144 template <typename ValT, typename LogT> 00145 inline Sacado::LFad::LogicalSparseImp<ValT,LogT>& 00146 Sacado::LFad::LogicalSparseImp<ValT,LogT>:: 00147 operator += (const ValT& v) 00148 { 00149 this->val() += v; 00150 00151 return *this; 00152 } 00153 00154 template <typename ValT, typename LogT> 00155 inline Sacado::LFad::LogicalSparseImp<ValT,LogT>& 00156 Sacado::LFad::LogicalSparseImp<ValT,LogT>:: 00157 operator -= (const ValT& v) 00158 { 00159 this->val() -= v; 00160 00161 return *this; 00162 } 00163 00164 template <typename ValT, typename LogT> 00165 inline Sacado::LFad::LogicalSparseImp<ValT,LogT>& 00166 Sacado::LFad::LogicalSparseImp<ValT,LogT>:: 00167 operator *= (const ValT& v) 00168 { 00169 this->val() *= v; 00170 00171 return *this; 00172 } 00173 00174 template <typename ValT, typename LogT> 00175 inline Sacado::LFad::LogicalSparseImp<ValT,LogT>& 00176 Sacado::LFad::LogicalSparseImp<ValT,LogT>:: 00177 operator /= (const ValT& v) 00178 { 00179 this->val() /= v; 00180 00181 return *this; 00182 } 00183 00184 template <typename ValT, typename LogT> 00185 template <typename S> 00186 inline Sacado::LFad::LogicalSparseImp<ValT,LogT>& 00187 Sacado::LFad::LogicalSparseImp<ValT,LogT>:: 00188 operator += (const Sacado::LFad::Expr<S>& x) 00189 { 00190 int xsz = x.size(), sz = this->size(); 00191 00192 #ifdef SACADO_DEBUG 00193 if ((xsz != sz) && (xsz != 0) && (sz != 0)) 00194 throw "LFad Error: Attempt to assign with incompatible sizes"; 00195 #endif 00196 00197 if (xsz) { 00198 if (sz) { 00199 if (x.hasFastAccess()) 00200 for (int i=0; i<xsz; ++i) 00201 this->fastAccessDx(i) = this->fastAccessDx(i) || x.fastAccessDx(i); 00202 else 00203 for (int i=0; i<xsz; ++i) 00204 this->fastAccessDx(i) = this->fastAccessDx(i) || x.dx(i); 00205 } 00206 else { 00207 this->resize(xsz); 00208 if (x.hasFastAccess()) 00209 for (int i=0; i<xsz; ++i) 00210 this->fastAccessDx(i) = x.fastAccessDx(i); 00211 else 00212 for (int i=0; i<xsz; ++i) 00213 this->fastAccessDx(i) = x.dx(i); 00214 } 00215 } 00216 00217 this->val() += x.val(); 00218 00219 return *this; 00220 } 00221 00222 template <typename ValT, typename LogT> 00223 template <typename S> 00224 inline Sacado::LFad::LogicalSparseImp<ValT,LogT>& 00225 Sacado::LFad::LogicalSparseImp<ValT,LogT>:: 00226 operator -= (const Sacado::LFad::Expr<S>& x) 00227 { 00228 int xsz = x.size(), sz = this->size(); 00229 00230 #ifdef SACADO_DEBUG 00231 if ((xsz != sz) && (xsz != 0) && (sz != 0)) 00232 throw "LFad Error: Attempt to assign with incompatible sizes"; 00233 #endif 00234 00235 if (xsz) { 00236 if (sz) { 00237 if (x.hasFastAccess()) 00238 for (int i=0; i<xsz; ++i) 00239 this->fastAccessDx(i) = this->fastAccessDx(i) || x.fastAccessDx(i); 00240 else 00241 for (int i=0; i<xsz; ++i) 00242 this->fastAccessDx(i) = this->fastAccessDx(i) || x.dx(i); 00243 } 00244 else { 00245 this->resize(xsz); 00246 if (x.hasFastAccess()) 00247 for (int i=0; i<xsz; ++i) 00248 this->fastAccessDx(i) = x.fastAccessDx(i); 00249 else 00250 for (int i=0; i<xsz; ++i) 00251 this->fastAccessDx(i) = x.dx(i); 00252 } 00253 } 00254 00255 this->val() -= x.val(); 00256 00257 00258 return *this; 00259 } 00260 00261 template <typename ValT, typename LogT> 00262 template <typename S> 00263 inline Sacado::LFad::LogicalSparseImp<ValT,LogT>& 00264 Sacado::LFad::LogicalSparseImp<ValT,LogT>:: 00265 operator *= (const Sacado::LFad::Expr<S>& x) 00266 { 00267 int xsz = x.size(), sz = this->size(); 00268 00269 #ifdef SACADO_DEBUG 00270 if ((xsz != sz) && (xsz != 0) && (sz != 0)) 00271 throw "LFad Error: Attempt to assign with incompatible sizes"; 00272 #endif 00273 00274 if (xsz) { 00275 if (sz) { 00276 if (x.hasFastAccess()) 00277 for (int i=0; i<xsz; ++i) 00278 this->fastAccessDx(i) = this->fastAccessDx(i) || x.fastAccessDx(i); 00279 else 00280 for (int i=0; i<xsz; ++i) 00281 this->fastAccessDx(i) = this->fastAccessDx(i) || x.dx(i); 00282 } 00283 else { 00284 this->resize(xsz); 00285 if (x.hasFastAccess()) 00286 for (int i=0; i<xsz; ++i) 00287 this->fastAccessDx(i) = x.fastAccessDx(i); 00288 else 00289 for (int i=0; i<xsz; ++i) 00290 this->fastAccessDx(i) = x.dx(i); 00291 } 00292 } 00293 00294 this->val() *= x.val(); 00295 00296 return *this; 00297 } 00298 00299 template <typename ValT, typename LogT> 00300 template <typename S> 00301 inline Sacado::LFad::LogicalSparseImp<ValT,LogT>& 00302 Sacado::LFad::LogicalSparseImp<ValT,LogT>:: 00303 operator /= (const Sacado::LFad::Expr<S>& x) 00304 { 00305 int xsz = x.size(), sz = this->size(); 00306 00307 #ifdef SACADO_DEBUG 00308 if ((xsz != sz) && (xsz != 0) && (sz != 0)) 00309 throw "LFad Error: Attempt to assign with incompatible sizes"; 00310 #endif 00311 00312 if (xsz) { 00313 if (sz) { 00314 if (x.hasFastAccess()) 00315 for (int i=0; i<xsz; ++i) 00316 this->fastAccessDx(i) = this->fastAccessDx(i) || x.fastAccessDx(i); 00317 else 00318 for (int i=0; i<xsz; ++i) 00319 this->fastAccessDx(i) = this->fastAccessDx(i) || x.dx(i); 00320 } 00321 else { 00322 this->resize(xsz); 00323 if (x.hasFastAccess()) 00324 for (int i=0; i<xsz; ++i) 00325 this->fastAccessDx(i) = x.fastAccessDx(i); 00326 else 00327 for (int i=0; i<xsz; ++i) 00328 this->fastAccessDx(i) = x.dx(i); 00329 } 00330 } 00331 00332 this->val() /= x.val(); 00333 00334 return *this; 00335 } 00336
1.7.4