deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+00:00
\(\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\}}\)
Loading...
Searching...
No Matches
numbers.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2006 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_numbers_h
16#define dealii_numbers_h
17
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/base/types.h>
22
23#include <Kokkos_Macros.hpp>
24
25#include <cmath>
26#include <complex>
27#include <cstddef>
28#include <type_traits>
29
30// Forward-declare the automatic differentiation types so we can add prototypes
31// for our own wrappers.
32#ifdef DEAL_II_WITH_ADOLC
33class adouble;
34namespace adtl
35{
36 class adouble;
37}
38#endif
39
41
42namespace internal
43{
60 template <typename Number>
62 {
66 constexpr static unsigned int max_width = 1;
67 };
68
75 template <>
77 {
81 constexpr static unsigned int max_width =
82#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512
83 8;
84#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256
85 4;
86#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128
87 2;
88#else
89 1;
90#endif
91 };
92
99 template <>
101 {
105 constexpr static unsigned int max_width =
106#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ALTIVEC__)
107 4;
108#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && defined(__AVX512F__)
109 16;
110#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)
111 8;
112#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__)
113 4;
114#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ARM_NEON)
115 4;
116#else
117 1;
118#endif
119 };
120
121
122} // namespace internal
123
124// forward declarations to support abs or sqrt operations on VectorizedArray
125#ifndef DOXYGEN
126template <typename Number,
127 std::size_t width =
129class VectorizedArray;
130template <typename T>
131struct EnableIfScalar;
132#endif
133
134#ifdef DEAL_II_WITH_ADOLC
135# ifndef DOXYGEN
136// Prototype some inline functions present in adolc_math.h for use in
137// NumberTraits.
138//
139// ADOL-C uses fabs(), but for genericity we want to use abs(). Simultaneously,
140// though, we don't want to include ADOL-C headers in this header since
141// numbers.h is in everything. To get around this: use C++ rules which permit
142// the use of forward-declared classes in function prototypes to declare some
143// functions which are defined in adolc_math.h. This permits us to write "using
144// ::abs;" in NumberTraits which will allow us to select the correct
145// overload (the one in ::) when instantiating NumberTraits for ADOL-C
146// types.
147
149abs(const adouble &x);
150
151adtl::adouble
152abs(const adtl::adouble &x);
153# endif
154#endif
155
157
158namespace std
159{
160 template <typename Number, std::size_t width>
163 template <typename Number, std::size_t width>
166 template <typename Number, std::size_t width>
170 template <typename Number, std::size_t width>
174 template <typename Number, size_t width>
177 template <typename Number, size_t width>
180 template <typename Number, size_t width>
183 template <typename Number, size_t width>
186 template <typename Number, size_t width>
189 template <typename Number, size_t width>
192} // namespace std
193
195
211namespace numbers
212{
216 constexpr double E = 2.7182818284590452354;
217
221 constexpr double LOG2E = 1.4426950408889634074;
222
226 constexpr double LOG10E = 0.43429448190325182765;
227
231 constexpr double LN2 = 0.69314718055994530942;
232
236 constexpr double LN10 = 2.30258509299404568402;
237
241 constexpr double PI = 3.14159265358979323846;
242
246 constexpr double PI_2 = 1.57079632679489661923;
247
251 constexpr double PI_4 = 0.78539816339744830962;
252
256 constexpr double SQRT2 = 1.41421356237309504880;
257
261 constexpr double SQRT1_2 = 0.70710678118654752440;
262
272 bool
273 is_finite(const double x);
274
279 bool
280 is_finite(const std::complex<double> &x);
281
286 bool
287 is_finite(const std::complex<float> &x);
288
297 bool
298 is_finite(const std::complex<long double> &x);
299
310 template <typename Number1, typename Number2>
311 constexpr DEAL_II_HOST_DEVICE bool
312 values_are_equal(const Number1 &value_1, const Number2 &value_2);
313
324 template <typename Number1, typename Number2>
325 constexpr bool
326 values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
327
335 template <typename Number>
336 constexpr DEAL_II_HOST_DEVICE bool
337 value_is_zero(const Number &value);
338
349 template <typename Number1, typename Number2>
350 bool
351 value_is_less_than(const Number1 &value_1, const Number2 &value_2);
352
363 template <typename Number1, typename Number2>
364 bool
366 const Number2 &value_2);
367
368
369
380 template <typename Number1, typename Number2>
381 bool
382 value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
383
394 template <typename Number1, typename Number2>
395 bool
397 const Number2 &value_2);
398
407 template <typename number>
409 {
415 static constexpr bool is_complex = false;
416
423 using real_type = number;
424
428 using double_type = double;
429
437 static constexpr DEAL_II_HOST_DEVICE const number &
438 conjugate(const number &x);
439
448 static constexpr DEAL_II_HOST_DEVICE real_type
449 abs_square(const number &x);
450
454 static real_type
455 abs(const number &x);
456 };
457
458
463 template <typename number>
464 struct NumberTraits<std::complex<number>>
465 {
471 static constexpr bool is_complex = true;
472
479 using real_type = number;
480
484 using double_type = std::complex<double>;
485
489 static constexpr std::complex<number>
490 conjugate(const std::complex<number> &x);
491
498 static constexpr real_type
499 abs_square(const std::complex<number> &x);
500
501
505 static real_type
506 abs(const std::complex<number> &x);
507 };
508
509 // --------------- inline and template functions ---------------- //
510
511 inline bool
512 is_nan(const double x)
513 {
514 return std::isnan(x);
515 }
516
517
518
519 inline bool
520 is_finite(const double x)
521 {
522 return std::isfinite(x);
523 }
524
525
526
527 inline bool
528 is_finite(const std::complex<double> &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<float> &x)
539 {
540 // Check complex numbers for infinity
541 // by testing real and imaginary part
542 return (is_finite(x.real()) && is_finite(x.imag()));
543 }
544
545
546
547 inline bool
548 is_finite(const std::complex<long double> &x)
549 {
550 // Same for std::complex<long double>
551 return (is_finite(x.real()) && is_finite(x.imag()));
552 }
553
554
555 template <typename number>
556 constexpr DEAL_II_HOST_DEVICE const number &
558 {
559 return x;
560 }
561
562
563
564 template <typename number>
567 {
568 return x * x;
569 }
570
571
572
573 template <typename number>
576 {
577 // Make things work with AD types
578 using std::abs;
579#ifdef DEAL_II_WITH_ADOLC
580 // This one is a little tricky - we have our own abs function in ::,
581 // prototyped with forward-declared types in this file, but it only exists
582 // if we have ADOL-C: hence we only add this using statement in that
583 // situation
584 using ::abs;
585#endif
586 return abs(x);
587 }
588
589
590
591 template <typename number>
592 constexpr std::complex<number>
593 NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
594 {
595 return std::conj(x);
596 }
597
598
599
600 template <typename number>
601 typename NumberTraits<std::complex<number>>::real_type
602 NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
603 {
604 // Make things work with AD types
605 using std::abs;
606#ifdef DEAL_II_WITH_ADOLC
607 // Same comment as the non-complex case holds here
608 using ::abs;
609#endif
610 return abs(x);
611 }
612
613
614
615 template <typename number>
616 constexpr typename NumberTraits<std::complex<number>>::real_type
617 NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
618 {
619 return std::norm(x);
620 }
621
622} // namespace numbers
623
624
625// Forward declarations
627{
628 namespace AD
629 {
630 namespace internal
631 {
632 // Defined in differentiation/ad/ad_number_traits.h
633 template <typename T>
635 } // namespace internal
636
637 // Defined in differentiation/ad/ad_number_traits.h
638 template <typename NumberType>
640 } // namespace AD
641} // namespace Differentiation
642
643
644namespace internal
645{
650 template <typename From, typename To>
652 {
653 // Source: https://stackoverflow.com/a/16944130
654 private:
655 template <typename T>
656 static void f(T);
657
658 template <typename F, typename T>
659 static constexpr auto
660 test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
661 {
662 return true;
663 }
664
665 template <typename F, typename T>
666 static constexpr auto
667 test(...) -> bool
668 {
669 return false;
670 }
671
672 public:
673 static const bool value = test<From, To>(0);
674 };
675
676 /*
677 * The structs below are needed to convert between some special number types.
678 * Also see tensor.h for another specialization.
679 */
680 template <typename T>
682 {
683 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T &
684 value(const T &t)
685 {
686 return t;
687 }
688
689 // Below are generic functions that allows an overload for any
690 // type U that is transformable to type T. This is particularly
691 // useful when needing to cast exotic number types
692 // (e.g. auto-differentiable or symbolic numbers) to a floating
693 // point one, such as might happen when converting between tensor
694 // types.
695
696 // Type T is constructible from F.
697 template <typename F>
698 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T
699 value(const F &f,
700 std::enable_if_t<!std::is_same_v<std::decay_t<T>, std::decay_t<F>> &&
701 std::is_constructible_v<T, F>> * = nullptr)
702 {
703 return T(f);
704 }
705
706 // Type T is explicitly convertible (but not constructible) from F.
707 template <typename F>
708 static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T
709 value(const F &f,
710 std::enable_if_t<!std::is_same_v<std::decay_t<T>, std::decay_t<F>> &&
711 !std::is_constructible_v<T, F> &&
713 nullptr)
714 {
715 return static_cast<T>(f);
716 }
717
718 // Sacado doesn't provide any conversion operators, so we have
719 // to extract the value and perform further conversions from there.
720 // To be safe, we extend this to other possible AD numbers that
721 // might fall into the same category.
722 template <typename F>
723 static T
725 const F &f,
726 std::enable_if_t<!std::is_same_v<std::decay_t<T>, std::decay_t<F>> &&
727 !std::is_constructible_v<T, F> &&
730 {
732 }
733 };
734
735 template <typename T>
736 struct NumberType<std::complex<T>>
737 {
738 static constexpr const std::complex<T> &
739 value(const std::complex<T> &t)
740 {
741 return t;
742 }
743
744 static constexpr std::complex<T>
745 value(const T &t)
746 {
747 return std::complex<T>(t);
748 }
749
750 // Facilitate cast from complex<double> to complex<float>
751 template <typename U>
752 static constexpr std::complex<T>
753 value(const std::complex<U> &t)
754 {
755 return std::complex<T>(NumberType<T>::value(t.real()),
756 NumberType<T>::value(t.imag()));
757 }
758 };
759
760} // namespace internal
761
762namespace numbers
763{
764#ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
765
776 // Defined in differentiation/ad/adolc_number_types.cc
777 bool
779
780
791 template <typename Number>
792 bool
794 {
795 // Use the specialized definition for two ADOL-C taped types
796 return values_are_equal(
798 }
799
800
811 template <typename Number>
812 bool
814 {
815 // Use the above definition
817 }
818
830 // Defined in differentiation/ad/adolc_number_types.cc
831 bool
833
834
846 template <typename Number>
847 bool
849 {
850 // Use the specialized definition for two ADOL-C taped types
851 return value_is_less_than(
853 }
854
855
867 template <typename Number>
868 bool
870 {
871 // Use the specialized definition for two ADOL-C taped types
872 return value_is_less_than(
874 }
875
876#endif
877
878
879 template <typename Number1, typename Number2>
880 constexpr DEAL_II_HOST_DEVICE bool
881 values_are_equal(const Number1 &value_1, const Number2 &value_2)
882 {
884 }
885
886
887 template <typename Number1, typename Number2>
888 inline constexpr bool
890 {
891 return !(values_are_equal(value_1, value_2));
892 }
893
894
895 template <typename Number>
896 constexpr DEAL_II_HOST_DEVICE bool
897 value_is_zero(const Number &value)
898 {
899 return values_are_equal(value, 0.0);
900 }
901
902
903 template <typename Number1, typename Number2>
904 inline bool
909
910
911 template <typename Number1, typename Number2>
912 inline bool
918
919
920 template <typename Number1, typename Number2>
921 bool
923 {
925 }
926
927
928 template <typename Number1, typename Number2>
929 inline bool
931 const Number2 &value_2)
932 {
934 }
935} // namespace numbers
936
938
939#endif
friend class Tensor
Definition tensor.h:865
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_HOST_DEVICE
Definition config.h:114
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
#define DEAL_II_HOST_DEVICE_ALWAYS_INLINE
Definition config.h:115
Definition numbers.h:35
constexpr double LOG10E
Definition numbers.h:226
constexpr double PI_2
Definition numbers.h:246
bool value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:913
constexpr double E
Definition numbers.h:216
constexpr double PI
Definition numbers.h:241
bool value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:922
constexpr bool values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:889
constexpr double SQRT2
Definition numbers.h:256
constexpr bool value_is_zero(const Number &value)
Definition numbers.h:897
constexpr double SQRT1_2
Definition numbers.h:261
constexpr bool values_are_equal(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:881
constexpr double PI_4
Definition numbers.h:251
constexpr double LN10
Definition numbers.h:236
constexpr double LN2
Definition numbers.h:231
bool value_is_less_than(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:905
bool is_finite(const double x)
Definition numbers.h:520
constexpr double LOG2E
Definition numbers.h:221
bool value_is_greater_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
Definition numbers.h:930
bool is_nan(const double x)
Definition numbers.h:512
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 > &)
static constexpr std::complex< T > value(const std::complex< U > &t)
Definition numbers.h:753
static constexpr std::complex< T > value(const T &t)
Definition numbers.h:745
static constexpr const std::complex< T > & value(const std::complex< T > &t)
Definition numbers.h:739
static constexpr T value(const F &f, std::enable_if_t<!std::is_same_v< std::decay_t< T >, std::decay_t< F > > &&std::is_constructible_v< T, F > > *=nullptr)
Definition numbers.h:699
static constexpr const T & value(const T &t)
Definition numbers.h:684
static constexpr T value(const F &f, std::enable_if_t<!std::is_same_v< std::decay_t< T >, std::decay_t< F > > &&!std::is_constructible_v< T, F > &&is_explicitly_convertible< const F, T >::value > *=nullptr)
Definition numbers.h:709
static T value(const F &f, std::enable_if_t<!std::is_same_v< std::decay_t< T >, std::decay_t< F > > &&!std::is_constructible_v< T, F > &&!is_explicitly_convertible< const F, T >::value &&Differentiation::AD::is_ad_number< F >::value > *=nullptr)
Definition numbers.h:724
static constexpr unsigned int max_width
Definition numbers.h:66
static constexpr auto test(...) -> bool
Definition numbers.h:667
static constexpr auto test(int) -> decltype(f(static_cast< T >(std::declval< F >())), true)
Definition numbers.h:660
static constexpr const number & conjugate(const number &x)
Definition numbers.h:557
static constexpr bool is_complex
Definition numbers.h:415
static real_type abs(const number &x)
Definition numbers.h:575
static constexpr real_type abs_square(const number &x)
Definition numbers.h:566