MFEM  v4.4.0 Finite element discretization library
fdual.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #ifndef FDUAL_H
13 #define FDUAL_H
14
15 #include <cmath>
16 #include <type_traits>
17
18 namespace mfem
19 {
21 {
22 /** The FDualNumber template class provides forward automatic differentiation
23  (see https://en.wikipedia.org/wiki/Automatic_differentiation) implementation
24  based on dual numbers.
25
26
27  The derivative of an arbitrary function double f(double a) can be obtained
28  by replacing the double type for the return value and the argument a with
29  FDualNumber<double>, i.e., FDualNumber<double> f(FDualNumber<double> a). The
30  derivative is evaluated automatically by calling the function r=f(a). The
31  value of the function is stored in r.pr and the derivative in r.du. These
32  can be extracted by the corresponding methods real()/prim() and dual().
33
34  Internally, the function f can be composed of standard functions predefined
35  for FDualNumber type. These consist of a large set of functions replicating
36  the functionality of the standard math library, i.e., sin, cos, exp, log,
37  etc. New functions (non-member) can be equally added to the class. Example:
38
39  \code{.cpp}
40  template<typename tbase>
41  inline FDualNumber<tbase> cos(const FDualNumber<tbase> &f)
42  {
43  return FDualNumber<tbase>(cos(f.real()), -f.dual() * sin(f.real()));
44  }
45  \endcode
46
47  The real part of the return value consists of the standard real value of the
48  function, i.e., cos(f.real()).
49
50  The dual part of the return value consists of the first derivative of the
51  function with respect to the real part of the argument -sin(f.reaf)
52  multiplied with the dual part of the argument f.dual().
53 */
54 template<typename tbase>
56 {
57 private:
58  /// Real value
59  tbase pr;
60  /// Dual value holding derivative information
61  tbase du;
62
63 public:
64  /// Standard constructor - both values are set to zero.
65  FDualNumber() : pr(0), du(0) {}
66
67  /// The constructor utilized in nested definition of dual numbers. It is
68  /// used for second and higher order derivatives.
69  template<class fltyp,
70  class = typename std::enable_if<std::is_arithmetic<fltyp>::value>::type>
71  FDualNumber(fltyp &f) : pr(f), du(0)
72  {}
73
74  /// The constructor utilized in nested definition of dual numbers. It is
75  /// used for second and higher order derivatives.
76  template<class fltyp,
77  class = typename std::enable_if<std::is_arithmetic<fltyp>::value>::type>
78  FDualNumber(const fltyp &f) : pr(f), du(0)
79  {}
80
81  /// Standard constructor with user supplied input for both parts of the dual
82  /// number.
83  FDualNumber(tbase &pr_, tbase &du_) : pr(pr_), du(du_) {}
84
85  /// Standard constructor with user supplied input for both parts of the dual
86  /// number.
87  FDualNumber(const tbase &pr_, const tbase &du_) : pr(pr_), du(du_) {}
88
89  /// Standard constructor with user supplied dual number.
90  FDualNumber(FDualNumber<tbase> &nm) : pr(nm.pr), du(nm.du) {}
91
92  /// Standard constructor with user supplied dual number.
93  FDualNumber(const FDualNumber<tbase> &nm) : pr(nm.pr), du(nm.du) {}
94
95  /// Return the real value of the dual number.
96  tbase prim() const { return pr; }
97
98  /// Same as prim(). Return the real value of the dual number.
99  tbase real() const { return pr; }
100
101  /// Return the dual value of the dual number.
102  tbase dual() const { return du; }
103
104  /// Set the primal and the dual values.
105  void set(const tbase &pr_, const tbase &du_)
106  {
107  pr = pr_;
108  du = du_;
109  }
110
111  /// Set the primal value.
112  void prim(const tbase &pr_) { pr = pr_; }
113
114  /// Set the primal value.
115  void real(const tbase &pr_) { pr = pr_; }
116
117  /// Set the dual value.
118  void dual(const tbase &du_) { du = du_; }
119
120  /// Set the primal value.
121  void setReal(const tbase &pr_) { pr = pr_; }
122
123  /// Set the dual value.
124  void setDual(const tbase &du_) { du = du_; }
125
126  /// operator =
128  {
129  pr = sc_;
130  du = tbase(0);
131  return *this;
132  }
133
134  /// operator +=
136  {
137  pr = pr + sc_;
138  return *this;
139  }
140
141  /// operator -=
143  {
144  pr = pr - sc_;
145  return *this;
146  }
147
148  /// operator *=
150  {
151  pr = pr * sc_;
152  du = du * sc_;
153  return *this;
154  }
155
156  /// operator /=
158  {
159  pr = pr / sc_;
160  du = du / sc_;
161  return *this;
162  }
163
164  /// operator =
166  {
167  pr = f.real();
168  du = f.dual();
169  return *this;
170  }
171
172  /// operator +=
174  {
175  pr += f.real();
176  du += f.dual();
177  return *this;
178  }
179
180  /// operator -=
182  {
183  pr -= f.real();
184  du -= f.dual();
185  return *this;
186  }
187
188  /// operator *=
190  {
191  du = du * f.real();
192  du = du + pr * f.dual();
193  pr = pr * f.real();
194  return *this;
195  }
196
197  /// operator /=
199  {
200  pr = pr / f_.real();
201  du = du - pr * f_.dual();
202  du = du / f_.real();
203  return *this;
204  }
205 };
206
207 /// non-member functions
208 /// boolean operation ==
209 template<typename tbase>
210 inline bool operator==(const FDualNumber<tbase> &a1,
211  const FDualNumber<tbase> &a2)
212 {
213  return a1.real() == a2.real();
214 }
215
216 /// boolean operation ==
217 template<typename tbase>
218 inline bool operator==(tbase a, const FDualNumber<tbase> &f_)
219 {
220  return a == f_.real();
221 }
222
223 /// boolean operation ==
224 template<typename tbase>
225 inline bool operator==(const FDualNumber<tbase> &a, tbase b)
226 {
227  return a.real() == b;
228 }
229
230 /// boolean operation <
231 template<typename tbase>
232 inline bool operator<(const FDualNumber<tbase> &f1,
233  const FDualNumber<tbase> &f2)
234 {
235  return f1.real() < f2.real();
236 }
237
238 /// boolean operation <
239 template<typename tbase>
240 inline bool operator<(const FDualNumber<tbase> &f, tbase a)
241 {
242  return f.real() < a;
243 }
244
245 /// boolean operation <
246 template<typename tbase>
247 inline bool operator<(tbase a, const FDualNumber<tbase> &f)
248 {
249  return a < f.real();
250 }
251
252 /// boolean operation >
253 template<typename tbase>
254 inline bool operator>(const FDualNumber<tbase> &f1,
255  const FDualNumber<tbase> &f2)
256 {
257  return f1.real() > f2.real();
258 }
259
260 /// boolean operation >
261 template<typename tbase>
262 inline bool operator>(const FDualNumber<tbase> &f, tbase a)
263 {
264  return f.real() > a;
265 }
266
267 /// boolean operation >
268 template<typename tbase>
269 inline bool operator>(tbase a, const FDualNumber<tbase> &f)
270 {
271  return (a > f.real());
272 }
273
274 /// Negate the real and the dual parts.
275 template<typename tbase>
277 {
278  return FDualNumber<tbase>(-f.real(), -f.dual());
279 }
280
281 /// [dual number] - [base number]
282 template<typename tbase>
284 {
285  return FDualNumber<tbase>(f.real() - a, f.dual());
286 }
287
288 /// [dual number<dual number>] - [base number]
289 template<typename tbase>
292 {
293  return FDualNumber<FDualNumber<tbase>>(f.real() - a, f.dual());
294 }
295
296 /// [dual number] + [base number]
297 template<typename tbase>
299 {
300  return FDualNumber<tbase>(f.real() + a, f.dual());
301 }
302
303 /// [dual number<dual number>] + [base number]
304 template<typename tbase>
307 {
308  return FDualNumber<FDualNumber<tbase>>(f.real() + a, f.dual());
309 }
310
311 /// [dual number] * [base number]
312 template<typename tbase>
314 {
315  return FDualNumber<tbase>(f.real() * a, f.dual() * a);
316 }
317
318 /// [dual number] / [base number]
319 template<typename tbase>
321 {
322  return FDualNumber<tbase>(f.real() / a, f.dual() / a);
323 }
324
325 /// [dual number<dual number>] / [base number]
326 template<typename tbase>
329 {
330  return FDualNumber<FDualNumber<tbase>>(f.real() / a, f.dual() / a);
331 }
332
333 /// [base number] + [dual number]
334 template<typename tbase>
336 {
337  return FDualNumber<tbase>(a + f.real(), f.dual());
338 }
339
340 /// [base number] + [dual number<dual number>]
341 template<typename tbase>
344 {
345  return FDualNumber<FDualNumber<tbase>>(a + f.real(), f.dual());
346 }
347
348 /// [base number] - [dual number]
349 template<typename tbase>
351 {
352  return FDualNumber<tbase>(a - f.real(), -f.dual());
353 }
354
355 /// [base number] - [dual number<dual number>]
356 template<typename tbase>
359 {
360  return FDualNumber<FDualNumber<tbase>>(a - f.real(), -f.dual());
361 }
362
363 /// [base number] * [dual number]
364 template<typename tbase>
366 {
367  return FDualNumber<tbase>(f.real() * a, f.dual() * a);
368 }
369
370 /// [base number] * [dual number<dual number>]
371 template<typename tbase>
374 {
375  return FDualNumber<FDualNumber<tbase>>(f.real() * a, f.dual() * a);
376 }
377
378 /// [base number] / [dual number]
379 template<typename tbase>
381 {
382  a = a / f.real();
383  return FDualNumber<tbase>(a, -a * f.dual() / f.real());
384 }
385
386 /// [dual number] + [dual number]
387 template<typename tbase>
389  const FDualNumber<tbase> &f2)
390 {
391  return FDualNumber<tbase>(f1.real() + f2.real(), f1.dual() + f2.dual());
392 }
393
394 /// [dual number] - [dual number]
395 template<typename tbase>
397  const FDualNumber<tbase> &f2)
398 {
399  return FDualNumber<tbase>(f1.real() - f2.real(), f1.dual() - f2.dual());
400 }
401
402 /// [dual number] * [dual number]
403 template<typename tbase>
405  const FDualNumber<tbase> &f2)
406 {
407  return FDualNumber<tbase>(f1.real() * f2.real(),
408  f1.real() * f2.dual() + f1.dual() * f2.real());
409 }
410
411 /// [dual number] / [dual number]
412 template<typename tbase>
414  const FDualNumber<tbase> &f2)
415 {
416  tbase a = tbase(1) / f2.real();
417  tbase b = f1.real() * a;
418  return FDualNumber<tbase>(b, (f1.dual() - f2.dual() * b) * a);
419 }
420
421 /// acos([dual number])
422 template<typename tbase>
424 {
425  return FDualNumber<tbase>(acos(f.real()),
426  -f.dual() / sqrt(tbase(1) - f.real() * f.real()));
427 }
428
429 /// acos([dual number<double>])
430 template<>
432 {
433  return FDualNumber<double>(std::acos(f.real()),
434  -f.dual() / std::sqrt(double(1) - f.real() * f.real()));
435 }
436
437 /// asin([dual number])
438 template<typename tbase>
440 {
441  return FDualNumber<tbase>(asin(f.real()),
442  f.dual() / sqrt(tbase(1) - f.real() * f.real()));
443 }
444
445 /// asin([dual number<double>])
446 template<>
448 {
449  return FDualNumber<double>(std::asin(f.real()),
450  f.dual() / std::sqrt(double(1) - f.real() * f.real()));
451 }
452
453 /// atan([dual number])
454 template<typename tbase>
456 {
457  return FDualNumber<tbase>(atan(f.real()),
458  f.dual() / (tbase(1) + f.real() * f.real()));
459 }
460
461 /// atan([dual number<double>])
462 template<>
464 {
465  return FDualNumber<double>(std::atan(f.real()),
466  f.dual() / (double(1) + f.real() * f.real()));
467 }
468
469 /// cos([dual number])
470 template<typename tbase>
472 {
473  return FDualNumber<tbase>(cos(f.real()), -f.dual() * sin(f.real()));
474 }
475
476 /// cos([dual number<double>])
477 template<>
479 {
480  return FDualNumber<double>(std::cos(f.real()), -f.dual() * std::sin(f.real()));
481 }
482
483 /// cosh([dual number])
484 template<typename tbase>
486 {
487  return FDualNumber<tbase>(cosh(f.real()), f.dual() * sinh(f.real()));
488 }
489
490 /// cosh([dual number<double>])
491 template<>
493 {
494  return FDualNumber<double>(std::cosh(f.real()), f.dual() * std::sinh(f.real()));
495 }
496
497 /// exp([dual number])
498 template<typename tbase>
500 {
501  tbase x = exp(f.real());
502  return FDualNumber<tbase>(x, f.dual() * x);
503 }
504
505 /// exp([dual number<double>])
506 template<>
508 {
509  double x = std::exp(f.real());
510  return FDualNumber<double>(x, f.dual() * x);
511 }
512
513 /// log([dual number])
514 template<typename tbase>
516 {
517  return FDualNumber<tbase>(log(f.real()), f.dual() / f.real());
518 }
519
520 /// log([dual number<double>])
521 template<>
523 {
524  return FDualNumber<double>(std::log(f.real()), f.dual() / f.real());
525 }
526
527 /// log10([dual number])
528 template<typename tbase>
530 {
531  return log(f) / log(tbase(10));
532 }
533
534 /// log10([dual number<double>])
535 template<>
537 {
538  return log(f) / std::log(double(10));
539 }
540
541 /// pow([dual number],[dual number])
542 template<typename tbase>
544  const FDualNumber<tbase> &b)
545 {
546  return exp(log(a) * b);
547 }
548
549 /// pow([dual number], [base number])
550 template<typename tbase, typename tbase1>
551 inline FDualNumber<tbase> pow(const FDualNumber<tbase> &a, const tbase1 &b)
552 {
553  return exp(log(a) * tbase(b));
554 }
555
556 /// pow([base number], [dual number])
557 template<typename tbase, typename tbase1>
558 inline FDualNumber<tbase> pow(const tbase1 &a, const FDualNumber<tbase> &b)
559 {
560  return exp(log(tbase(a)) * b);
561 }
562
563 /// pow([base number], [dual number<double>])
564 template<>
565 inline FDualNumber<double> pow(const double &a, const FDualNumber<double> &b)
566 {
567  return exp(std::log(a) * b);
568 }
569
570 /// sin([dual number])
571 template<typename tbase>
573 {
574  return FDualNumber<tbase>(sin(f.real()), f.dual() * cos(f.real()));
575 }
576
577 /// sin([dual number<double>])
578 template<>
580 {
581  return FDualNumber<double>(std::sin(f.real()), f.dual() * std::cos(f.real()));
582 }
583
584 /// sinh([dual number])
585 template<typename tbase>
587 {
588  return FDualNumber<tbase>(sinh(f.real()), f.dual() * cosh(f.real()));
589 }
590
591 /// sinh([dual number<double>])
592 template<>
594 {
595  return FDualNumber<double>(std::sinh(f.real()), f.dual() * std::cosh(f.real()));
596 }
597
598 /// sqrt([dual number])
599 template<typename tbase>
601 {
602  tbase a = sqrt(f.real());
603  return FDualNumber<tbase>(a, f.dual() / (tbase(2) * a));
604 }
605
606 /// sqrt([dual number<double>])
607 template<>
609 {
610  double a = std::sqrt(f.real());
611  return FDualNumber<double>(a, f.dual() / (double(2) * a));
612 }
613
614 /// tan([dual number])
615 template<typename tbase>
617 {
618  tbase a = tan(f.real());
619  return FDualNumber<tbase>(a, f.dual() * (tbase(1) + a * a));
620 }
621
622 /// tan([dual number<double>])
623 template<>
625 {
626  double a = std::tan(f.real());
627  return FDualNumber<double>(a, f.dual() * (double(1) + a * a));
628 }
629
630 /// tanh([dual number])
631 template<typename tbase>
633 {
634  tbase a = tanh(f.real());
635  return FDualNumber<tbase>(a, f.dual() * (tbase(1) - a * a));
636 }
637
638 /// tanh([dual number<double>])
639 template<>
641 {
642  double a = std::tanh(f.real());
643  return FDualNumber<double>(a, f.dual() * (double(1) - a * a));
644 }
645
647
648 } // namespace mfem
649
650 #endif
FDualNumber< tbase > asin(const FDualNumber< tbase > &f)
asin([dual number])
Definition: fdual.hpp:439
FDualNumber< tbase > tan(const FDualNumber< tbase > &f)
tan([dual number])
Definition: fdual.hpp:616
FDualNumber(tbase &pr_, tbase &du_)
Definition: fdual.hpp:83
void setDual(const tbase &du_)
Set the dual value.
Definition: fdual.hpp:124
FDualNumber< tbase > & operator/=(const FDualNumber< tbase > &f_)
operator /=
Definition: fdual.hpp:198
FDualNumber< tbase > & operator*=(tbase sc_)
operator *=
Definition: fdual.hpp:149
FDualNumber< tbase > operator-(const FDualNumber< tbase > &f)
Negate the real and the dual parts.
Definition: fdual.hpp:276
FDualNumber< tbase > tanh(const FDualNumber< tbase > &f)
tanh([dual number])
Definition: fdual.hpp:632
FDualNumber< tbase > & operator+=(const FDualNumber< tbase > &f)
operator +=
Definition: fdual.hpp:173
tbase real() const
Same as prim(). Return the real value of the dual number.
Definition: fdual.hpp:99
void prim(const tbase &pr_)
Set the primal value.
Definition: fdual.hpp:112
FDualNumber< tbase > operator/(const FDualNumber< tbase > &f, tbase a)
[dual number] / [base number]
Definition: fdual.hpp:320
FDualNumber(fltyp &f)
Definition: fdual.hpp:71
FDualNumber(const fltyp &f)
Definition: fdual.hpp:78
FDualNumber< tbase > & operator-=(const FDualNumber< tbase > &f)
operator -=
Definition: fdual.hpp:181
FDualNumber< tbase > & operator=(tbase sc_)
operator =
Definition: fdual.hpp:127
FDualNumber()
Standard constructor - both values are set to zero.
Definition: fdual.hpp:65
FDualNumber< tbase > acos(const FDualNumber< tbase > &f)
acos([dual number])
Definition: fdual.hpp:423
FDualNumber(FDualNumber< tbase > &nm)
Standard constructor with user supplied dual number.
Definition: fdual.hpp:90
void dual(const tbase &du_)
Set the dual value.
Definition: fdual.hpp:118
FDualNumber< tbase > exp(const FDualNumber< tbase > &f)
exp([dual number])
Definition: fdual.hpp:499
void real(const tbase &pr_)
Set the primal value.
Definition: fdual.hpp:115
tbase prim() const
Return the real value of the dual number.
Definition: fdual.hpp:96
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
FDualNumber< tbase > operator*(const FDualNumber< tbase > &f, tbase a)
[dual number] * [base number]
Definition: fdual.hpp:313
double b
Definition: lissajous.cpp:42
FDualNumber< tbase > & operator-=(tbase sc_)
operator -=
Definition: fdual.hpp:142
FDualNumber< tbase > & operator+=(tbase sc_)
operator +=
Definition: fdual.hpp:135
FDualNumber< tbase > & operator=(const FDualNumber< tbase > &f)
operator =
Definition: fdual.hpp:165
FDualNumber< tbase > cos(const FDualNumber< tbase > &f)
cos([dual number])
Definition: fdual.hpp:471
FDualNumber(const tbase &pr_, const tbase &du_)
Definition: fdual.hpp:87
FDualNumber< tbase > pow(const FDualNumber< tbase > &a, const FDualNumber< tbase > &b)
pow([dual number],[dual number])
Definition: fdual.hpp:543
FDualNumber< tbase > sin(const FDualNumber< tbase > &f)
sin([dual number])
Definition: fdual.hpp:572
void setReal(const tbase &pr_)
Set the primal value.
Definition: fdual.hpp:121
FDualNumber< tbase > & operator*=(const FDualNumber< tbase > &f)
operator *=
Definition: fdual.hpp:189
FDualNumber< tbase > log(const FDualNumber< tbase > &f)
log([dual number])
Definition: fdual.hpp:515
FDualNumber< tbase > operator+(const FDualNumber< tbase > &f, tbase a)
[dual number] + [base number]
Definition: fdual.hpp:298
double a
Definition: lissajous.cpp:41
FDualNumber< tbase > sqrt(const FDualNumber< tbase > &f)
sqrt([dual number])
Definition: fdual.hpp:600
tbase dual() const
Return the dual value of the dual number.
Definition: fdual.hpp:102
bool operator>(const FDualNumber< tbase > &f1, const FDualNumber< tbase > &f2)
boolean operation &gt;
Definition: fdual.hpp:254
FDualNumber< tbase > & operator/=(tbase sc_)
operator /=
Definition: fdual.hpp:157
FDualNumber< tbase > log10(const FDualNumber< tbase > &f)
log10([dual number])
Definition: fdual.hpp:529
bool operator==(const FDualNumber< tbase > &a1, const FDualNumber< tbase > &a2)
Definition: fdual.hpp:210
FDualNumber< tbase > cosh(const FDualNumber< tbase > &f)
cosh([dual number])
Definition: fdual.hpp:485
void set(const tbase &pr_, const tbase &du_)
Set the primal and the dual values.
Definition: fdual.hpp:105
FDualNumber< tbase > atan(const FDualNumber< tbase > &f)
atan([dual number])
Definition: fdual.hpp:455
FDualNumber< tbase > sinh(const FDualNumber< tbase > &f)
sinh([dual number])
Definition: fdual.hpp:586
FDualNumber(const FDualNumber< tbase > &nm)
Standard constructor with user supplied dual number.
Definition: fdual.hpp:93