Reference documentation for deal.II version Git 2c472e7589 2019-11-13 14:15:54 +0100
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
numbers.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2006 - 2019 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 <cstdlib>
31 
32 #ifdef DEAL_II_COMPILER_CUDA_AWARE
33 # define DEAL_II_CUDA_HOST_DEV __host__ __device__
34 #else
35 # define DEAL_II_CUDA_HOST_DEV
36 #endif
37 
38 DEAL_II_NAMESPACE_OPEN
39 
40 namespace internal
41 {
60  template <typename Number>
62  {
66  constexpr static unsigned int max_width = 1;
67  };
68 
77  template <>
79  {
83  constexpr static unsigned int max_width =
84 #if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__ALTIVEC__)
85  2;
86 #elif DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 3 && defined(__AVX512F__)
87  8;
88 #elif DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 2 && defined(__AVX__)
89  4;
90 #elif DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__SSE2__)
91  2;
92 #else
93  1;
94 #endif
95  };
96 
105  template <>
107  {
111  constexpr static unsigned int max_width =
112 #if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__ALTIVEC__)
113  4;
114 #elif DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 3 && defined(__AVX512F__)
115  16;
116 #elif DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 2 && defined(__AVX__)
117  8;
118 #elif DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__SSE2__)
119  4;
120 #else
121  1;
122 #endif
123  };
124 
125 
126 } // namespace internal
127 
128 // forward declarations to support abs or sqrt operations on VectorizedArray
129 #ifndef DOXYGEN
130 template <typename Number,
131  int width =
133 class VectorizedArray;
134 template <typename T>
135 struct EnableIfScalar;
136 #endif
137 
138 DEAL_II_NAMESPACE_CLOSE
139 
140 // Declare / Import auto-differentiable math functions in(to) standard
141 // namespace before numbers::NumberTraits is defined
142 #ifdef DEAL_II_WITH_ADOLC
143 # include <deal.II/differentiation/ad/adolc_math.h>
144 
145 # include <adolc/adouble.h> // Taped double
146 #endif
147 // Ideally we'd like to #include <deal.II/differentiation/ad/sacado_math.h>
148 // but header indirectly references numbers.h. We therefore simply
149 // import the whole Sacado header at this point to get the math
150 // functions imported into the standard namespace.
151 #ifdef DEAL_II_TRILINOS_WITH_SACADO
152 # include <Sacado.hpp>
153 #endif
154 
155 namespace std
156 {
157  template <typename Number, int width>
158  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
159  sqrt(const ::VectorizedArray<Number, width> &);
160  template <typename Number, int width>
161  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
162  abs(const ::VectorizedArray<Number, width> &);
163  template <typename Number, int width>
164  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
165  max(const ::VectorizedArray<Number, width> &,
166  const ::VectorizedArray<Number, width> &);
167  template <typename Number, int width>
168  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
169  min(const ::VectorizedArray<Number, width> &,
170  const ::VectorizedArray<Number, width> &);
171  template <typename Number, int width>
173  pow(const ::VectorizedArray<Number, width> &, const Number p);
174  template <typename Number, int width>
176  sin(const ::VectorizedArray<Number, width> &);
177  template <typename Number, int width>
179  cos(const ::VectorizedArray<Number, width> &);
180  template <typename Number, int width>
182  tan(const ::VectorizedArray<Number, width> &);
183  template <typename Number, int width>
185  exp(const ::VectorizedArray<Number, width> &);
186  template <typename Number, int width>
188  log(const ::VectorizedArray<Number, width> &);
189 } // namespace std
190 
191 DEAL_II_NAMESPACE_OPEN
192 
208 namespace numbers
209 {
213  static constexpr double E = 2.7182818284590452354;
214 
218  static constexpr double LOG2E = 1.4426950408889634074;
219 
223  static constexpr double LOG10E = 0.43429448190325182765;
224 
228  static constexpr double LN2 = 0.69314718055994530942;
229 
233  static constexpr double LN10 = 2.30258509299404568402;
234 
238  static constexpr double PI = 3.14159265358979323846;
239 
243  static constexpr double PI_2 = 1.57079632679489661923;
244 
248  static constexpr double PI_4 = 0.78539816339744830962;
249 
253  static constexpr double SQRT2 = 1.41421356237309504880;
254 
258  static constexpr double SQRT1_2 = 0.70710678118654752440;
259 
265  template <typename Number, typename = void>
266  struct is_cuda_compatible : std::true_type
267  {};
268 
272  template <typename Number>
273  struct is_cuda_compatible<std::complex<Number>, void> : std::false_type
274  {};
275 
289  DEAL_II_DEPRECATED
290  bool
291  is_nan(const double x);
292 
302  bool
303  is_finite(const double x);
304 
309  bool
310  is_finite(const std::complex<double> &x);
311 
316  bool
317  is_finite(const std::complex<float> &x);
318 
327  bool
328  is_finite(const std::complex<long double> &x);
329 
340  template <typename Number1, typename Number2>
341  constexpr bool
342  values_are_equal(const Number1 &value_1, const Number2 &value_2);
343 
354  template <typename Number1, typename Number2>
355  bool
356  values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
357 
365  template <typename Number>
366  constexpr bool
367  value_is_zero(const Number &value);
368 
379  template <typename Number1, typename Number2>
380  bool
381  value_is_less_than(const Number1 &value_1, const Number2 &value_2);
382 
393  template <typename Number1, typename Number2>
394  bool
395  value_is_less_than_or_equal_to(const Number1 &value_1,
396  const Number2 &value_2);
397 
398 
399 
410  template <typename Number1, typename Number2>
411  bool
412  value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
413 
424  template <typename Number1, typename Number2>
425  bool
426  value_is_greater_than_or_equal_to(const Number1 &value_1,
427  const Number2 &value_2);
428 
439  template <typename number>
441  {
447  static constexpr bool is_complex = false;
448 
455  using real_type = number;
456 
464  static constexpr DEAL_II_CUDA_HOST_DEV const number &
465  conjugate(const number &x);
466 
475  template <typename Dummy = number>
476  static constexpr DEAL_II_CUDA_HOST_DEV
477  typename std::enable_if<std::is_same<Dummy, number>::value &&
479  real_type>::type
480  abs_square(const number &x);
481 
482  template <typename Dummy = number>
483  static constexpr
484  typename std::enable_if<std::is_same<Dummy, number>::value &&
485  !is_cuda_compatible<Dummy>::value,
486  real_type>::type
487  abs_square(const number &x);
488 
492  static real_type
493  abs(const number &x);
494  };
495 
496 
503  template <typename number>
504  struct NumberTraits<std::complex<number>>
505  {
511  static constexpr bool is_complex = true;
512 
519  using real_type = number;
520 
524  static constexpr std::complex<number>
525  conjugate(const std::complex<number> &x);
526 
533  static constexpr real_type
534  abs_square(const std::complex<number> &x);
535 
536 
540  static real_type
541  abs(const std::complex<number> &x);
542  };
543 
544  // --------------- inline and template functions ---------------- //
545 
546  inline bool
547  is_nan(const double x)
548  {
549  return std::isnan(x);
550  }
551 
552 
553 
554  inline bool
555  is_finite(const double x)
556  {
557  return std::isfinite(x);
558  }
559 
560 
561 
562  inline bool
563  is_finite(const std::complex<double> &x)
564  {
565  // Check complex numbers for infinity
566  // by testing real and imaginary part
567  return (is_finite(x.real()) && is_finite(x.imag()));
568  }
569 
570 
571 
572  inline bool
573  is_finite(const std::complex<float> &x)
574  {
575  // Check complex numbers for infinity
576  // by testing real and imaginary part
577  return (is_finite(x.real()) && is_finite(x.imag()));
578  }
579 
580 
581 
582  inline bool
583  is_finite(const std::complex<long double> &x)
584  {
585  // Same for std::complex<long double>
586  return (is_finite(x.real()) && is_finite(x.imag()));
587  }
588 
589 
590  template <typename number>
591  constexpr DEAL_II_CUDA_HOST_DEV const number &
593  {
594  return x;
595  }
596 
597 
598 
599  template <typename number>
600  template <typename Dummy>
601  constexpr DEAL_II_CUDA_HOST_DEV
602  typename std::enable_if<std::is_same<Dummy, number>::value &&
604  typename NumberTraits<number>::real_type>::type
606  {
607  return x * x;
608  }
609 
610 
611 
612  template <typename number>
613  template <typename Dummy>
614  constexpr
615  typename std::enable_if<std::is_same<Dummy, number>::value &&
616  !is_cuda_compatible<Dummy>::value,
617  typename NumberTraits<number>::real_type>::type
618  NumberTraits<number>::abs_square(const number &x)
619  {
620  return x * x;
621  }
622 
623 
624 
625  template <typename number>
627  NumberTraits<number>::abs(const number &x)
628  {
629  return std::abs(x);
630  }
631 
632 
633 
634  template <typename number>
635  constexpr std::complex<number>
636  NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
637  {
638  return std::conj(x);
639  }
640 
641 
642 
643  template <typename number>
645  NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
646  {
647  return std::abs(x);
648  }
649 
650 
651 
652  template <typename number>
654  NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
655  {
656  return std::norm(x);
657  }
658 
659 } // namespace numbers
660 
661 
662 // Forward declarations
664 {
665  namespace AD
666  {
667  namespace internal
668  {
669  // Defined in differentiation/ad/ad_number_traits.h
670  template <typename T>
671  struct NumberType;
672  } // namespace internal
673 
674  // Defined in differentiation/ad/ad_number_traits.h
675  template <typename NumberType>
676  struct is_ad_number;
677  } // namespace AD
678 } // namespace Differentiation
679 
680 
681 namespace internal
682 {
687  template <typename From, typename To>
689  {
690  // Source: https://stackoverflow.com/a/16944130
691  private:
692  template <typename T>
693  static void f(T);
694 
695  template <typename F, typename T>
696  static constexpr auto
697  test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
698  {
699  return true;
700  }
701 
702  template <typename F, typename T>
703  static constexpr auto
704  test(...) -> bool
705  {
706  return false;
707  }
708 
709  public:
710  static bool const value = test<From, To>(0);
711  };
712 
713  /*
714  * The structs below are needed to convert between some special number types.
715  * Also see tensor.h for another specialization.
716  */
717  template <typename T>
718  struct NumberType
719  {
720  static constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV const T &
721  value(const T &t)
722  {
723  return t;
724  }
725 
726  // Below are generic functions that allows an overload for any
727  // type U that is transformable to type T. This is particularly
728  // useful when needing to cast exotic number types
729  // (e.g. auto-differentiable or symbolic numbers) to a floating
730  // point one, such as might happen when converting between tensor
731  // types.
732 
733  // Type T is constructible from F.
734  template <typename F>
735  static constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV 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>::type * = nullptr)
741  {
742  return T(f);
743  }
744 
745  // Type T is explicitly convertible (but not constructible) from F.
746  template <typename F>
747  static constexpr DEAL_II_ALWAYS_INLINE 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 &&
754  {
755  return static_cast<T>(f);
756  }
757 
758  // Sacado doesn't provide any conversion operators, so we have
759  // to extract the value and perform further conversions from there.
760  // To be safe, we extend this to other possible AD numbers that
761  // might fall into the same category.
762  template <typename F>
763  static T
764  value(const F &f,
765  typename std::enable_if<
766  !std::is_same<typename std::decay<T>::type,
767  typename std::decay<F>::type>::value &&
768  !std::is_constructible<T, F>::value &&
771  {
773  }
774  };
775 
776  template <typename T>
777  struct NumberType<std::complex<T>>
778  {
779  static constexpr const std::complex<T> &
780  value(const std::complex<T> &t)
781  {
782  return t;
783  }
784 
785  static constexpr std::complex<T>
786  value(const T &t)
787  {
788  return std::complex<T>(t);
789  }
790 
791  // Facilitate cast from complex<double> to complex<float>
792  template <typename U>
793  static constexpr std::complex<T>
794  value(const std::complex<U> &t)
795  {
796  return std::complex<T>(NumberType<T>::value(t.real()),
797  NumberType<T>::value(t.imag()));
798  }
799  };
800 
801 #ifdef DEAL_II_COMPILER_CUDA_AWARE
802  template <>
803  struct NumberType<cuComplex>
804  {
805  static cuComplex
806  value(const float t)
807  {
808  return make_cuComplex(t, 0.f);
809  }
810  };
811 
812  template <>
813  struct NumberType<cuDoubleComplex>
814  {
815  static cuDoubleComplex
816  value(const double t)
817  {
818  return make_cuDoubleComplex(t, 0.);
819  }
820  };
821 #endif
822 } // namespace internal
823 
824 namespace numbers
825 {
826 #ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
827 
838  // Defined in differentiation/ad/adolc_number_types.cc
839  bool
840  values_are_equal(const adouble &value_1, const adouble &value_2);
841 
842 
853  template <typename Number>
854  bool
855  values_are_equal(const adouble &value_1, const Number &value_2)
856  {
857  // Use the specialized definition for two ADOL-C taped types
858  return values_are_equal(value_1,
859  internal::NumberType<adouble>::value(value_2));
860  }
861 
862 
873  template <typename Number>
874  bool
875  values_are_equal(const Number &value_1, const adouble &value_2)
876  {
877  // Use the above definition
878  return values_are_equal(value_2, value_1);
879  }
880 
892  // Defined in differentiation/ad/adolc_number_types.cc
893  bool
894  value_is_less_than(const adouble &value_1, const adouble &value_2);
895 
896 
908  template <typename Number>
909  bool
910  value_is_less_than(const adouble &value_1, const Number &value_2)
911  {
912  // Use the specialized definition for two ADOL-C taped types
913  return value_is_less_than(value_1,
914  internal::NumberType<adouble>::value(value_2));
915  }
916 
917 
929  template <typename Number>
930  bool
931  value_is_less_than(const Number &value_1, const adouble &value_2)
932  {
933  // Use the specialized definition for two ADOL-C taped types
934  return value_is_less_than(internal::NumberType<adouble>::value(value_1),
935  value_2);
936  }
937 
938 #endif
939 
940 
941  template <typename Number1, typename Number2>
942  constexpr bool
943  values_are_equal(const Number1 &value_1, const Number2 &value_2)
944  {
945  return (value_1 == internal::NumberType<Number1>::value(value_2));
946  }
947 
948 
949  template <typename Number1, typename Number2>
950  inline bool
951  values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
952  {
953  return !(values_are_equal(value_1, value_2));
954  }
955 
956 
957  template <typename Number>
958  constexpr bool
959  value_is_zero(const Number &value)
960  {
961  return values_are_equal(value, 0.0);
962  }
963 
964 
965  template <typename Number1, typename Number2>
966  inline bool
967  value_is_less_than(const Number1 &value_1, const Number2 &value_2)
968  {
969  return (value_1 < internal::NumberType<Number1>::value(value_2));
970  }
971 
972 
973  template <typename Number1, typename Number2>
974  inline bool
975  value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
976  {
977  return (value_is_less_than(value_1, value_2) ||
978  values_are_equal(value_1, value_2));
979  }
980 
981 
982  template <typename Number1, typename Number2>
983  bool
984  value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
985  {
986  return !(value_is_less_than_or_equal_to(value_1, value_2));
987  }
988 
989 
990  template <typename Number1, typename Number2>
991  inline bool
992  value_is_greater_than_or_equal_to(const Number1 &value_1,
993  const Number2 &value_2)
994  {
995  return !(value_is_less_than(value_1, value_2));
996  }
997 } // namespace numbers
998 
999 DEAL_II_NAMESPACE_CLOSE
1000 
1001 #endif
static constexpr double LOG2E
Definition: numbers.h:218
static constexpr unsigned int max_width
Definition: numbers.h:66
bool value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:975
static constexpr double LN10
Definition: numbers.h:233
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:943
STL namespace.
static constexpr double LOG10E
Definition: numbers.h:223
bool is_finite(const double x)
Definition: numbers.h:555
static constexpr double SQRT1_2
Definition: numbers.h:258
static constexpr double E
Definition: numbers.h:213
static constexpr double SQRT2
Definition: numbers.h:253
bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:951
static constexpr double PI_2
Definition: numbers.h:243
bool is_nan(const double x)
Definition: numbers.h:547
bool value_is_greater_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:992
static constexpr double LN2
Definition: numbers.h:228
bool value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:984
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:967
static constexpr double PI
Definition: numbers.h:238
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:959
static constexpr double PI_4
Definition: numbers.h:248