NGSolve  4.9
ngstd/autodiff.hpp
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