Reference documentation for deal.II version GIT eef19498cf 2022-10-05 14:45: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\}}\)
numbers.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2006 - 2022 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_COMPILER_CUDA_AWARE
25 # include <cuComplex.h>
26 #endif
27 
28 #include <cmath>
29 #include <complex>
30 #include <cstddef>
31 #include <type_traits>
32 
33 #ifdef DEAL_II_COMPILER_CUDA_AWARE
34 # define DEAL_II_CUDA_HOST_DEV __host__ __device__
35 #else
36 # define DEAL_II_CUDA_HOST_DEV
37 #endif
38 
39 // Forward-declare the automatic differentiation types so we can add prototypes
40 // for our own wrappers.
41 #ifdef DEAL_II_WITH_ADOLC
42 class adouble;
43 namespace adtl
44 {
45  class adouble;
46 }
47 #endif
48 
50 
51 namespace internal
52 {
69  template <typename Number>
71  {
75  constexpr static unsigned int max_width = 1;
76  };
77 
84  template <>
86  {
90  constexpr static unsigned int max_width =
91 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512
92  8;
93 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256
94  4;
95 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128
96  2;
97 #else
98  1;
99 #endif
100  };
101 
108  template <>
110  {
114  constexpr static unsigned int max_width =
115 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ALTIVEC__)
116  4;
117 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && defined(__AVX512F__)
118  16;
119 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)
120  8;
121 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__)
122  4;
123 #else
124  1;
125 #endif
126  };
127 
128 
129 } // namespace internal
130 
131 // forward declarations to support abs or sqrt operations on VectorizedArray
132 #ifndef DOXYGEN
133 template <typename Number,
134  std::size_t width =
136 class VectorizedArray;
137 template <typename T>
138 struct EnableIfScalar;
139 #endif
140 
141 #ifdef DEAL_II_WITH_ADOLC
142 # ifndef DOXYGEN
143 // Prototype some inline functions present in adolc_math.h for use in
144 // NumberTraits.
145 //
146 // ADOL-C uses fabs(), but for genericity we want to use abs(). Simultaneously,
147 // though, we don't want to include ADOL-C headers in this header since
148 // numbers.h is in everything. To get around this: use C++ rules which permit
149 // the use of forward-declared classes in function prototypes to declare some
150 // functions which are defined in adolc_math.h. This permits us to write "using
151 // ::abs;" in NumberTraits which will allow us to select the correct
152 // overload (the one in ::) when instantiating NumberTraits for ADOL-C
153 // types.
154 
155 adouble
156 abs(const adouble &x);
157 
158 adtl::adouble
159 abs(const adtl::adouble &x);
160 # endif
161 #endif
162 
164 
165 namespace std
166 {
167  template <typename Number, std::size_t width>
168  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
169  sqrt(const ::VectorizedArray<Number, width> &);
170  template <typename Number, std::size_t width>
171  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
172  abs(const ::VectorizedArray<Number, width> &);
173  template <typename Number, std::size_t width>
174  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
175  max(const ::VectorizedArray<Number, width> &,
176  const ::VectorizedArray<Number, width> &);
177  template <typename Number, std::size_t width>
178  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
179  min(const ::VectorizedArray<Number, width> &,
180  const ::VectorizedArray<Number, width> &);
181  template <typename Number, size_t width>
183  pow(const ::VectorizedArray<Number, width> &, const Number p);
184  template <typename Number, size_t width>
186  sin(const ::VectorizedArray<Number, width> &);
187  template <typename Number, size_t width>
189  cos(const ::VectorizedArray<Number, width> &);
190  template <typename Number, size_t width>
192  tan(const ::VectorizedArray<Number, width> &);
193  template <typename Number, size_t width>
195  exp(const ::VectorizedArray<Number, width> &);
196  template <typename Number, size_t width>
198  log(const ::VectorizedArray<Number, width> &);
199 } // namespace std
200 
202 
218 namespace numbers
219 {
223  static constexpr double E = 2.7182818284590452354;
224 
228  static constexpr double LOG2E = 1.4426950408889634074;
229 
233  static constexpr double LOG10E = 0.43429448190325182765;
234 
238  static constexpr double LN2 = 0.69314718055994530942;
239 
243  static constexpr double LN10 = 2.30258509299404568402;
244 
248  static constexpr double PI = 3.14159265358979323846;
249 
253  static constexpr double PI_2 = 1.57079632679489661923;
254 
258  static constexpr double PI_4 = 0.78539816339744830962;
259 
263  static constexpr double SQRT2 = 1.41421356237309504880;
264 
268  static constexpr double SQRT1_2 = 0.70710678118654752440;
269 
275  template <typename Number, typename = void>
276  struct is_cuda_compatible : std::true_type
277  {};
278 
282  template <typename Number>
283  struct is_cuda_compatible<std::complex<Number>, void> : std::false_type
284  {};
285 
295  bool
296  is_finite(const double x);
297 
302  bool
303  is_finite(const std::complex<double> &x);
304 
309  bool
310  is_finite(const std::complex<float> &x);
311 
320  bool
321  is_finite(const std::complex<long double> &x);
322 
333  template <typename Number1, typename Number2>
334  constexpr bool
335  values_are_equal(const Number1 &value_1, const Number2 &value_2);
336 
347  template <typename Number1, typename Number2>
348  bool
349  values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
350 
358  template <typename Number>
359  constexpr bool
360  value_is_zero(const Number &value);
361 
372  template <typename Number1, typename Number2>
373  bool
374  value_is_less_than(const Number1 &value_1, const Number2 &value_2);
375 
386  template <typename Number1, typename Number2>
387  bool
388  value_is_less_than_or_equal_to(const Number1 &value_1,
389  const Number2 &value_2);
390 
391 
392 
403  template <typename Number1, typename Number2>
404  bool
405  value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
406 
417  template <typename Number1, typename Number2>
418  bool
419  value_is_greater_than_or_equal_to(const Number1 &value_1,
420  const Number2 &value_2);
421 
430  template <typename number>
432  {
438  static constexpr bool is_complex = false;
439 
446  using real_type = number;
447 
452 
460  static constexpr DEAL_II_CUDA_HOST_DEV const number &
461  conjugate(const number &x);
462 
471  template <typename Dummy = number>
472  static constexpr DEAL_II_CUDA_HOST_DEV
473  std::enable_if_t<std::is_same<Dummy, number>::value &&
475  real_type>
476  abs_square(const number &x);
477 
478  template <typename Dummy = number>
479  static constexpr std::enable_if_t<std::is_same<Dummy, number>::value &&
481  real_type>
482  abs_square(const number &x);
483 
487  static real_type
488  abs(const number &x);
489  };
490 
491 
496  template <typename number>
497  struct NumberTraits<std::complex<number>>
498  {
504  static constexpr bool is_complex = true;
505 
512  using real_type = number;
513 
517  using double_type = std::complex<double>;
518 
522  static constexpr std::complex<number>
523  conjugate(const std::complex<number> &x);
524 
531  static constexpr real_type
532  abs_square(const std::complex<number> &x);
533 
534 
538  static real_type
539  abs(const std::complex<number> &x);
540  };
541 
542  // --------------- inline and template functions ---------------- //
543 
544  inline bool
545  is_nan(const double x)
546  {
547  return std::isnan(x);
548  }
549 
550 
551 
552  inline bool
553  is_finite(const double x)
554  {
555  return std::isfinite(x);
556  }
557 
558 
559 
560  inline bool
561  is_finite(const std::complex<double> &x)
562  {
563  // Check complex numbers for infinity
564  // by testing real and imaginary part
565  return (is_finite(x.real()) && is_finite(x.imag()));
566  }
567 
568 
569 
570  inline bool
571  is_finite(const std::complex<float> &x)
572  {
573  // Check complex numbers for infinity
574  // by testing real and imaginary part
575  return (is_finite(x.real()) && is_finite(x.imag()));
576  }
577 
578 
579 
580  inline bool
581  is_finite(const std::complex<long double> &x)
582  {
583  // Same for std::complex<long double>
584  return (is_finite(x.real()) && is_finite(x.imag()));
585  }
586 
587 
588  template <typename number>
589  constexpr DEAL_II_CUDA_HOST_DEV const number &
591  {
592  return x;
593  }
594 
595 
596 
597  template <typename number>
598  template <typename Dummy>
599  constexpr DEAL_II_CUDA_HOST_DEV
600  std::enable_if_t<std::is_same<Dummy, number>::value &&
604  {
605  return x * x;
606  }
607 
608 
609 
610  template <typename number>
611  template <typename Dummy>
612  constexpr std::enable_if_t<std::is_same<Dummy, number>::value &&
616  {
617  return x * x;
618  }
619 
620 
621 
622  template <typename number>
624  NumberTraits<number>::abs(const number &x)
625  {
626  // Make things work with AD types
627  using std::abs;
628 #ifdef DEAL_II_WITH_ADOLC
629  // This one is a little tricky - we have our own abs function in ::,
630  // prototyped with forward-declared types in this file, but it only exists
631  // if we have ADOL-C: hence we only add this using statement in that
632  // situation
633  using ::abs;
634 #endif
635  return abs(x);
636  }
637 
638 
639 
640  template <typename number>
641  constexpr std::complex<number>
642  NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
643  {
644  return std::conj(x);
645  }
646 
647 
648 
649  template <typename number>
650  typename NumberTraits<std::complex<number>>::real_type
651  NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
652  {
653  // Make things work with AD types
654  using std::abs;
655 #ifdef DEAL_II_WITH_ADOLC
656  // Same comment as the non-complex case holds here
657  using ::abs;
658 #endif
659  return abs(x);
660  }
661 
662 
663 
664  template <typename number>
665  constexpr typename NumberTraits<std::complex<number>>::real_type
666  NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
667  {
668  return std::norm(x);
669  }
670 
671 } // namespace numbers
672 
673 
674 // Forward declarations
676 {
677  namespace AD
678  {
679  namespace internal
680  {
681  // Defined in differentiation/ad/ad_number_traits.h
682  template <typename T>
683  struct NumberType;
684  } // namespace internal
685 
686  // Defined in differentiation/ad/ad_number_traits.h
687  template <typename NumberType>
688  struct is_ad_number;
689  } // namespace AD
690 } // namespace Differentiation
691 
692 
693 namespace internal
694 {
699  template <typename From, typename To>
701  {
702  // Source: https://stackoverflow.com/a/16944130
703  private:
704  template <typename T>
705  static void f(T);
706 
707  template <typename F, typename T>
708  static constexpr auto
709  test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
710  {
711  return true;
712  }
713 
714  template <typename F, typename T>
715  static constexpr auto
716  test(...) -> bool
717  {
718  return false;
719  }
720 
721  public:
722  static bool const value = test<From, To>(0);
723  };
724 
725  /*
726  * The structs below are needed to convert between some special number types.
727  * Also see tensor.h for another specialization.
728  */
729  template <typename T>
730  struct NumberType
731  {
732  static constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV const T &
733  value(const T &t)
734  {
735  return t;
736  }
737 
738  // Below are generic functions that allows an overload for any
739  // type U that is transformable to type T. This is particularly
740  // useful when needing to cast exotic number types
741  // (e.g. auto-differentiable or symbolic numbers) to a floating
742  // point one, such as might happen when converting between tensor
743  // types.
744 
745  // Type T is constructible from F.
746  template <typename F>
748  value(const F &f,
749  std::enable_if_t<!std::is_same<typename std::decay<T>::type,
750  typename std::decay<F>::type>::value &&
751  std::is_constructible<T, F>::value> * = nullptr)
752  {
753  return T(f);
754  }
755 
756  // Type T is explicitly convertible (but not constructible) from F.
757  template <typename F>
758  static constexpr DEAL_II_ALWAYS_INLINE T
759  value(const F &f,
760  std::enable_if_t<!std::is_same<typename std::decay<T>::type,
761  typename std::decay<F>::type>::value &&
762  !std::is_constructible<T, F>::value &&
764  nullptr)
765  {
766  return static_cast<T>(f);
767  }
768 
769  // Sacado doesn't provide any conversion operators, so we have
770  // to extract the value and perform further conversions from there.
771  // To be safe, we extend this to other possible AD numbers that
772  // might fall into the same category.
773  template <typename F>
774  static T
776  const F &f,
777  std::enable_if_t<!std::is_same<typename std::decay<T>::type,
778  typename std::decay<F>::type>::value &&
779  !std::is_constructible<T, F>::value &&
782  {
784  }
785  };
786 
787  template <typename T>
788  struct NumberType<std::complex<T>>
789  {
790  static constexpr const std::complex<T> &
791  value(const std::complex<T> &t)
792  {
793  return t;
794  }
795 
796  static constexpr std::complex<T>
797  value(const T &t)
798  {
799  return std::complex<T>(t);
800  }
801 
802  // Facilitate cast from complex<double> to complex<float>
803  template <typename U>
804  static constexpr std::complex<T>
805  value(const std::complex<U> &t)
806  {
807  return std::complex<T>(NumberType<T>::value(t.real()),
808  NumberType<T>::value(t.imag()));
809  }
810  };
811 
812 #ifdef DEAL_II_COMPILER_CUDA_AWARE
813  template <>
814  struct NumberType<cuComplex>
815  {
816  static cuComplex
817  value(const float t)
818  {
819  return make_cuComplex(t, 0.f);
820  }
821  };
822 
823  template <>
824  struct NumberType<cuDoubleComplex>
825  {
826  static cuDoubleComplex
827  value(const double t)
828  {
829  return make_cuDoubleComplex(t, 0.);
830  }
831  };
832 #endif
833 } // namespace internal
834 
835 namespace numbers
836 {
837 #ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
838 
849  // Defined in differentiation/ad/adolc_number_types.cc
850  bool
851  values_are_equal(const adouble &value_1, const adouble &value_2);
852 
853 
864  template <typename Number>
865  bool
866  values_are_equal(const adouble &value_1, const Number &value_2)
867  {
868  // Use the specialized definition for two ADOL-C taped types
869  return values_are_equal(
870  value_1, ::internal::NumberType<adouble>::value(value_2));
871  }
872 
873 
884  template <typename Number>
885  bool
886  values_are_equal(const Number &value_1, const adouble &value_2)
887  {
888  // Use the above definition
889  return values_are_equal(value_2, value_1);
890  }
891 
903  // Defined in differentiation/ad/adolc_number_types.cc
904  bool
905  value_is_less_than(const adouble &value_1, const adouble &value_2);
906 
907 
919  template <typename Number>
920  bool
921  value_is_less_than(const adouble &value_1, const Number &value_2)
922  {
923  // Use the specialized definition for two ADOL-C taped types
924  return value_is_less_than(
925  value_1, ::internal::NumberType<adouble>::value(value_2));
926  }
927 
928 
940  template <typename Number>
941  bool
942  value_is_less_than(const Number &value_1, const adouble &value_2)
943  {
944  // Use the specialized definition for two ADOL-C taped types
945  return value_is_less_than(
946  ::internal::NumberType<adouble>::value(value_1), value_2);
947  }
948 
949 #endif
950 
951 
952  template <typename Number1, typename Number2>
953  constexpr bool
954  values_are_equal(const Number1 &value_1, const Number2 &value_2)
955  {
956  return (value_1 == ::internal::NumberType<Number1>::value(value_2));
957  }
958 
959 
960  template <typename Number1, typename Number2>
961  inline bool
962  values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
963  {
964  return !(values_are_equal(value_1, value_2));
965  }
966 
967 
968  template <typename Number>
969  constexpr bool
970  value_is_zero(const Number &value)
971  {
972  return values_are_equal(value, 0.0);
973  }
974 
975 
976  template <typename Number1, typename Number2>
977  inline bool
978  value_is_less_than(const Number1 &value_1, const Number2 &value_2)
979  {
980  return (value_1 < ::internal::NumberType<Number1>::value(value_2));
981  }
982 
983 
984  template <typename Number1, typename Number2>
985  inline bool
986  value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
987  {
988  return (value_is_less_than(value_1, value_2) ||
989  values_are_equal(value_1, value_2));
990  }
991 
992 
993  template <typename Number1, typename Number2>
994  bool
995  value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
996  {
997  return !(value_is_less_than_or_equal_to(value_1, value_2));
998  }
999 
1000 
1001  template <typename Number1, typename Number2>
1002  inline bool
1003  value_is_greater_than_or_equal_to(const Number1 &value_1,
1004  const Number2 &value_2)
1005  {
1006  return !(value_is_less_than(value_1, value_2));
1007  }
1008 } // namespace numbers
1009 
1011 
1012 #endif
VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &x)
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:102
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:457
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:458
static const char T
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
Definition: numbers.h:44
static constexpr double LOG10E
Definition: numbers.h:233
static constexpr double PI_2
Definition: numbers.h:253
bool value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:986
static constexpr double E
Definition: numbers.h:223
static constexpr double PI
Definition: numbers.h:248
bool value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:995
static constexpr double SQRT2
Definition: numbers.h:263
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:970
bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:962
static constexpr double SQRT1_2
Definition: numbers.h:268
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:954
static constexpr double PI_4
Definition: numbers.h:258
static constexpr double LN10
Definition: numbers.h:243
static constexpr double LN2
Definition: numbers.h:238
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:978
bool is_finite(const double x)
Definition: numbers.h:553
static constexpr double LOG2E
Definition: numbers.h:228
bool value_is_greater_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:1003
bool is_nan(const double x)
Definition: numbers.h:545
#define DEAL_II_CUDA_HOST_DEV
Definition: numbers.h:34
static cuComplex value(const float t)
Definition: numbers.h:817
static cuDoubleComplex value(const double t)
Definition: numbers.h:827
static constexpr std::complex< T > value(const std::complex< U > &t)
Definition: numbers.h:805
static constexpr std::complex< T > value(const T &t)
Definition: numbers.h:797
static constexpr const std::complex< T > & value(const std::complex< T > &t)
Definition: numbers.h:791
static constexpr 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:748
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:775
static constexpr const T & value(const T &t)
Definition: numbers.h:733
static constexpr 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:759
constexpr static unsigned int max_width
Definition: numbers.h:75
static constexpr auto test(...) -> bool
Definition: numbers.h:716
static constexpr auto test(int) -> decltype(f(static_cast< T >(std::declval< F >())), true)
Definition: numbers.h:709
static constexpr const number & conjugate(const number &x)
Definition: numbers.h:590
static constexpr std::enable_if_t< std::is_same< Dummy, number >::value &&!is_cuda_compatible< Dummy >::value, real_type > abs_square(const number &x)
static constexpr bool is_complex
Definition: numbers.h:438
static real_type abs(const number &x)
Definition: numbers.h:624
static constexpr std::enable_if_t< std::is_same< Dummy, number >::value &&is_cuda_compatible< Dummy >::value, real_type > abs_square(const number &x)