NGSolve
5.3
|
00001 #ifndef FILE_AUTODIFFDIFF 00002 #define FILE_AUTODIFFDIFF 00003 00004 /**************************************************************************/ 00005 /* File: autodiffdiff.hpp */ 00006 /* Author: Joachim Schoeberl */ 00007 /* Date: 13. June. 05 */ 00008 /**************************************************************************/ 00009 00010 namespace ngstd 00011 { 00012 00013 // Automatic second differentiation datatype 00014 00015 00021 template <int D> 00022 class AutoDiffDiff 00023 { 00024 double val; 00025 double dval[D]; 00026 double ddval[D*D]; 00027 public: 00028 00029 typedef AutoDiffDiff<D> TELEM; 00030 typedef double TSCAL; 00031 00032 00034 AutoDiffDiff () throw() { ; } 00035 00037 AutoDiffDiff (const AutoDiffDiff & ad2) throw() 00038 { 00039 val = ad2.val; 00040 for (int i = 0; i < D; i++) 00041 dval[i] = ad2.dval[i]; 00042 for (int i = 0; i < D*D; i++) 00043 ddval[i] = ad2.ddval[i]; 00044 } 00045 00047 AutoDiffDiff (double aval) throw() 00048 { 00049 val = aval; 00050 for (int i = 0; i < D; i++) 00051 dval[i] = 0; 00052 for (int i = 0; i < D*D; i++) 00053 ddval[i] = 0; 00054 } 00055 00057 AutoDiffDiff (const AutoDiff<D> & ad2) throw() 00058 { 00059 val = ad2.Value(); 00060 for (int i = 0; i < D; i++) 00061 dval[i] = ad2.DValue(i); 00062 for (int i = 0; i < D*D; i++) 00063 ddval[i] = 0; 00064 } 00065 00067 AutoDiffDiff (double aval, int diffindex) throw() 00068 { 00069 val = aval; 00070 for (int i = 0; i < D; i++) 00071 dval[i] = 0; 00072 for (int i = 0; i < D*D; i++) 00073 ddval[i] = 0; 00074 dval[diffindex] = 1; 00075 } 00076 00078 AutoDiffDiff & operator= (double aval) throw() 00079 { 00080 val = aval; 00081 for (int i = 0; i < D; i++) 00082 dval[i] = 0; 00083 for (int i = 0; i < D*D; i++) 00084 ddval[i] = 0; 00085 return *this; 00086 } 00087 00089 double Value() const throw() { return val; } 00090 00092 double DValue (int i) const throw() { return dval[i]; } 00093 00095 double DDValue (int i) const throw() { return ddval[i]; } 00096 00098 double DDValue (int i, int j) const throw() { return ddval[i*D+j]; } 00099 00101 double & Value() throw() { return val; } 00102 00104 double & DValue (int i) throw() { return dval[i]; } 00105 00107 double & DDValue (int i) throw() { return ddval[i]; } 00108 00110 double & DDValue (int i, int j) throw() { return ddval[i*D+j]; } 00111 00113 AutoDiffDiff<D> & operator+= (const AutoDiffDiff<D> & y) throw() 00114 { 00115 val += y.val; 00116 for (int i = 0; i < D; i++) 00117 dval[i] += y.dval[i]; 00118 for (int i = 0; i < D*D; i++) 00119 ddval[i] += y.ddval[i]; 00120 return *this; 00121 } 00122 00124 AutoDiffDiff<D> & operator-= (const AutoDiffDiff<D> & y) throw() 00125 { 00126 val -= y.val; 00127 for (int i = 0; i < D; i++) 00128 dval[i] -= y.dval[i]; 00129 for (int i = 0; i < D*D; i++) 00130 ddval[i] -= y.ddval[i]; 00131 return *this; 00132 } 00133 00135 AutoDiffDiff<D> & operator*= (const AutoDiffDiff<D> & y) throw() 00136 { 00137 for (int i = 0; i < D*D; i++) 00138 ddval[i] = val * y.ddval[i] + y.val * ddval[i]; 00139 00140 for (int i = 0; i < D; i++) 00141 for (int j = 0; j < D; j++) 00142 ddval[i*D+j] += dval[i] * y.dval[j] + dval[j] * y.dval[i]; 00143 00144 for (int i = 0; i < D; i++) 00145 { 00146 dval[i] *= y.val; 00147 dval[i] += val * y.dval[i]; 00148 } 00149 val *= y.val; 00150 return *this; 00151 } 00152 00154 AutoDiffDiff<D> & operator*= (const double & y) throw() 00155 { 00156 for ( int i = 0; i < D*D; i++ ) 00157 ddval[i] *= y; 00158 for (int i = 0; i < D; i++) 00159 dval[i] *= y; 00160 val *= y; 00161 return *this; 00162 } 00163 00165 AutoDiffDiff<D> & operator/= (const double & y) throw() 00166 { 00167 double iy = 1.0 / y; 00168 for ( int i = 0; i < D*D; i++ ) 00169 ddval[i] *= iy; 00170 for (int i = 0; i < D; i++) 00171 dval[i] *= iy; 00172 val *= iy; 00173 return *this; 00174 } 00175 00177 bool operator== (double val2) throw() 00178 { 00179 return val == val2; 00180 } 00181 00183 bool operator!= (double val2) throw() 00184 { 00185 return val != val2; 00186 } 00187 00189 bool operator< (double val2) throw() 00190 { 00191 return val < val2; 00192 } 00193 00195 bool operator> (double val2) throw() 00196 { 00197 return val > val2; 00198 } 00199 }; 00200 00201 00203 00205 template<int D> 00206 inline ostream & operator<< (ostream & ost, const AutoDiffDiff<D> & x) 00207 { 00208 ost << x.Value() << ", D = "; 00209 for (int i = 0; i < D; i++) 00210 ost << x.DValue(i) << " "; 00211 ost << ", DD = "; 00212 for (int i = 0; i < D*D; i++) 00213 ost << x.DDValue(i) << " "; 00214 return ost; 00215 } 00216 00218 template<int D> 00219 inline AutoDiffDiff<D> operator+ (const AutoDiffDiff<D> & x, const AutoDiffDiff<D> & y) throw() 00220 { 00221 AutoDiffDiff<D> res; 00222 res.Value () = x.Value()+y.Value(); 00223 for (int i = 0; i < D; i++) 00224 res.DValue(i) = x.DValue(i) + y.DValue(i); 00225 for (int i = 0; i < D*D; i++) 00226 res.DDValue(i) = x.DDValue(i) + y.DDValue(i); 00227 return res; 00228 } 00229 00230 00232 template<int D> 00233 inline AutoDiffDiff<D> operator- (const AutoDiffDiff<D> & x, const AutoDiffDiff<D> & y) throw() 00234 { 00235 AutoDiffDiff<D> res; 00236 res.Value() = x.Value()-y.Value(); 00237 for (int i = 0; i < D; i++) 00238 res.DValue(i) = x.DValue(i) - y.DValue(i); 00239 for (int i = 0; i < D*D; i++) 00240 res.DDValue(i) = x.DDValue(i) - y.DDValue(i); 00241 return res; 00242 } 00243 00244 00246 template<int D> 00247 inline AutoDiffDiff<D> operator+ (double x, const AutoDiffDiff<D> & y) throw() 00248 { 00249 AutoDiffDiff<D> res; 00250 res.Value() = x+y.Value(); 00251 for (int i = 0; i < D; i++) 00252 res.DValue(i) = y.DValue(i); 00253 for (int i = 0; i < D*D; i++) 00254 res.DDValue(i) = y.DDValue(i); 00255 return res; 00256 } 00257 00259 template<int D> 00260 inline AutoDiffDiff<D> operator+ (const AutoDiffDiff<D> & y, double x) throw() 00261 { 00262 AutoDiffDiff<D> res; 00263 res.Value() = x+y.Value(); 00264 for (int i = 0; i < D; i++) 00265 res.DValue(i) = y.DValue(i); 00266 for (int i = 0; i < D*D; i++) 00267 res.DDValue(i) = y.DDValue(i); 00268 return res; 00269 } 00270 00271 00273 template<int D> 00274 inline AutoDiffDiff<D> operator- (const AutoDiffDiff<D> & x) throw() 00275 { 00276 AutoDiffDiff<D> res; 00277 res.Value() = -x.Value(); 00278 for (int i = 0; i < D; i++) 00279 res.DValue(i) = -x.DValue(i); 00280 for (int i = 0; i < D*D; i++) 00281 res.DDValue(i) = -x.DDValue(i); 00282 return res; 00283 } 00284 00286 template<int D> 00287 inline AutoDiffDiff<D> operator- (const AutoDiffDiff<D> & x, double y) throw() 00288 { 00289 AutoDiffDiff<D> res; 00290 res.Value() = x.Value()-y; 00291 for (int i = 0; i < D; i++) 00292 res.DValue(i) = x.DValue(i); 00293 for (int i = 0; i < D*D; i++) 00294 res.DDValue(i) = x.DDValue(i); 00295 return res; 00296 } 00297 00299 template<int D> 00300 inline AutoDiffDiff<D> operator- (double x, const AutoDiffDiff<D> & y) throw() 00301 { 00302 AutoDiffDiff<D> res; 00303 res.Value() = x-y.Value(); 00304 for (int i = 0; i < D; i++) 00305 res.DValue(i) = -y.DValue(i); 00306 for (int i = 0; i < D*D; i++) 00307 res.DDValue(i) = -y.DDValue(i); 00308 return res; 00309 } 00310 00311 00313 template<int D> 00314 inline AutoDiffDiff<D> operator* (double x, const AutoDiffDiff<D> & y) throw() 00315 { 00316 AutoDiffDiff<D> res; 00317 res.Value() = x*y.Value(); 00318 for (int i = 0; i < D; i++) 00319 res.DValue(i) = x*y.DValue(i); 00320 for (int i = 0; i < D*D; i++) 00321 res.DDValue(i) = x*y.DDValue(i); 00322 return res; 00323 } 00324 00326 template<int D> 00327 inline AutoDiffDiff<D> operator* (const AutoDiffDiff<D> & y, double x) throw() 00328 { 00329 AutoDiffDiff<D> res; 00330 res.Value() = x*y.Value(); 00331 for (int i = 0; i < D; i++) 00332 res.DValue(i) = x*y.DValue(i); 00333 for (int i = 0; i < D*D; i++) 00334 res.DDValue(i) = x*y.DDValue(i); 00335 return res; 00336 } 00337 00339 template<int D> 00340 inline AutoDiffDiff<D> operator* (const AutoDiffDiff<D> & x, const AutoDiffDiff<D> & y) throw() 00341 { 00342 AutoDiffDiff<D> res; 00343 double hx = x.Value(); 00344 double hy = y.Value(); 00345 00346 res.Value() = hx*hy; 00347 for (int i = 0; i < D; i++) 00348 res.DValue(i) = hx*y.DValue(i) + hy*x.DValue(i); 00349 00350 for (int i = 0; i < D; i++) 00351 for (int j = 0; j < D; j++) 00352 res.DDValue(i,j) = hx * y.DDValue(i,j) + hy * x.DDValue(i,j) 00353 + x.DValue(i) * y.DValue(j) + x.DValue(j) * y.DValue(i); 00354 00355 return res; 00356 } 00357 00358 00359 00360 template<int D> 00361 inline AutoDiffDiff<D> Inv (const AutoDiffDiff<D> & x) 00362 { 00363 AutoDiffDiff<D> res(1.0 / x.Value()); 00364 for (int i = 0; i < D; i++) 00365 res.DValue(i) = -x.DValue(i) / (x.Value() * x.Value()); 00366 cout << "ADD::Inv not implemented" << endl; 00367 return res; 00368 } 00369 00370 00371 template<int D> 00372 inline AutoDiffDiff<D> operator/ (const AutoDiffDiff<D> & x, const AutoDiffDiff<D> & y) 00373 { 00374 return x * Inv (y); 00375 } 00376 00377 template<int D> 00378 inline AutoDiffDiff<D> operator/ (const AutoDiffDiff<D> & x, double y) 00379 { 00380 return (1/y) * x; 00381 } 00382 00383 template<int D> 00384 inline AutoDiffDiff<D> operator/ (double x, const AutoDiffDiff<D> & y) 00385 { 00386 return x * Inv(y); 00387 } 00388 00389 00390 template<int D> 00391 inline AutoDiffDiff<D> sqrt (const AutoDiffDiff<D> & x) 00392 { 00393 AutoDiffDiff<D> res; 00394 res.Value() = std::sqrt(x.Value()); 00395 for (int j = 0; j < D; j++) 00396 res.DValue(j) = 0.5 / res.Value() * x.DValue(j); 00397 00398 00399 for (int i = 0; i < D; i++) 00400 for (int j = 0; j < D; j++) 00401 res.DDValue(i,j) = 0.5/res.Value() * x.DDValue(i,j) - 0.25 / (x.Value()*res.Value()) * x.DValue(i) * x.DValue(j); 00402 00403 return res; 00404 } 00405 00406 00407 00408 00409 00411 00412 } 00413 00414 00415 #endif