Reference documentation for deal.II version Git f40be01994 2020-04-09 07:13:12 +0200
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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_VECTORIZATION_WIDTH_IN_BITS >= 512
85  8;
86 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256
87  4;
88 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128
89  2;
90 #else
91  1;
92 #endif
93  };
94 
103  template <>
105  {
109  constexpr static unsigned int max_width =
110 #if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ALTIVEC__)
111  4;
112 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && defined(__AVX512F__)
113  16;
114 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)
115  8;
116 #elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__)
117  4;
118 #else
119  1;
120 #endif
121  };
122 
123 
124 } // namespace internal
125 
126 // forward declarations to support abs or sqrt operations on VectorizedArray
127 #ifndef DOXYGEN
128 template <typename Number,
129  std::size_t width =
131 class VectorizedArray;
132 template <typename T>
133 struct EnableIfScalar;
134 #endif
135 
136 DEAL_II_NAMESPACE_CLOSE
137 
138 // Declare / Import auto-differentiable math functions in(to) standard
139 // namespace before numbers::NumberTraits is defined
140 #ifdef DEAL_II_WITH_ADOLC
141 # include <deal.II/differentiation/ad/adolc_math.h>
142 
143 # include <adolc/adouble.h> // Taped double
144 #endif
145 // Ideally we'd like to #include <deal.II/differentiation/ad/sacado_math.h>
146 // but header indirectly references numbers.h. We therefore simply
147 // import the whole Sacado header at this point to get the math
148 // functions imported into the standard namespace.
149 #ifdef DEAL_II_TRILINOS_WITH_SACADO
150 # include <Sacado.hpp>
151 #endif
152 
153 namespace std
154 {
155  template <typename Number, std::size_t width>
156  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
157  sqrt(const ::VectorizedArray<Number, width> &);
158  template <typename Number, std::size_t width>
159  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
160  abs(const ::VectorizedArray<Number, width> &);
161  template <typename Number, std::size_t width>
162  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
163  max(const ::VectorizedArray<Number, width> &,
164  const ::VectorizedArray<Number, width> &);
165  template <typename Number, std::size_t width>
166  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
167  min(const ::VectorizedArray<Number, width> &,
168  const ::VectorizedArray<Number, width> &);
169  template <typename Number, size_t width>
171  pow(const ::VectorizedArray<Number, width> &, const Number p);
172  template <typename Number, size_t width>
174  sin(const ::VectorizedArray<Number, width> &);
175  template <typename Number, size_t width>
177  cos(const ::VectorizedArray<Number, width> &);
178  template <typename Number, size_t width>
180  tan(const ::VectorizedArray<Number, width> &);
181  template <typename Number, size_t width>
183  exp(const ::VectorizedArray<Number, width> &);
184  template <typename Number, size_t width>
186  log(const ::VectorizedArray<Number, width> &);
187 } // namespace std
188 
189 DEAL_II_NAMESPACE_OPEN
190 
206 namespace numbers
207 {
211  static constexpr double E = 2.7182818284590452354;
212 
216  static constexpr double LOG2E = 1.4426950408889634074;
217 
221  static constexpr double LOG10E = 0.43429448190325182765;
222 
226  static constexpr double LN2 = 0.69314718055994530942;
227 
231  static constexpr double LN10 = 2.30258509299404568402;
232 
236  static constexpr double PI = 3.14159265358979323846;
237 
241  static constexpr double PI_2 = 1.57079632679489661923;
242 
246  static constexpr double PI_4 = 0.78539816339744830962;
247 
251  static constexpr double SQRT2 = 1.41421356237309504880;
252 
256  static constexpr double SQRT1_2 = 0.70710678118654752440;
257 
263  template <typename Number, typename = void>
264  struct is_cuda_compatible : std::true_type
265  {};
266 
270  template <typename Number>
271  struct is_cuda_compatible<std::complex<Number>, void> : std::false_type
272  {};
273 
287  DEAL_II_DEPRECATED
288  bool
289  is_nan(const double x);
290 
300  bool
301  is_finite(const double x);
302 
307  bool
308  is_finite(const std::complex<double> &x);
309 
314  bool
315  is_finite(const std::complex<float> &x);
316 
325  bool
326  is_finite(const std::complex<long double> &x);
327 
338  template <typename Number1, typename Number2>
339  constexpr bool
340  values_are_equal(const Number1 &value_1, const Number2 &value_2);
341 
352  template <typename Number1, typename Number2>
353  bool
354  values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
355 
363  template <typename Number>
364  constexpr bool
365  value_is_zero(const Number &value);
366 
377  template <typename Number1, typename Number2>
378  bool
379  value_is_less_than(const Number1 &value_1, const Number2 &value_2);
380 
391  template <typename Number1, typename Number2>
392  bool
393  value_is_less_than_or_equal_to(const Number1 &value_1,
394  const Number2 &value_2);
395 
396 
397 
408  template <typename Number1, typename Number2>
409  bool
410  value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
411 
422  template <typename Number1, typename Number2>
423  bool
424  value_is_greater_than_or_equal_to(const Number1 &value_1,
425  const Number2 &value_2);
426 
437  template <typename number>
439  {
445  static constexpr bool is_complex = false;
446 
453  using real_type = number;
454 
462  static constexpr DEAL_II_CUDA_HOST_DEV const number &
463  conjugate(const number &x);
464 
473  template <typename Dummy = number>
474  static constexpr DEAL_II_CUDA_HOST_DEV
475  typename std::enable_if<std::is_same<Dummy, number>::value &&
477  real_type>::type
478  abs_square(const number &x);
479 
480  template <typename Dummy = number>
481  static constexpr
482  typename std::enable_if<std::is_same<Dummy, number>::value &&
483  !is_cuda_compatible<Dummy>::value,
484  real_type>::type
485  abs_square(const number &x);
486 
490  static real_type
491  abs(const number &x);
492  };
493 
494 
501  template <typename number>
502  struct NumberTraits<std::complex<number>>
503  {
509  static constexpr bool is_complex = true;
510 
517  using real_type = number;
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  typename std::enable_if<std::is_same<Dummy, number>::value &&
602  typename NumberTraits<number>::real_type>::type
604  {
605  return x * x;
606  }
607 
608 
609 
610  template <typename number>
611  template <typename Dummy>
612  constexpr
613  typename std::enable_if<std::is_same<Dummy, number>::value &&
614  !is_cuda_compatible<Dummy>::value,
615  typename NumberTraits<number>::real_type>::type
616  NumberTraits<number>::abs_square(const number &x)
617  {
618  return x * x;
619  }
620 
621 
622 
623  template <typename number>
625  NumberTraits<number>::abs(const number &x)
626  {
627  return std::abs(x);
628  }
629 
630 
631 
632  template <typename number>
633  constexpr std::complex<number>
634  NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
635  {
636  return std::conj(x);
637  }
638 
639 
640 
641  template <typename number>
643  NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
644  {
645  return std::abs(x);
646  }
647 
648 
649 
650  template <typename number>
652  NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
653  {
654  return std::norm(x);
655  }
656 
657 } // namespace numbers
658 
659 
660 // Forward declarations
662 {
663  namespace AD
664  {
665  namespace internal
666  {
667  // Defined in differentiation/ad/ad_number_traits.h
668  template <typename T>
669  struct NumberType;
670  } // namespace internal
671 
672  // Defined in differentiation/ad/ad_number_traits.h
673  template <typename NumberType>
674  struct is_ad_number;
675  } // namespace AD
676 } // namespace Differentiation
677 
678 
679 namespace internal
680 {
685  template <typename From, typename To>
687  {
688  // Source: https://stackoverflow.com/a/16944130
689  private:
690  template <typename T>
691  static void f(T);
692 
693  template <typename F, typename T>
694  static constexpr auto
695  test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
696  {
697  return true;
698  }
699 
700  template <typename F, typename T>
701  static constexpr auto
702  test(...) -> bool
703  {
704  return false;
705  }
706 
707  public:
708  static bool const value = test<From, To>(0);
709  };
710 
711  /*
712  * The structs below are needed to convert between some special number types.
713  * Also see tensor.h for another specialization.
714  */
715  template <typename T>
716  struct NumberType
717  {
718  static constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV const T &
719  value(const T &t)
720  {
721  return t;
722  }
723 
724  // Below are generic functions that allows an overload for any
725  // type U that is transformable to type T. This is particularly
726  // useful when needing to cast exotic number types
727  // (e.g. auto-differentiable or symbolic numbers) to a floating
728  // point one, such as might happen when converting between tensor
729  // types.
730 
731  // Type T is constructible from F.
732  template <typename F>
733  static constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV T
734  value(const F &f,
735  typename std::enable_if<
736  !std::is_same<typename std::decay<T>::type,
737  typename std::decay<F>::type>::value &&
738  std::is_constructible<T, F>::value>::type * = nullptr)
739  {
740  return T(f);
741  }
742 
743  // Type T is explicitly convertible (but not constructible) from F.
744  template <typename F>
745  static constexpr DEAL_II_ALWAYS_INLINE 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 &&
750  !std::is_constructible<T, F>::value &&
752  {
753  return static_cast<T>(f);
754  }
755 
756  // Sacado doesn't provide any conversion operators, so we have
757  // to extract the value and perform further conversions from there.
758  // To be safe, we extend this to other possible AD numbers that
759  // might fall into the same category.
760  template <typename F>
761  static T
762  value(const F &f,
763  typename std::enable_if<
764  !std::is_same<typename std::decay<T>::type,
765  typename std::decay<F>::type>::value &&
766  !std::is_constructible<T, F>::value &&
769  {
771  }
772  };
773 
774  template <typename T>
775  struct NumberType<std::complex<T>>
776  {
777  static constexpr const std::complex<T> &
778  value(const std::complex<T> &t)
779  {
780  return t;
781  }
782 
783  static constexpr std::complex<T>
784  value(const T &t)
785  {
786  return std::complex<T>(t);
787  }
788 
789  // Facilitate cast from complex<double> to complex<float>
790  template <typename U>
791  static constexpr std::complex<T>
792  value(const std::complex<U> &t)
793  {
794  return std::complex<T>(NumberType<T>::value(t.real()),
795  NumberType<T>::value(t.imag()));
796  }
797  };
798 
799 #ifdef DEAL_II_COMPILER_CUDA_AWARE
800  template <>
801  struct NumberType<cuComplex>
802  {
803  static cuComplex
804  value(const float t)
805  {
806  return make_cuComplex(t, 0.f);
807  }
808  };
809 
810  template <>
811  struct NumberType<cuDoubleComplex>
812  {
813  static cuDoubleComplex
814  value(const double t)
815  {
816  return make_cuDoubleComplex(t, 0.);
817  }
818  };
819 #endif
820 } // namespace internal
821 
822 namespace numbers
823 {
824 #ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
825 
836  // Defined in differentiation/ad/adolc_number_types.cc
837  bool
838  values_are_equal(const adouble &value_1, const adouble &value_2);
839 
840 
851  template <typename Number>
852  bool
853  values_are_equal(const adouble &value_1, const Number &value_2)
854  {
855  // Use the specialized definition for two ADOL-C taped types
856  return values_are_equal(value_1,
857  internal::NumberType<adouble>::value(value_2));
858  }
859 
860 
871  template <typename Number>
872  bool
873  values_are_equal(const Number &value_1, const adouble &value_2)
874  {
875  // Use the above definition
876  return values_are_equal(value_2, value_1);
877  }
878 
890  // Defined in differentiation/ad/adolc_number_types.cc
891  bool
892  value_is_less_than(const adouble &value_1, const adouble &value_2);
893 
894 
906  template <typename Number>
907  bool
908  value_is_less_than(const adouble &value_1, const Number &value_2)
909  {
910  // Use the specialized definition for two ADOL-C taped types
911  return value_is_less_than(value_1,
912  internal::NumberType<adouble>::value(value_2));
913  }
914 
915 
927  template <typename Number>
928  bool
929  value_is_less_than(const Number &value_1, const adouble &value_2)
930  {
931  // Use the specialized definition for two ADOL-C taped types
932  return value_is_less_than(internal::NumberType<adouble>::value(value_1),
933  value_2);
934  }
935 
936 #endif
937 
938 
939  template <typename Number1, typename Number2>
940  constexpr bool
941  values_are_equal(const Number1 &value_1, const Number2 &value_2)
942  {
943  return (value_1 == internal::NumberType<Number1>::value(value_2));
944  }
945 
946 
947  template <typename Number1, typename Number2>
948  inline bool
949  values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
950  {
951  return !(values_are_equal(value_1, value_2));
952  }
953 
954 
955  template <typename Number>
956  constexpr bool
957  value_is_zero(const Number &value)
958  {
959  return values_are_equal(value, 0.0);
960  }
961 
962 
963  template <typename Number1, typename Number2>
964  inline bool
965  value_is_less_than(const Number1 &value_1, const Number2 &value_2)
966  {
967  return (value_1 < internal::NumberType<Number1>::value(value_2));
968  }
969 
970 
971  template <typename Number1, typename Number2>
972  inline bool
973  value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
974  {
975  return (value_is_less_than(value_1, value_2) ||
976  values_are_equal(value_1, value_2));
977  }
978 
979 
980  template <typename Number1, typename Number2>
981  bool
982  value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
983  {
984  return !(value_is_less_than_or_equal_to(value_1, value_2));
985  }
986 
987 
988  template <typename Number1, typename Number2>
989  inline bool
990  value_is_greater_than_or_equal_to(const Number1 &value_1,
991  const Number2 &value_2)
992  {
993  return !(value_is_less_than(value_1, value_2));
994  }
995 } // namespace numbers
996 
997 DEAL_II_NAMESPACE_CLOSE
998 
999 #endif
static constexpr double LOG2E
Definition: numbers.h:216
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:973
static constexpr double LN10
Definition: numbers.h:231
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:941
STL namespace.
static constexpr double LOG10E
Definition: numbers.h:221
bool is_finite(const double x)
Definition: numbers.h:553
static constexpr double SQRT1_2
Definition: numbers.h:256
static constexpr double E
Definition: numbers.h:211
static constexpr double SQRT2
Definition: numbers.h:251
bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:949
static constexpr double PI_2
Definition: numbers.h:241
bool is_nan(const double x)
Definition: numbers.h:545
bool value_is_greater_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:990
static constexpr double LN2
Definition: numbers.h:226
bool value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:982
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:965
static constexpr double PI
Definition: numbers.h:236
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:957
static constexpr double PI_4
Definition: numbers.h:246