Reference documentation for deal.II version GIT fac6d58a2d 2022-05-23 17:55: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 
40 
41 namespace internal
42 {
59  template <typename Number>
61  {
65  constexpr static unsigned int max_width = 1;
66  };
67 
74  template <>
76  {
80  constexpr static unsigned int max_width =
81 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512
82  8;
83 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256
84  4;
85 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128
86  2;
87 #else
88  1;
89 #endif
90  };
91 
98  template <>
100  {
104  constexpr static unsigned int max_width =
105 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ALTIVEC__)
106  4;
107 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && defined(__AVX512F__)
108  16;
109 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)
110  8;
111 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__)
112  4;
113 #else
114  1;
115 #endif
116  };
117 
118 
119 } // namespace internal
120 
121 // forward declarations to support abs or sqrt operations on VectorizedArray
122 #ifndef DOXYGEN
123 template <typename Number,
124  std::size_t width =
126 class VectorizedArray;
127 template <typename T>
128 struct EnableIfScalar;
129 #endif
130 
132 
133 // Declare / Import auto-differentiable math functions in(to) standard
134 // namespace before numbers::NumberTraits is defined
135 #ifdef DEAL_II_WITH_ADOLC
137 
138 # include <adolc/adouble.h> // Taped double
139 #endif
140 // Ideally we'd like to #include <deal.II/differentiation/ad/sacado_math.h>
141 // but header indirectly references numbers.h. We therefore simply
142 // import the whole Sacado header at this point to get the math
143 // functions imported into the standard namespace.
144 #ifdef DEAL_II_TRILINOS_WITH_SACADO
146 # include <Sacado.hpp>
148 #endif
149 
150 namespace std
151 {
152  template <typename Number, std::size_t width>
153  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
154  sqrt(const ::VectorizedArray<Number, width> &);
155  template <typename Number, std::size_t width>
156  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
157  abs(const ::VectorizedArray<Number, width> &);
158  template <typename Number, std::size_t width>
159  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
160  max(const ::VectorizedArray<Number, width> &,
161  const ::VectorizedArray<Number, width> &);
162  template <typename Number, std::size_t width>
163  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
164  min(const ::VectorizedArray<Number, width> &,
165  const ::VectorizedArray<Number, width> &);
166  template <typename Number, size_t width>
168  pow(const ::VectorizedArray<Number, width> &, const Number p);
169  template <typename Number, size_t width>
171  sin(const ::VectorizedArray<Number, width> &);
172  template <typename Number, size_t width>
174  cos(const ::VectorizedArray<Number, width> &);
175  template <typename Number, size_t width>
177  tan(const ::VectorizedArray<Number, width> &);
178  template <typename Number, size_t width>
180  exp(const ::VectorizedArray<Number, width> &);
181  template <typename Number, size_t width>
183  log(const ::VectorizedArray<Number, width> &);
184 } // namespace std
185 
187 
203 namespace numbers
204 {
208  static constexpr double E = 2.7182818284590452354;
209 
213  static constexpr double LOG2E = 1.4426950408889634074;
214 
218  static constexpr double LOG10E = 0.43429448190325182765;
219 
223  static constexpr double LN2 = 0.69314718055994530942;
224 
228  static constexpr double LN10 = 2.30258509299404568402;
229 
233  static constexpr double PI = 3.14159265358979323846;
234 
238  static constexpr double PI_2 = 1.57079632679489661923;
239 
243  static constexpr double PI_4 = 0.78539816339744830962;
244 
248  static constexpr double SQRT2 = 1.41421356237309504880;
249 
253  static constexpr double SQRT1_2 = 0.70710678118654752440;
254 
260  template <typename Number, typename = void>
261  struct is_cuda_compatible : std::true_type
262  {};
263 
267  template <typename Number>
268  struct is_cuda_compatible<std::complex<Number>, void> : std::false_type
269  {};
270 
280  bool
281  is_finite(const double x);
282 
287  bool
288  is_finite(const std::complex<double> &x);
289 
294  bool
295  is_finite(const std::complex<float> &x);
296 
305  bool
306  is_finite(const std::complex<long double> &x);
307 
318  template <typename Number1, typename Number2>
319  constexpr bool
320  values_are_equal(const Number1 &value_1, const Number2 &value_2);
321 
332  template <typename Number1, typename Number2>
333  bool
334  values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
335 
343  template <typename Number>
344  constexpr bool
345  value_is_zero(const Number &value);
346 
357  template <typename Number1, typename Number2>
358  bool
359  value_is_less_than(const Number1 &value_1, const Number2 &value_2);
360 
371  template <typename Number1, typename Number2>
372  bool
373  value_is_less_than_or_equal_to(const Number1 &value_1,
374  const Number2 &value_2);
375 
376 
377 
388  template <typename Number1, typename Number2>
389  bool
390  value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
391 
402  template <typename Number1, typename Number2>
403  bool
404  value_is_greater_than_or_equal_to(const Number1 &value_1,
405  const Number2 &value_2);
406 
415  template <typename number>
417  {
423  static constexpr bool is_complex = false;
424 
431  using real_type = number;
432 
437 
445  static constexpr DEAL_II_CUDA_HOST_DEV const number &
446  conjugate(const number &x);
447 
456  template <typename Dummy = number>
457  static constexpr DEAL_II_CUDA_HOST_DEV
458  typename std::enable_if<std::is_same<Dummy, number>::value &&
460  real_type>::type
461  abs_square(const number &x);
462 
463  template <typename Dummy = number>
464  static constexpr
465  typename std::enable_if<std::is_same<Dummy, number>::value &&
467  real_type>::type
468  abs_square(const number &x);
469 
473  static real_type
474  abs(const number &x);
475  };
476 
477 
482  template <typename number>
483  struct NumberTraits<std::complex<number>>
484  {
490  static constexpr bool is_complex = true;
491 
498  using real_type = number;
499 
503  using double_type = std::complex<double>;
504 
508  static constexpr std::complex<number>
509  conjugate(const std::complex<number> &x);
510 
517  static constexpr real_type
518  abs_square(const std::complex<number> &x);
519 
520 
524  static real_type
525  abs(const std::complex<number> &x);
526  };
527 
528  // --------------- inline and template functions ---------------- //
529 
530  inline bool
531  is_nan(const double x)
532  {
533  return std::isnan(x);
534  }
535 
536 
537 
538  inline bool
539  is_finite(const double x)
540  {
541  return std::isfinite(x);
542  }
543 
544 
545 
546  inline bool
547  is_finite(const std::complex<double> &x)
548  {
549  // Check complex numbers for infinity
550  // by testing real and imaginary part
551  return (is_finite(x.real()) && is_finite(x.imag()));
552  }
553 
554 
555 
556  inline bool
557  is_finite(const std::complex<float> &x)
558  {
559  // Check complex numbers for infinity
560  // by testing real and imaginary part
561  return (is_finite(x.real()) && is_finite(x.imag()));
562  }
563 
564 
565 
566  inline bool
567  is_finite(const std::complex<long double> &x)
568  {
569  // Same for std::complex<long double>
570  return (is_finite(x.real()) && is_finite(x.imag()));
571  }
572 
573 
574  template <typename number>
575  constexpr DEAL_II_CUDA_HOST_DEV const number &
577  {
578  return x;
579  }
580 
581 
582 
583  template <typename number>
584  template <typename Dummy>
585  constexpr DEAL_II_CUDA_HOST_DEV
586  typename std::enable_if<std::is_same<Dummy, number>::value &&
588  typename NumberTraits<number>::real_type>::type
590  {
591  return x * x;
592  }
593 
594 
595 
596  template <typename number>
597  template <typename Dummy>
598  constexpr
599  typename std::enable_if<std::is_same<Dummy, number>::value &&
601  typename NumberTraits<number>::real_type>::type
602  NumberTraits<number>::abs_square(const number &x)
603  {
604  return x * x;
605  }
606 
607 
608 
609  template <typename number>
611  NumberTraits<number>::abs(const number &x)
612  {
613  return std::abs(x);
614  }
615 
616 
617 
618  template <typename number>
619  constexpr std::complex<number>
620  NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
621  {
622  return std::conj(x);
623  }
624 
625 
626 
627  template <typename number>
628  typename NumberTraits<std::complex<number>>::real_type
629  NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
630  {
631  return std::abs(x);
632  }
633 
634 
635 
636  template <typename number>
637  constexpr typename NumberTraits<std::complex<number>>::real_type
638  NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
639  {
640  return std::norm(x);
641  }
642 
643 } // namespace numbers
644 
645 
646 // Forward declarations
648 {
649  namespace AD
650  {
651  namespace internal
652  {
653  // Defined in differentiation/ad/ad_number_traits.h
654  template <typename T>
655  struct NumberType;
656  } // namespace internal
657 
658  // Defined in differentiation/ad/ad_number_traits.h
659  template <typename NumberType>
660  struct is_ad_number;
661  } // namespace AD
662 } // namespace Differentiation
663 
664 
665 namespace internal
666 {
671  template <typename From, typename To>
673  {
674  // Source: https://stackoverflow.com/a/16944130
675  private:
676  template <typename T>
677  static void f(T);
678 
679  template <typename F, typename T>
680  static constexpr auto
681  test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
682  {
683  return true;
684  }
685 
686  template <typename F, typename T>
687  static constexpr auto
688  test(...) -> bool
689  {
690  return false;
691  }
692 
693  public:
694  static bool const value = test<From, To>(0);
695  };
696 
697  /*
698  * The structs below are needed to convert between some special number types.
699  * Also see tensor.h for another specialization.
700  */
701  template <typename T>
702  struct NumberType
703  {
704  static constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV const T &
705  value(const T &t)
706  {
707  return t;
708  }
709 
710  // Below are generic functions that allows an overload for any
711  // type U that is transformable to type T. This is particularly
712  // useful when needing to cast exotic number types
713  // (e.g. auto-differentiable or symbolic numbers) to a floating
714  // point one, such as might happen when converting between tensor
715  // types.
716 
717  // Type T is constructible from F.
718  template <typename F>
720  value(const F &f,
721  typename std::enable_if<
722  !std::is_same<typename std::decay<T>::type,
723  typename std::decay<F>::type>::value &&
724  std::is_constructible<T, F>::value>::type * = nullptr)
725  {
726  return T(f);
727  }
728 
729  // Type T is explicitly convertible (but not constructible) from F.
730  template <typename F>
731  static constexpr DEAL_II_ALWAYS_INLINE T
732  value(const F &f,
733  typename std::enable_if<
734  !std::is_same<typename std::decay<T>::type,
735  typename std::decay<F>::type>::value &&
736  !std::is_constructible<T, F>::value &&
738  {
739  return static_cast<T>(f);
740  }
741 
742  // Sacado doesn't provide any conversion operators, so we have
743  // to extract the value and perform further conversions from there.
744  // To be safe, we extend this to other possible AD numbers that
745  // might fall into the same category.
746  template <typename F>
747  static T
748  value(const F &f,
749  typename std::enable_if<
750  !std::is_same<typename std::decay<T>::type,
751  typename std::decay<F>::type>::value &&
752  !std::is_constructible<T, F>::value &&
755  {
757  }
758  };
759 
760  template <typename T>
761  struct NumberType<std::complex<T>>
762  {
763  static constexpr const std::complex<T> &
764  value(const std::complex<T> &t)
765  {
766  return t;
767  }
768 
769  static constexpr std::complex<T>
770  value(const T &t)
771  {
772  return std::complex<T>(t);
773  }
774 
775  // Facilitate cast from complex<double> to complex<float>
776  template <typename U>
777  static constexpr std::complex<T>
778  value(const std::complex<U> &t)
779  {
780  return std::complex<T>(NumberType<T>::value(t.real()),
781  NumberType<T>::value(t.imag()));
782  }
783  };
784 
785 #ifdef DEAL_II_COMPILER_CUDA_AWARE
786  template <>
787  struct NumberType<cuComplex>
788  {
789  static cuComplex
790  value(const float t)
791  {
792  return make_cuComplex(t, 0.f);
793  }
794  };
795 
796  template <>
797  struct NumberType<cuDoubleComplex>
798  {
799  static cuDoubleComplex
800  value(const double t)
801  {
802  return make_cuDoubleComplex(t, 0.);
803  }
804  };
805 #endif
806 } // namespace internal
807 
808 namespace numbers
809 {
810 #ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
811 
822  // Defined in differentiation/ad/adolc_number_types.cc
823  bool
824  values_are_equal(const adouble &value_1, const adouble &value_2);
825 
826 
837  template <typename Number>
838  bool
839  values_are_equal(const adouble &value_1, const Number &value_2)
840  {
841  // Use the specialized definition for two ADOL-C taped types
842  return values_are_equal(value_1,
844  }
845 
846 
857  template <typename Number>
858  bool
859  values_are_equal(const Number &value_1, const adouble &value_2)
860  {
861  // Use the above definition
862  return values_are_equal(value_2, value_1);
863  }
864 
876  // Defined in differentiation/ad/adolc_number_types.cc
877  bool
878  value_is_less_than(const adouble &value_1, const adouble &value_2);
879 
880 
892  template <typename Number>
893  bool
894  value_is_less_than(const adouble &value_1, const Number &value_2)
895  {
896  // Use the specialized definition for two ADOL-C taped types
897  return value_is_less_than(value_1,
899  }
900 
901 
913  template <typename Number>
914  bool
915  value_is_less_than(const Number &value_1, const adouble &value_2)
916  {
917  // Use the specialized definition for two ADOL-C taped types
919  value_2);
920  }
921 
922 #endif
923 
924 
925  template <typename Number1, typename Number2>
926  constexpr bool
927  values_are_equal(const Number1 &value_1, const Number2 &value_2)
928  {
929  return (value_1 == internal::NumberType<Number1>::value(value_2));
930  }
931 
932 
933  template <typename Number1, typename Number2>
934  inline bool
935  values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
936  {
937  return !(values_are_equal(value_1, value_2));
938  }
939 
940 
941  template <typename Number>
942  constexpr bool
943  value_is_zero(const Number &value)
944  {
945  return values_are_equal(value, 0.0);
946  }
947 
948 
949  template <typename Number1, typename Number2>
950  inline bool
951  value_is_less_than(const Number1 &value_1, const Number2 &value_2)
952  {
953  return (value_1 < internal::NumberType<Number1>::value(value_2));
954  }
955 
956 
957  template <typename Number1, typename Number2>
958  inline bool
959  value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
960  {
961  return (value_is_less_than(value_1, value_2) ||
962  values_are_equal(value_1, value_2));
963  }
964 
965 
966  template <typename Number1, typename Number2>
967  bool
968  value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
969  {
970  return !(value_is_less_than_or_equal_to(value_1, value_2));
971  }
972 
973 
974  template <typename Number1, typename Number2>
975  inline bool
976  value_is_greater_than_or_equal_to(const Number1 &value_1,
977  const Number2 &value_2)
978  {
979  return !(value_is_less_than(value_1, value_2));
980  }
981 } // namespace numbers
982 
984 
985 #endif
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:102
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:416
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:430
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:417
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:469
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)
static constexpr double LOG10E
Definition: numbers.h:218
static constexpr double PI_2
Definition: numbers.h:238
bool value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:959
static constexpr double E
Definition: numbers.h:208
static constexpr double PI
Definition: numbers.h:233
bool value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:968
static constexpr double SQRT2
Definition: numbers.h:248
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:943
bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:935
static constexpr double SQRT1_2
Definition: numbers.h:253
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:927
static constexpr double PI_4
Definition: numbers.h:243
static constexpr double LN10
Definition: numbers.h:228
static constexpr double LN2
Definition: numbers.h:223
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:951
bool is_finite(const double x)
Definition: numbers.h:539
static constexpr double LOG2E
Definition: numbers.h:213
bool value_is_greater_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:976
bool is_nan(const double x)
Definition: numbers.h:531
#define DEAL_II_CUDA_HOST_DEV
Definition: numbers.h:34
static cuComplex value(const float t)
Definition: numbers.h:790
static cuDoubleComplex value(const double t)
Definition: numbers.h:800
static constexpr std::complex< T > value(const std::complex< U > &t)
Definition: numbers.h:778
static constexpr std::complex< T > value(const T &t)
Definition: numbers.h:770
static constexpr const std::complex< T > & value(const std::complex< T > &t)
Definition: numbers.h:764
static constexpr T value(const F &f, typename std::enable_if< !std::is_same< typename std::decay< T >::type, typename std::decay< F >::type >::value &&std::is_constructible< T, F >::value >::type *=nullptr)
Definition: numbers.h:720
static constexpr const T & value(const T &t)
Definition: numbers.h:705
static T value(const F &f, typename std::enable_if< !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 >::type *=nullptr)
Definition: numbers.h:748
static constexpr T value(const F &f, typename std::enable_if< !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 >::type *=nullptr)
Definition: numbers.h:732
constexpr static unsigned int max_width
Definition: numbers.h:65
static constexpr auto test(...) -> bool
Definition: numbers.h:688
static constexpr auto test(int) -> decltype(f(static_cast< T >(std::declval< F >())), true)
Definition: numbers.h:681
static constexpr const number & conjugate(const number &x)
Definition: numbers.h:576
static constexpr bool is_complex
Definition: numbers.h:423
static constexpr std::enable_if< std::is_same< Dummy, number >::value &&is_cuda_compatible< Dummy >::value, real_type >::type abs_square(const number &x)
Definition: numbers.h:589
static real_type abs(const number &x)
Definition: numbers.h:611
static constexpr std::enable_if< std::is_same< Dummy, number >::value &&!is_cuda_compatible< Dummy >::value, real_type >::type abs_square(const number &x)