Reference documentation for deal.II version 9.1.1
\(\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 // forward declarations to support abs or sqrt operations on VectorizedArray
41 template <typename Number>
42 class VectorizedArray;
43 template <typename T>
44 struct EnableIfScalar;
45 
46 DEAL_II_NAMESPACE_CLOSE
47 
48 // Declare / Import auto-differentiable math functions in(to) standard
49 // namespace before numbers::NumberTraits is defined
50 #ifdef DEAL_II_WITH_ADOLC
51 # include <deal.II/differentiation/ad/adolc_math.h>
52 
53 # include <adolc/adouble.h> // Taped double
54 #endif
55 // Ideally we'd like to #include <deal.II/differentiation/ad/sacado_math.h>
56 // but header indirectly references numbers.h. We therefore simply
57 // import the whole Sacado header at this point to get the math
58 // functions imported into the standard namespace.
59 #ifdef DEAL_II_TRILINOS_WITH_SACADO
60 # include <Sacado.hpp>
61 #endif
62 
63 namespace std
64 {
65  template <typename Number>
66  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number>
67  sqrt(const ::VectorizedArray<Number> &);
68  template <typename Number>
69  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number>
70  abs(const ::VectorizedArray<Number> &);
71  template <typename Number>
72  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number>
73  max(const ::VectorizedArray<Number> &,
74  const ::VectorizedArray<Number> &);
75  template <typename Number>
76  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number>
77  min(const ::VectorizedArray<Number> &,
78  const ::VectorizedArray<Number> &);
79  template <typename Number>
81  pow(const ::VectorizedArray<Number> &, const Number p);
82  template <typename Number>
84  sin(const ::VectorizedArray<Number> &);
85  template <typename Number>
87  cos(const ::VectorizedArray<Number> &);
88  template <typename Number>
90  tan(const ::VectorizedArray<Number> &);
91  template <typename Number>
93  exp(const ::VectorizedArray<Number> &);
94  template <typename Number>
96  log(const ::VectorizedArray<Number> &);
97 } // namespace std
98 
99 DEAL_II_NAMESPACE_OPEN
100 
116 namespace numbers
117 {
121  static constexpr double E = 2.7182818284590452354;
122 
126  static constexpr double LOG2E = 1.4426950408889634074;
127 
131  static constexpr double LOG10E = 0.43429448190325182765;
132 
136  static constexpr double LN2 = 0.69314718055994530942;
137 
141  static constexpr double LN10 = 2.30258509299404568402;
142 
146  static constexpr double PI = 3.14159265358979323846;
147 
151  static constexpr double PI_2 = 1.57079632679489661923;
152 
156  static constexpr double PI_4 = 0.78539816339744830962;
157 
161  static constexpr double SQRT2 = 1.41421356237309504880;
162 
166  static constexpr double SQRT1_2 = 0.70710678118654752440;
167 
173  template <typename Number, typename = void>
174  struct is_cuda_compatible : std::true_type
175  {};
176 
180  template <typename Number>
181  struct is_cuda_compatible<std::complex<Number>, void> : std::false_type
182  {};
183 
197  DEAL_II_DEPRECATED
198  bool
199  is_nan(const double x);
200 
210  bool
211  is_finite(const double x);
212 
217  bool
218  is_finite(const std::complex<double> &x);
219 
224  bool
225  is_finite(const std::complex<float> &x);
226 
235  bool
236  is_finite(const std::complex<long double> &x);
237 
248  template <typename Number1, typename Number2>
249  bool
250  values_are_equal(const Number1 &value_1, const Number2 &value_2);
251 
262  template <typename Number1, typename Number2>
263  bool
264  values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
265 
273  template <typename Number>
274  bool
275  value_is_zero(const Number &value);
276 
287  template <typename Number1, typename Number2>
288  bool
289  value_is_less_than(const Number1 &value_1, const Number2 &value_2);
290 
301  template <typename Number1, typename Number2>
302  bool
303  value_is_less_than_or_equal_to(const Number1 &value_1,
304  const Number2 &value_2);
305 
306 
307 
318  template <typename Number1, typename Number2>
319  bool
320  value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
321 
332  template <typename Number1, typename Number2>
333  bool
334  value_is_greater_than_or_equal_to(const Number1 &value_1,
335  const Number2 &value_2);
336 
347  template <typename number>
349  {
355  static constexpr bool is_complex = false;
356 
363  using real_type = number;
364 
372  static constexpr DEAL_II_CUDA_HOST_DEV const number &
373  conjugate(const number &x);
374 
383  template <typename Dummy = number>
384  static constexpr DEAL_II_CUDA_HOST_DEV
385  typename std::enable_if<std::is_same<Dummy, number>::value &&
387  real_type>::type
388  abs_square(const number &x);
389 
390  template <typename Dummy = number>
391  static constexpr
392  typename std::enable_if<std::is_same<Dummy, number>::value &&
394  real_type>::type
395  abs_square(const number &x);
396 
400  static real_type
401  abs(const number &x);
402  };
403 
404 
411  template <typename number>
412  struct NumberTraits<std::complex<number>>
413  {
419  static constexpr bool is_complex = true;
420 
427  using real_type = number;
428 
432  static constexpr std::complex<number>
433  conjugate(const std::complex<number> &x);
434 
441  static constexpr real_type
442  abs_square(const std::complex<number> &x);
443 
444 
448  static real_type
449  abs(const std::complex<number> &x);
450  };
451 
452  // --------------- inline and template functions ---------------- //
453 
454  inline bool
455  is_nan(const double x)
456  {
457  return std::isnan(x);
458  }
459 
460 
461 
462  inline bool
463  is_finite(const double x)
464  {
465  return std::isfinite(x);
466  }
467 
468 
469 
470  inline bool
471  is_finite(const std::complex<double> &x)
472  {
473  // Check complex numbers for infinity
474  // by testing real and imaginary part
475  return (is_finite(x.real()) && is_finite(x.imag()));
476  }
477 
478 
479 
480  inline bool
481  is_finite(const std::complex<float> &x)
482  {
483  // Check complex numbers for infinity
484  // by testing real and imaginary part
485  return (is_finite(x.real()) && is_finite(x.imag()));
486  }
487 
488 
489 
490  inline bool
491  is_finite(const std::complex<long double> &x)
492  {
493  // Same for std::complex<long double>
494  return (is_finite(x.real()) && is_finite(x.imag()));
495  }
496 
497 
498  template <typename number>
499  constexpr DEAL_II_CUDA_HOST_DEV const number &
501  {
502  return x;
503  }
504 
505 
506 
507  template <typename number>
508  template <typename Dummy>
509  constexpr DEAL_II_CUDA_HOST_DEV
510  typename std::enable_if<std::is_same<Dummy, number>::value &&
512  typename NumberTraits<number>::real_type>::type
514  {
515  return x * x;
516  }
517 
518 
519 
520  template <typename number>
521  template <typename Dummy>
522  constexpr
523  typename std::enable_if<std::is_same<Dummy, number>::value &&
525  typename NumberTraits<number>::real_type>::type
526  NumberTraits<number>::abs_square(const number &x)
527  {
528  return x * x;
529  }
530 
531 
532 
533  template <typename number>
535  NumberTraits<number>::abs(const number &x)
536  {
537  return std::abs(x);
538  }
539 
540 
541 
542  template <typename number>
543  constexpr std::complex<number>
544  NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
545  {
546  return std::conj(x);
547  }
548 
549 
550 
551  template <typename number>
552  typename NumberTraits<std::complex<number>>::real_type
553  NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
554  {
555  return std::abs(x);
556  }
557 
558 
559 
560  template <typename number>
561  constexpr typename NumberTraits<std::complex<number>>::real_type
562  NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
563  {
564  return std::norm(x);
565  }
566 
567 } // namespace numbers
568 
569 
570 // Forward declarations
572 {
573  namespace AD
574  {
575  namespace internal
576  {
577  // Defined in differentiation/ad/ad_number_traits.h
578  template <typename T>
579  struct NumberType;
580  } // namespace internal
581 
582  // Defined in differentiation/ad/ad_number_traits.h
583  template <typename NumberType>
584  struct is_ad_number;
585  } // namespace AD
586 } // namespace Differentiation
587 
588 
589 namespace internal
590 {
595  template <typename From, typename To>
597  {
598  // Source: https://stackoverflow.com/a/16944130
599  private:
600  template <typename T>
601  static void f(T);
602 
603  template <typename F, typename T>
604  static constexpr auto
605  test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
606  {
607  return true;
608  }
609 
610  template <typename F, typename T>
611  static constexpr auto
612  test(...) -> bool
613  {
614  return false;
615  }
616 
617  public:
618  static bool const value = test<From, To>(0);
619  };
620 
634  template <typename T>
635  struct NumberType
636  {
637  static constexpr DEAL_II_CUDA_HOST_DEV const T &
638  value(const T &t)
639  {
640  return t;
641  }
642 
643  // Below are generic functions that allows an overload for any
644  // type U that is transformable to type T. This is particularly
645  // useful when needing to cast exotic number types
646  // (e.g. auto-differentiable or symbolic numbers) to a floating
647  // point one, such as might happen when converting between tensor
648  // types.
649 
650  // Type T is constructible from F.
651  template <typename F>
652  static constexpr DEAL_II_CUDA_HOST_DEV T
653  value(const F &f,
654  typename std::enable_if<
655  !std::is_same<typename std::decay<T>::type,
656  typename std::decay<F>::type>::value &&
657  std::is_constructible<T, F>::value>::type * = nullptr)
658  {
659  return T(f);
660  }
661 
662  // Type T is explicitly convertible (but not constructible) from F.
663  template <typename F>
664  static constexpr T
665  value(const F &f,
666  typename std::enable_if<
667  !std::is_same<typename std::decay<T>::type,
668  typename std::decay<F>::type>::value &&
669  !std::is_constructible<T, F>::value &&
671  {
672  return static_cast<T>(f);
673  }
674 
675  // Sacado doesn't provide any conversion operators, so we have
676  // to extract the value and perform further conversions from there.
677  // To be safe, we extend this to other possible AD numbers that
678  // might fall into the same category.
679  template <typename F>
680  static T
681  value(const F &f,
682  typename std::enable_if<
683  !std::is_same<typename std::decay<T>::type,
684  typename std::decay<F>::type>::value &&
685  !std::is_constructible<T, F>::value &&
688  {
690  }
691  };
692 
693  template <typename T>
694  struct NumberType<std::complex<T>>
695  {
696  static constexpr const std::complex<T> &
697  value(const std::complex<T> &t)
698  {
699  return t;
700  }
701 
702  static constexpr std::complex<T>
703  value(const T &t)
704  {
705  return std::complex<T>(t);
706  }
707 
708  // Facilitate cast from complex<double> to complex<float>
709  template <typename U>
710  static constexpr std::complex<T>
711  value(const std::complex<U> &t)
712  {
713  return std::complex<T>(NumberType<T>::value(t.real()),
714  NumberType<T>::value(t.imag()));
715  }
716  };
717 
718 #ifdef DEAL_II_COMPILER_CUDA_AWARE
719  template <>
720  struct NumberType<cuComplex>
721  {
722  static cuComplex
723  value(const float t)
724  {
725  return make_cuComplex(t, 0.f);
726  }
727  };
728 
729  template <>
730  struct NumberType<cuDoubleComplex>
731  {
732  static cuDoubleComplex
733  value(const double t)
734  {
735  return make_cuDoubleComplex(t, 0.);
736  }
737  };
738 #endif
739 } // namespace internal
740 
741 namespace numbers
742 {
743 #ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
744 
755  // Defined in differentiation/ad/adolc_number_types.cc
756  bool
757  values_are_equal(const adouble &value_1, const adouble &value_2);
758 
759 
770  template <typename Number>
771  bool
772  values_are_equal(const adouble &value_1, const Number &value_2)
773  {
774  // Use the specialized definition for two ADOL-C taped types
775  return values_are_equal(value_1,
777  }
778 
779 
790  template <typename Number>
791  bool
792  values_are_equal(const Number &value_1, const adouble &value_2)
793  {
794  // Use the above definition
795  return values_are_equal(value_2, value_1);
796  }
797 
809  // Defined in differentiation/ad/adolc_number_types.cc
810  bool
811  value_is_less_than(const adouble &value_1, const adouble &value_2);
812 
813 
825  template <typename Number>
826  bool
827  value_is_less_than(const adouble &value_1, const Number &value_2)
828  {
829  // Use the specialized definition for two ADOL-C taped types
830  return value_is_less_than(value_1,
832  }
833 
834 
846  template <typename Number>
847  bool
848  value_is_less_than(const Number &value_1, const adouble &value_2)
849  {
850  // Use the specialized definition for two ADOL-C taped types
852  value_2);
853  }
854 
855 #endif
856 
857 
858  template <typename Number1, typename Number2>
859  inline bool
860  values_are_equal(const Number1 &value_1, const Number2 &value_2)
861  {
862  return (value_1 == internal::NumberType<Number1>::value(value_2));
863  }
864 
865 
866  template <typename Number1, typename Number2>
867  inline bool
868  values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
869  {
870  return !(values_are_equal(value_1, value_2));
871  }
872 
873 
874  template <typename Number>
875  inline bool
876  value_is_zero(const Number &value)
877  {
878  return values_are_equal(value, 0.0);
879  }
880 
881 
882  template <typename Number1, typename Number2>
883  inline bool
884  value_is_less_than(const Number1 &value_1, const Number2 &value_2)
885  {
886  return (value_1 < internal::NumberType<Number1>::value(value_2));
887  }
888 
889 
890  template <typename Number1, typename Number2>
891  inline bool
892  value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
893  {
894  return (value_is_less_than(value_1, value_2) ||
895  values_are_equal(value_1, value_2));
896  }
897 
898 
899  template <typename Number1, typename Number2>
900  bool
901  value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
902  {
903  return !(value_is_less_than_or_equal_to(value_1, value_2));
904  }
905 
906 
907  template <typename Number1, typename Number2>
908  inline bool
909  value_is_greater_than_or_equal_to(const Number1 &value_1,
910  const Number2 &value_2)
911  {
912  return !(value_is_less_than(value_1, value_2));
913  }
914 } // namespace numbers
915 
916 DEAL_II_NAMESPACE_CLOSE
917 
918 #endif
static constexpr double LOG2E
Definition: numbers.h:126
bool value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:892
static constexpr double LN10
Definition: numbers.h:141
bool value_is_zero(const Number &value)
Definition: numbers.h:876
static constexpr std::enable_if< std::is_same< Dummy, number >::value &&is_cuda_compatible< Dummy >::value, real_type >::type abs_square(const number &x)
Definition: numbers.h:513
STL namespace.
static real_type abs(const number &x)
Definition: numbers.h:535
static constexpr double LOG10E
Definition: numbers.h:131
bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:860
bool is_finite(const double x)
Definition: numbers.h:463
static constexpr double SQRT1_2
Definition: numbers.h:166
static constexpr double E
Definition: numbers.h:121
static constexpr double SQRT2
Definition: numbers.h:161
bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:868
static constexpr double PI_2
Definition: numbers.h:151
bool is_nan(const double x)
Definition: numbers.h:455
bool value_is_greater_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:909
static constexpr double LN2
Definition: numbers.h:136
bool value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:901
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:884
static constexpr double PI
Definition: numbers.h:146
T min(const T &t, const MPI_Comm &mpi_communicator)
static constexpr const number & conjugate(const number &x)
Definition: numbers.h:500
static constexpr double PI_4
Definition: numbers.h:156
static constexpr bool is_complex
Definition: numbers.h:355
T max(const T &t, const MPI_Comm &mpi_communicator)