Rheolef  7.1
an efficient C++ finite element environment
ad3.h
Go to the documentation of this file.
1 #ifndef _RHEO_AD_POINT_H
2 #define _RHEO_AD_POINT_H
23 /*Class:ad3
24 NAME: @code{ad3} - automatic differentiation with respect to variables in R^3
25 @clindex ad3
26 DESCRIPTION:
27  @noindent
28  The @code{ad3} class defines a forward automatic differentiation
29  in R^3 for computing automatically the gradient of a function from R^3 to R.
30  The implementation uses a simple forward automatic differentiation method.
31 EXAMPLE:
32  @example
33  template<class T>
34  T f (const point_basic<T>& x) @{ return sqr(x[0]) + x[0]*x[1] + 1/x[2]; @}
35  ...
36  point_basic<ad3> x_ad = ad3::point (x);
37  ad3 y_ad = f(x_ad);
38  cout << "f ="<<y_ad <<endl
39  << "gradq(f)="<<y_ad.grad() <<endl;
40  @end example
41 AUTHOR: Pierre.Saramito@imag.fr
42 DATE: 26 september 2017
43 End:
44 */
45 #include "rheolef/point.h"
46 
47 namespace rheolef {
48 
49 template <class T>
50 class ad3_basic {
51 public:
52 
53 // typedefs:
54 
56  typedef T value_type;
57 
58 // allocators:
59 
60  ad3_basic ();
61  ad3_basic (const ad3_basic&);
63  template <class U>
64  ad3_basic(const U& x);
65 
66 // helper:
67 
68  static point_basic<ad3_basic<T> >
69  point (const point_basic<T>& x0);
70 
71 // accessors:
72 
73  const T& value() const;
74  const point_basic<T>& grad() const;
75 
76 // algebra:
77 
78  ad3_basic<T> operator+ () const;
79  ad3_basic<T> operator- () const;
80  ad3_basic<T> operator+ (const ad3_basic<T>& b) const;
81  ad3_basic<T> operator- (const ad3_basic<T>& b) const;
82  ad3_basic<T> operator* (const ad3_basic<T>& b) const;
83  ad3_basic<T> operator/ (const ad3_basic<T>& b) const;
84 
89 
90 // data:
91 public:
92  T _v;
94 };
96 
97 // output:
98 template <class T>
99 std::ostream& operator<< (std::ostream& os, const ad3_basic<T>& a);
100 
101 // get the underlying float
102 template<class T> struct float_traits<ad3_basic<T> > { typedef typename float_traits<T>::type type; };
103 
104 // ----------------------------------------------------------------------------
105 // inlined
106 // ----------------------------------------------------------------------------
107 template <class T>
108 inline
110  : _v(0),
111  _g()
112 {
113 }
114 template <class T>
115 inline
117  : _v(a._v),
118  _g(a._g)
119 {
120 }
121 template <class T>
122 inline
125 {
126  _v = a._v;
127  _g = a._g;
128  return *this;
129 }
130 template <class T>
131 template <class U>
132 inline
134  : _v(v),
135  _g()
136 {
137 }
138 template <class T>
139 inline
142 {
144  for (size_type i = 0; i < 3; ++i) {
145  x[i]._v = x0[i];
146  x[i]._g[i] = 1;
147  }
148  return x;
149 }
150 template <class T>
151 inline
152 const T&
154 {
155  return _v;
156 }
157 template <class T>
158 inline
159 const point_basic<T>&
161 {
162  return _g;
163 }
164 template <class T>
165 inline
166 std::ostream&
167 operator<< (std::ostream& os, const ad3_basic<T>& a)
168 {
169  os << a.value();
170  return os;
171 }
172 // ----------------------------------------------------------------------------
173 // op+-
174 // ----------------------------------------------------------------------------
175 template <class T>
176 inline
177 ad3_basic<T>
179 {
180  ad3_basic<T> c;
181  c._v = _v;
182  c._g = _g;
183  return c;
184 }
185 template <class T>
186 inline
189 {
190  ad3_basic<T> c;
191  c._v = - _v;
192  c._g = - _g;
193  return c;
194 }
195 template <class T>
196 inline
199 {
200  ad3_basic<T> c;
201  c._v = _v + b._v;
202  c._g = _g + b._g;
203  return c;
204 }
205 template <class T>
206 inline
209 {
210  ad3_basic<T> c;
211  c._v = _v - b._v;
212  c._g = _g - b._g;
213  return c;
214 }
215 template <class T, class U>
216 inline
217 typename
218 std::enable_if<
219  details::is_rheolef_arithmetic<U>::value
220  ,ad3_basic<T>
221 >::type
222 operator+ (const U& a, const ad3_basic<T>& b)
223 {
224  ad3_basic<T> c;
225  c._v = a + b._v;
226  c._g = b._g;
227  return c;
228 }
229 template <class T, class U>
230 inline
231 typename
232 std::enable_if<
233  details::is_rheolef_arithmetic<U>::value
234  ,ad3_basic<T>
235 >::type
236 operator+ (const ad3_basic<T>& a, const U& b)
237 {
238  ad3_basic<T> c;
239  c._v = a._v + b;
240  c._g = a._g;
241  return c;
242 }
243 template <class T, class U>
244 inline
245 typename
246 std::enable_if<
247  details::is_rheolef_arithmetic<U>::value
248  ,ad3_basic<T>
249 >::type
250 operator- (const U& a, const ad3_basic<T>& b)
251 {
252  ad3_basic<T> c;
253  c._v = a - b._v;
254  c._g = - b._g;
255  return c;
256 }
257 template <class T, class U>
258 inline
259 typename
260 std::enable_if<
261  details::is_rheolef_arithmetic<U>::value
262  ,ad3_basic<T>
263 >::type
264 operator- (const ad3_basic<T>& a, const U& b)
265 {
266  ad3_basic<T> c;
267  c._v = a._v - b;
268  c._g = a._g;
269  return c;
270 }
271 template <class T>
272 inline
273 ad3_basic<T>&
275 {
276  *this = (*this) + b;
277  return *this;
278 }
279 template <class T, class U>
280 inline
281 typename
282 std::enable_if<
283  details::is_rheolef_arithmetic<U>::value
284  ,ad3_basic<T>&
285 >::type
287 {
288  a = a + b;
289  return a;
290 }
291 template <class T>
292 inline
293 ad3_basic<T>&
295 {
296  *this = (*this) - b;
297  return *this;
298 }
299 template <class T, class U>
300 inline
301 typename
302 std::enable_if<
303  details::is_rheolef_arithmetic<U>::value
304  ,ad3_basic<T>&
305 >::type
307 {
308  a = a - b;
309  return a;
310 }
311 // ----------------------------------------------------------------------------
312 // op*
313 // ----------------------------------------------------------------------------
314 template <class T>
315 inline
316 ad3_basic<T>
318 {
319  ad3_basic<T> c;
320  c._v = _v*b._v;
321  c._g = _v*b._g + _g*b._v;
322  return c;
323 }
324 template <class T, class U>
325 inline
326 typename
327 std::enable_if<
328  details::is_rheolef_arithmetic<U>::value
329  ,ad3_basic<T>
330 >::type
331 operator* (const U& a, const ad3_basic<T>& b)
332 {
333  ad3_basic<T> c;
334  c._v = a*b._v;
335  c._g = a*b._g;
336  return c;
337 }
338 template <class T, class U>
339 inline
340 typename
341 std::enable_if<
342  details::is_rheolef_arithmetic<U>::value
343  ,ad3_basic<T>
344 >::type
345 operator* (const ad3_basic<T>& a, const U& b)
346 {
347  ad3_basic<T> c;
348  c._v = a._v*b;
349  c._g = a._g*b;
350  return c;
351 }
352 template <class T>
353 inline
354 ad3_basic<T>&
356 {
357  *this = (*this) * b;
358  return *this;
359 }
360 template <class T, class U>
361 inline
362 typename
363 std::enable_if<
364  details::is_rheolef_arithmetic<U>::value
365  ,ad3_basic<T>&
366 >::type
368 {
369  a = a * b;
370  return a;
371 }
372 // ----------------------------------------------------------------------------
373 // op/
374 // ----------------------------------------------------------------------------
375 template <class T>
376 inline
377 ad3_basic<T>
379 {
380  ad3_basic<T> c;
381  c._v = _v/b._v;
382  c._g = _g/b._v - _v*b._g/b._v;
383  return c;
384 }
385 template <class T, class U>
386 inline
387 typename
388 std::enable_if<
389  details::is_rheolef_arithmetic<U>::value
390  ,ad3_basic<T>
391 >::type
392 operator/ (const U& a, const ad3_basic<T>& b)
393 {
394  ad3_basic<T> c;
395  c._v = a/b._v;
396  c._g = - a*b._g/b._v;
397  return c;
398 }
399 template <class T, class U>
400 inline
401 typename
402 std::enable_if<
403  details::is_rheolef_arithmetic<U>::value
404  ,ad3_basic<T>
405 >::type
406 operator/ (const ad3_basic<T>& a, const U& b)
407 {
408  ad3_basic<T> c;
409  c._v = a._v/b;
410  c._g = a._g/b;
411  return c;
412 }
413 template <class T>
414 inline
415 ad3_basic<T>&
417 {
418  *this = (*this) / b;
419  return *this;
420 }
421 template <class T, class U>
422 inline
423 typename
424 std::enable_if<
425  details::is_rheolef_arithmetic<U>::value
426  ,ad3_basic<T>&
427 >::type
429 {
430  a = a / b;
431  return a;
432 }
433 // -----------------------------------------------------------------------------
434 // comparators
435 // -----------------------------------------------------------------------------
436 #define _RHEOLEF_ad3_comparator(OP) \
437 template <class T> \
438 inline \
439 bool \
440 operator OP (const ad3_basic<T> &a, const ad3_basic<T> &b) \
441 { \
442  return a._v OP b._v; \
443 } \
444 template <class T, class U> \
445 inline \
446 typename \
447 std::enable_if< \
448  details::is_rheolef_arithmetic<U>::value \
449  ,bool \
450 >::type \
451 operator OP (const ad3_basic<T> &a, const U &b) \
452 { \
453  return a._v OP b; \
454 } \
455 template <class T, class U> \
456 inline \
457 typename \
458 std::enable_if< \
459  details::is_rheolef_arithmetic<U>::value \
460  ,bool \
461 >::type \
462 operator OP (const U &a, const ad3_basic<T> &b) \
463 { \
464  return a OP b._v; \
465 }
466 
467 _RHEOLEF_ad3_comparator(==)
468 _RHEOLEF_ad3_comparator(!=)
469 _RHEOLEF_ad3_comparator(>)
470 _RHEOLEF_ad3_comparator(>=)
471 _RHEOLEF_ad3_comparator(<)
472 _RHEOLEF_ad3_comparator(<=)
473 #undef _RHEOLEF_ad3_comparator
474 // ----------------------------------------------------------------------------
475 // std math functions: sqr, pow, exp...
476 // ----------------------------------------------------------------------------
477 // unary functions
478 // ----------------------------------------------------------------------------
479 template <class T>
480 inline
481 ad3_basic<T>
482 sqr (const ad3_basic<T>& a)
483 {
484  ad3_basic<T> c;
485  c._v = sqr(a._v);
486  c._g = (2*a._v)*a._g;
487  return c;
488 }
489 template <class T>
490 inline
491 ad3_basic<T>
492 fabs (const ad3_basic<T>& a)
493 {
494  ad3_basic<T> c;
495  c._v = fabs(a._v);
496  c._g = (a._v >= 0 ? 1 : -1)*a._g;
497  return c;
498 }
499 #ifdef TODO
500 // ----------------------------------------------------------------------------
501 // unary functions (cont.)
502 // ----------------------------------------------------------------------------
503 template <class T>
504 INLINE2 ad3_basic<T> exp (const ad3_basic<T>& a)
505 {
506  ad3_basic<T> c(exp(a.v));
507  c.touchg(a.gsize);
508  for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*c.v;
509  return c;
510 }
511 
512 template <class T>
513 INLINE2 ad3_basic<T> log (const ad3_basic<T>& a)
514 {
515  ad3_basic<T> c(log(a.v));
516  c.touchg(a.gsize);
517  for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]/a.v;
518  return c;
519 }
520 
521 template <class T>
522 INLINE2 ad3_basic<T> sqrt (const ad3_basic<T>& a)
523 {
524  ad3_basic<T> c(sqrt(a.v));
525  T tmp(c.v*FADBAD_TWO);
526  c.touchg(a.gsize);
527  for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]/tmp;
528  return c;
529 }
530 
531 template <class T>
532 INLINE2 ad3_basic<T> sin (const ad3_basic<T>& a)
533 {
534  ad3_basic<T> c(sin(a.v));
535  T tmp(cos(a.v));
536  c.touchg(a.gsize);
537  for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*tmp;
538  return c;
539 }
540 
541 template <class T>
542 INLINE2 ad3_basic<T> cos (const ad3_basic<T>& a)
543 {
544  ad3_basic<T> c(cos(a.v));
545  T tmp(-sin(a.v));
546  c.touchg(a.gsize);
547  for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*tmp;
548  return c;
549 }
550 
551 template <class T>
552 INLINE2 ad3_basic<T> tan (const ad3_basic<T>& a)
553 {
554  ad3_basic<T> c(tan(a.v));
555  c.touchg(a.gsize);
556  T tmp(FADBAD_ONE+_sqr(c.v));
557  for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*tmp;
558  return c;
559 }
560 
561 template <class T>
562 INLINE2 ad3_basic<T> asin (const ad3_basic<T>& a)
563 {
564  ad3_basic<T> c(asin(a.v));
565  c.touchg(a.gsize);
566  T tmp(FADBAD_ONE/sqrt(FADBAD_ONE-_sqr(a.v)));
567  for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*tmp;
568  return c;
569 }
570 
571 template <class T>
572 INLINE2 ad3_basic<T> acos (const ad3_basic<T>& a)
573 {
574  ad3_basic<T> c(acos(a.v));
575  c.touchg(a.gsize);
576  T tmp(-FADBAD_ONE/sqrt(FADBAD_ONE-_sqr(a.v)));
577  for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*tmp;
578  return c;
579 }
580 
581 template <class T>
582 INLINE2 ad3_basic<T> atan (const ad3_basic<T>& a)
583 {
584  ad3_basic<T> c(atan(a.v));
585  c.touchg(a.gsize);
586  T tmp(FADBAD_ONE/(FADBAD_ONE+_sqr(a.v)));
587  for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*tmp;
588  return c;
589 }
590 // ----------------------------------------------------------------------------
591 // binary functions
592 // ----------------------------------------------------------------------------
593 template <class T>
594 INLINE2 ad3_basic<T> pow (const BaseType& a, const ad3_basic<T>& b)
595 {
596  ad3_basic<T> c(pow(a,b.v));
597  T tmp(c.v*log(a));
598  c.touchg(b.gsize);
599  for (int i=0;i<b.gsize;i++) c.g[i]=tmp*b.g[i];
600  return c;
601 }
602 
603 template <class T>
604 INLINE2 ad3_basic<T> pow (const ad3_basic<T>& a,const BaseType& b)
605 {
606  ad3_basic<T> c(pow(a.v,b));
607  T tmp(b*pow(a.v,b-FADBAD_ONE));
608  c.touchg(a.gsize);
609  for (int i=0;i<a.gsize;i++) c.g[i]=tmp*a.g[i];
610  return c;
611 }
612 template <class T>
613 INLINE2 ad3_basic<T> pow (const ad3_basic<T>& a, const ad3_basic<T>& b)
614 {
615  if (a.gsize==0) return pow1(a.v,b);
616  if (b.gsize==0) return pow2(a,b.v);
617  ad3_basic<T> c(pow(a.v,b.v));
618  USER_ASSERT(a.gsize==b.gsize,"derivative vectors not of same size in pow");
619  T tmp(b.v*pow(a.v,b.v-FADBAD_ONE)),tmp1(c.v*log(a.v));
620  c.touchg(a.gsize);
621  for (int i=0;i<a.gsize;i++)
622  c.g[i]=tmp*a.g[i]+tmp1*b.g[i];
623  return c;
624 }
625 template <class T>
626 INLINE2 ad3_basic<T> pow (const ad3_basic<T>& a,int b)
627 {
628  ad3_basic<T> c(pow(a.v,b));
629  c.touchg(a.gsize);
630  T tmp(b*pow(a.v,b-1));
631  for (int i=0;i<a.gsize;i++) c.g[i]=a.g[i]*tmp;
632  return c;
633 }
634 #endif // TODO
635 
636 } // namespace rheolef
637 #endif // _RHEO_AD_POINT_H
ad3_basic< T > & operator*=(const ad3_basic< T > &b)
Definition: ad3.h:355
point_basic< T > _g
Definition: ad3.h:93
ad3_basic< T > & operator-=(const ad3_basic< T > &b)
Definition: ad3.h:294
ad3_basic< T > operator/(const ad3_basic< T > &b) const
Definition: ad3.h:378
ad3_basic< T > operator*(const ad3_basic< T > &b) const
Definition: ad3.h:317
ad3_basic< T > operator-() const
Definition: ad3.h:188
const point_basic< T > & grad() const
Definition: ad3.h:160
static point_basic< ad3_basic< T > > point(const point_basic< T > &x0)
Definition: ad3.h:141
ad3_basic< T > & operator/=(const ad3_basic< T > &b)
Definition: ad3.h:416
point_basic< T >::size_type size_type
Definition: ad3.h:55
const T & value() const
Definition: ad3.h:153
ad3_basic< T > & operator+=(const ad3_basic< T > &b)
Definition: ad3.h:274
ad3_basic(const ad3_basic &)
ad3_basic< T > operator+() const
Definition: ad3.h:178
ad3_basic< T > & operator=(const ad3_basic< T > &)
Definition: ad3.h:124
size_t size_type
Definition: point.h:92
rheolef::std type
ad3_basic< Float > ad3
Definition: ad3.h:95
Expr1::float_type T
Definition: field_expr.h:261
This file is part of Rheolef.
std::enable_if< details::is_rheolef_arithmetic< U >::value,ad3_basic< T > & >::type operator*=(ad3_basic< T > &a, const U &b)
Definition: ad3.h:367
dia< T, M > operator/(const T &lambda, const dia< T, M > &d)
Definition: dia.h:145
std::enable_if< details::is_rheolef_arithmetic< U >::value,ad3_basic< T > & >::type operator-=(ad3_basic< T > &a, const U &b)
Definition: ad3.h:306
std::enable_if< details::is_rheolef_arithmetic< U >::value,ad3_basic< T >>::type operator+(const U &a, const ad3_basic< T > &b)
Definition: ad3.h:222
csr< T, sequential > operator-(const csr< T, sequential > &a)
Definition: csr.h:447
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition: csr.h:437
std::enable_if< details::is_rheolef_arithmetic< U >::value,ad3_basic< T > & >::type operator/=(ad3_basic< T > &a, const U &b)
Definition: ad3.h:428
std::ostream & operator<<(std::ostream &os, const catchmark &m)
Definition: catchmark.h:99
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition: space_mult.h:120
std::enable_if< details::is_rheolef_arithmetic< U >::value,ad3_basic< T > & >::type operator+=(ad3_basic< T > &a, const U &b)
Definition: ad3.h:286
tensor_basic< T > exp(const tensor_basic< T > &a, size_t d)
Definition: tensor-exp.cc:92
float_traits< T >::type type
Definition: ad3.h:102
helper for std::complex<T>: get basic T type
Definition: Float.h:93