Reference documentation for deal.II version GIT relicensing-136-gb80d0be4af 2024-03-18 08:20:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
numbers.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2006 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_numbers_h
16#define dealii_numbers_h
17
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/base/types.h>
22
23#ifdef DEAL_II_WITH_CUDA
24# include <cuComplex.h>
25#endif
26
27#include <Kokkos_Macros.hpp>
28
29#include <cmath>
30#include <complex>
31#include <cstddef>
32#include <type_traits>
33
34#define DEAL_II_HOST_DEVICE KOKKOS_FUNCTION
35#define DEAL_II_CUDA_HOST_DEV DEAL_II_HOST_DEVICE
36#define DEAL_II_HOST_DEVICE_ALWAYS_INLINE KOKKOS_FORCEINLINE_FUNCTION
37
38// clang++ assumes that all constexpr functions are __host__ __device__ when
39// Kokkos was configured with CUDA or HIP support. This is problematic
40// when calling non-constexpr functions in constexpr functions. Hence, we
41// need a way to annotate functions explicitly as host-only.
42#if (defined(__clang__) && defined(__CUDA__)) || defined(KOKKOS_ENABLE_HIP)
43# define DEAL_II_HOST __host__
44#else
45# define DEAL_II_HOST
46#endif
47
48// Forward-declare the automatic differentiation types so we can add prototypes
49// for our own wrappers.
50#ifdef DEAL_II_WITH_ADOLC
51class adouble;
52namespace adtl
53{
54 class adouble;
55}
56#endif
57
59
60namespace internal
61{
78 template <typename Number>
80 {
84 constexpr static unsigned int max_width = 1;
85 };
86
93 template <>
95 {
99 constexpr static unsigned int max_width =
100#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512
101 8;
102#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256
103 4;
104#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128
105 2;
106#else
107 1;
108#endif
109 };
110
117 template <>
119 {
123 constexpr static unsigned int max_width =
124#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ALTIVEC__)
125 4;
126#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && defined(__AVX512F__)
127 16;
128#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)
129 8;
130#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__)
131 4;
132#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ARM_NEON)
133 4;
134#else
135 1;
136#endif
137 };
138
139
140} // namespace internal
141
142// forward declarations to support abs or sqrt operations on VectorizedArray
143#ifndef DOXYGEN
144template <typename Number,
145 std::size_t width =
147class VectorizedArray;
148template <typename T>
149struct EnableIfScalar;
150#endif
151
152#ifdef DEAL_II_WITH_ADOLC
153# ifndef DOXYGEN
154// Prototype some inline functions present in adolc_math.h for use in
155// NumberTraits.
156//
157// ADOL-C uses fabs(), but for genericity we want to use abs(). Simultaneously,
158// though, we don't want to include ADOL-C headers in this header since
159// numbers.h is in everything. To get around this: use C++ rules which permit
160// the use of forward-declared classes in function prototypes to declare some
161// functions which are defined in adolc_math.h. This permits us to write "using
162// ::abs;" in NumberTraits which will allow us to select the correct
163// overload (the one in ::) when instantiating NumberTraits for ADOL-C
164// types.
165
166adouble
167abs(const adouble &x);
168
169adtl::adouble
170abs(const adtl::adouble &x);
171# endif
172#endif
173
175
176namespace std
177{
178 template <typename Number, std::size_t width>
179 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
180 sqrt(const ::VectorizedArray<Number, width> &);
181 template <typename Number, std::size_t width>
182 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
183 abs(const ::VectorizedArray<Number, width> &);
184 template <typename Number, std::size_t width>
185 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
186 max(const ::VectorizedArray<Number, width> &,
187 const ::VectorizedArray<Number, width> &);
188 template <typename Number, std::size_t width>
189 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
190 min(const ::VectorizedArray<Number, width> &,
191 const ::VectorizedArray<Number, width> &);
192 template <typename Number, size_t width>
194 pow(const ::VectorizedArray<Number, width> &, const Number p);
195 template <typename Number, size_t width>
197 sin(const ::VectorizedArray<Number, width> &);
198 template <typename Number, size_t width>
200 cos(const ::VectorizedArray<Number, width> &);
201 template <typename Number, size_t width>
203 tan(const ::VectorizedArray<Number, width> &);
204 template <typename Number, size_t width>
206 exp(const ::VectorizedArray<Number, width> &);
207 template <typename Number, size_t width>
209 log(const ::VectorizedArray<Number, width> &);
210} // namespace std
211
213
229namespace numbers
230{
234 static constexpr double E = 2.7182818284590452354;
235
239 static constexpr double LOG2E = 1.4426950408889634074;
240
244 static constexpr double LOG10E = 0.43429448190325182765;
245
249 static constexpr double LN2 = 0.69314718055994530942;
250
254 static constexpr double LN10 = 2.30258509299404568402;
255
259 static constexpr double PI = 3.14159265358979323846;
260
264 static constexpr double PI_2 = 1.57079632679489661923;
265
269 static constexpr double PI_4 = 0.78539816339744830962;
270
274 static constexpr double SQRT2 = 1.41421356237309504880;
275
279 static constexpr double SQRT1_2 = 0.70710678118654752440;
280
290 bool
291 is_finite(const double x);
292
297 bool
298 is_finite(const std::complex<double> &x);
299
304 bool
305 is_finite(const std::complex<float> &x);
306
315 bool
316 is_finite(const std::complex<long double> &x);
317
328 template <typename Number1, typename Number2>
329 constexpr DEAL_II_HOST_DEVICE bool
330 values_are_equal(const Number1 &value_1, const Number2 &value_2);
331
342 template <typename Number1, typename Number2>
343 constexpr bool
344 values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
345
353 template <typename Number>
354 constexpr DEAL_II_HOST_DEVICE bool
355 value_is_zero(const Number &value);
356
367 template <typename Number1, typename Number2>
368 bool
369 value_is_less_than(const Number1 &value_1, const Number2 &value_2);
370
381 template <typename Number1, typename Number2>
382 bool
383 value_is_less_than_or_equal_to(const Number1 &value_1,
384 const Number2 &value_2);
385
386
387
398 template <typename Number1, typename Number2>
399 bool
400 value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
401
412 template <typename Number1, typename Number2>
413 bool
414 value_is_greater_than_or_equal_to(const Number1 &value_1,
415 const Number2 &value_2);
416
425 template <typename number>
427 {
433 static constexpr bool is_complex = false;
434
441 using real_type = number;
442
446 using double_type = double;
447
455 static constexpr DEAL_II_HOST_DEVICE const number &
456 conjugate(const number &x);
457
466 static constexpr DEAL_II_HOST_DEVICE real_type
467 abs_square(const number &x);
468
472 static real_type
473 abs(const number &x);
474 };
475
476
481 template <typename number>
482 struct NumberTraits<std::complex<number>>
483 {
489 static constexpr bool is_complex = true;
490
497 using real_type = number;
498
502 using double_type = std::complex<double>;
503
507 static constexpr std::complex<number>
508 conjugate(const std::complex<number> &x);
509
516 static constexpr real_type
517 abs_square(const std::complex<number> &x);
518
519
523 static real_type
524 abs(const std::complex<number> &x);
525 };
526
527 // --------------- inline and template functions ---------------- //
528
529 inline bool
530 is_nan(const double x)
531 {
532 return std::isnan(x);
533 }
534
535
536
537 inline bool
538 is_finite(const double x)
539 {
540 return std::isfinite(x);
541 }
542
543
544
545 inline bool
546 is_finite(const std::complex<double> &x)
547 {
548 // Check complex numbers for infinity
549 // by testing real and imaginary part
550 return (is_finite(x.real()) && is_finite(x.imag()));
551 }
552
553
554
555 inline bool
556 is_finite(const std::complex<float> &x)
557 {
558 // Check complex numbers for infinity
559 // by testing real and imaginary part
560 return (is_finite(x.real()) && is_finite(x.imag()));
561 }
562
563
564
565 inline bool
566 is_finite(const std::complex<long double> &x)
567 {
568 // Same for std::complex<long double>
569 return (is_finite(x.real()) && is_finite(x.imag()));
570 }
571
572
573 template <typename number>
574 constexpr DEAL_II_HOST_DEVICE const number &
576 {
577 return x;
578 }
579
580
581
582 template <typename number>
585 {
586 return x * x;
587 }
588
589
590
591 template <typename number>
594 {
595 // Make things work with AD types
596 using std::abs;
597#ifdef DEAL_II_WITH_ADOLC
598 // This one is a little tricky - we have our own abs function in ::,
599 // prototyped with forward-declared types in this file, but it only exists
600 // if we have ADOL-C: hence we only add this using statement in that
601 // situation
602 using ::abs;
603#endif
604 return abs(x);
605 }
606
607
608
609 template <typename number>
610 constexpr std::complex<number>
611 NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
612 {
613 return std::conj(x);
614 }
615
616
617
618 template <typename number>
619 typename NumberTraits<std::complex<number>>::real_type
620 NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
621 {
622 // Make things work with AD types
623 using std::abs;
624#ifdef DEAL_II_WITH_ADOLC
625 // Same comment as the non-complex case holds here
626 using ::abs;
627#endif
628 return abs(x);
629 }
630
631
632
633 template <typename number>
634 constexpr typename NumberTraits<std::complex<number>>::real_type
635 NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
636 {
637 return std::norm(x);
638 }
639
640} // namespace numbers
641
642
643// Forward declarations
645{
646 namespace AD
647 {
648 namespace internal
649 {
650 // Defined in differentiation/ad/ad_number_traits.h
651 template <typename T>
653 } // namespace internal
654
655 // Defined in differentiation/ad/ad_number_traits.h
656 template <typename NumberType>
658 } // namespace AD
659} // namespace Differentiation
660
661
662namespace internal
663{
668 template <typename From, typename To>
670 {
671 // Source: https://stackoverflow.com/a/16944130
672 private:
673 template <typename T>
674 static void f(T);
675
676 template <typename F, typename T>
677 static constexpr auto
678 test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
679 {
680 return true;
681 }
682
683 template <typename F, typename T>
684 static constexpr auto
685 test(...) -> bool
686 {
687 return false;
688 }
689
690 public:
691 static const bool value = test<From, To>(0);
692 };
693
694 /*
695 * The structs below are needed to convert between some special number types.
696 * Also see tensor.h for another specialization.
697 */
698 template <typename T>
700 {
701 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T &
702 value(const T &t)
703 {
704 return t;
705 }
706
707 // Below are generic functions that allows an overload for any
708 // type U that is transformable to type T. This is particularly
709 // useful when needing to cast exotic number types
710 // (e.g. auto-differentiable or symbolic numbers) to a floating
711 // point one, such as might happen when converting between tensor
712 // types.
713
714 // Type T is constructible from F.
715 template <typename F>
716 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T
717 value(const F &f,
718 std::enable_if_t<!std::is_same_v<std::decay_t<T>, std::decay_t<F>> &&
719 std::is_constructible_v<T, F>> * = nullptr)
720 {
721 return T(f);
722 }
723
724 // Type T is explicitly convertible (but not constructible) from F.
725 template <typename F>
726 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T
727 value(const F &f,
728 std::enable_if_t<!std::is_same_v<std::decay_t<T>, std::decay_t<F>> &&
729 !std::is_constructible_v<T, F> &&
731 nullptr)
732 {
733 return static_cast<T>(f);
734 }
735
736 // Sacado doesn't provide any conversion operators, so we have
737 // to extract the value and perform further conversions from there.
738 // To be safe, we extend this to other possible AD numbers that
739 // might fall into the same category.
740 template <typename F>
741 static T
743 const F &f,
744 std::enable_if_t<!std::is_same_v<std::decay_t<T>, std::decay_t<F>> &&
745 !std::is_constructible_v<T, F> &&
748 {
750 }
751 };
752
753 template <typename T>
754 struct NumberType<std::complex<T>>
755 {
756 static constexpr const std::complex<T> &
757 value(const std::complex<T> &t)
758 {
759 return t;
760 }
761
762 static constexpr std::complex<T>
763 value(const T &t)
764 {
765 return std::complex<T>(t);
766 }
767
768 // Facilitate cast from complex<double> to complex<float>
769 template <typename U>
770 static constexpr std::complex<T>
771 value(const std::complex<U> &t)
772 {
773 return std::complex<T>(NumberType<T>::value(t.real()),
774 NumberType<T>::value(t.imag()));
775 }
776 };
777
778#ifdef DEAL_II_WITH_CUDA
779 template <>
780 struct NumberType<cuComplex>
781 {
782 static cuComplex
783 value(const float t)
784 {
785 return make_cuComplex(t, 0.f);
786 }
787 };
788
789 template <>
790 struct NumberType<cuDoubleComplex>
791 {
792 static cuDoubleComplex
793 value(const double t)
794 {
795 return make_cuDoubleComplex(t, 0.);
796 }
797 };
798#endif
799} // namespace internal
800
801namespace numbers
802{
803#ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
804
815 // Defined in differentiation/ad/adolc_number_types.cc
816 bool
817 values_are_equal(const adouble &value_1, const adouble &value_2);
818
819
830 template <typename Number>
831 bool
832 values_are_equal(const adouble &value_1, const Number &value_2)
833 {
834 // Use the specialized definition for two ADOL-C taped types
835 return values_are_equal(
837 }
838
839
850 template <typename Number>
851 bool
852 values_are_equal(const Number &value_1, const adouble &value_2)
853 {
854 // Use the above definition
855 return values_are_equal(value_2, value_1);
856 }
857
869 // Defined in differentiation/ad/adolc_number_types.cc
870 bool
871 value_is_less_than(const adouble &value_1, const adouble &value_2);
872
873
885 template <typename Number>
886 bool
887 value_is_less_than(const adouble &value_1, const Number &value_2)
888 {
889 // Use the specialized definition for two ADOL-C taped types
890 return value_is_less_than(
892 }
893
894
906 template <typename Number>
907 bool
908 value_is_less_than(const Number &value_1, const adouble &value_2)
909 {
910 // Use the specialized definition for two ADOL-C taped types
911 return value_is_less_than(
913 }
914
915#endif
916
917
918 template <typename Number1, typename Number2>
919 constexpr DEAL_II_HOST_DEVICE bool
920 values_are_equal(const Number1 &value_1, const Number2 &value_2)
921 {
922 return (value_1 == ::internal::NumberType<Number1>::value(value_2));
923 }
924
925
926 template <typename Number1, typename Number2>
927 inline constexpr bool
928 values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
929 {
930 return !(values_are_equal(value_1, value_2));
931 }
932
933
934 template <typename Number>
935 constexpr DEAL_II_HOST_DEVICE bool
936 value_is_zero(const Number &value)
937 {
938 return values_are_equal(value, 0.0);
939 }
940
941
942 template <typename Number1, typename Number2>
943 inline bool
944 value_is_less_than(const Number1 &value_1, const Number2 &value_2)
945 {
946 return (value_1 < ::internal::NumberType<Number1>::value(value_2));
947 }
948
949
950 template <typename Number1, typename Number2>
951 inline bool
952 value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
953 {
954 return (value_is_less_than(value_1, value_2) ||
955 values_are_equal(value_1, value_2));
956 }
957
958
959 template <typename Number1, typename Number2>
960 bool
961 value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
962 {
963 return !(value_is_less_than_or_equal_to(value_1, value_2));
964 }
965
966
967 template <typename Number1, typename Number2>
968 inline bool
969 value_is_greater_than_or_equal_to(const Number1 &value_1,
970 const Number2 &value_2)
971 {
972 return !(value_is_less_than(value_1, value_2));
973 }
974} // namespace numbers
975
977
978#endif
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
Definition numbers.h:53
static constexpr double LOG10E
Definition numbers.h:244
static constexpr double PI_2
Definition numbers.h:264
bool value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:952
static constexpr double E
Definition numbers.h:234
static constexpr double PI
Definition numbers.h:259
bool value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:961
constexpr bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:928
static constexpr double SQRT2
Definition numbers.h:274
constexpr bool value_is_zero(const Number &value)
Definition numbers.h:936
static constexpr double SQRT1_2
Definition numbers.h:279
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:920
static constexpr double PI_4
Definition numbers.h:269
static constexpr double LN10
Definition numbers.h:254
static constexpr double LN2
Definition numbers.h:249
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:944
bool is_finite(const double x)
Definition numbers.h:538
static constexpr double LOG2E
Definition numbers.h:239
bool value_is_greater_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:969
bool is_nan(const double x)
Definition numbers.h:530
STL namespace.
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > tan(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
#define DEAL_II_HOST_DEVICE
Definition numbers.h:34
#define DEAL_II_HOST_DEVICE_ALWAYS_INLINE
Definition numbers.h:36
static cuComplex value(const float t)
Definition numbers.h:783
static cuDoubleComplex value(const double t)
Definition numbers.h:793
static constexpr std::complex< T > value(const std::complex< U > &t)
Definition numbers.h:771
static constexpr std::complex< T > value(const T &t)
Definition numbers.h:763
static constexpr const std::complex< T > & value(const std::complex< T > &t)
Definition numbers.h:757
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T value(const F &f, std::enable_if_t<!std::is_same_v< std::decay_t< T >, std::decay_t< F > > &&std::is_constructible_v< T, F > > *=nullptr)
Definition numbers.h:717
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T value(const F &f, std::enable_if_t<!std::is_same_v< std::decay_t< T >, std::decay_t< F > > &&!std::is_constructible_v< T, F > &&is_explicitly_convertible< const F, T >::value > *=nullptr)
Definition numbers.h:727
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition numbers.h:702
static T value(const F &f, std::enable_if_t<!std::is_same_v< std::decay_t< T >, std::decay_t< F > > &&!std::is_constructible_v< T, F > &&!is_explicitly_convertible< const F, T >::value &&Differentiation::AD::is_ad_number< F >::value > *=nullptr)
Definition numbers.h:742
static constexpr unsigned int max_width
Definition numbers.h:84
static constexpr auto test(...) -> bool
Definition numbers.h:685
static constexpr auto test(int) -> decltype(f(static_cast< T >(std::declval< F >())), true)
Definition numbers.h:678
static constexpr const number & conjugate(const number &x)
Definition numbers.h:575
static constexpr bool is_complex
Definition numbers.h:433
static real_type abs(const number &x)
Definition numbers.h:593
static constexpr real_type abs_square(const number &x)
Definition numbers.h:584