Reference documentation for deal.II version 9.4.0
\(\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 - 2022 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
41namespace 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
123template <typename Number,
124 std::size_t width =
126class VectorizedArray;
127template <typename T>
128struct 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
146# include <Sacado.hpp>
148#endif
149
150namespace std
151{
152 template <typename Number, std::size_t width>
153 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
154 sqrt(const ::VectorizedArray<Number, width> &);
155 template <typename Number, std::size_t width>
156 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
157 abs(const ::VectorizedArray<Number, width> &);
158 template <typename Number, std::size_t width>
159 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
160 max(const ::VectorizedArray<Number, width> &,
161 const ::VectorizedArray<Number, width> &);
162 template <typename Number, std::size_t width>
163 DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
164 min(const ::VectorizedArray<Number, width> &,
165 const ::VectorizedArray<Number, width> &);
166 template <typename Number, size_t width>
168 pow(const ::VectorizedArray<Number, width> &, const Number p);
169 template <typename Number, size_t width>
171 sin(const ::VectorizedArray<Number, width> &);
172 template <typename Number, size_t width>
174 cos(const ::VectorizedArray<Number, width> &);
175 template <typename Number, size_t width>
177 tan(const ::VectorizedArray<Number, width> &);
178 template <typename Number, size_t width>
180 exp(const ::VectorizedArray<Number, width> &);
181 template <typename Number, size_t width>
183 log(const ::VectorizedArray<Number, width> &);
184} // namespace std
185
187
203namespace numbers
204{
208 static constexpr double E = 2.7182818284590452354;
209
213 static constexpr double LOG2E = 1.4426950408889634074;
214
218 static constexpr double LOG10E = 0.43429448190325182765;
219
223 static constexpr double LN2 = 0.69314718055994530942;
224
228 static constexpr double LN10 = 2.30258509299404568402;
229
233 static constexpr double PI = 3.14159265358979323846;
234
238 static constexpr double PI_2 = 1.57079632679489661923;
239
243 static constexpr double PI_4 = 0.78539816339744830962;
244
248 static constexpr double SQRT2 = 1.41421356237309504880;
249
253 static constexpr double SQRT1_2 = 0.70710678118654752440;
254
260 template <typename Number, typename = void>
261 struct is_cuda_compatible : std::true_type
262 {};
263
267 template <typename Number>
268 struct is_cuda_compatible<std::complex<Number>, void> : std::false_type
269 {};
270
280 bool
281 is_finite(const double x);
282
287 bool
288 is_finite(const std::complex<double> &x);
289
294 bool
295 is_finite(const std::complex<float> &x);
296
305 bool
306 is_finite(const std::complex<long double> &x);
307
318 template <typename Number1, typename Number2>
319 constexpr bool
320 values_are_equal(const Number1 &value_1, const Number2 &value_2);
321
332 template <typename Number1, typename Number2>
333 bool
334 values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
335
343 template <typename Number>
344 constexpr bool
345 value_is_zero(const Number &value);
346
357 template <typename Number1, typename Number2>
358 bool
359 value_is_less_than(const Number1 &value_1, const Number2 &value_2);
360
371 template <typename Number1, typename Number2>
372 bool
373 value_is_less_than_or_equal_to(const Number1 &value_1,
374 const Number2 &value_2);
375
376
377
388 template <typename Number1, typename Number2>
389 bool
390 value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
391
402 template <typename Number1, typename Number2>
403 bool
404 value_is_greater_than_or_equal_to(const Number1 &value_1,
405 const Number2 &value_2);
406
415 template <typename number>
417 {
423 static constexpr bool is_complex = false;
424
431 using real_type = number;
432
436 using double_type = double;
437
445 static constexpr DEAL_II_CUDA_HOST_DEV const number &
446 conjugate(const number &x);
447
456 template <typename Dummy = number>
457 static constexpr DEAL_II_CUDA_HOST_DEV
458 typename std::enable_if<std::is_same<Dummy, number>::value &&
460 real_type>::type
461 abs_square(const number &x);
462
463 template <typename Dummy = number>
464 static constexpr
465 typename std::enable_if<std::is_same<Dummy, number>::value &&
467 real_type>::type
468 abs_square(const number &x);
469
473 static real_type
474 abs(const number &x);
475 };
476
477
482 template <typename number>
483 struct NumberTraits<std::complex<number>>
484 {
490 static constexpr bool is_complex = true;
491
498 using real_type = number;
499
503 using double_type = std::complex<double>;
504
508 static constexpr std::complex<number>
509 conjugate(const std::complex<number> &x);
510
517 static constexpr real_type
518 abs_square(const std::complex<number> &x);
519
520
524 static real_type
525 abs(const std::complex<number> &x);
526 };
527
528 // --------------- inline and template functions ---------------- //
529
530 inline bool
531 is_nan(const double x)
532 {
533 return std::isnan(x);
534 }
535
536
537
538 inline bool
539 is_finite(const double x)
540 {
541 return std::isfinite(x);
542 }
543
544
545
546 inline bool
547 is_finite(const std::complex<double> &x)
548 {
549 // Check complex numbers for infinity
550 // by testing real and imaginary part
551 return (is_finite(x.real()) && is_finite(x.imag()));
552 }
553
554
555
556 inline bool
557 is_finite(const std::complex<float> &x)
558 {
559 // Check complex numbers for infinity
560 // by testing real and imaginary part
561 return (is_finite(x.real()) && is_finite(x.imag()));
562 }
563
564
565
566 inline bool
567 is_finite(const std::complex<long double> &x)
568 {
569 // Same for std::complex<long double>
570 return (is_finite(x.real()) && is_finite(x.imag()));
571 }
572
573
574 template <typename number>
575 constexpr DEAL_II_CUDA_HOST_DEV const number &
577 {
578 return x;
579 }
580
581
582
583 template <typename number>
584 template <typename Dummy>
585 constexpr DEAL_II_CUDA_HOST_DEV
586 typename std::enable_if<std::is_same<Dummy, number>::value &&
588 typename NumberTraits<number>::real_type>::type
590 {
591 return x * x;
592 }
593
594
595
596 template <typename number>
597 template <typename Dummy>
598 constexpr
599 typename std::enable_if<std::is_same<Dummy, number>::value &&
601 typename NumberTraits<number>::real_type>::type
602 NumberTraits<number>::abs_square(const number &x)
603 {
604 return x * x;
605 }
606
607
608
609 template <typename number>
612 {
613 return std::abs(x);
614 }
615
616
617
618 template <typename number>
619 constexpr std::complex<number>
620 NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
621 {
622 return std::conj(x);
623 }
624
625
626
627 template <typename number>
628 typename NumberTraits<std::complex<number>>::real_type
629 NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
630 {
631 return std::abs(x);
632 }
633
634
635
636 template <typename number>
637 constexpr typename NumberTraits<std::complex<number>>::real_type
638 NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
639 {
640 return std::norm(x);
641 }
642
643} // namespace numbers
644
645
646// Forward declarations
648{
649 namespace AD
650 {
651 namespace internal
652 {
653 // Defined in differentiation/ad/ad_number_traits.h
654 template <typename T>
656 } // namespace internal
657
658 // Defined in differentiation/ad/ad_number_traits.h
659 template <typename NumberType>
661 } // namespace AD
662} // namespace Differentiation
663
664
665namespace internal
666{
671 template <typename From, typename To>
673 {
674 // Source: https://stackoverflow.com/a/16944130
675 private:
676 template <typename T>
677 static void f(T);
678
679 template <typename F, typename T>
680 static constexpr auto
681 test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
682 {
683 return true;
684 }
685
686 template <typename F, typename T>
687 static constexpr auto
688 test(...) -> bool
689 {
690 return false;
691 }
692
693 public:
694 static bool const value = test<From, To>(0);
695 };
696
697 /*
698 * The structs below are needed to convert between some special number types.
699 * Also see tensor.h for another specialization.
700 */
701 template <typename T>
703 {
704 static constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV const T &
705 value(const T &t)
706 {
707 return t;
708 }
709
710 // Below are generic functions that allows an overload for any
711 // type U that is transformable to type T. This is particularly
712 // useful when needing to cast exotic number types
713 // (e.g. auto-differentiable or symbolic numbers) to a floating
714 // point one, such as might happen when converting between tensor
715 // types.
716
717 // Type T is constructible from F.
718 template <typename F>
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>::type * = nullptr)
725 {
726 return T(f);
727 }
728
729 // Type T is explicitly convertible (but not constructible) from F.
730 template <typename F>
731 static constexpr DEAL_II_ALWAYS_INLINE T
732 value(const F &f,
733 typename std::enable_if<
734 !std::is_same<typename std::decay<T>::type,
735 typename std::decay<F>::type>::value &&
736 !std::is_constructible<T, F>::value &&
738 {
739 return static_cast<T>(f);
740 }
741
742 // Sacado doesn't provide any conversion operators, so we have
743 // to extract the value and perform further conversions from there.
744 // To be safe, we extend this to other possible AD numbers that
745 // might fall into the same category.
746 template <typename F>
747 static 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 &&
755 {
757 }
758 };
759
760 template <typename T>
761 struct NumberType<std::complex<T>>
762 {
763 static constexpr const std::complex<T> &
764 value(const std::complex<T> &t)
765 {
766 return t;
767 }
768
769 static constexpr std::complex<T>
770 value(const T &t)
771 {
772 return std::complex<T>(t);
773 }
774
775 // Facilitate cast from complex<double> to complex<float>
776 template <typename U>
777 static constexpr std::complex<T>
778 value(const std::complex<U> &t)
779 {
780 return std::complex<T>(NumberType<T>::value(t.real()),
781 NumberType<T>::value(t.imag()));
782 }
783 };
784
785#ifdef DEAL_II_COMPILER_CUDA_AWARE
786 template <>
787 struct NumberType<cuComplex>
788 {
789 static cuComplex
790 value(const float t)
791 {
792 return make_cuComplex(t, 0.f);
793 }
794 };
795
796 template <>
797 struct NumberType<cuDoubleComplex>
798 {
799 static cuDoubleComplex
800 value(const double t)
801 {
802 return make_cuDoubleComplex(t, 0.);
803 }
804 };
805#endif
806} // namespace internal
807
808namespace numbers
809{
810#ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
811
822 // Defined in differentiation/ad/adolc_number_types.cc
823 bool
824 values_are_equal(const adouble &value_1, const adouble &value_2);
825
826
837 template <typename Number>
838 bool
839 values_are_equal(const adouble &value_1, const Number &value_2)
840 {
841 // Use the specialized definition for two ADOL-C taped types
842 return values_are_equal(value_1,
844 }
845
846
857 template <typename Number>
858 bool
859 values_are_equal(const Number &value_1, const adouble &value_2)
860 {
861 // Use the above definition
862 return values_are_equal(value_2, value_1);
863 }
864
876 // Defined in differentiation/ad/adolc_number_types.cc
877 bool
878 value_is_less_than(const adouble &value_1, const adouble &value_2);
879
880
892 template <typename Number>
893 bool
894 value_is_less_than(const adouble &value_1, const Number &value_2)
895 {
896 // Use the specialized definition for two ADOL-C taped types
897 return value_is_less_than(value_1,
899 }
900
901
913 template <typename Number>
914 bool
915 value_is_less_than(const Number &value_1, const adouble &value_2)
916 {
917 // Use the specialized definition for two ADOL-C taped types
919 value_2);
920 }
921
922#endif
923
924
925 template <typename Number1, typename Number2>
926 constexpr bool
927 values_are_equal(const Number1 &value_1, const Number2 &value_2)
928 {
929 return (value_1 == internal::NumberType<Number1>::value(value_2));
930 }
931
932
933 template <typename Number1, typename Number2>
934 inline bool
935 values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
936 {
937 return !(values_are_equal(value_1, value_2));
938 }
939
940
941 template <typename Number>
942 constexpr bool
943 value_is_zero(const Number &value)
944 {
945 return values_are_equal(value, 0.0);
946 }
947
948
949 template <typename Number1, typename Number2>
950 inline bool
951 value_is_less_than(const Number1 &value_1, const Number2 &value_2)
952 {
953 return (value_1 < internal::NumberType<Number1>::value(value_2));
954 }
955
956
957 template <typename Number1, typename Number2>
958 inline bool
959 value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
960 {
961 return (value_is_less_than(value_1, value_2) ||
962 values_are_equal(value_1, value_2));
963 }
964
965
966 template <typename Number1, typename Number2>
967 bool
968 value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
969 {
970 return !(value_is_less_than_or_equal_to(value_1, value_2));
971 }
972
973
974 template <typename Number1, typename Number2>
975 inline bool
976 value_is_greater_than_or_equal_to(const Number1 &value_1,
977 const Number2 &value_2)
978 {
979 return !(value_is_less_than(value_1, value_2));
980 }
981} // namespace numbers
982
984
985#endif
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:102
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_DISABLE_EXTRA_DIAGNOSTICS
Definition: config.h:456
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define DEAL_II_ENABLE_EXTRA_DIAGNOSTICS
Definition: config.h:495
static const char T
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition: divergence.h:472
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
static constexpr double LOG10E
Definition: numbers.h:218
static constexpr double PI_2
Definition: numbers.h:238
bool value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:959
static constexpr double E
Definition: numbers.h:208
static constexpr double PI
Definition: numbers.h:233
bool value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:968
static constexpr double SQRT2
Definition: numbers.h:248
constexpr bool value_is_zero(const Number &value)
Definition: numbers.h:943
bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:935
static constexpr double SQRT1_2
Definition: numbers.h:253
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:927
static constexpr double PI_4
Definition: numbers.h:243
static constexpr double LN10
Definition: numbers.h:228
static constexpr double LN2
Definition: numbers.h:223
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:951
bool is_finite(const double x)
Definition: numbers.h:539
static constexpr double LOG2E
Definition: numbers.h:213
bool value_is_greater_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition: numbers.h:976
bool is_nan(const double x)
Definition: numbers.h:531
STL namespace.
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > tan(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
#define DEAL_II_CUDA_HOST_DEV
Definition: numbers.h:34
static cuComplex value(const float t)
Definition: numbers.h:790
static cuDoubleComplex value(const double t)
Definition: numbers.h:800
static constexpr std::complex< T > value(const std::complex< U > &t)
Definition: numbers.h:778
static constexpr std::complex< T > value(const T &t)
Definition: numbers.h:770
static constexpr const std::complex< T > & value(const std::complex< T > &t)
Definition: numbers.h:764
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:720
static constexpr const T & value(const T &t)
Definition: numbers.h:705
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:748
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:732
static constexpr unsigned int max_width
Definition: numbers.h:65
static constexpr auto test(...) -> bool
Definition: numbers.h:688
static constexpr auto test(int) -> decltype(f(static_cast< T >(std::declval< F >())), true)
Definition: numbers.h:681
static constexpr std::enable_if< std::is_same< Dummy, number >::value &&!is_cuda_compatible< Dummy >::value, real_type >::type abs_square(const number &x)
static constexpr const number & conjugate(const number &x)
Definition: numbers.h:576
static constexpr bool is_complex
Definition: numbers.h:423
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:589
static real_type abs(const number &x)
Definition: numbers.h:611