Reference documentation for deal.II version Git 423cd11810 2020-05-28 01:07:25 +0200
\(\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 - 2020 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 {
61  template <typename Number>
63  {
67  constexpr static unsigned int max_width = 1;
68  };
69 
78  template <>
80  {
84  constexpr static unsigned int max_width =
85 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512
86  8;
87 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256
88  4;
89 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128
90  2;
91 #else
92  1;
93 #endif
94  };
95 
104  template <>
106  {
110  constexpr static unsigned int max_width =
111 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ALTIVEC__)
112  4;
113 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && defined(__AVX512F__)
114  16;
115 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)
116  8;
117 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__)
118  4;
119 #else
120  1;
121 #endif
122  };
123 
124 
125 } // namespace internal
126 
127 // forward declarations to support abs or sqrt operations on VectorizedArray
128 #ifndef DOXYGEN
129 template <typename Number,
130  std::size_t width =
132 class VectorizedArray;
133 template <typename T>
134 struct EnableIfScalar;
135 #endif
136 
138 
139 // Declare / Import auto-differentiable math functions in(to) standard
140 // namespace before numbers::NumberTraits is defined
141 #ifdef DEAL_II_WITH_ADOLC
143 
144 # include <adolc/adouble.h> // Taped double
145 #endif
146 // Ideally we'd like to #include <deal.II/differentiation/ad/sacado_math.h>
147 // but header indirectly references numbers.h. We therefore simply
148 // import the whole Sacado header at this point to get the math
149 // functions imported into the standard namespace.
150 #ifdef DEAL_II_TRILINOS_WITH_SACADO
151 # include <Sacado.hpp>
152 #endif
153 
154 namespace std
155 {
156  template <typename Number, std::size_t width>
157  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
158  sqrt(const ::VectorizedArray<Number, width> &);
159  template <typename Number, std::size_t width>
160  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
161  abs(const ::VectorizedArray<Number, width> &);
162  template <typename Number, std::size_t width>
163  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
164  max(const ::VectorizedArray<Number, width> &,
165  const ::VectorizedArray<Number, width> &);
166  template <typename Number, std::size_t width>
167  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
168  min(const ::VectorizedArray<Number, width> &,
169  const ::VectorizedArray<Number, width> &);
170  template <typename Number, size_t width>
172  pow(const ::VectorizedArray<Number, width> &, const Number p);
173  template <typename Number, size_t width>
175  sin(const ::VectorizedArray<Number, width> &);
176  template <typename Number, size_t width>
178  cos(const ::VectorizedArray<Number, width> &);
179  template <typename Number, size_t width>
181  tan(const ::VectorizedArray<Number, width> &);
182  template <typename Number, size_t width>
184  exp(const ::VectorizedArray<Number, width> &);
185  template <typename Number, size_t width>
187  log(const ::VectorizedArray<Number, width> &);
188 } // namespace std
189 
191 
207 namespace numbers
208 {
212  static constexpr double E = 2.7182818284590452354;
213 
217  static constexpr double LOG2E = 1.4426950408889634074;
218 
222  static constexpr double LOG10E = 0.43429448190325182765;
223 
227  static constexpr double LN2 = 0.69314718055994530942;
228 
232  static constexpr double LN10 = 2.30258509299404568402;
233 
237  static constexpr double PI = 3.14159265358979323846;
238 
242  static constexpr double PI_2 = 1.57079632679489661923;
243 
247  static constexpr double PI_4 = 0.78539816339744830962;
248 
252  static constexpr double SQRT2 = 1.41421356237309504880;
253 
257  static constexpr double SQRT1_2 = 0.70710678118654752440;
258 
264  template <typename Number, typename = void>
265  struct is_cuda_compatible : std::true_type
266  {};
267 
271  template <typename Number>
272  struct is_cuda_compatible<std::complex<Number>, void> : std::false_type
273  {};
274 
284  bool
285  is_finite(const double x);
286 
291  bool
292  is_finite(const std::complex<double> &x);
293 
298  bool
299  is_finite(const std::complex<float> &x);
300 
309  bool
310  is_finite(const std::complex<long double> &x);
311 
322  template <typename Number1, typename Number2>
323  constexpr bool
324  values_are_equal(const Number1 &value_1, const Number2 &value_2);
325 
336  template <typename Number1, typename Number2>
337  bool
338  values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
339 
347  template <typename Number>
348  constexpr bool
349  value_is_zero(const Number &value);
350 
361  template <typename Number1, typename Number2>
362  bool
363  value_is_less_than(const Number1 &value_1, const Number2 &value_2);
364 
375  template <typename Number1, typename Number2>
376  bool
377  value_is_less_than_or_equal_to(const Number1 &value_1,
378  const Number2 &value_2);
379 
380 
381 
392  template <typename Number1, typename Number2>
393  bool
394  value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
395 
406  template <typename Number1, typename Number2>
407  bool
408  value_is_greater_than_or_equal_to(const Number1 &value_1,
409  const Number2 &value_2);
410 
421  template <typename number>
423  {
429  static constexpr bool is_complex = false;
430 
437  using real_type = number;
438 
446  static constexpr DEAL_II_CUDA_HOST_DEV const number &
447  conjugate(const number &x);
448 
457  template <typename Dummy = number>
458  static constexpr DEAL_II_CUDA_HOST_DEV
461  real_type>::type
462  abs_square(const number &x);
463 
464  template <typename Dummy = number>
465  static constexpr
467  !is_cuda_compatible<Dummy>::value,
468  real_type>::type
469  abs_square(const number &x);
470 
474  static real_type
475  abs(const number &x);
476  };
477 
478 
485  template <typename number>
486  struct NumberTraits<std::complex<number>>
487  {
493  static constexpr bool is_complex = true;
494 
501  using real_type = number;
502 
506  static constexpr std::complex<number>
507  conjugate(const std::complex<number> &x);
508 
515  static constexpr real_type
516  abs_square(const std::complex<number> &x);
517 
518 
522  static real_type
523  abs(const std::complex<number> &x);
524  };
525 
526  // --------------- inline and template functions ---------------- //
527 
528  inline bool
529  is_nan(const double x)
530  {
531  return std::isnan(x);
532  }
533 
534 
535 
536  inline bool
537  is_finite(const double x)
538  {
539  return std::isfinite(x);
540  }
541 
542 
543 
544  inline bool
545  is_finite(const std::complex<double> &x)
546  {
547  // Check complex numbers for infinity
548  // by testing real and imaginary part
549  return (is_finite(x.real()) && is_finite(x.imag()));
550  }
551 
552 
553 
554  inline bool
555  is_finite(const std::complex<float> &x)
556  {
557  // Check complex numbers for infinity
558  // by testing real and imaginary part
559  return (is_finite(x.real()) && is_finite(x.imag()));
560  }
561 
562 
563 
564  inline bool
565  is_finite(const std::complex<long double> &x)
566  {
567  // Same for std::complex<long double>
568  return (is_finite(x.real()) && is_finite(x.imag()));
569  }
570 
571 
572  template <typename number>
573  constexpr DEAL_II_CUDA_HOST_DEV const number &
575  {
576  return x;
577  }
578 
579 
580 
581  template <typename number>
582  template <typename Dummy>
583  constexpr DEAL_II_CUDA_HOST_DEV
586  typename NumberTraits<number>::real_type>::type
588  {
589  return x * x;
590  }
591 
592 
593 
594  template <typename number>
595  template <typename Dummy>
596  constexpr
598  !is_cuda_compatible<Dummy>::value,
599  typename NumberTraits<number>::real_type>::type
600  NumberTraits<number>::abs_square(const number &x)
601  {
602  return x * x;
603  }
604 
605 
606 
607  template <typename number>
609  NumberTraits<number>::abs(const number &x)
610  {
611  return std::abs(x);
612  }
613 
614 
615 
616  template <typename number>
617  constexpr std::complex<number>
618  NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
619  {
620  return std::conj(x);
621  }
622 
623 
624 
625  template <typename number>
627  NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
628  {
629  return std::abs(x);
630  }
631 
632 
633 
634  template <typename number>
636  NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
637  {
638  return std::norm(x);
639  }
640 
641 } // namespace numbers
642 
643 
644 // Forward declarations
646 {
647  namespace AD
648  {
649  namespace internal
650  {
651  // Defined in differentiation/ad/ad_number_traits.h
652  template <typename T>
653  struct NumberType;
654  } // namespace internal
655 
656  // Defined in differentiation/ad/ad_number_traits.h
657  template <typename NumberType>
658  struct is_ad_number;
659  } // namespace AD
660 } // namespace Differentiation
661 
662 
663 namespace internal
664 {
669  template <typename From, typename To>
671  {
672  // Source: https://stackoverflow.com/a/16944130
673  private:
674  template <typename T>
675  static void f(T);
676 
677  template <typename F, typename T>
678  static constexpr auto
679  test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
680  {
681  return true;
682  }
683 
684  template <typename F, typename T>
685  static constexpr auto
686  test(...) -> bool
687  {
688  return false;
689  }
690 
691  public:
692  static bool const value = test<From, To>(0);
693  };
694 
695  /*
696  * The structs below are needed to convert between some special number types.
697  * Also see tensor.h for another specialization.
698  */
699  template <typename T>
700  struct NumberType
701  {
702  static constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV const T &
703  value(const T &t)
704  {
705  return t;
706  }
707 
708  // Below are generic functions that allows an overload for any
709  // type U that is transformable to type T. This is particularly
710  // useful when needing to cast exotic number types
711  // (e.g. auto-differentiable or symbolic numbers) to a floating
712  // point one, such as might happen when converting between tensor
713  // types.
714 
715  // Type T is constructible from F.
716  template <typename F>
718  value(const F &f,
719  typename std::enable_if<
720  !std::is_same<typename std::decay<T>::type,
721  typename std::decay<F>::type>::value &&
722  std::is_constructible<T, F>::value>::type * = nullptr)
723  {
724  return T(f);
725  }
726 
727  // Type T is explicitly convertible (but not constructible) from F.
728  template <typename F>
729  static constexpr DEAL_II_ALWAYS_INLINE T
730  value(const F &f,
731  typename std::enable_if<
732  !std::is_same<typename std::decay<T>::type,
733  typename std::decay<F>::type>::value &&
736  {
737  return static_cast<T>(f);
738  }
739 
740  // Sacado doesn't provide any conversion operators, so we have
741  // to extract the value and perform further conversions from there.
742  // To be safe, we extend this to other possible AD numbers that
743  // might fall into the same category.
744  template <typename F>
745  static T
746  value(const F &f,
747  typename std::enable_if<
748  !std::is_same<typename std::decay<T>::type,
749  typename std::decay<F>::type>::value &&
753  {
755  }
756  };
757 
758  template <typename T>
759  struct NumberType<std::complex<T>>
760  {
761  static constexpr const std::complex<T> &
762  value(const std::complex<T> &t)
763  {
764  return t;
765  }
766 
767  static constexpr std::complex<T>
768  value(const T &t)
769  {
770  return std::complex<T>(t);
771  }
772 
773  // Facilitate cast from complex<double> to complex<float>
774  template <typename U>
775  static constexpr std::complex<T>
776  value(const std::complex<U> &t)
777  {
778  return std::complex<T>(NumberType<T>::value(t.real()),
779  NumberType<T>::value(t.imag()));
780  }
781  };
782 
783 #ifdef DEAL_II_COMPILER_CUDA_AWARE
784  template <>
785  struct NumberType<cuComplex>
786  {
787  static cuComplex
788  value(const float t)
789  {
790  return make_cuComplex(t, 0.f);
791  }
792  };
793 
794  template <>
795  struct NumberType<cuDoubleComplex>
796  {
797  static cuDoubleComplex
798  value(const double t)
799  {
800  return make_cuDoubleComplex(t, 0.);
801  }
802  };
803 #endif
804 } // namespace internal
805 
806 namespace numbers
807 {
808 #ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
809 
820  // Defined in differentiation/ad/adolc_number_types.cc
821  bool
822  values_are_equal(const adouble &value_1, const adouble &value_2);
823 
824 
835  template <typename Number>
836  bool
837  values_are_equal(const adouble &value_1, const Number &value_2)
838  {
839  // Use the specialized definition for two ADOL-C taped types
840  return values_are_equal(value_1,
842  }
843 
844 
855  template <typename Number>
856  bool
857  values_are_equal(const Number &value_1, const adouble &value_2)
858  {
859  // Use the above definition
860  return values_are_equal(value_2, value_1);
861  }
862 
874  // Defined in differentiation/ad/adolc_number_types.cc
875  bool
876  value_is_less_than(const adouble &value_1, const adouble &value_2);
877 
878 
890  template <typename Number>
891  bool
892  value_is_less_than(const adouble &value_1, const Number &value_2)
893  {
894  // Use the specialized definition for two ADOL-C taped types
895  return value_is_less_than(value_1,
897  }
898 
899 
911  template <typename Number>
912  bool
913  value_is_less_than(const Number &value_1, const adouble &value_2)
914  {
915  // Use the specialized definition for two ADOL-C taped types
917  value_2);
918  }
919 
920 #endif
921 
922 
923  template <typename Number1, typename Number2>
924  constexpr bool
925  values_are_equal(const Number1 &value_1, const Number2 &value_2)
926  {
927  return (value_1 == internal::NumberType<Number1>::value(value_2));
928  }
929 
930 
931  template <typename Number1, typename Number2>
932  inline bool
933  values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
934  {
935  return !(values_are_equal(value_1, value_2));
936  }
937 
938 
939  template <typename Number>
940  constexpr bool
941  value_is_zero(const Number &value)
942  {
943  return values_are_equal(value, 0.0);
944  }
945 
946 
947  template <typename Number1, typename Number2>
948  inline bool
949  value_is_less_than(const Number1 &value_1, const Number2 &value_2)
950  {
951  return (value_1 < internal::NumberType<Number1>::value(value_2));
952  }
953 
954 
955  template <typename Number1, typename Number2>
956  inline bool
957  value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
958  {
959  return (value_is_less_than(value_1, value_2) ||
960  values_are_equal(value_1, value_2));
961  }
962 
963 
964  template <typename Number1, typename Number2>
965  bool
966  value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
967  {
968  return !(value_is_less_than_or_equal_to(value_1, value_2));
969  }
970 
971 
972  template <typename Number1, typename Number2>
973  inline bool
974  value_is_greater_than_or_equal_to(const Number1 &value_1,
975  const Number2 &value_2)
976  {
977  return !(value_is_less_than(value_1, value_2));
978  }
979 } // namespace numbers
980 
982 
983 #endif
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:746
static constexpr double LOG2E
Definition: numbers.h:217
static constexpr unsigned int max_width
Definition: numbers.h:67
bool value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:957
static constexpr const T & value(const T &t)
Definition: numbers.h:703
static constexpr double LN10
Definition: numbers.h:232
static constexpr std::complex< T > value(const T &t)
Definition: numbers.h:768
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:730
bool values_are_equal(const Number &value_1, const adouble &value_2)
Definition: numbers.h:857
static cuDoubleComplex value(const double t)
Definition: numbers.h:798
static constexpr std::complex< T > value(const std::complex< U > &t)
Definition: numbers.h:776
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
STL namespace.
static constexpr double LOG10E
Definition: numbers.h:222
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
bool value_is_less_than(const Number &value_1, const adouble &value_2)
Definition: numbers.h:913
bool is_finite(const std::complex< long double > &x)
Definition: numbers.h:565
static constexpr double SQRT1_2
Definition: numbers.h:257
static constexpr double E
Definition: numbers.h:212
static const char T
static constexpr double SQRT2
Definition: numbers.h:252
static constexpr auto test(int) -> decltype(f(static_cast< T >(std::declval< F >())), true)
Definition: numbers.h:679
bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:933
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:362
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:93
static constexpr auto test(...) -> bool
Definition: numbers.h:686
static cuComplex value(const float t)
Definition: numbers.h:788
static constexpr double PI_2
Definition: numbers.h:242
bool is_nan(const double x)
Definition: numbers.h:529
bool value_is_greater_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:974
static constexpr double LN2
Definition: numbers.h:227
bool value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:966
static constexpr double PI
Definition: numbers.h:237
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:361
static constexpr const std::complex< T > & value(const std::complex< T > &t)
Definition: numbers.h:762
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:941
#define DEAL_II_CUDA_HOST_DEV
Definition: numbers.h:34
static constexpr double PI_4
Definition: numbers.h:247
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:718
static const bool value