Rheolef  7.2
an efficient C++ finite element environment
Loading...
Searching...
No Matches
ad3.h
Go to the documentation of this file.
1#ifndef _RHEO_AD_POINT_H
2#define _RHEO_AD_POINT_H
23/*Class:ad3
24NAME: @code{ad3} - automatic differentiation with respect to variables in R^3
25@clindex ad3
26DESCRIPTION:
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.
31EXAMPLE:
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
41AUTHOR: Pierre.Saramito@imag.fr
42DATE: 26 september 2017
43End:
44*/
45#include "rheolef/point.h"
46
47namespace rheolef {
48
49template <class T>
50class ad3_basic {
51public:
52
53// typedefs:
54
55 typedef typename point_basic<T>::size_type size_type;
56 typedef T value_type;
57
58// allocators:
59
60 ad3_basic ();
61 ad3_basic (const ad3_basic&);
62 ad3_basic<T>& operator= (const ad3_basic<T>&);
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
85 ad3_basic<T>& operator += (const ad3_basic<T>& b);
86 ad3_basic<T>& operator -= (const ad3_basic<T>& b);
87 ad3_basic<T>& operator *= (const ad3_basic<T>& b);
88 ad3_basic<T>& operator /= (const ad3_basic<T>& b);
89
90// data:
91public:
92 T _v;
93 point_basic<T> _g;
94};
95typedef ad3_basic<Float> ad3;
96
97// output:
98template <class T>
99std::ostream& operator<< (std::ostream& os, const ad3_basic<T>& a);
100
101// get the underlying float
102template<class T> struct float_traits<ad3_basic<T> > { typedef typename float_traits<T>::type type; };
103
104// ----------------------------------------------------------------------------
105// inlined
106// ----------------------------------------------------------------------------
107template <class T>
108inline
109ad3_basic<T>::ad3_basic()
110 : _v(0),
111 _g()
112{
113}
114template <class T>
115inline
116ad3_basic<T>::ad3_basic(const ad3_basic<T>& a)
117 : _v(a._v),
118 _g(a._g)
119{
120}
121template <class T>
122inline
123ad3_basic<T>&
124ad3_basic<T>::operator= (const ad3_basic<T>& a)
125{
126 _v = a._v;
127 _g = a._g;
128 return *this;
129}
130template <class T>
131template <class U>
132inline
133ad3_basic<T>::ad3_basic(const U& v)
134 : _v(v),
135 _g()
136{
137}
138template <class T>
139inline
140point_basic<ad3_basic<T> >
141ad3_basic<T>::point (const point_basic<T>& x0)
142{
143 point_basic<ad3_basic<T> > x;
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}
150template <class T>
151inline
152const T&
153ad3_basic<T>::value () const
154{
155 return _v;
156}
157template <class T>
158inline
159const point_basic<T>&
160ad3_basic<T>::grad () const
161{
162 return _g;
163}
164template <class T>
165inline
166std::ostream&
167operator<< (std::ostream& os, const ad3_basic<T>& a)
168{
169 os << a.value();
170 return os;
171}
172// ----------------------------------------------------------------------------
173// op+-
174// ----------------------------------------------------------------------------
175template <class T>
176inline
177ad3_basic<T>
178ad3_basic<T>::operator+ () const
179{
180 ad3_basic<T> c;
181 c._v = _v;
182 c._g = _g;
183 return c;
184}
185template <class T>
186inline
187ad3_basic<T>
188ad3_basic<T>::operator- () const
189{
190 ad3_basic<T> c;
191 c._v = - _v;
192 c._g = - _g;
193 return c;
194}
195template <class T>
196inline
197ad3_basic<T>
198ad3_basic<T>::operator+ (const ad3_basic<T>& b) const
199{
200 ad3_basic<T> c;
201 c._v = _v + b._v;
202 c._g = _g + b._g;
203 return c;
204}
205template <class T>
206inline
207ad3_basic<T>
208ad3_basic<T>::operator- (const ad3_basic<T>& b) const
209{
210 ad3_basic<T> c;
211 c._v = _v - b._v;
212 c._g = _g - b._g;
213 return c;
214}
215template <class T, class U>
216inline
217typename
218std::enable_if<
219 details::is_rheolef_arithmetic<U>::value
220 ,ad3_basic<T>
221>::type
222operator+ (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}
229template <class T, class U>
230inline
231typename
232std::enable_if<
233 details::is_rheolef_arithmetic<U>::value
234 ,ad3_basic<T>
235>::type
236operator+ (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}
243template <class T, class U>
244inline
245typename
246std::enable_if<
247 details::is_rheolef_arithmetic<U>::value
248 ,ad3_basic<T>
249>::type
250operator- (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}
257template <class T, class U>
258inline
259typename
260std::enable_if<
261 details::is_rheolef_arithmetic<U>::value
262 ,ad3_basic<T>
263>::type
264operator- (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}
271template <class T>
272inline
273ad3_basic<T>&
274ad3_basic<T>::operator+= (const ad3_basic<T>& b)
275{
276 *this = (*this) + b;
277 return *this;
278}
279template <class T, class U>
280inline
281typename
282std::enable_if<
283 details::is_rheolef_arithmetic<U>::value
284 ,ad3_basic<T>&
285>::type
286operator+= (ad3_basic<T>& a, const U& b)
287{
288 a = a + b;
289 return a;
290}
291template <class T>
292inline
293ad3_basic<T>&
294ad3_basic<T>::operator-= (const ad3_basic<T>& b)
295{
296 *this = (*this) - b;
297 return *this;
298}
299template <class T, class U>
300inline
301typename
302std::enable_if<
303 details::is_rheolef_arithmetic<U>::value
304 ,ad3_basic<T>&
305>::type
306operator-= (ad3_basic<T>& a, const U& b)
307{
308 a = a - b;
309 return a;
310}
311// ----------------------------------------------------------------------------
312// op*
313// ----------------------------------------------------------------------------
314template <class T>
315inline
316ad3_basic<T>
317ad3_basic<T>::operator* (const ad3_basic<T>& b) const
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}
324template <class T, class U>
325inline
326typename
327std::enable_if<
328 details::is_rheolef_arithmetic<U>::value
329 ,ad3_basic<T>
330>::type
331operator* (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}
338template <class T, class U>
339inline
340typename
341std::enable_if<
342 details::is_rheolef_arithmetic<U>::value
343 ,ad3_basic<T>
344>::type
345operator* (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}
352template <class T>
353inline
354ad3_basic<T>&
355ad3_basic<T>::operator*= (const ad3_basic<T>& b)
356{
357 *this = (*this) * b;
358 return *this;
359}
360template <class T, class U>
361inline
362typename
363std::enable_if<
364 details::is_rheolef_arithmetic<U>::value
365 ,ad3_basic<T>&
366>::type
367operator*= (ad3_basic<T>& a, const U& b)
368{
369 a = a * b;
370 return a;
371}
372// ----------------------------------------------------------------------------
373// op/
374// ----------------------------------------------------------------------------
375template <class T>
376inline
377ad3_basic<T>
378ad3_basic<T>::operator/ (const ad3_basic<T>& b) const
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}
385template <class T, class U>
386inline
387typename
388std::enable_if<
389 details::is_rheolef_arithmetic<U>::value
390 ,ad3_basic<T>
391>::type
392operator/ (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}
399template <class T, class U>
400inline
401typename
402std::enable_if<
403 details::is_rheolef_arithmetic<U>::value
404 ,ad3_basic<T>
405>::type
406operator/ (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}
413template <class T>
414inline
415ad3_basic<T>&
416ad3_basic<T>::operator/= (const ad3_basic<T>& b)
417{
418 *this = (*this) / b;
419 return *this;
420}
421template <class T, class U>
422inline
423typename
424std::enable_if<
425 details::is_rheolef_arithmetic<U>::value
426 ,ad3_basic<T>&
427>::type
428operator/= (ad3_basic<T>& a, const U& b)
429{
430 a = a / b;
431 return a;
432}
433// -----------------------------------------------------------------------------
434// comparators
435// -----------------------------------------------------------------------------
436#define _RHEOLEF_ad3_comparator(OP) \
437template <class T> \
438inline \
439bool \
440operator OP (const ad3_basic<T> &a, const ad3_basic<T> &b) \
441{ \
442 return a._v OP b._v; \
443} \
444template <class T, class U> \
445inline \
446typename \
447std::enable_if< \
448 details::is_rheolef_arithmetic<U>::value \
449 ,bool \
450>::type \
451operator OP (const ad3_basic<T> &a, const U &b) \
452{ \
453 return a._v OP b; \
454} \
455template <class T, class U> \
456inline \
457typename \
458std::enable_if< \
459 details::is_rheolef_arithmetic<U>::value \
460 ,bool \
461>::type \
462operator 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// ----------------------------------------------------------------------------
479template <class T>
480inline
481ad3_basic<T>
482sqr (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}
489template <class T>
490inline
491ad3_basic<T>
492fabs (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// ----------------------------------------------------------------------------
503template <class T>
504INLINE2 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
512template <class T>
513INLINE2 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
521template <class T>
522INLINE2 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
531template <class T>
532INLINE2 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
541template <class T>
542INLINE2 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
551template <class T>
552INLINE2 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
561template <class T>
562INLINE2 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
571template <class T>
572INLINE2 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
581template <class T>
582INLINE2 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// ----------------------------------------------------------------------------
593template <class T>
594INLINE2 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
603template <class T>
604INLINE2 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}
612template <class T>
613INLINE2 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}
625template <class T>
626INLINE2 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
rheolef::ad3_basic
Definition ad3.h:50
rheolef::ad3_basic::operator*=
ad3_basic< T > & operator*=(const ad3_basic< T > &b)
Definition ad3.h:355
rheolef::ad3_basic::_g
point_basic< T > _g
Definition ad3.h:93
rheolef::ad3_basic::value_type
T value_type
Definition ad3.h:56
rheolef::ad3_basic::operator-=
ad3_basic< T > & operator-=(const ad3_basic< T > &b)
Definition ad3.h:294
rheolef::ad3_basic::operator/
ad3_basic< T > operator/(const ad3_basic< T > &b) const
Definition ad3.h:378
rheolef::ad3_basic::operator*
ad3_basic< T > operator*(const ad3_basic< T > &b) const
Definition ad3.h:317
rheolef::ad3_basic::operator-
ad3_basic< T > operator-() const
Definition ad3.h:188
rheolef::ad3_basic::grad
const point_basic< T > & grad() const
Definition ad3.h:160
rheolef::ad3_basic::ad3_basic
ad3_basic()
Definition ad3.h:109
rheolef::ad3_basic::point
static point_basic< ad3_basic< T > > point(const point_basic< T > &x0)
Definition ad3.h:141
rheolef::ad3_basic::operator/=
ad3_basic< T > & operator/=(const ad3_basic< T > &b)
Definition ad3.h:416
rheolef::ad3_basic::size_type
point_basic< T >::size_type size_type
Definition ad3.h:55
rheolef::ad3_basic::value
const T & value() const
Definition ad3.h:153
rheolef::ad3_basic::operator+=
ad3_basic< T > & operator+=(const ad3_basic< T > &b)
Definition ad3.h:274
rheolef::ad3_basic::ad3_basic
ad3_basic(const ad3_basic &)
rheolef::ad3_basic::operator+
ad3_basic< T > operator+() const
Definition ad3.h:178
rheolef::ad3_basic::operator=
ad3_basic< T > & operator=(const ad3_basic< T > &)
Definition ad3.h:124
rheolef::ad3_basic::_v
T _v
Definition ad3.h:92
rheolef::point_basic
Definition point.h:86
rheolef::point_basic::size_type
size_t size_type
Definition point.h:91
rheolef::ad3
ad3_basic< Float > ad3
Definition ad3.h:95
rheolef::type
rheolef::std::enable_if< details::is_field_expr_v2_nonlinear_arg< Expr1 >::value &&! details::is_field_expr_v2_constant< Expr1 >::value &&details::is_field_expr_v2_nonlinear_arg< Expr2 >::value &&! details::is_field_expr_v2_constant< Expr2 >::value, details::field_expr_v2_nonlinear_node_binary< details::function_traits< Function >::functor_type, details::field_expr_v2_nonlinear_terminal_wrapper_traits< Expr1 >::type, details::field_expr_v2_nonlinear_terminal_wrapper_traits< Expr2 >::type > > type
T
Expr1::float_type T
Definition field_expr.h:230
mkgeo_ball.b
int b
Definition mkgeo_ball.sh:152
mkgeo_ball.c
int c
Definition mkgeo_ball.sh:153
mkgeo_ball.tmp
str tmp
Definition mkgeo_ball.sh:380
mkgeo_ball.a
int a
Definition mkgeo_ball.sh:151
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::operator-
csr< T, sequential > operator-(const csr< T, sequential > &a)
Definition csr.h:447
rheolef::exp
tensor_basic< T > exp(const tensor_basic< T > &a, size_t d)
Definition tensor-exp.cc:92
rheolef::operator<<
std::ostream & operator<<(std::ostream &os, const catchmark &m)
Definition catchmark.h:99
rheolef::operator/
dia< T, M > operator/(const T &lambda, const dia< T, M > &d)
Definition dia.h:145
rheolef::pow
space_mult_list< T, M > pow(const space_basic< T, M > &X, size_t n)
Definition space_mult.h:120
rheolef::operator*=
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
rheolef::operator/=
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
rheolef::operator-=
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
rheolef::operator+=
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
rheolef::operator+
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
rheolef::operator*
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition csr.h:437
rheolef::float_traits< ad3_basic< T > >::type
float_traits< T >::type type
Definition ad3.h:102
rheolef::float_traits
helper for std::complex<T>: get basic T type
Definition Float.h:93
rheolef::float_traits::type
T type
Definition Float.h:94