NGSolve
4.9
|
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 ALWAYS_INLINE AutoDiff () throw() { }; 00034 // { val = 0; for (int i = 0; i < D; i++) dval[i] = 0; } // ! 00035 00037 ALWAYS_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 ALWAYS_INLINE AutoDiff (SCAL aval) throw() 00046 { 00047 val = aval; 00048 for (int i = 0; i < D; i++) 00049 dval[i] = 0; 00050 } 00051 00053 ALWAYS_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 00062 ALWAYS_INLINE AutoDiff & operator= (SCAL aval) throw() 00063 { 00064 val = aval; 00065 for (int i = 0; i < D; i++) 00066 dval[i] = 0; 00067 return *this; 00068 } 00069 00071 SCAL Value() const throw() { return val; } 00072 00074 SCAL DValue (int i) const throw() { return dval[i]; } 00075 00077 SCAL & Value() throw() { return val; } 00078 00080 SCAL & DValue (int i) throw() { return dval[i]; } 00081 00082 00083 ALWAYS_INLINE AutoDiff<D,SCAL> & operator+= (const SCAL & y) throw() 00084 { 00085 val += y; 00086 return *this; 00087 } 00088 00089 00091 ALWAYS_INLINE AutoDiff<D,SCAL> & operator+= (const AutoDiff<D,SCAL> & y) throw() 00092 { 00093 val += y.val; 00094 for (int i = 0; i < D; i++) 00095 dval[i] += y.dval[i]; 00096 return *this; 00097 } 00098 00100 ALWAYS_INLINE AutoDiff<D,SCAL> & operator-= (const AutoDiff<D,SCAL> & y) throw() 00101 { 00102 val -= y.val; 00103 for (int i = 0; i < D; i++) 00104 dval[i] -= y.dval[i]; 00105 return *this; 00106 00107 } 00108 00110 ALWAYS_INLINE AutoDiff<D,SCAL> & operator*= (const AutoDiff<D,SCAL> & y) throw() 00111 { 00112 for (int i = 0; i < D; i++) 00113 { 00114 // dval[i] *= y.val; 00115 // dval[i] += val * y.dval[i]; 00116 dval[i] = dval[i] * y.val + val * y.dval[i]; 00117 } 00118 val *= y.val; 00119 return *this; 00120 } 00121 00123 ALWAYS_INLINE AutoDiff<D,SCAL> & operator*= (const SCAL & y) throw() 00124 { 00125 val *= y; 00126 for (int i = 0; i < D; i++) 00127 dval[i] *= y; 00128 return *this; 00129 } 00130 00132 ALWAYS_INLINE AutoDiff<D,SCAL> & operator/= (const SCAL & y) throw() 00133 { 00134 SCAL iy = 1.0 / y; 00135 val *= iy; 00136 for (int i = 0; i < D; i++) 00137 dval[i] *= iy; 00138 return *this; 00139 } 00140 00142 ALWAYS_INLINE bool operator== (SCAL val2) throw() 00143 { 00144 return val == val2; 00145 } 00146 00148 bool operator!= (SCAL val2) throw() 00149 { 00150 return val != val2; 00151 } 00152 00154 bool operator< (SCAL val2) throw() 00155 { 00156 return val < val2; 00157 } 00158 00160 bool operator> (SCAL val2) throw() 00161 { 00162 return val > val2; 00163 } 00164 }; 00165 00166 00168 00170 template<int D, typename SCAL> 00171 inline ostream & operator<< (ostream & ost, const AutoDiff<D,SCAL> & x) 00172 { 00173 ost << x.Value() << ", D = "; 00174 for (int i = 0; i < D; i++) 00175 ost << x.DValue(i) << " "; 00176 return ost; 00177 } 00178 00180 template<int D, typename SCAL> 00181 ALWAYS_INLINE inline AutoDiff<D,SCAL> operator+ (const AutoDiff<D,SCAL> & x, const AutoDiff<D,SCAL> & y) throw() 00182 { 00183 AutoDiff<D,SCAL> res; 00184 res.Value () = x.Value()+y.Value(); 00185 // AutoDiff<D,SCAL> res(x.Value()+y.Value()); 00186 for (int i = 0; i < D; i++) 00187 res.DValue(i) = x.DValue(i) + y.DValue(i); 00188 return res; 00189 } 00190 00191 00193 template<int D, typename SCAL> 00194 ALWAYS_INLINE inline AutoDiff<D,SCAL> operator- (const AutoDiff<D,SCAL> & x, const AutoDiff<D,SCAL> & y) throw() 00195 { 00196 AutoDiff<D,SCAL> res; 00197 res.Value() = x.Value()-y.Value(); 00198 // AutoDiff<D,SCAL> res (x.Value()-y.Value()); 00199 for (int i = 0; i < D; i++) 00200 res.DValue(i) = x.DValue(i) - y.DValue(i); 00201 return res; 00202 } 00203 00205 template<int D, typename SCAL> 00206 ALWAYS_INLINE inline AutoDiff<D,SCAL> operator+ (double x, const AutoDiff<D,SCAL> & y) throw() 00207 { 00208 AutoDiff<D,SCAL> res; 00209 res.Value() = x+y.Value(); 00210 for (int i = 0; i < D; i++) 00211 res.DValue(i) = y.DValue(i); 00212 return res; 00213 } 00214 00216 template<int D, typename SCAL> 00217 ALWAYS_INLINE inline AutoDiff<D,SCAL> operator+ (const AutoDiff<D,SCAL> & y, double x) throw() 00218 { 00219 AutoDiff<D,SCAL> res; 00220 res.Value() = x+y.Value(); 00221 for (int i = 0; i < D; i++) 00222 res.DValue(i) = y.DValue(i); 00223 return res; 00224 } 00225 00226 00228 template<int D, typename SCAL> 00229 ALWAYS_INLINE inline AutoDiff<D,SCAL> operator- (const AutoDiff<D,SCAL> & x) throw() 00230 { 00231 AutoDiff<D,SCAL> res; 00232 res.Value() = -x.Value(); 00233 for (int i = 0; i < D; i++) 00234 res.DValue(i) = -x.DValue(i); 00235 return res; 00236 } 00237 00239 template<int D, typename SCAL> 00240 ALWAYS_INLINE inline AutoDiff<D,SCAL> operator- (const AutoDiff<D,SCAL> & x, double y) throw() 00241 { 00242 AutoDiff<D,SCAL> res; 00243 res.Value() = x.Value()-y; 00244 for (int i = 0; i < D; i++) 00245 res.DValue(i) = x.DValue(i); 00246 return res; 00247 } 00248 00250 template<int D, typename SCAL> 00251 ALWAYS_INLINE inline AutoDiff<D,SCAL> operator- (double x, const AutoDiff<D,SCAL> & y) throw() 00252 { 00253 AutoDiff<D,SCAL> res; 00254 res.Value() = x-y.Value(); 00255 for (int i = 0; i < D; i++) 00256 res.DValue(i) = -y.DValue(i); 00257 return res; 00258 } 00259 00260 00262 template<int D, typename SCAL> 00263 ALWAYS_INLINE inline AutoDiff<D,SCAL> operator* (double x, const AutoDiff<D,SCAL> & y) throw() 00264 { 00265 AutoDiff<D,SCAL> res; 00266 res.Value() = x*y.Value(); 00267 for (int i = 0; i < D; i++) 00268 res.DValue(i) = x*y.DValue(i); 00269 return res; 00270 } 00271 00273 template<int D, typename SCAL> 00274 ALWAYS_INLINE inline AutoDiff<D,SCAL> operator* (const AutoDiff<D,SCAL> & y, double x) throw() 00275 { 00276 AutoDiff<D,SCAL> res; 00277 res.Value() = x*y.Value(); 00278 for (int i = 0; i < D; i++) 00279 res.DValue(i) = x*y.DValue(i); 00280 return res; 00281 } 00282 00284 template<int D, typename SCAL> 00285 ALWAYS_INLINE inline AutoDiff<D,SCAL> operator* (const AutoDiff<D,SCAL> & x, const AutoDiff<D,SCAL> & y) throw() 00286 { 00287 AutoDiff<D,SCAL> res; 00288 SCAL hx = x.Value(); 00289 SCAL hy = y.Value(); 00290 00291 res.Value() = hx*hy; 00292 for (int i = 0; i < D; i++) 00293 res.DValue(i) = hx*y.DValue(i) + hy*x.DValue(i); 00294 00295 return res; 00296 } 00297 00299 template<int D, typename SCAL> 00300 inline AutoDiff<D,SCAL> sqr (const AutoDiff<D,SCAL> & x) throw() 00301 { 00302 AutoDiff<D,SCAL> res; 00303 SCAL hx = x.Value(); 00304 res.Value() = hx*hx; 00305 hx *= 2; 00306 for (int i = 0; i < D; i++) 00307 res.DValue(i) = hx*x.DValue(i); 00308 return res; 00309 } 00310 00312 template<int D, typename SCAL> 00313 inline AutoDiff<D,SCAL> Inv (const AutoDiff<D,SCAL> & x) 00314 { 00315 AutoDiff<D,SCAL> res(1.0 / x.Value()); 00316 for (int i = 0; i < D; i++) 00317 res.DValue(i) = -x.DValue(i) / (x.Value() * x.Value()); 00318 return res; 00319 } 00320 00321 00323 template<int D, typename SCAL> 00324 inline AutoDiff<D,SCAL> operator/ (const AutoDiff<D,SCAL> & x, const AutoDiff<D,SCAL> & y) 00325 { 00326 return x * Inv (y); 00327 } 00328 00330 template<int D, typename SCAL> 00331 inline AutoDiff<D,SCAL> operator/ (const AutoDiff<D,SCAL> & x, double y) 00332 { 00333 return (1/y) * x; 00334 } 00335 00337 template<int D, typename SCAL> 00338 inline AutoDiff<D,SCAL> operator/ (double x, const AutoDiff<D,SCAL> & y) 00339 { 00340 return x * Inv(y); 00341 } 00342 00343 00344 00345 00346 template<int D, typename SCAL> 00347 inline AutoDiff<D,SCAL> fabs (const AutoDiff<D,SCAL> & x) 00348 { 00349 double abs = std::fabs (x.Value()); 00350 AutoDiff<D,SCAL> res( abs ); 00351 if (abs != 0.0) 00352 for (int i = 0; i < D; i++) 00353 res.DValue(i) = x.DValue(i) / abs; 00354 else 00355 for (int i = 0; i < D; i++) 00356 res.DValue(i) = 0.0; 00357 return res; 00358 } 00359 00360 template<int D, typename SCAL> 00361 inline AutoDiff<D,SCAL> sqrt (const AutoDiff<D,SCAL> & x) 00362 { 00363 AutoDiff<D,SCAL> res; 00364 res.Value() = std::sqrt(x.Value()); 00365 for (int j = 0; j < D; j++) 00366 res.DValue(j) = 0.5 / res.Value() * x.DValue(j); 00367 return res; 00368 } 00369 00370 00371 00373 00374 } 00375 00376 #endif