NGSolve  5.3
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   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