NGSolve  5.3
ngstd/autodiffdiff.hpp
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