Reference documentation for deal.II version 9.5.0
\(\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// Copyright (C) 2006 - 2023 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_numbers_h
17#define dealii_numbers_h
18
19
20#include <deal.II/base/config.h>
21
22#include <deal.II/base/types.h>
23
24#ifdef DEAL_II_WITH_CUDA
25# include <cuComplex.h>
26#endif
27
28#include <Kokkos_Macros.hpp>
29
30#include <cmath>
31#include <complex>
32#include <cstddef>
33#include <type_traits>
34
35#define DEAL_II_HOST_DEVICE KOKKOS_FUNCTION
36#define DEAL_II_CUDA_HOST_DEV DEAL_II_HOST_DEVICE
37#define DEAL_II_HOST_DEVICE_ALWAYS_INLINE KOKKOS_FORCEINLINE_FUNCTION
38
39// clang++ assumes that all constexpr functions are __host__ __device__ when
40// compiling CUDA code, i.e, when Kokkos was configured with CUDA support.
41// This is problematic when calling non-constexpr functions in constexpr
42// functions. Hence, we need a way to annotate functions explicitly as
43// host-only.
44#if defined(__clang__) && defined(__CUDA__)
45# define DEAL_II_HOST __host__
46#else
47# define DEAL_II_HOST
48#endif
49
50// Forward-declare the automatic differentiation types so we can add prototypes
51// for our own wrappers.
52#ifdef DEAL_II_WITH_ADOLC
53class adouble;
54namespace adtl
55{
56 class adouble;
57}
58#endif
59
61
62namespace internal
63{
80 template <typename Number>
82 {
86 constexpr static unsigned int max_width = 1;
87 };
88
95 template <>
97 {
101 constexpr static unsigned int max_width =
102#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512
103 8;
104#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256
105 4;
106#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128
107 2;
108#else
109 1;
110#endif
111 };
112
119 template <>
121 {
125 constexpr static unsigned int max_width =
126#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ALTIVEC__)
127 4;
128#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && defined(__AVX512F__)
129 16;
130#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)
131 8;
132#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__)
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 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 bool const 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<typename std::decay<T>::type,
719 typename std::decay<F>::type>::value &&
720 std::is_constructible<T, F>::value> * = nullptr)
721 {
722 return T(f);
723 }
724
725 // Type T is explicitly convertible (but not constructible) from F.
726 template <typename F>
727 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T
728 value(const F &f,
729 std::enable_if_t<!std::is_same<typename std::decay<T>::type,
730 typename std::decay<F>::type>::value &&
731 !std::is_constructible<T, F>::value &&
733 nullptr)
734 {
735 return static_cast<T>(f);
736 }
737
738 // Sacado doesn't provide any conversion operators, so we have
739 // to extract the value and perform further conversions from there.
740 // To be safe, we extend this to other possible AD numbers that
741 // might fall into the same category.
742 template <typename F>
743 static T
745 const F &f,
746 std::enable_if_t<!std::is_same<typename std::decay<T>::type,
747 typename std::decay<F>::type>::value &&
748 !std::is_constructible<T, F>::value &&
751 {
753 }
754 };
755
756 template <typename T>
757 struct NumberType<std::complex<T>>
758 {
759 static constexpr const std::complex<T> &
760 value(const std::complex<T> &t)
761 {
762 return t;
763 }
764
765 static constexpr std::complex<T>
766 value(const T &t)
767 {
768 return std::complex<T>(t);
769 }
770
771 // Facilitate cast from complex<double> to complex<float>
772 template <typename U>
773 static constexpr std::complex<T>
774 value(const std::complex<U> &t)
775 {
776 return std::complex<T>(NumberType<T>::value(t.real()),
777 NumberType<T>::value(t.imag()));
778 }
779 };
780
781#ifdef DEAL_II_WITH_CUDA
782 template <>
783 struct NumberType<cuComplex>
784 {
785 static cuComplex
786 value(const float t)
787 {
788 return make_cuComplex(t, 0.f);
789 }
790 };
791
792 template <>
793 struct NumberType<cuDoubleComplex>
794 {
795 static cuDoubleComplex
796 value(const double t)
797 {
798 return make_cuDoubleComplex(t, 0.);
799 }
800 };
801#endif
802} // namespace internal
803
804namespace numbers
805{
806#ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
807
818 // Defined in differentiation/ad/adolc_number_types.cc
819 bool
820 values_are_equal(const adouble &value_1, const adouble &value_2);
821
822
833 template <typename Number>
834 bool
835 values_are_equal(const adouble &value_1, const Number &value_2)
836 {
837 // Use the specialized definition for two ADOL-C taped types
838 return values_are_equal(
840 }
841
842
853 template <typename Number>
854 bool
855 values_are_equal(const Number &value_1, const adouble &value_2)
856 {
857 // Use the above definition
858 return values_are_equal(value_2, value_1);
859 }
860
872 // Defined in differentiation/ad/adolc_number_types.cc
873 bool
874 value_is_less_than(const adouble &value_1, const adouble &value_2);
875
876
888 template <typename Number>
889 bool
890 value_is_less_than(const adouble &value_1, const Number &value_2)
891 {
892 // Use the specialized definition for two ADOL-C taped types
893 return value_is_less_than(
895 }
896
897
909 template <typename Number>
910 bool
911 value_is_less_than(const Number &value_1, const adouble &value_2)
912 {
913 // Use the specialized definition for two ADOL-C taped types
914 return value_is_less_than(
916 }
917
918#endif
919
920
921 template <typename Number1, typename Number2>
922 constexpr DEAL_II_HOST_DEVICE bool
923 values_are_equal(const Number1 &value_1, const Number2 &value_2)
924 {
925 return (value_1 == ::internal::NumberType<Number1>::value(value_2));
926 }
927
928
929 template <typename Number1, typename Number2>
930 inline bool
931 values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
932 {
933 return !(values_are_equal(value_1, value_2));
934 }
935
936
937 template <typename Number>
938 constexpr DEAL_II_HOST_DEVICE bool
939 value_is_zero(const Number &value)
940 {
941 return values_are_equal(value, 0.0);
942 }
943
944
945 template <typename Number1, typename Number2>
946 inline bool
947 value_is_less_than(const Number1 &value_1, const Number2 &value_2)
948 {
949 return (value_1 < ::internal::NumberType<Number1>::value(value_2));
950 }
951
952
953 template <typename Number1, typename Number2>
954 inline bool
955 value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
956 {
957 return (value_is_less_than(value_1, value_2) ||
958 values_are_equal(value_1, value_2));
959 }
960
961
962 template <typename Number1, typename Number2>
963 bool
964 value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
965 {
966 return !(value_is_less_than_or_equal_to(value_1, value_2));
967 }
968
969
970 template <typename Number1, typename Number2>
971 inline bool
972 value_is_greater_than_or_equal_to(const Number1 &value_1,
973 const Number2 &value_2)
974 {
975 return !(value_is_less_than(value_1, value_2));
976 }
977} // namespace numbers
978
980
981#endif
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Definition numbers.h:55
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:955
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:964
static constexpr double SQRT2
Definition numbers.h:274
constexpr bool value_is_zero(const Number &value)
Definition numbers.h:939
bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:931
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:923
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:947
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:972
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 > 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:35
#define DEAL_II_HOST_DEVICE_ALWAYS_INLINE
Definition numbers.h:37
static cuComplex value(const float t)
Definition numbers.h:786
static cuDoubleComplex value(const double t)
Definition numbers.h:796
static constexpr std::complex< T > value(const std::complex< U > &t)
Definition numbers.h:774
static constexpr std::complex< T > value(const T &t)
Definition numbers.h:766
static constexpr const std::complex< T > & value(const std::complex< T > &t)
Definition numbers.h:760
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T value(const F &f, std::enable_if_t<!std::is_same< typename std::decay< T >::type, typename std::decay< F >::type >::value &&!std::is_constructible< T, F >::value &&is_explicitly_convertible< const F, T >::value > *=nullptr)
Definition numbers.h:728
static T value(const F &f, std::enable_if_t<!std::is_same< typename std::decay< T >::type, typename std::decay< F >::type >::value &&!std::is_constructible< T, F >::value &&!is_explicitly_convertible< const F, T >::value &&Differentiation::AD::is_ad_number< F >::value > *=nullptr)
Definition numbers.h:744
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T value(const F &f, std::enable_if_t<!std::is_same< typename std::decay< T >::type, typename std::decay< F >::type >::value &&std::is_constructible< T, F >::value > *=nullptr)
Definition numbers.h:717
static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T & value(const T &t)
Definition numbers.h:702
static constexpr unsigned int max_width
Definition numbers.h:86
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