MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
dual.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
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 /**
13  * @file dual.hpp
14  *
15  * @brief This file contains the declaration of a dual number class
16  */
17 
18 #ifndef MFEM_INTERNAL_DUAL_HPP
19 #define MFEM_INTERNAL_DUAL_HPP
20 
21 #include <type_traits> // for is_arithmetic
22 #include <cmath>
23 #include "../general/backends.hpp"
24 
25 namespace mfem
26 {
27 namespace internal
28 {
29 
30 /**
31  * @brief Dual number struct (value plus gradient)
32  * @tparam gradient_type The type of the gradient (should support addition,
33  * scalar multiplication/division, and unary negation operators)
34  */
35 template <typename value_type, typename gradient_type>
36 struct dual
37 {
38  /// the actual numerical value
39  value_type value;
40  /// the partial derivatives of value w.r.t. some other quantity
41  gradient_type gradient;
42 
43  /** @brief assignment of a double to a value of a dual. Promotes a double to
44  * a dual with a zero gradient value. */
45  auto operator=(double a) -> dual<value_type, gradient_type>&
46  {
47  value = a;
48  gradient = {};
49  return *this;
50  }
51 };
52 
53 /** @brief class for checking if a type is a dual number or not */
54 template <typename T>
55 struct is_dual_number
56 {
57  /// whether or not type T is a dual number
58  static constexpr bool value = false;
59 };
60 
61 /** @brief class for checking if a type is a dual number or not */
62 template <typename value_type, typename gradient_type>
63 struct is_dual_number<dual<value_type, gradient_type> >
64 {
65  static constexpr bool value = true; ///< whether or not type T is a dual number
66 };
67 
68 /** @brief addition of a dual number and a non-dual number */
69 template <typename other_type, typename value_type, typename gradient_type,
70  typename = typename std::enable_if<
71  std::is_arithmetic<other_type>::value ||
72  is_dual_number<other_type>::value>::type>
73 MFEM_HOST_DEVICE
74 constexpr auto operator+(dual<value_type, gradient_type> a,
75  other_type b) -> dual<value_type, gradient_type>
76 {
77  return {a.value + b, a.gradient};
78 }
79 
80 // C++17 version of the above
81 //
82 // template <typename value_type, typename gradient_type>
83 // constexpr auto operator+(dual<value_type, gradient_type> a, value_type b)
84 // {
85 // return dual{a.value + b, a.gradient};
86 // }
87 
88 /** @brief addition of a dual number and a non-dual number */
89 template <typename other_type, typename value_type, typename gradient_type,
90  typename = typename std::enable_if<
91  std::is_arithmetic<other_type>::value ||
92  is_dual_number<other_type>::value>::type>
93 MFEM_HOST_DEVICE
94 constexpr auto operator+(other_type a,
95  dual<value_type, gradient_type> b) -> dual<value_type, gradient_type>
96 {
97  return {a + b.value, b.gradient};
98 }
99 
100 /** @brief addition of two dual numbers */
101 template <typename value_type_a, typename gradient_type_a, typename value_type_b, typename gradient_type_b>
102 MFEM_HOST_DEVICE
103 constexpr auto operator+(dual<value_type_a, gradient_type_a> a,
104  dual<value_type_b, gradient_type_b> b) -> dual<decltype(a.value + b.value),
105  decltype(a.gradient + b.gradient)>
106 {
107  return {a.value + b.value, a.gradient + b.gradient};
108 }
109 
110 /** @brief unary negation of a dual number */
111 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
112 constexpr auto operator-(dual<value_type, gradient_type> x) ->
113 dual<value_type, gradient_type>
114 {
115  return {-x.value, -x.gradient};
116 }
117 
118 /** @brief subtraction of a non-dual number from a dual number */
119 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
120 constexpr auto operator-(dual<value_type, gradient_type> a,
121  double b) -> dual<value_type, gradient_type>
122 {
123  return {a.value - b, a.gradient};
124 }
125 
126 /** @brief subtraction of a dual number from a non-dual number */
127 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
128 constexpr auto operator-(double a,
129  dual<value_type, gradient_type> b) -> dual<value_type, gradient_type>
130 {
131  return {a - b.value, -b.gradient};
132 }
133 
134 /** @brief subtraction of two dual numbers */
135 template <typename value_type_a, typename gradient_type_a, typename value_type_b, typename gradient_type_b>
136 MFEM_HOST_DEVICE
137 constexpr auto operator-(dual<value_type_a, gradient_type_a> a,
138  dual<value_type_b, gradient_type_b> b) -> dual<decltype(a.value - b.value),
139  decltype(a.gradient - b.gradient)>
140 {
141  return {a.value - b.value, a.gradient - b.gradient};
142 }
143 
144 /** @brief multiplication of a dual number and a non-dual number */
145 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
146 constexpr auto operator*(const dual<value_type, gradient_type>& a,
147  double b) -> dual<decltype(a.value * b), decltype(a.gradient * b)>
148 {
149  return {a.value * b, a.gradient * b};
150 }
151 
152 /** @brief multiplication of a dual number and a non-dual number */
153 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
154 constexpr auto operator*(double a,
155  const dual<value_type, gradient_type>& b) ->
156 dual<decltype(a * b.value), decltype(a * b.gradient)>
157 {
158  return {a * b.value, a * b.gradient};
159 }
160 
161 /** @brief multiplication of two dual numbers */
162 template <typename value_type_a, typename gradient_type_a, typename value_type_b, typename gradient_type_b>
163 MFEM_HOST_DEVICE
164 constexpr auto operator*(dual<value_type_a, gradient_type_a> a,
165  dual<value_type_b, gradient_type_b> b) -> dual<decltype(a.value * b.value),
166  decltype(b.value * a.gradient + a.value * b.gradient)>
167 {
168  return {a.value * b.value, b.value * a.gradient + a.value * b.gradient};
169 }
170 
171 /** @brief division of a dual number by a non-dual number */
172 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
173 constexpr auto operator/(const dual<value_type, gradient_type>& a,
174  double b) -> dual<decltype(a.value / b), decltype(a.gradient / b)>
175 {
176  return {a.value / b, a.gradient / b};
177 }
178 
179 /** @brief division of a non-dual number by a dual number */
180 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
181 constexpr auto operator/(double a,
182  const dual<value_type, gradient_type>& b) -> dual<decltype(a / b.value),
183  decltype(-(a / (b.value * b.value)) * b.gradient)>
184 {
185  return {a / b.value, -(a / (b.value * b.value)) * b.gradient};
186 }
187 
188 /** @brief division of two dual numbers */
189 template <typename value_type_a, typename gradient_type_a, typename value_type_b, typename gradient_type_b>
190 MFEM_HOST_DEVICE
191 constexpr auto operator/(dual<value_type_a, gradient_type_a> a,
192  dual<value_type_b, gradient_type_b> b) -> dual<decltype(a.value / b.value),
193  decltype((a.gradient / b.value) -
194  (a.value * b.gradient) /
195  (b.value * b.value))>
196 {
197  return {a.value / b.value, (a.gradient / b.value) - (a.value * b.gradient) / (b.value * b.value)};
198 }
199 
200 /**
201  * @brief Generates const + non-const overloads for a binary comparison operator
202  * Comparisons are conducted against the "value" part of the dual number
203  * @param[in] x The comparison operator to overload
204  */
205 #define mfem_binary_comparator_overload(x) \
206  template <typename value_type, typename gradient_type> \
207  MFEM_HOST_DEVICE constexpr bool operator x( \
208  const dual<value_type, gradient_type>& a, \
209  double b) \
210  { \
211  return a.value x b; \
212  } \
213  \
214  template <typename value_type, typename gradient_type> \
215  MFEM_HOST_DEVICE constexpr bool operator x( \
216  double a, \
217  const dual<value_type, gradient_type>& b) \
218  { \
219  return a x b.value; \
220  } \
221  \
222  template <typename value_type_a, \
223  typename gradient_type_a, \
224  typename value_type_b, \
225  typename gradient_type_b> MFEM_HOST_DEVICE \
226  constexpr bool operator x( \
227  const dual<value_type_a, gradient_type_a>& a, \
228  const dual<value_type_b, gradient_type_b>& b) \
229  { \
230  return a.value x b.value; \
231  }
232 
233 mfem_binary_comparator_overload(<) ///< implement operator< for dual numbers
234 mfem_binary_comparator_overload(<=) ///< implement operator<= for dual numbers
235 mfem_binary_comparator_overload(==) ///< implement operator== for dual numbers
236 mfem_binary_comparator_overload(>=) ///< implement operator>= for dual numbers
237 mfem_binary_comparator_overload(>) ///< implement operator> for dual numbers
238 
239 #undef mfem_binary_comparator_overload
240 
241 /** @brief compound assignment (+) for dual numbers */
242 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
243 dual<value_type, gradient_type>& operator+=(dual<value_type, gradient_type>& a,
244  const dual<value_type, gradient_type>& b)
245 {
246  a.value += b.value;
247  a.gradient += b.gradient;
248  return a;
249 }
250 
251 /** @brief compound assignment (-) for dual numbers */
252 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
253 dual<value_type, gradient_type>& operator-=(dual<value_type, gradient_type>& a,
254  const dual<value_type, gradient_type>& b)
255 {
256  a.value -= b.value;
257  a.gradient -= b.gradient;
258  return a;
259 }
260 
261 /** @brief compound assignment (+) for dual numbers with `double` righthand side */
262 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
263 dual<value_type, gradient_type>& operator+=(dual<value_type, gradient_type>& a,
264  double b)
265 {
266  a.value += b;
267  return a;
268 }
269 
270 /** @brief compound assignment (-) for dual numbers with `double` righthand side */
271 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
272 dual<value_type, gradient_type>& operator-=(dual<value_type, gradient_type>& a,
273  double b)
274 {
275  a.value -= b;
276  return a;
277 }
278 
279 /** @brief implementation of absolute value function for dual numbers */
280 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
281 dual<value_type, gradient_type> abs(dual<value_type, gradient_type> x)
282 {
283  return (x.value >= 0) ? x : -x;
284 }
285 
286 /** @brief implementation of square root for dual numbers */
287 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
288 dual<value_type, gradient_type> sqrt(dual<value_type, gradient_type> x)
289 {
290  return {std::sqrt(x.value), x.gradient / (2.0 * std::sqrt(x.value))};
291 }
292 
293 /** @brief implementation of cosine for dual numbers */
294 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
295 dual<value_type, gradient_type> cos(dual<value_type, gradient_type> a)
296 {
297  return {std::cos(a.value), -a.gradient * std::sin(a.value)};
298 }
299 
300 /** @brief implementation of sine for dual numbers */
301 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
302 dual<value_type, gradient_type> sin(dual<value_type, gradient_type> a)
303 {
304  return {std::sin(a.value), a.gradient * std::cos(a.value)};
305 }
306 
307 /** @brief implementation of sinh for dual numbers */
308 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
309 dual<value_type, gradient_type> sinh(dual<value_type, gradient_type> a)
310 {
311  return {std::sinh(a.value), a.gradient * std::cosh(a.value)};
312 }
313 
314 /** @brief implementation of acos for dual numbers */
315 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
316 dual<value_type, gradient_type> acos(dual<value_type, gradient_type> a)
317 {
318  using std::sqrt;
319  using std::acos;
320  return {acos(a.value), -a.gradient / sqrt(value_type{1} - a.value * a.value)};
321 }
322 
323 /** @brief implementation of asin for dual numbers */
324 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
325 dual<value_type, gradient_type> asin(dual<value_type, gradient_type> a)
326 {
327  using std::sqrt;
328  using std::asin;
329  return {asin(a.value), a.gradient / sqrt(value_type{1} - a.value * a.value)};
330 }
331 
332 /** @brief implementation of tan for dual numbers */
333 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
334 dual<value_type, gradient_type> tan(dual<value_type, gradient_type> a)
335 {
336  using std::tan;
337  value_type f = tan(a.value);
338  return {f, a.gradient * (value_type{1} + f * f)};
339 }
340 
341 /** @brief implementation of atan for dual numbers */
342 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
343 dual<value_type, gradient_type> atan(dual<value_type, gradient_type> a)
344 {
345  return {atan(a.value), a.gradient / (value_type{1} + a.value * a.value)};
346 }
347 
348 /** @brief implementation of exponential function for dual numbers */
349 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
350 dual<value_type, gradient_type> exp(dual<value_type, gradient_type> a)
351 {
352  return {std::exp(a.value), std::exp(a.value) * a.gradient};
353 }
354 
355 /** @brief implementation of the natural logarithm function for dual numbers */
356 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
357 dual<value_type, gradient_type> log(dual<value_type, gradient_type> a)
358 {
359  return {std::log(a.value), a.gradient / a.value};
360 }
361 
362 /** @brief implementation of `a` (dual) raised to the `b` (dual) power */
363 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
364 dual<value_type, gradient_type> pow(dual<value_type, gradient_type> a,
365  dual<value_type, gradient_type> b)
366 {
367  value_type value = pow(a.value, b.value);
368  return {value, value * (a.gradient * (b.value / a.value) + b.gradient * std::log(a.value))};
369 }
370 
371 /** @brief implementation of `a` (non-dual) raised to the `b` (dual) power */
372 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
373 dual<value_type, gradient_type> pow(double a, dual<value_type, gradient_type> b)
374 {
375  value_type value = pow(a, b.value);
376  return {value, value * b.gradient * std::log(a)};
377 }
378 
379 /** @brief implementation of `a` (non-dual) raised to the `b` (non-dual) power */
380 template <typename value_type > MFEM_HOST_DEVICE
381 value_type pow(value_type a, value_type b) { return std::pow(a, b); }
382 
383 /** @brief implementation of `a` (dual) raised to the `b` (non-dual) power */
384 template <typename value_type, typename gradient_type> MFEM_HOST_DEVICE
385 dual<value_type, gradient_type> pow(dual<value_type, gradient_type> a, double b)
386 {
387  value_type value = pow(a.value, b);
388  return {value, value * a.gradient * b / a.value};
389 }
390 
391 /** @brief overload of operator<< for `dual` to work with work with standard output streams */
392 template <typename value_type, typename gradient_type, int... n>
393 std::ostream& operator<<(std::ostream& os, dual<value_type, gradient_type> A)
394 {
395  os << '(' << A.value << ' ' << A.gradient << ')';
396  return os;
397 }
398 
399 /** @brief promote a value to a dual number of the appropriate type */
400 MFEM_HOST_DEVICE constexpr dual<double, double> make_dual(double x) { return {x, 1.0}; }
401 
402 /** @brief return the "value" part from a given type. For non-dual types, this is just the identity function */
403 template <typename T> MFEM_HOST_DEVICE T get_value(const T& arg) { return arg; }
404 
405 /** @brief return the "value" part from a dual number type */
406 template <typename value_type, typename gradient_type>
407 MFEM_HOST_DEVICE gradient_type get_value(dual<value_type, gradient_type> arg)
408 {
409  return arg.value;
410 }
411 
412 /** @brief return the "gradient" part from a dual number type */
413 template <typename value_type, typename gradient_type>
414 MFEM_HOST_DEVICE gradient_type get_gradient(dual<value_type, gradient_type> arg)
415 {
416  return arg.gradient;
417 }
418 
419 } // namespace internal
420 } // namespace mfem
421 
422 #endif
MFEM_ALWAYS_INLINE AutoSIMD< scalar_t, S, A > operator/(const scalar_t &e, const AutoSIMD< scalar_t, S, A > &v)
Definition: auto.hpp:271
MFEM_ALWAYS_INLINE AutoSIMD< scalar_t, S, A > operator+(const scalar_t &e, const AutoSIMD< scalar_t, S, A > &v)
Definition: auto.hpp:238
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
double b
Definition: lissajous.cpp:42
double a
Definition: lissajous.cpp:41
MemoryClass operator*(MemoryClass mc1, MemoryClass mc2)
Return a suitable MemoryClass from a pair of MemoryClasses.
MFEM_ALWAYS_INLINE AutoSIMD< scalar_t, S, A > operator-(const scalar_t &e, const AutoSIMD< scalar_t, S, A > &v)
Definition: auto.hpp:249