NGSolve
5.3
|
00001 #ifndef FILE_AUTODIFF 00002 #define FILE_AUTODIFF 00003 00004 /**************************************************************************/ 00005 /* File: autodiff.hpp */ 00006 /* Author: Joachim Schoeberl */ 00007 /* Date: 24. Oct. 02 */ 00008 /**************************************************************************/ 00009 00010 namespace ngstd 00011 { 00012 00013 // Automatic differentiation datatype 00014 00015 00021 template <int D, typename SCAL = double> 00022 class AutoDiff 00023 { 00024 SCAL val; 00025 SCAL dval[D?D:1]; 00026 public: 00027 00028 typedef AutoDiff<D,SCAL> TELEM; 00029 typedef SCAL TSCAL; 00030 00031 00033 INLINE AutoDiff () throw() { }; 00034 // { val = 0; for (int i = 0; i < D; i++) dval[i] = 0; } // ! 00035 00037 INLINE AutoDiff (const AutoDiff & ad2) throw() 00038 { 00039 val = ad2.val; 00040 for (int i = 0; i < D; i++) 00041 dval[i] = ad2.dval[i]; 00042 } 00043 00045 INLINE AutoDiff (SCAL aval) throw() 00046 { 00047 val = aval; 00048 for (int i = 0; i < D; i++) 00049 dval[i] = 0; 00050 } 00051 00053 INLINE AutoDiff (SCAL aval, int diffindex) throw() 00054 { 00055 val = aval; 00056 for (int i = 0; i < D; i++) 00057 dval[i] = 0; 00058 dval[diffindex] = 1; 00059 } 00060 00061 INLINE AutoDiff (SCAL aval, const SCAL * grad) 00062 { 00063 val = aval; 00064 LoadGradient (grad); 00065 } 00066 00068 INLINE AutoDiff & operator= (SCAL aval) throw() 00069 { 00070 val = aval; 00071 for (int i = 0; i < D; i++) 00072 dval[i] = 0; 00073 return *this; 00074 } 00075 00077 INLINE SCAL Value() const throw() { return val; } 00078 00080 INLINE SCAL DValue (int i) const throw() { return dval[i]; } 00081 00083 INLINE void StoreGradient (SCAL * p) const 00084 { 00085 for (int i = 0; i < D; i++) 00086 p[i] = dval[i]; 00087 } 00088 00089 INLINE void LoadGradient (const SCAL * p) 00090 { 00091 for (int i = 0; i < D; i++) 00092 dval[i] = p[i]; 00093 } 00094 00096 INLINE SCAL & Value() throw() { return val; } 00097 00099 INLINE SCAL & DValue (int i) throw() { return dval[i]; } 00100 }; 00101 00102 00104 00106 template<int D, typename SCAL> 00107 inline ostream & operator<< (ostream & ost, const AutoDiff<D,SCAL> & x) 00108 { 00109 ost << x.Value() << ", D = "; 00110 for (int i = 0; i < D; i++) 00111 ost << x.DValue(i) << " "; 00112 return ost; 00113 } 00114 00116 template<int D, typename SCAL> 00117 INLINE AutoDiff<D,SCAL> operator+ (const AutoDiff<D,SCAL> & x, const AutoDiff<D,SCAL> & y) throw() 00118 { 00119 AutoDiff<D,SCAL> res; 00120 res.Value () = x.Value()+y.Value(); 00121 // AutoDiff<D,SCAL> res(x.Value()+y.Value()); 00122 for (int i = 0; i < D; i++) 00123 res.DValue(i) = x.DValue(i) + y.DValue(i); 00124 return res; 00125 } 00126 00127 00129 template<int D, typename SCAL> 00130 INLINE AutoDiff<D,SCAL> operator- (const AutoDiff<D,SCAL> & x, const AutoDiff<D,SCAL> & y) throw() 00131 { 00132 AutoDiff<D,SCAL> res; 00133 res.Value() = x.Value()-y.Value(); 00134 // AutoDiff<D,SCAL> res (x.Value()-y.Value()); 00135 for (int i = 0; i < D; i++) 00136 res.DValue(i) = x.DValue(i) - y.DValue(i); 00137 return res; 00138 } 00139 00141 template<int D, typename SCAL> 00142 INLINE AutoDiff<D,SCAL> operator+ (double x, const AutoDiff<D,SCAL> & y) throw() 00143 { 00144 AutoDiff<D,SCAL> res; 00145 res.Value() = x+y.Value(); 00146 for (int i = 0; i < D; i++) 00147 res.DValue(i) = y.DValue(i); 00148 return res; 00149 } 00150 00152 template<int D, typename SCAL> 00153 INLINE AutoDiff<D,SCAL> operator+ (const AutoDiff<D,SCAL> & y, double x) throw() 00154 { 00155 AutoDiff<D,SCAL> res; 00156 res.Value() = x+y.Value(); 00157 for (int i = 0; i < D; i++) 00158 res.DValue(i) = y.DValue(i); 00159 return res; 00160 } 00161 00162 00164 template<int D, typename SCAL> 00165 INLINE AutoDiff<D,SCAL> operator- (const AutoDiff<D,SCAL> & x) throw() 00166 { 00167 AutoDiff<D,SCAL> res; 00168 res.Value() = -x.Value(); 00169 for (int i = 0; i < D; i++) 00170 res.DValue(i) = -x.DValue(i); 00171 return res; 00172 } 00173 00175 template<int D, typename SCAL> 00176 INLINE AutoDiff<D,SCAL> operator- (const AutoDiff<D,SCAL> & x, double y) throw() 00177 { 00178 AutoDiff<D,SCAL> res; 00179 res.Value() = x.Value()-y; 00180 for (int i = 0; i < D; i++) 00181 res.DValue(i) = x.DValue(i); 00182 return res; 00183 } 00184 00186 template<int D, typename SCAL> 00187 INLINE AutoDiff<D,SCAL> operator- (double x, const AutoDiff<D,SCAL> & y) throw() 00188 { 00189 AutoDiff<D,SCAL> res; 00190 res.Value() = x-y.Value(); 00191 for (int i = 0; i < D; i++) 00192 res.DValue(i) = -y.DValue(i); 00193 return res; 00194 } 00195 00196 00198 template<int D, typename SCAL> 00199 INLINE AutoDiff<D,SCAL> operator* (double x, const AutoDiff<D,SCAL> & y) throw() 00200 { 00201 AutoDiff<D,SCAL> res; 00202 res.Value() = x*y.Value(); 00203 for (int i = 0; i < D; i++) 00204 res.DValue(i) = x*y.DValue(i); 00205 return res; 00206 } 00207 00209 template<int D, typename SCAL> 00210 INLINE AutoDiff<D,SCAL> operator* (const AutoDiff<D,SCAL> & y, double x) throw() 00211 { 00212 AutoDiff<D,SCAL> res; 00213 res.Value() = x*y.Value(); 00214 for (int i = 0; i < D; i++) 00215 res.DValue(i) = x*y.DValue(i); 00216 return res; 00217 } 00218 00220 template<int D, typename SCAL> 00221 INLINE AutoDiff<D,SCAL> operator* (const AutoDiff<D,SCAL> & x, const AutoDiff<D,SCAL> & y) throw() 00222 { 00223 AutoDiff<D,SCAL> res; 00224 SCAL hx = x.Value(); 00225 SCAL hy = y.Value(); 00226 00227 res.Value() = hx*hy; 00228 for (int i = 0; i < D; i++) 00229 res.DValue(i) = hx*y.DValue(i) + hy*x.DValue(i); 00230 00231 return res; 00232 } 00233 00235 template<int D, typename SCAL> 00236 INLINE AutoDiff<D,SCAL> sqr (const AutoDiff<D,SCAL> & x) throw() 00237 { 00238 AutoDiff<D,SCAL> res; 00239 SCAL hx = x.Value(); 00240 res.Value() = hx*hx; 00241 hx *= 2; 00242 for (int i = 0; i < D; i++) 00243 res.DValue(i) = hx*x.DValue(i); 00244 return res; 00245 } 00246 00248 template<int D, typename SCAL> 00249 INLINE AutoDiff<D,SCAL> Inv (const AutoDiff<D,SCAL> & x) 00250 { 00251 AutoDiff<D,SCAL> res(1.0 / x.Value()); 00252 for (int i = 0; i < D; i++) 00253 res.DValue(i) = -x.DValue(i) / (x.Value() * x.Value()); 00254 return res; 00255 } 00256 00257 00259 template<int D, typename SCAL> 00260 INLINE AutoDiff<D,SCAL> operator/ (const AutoDiff<D,SCAL> & x, const AutoDiff<D,SCAL> & y) 00261 { 00262 return x * Inv (y); 00263 } 00264 00266 template<int D, typename SCAL> 00267 INLINE AutoDiff<D,SCAL> operator/ (const AutoDiff<D,SCAL> & x, double y) 00268 { 00269 return (1/y) * x; 00270 } 00271 00273 template<int D, typename SCAL> 00274 INLINE AutoDiff<D,SCAL> operator/ (double x, const AutoDiff<D,SCAL> & y) 00275 { 00276 return x * Inv(y); 00277 } 00278 00279 00280 00281 00282 template <int D, typename SCAL, typename SCAL2> 00283 INLINE AutoDiff<D,SCAL> & operator+= (AutoDiff<D,SCAL> & x, SCAL2 y) throw() 00284 { 00285 x.Value() += y; 00286 return x; 00287 } 00288 00289 00291 template <int D, typename SCAL> 00292 INLINE AutoDiff<D,SCAL> & operator+= (AutoDiff<D,SCAL> & x, AutoDiff<D,SCAL> y) 00293 { 00294 x.Value() += y.Value(); 00295 for (int i = 0; i < D; i++) 00296 x.DValue(i) += y.DValue(i); 00297 return x; 00298 } 00299 00301 template <int D, typename SCAL> 00302 INLINE AutoDiff<D,SCAL> & operator-= (AutoDiff<D,SCAL> & x, AutoDiff<D,SCAL> y) 00303 { 00304 x.Value() -= y.Value(); 00305 for (int i = 0; i < D; i++) 00306 x.DValue(i) -= y.DValue(i); 00307 return x; 00308 00309 } 00310 00311 template <int D, typename SCAL, typename SCAL2> 00312 INLINE AutoDiff<D,SCAL> & operator-= (AutoDiff<D,SCAL> & x, SCAL2 y) 00313 { 00314 x.Value() -= y; 00315 return x; 00316 } 00317 00319 template <int D, typename SCAL> 00320 INLINE AutoDiff<D,SCAL> & operator*= (AutoDiff<D,SCAL> & x, AutoDiff<D,SCAL> y) 00321 { 00322 for (int i = 0; i < D; i++) 00323 x.DValue(i) = x.DValue(i)*y.Value() + x.Value() * y.DValue(i); 00324 x.Value() *= y.Value(); 00325 return x; 00326 } 00327 00329 template <int D, typename SCAL, typename SCAL2> 00330 INLINE AutoDiff<D,SCAL> & operator*= (AutoDiff<D,SCAL> & x, SCAL2 y) 00331 { 00332 x.Value() *= y; 00333 for (int i = 0; i < D; i++) 00334 x.DValue(i) *= y; 00335 return x; 00336 } 00337 00339 template <int D, typename SCAL> 00340 INLINE AutoDiff<D,SCAL> & operator/= (AutoDiff<D,SCAL> & x, SCAL y) 00341 { 00342 SCAL iy = 1.0 / y; 00343 x.Value() *= iy; 00344 for (int i = 0; i < D; i++) 00345 x.DValue(i) *= iy; 00346 return x; 00347 } 00348 00349 00350 00351 00353 template <int D, typename SCAL> 00354 INLINE bool operator== (AutoDiff<D,SCAL> x, SCAL val2) 00355 { 00356 return x.Value() == val2; 00357 } 00358 00360 template <int D, typename SCAL> 00361 INLINE bool operator!= (AutoDiff<D,SCAL> x, SCAL val2) throw() 00362 { 00363 return x.Value() != val2; 00364 } 00365 00367 template <int D, typename SCAL> 00368 INLINE bool operator< (AutoDiff<D,SCAL> x, SCAL val2) throw() 00369 { 00370 return x.Value() < val2; 00371 } 00372 00374 template <int D, typename SCAL> 00375 INLINE bool operator> (AutoDiff<D,SCAL> x, SCAL val2) throw() 00376 { 00377 return x.Value() > val2; 00378 } 00379 00380 00381 00382 00383 template<int D, typename SCAL> 00384 INLINE AutoDiff<D,SCAL> fabs (const AutoDiff<D,SCAL> & x) 00385 { 00386 double abs = std::fabs (x.Value()); 00387 AutoDiff<D,SCAL> res( abs ); 00388 if (abs != 0.0) 00389 for (int i = 0; i < D; i++) 00390 res.DValue(i) = x.DValue(i) / abs; 00391 else 00392 for (int i = 0; i < D; i++) 00393 res.DValue(i) = 0.0; 00394 return res; 00395 } 00396 00397 template<int D, typename SCAL> 00398 INLINE AutoDiff<D,SCAL> sqrt (const AutoDiff<D,SCAL> & x) 00399 { 00400 AutoDiff<D,SCAL> res; 00401 res.Value() = std::sqrt(x.Value()); 00402 for (int j = 0; j < D; j++) 00403 res.DValue(j) = 0.5 / res.Value() * x.DValue(j); 00404 return res; 00405 } 00406 00407 using std::log; 00408 template <int D, typename SCAL> 00409 AutoDiff<D,SCAL> log (AutoDiff<D,SCAL> x) 00410 { 00411 AutoDiff<D,SCAL> res; 00412 res.Value() = std::log(x.Value()); 00413 for (int k = 0; k < D; k++) 00414 res.DValue(k) = x.DValue(k) / x.Value(); 00415 return res; 00416 } 00417 00418 00419 00421 00422 } 00423 00424 #endif 00425 00426