Reference documentation for deal.II version Git d199cacc3e 2019-07-18 17:13:56 -0400
\(\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 {
42  template <typename Number>
43  struct VectorizedArrayWidthSpecifier
44  {
45  static const unsigned int max_width = 1;
46  };
47 
48  template <>
49  struct VectorizedArrayWidthSpecifier<double>
50  {
51  static const unsigned int max_width =
52 #if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__ALTIVEC__)
53  2;
54 #elif DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 3 && defined(__AVX512F__)
55  8;
56 #elif DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 2 && defined(__AVX__)
57  4;
58 #elif DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__SSE2__)
59  2;
60 #else
61  1;
62 #endif
63  };
64 
65  template <>
66  struct VectorizedArrayWidthSpecifier<float>
67  {
68  static const unsigned int max_width =
69 #if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__ALTIVEC__)
70  4;
71 #elif DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 3 && defined(__AVX512F__)
72  16;
73 #elif DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 2 && defined(__AVX__)
74  8;
75 #elif DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1 && defined(__SSE2__)
76  4;
77 #else
78  1;
79 #endif
80  };
81 
82 
83 } // namespace internal
84 
85 // forward declarations to support abs or sqrt operations on VectorizedArray
86 template <typename Number,
87  int width =
88  internal::VectorizedArrayWidthSpecifier<Number>::max_width>
90 template <typename T>
91 struct EnableIfScalar;
92 
93 DEAL_II_NAMESPACE_CLOSE
94 
95 // Declare / Import auto-differentiable math functions in(to) standard
96 // namespace before numbers::NumberTraits is defined
97 #ifdef DEAL_II_WITH_ADOLC
98 # include <deal.II/differentiation/ad/adolc_math.h>
99 
100 # include <adolc/adouble.h> // Taped double
101 #endif
102 // Ideally we'd like to #include <deal.II/differentiation/ad/sacado_math.h>
103 // but header indirectly references numbers.h. We therefore simply
104 // import the whole Sacado header at this point to get the math
105 // functions imported into the standard namespace.
106 #ifdef DEAL_II_TRILINOS_WITH_SACADO
107 # include <Sacado.hpp>
108 #endif
109 
110 namespace std
111 {
112  template <typename Number, int width>
113  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
114  sqrt(const ::VectorizedArray<Number, width> &);
115  template <typename Number, int width>
116  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
117  abs(const ::VectorizedArray<Number, width> &);
118  template <typename Number, int width>
119  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
120  max(const ::VectorizedArray<Number, width> &,
121  const ::VectorizedArray<Number, width> &);
122  template <typename Number, int width>
123  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
124  min(const ::VectorizedArray<Number, width> &,
125  const ::VectorizedArray<Number, width> &);
126  template <typename Number, int width>
128  pow(const ::VectorizedArray<Number, width> &, const Number p);
129  template <typename Number, int width>
131  sin(const ::VectorizedArray<Number, width> &);
132  template <typename Number, int width>
134  cos(const ::VectorizedArray<Number, width> &);
135  template <typename Number, int width>
137  tan(const ::VectorizedArray<Number, width> &);
138  template <typename Number, int width>
140  exp(const ::VectorizedArray<Number, width> &);
141  template <typename Number, int width>
143  log(const ::VectorizedArray<Number, width> &);
144 } // namespace std
145 
146 DEAL_II_NAMESPACE_OPEN
147 
163 namespace numbers
164 {
168  static constexpr double E = 2.7182818284590452354;
169 
173  static constexpr double LOG2E = 1.4426950408889634074;
174 
178  static constexpr double LOG10E = 0.43429448190325182765;
179 
183  static constexpr double LN2 = 0.69314718055994530942;
184 
188  static constexpr double LN10 = 2.30258509299404568402;
189 
193  static constexpr double PI = 3.14159265358979323846;
194 
198  static constexpr double PI_2 = 1.57079632679489661923;
199 
203  static constexpr double PI_4 = 0.78539816339744830962;
204 
208  static constexpr double SQRT2 = 1.41421356237309504880;
209 
213  static constexpr double SQRT1_2 = 0.70710678118654752440;
214 
220  template <typename Number, typename = void>
221  struct is_cuda_compatible : std::true_type
222  {};
223 
227  template <typename Number>
228  struct is_cuda_compatible<std::complex<Number>, void> : std::false_type
229  {};
230 
244  DEAL_II_DEPRECATED
245  bool
246  is_nan(const double x);
247 
257  bool
258  is_finite(const double x);
259 
264  bool
265  is_finite(const std::complex<double> &x);
266 
271  bool
272  is_finite(const std::complex<float> &x);
273 
282  bool
283  is_finite(const std::complex<long double> &x);
284 
295  template <typename Number1, typename Number2>
296  bool
297  values_are_equal(const Number1 &value_1, const Number2 &value_2);
298 
309  template <typename Number1, typename Number2>
310  bool
311  values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
312 
320  template <typename Number>
321  bool
322  value_is_zero(const Number &value);
323 
334  template <typename Number1, typename Number2>
335  bool
336  value_is_less_than(const Number1 &value_1, const Number2 &value_2);
337 
348  template <typename Number1, typename Number2>
349  bool
350  value_is_less_than_or_equal_to(const Number1 &value_1,
351  const Number2 &value_2);
352 
353 
354 
365  template <typename Number1, typename Number2>
366  bool
367  value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
368 
379  template <typename Number1, typename Number2>
380  bool
381  value_is_greater_than_or_equal_to(const Number1 &value_1,
382  const Number2 &value_2);
383 
394  template <typename number>
396  {
402  static constexpr bool is_complex = false;
403 
410  using real_type = number;
411 
419  static constexpr DEAL_II_CUDA_HOST_DEV const number &
420  conjugate(const number &x);
421 
430  template <typename Dummy = number>
431  static constexpr DEAL_II_CUDA_HOST_DEV
432  typename std::enable_if<std::is_same<Dummy, number>::value &&
434  real_type>::type
435  abs_square(const number &x);
436 
437  template <typename Dummy = number>
438  static constexpr
439  typename std::enable_if<std::is_same<Dummy, number>::value &&
440  !is_cuda_compatible<Dummy>::value,
441  real_type>::type
442  abs_square(const number &x);
443 
447  static real_type
448  abs(const number &x);
449  };
450 
451 
458  template <typename number>
459  struct NumberTraits<std::complex<number>>
460  {
466  static constexpr bool is_complex = true;
467 
474  using real_type = number;
475 
479  static constexpr std::complex<number>
480  conjugate(const std::complex<number> &x);
481 
488  static constexpr real_type
489  abs_square(const std::complex<number> &x);
490 
491 
495  static real_type
496  abs(const std::complex<number> &x);
497  };
498 
499  // --------------- inline and template functions ---------------- //
500 
501  inline bool
502  is_nan(const double x)
503  {
504  return std::isnan(x);
505  }
506 
507 
508 
509  inline bool
510  is_finite(const double x)
511  {
512  return std::isfinite(x);
513  }
514 
515 
516 
517  inline bool
518  is_finite(const std::complex<double> &x)
519  {
520  // Check complex numbers for infinity
521  // by testing real and imaginary part
522  return (is_finite(x.real()) && is_finite(x.imag()));
523  }
524 
525 
526 
527  inline bool
528  is_finite(const std::complex<float> &x)
529  {
530  // Check complex numbers for infinity
531  // by testing real and imaginary part
532  return (is_finite(x.real()) && is_finite(x.imag()));
533  }
534 
535 
536 
537  inline bool
538  is_finite(const std::complex<long double> &x)
539  {
540  // Same for std::complex<long double>
541  return (is_finite(x.real()) && is_finite(x.imag()));
542  }
543 
544 
545  template <typename number>
546  constexpr DEAL_II_CUDA_HOST_DEV const number &
548  {
549  return x;
550  }
551 
552 
553 
554  template <typename number>
555  template <typename Dummy>
556  constexpr DEAL_II_CUDA_HOST_DEV
557  typename std::enable_if<std::is_same<Dummy, number>::value &&
559  typename NumberTraits<number>::real_type>::type
561  {
562  return x * x;
563  }
564 
565 
566 
567  template <typename number>
568  template <typename Dummy>
569  constexpr
570  typename std::enable_if<std::is_same<Dummy, number>::value &&
571  !is_cuda_compatible<Dummy>::value,
572  typename NumberTraits<number>::real_type>::type
573  NumberTraits<number>::abs_square(const number &x)
574  {
575  return x * x;
576  }
577 
578 
579 
580  template <typename number>
582  NumberTraits<number>::abs(const number &x)
583  {
584  return std::abs(x);
585  }
586 
587 
588 
589  template <typename number>
590  constexpr std::complex<number>
591  NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
592  {
593  return std::conj(x);
594  }
595 
596 
597 
598  template <typename number>
600  NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
601  {
602  return std::abs(x);
603  }
604 
605 
606 
607  template <typename number>
609  NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
610  {
611  return std::norm(x);
612  }
613 
614 } // namespace numbers
615 
616 
617 // Forward declarations
619 {
620  namespace AD
621  {
622  namespace internal
623  {
624  // Defined in differentiation/ad/ad_number_traits.h
625  template <typename T>
626  struct NumberType;
627  } // namespace internal
628 
629  // Defined in differentiation/ad/ad_number_traits.h
630  template <typename NumberType>
631  struct is_ad_number;
632  } // namespace AD
633 } // namespace Differentiation
634 
635 
636 namespace internal
637 {
642  template <typename From, typename To>
644  {
645  // Source: https://stackoverflow.com/a/16944130
646  private:
647  template <typename T>
648  static void f(T);
649 
650  template <typename F, typename T>
651  static constexpr auto
652  test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
653  {
654  return true;
655  }
656 
657  template <typename F, typename T>
658  static constexpr auto
659  test(...) -> bool
660  {
661  return false;
662  }
663 
664  public:
665  static bool const value = test<From, To>(0);
666  };
667 
681  template <typename T>
682  struct NumberType
683  {
684  static constexpr DEAL_II_CUDA_HOST_DEV const T &
685  value(const T &t)
686  {
687  return t;
688  }
689 
690  // Below are generic functions that allows an overload for any
691  // type U that is transformable to type T. This is particularly
692  // useful when needing to cast exotic number types
693  // (e.g. auto-differentiable or symbolic numbers) to a floating
694  // point one, such as might happen when converting between tensor
695  // types.
696 
697  // Type T is constructible from F.
698  template <typename F>
699  static constexpr DEAL_II_CUDA_HOST_DEV T
700  value(const F &f,
701  typename std::enable_if<
702  !std::is_same<typename std::decay<T>::type,
703  typename std::decay<F>::type>::value &&
704  std::is_constructible<T, F>::value>::type * = nullptr)
705  {
706  return T(f);
707  }
708 
709  // Type T is explicitly convertible (but not constructible) from F.
710  template <typename F>
711  static constexpr T
712  value(const F &f,
713  typename std::enable_if<
714  !std::is_same<typename std::decay<T>::type,
715  typename std::decay<F>::type>::value &&
716  !std::is_constructible<T, F>::value &&
718  {
719  return static_cast<T>(f);
720  }
721 
722  // Sacado doesn't provide any conversion operators, so we have
723  // to extract the value and perform further conversions from there.
724  // To be safe, we extend this to other possible AD numbers that
725  // might fall into the same category.
726  template <typename F>
727  static T
728  value(const F &f,
729  typename std::enable_if<
730  !std::is_same<typename std::decay<T>::type,
731  typename std::decay<F>::type>::value &&
732  !std::is_constructible<T, F>::value &&
735  {
737  }
738  };
739 
740  template <typename T>
741  struct NumberType<std::complex<T>>
742  {
743  static constexpr const std::complex<T> &
744  value(const std::complex<T> &t)
745  {
746  return t;
747  }
748 
749  static constexpr std::complex<T>
750  value(const T &t)
751  {
752  return std::complex<T>(t);
753  }
754 
755  // Facilitate cast from complex<double> to complex<float>
756  template <typename U>
757  static constexpr std::complex<T>
758  value(const std::complex<U> &t)
759  {
760  return std::complex<T>(NumberType<T>::value(t.real()),
761  NumberType<T>::value(t.imag()));
762  }
763  };
764 
765 #ifdef DEAL_II_COMPILER_CUDA_AWARE
766  template <>
767  struct NumberType<cuComplex>
768  {
769  static cuComplex
770  value(const float t)
771  {
772  return make_cuComplex(t, 0.f);
773  }
774  };
775 
776  template <>
777  struct NumberType<cuDoubleComplex>
778  {
779  static cuDoubleComplex
780  value(const double t)
781  {
782  return make_cuDoubleComplex(t, 0.);
783  }
784  };
785 #endif
786 } // namespace internal
787 
788 namespace numbers
789 {
790 #ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
791 
802  // Defined in differentiation/ad/adolc_number_types.cc
803  bool
804  values_are_equal(const adouble &value_1, const adouble &value_2);
805 
806 
817  template <typename Number>
818  bool
819  values_are_equal(const adouble &value_1, const Number &value_2)
820  {
821  // Use the specialized definition for two ADOL-C taped types
822  return values_are_equal(value_1,
824  }
825 
826 
837  template <typename Number>
838  bool
839  values_are_equal(const Number &value_1, const adouble &value_2)
840  {
841  // Use the above definition
842  return values_are_equal(value_2, value_1);
843  }
844 
856  // Defined in differentiation/ad/adolc_number_types.cc
857  bool
858  value_is_less_than(const adouble &value_1, const adouble &value_2);
859 
860 
872  template <typename Number>
873  bool
874  value_is_less_than(const adouble &value_1, const Number &value_2)
875  {
876  // Use the specialized definition for two ADOL-C taped types
877  return value_is_less_than(value_1,
879  }
880 
881 
893  template <typename Number>
894  bool
895  value_is_less_than(const Number &value_1, const adouble &value_2)
896  {
897  // Use the specialized definition for two ADOL-C taped types
898  return value_is_less_than(internal::NumberType<adouble>::value(value_1),
899  value_2);
900  }
901 
902 #endif
903 
904 
905  template <typename Number1, typename Number2>
906  inline bool
907  values_are_equal(const Number1 &value_1, const Number2 &value_2)
908  {
909  return (value_1 == internal::NumberType<Number1>::value(value_2));
910  }
911 
912 
913  template <typename Number1, typename Number2>
914  inline bool
915  values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
916  {
917  return !(values_are_equal(value_1, value_2));
918  }
919 
920 
921  template <typename Number>
922  inline bool
923  value_is_zero(const Number &value)
924  {
925  return values_are_equal(value, 0.0);
926  }
927 
928 
929  template <typename Number1, typename Number2>
930  inline bool
931  value_is_less_than(const Number1 &value_1, const Number2 &value_2)
932  {
933  return (value_1 < internal::NumberType<Number1>::value(value_2));
934  }
935 
936 
937  template <typename Number1, typename Number2>
938  inline bool
939  value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
940  {
941  return (value_is_less_than(value_1, value_2) ||
942  values_are_equal(value_1, value_2));
943  }
944 
945 
946  template <typename Number1, typename Number2>
947  bool
948  value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
949  {
950  return !(value_is_less_than_or_equal_to(value_1, value_2));
951  }
952 
953 
954  template <typename Number1, typename Number2>
955  inline bool
956  value_is_greater_than_or_equal_to(const Number1 &value_1,
957  const Number2 &value_2)
958  {
959  return !(value_is_less_than(value_1, value_2));
960  }
961 } // namespace numbers
962 
963 DEAL_II_NAMESPACE_CLOSE
964 
965 #endif
static constexpr double LOG2E
Definition: numbers.h:173
bool value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:939
static constexpr double LN10
Definition: numbers.h:188
bool value_is_zero(const Number &value)
Definition: numbers.h:923
STL namespace.
static constexpr double LOG10E
Definition: numbers.h:178
bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:907
bool is_finite(const double x)
Definition: numbers.h:510
static constexpr double SQRT1_2
Definition: numbers.h:213
static constexpr double E
Definition: numbers.h:168
static constexpr double SQRT2
Definition: numbers.h:208
bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:915
static constexpr double PI_2
Definition: numbers.h:198
bool is_nan(const double x)
Definition: numbers.h:502
bool value_is_greater_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:956
static constexpr double LN2
Definition: numbers.h:183
bool value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:948
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:931
static constexpr double PI
Definition: numbers.h:193
static constexpr double PI_4
Definition: numbers.h:203