Reference documentation for deal.II version Git ac8d010384 2020-11-27 19:49:05 +0100
\(\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 {
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
145 # include <Sacado.hpp>
146 #endif
147 
148 namespace std
149 {
150  template <typename Number, std::size_t width>
151  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
152  sqrt(const ::VectorizedArray<Number, width> &);
153  template <typename Number, std::size_t width>
154  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
155  abs(const ::VectorizedArray<Number, width> &);
156  template <typename Number, std::size_t width>
157  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
158  max(const ::VectorizedArray<Number, width> &,
159  const ::VectorizedArray<Number, width> &);
160  template <typename Number, std::size_t width>
161  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
162  min(const ::VectorizedArray<Number, width> &,
163  const ::VectorizedArray<Number, width> &);
164  template <typename Number, size_t width>
166  pow(const ::VectorizedArray<Number, width> &, const Number p);
167  template <typename Number, size_t width>
169  sin(const ::VectorizedArray<Number, width> &);
170  template <typename Number, size_t width>
172  cos(const ::VectorizedArray<Number, width> &);
173  template <typename Number, size_t width>
175  tan(const ::VectorizedArray<Number, width> &);
176  template <typename Number, size_t width>
178  exp(const ::VectorizedArray<Number, width> &);
179  template <typename Number, size_t width>
181  log(const ::VectorizedArray<Number, width> &);
182 } // namespace std
183 
185 
201 namespace numbers
202 {
206  static constexpr double E = 2.7182818284590452354;
207 
211  static constexpr double LOG2E = 1.4426950408889634074;
212 
216  static constexpr double LOG10E = 0.43429448190325182765;
217 
221  static constexpr double LN2 = 0.69314718055994530942;
222 
226  static constexpr double LN10 = 2.30258509299404568402;
227 
231  static constexpr double PI = 3.14159265358979323846;
232 
236  static constexpr double PI_2 = 1.57079632679489661923;
237 
241  static constexpr double PI_4 = 0.78539816339744830962;
242 
246  static constexpr double SQRT2 = 1.41421356237309504880;
247 
251  static constexpr double SQRT1_2 = 0.70710678118654752440;
252 
258  template <typename Number, typename = void>
259  struct is_cuda_compatible : std::true_type
260  {};
261 
265  template <typename Number>
266  struct is_cuda_compatible<std::complex<Number>, void> : std::false_type
267  {};
268 
278  bool
279  is_finite(const double x);
280 
285  bool
286  is_finite(const std::complex<double> &x);
287 
292  bool
293  is_finite(const std::complex<float> &x);
294 
303  bool
304  is_finite(const std::complex<long double> &x);
305 
316  template <typename Number1, typename Number2>
317  constexpr bool
318  values_are_equal(const Number1 &value_1, const Number2 &value_2);
319 
330  template <typename Number1, typename Number2>
331  bool
332  values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
333 
341  template <typename Number>
342  constexpr bool
343  value_is_zero(const Number &value);
344 
355  template <typename Number1, typename Number2>
356  bool
357  value_is_less_than(const Number1 &value_1, const Number2 &value_2);
358 
369  template <typename Number1, typename Number2>
370  bool
371  value_is_less_than_or_equal_to(const Number1 &value_1,
372  const Number2 &value_2);
373 
374 
375 
386  template <typename Number1, typename Number2>
387  bool
388  value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
389 
400  template <typename Number1, typename Number2>
401  bool
402  value_is_greater_than_or_equal_to(const Number1 &value_1,
403  const Number2 &value_2);
404 
413  template <typename number>
415  {
421  static constexpr bool is_complex = false;
422 
429  using real_type = number;
430 
438  static constexpr DEAL_II_CUDA_HOST_DEV const number &
439  conjugate(const number &x);
440 
449  template <typename Dummy = number>
450  static constexpr DEAL_II_CUDA_HOST_DEV
451  typename std::enable_if<std::is_same<Dummy, number>::value &&
453  real_type>::type
454  abs_square(const number &x);
455 
456  template <typename Dummy = number>
457  static constexpr
458  typename std::enable_if<std::is_same<Dummy, number>::value &&
459  !is_cuda_compatible<Dummy>::value,
460  real_type>::type
461  abs_square(const number &x);
462 
466  static real_type
467  abs(const number &x);
468  };
469 
470 
475  template <typename number>
476  struct NumberTraits<std::complex<number>>
477  {
483  static constexpr bool is_complex = true;
484 
491  using real_type = number;
492 
496  static constexpr std::complex<number>
497  conjugate(const std::complex<number> &x);
498 
505  static constexpr real_type
506  abs_square(const std::complex<number> &x);
507 
508 
512  static real_type
513  abs(const std::complex<number> &x);
514  };
515 
516  // --------------- inline and template functions ---------------- //
517 
518  inline bool
519  is_nan(const double x)
520  {
521  return std::isnan(x);
522  }
523 
524 
525 
526  inline bool
527  is_finite(const double x)
528  {
529  return std::isfinite(x);
530  }
531 
532 
533 
534  inline bool
535  is_finite(const std::complex<double> &x)
536  {
537  // Check complex numbers for infinity
538  // by testing real and imaginary part
539  return (is_finite(x.real()) && is_finite(x.imag()));
540  }
541 
542 
543 
544  inline bool
545  is_finite(const std::complex<float> &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<long double> &x)
556  {
557  // Same for std::complex<long double>
558  return (is_finite(x.real()) && is_finite(x.imag()));
559  }
560 
561 
562  template <typename number>
563  constexpr DEAL_II_CUDA_HOST_DEV const number &
565  {
566  return x;
567  }
568 
569 
570 
571  template <typename number>
572  template <typename Dummy>
573  constexpr DEAL_II_CUDA_HOST_DEV
574  typename std::enable_if<std::is_same<Dummy, number>::value &&
576  typename NumberTraits<number>::real_type>::type
578  {
579  return x * x;
580  }
581 
582 
583 
584  template <typename number>
585  template <typename Dummy>
586  constexpr
587  typename std::enable_if<std::is_same<Dummy, number>::value &&
588  !is_cuda_compatible<Dummy>::value,
589  typename NumberTraits<number>::real_type>::type
590  NumberTraits<number>::abs_square(const number &x)
591  {
592  return x * x;
593  }
594 
595 
596 
597  template <typename number>
599  NumberTraits<number>::abs(const number &x)
600  {
601  return std::abs(x);
602  }
603 
604 
605 
606  template <typename number>
607  constexpr std::complex<number>
608  NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
609  {
610  return std::conj(x);
611  }
612 
613 
614 
615  template <typename number>
617  NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
618  {
619  return std::abs(x);
620  }
621 
622 
623 
624  template <typename number>
626  NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
627  {
628  return std::norm(x);
629  }
630 
631 } // namespace numbers
632 
633 
634 // Forward declarations
636 {
637  namespace AD
638  {
639  namespace internal
640  {
641  // Defined in differentiation/ad/ad_number_traits.h
642  template <typename T>
643  struct NumberType;
644  } // namespace internal
645 
646  // Defined in differentiation/ad/ad_number_traits.h
647  template <typename NumberType>
648  struct is_ad_number;
649  } // namespace AD
650 } // namespace Differentiation
651 
652 
653 namespace internal
654 {
659  template <typename From, typename To>
661  {
662  // Source: https://stackoverflow.com/a/16944130
663  private:
664  template <typename T>
665  static void f(T);
666 
667  template <typename F, typename T>
668  static constexpr auto
669  test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
670  {
671  return true;
672  }
673 
674  template <typename F, typename T>
675  static constexpr auto
676  test(...) -> bool
677  {
678  return false;
679  }
680 
681  public:
682  static bool const value = test<From, To>(0);
683  };
684 
685  /*
686  * The structs below are needed to convert between some special number types.
687  * Also see tensor.h for another specialization.
688  */
689  template <typename T>
690  struct NumberType
691  {
692  static constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV const T &
693  value(const T &t)
694  {
695  return t;
696  }
697 
698  // Below are generic functions that allows an overload for any
699  // type U that is transformable to type T. This is particularly
700  // useful when needing to cast exotic number types
701  // (e.g. auto-differentiable or symbolic numbers) to a floating
702  // point one, such as might happen when converting between tensor
703  // types.
704 
705  // Type T is constructible from F.
706  template <typename F>
708  value(const F &f,
709  typename std::enable_if<
710  !std::is_same<typename std::decay<T>::type,
711  typename std::decay<F>::type>::value &&
712  std::is_constructible<T, F>::value>::type * = nullptr)
713  {
714  return T(f);
715  }
716 
717  // Type T is explicitly convertible (but not constructible) from F.
718  template <typename F>
719  static constexpr DEAL_II_ALWAYS_INLINE T
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 &&
726  {
727  return static_cast<T>(f);
728  }
729 
730  // Sacado doesn't provide any conversion operators, so we have
731  // to extract the value and perform further conversions from there.
732  // To be safe, we extend this to other possible AD numbers that
733  // might fall into the same category.
734  template <typename F>
735  static T
736  value(const F &f,
737  typename std::enable_if<
738  !std::is_same<typename std::decay<T>::type,
739  typename std::decay<F>::type>::value &&
740  !std::is_constructible<T, F>::value &&
743  {
745  }
746  };
747 
748  template <typename T>
749  struct NumberType<std::complex<T>>
750  {
751  static constexpr const std::complex<T> &
752  value(const std::complex<T> &t)
753  {
754  return t;
755  }
756 
757  static constexpr std::complex<T>
758  value(const T &t)
759  {
760  return std::complex<T>(t);
761  }
762 
763  // Facilitate cast from complex<double> to complex<float>
764  template <typename U>
765  static constexpr std::complex<T>
766  value(const std::complex<U> &t)
767  {
768  return std::complex<T>(NumberType<T>::value(t.real()),
769  NumberType<T>::value(t.imag()));
770  }
771  };
772 
773 #ifdef DEAL_II_COMPILER_CUDA_AWARE
774  template <>
775  struct NumberType<cuComplex>
776  {
777  static cuComplex
778  value(const float t)
779  {
780  return make_cuComplex(t, 0.f);
781  }
782  };
783 
784  template <>
785  struct NumberType<cuDoubleComplex>
786  {
787  static cuDoubleComplex
788  value(const double t)
789  {
790  return make_cuDoubleComplex(t, 0.);
791  }
792  };
793 #endif
794 } // namespace internal
795 
796 namespace numbers
797 {
798 #ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
799 
810  // Defined in differentiation/ad/adolc_number_types.cc
811  bool
812  values_are_equal(const adouble &value_1, const adouble &value_2);
813 
814 
825  template <typename Number>
826  bool
827  values_are_equal(const adouble &value_1, const Number &value_2)
828  {
829  // Use the specialized definition for two ADOL-C taped types
830  return values_are_equal(value_1,
832  }
833 
834 
845  template <typename Number>
846  bool
847  values_are_equal(const Number &value_1, const adouble &value_2)
848  {
849  // Use the above definition
850  return values_are_equal(value_2, value_1);
851  }
852 
864  // Defined in differentiation/ad/adolc_number_types.cc
865  bool
866  value_is_less_than(const adouble &value_1, const adouble &value_2);
867 
868 
880  template <typename Number>
881  bool
882  value_is_less_than(const adouble &value_1, const Number &value_2)
883  {
884  // Use the specialized definition for two ADOL-C taped types
885  return value_is_less_than(value_1,
887  }
888 
889 
901  template <typename Number>
902  bool
903  value_is_less_than(const Number &value_1, const adouble &value_2)
904  {
905  // Use the specialized definition for two ADOL-C taped types
907  value_2);
908  }
909 
910 #endif
911 
912 
913  template <typename Number1, typename Number2>
914  constexpr bool
915  values_are_equal(const Number1 &value_1, const Number2 &value_2)
916  {
917  return (value_1 == internal::NumberType<Number1>::value(value_2));
918  }
919 
920 
921  template <typename Number1, typename Number2>
922  inline bool
923  values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
924  {
925  return !(values_are_equal(value_1, value_2));
926  }
927 
928 
929  template <typename Number>
930  constexpr bool
931  value_is_zero(const Number &value)
932  {
933  return values_are_equal(value, 0.0);
934  }
935 
936 
937  template <typename Number1, typename Number2>
938  inline bool
939  value_is_less_than(const Number1 &value_1, const Number2 &value_2)
940  {
941  return (value_1 < internal::NumberType<Number1>::value(value_2));
942  }
943 
944 
945  template <typename Number1, typename Number2>
946  inline bool
947  value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
948  {
949  return (value_is_less_than(value_1, value_2) ||
950  values_are_equal(value_1, value_2));
951  }
952 
953 
954  template <typename Number1, typename Number2>
955  bool
956  value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
957  {
958  return !(value_is_less_than_or_equal_to(value_1, value_2));
959  }
960 
961 
962  template <typename Number1, typename Number2>
963  inline bool
964  value_is_greater_than_or_equal_to(const Number1 &value_1,
965  const Number2 &value_2)
966  {
967  return !(value_is_less_than(value_1, value_2));
968  }
969 } // namespace numbers
970 
972 
973 #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:736
static constexpr double LOG2E
Definition: numbers.h:211
static constexpr unsigned int max_width
Definition: numbers.h:65
bool value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:947
static constexpr const T & value(const T &t)
Definition: numbers.h:693
static constexpr double LN10
Definition: numbers.h:226
static constexpr std::complex< T > value(const T &t)
Definition: numbers.h:758
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:720
bool values_are_equal(const Number &value_1, const adouble &value_2)
Definition: numbers.h:847
static cuDoubleComplex value(const double t)
Definition: numbers.h:788
static constexpr std::complex< T > value(const std::complex< U > &t)
Definition: numbers.h:766
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
STL namespace.
static constexpr double LOG10E
Definition: numbers.h:216
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:903
bool is_finite(const std::complex< long double > &x)
Definition: numbers.h:555
static constexpr double SQRT1_2
Definition: numbers.h:251
static constexpr double E
Definition: numbers.h:206
static const char T
static constexpr double SQRT2
Definition: numbers.h:246
static constexpr auto test(int) -> decltype(f(static_cast< T >(std::declval< F >())), true)
Definition: numbers.h:669
bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:923
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:94
static constexpr auto test(...) -> bool
Definition: numbers.h:676
static cuComplex value(const float t)
Definition: numbers.h:778
static constexpr double PI_2
Definition: numbers.h:236
bool is_nan(const double x)
Definition: numbers.h:519
bool value_is_greater_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:964
static constexpr double LN2
Definition: numbers.h:221
bool value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:956
static constexpr double PI
Definition: numbers.h:231
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
static constexpr const std::complex< T > & value(const std::complex< T > &t)
Definition: numbers.h:752
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:931
#define DEAL_II_CUDA_HOST_DEV
Definition: numbers.h:34
static constexpr double PI_4
Definition: numbers.h:241
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:708