NGSolve 5.3
autodiff.hpp
1#ifndef FILE_AUTODIFF
2#define FILE_AUTODIFF
3
4/**************************************************************************/
5/* File: autodiff.hpp */
6/* Author: Joachim Schoeberl */
7/* Date: 24. Oct. 02 */
8/**************************************************************************/
9
10namespace ngstd
11{
12 using ngcore::IfPos;
13
14// Automatic differentiation datatype
15
16 template <int D, typename SCAL = double> class AutoDiffRec;
17
18
24template <int D, typename SCAL = double>
26{
27 SCAL val;
28 SCAL dval[D?D:1];
29public:
30
32 typedef SCAL TSCAL;
33
34
36 // INLINE AutoDiffVec () throw() { };
37 AutoDiffVec() = default;
38 // { val = 0; for (int i = 0; i < D; i++) dval[i] = 0; } // !
39
41 AutoDiffVec (const AutoDiffVec & ad2) = default;
42 /*
43 INLINE AutoDiffVec (const AutoDiffVec & ad2) throw()
44 {
45 val = ad2.val;
46 for (int i = 0; i < D; i++)
47 dval[i] = ad2.dval[i];
48 }
49 */
51 INLINE AutoDiffVec (SCAL aval) throw()
52 {
53 val = aval;
54 for (int i = 0; i < D; i++)
55 dval[i] = 0;
56 }
57
59 INLINE AutoDiffVec (SCAL aval, int diffindex) throw()
60 {
61 val = aval;
62 for (int i = 0; i < D; i++)
63 dval[i] = 0;
64 dval[diffindex] = 1;
65 }
66
67 INLINE AutoDiffVec (SCAL aval, const SCAL * grad)
68 {
69 val = aval;
70 LoadGradient (grad);
71 }
72
74 INLINE AutoDiffVec & operator= (SCAL aval) throw()
75 {
76 val = aval;
77 for (int i = 0; i < D; i++)
78 dval[i] = 0;
79 return *this;
80 }
81
82 AutoDiffVec & operator= (const AutoDiffVec & ad2) = default;
83
85 INLINE SCAL Value() const throw() { return val; }
86
88 INLINE SCAL DValue (int i) const throw() { return dval[i]; }
89
91 INLINE void StoreGradient (SCAL * p) const
92 {
93 for (int i = 0; i < D; i++)
94 p[i] = dval[i];
95 }
96
97 INLINE void LoadGradient (const SCAL * p)
98 {
99 for (int i = 0; i < D; i++)
100 dval[i] = p[i];
101 }
102
104 INLINE SCAL & Value() throw() { return val; }
105
107 INLINE SCAL & DValue (int i) throw() { return dval[i]; }
108};
109
110
112
114template<int D, typename SCAL>
115inline ostream & operator<< (ostream & ost, const AutoDiffVec<D,SCAL> & x)
116{
117 ost << x.Value() << ", D = ";
118 for (int i = 0; i < D; i++)
119 ost << x.DValue(i) << " ";
120 return ost;
121}
122
124template<int D, typename SCAL>
125INLINE AutoDiffVec<D,SCAL> operator+ (const AutoDiffVec<D,SCAL> & x, const AutoDiffVec<D,SCAL> & y) throw()
126{
128 res.Value () = x.Value()+y.Value();
129 // AutoDiffVec<D,SCAL> res(x.Value()+y.Value());
130 for (int i = 0; i < D; i++)
131 res.DValue(i) = x.DValue(i) + y.DValue(i);
132 return res;
133}
134
135
137template<int D, typename SCAL>
138INLINE AutoDiffVec<D,SCAL> operator- (const AutoDiffVec<D,SCAL> & x, const AutoDiffVec<D,SCAL> & y) throw()
139{
141 res.Value() = x.Value()-y.Value();
142 // AutoDiffVec<D,SCAL> res (x.Value()-y.Value());
143 for (int i = 0; i < D; i++)
144 res.DValue(i) = x.DValue(i) - y.DValue(i);
145 return res;
146}
147
149 template<int D, typename SCAL, typename SCAL2,
150 typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
151INLINE AutoDiffVec<D,SCAL> operator+ (SCAL2 x, const AutoDiffVec<D,SCAL> & y) throw()
152{
154 res.Value() = x+y.Value();
155 for (int i = 0; i < D; i++)
156 res.DValue(i) = y.DValue(i);
157 return res;
158}
159
161 template<int D, typename SCAL, typename SCAL2,
162 typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
163INLINE AutoDiffVec<D,SCAL> operator+ (const AutoDiffVec<D,SCAL> & y, SCAL2 x) throw()
164{
166 res.Value() = x+y.Value();
167 for (int i = 0; i < D; i++)
168 res.DValue(i) = y.DValue(i);
169 return res;
170}
171
172
174template<int D, typename SCAL>
175INLINE AutoDiffVec<D,SCAL> operator- (const AutoDiffVec<D,SCAL> & x) throw()
176{
178 res.Value() = -x.Value();
179 for (int i = 0; i < D; i++)
180 res.DValue(i) = -x.DValue(i);
181 return res;
182}
183
185 template<int D, typename SCAL, typename SCAL2,
186 typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
187INLINE AutoDiffVec<D,SCAL> operator- (const AutoDiffVec<D,SCAL> & x, SCAL2 y) throw()
188{
190 res.Value() = x.Value()-y;
191 for (int i = 0; i < D; i++)
192 res.DValue(i) = x.DValue(i);
193 return res;
194}
195
197 template<int D, typename SCAL, typename SCAL2,
198 typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
199INLINE AutoDiffVec<D,SCAL> operator- (SCAL2 x, const AutoDiffVec<D,SCAL> & y) throw()
200{
201 AutoDiffVec<D,SCAL> res;
202 res.Value() = x-y.Value();
203 for (int i = 0; i < D; i++)
204 res.DValue(i) = -y.DValue(i);
205 return res;
206}
207
208
210 template<int D, typename SCAL, typename SCAL2,
211 typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
213{
215 res.Value() = x*y.Value();
216 for (int i = 0; i < D; i++)
217 res.DValue(i) = x*y.DValue(i);
218 return res;
219}
220
222 template<int D, typename SCAL, typename SCAL2,
223 typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
224
226{
228 res.Value() = x*y.Value();
229 for (int i = 0; i < D; i++)
230 res.DValue(i) = x*y.DValue(i);
231 return res;
232}
233
235template<int D, typename SCAL>
237{
239 SCAL hx = x.Value();
240 SCAL hy = y.Value();
241
242 res.Value() = hx*hy;
243 for (int i = 0; i < D; i++)
244 res.DValue(i) = hx*y.DValue(i) + hy*x.DValue(i);
245
246 return res;
247}
248
250using ngcore::sqr;
251template<int D, typename SCAL>
252INLINE AutoDiffVec<D,SCAL> sqr (const AutoDiffVec<D,SCAL> & x) throw()
253{
254 AutoDiffVec<D,SCAL> res;
255 SCAL hx = x.Value();
256 res.Value() = hx*hx;
257 hx *= 2;
258 for (int i = 0; i < D; i++)
259 res.DValue(i) = hx*x.DValue(i);
260 return res;
261}
262
264template<int D, typename SCAL>
266{
268 for (int i = 0; i < D; i++)
269 res.DValue(i) = -x.DValue(i) / (x.Value() * x.Value());
270 return res;
271}
272
273
275template<int D, typename SCAL>
277{
278 return x * Inv (y);
279}
280
282template<int D, typename SCAL, typename SCAL2,
283 typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
285{
286 return (1.0/y) * x;
287}
288
290template<int D, typename SCAL, typename SCAL2,
291 typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
293 {
294 return x * Inv(y);
295 }
296
297
298
299
300 template <int D, typename SCAL, typename SCAL2>
301 INLINE AutoDiffVec<D,SCAL> & operator+= (AutoDiffVec<D,SCAL> & x, SCAL2 y) throw()
302 {
303 x.Value() += y;
304 return x;
305 }
306
307
309 template <int D, typename SCAL>
310 INLINE AutoDiffVec<D,SCAL> & operator+= (AutoDiffVec<D,SCAL> & x, AutoDiffVec<D,SCAL> y)
311 {
312 x.Value() += y.Value();
313 for (int i = 0; i < D; i++)
314 x.DValue(i) += y.DValue(i);
315 return x;
316 }
317
319 template <int D, typename SCAL>
320 INLINE AutoDiffVec<D,SCAL> & operator-= (AutoDiffVec<D,SCAL> & x, AutoDiffVec<D,SCAL> y)
321 {
322 x.Value() -= y.Value();
323 for (int i = 0; i < D; i++)
324 x.DValue(i) -= y.DValue(i);
325 return x;
326
327 }
328
329 template <int D, typename SCAL, typename SCAL2>
330 INLINE AutoDiffVec<D,SCAL> & operator-= (AutoDiffVec<D,SCAL> & x, SCAL2 y)
331 {
332 x.Value() -= y;
333 return x;
334 }
335
337 template <int D, typename SCAL>
338 INLINE AutoDiffVec<D,SCAL> & operator*= (AutoDiffVec<D,SCAL> & x, AutoDiffVec<D,SCAL> y)
339 {
340 for (int i = 0; i < D; i++)
341 x.DValue(i) = x.DValue(i)*y.Value() + x.Value() * y.DValue(i);
342 x.Value() *= y.Value();
343 return x;
344 }
345
347 template <int D, typename SCAL, typename SCAL2>
348 INLINE AutoDiffVec<D,SCAL> & operator*= (AutoDiffVec<D,SCAL> & x, SCAL2 y)
349 {
350 x.Value() *= y;
351 for (int i = 0; i < D; i++)
352 x.DValue(i) *= y;
353 return x;
354 }
355
357 template <int D, typename SCAL>
358 INLINE AutoDiffVec<D,SCAL> & operator/= (AutoDiffVec<D,SCAL> & x, SCAL y)
359 {
360 SCAL iy = 1.0 / y;
361 x.Value() *= iy;
362 for (int i = 0; i < D; i++)
363 x.DValue(i) *= iy;
364 return x;
365 }
366
367
368
369
371 template <int D, typename SCAL>
372 INLINE bool operator== (AutoDiffVec<D,SCAL> x, SCAL val2)
373 {
374 return x.Value() == val2;
375 }
376
378 template <int D, typename SCAL>
379 INLINE bool operator!= (AutoDiffVec<D,SCAL> x, SCAL val2) throw()
380 {
381 return x.Value() != val2;
382 }
383
385 template <int D, typename SCAL>
386 INLINE bool operator< (AutoDiffVec<D,SCAL> x, SCAL val2) throw()
387 {
388 return x.Value() < val2;
389 }
390
392 template <int D, typename SCAL>
393 INLINE bool operator> (AutoDiffVec<D,SCAL> x, SCAL val2) throw()
394 {
395 return x.Value() > val2;
396 }
397
398
399
400
401template<int D, typename SCAL>
402INLINE AutoDiffVec<D,SCAL> fabs (const AutoDiffVec<D,SCAL> & x)
403{
404 double abs = fabs (x.Value());
405 AutoDiffVec<D,SCAL> res( abs );
406 if (abs != 0.0)
407 for (int i = 0; i < D; i++)
408 res.DValue(i) = x.Value()*x.DValue(i) / abs;
409 else
410 for (int i = 0; i < D; i++)
411 res.DValue(i) = 0.0;
412 return res;
413}
414
415using std::sqrt;
416template<int D, typename SCAL>
417INLINE AutoDiffVec<D,SCAL> sqrt (const AutoDiffVec<D,SCAL> & x)
418{
419 AutoDiffVec<D,SCAL> res;
420 res.Value() = sqrt(x.Value());
421 for (int j = 0; j < D; j++)
422 res.DValue(j) = 0.5 / res.Value() * x.DValue(j);
423 return res;
424}
425
426using std::log;
427template <int D, typename SCAL>
428AutoDiffVec<D,SCAL> log (AutoDiffVec<D,SCAL> x)
429{
430 AutoDiffVec<D,SCAL> res;
431 res.Value() = log(x.Value());
432 for (int k = 0; k < D; k++)
433 res.DValue(k) = x.DValue(k) / x.Value();
434 return res;
435}
436
437using std::exp;
438template <int D, typename SCAL>
439INLINE AutoDiffVec<D,SCAL> exp (AutoDiffVec<D,SCAL> x)
440{
441 AutoDiffVec<D,SCAL> res;
442 res.Value() = exp(x.Value());
443 for (int k = 0; k < D; k++)
444 res.DValue(k) = x.DValue(k) * res.Value();
445 return res;
446}
447
448using std::pow;
449template <int D, typename SCAL>
450INLINE AutoDiffVec<D,SCAL> pow (AutoDiffVec<D,SCAL> x, AutoDiffVec<D,SCAL> y )
451{
452 return exp(log(x)*y);
453}
454
455using std::sin;
456/*
457template <int D, typename SCAL>
458INLINE AutoDiffVec<D,SCAL> sin (AutoDiffVec<D,SCAL> x)
459{
460 AutoDiffVec<D,SCAL> res;
461 res.Value() = sin(x.Value());
462 SCAL c = cos(x.Value());
463 for (int k = 0; k < D; k++)
464 res.DValue(k) = x.DValue(k) * c;
465 return res;
466}
467*/
468
469template <int D, typename SCAL>
470INLINE AutoDiffVec<D,SCAL> sin (AutoDiffVec<D,SCAL> x)
471{
472 return sin(AutoDiffRec<D,SCAL>(x));
473}
474
475using std::cos;
476/*
477template <int D, typename SCAL>
478INLINE AutoDiffVec<D,SCAL> cos (AutoDiffVec<D,SCAL> x)
479{
480 AutoDiffVec<D,SCAL> res;
481 res.Value() = cos(x.Value());
482 SCAL ms = -sin(x.Value());
483 for (int k = 0; k < D; k++)
484 res.DValue(k) = x.DValue(k) * ms;
485 return res;
486}
487*/
488template <int D, typename SCAL>
489INLINE AutoDiffVec<D,SCAL> cos (AutoDiffVec<D,SCAL> x)
490{
491 return cos(AutoDiffRec<D,SCAL>(x));
492}
493
494using std::tan;
495template <int D, typename SCAL>
496INLINE AutoDiffVec<D,SCAL> tan (AutoDiffVec<D,SCAL> x)
497{ return sin(x) / cos(x); }
498
499using std::sinh;
500template <int D, typename SCAL>
501INLINE AutoDiffVec<D,SCAL> sinh (AutoDiffVec<D,SCAL> x)
502{
503 AutoDiffVec<D,SCAL> res;
504 res.Value() = sinh(x.Value());
505 SCAL ch = cosh(x.Value());
506 for (int k = 0; k < D; k++)
507 res.DValue(k) = x.DValue(k) * ch;
508 return res;
509}
510
511using std::cosh;
512template <int D, typename SCAL>
513INLINE AutoDiffVec<D,SCAL> cosh (AutoDiffVec<D,SCAL> x)
514{
515 AutoDiffVec<D,SCAL> res;
516 res.Value() = cosh(x.Value());
517 SCAL sh = sinh(x.Value());
518 for (int k = 0; k < D; k++)
519 res.DValue(k) = x.DValue(k) * sh;
520 return res;
521}
522
523using std::erf;
524template <int D, typename SCAL>
525INLINE AutoDiffVec<D,SCAL> erf (AutoDiffVec<D,SCAL> x)
526{
527 return erf(AutoDiffRec<D,SCAL>(x));
528}
529
530using std::floor;
531template<int D, typename SCAL>
532INLINE AutoDiffVec<D,SCAL> floor (const AutoDiffVec<D,SCAL> & x)
533{
534 AutoDiffVec<D,SCAL> res;
535 res.Value() = floor(x.Value());
536 for (int j = 0; j < D; j++)
537 res.DValue(j) = 0.0;
538 return res;
539}
540
541using std::ceil;
542template<int D, typename SCAL>
543INLINE AutoDiffVec<D,SCAL> ceil (const AutoDiffVec<D,SCAL> & x)
544{
545 AutoDiffVec<D,SCAL> res;
546 res.Value() = ceil(x.Value());
547 for (int j = 0; j < D; j++)
548 res.DValue(j) = 0.0;
549 return res;
550}
551
552
553using std::atan;
554/*
555template <int D, typename SCAL>
556INLINE AutoDiffVec<D,SCAL> atan (AutoDiffVec<D,SCAL> x)
557{
558 AutoDiffVec<D,SCAL> res;
559 SCAL a = atan(x.Value());
560 res.Value() = a;
561 for (int k = 0; k < D; k++)
562 res.DValue(k) = x.DValue(k)/(1+x.Value()*x.Value()) ;
563 return res;
564}
565*/
566template <int D, typename SCAL>
567AutoDiffVec<D,SCAL> atan (AutoDiffVec<D,SCAL> x)
568{
569 return atan (AutoDiffRec<D,SCAL> (x));
570}
571
572using std::atan2;
573template <int D, typename SCAL>
574INLINE AutoDiffVec<D,SCAL> atan2 (AutoDiffVec<D,SCAL> x, AutoDiffVec<D,SCAL> y)
575{
576 AutoDiffVec<D,SCAL> res;
577 SCAL a = atan2(x.Value(), y.Value());
578 res.Value() = a;
579 for (int k = 0; k < D; k++)
580 res.DValue(k) = (x.Value()*y.DValue(k)-y.Value()*x.DValue(k))/(y.Value()*y.Value()+x.Value()*x.Value());
581 return res;
582}
583
584
585using std::acos;
586template <int D, typename SCAL>
587INLINE AutoDiffVec<D,SCAL> acos (AutoDiffVec<D,SCAL> x)
588{
589 AutoDiffVec<D,SCAL> res;
590 SCAL a = acos(x.Value());
591 res.Value() = a;
592 SCAL da = -1 / sqrt(1-x.Value()*x.Value());
593 for (int k = 0; k < D; k++)
594 res.DValue(k) = x.DValue(k)*da;
595 return res;
596}
597
598
599using std::asin;
600template <int D, typename SCAL>
601INLINE AutoDiffVec<D,SCAL> asin (AutoDiffVec<D,SCAL> x)
602{
603 AutoDiffVec<D,SCAL> res;
604 SCAL a = asin(x.Value());
605 res.Value() = a;
606 SCAL da = 1 / sqrt(1-x.Value()*x.Value());
607 for (int k = 0; k < D; k++)
608 res.DValue(k) = x.DValue(k)*da;
609 return res;
610}
611
612
613
614
615 template <int D, typename SCAL, typename TB, typename TC>
616 auto IfPos (AutoDiffVec<D,SCAL> a, TB b, TC c) // -> decltype(IfPos (a.Value(), b, c))
617 {
618 return IfPos (a.Value(), b, c);
619 }
620
621 template <int D, typename SCAL>
622 INLINE AutoDiffVec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffVec<D,SCAL> b, AutoDiffVec<D,SCAL> c)
623 {
624 AutoDiffVec<D,SCAL> res;
625 res.Value() = IfPos (a, b.Value(), c.Value());
626 for (int j = 0; j < D; j++)
627 res.DValue(j) = IfPos (a, b.DValue(j), c.DValue(j));
628 return res;
629 }
630
631 template <int D, typename SCAL, typename TC>
632 INLINE AutoDiffVec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffVec<D,SCAL> b, TC c)
633 {
634 return IfPos (a, b, AutoDiffVec<D,SCAL> (c));
635 }
636
638
639
640
641 template <int D, typename SCAL>
643 {
644 AutoDiffRec<D-1, SCAL> rec;
645 SCAL last;
646
647 public:
648 INLINE AutoDiffRec () = default;
649 INLINE AutoDiffRec (const AutoDiffRec &) = default;
651 INLINE AutoDiffRec & operator= (const AutoDiffRec &) = default;
652
653 INLINE AutoDiffRec (SCAL aval) : rec(aval), last(0.0) { ; }
654 INLINE AutoDiffRec (SCAL aval, int diffindex) : rec(aval, diffindex), last((diffindex==D-1) ? 1.0 : 0.0) { ; }
656 : rec(aval, grad), last(grad[D-1]) { }
657
659 {
660 Value() = ad.Value();
661 for (int i = 0; i < D; i++)
662 DValue(i) = ad.DValue(i);
663 }
664
665 INLINE AutoDiffRec & operator= (SCAL aval) { rec = aval; last = 0.0; return *this; }
666 INLINE SCAL Value() const { return rec.Value(); }
667 INLINE SCAL DValue(int i) const { return (i == D-1) ? last : rec.DValue(i); }
668 INLINE SCAL & Value() { return rec.Value(); }
669 INLINE SCAL & DValue(int i) { return (i == D-1) ? last : rec.DValue(i); }
670 INLINE auto Rec() const { return rec; }
671 INLINE auto Last() const { return last; }
672 INLINE auto & Rec() { return rec; }
673 INLINE auto & Last() { return last; }
674 INLINE operator AutoDiffVec<D,SCAL> () const
675 {
676 AutoDiffVec<D,SCAL> res(Value());
677 for (int i = 0; i < D; i++)
678 res.DValue(i) = DValue(i);
679 return res;
680 }
681 };
682
683 template<int D, typename SCAL>
685 {
687 }
688
689 template <typename SCAL>
691 {
692 SCAL val;
693 public:
694 INLINE AutoDiffRec () = default;
695 INLINE AutoDiffRec (const AutoDiffRec &) = default;
696 INLINE AutoDiffRec (SCAL _val) : val(_val) { ; }
697 INLINE AutoDiffRec (SCAL _val, SCAL /* _dummylast */) : val(_val) { ; }
698 INLINE AutoDiffRec (SCAL aval, const SCAL * /* grad */)
699 : val(aval) { }
700
701 INLINE AutoDiffRec & operator= (const AutoDiffRec &) = default;
702 INLINE AutoDiffRec & operator= (SCAL aval) { val = aval; return *this; }
703
704 INLINE SCAL Value() const { return val; }
705 INLINE SCAL DValue(int /* i */) const { return SCAL(0); }
706 INLINE SCAL & Value() { return val; }
707 // SCAL & DValue(int i) { return val; }
708 INLINE auto Rec() const { return val; }
709 INLINE auto Last() const { return SCAL(0); }
710 INLINE auto & Rec() { return val; }
711 INLINE auto & Last() { return val; }
712 INLINE operator AutoDiffVec<0,SCAL> () const { return AutoDiffVec<0,SCAL>(); }
713 };
714
715
716 template <typename SCAL>
718 {
719 SCAL val;
720 SCAL last;
721 public:
722 INLINE AutoDiffRec () = default;
723 INLINE AutoDiffRec (const AutoDiffRec &) = default;
724 INLINE AutoDiffRec (SCAL _val) : val(_val), last(0.0) { ; }
725 INLINE AutoDiffRec (SCAL _val, SCAL _last) : val(_val), last(_last) { ; }
726 INLINE AutoDiffRec (SCAL aval, int diffindex) : val(aval), last((diffindex==0) ? 1.0 : 0.0) { ; }
728 : val(aval), last(grad[0]) { }
729
731 {
732 Value() = ad.Value();
733 DValue(0) = ad.DValue(0);
734 }
735
736 INLINE AutoDiffRec & operator= (const AutoDiffRec &) = default;
737 INLINE AutoDiffRec & operator= (SCAL aval) { val = aval; last = 0.0; return *this; }
738
739 INLINE SCAL Value() const { return val; }
740 INLINE SCAL DValue(int /* i */) const { return last; }
741 INLINE SCAL & Value() { return val; }
742 INLINE SCAL & DValue(int /* i */) { return last; }
743 INLINE auto Rec() const { return val; }
744 INLINE auto Last() const { return last; }
745 INLINE auto & Rec() { return val; }
746 INLINE auto & Last() { return last; }
747
748 INLINE operator AutoDiffVec<1,SCAL> () const
749 {
750 AutoDiffVec<1,SCAL> res(Value());
751 res.DValue(0) = DValue(0);
752 return res;
753 }
754 };
755
756 template <int D, typename SCAL, typename SCAL2,
757 typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
759 {
760 return AutoDiffRec<D,SCAL> (a+b.Rec(), b.Last());
761 }
762
763 template <int D, typename SCAL, typename SCAL2,
764 typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
765 INLINE AutoDiffRec<D,SCAL> operator+ (AutoDiffRec<D,SCAL> a, SCAL2 b)
766 {
767 return AutoDiffRec<D,SCAL> (a.Rec()+b, a.Last());
768 }
769
770 template <int D, typename SCAL>
771 INLINE AutoDiffRec<D,SCAL> operator+ (AutoDiffRec<D,SCAL> a, AutoDiffRec<D,SCAL> b)
772 {
773 return AutoDiffRec<D,SCAL> (a.Rec()+b.Rec(), a.Last()+b.Last());
774 }
775
776 template <int D, typename SCAL, typename SCAL2,
777 typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
778 INLINE AutoDiffRec<D,SCAL> operator- (SCAL2 b, AutoDiffRec<D,SCAL> a)
779 {
780 return AutoDiffRec<D,SCAL> (b-a.Rec(), -a.Last());
781 }
782
783 template <int D, typename SCAL, typename SCAL2,
784 typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
785 INLINE AutoDiffRec<D,SCAL> operator- (AutoDiffRec<D,SCAL> a, SCAL2 b)
786 {
787 return AutoDiffRec<D,SCAL> (a.Rec()-b, a.Last());
788 }
789
790 template <int D, typename SCAL>
791 INLINE AutoDiffRec<D,SCAL> operator- (AutoDiffRec<D,SCAL> a, AutoDiffRec<D,SCAL> b)
792 {
793 return AutoDiffRec<D,SCAL> (a.Rec()-b.Rec(), a.Last()-b.Last());
794 }
795
797 template<int D, typename SCAL>
802
803 template <int D, typename SCAL>
804 INLINE AutoDiffRec<D,SCAL> operator* (AutoDiffRec<D,SCAL> a, AutoDiffRec<D,SCAL> b)
805 {
806 return AutoDiffRec<D,SCAL> (a.Rec()*b.Rec(), a.Value()*b.Last()+b.Value()*a.Last());
807 }
808
809 template <int D, typename SCAL, typename SCAL1,
810 typename std::enable_if<std::is_convertible<SCAL1,SCAL>::value, int>::type = 0>
811 INLINE AutoDiffRec<D,SCAL> operator* (AutoDiffRec<D,SCAL> b, SCAL1 a)
812 {
813 return AutoDiffRec<D,SCAL> (a*b.Rec(), a*b.Last());
814 }
815
816 template <int D, typename SCAL, typename SCAL1,
817 typename std::enable_if<std::is_convertible<SCAL1,SCAL>::value, int>::type = 0>
818 INLINE AutoDiffRec<D,SCAL> operator* (SCAL1 a, AutoDiffRec<D,SCAL> b)
819 {
820 return AutoDiffRec<D,SCAL> (a*b.Rec(), a*b.Last());
821 }
822
823 template <int D, typename SCAL>
824 INLINE AutoDiffRec<D,SCAL> & operator+= (AutoDiffRec<D,SCAL> & a, AutoDiffRec<D,SCAL> b)
825 {
826 a.Rec() += b.Rec();
827 a.Last() += b.Last();
828 return a;
829 }
830
831 template <int D, typename SCAL>
832 INLINE AutoDiffRec<D,SCAL> & operator-= (AutoDiffRec<D,SCAL> & a, double b)
833 {
834 a.Rec() -= b;
835 return a;
836 }
837
838 template <int D, typename SCAL>
839 INLINE AutoDiffRec<D,SCAL> & operator-= (AutoDiffRec<D,SCAL> & a, AutoDiffRec<D,SCAL> b)
840 {
841 a.Rec() -= b.Rec();
842 a.Last() -= b.Last();
843 return a;
844 }
845
846
847 template <int D, typename SCAL>
848 INLINE AutoDiffRec<D,SCAL> & operator*= (AutoDiffRec<D,SCAL> & a, AutoDiffRec<D,SCAL> b)
849 {
850 a = a*b;
851 return a;
852 }
853
854
855 template <int D, typename SCAL, typename SCAL2>
856 INLINE AutoDiffRec<D,SCAL> & operator*= (AutoDiffRec<D,SCAL> & b, SCAL2 a)
857 {
858 b.Rec() *= a;
859 b.Last() *= a;
860 return b;
861 }
862
864
865 template <typename SCAL>
866 auto Inv1 (SCAL x) { return 1.0/x; }
867
868 template<int D, typename SCAL>
869 INLINE AutoDiffRec<D,SCAL> Inv1 (AutoDiffRec<D,SCAL> x)
870 {
871 return AutoDiffRec<D,SCAL> (Inv1(x.Rec()), (-sqr(1.0/x.Value())) * x.Last());
872 }
873
875 template<int D, typename SCAL>
877 {
878 return x * Inv1 (y);
879 }
880
881
883 template<int D, typename SCAL, typename SCAL2,
884 typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
886 {
887 return (1.0/y) * x;
888 }
889
890
892 template<int D, typename SCAL, typename SCAL2,
893 typename std::enable_if<std::is_convertible<SCAL2,SCAL>::value, int>::type = 0>
895 {
896 return x * Inv1(y);
897 }
898
899
900
901
902
903
904
906 template <int D, typename SCAL>
907 INLINE bool operator== (AutoDiffRec<D,SCAL> x, SCAL val2)
908 {
909 return x.Value() == val2;
910 }
911
913 template <int D, typename SCAL>
914 INLINE bool operator!= (AutoDiffRec<D,SCAL> x, SCAL val2) throw()
915 {
916 return x.Value() != val2;
917 }
918
920 template <int D, typename SCAL>
921 INLINE bool operator< (AutoDiffRec<D,SCAL> x, SCAL val2) throw()
922 {
923 return x.Value() < val2;
924 }
925
927 template <int D, typename SCAL>
928 INLINE bool operator> (AutoDiffRec<D,SCAL> x, SCAL val2) throw()
929 {
930 return x.Value() > val2;
931 }
932
933 using std::fabs;
934 template<int D, typename SCAL>
935 INLINE AutoDiffRec<D,SCAL> fabs (const AutoDiffRec<D,SCAL> & x)
936 {
937 auto sign = IfPos(x.Value(), SCAL(1.0), IfPos(-x.Value(), SCAL(-1.0), SCAL(0.0)));
938 return AutoDiffRec<D,SCAL> (fabs(x.Rec()), sign*x.Last());
939 // return fabs (AutoDiffVec<D,SCAL>(x));
940 /*
941 double abs = fabs (x.Value());
942 AutoDiffVec<D,SCAL> res( abs );
943 if (abs != 0.0)
944 for (int i = 0; i < D; i++)
945 res.DValue(i) = x.Value()*x.DValue(i) / abs;
946 else
947 for (int i = 0; i < D; i++)
948 res.DValue(i) = 0.0;
949 return res;
950 */
951 }
952
953
954 template<int D, typename SCAL>
955 INLINE auto sqrt (const AutoDiffRec<D,SCAL> & x)
956 {
957 return AutoDiffRec<D,SCAL> (sqrt(x.Rec()), (0.5/sqrt(x.Value()))*x.Last());
958 }
959
960
961
962 template <int D, typename SCAL>
963 auto log (AutoDiffRec<D,SCAL> x)
964 {
965 return AutoDiffRec<D,SCAL> (log(x.Rec()), (1.0/x.Value())*x.Last());
966 }
967
968 template <int D, typename SCAL>
969 auto exp (AutoDiffRec<D,SCAL> x)
970 {
971 return AutoDiffRec<D,SCAL> (exp(x.Rec()), exp(x.Value())*x.Last());
972 }
973
974 template <int D, typename SCAL>
975 INLINE AutoDiffRec<D,SCAL> pow (AutoDiffRec<D,SCAL> x, AutoDiffRec<D,SCAL> y )
976 {
977 return exp(log(x)*y);
978 }
979
980
981 template <int D, typename SCAL>
982 auto sin (AutoDiffRec<D,SCAL> x)
983 {
984 return AutoDiffRec<D,SCAL> (sin(x.Rec()), cos(x.Value())*x.Last());
985 }
986
987 template <int D, typename SCAL>
988 auto cos (AutoDiffRec<D,SCAL> x)
989 {
990 return AutoDiffRec<D,SCAL> (cos(x.Rec()), -sin(x.Value())*x.Last());
991 }
992
993 template <int D, typename SCAL>
994 auto tan (AutoDiffRec<D,SCAL> x)
995 {
996 return sin(x) / cos(x);
997 }
998
999 template <int D, typename SCAL>
1000 auto sinh (AutoDiffRec<D,SCAL> x)
1001 {
1002 return AutoDiffRec<D,SCAL> (sinh(x.Rec()), cosh(x.Value())*x.Last());
1003 }
1004
1005 template <int D, typename SCAL>
1006 auto cosh (AutoDiffRec<D,SCAL> x)
1007 {
1008 return AutoDiffRec<D,SCAL> (cosh(x.Rec()), sinh(x.Value())*x.Last());
1009 }
1010
1011 template <int D, typename SCAL>
1012 auto erf (AutoDiffRec<D,SCAL> x)
1013 {
1014 return AutoDiffRec<D,SCAL> (erf(x.Rec()), 2. / sqrt(M_PI) * exp(- x.Value() * x.Value())*x.Last());
1015 }
1016
1017 template <int D, typename SCAL>
1018 auto floor (AutoDiffRec<D,SCAL> x)
1019 {
1020 return AutoDiffRec<D,SCAL> (floor(x.Rec()), 0.0);
1021 }
1022
1023 template <int D, typename SCAL>
1024 auto ceil (AutoDiffRec<D,SCAL> x)
1025 {
1026 return AutoDiffRec<D,SCAL> (ceil(x.Rec()), 0.0);
1027 }
1028
1029
1030
1031 template <int D, typename SCAL>
1032 auto atan (AutoDiffRec<D,SCAL> x)
1033 {
1034 return AutoDiffRec<D,SCAL> (atan(x.Rec()), (1./(1.+x.Value()*x.Value()))*x.Last());
1035 }
1036
1037 template <int D, typename SCAL>
1038 auto atan2 (AutoDiffRec<D,SCAL> x, AutoDiffRec<D,SCAL> y)
1039 {
1040 return AutoDiffRec<D,SCAL> (atan2(x.Rec(), y.Rec()),
1041 (1./(x.Value()*x.Value()+y.Value()*y.Value()))*(x.Value()*y.Last()-y.Value()*x.Last()));
1042 }
1043
1044 template <int D, typename SCAL>
1045 auto acos (AutoDiffRec<D,SCAL> x)
1046 {
1047 return AutoDiffRec<D,SCAL> (acos(x.Rec()), (-1./sqrt(1.-x.Value()*x.Value()))*x.Last());
1048 }
1049
1050 template <int D, typename SCAL>
1051 auto asin (AutoDiffRec<D,SCAL> x)
1052 {
1053 return AutoDiffRec<D,SCAL> (asin(x.Rec()), (1./sqrt(1.-x.Value()*x.Value()))*x.Last());
1054 }
1055
1056
1057 template <int D, typename SCAL, typename TB, typename TC>
1058 auto IfPos (AutoDiffRec<D,SCAL> a, TB b, TC c) // -> decltype(IfPos (a.Value(), b, c))
1059 {
1060 return IfPos (a.Value(), b, c);
1061 }
1062
1063 template <int D, typename SCAL>
1064 INLINE AutoDiffRec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffRec<D,SCAL> b, AutoDiffRec<D,SCAL> c)
1065 {
1066 /*
1067 AutoDiffRec<D,SCAL> res;
1068 res.Value() = IfPos (a, b.Value(), c.Value());
1069 for (int j = 0; j < D; j++)
1070 res.DValue(j) = IfPos (a, b.DValue(j), c.DValue(j));
1071 return res;
1072 */
1073 return AutoDiffRec<D,SCAL> (IfPos(a, b.Rec(), c.Rec()), IfPos(a, b.Last(), c.Last()));
1074 }
1075
1076 template <int D, typename SCAL, typename TC>
1077 INLINE AutoDiffRec<D,SCAL> IfPos (SCAL /* SIMD<double> */ a, AutoDiffRec<D,SCAL> b, TC c)
1078 {
1079 return IfPos (a, b, AutoDiffRec<D,SCAL> (c));
1080 }
1081
1082
1083
1084template <int D, typename SCAL = double>
1085using AutoDiff = AutoDiffRec<D,SCAL>;
1086
1087}
1088
1089
1090
1091namespace ngbla
1092{
1093 template <typename T> struct is_scalar_type;
1094 template <int D, typename T>
1095 struct is_scalar_type<ngstd::AutoDiff<D,T>> { static constexpr bool value = true; };
1096
1097 // not meaningful for AutoDiff<D,Complex>, since this is
1098 // not (complex) differentiable anyway
1099 template<int D, typename SCAL>
1100 inline auto L2Norm2 (const ngstd::AutoDiff<D,SCAL> & x)
1101 {
1102 return x*x;
1103 }
1104
1105 template<int D, typename SCAL>
1106 inline auto L2Norm (const ngstd::AutoDiff<D,SCAL> & x) throw()
1107 {
1108 return IfPos(x.Value(), x, -x);
1109 }
1110
1111
1112
1113 template<int D, typename TAD>
1114 INLINE auto Conj (const ngstd::AutoDiff<D,TAD> & a)
1115 {
1117 b.Value() = conj(a.Value());
1118
1119 for(int i=0;i<D;i++)
1120 b.DValue(i) = conj(a.DValue(i));
1121
1122 return b;
1123 }
1124
1125
1126}
1127
1128
1129#endif
1130
1131
Definition autodiff.hpp:643
Datatype for automatic differentiation.
Definition autodiff.hpp:26
INLINE AutoDiffVec(SCAL aval)
initial object with constant value
Definition autodiff.hpp:51
INLINE SCAL & Value()
access value
Definition autodiff.hpp:104
AutoDiffVec(const AutoDiffVec &ad2)=default
copy constructor
AutoDiffVec()=default
elements are undefined
INLINE SCAL & DValue(int i)
accesses partial derivative
Definition autodiff.hpp:107
INLINE AutoDiffVec(SCAL aval, int diffindex)
init object with (val, e_diffindex)
Definition autodiff.hpp:59
INLINE AutoDiffVec & operator=(SCAL aval)
assign constant value
Definition autodiff.hpp:74
INLINE SCAL DValue(int i) const
returns partial derivative
Definition autodiff.hpp:88
INLINE SCAL Value() const
returns value
Definition autodiff.hpp:85
namespace for standard data types and algorithms.
Definition ngstd.hpp:42
ostream & operator<<(ostream &ost, const AutoDiffDiff< D, SCAL > &x)
Prints AudoDiffDiff.
Definition autodiffdiff.hpp:256
auto Inv1(SCAL x)
Inverse of AutoDiffRec.
Definition autodiff.hpp:866
Definition autodiffdiff.hpp:716