numbers.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2006 - 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------
 
#ifndef dealii_numbers_h
#define dealii_numbers_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/types.h>
 
#ifdef DEAL_II_WITH_CUDA
#  include <cuComplex.h>
#endif
 
#include <Kokkos_Macros.hpp>
 
#include <cmath>
#include <complex>
#include <cstddef>
#include <type_traits>
 
#define DEAL_II_HOST_DEVICE KOKKOS_FUNCTION
#define DEAL_II_CUDA_HOST_DEV DEAL_II_HOST_DEVICE
#define DEAL_II_HOST_DEVICE_ALWAYS_INLINE KOKKOS_FORCEINLINE_FUNCTION
 
// clang++ assumes that all constexpr functions are __host__ __device__ when
// compiling CUDA code, i.e, when Kokkos was configured with CUDA support.
// This is problematic when calling non-constexpr functions in constexpr
// functions. Hence, we need a way to annotate functions explicitly as
// host-only.
#if defined(__clang__) && defined(__CUDA__)
#  define DEAL_II_HOST __host__
#else
#  define DEAL_II_HOST
#endif
 
// Forward-declare the automatic differentiation types so we can add prototypes
// for our own wrappers.
#ifdef DEAL_II_WITH_ADOLC
class adouble;
namespace adtl
{
  class adouble;
}
#endif
 
DEAL_II_NAMESPACE_OPEN
 
namespace internal
{
  template <typename Number>
  struct VectorizedArrayWidthSpecifier
  {
    constexpr static unsigned int max_width = 1;
  };
 
  template <>
  struct VectorizedArrayWidthSpecifier<double>
  {
    constexpr static unsigned int max_width =
#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512
      8;
#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256
      4;
#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128
      2;
#else
      1;
#endif
  };
 
  template <>
  struct VectorizedArrayWidthSpecifier<float>
  {
    constexpr static unsigned int max_width =
#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ALTIVEC__)
      4;
#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && defined(__AVX512F__)
      16;
#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)
      8;
#elif DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__)
      4;
#else
      1;
#endif
  };
 
 
} // namespace internal
 
// forward declarations to support abs or sqrt operations on VectorizedArray
#ifndef DOXYGEN
template <typename Number,
          std::size_t width =
            internal::VectorizedArrayWidthSpecifier<Number>::max_width>
class VectorizedArray;
template <typename T>
struct EnableIfScalar;
#endif
 
#ifdef DEAL_II_WITH_ADOLC
#  ifndef DOXYGEN
// Prototype some inline functions present in adolc_math.h for use in
// NumberTraits.
//
// ADOL-C uses fabs(), but for genericity we want to use abs(). Simultaneously,
// though, we don't want to include ADOL-C headers in this header since
// numbers.h is in everything. To get around this: use C++ rules which permit
// the use of forward-declared classes in function prototypes to declare some
// functions which are defined in adolc_math.h. This permits us to write "using
// ::abs;" in NumberTraits which will allow us to select the correct
// overload (the one in ::) when instantiating NumberTraits for ADOL-C
// types.
 
adouble
abs(const adouble &x);
 
adtl::adouble
abs(const adtl::adouble &x);
#  endif
#endif
 
DEAL_II_NAMESPACE_CLOSE
 
namespace std
{
  template <typename Number, std::size_t width>
  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
  sqrt(const ::VectorizedArray<Number, width> &);
  template <typename Number, std::size_t width>
  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
  abs(const ::VectorizedArray<Number, width> &);
  template <typename Number, std::size_t width>
  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
  max(const ::VectorizedArray<Number, width> &,
      const ::VectorizedArray<Number, width> &);
  template <typename Number, std::size_t width>
  DEAL_II_ALWAYS_INLINE ::VectorizedArray<Number, width>
  min(const ::VectorizedArray<Number, width> &,
      const ::VectorizedArray<Number, width> &);
  template <typename Number, size_t width>
  ::VectorizedArray<Number, width>
  pow(const ::VectorizedArray<Number, width> &, const Number p);
  template <typename Number, size_t width>
  ::VectorizedArray<Number, width>
  sin(const ::VectorizedArray<Number, width> &);
  template <typename Number, size_t width>
  ::VectorizedArray<Number, width>
  cos(const ::VectorizedArray<Number, width> &);
  template <typename Number, size_t width>
  ::VectorizedArray<Number, width>
  tan(const ::VectorizedArray<Number, width> &);
  template <typename Number, size_t width>
  ::VectorizedArray<Number, width>
  exp(const ::VectorizedArray<Number, width> &);
  template <typename Number, size_t width>
  ::VectorizedArray<Number, width>
  log(const ::VectorizedArray<Number, width> &);
} // namespace std
 
DEAL_II_NAMESPACE_OPEN
 
namespace numbers
{
  static constexpr double E = 2.7182818284590452354;
 
  static constexpr double LOG2E = 1.4426950408889634074;
 
  static constexpr double LOG10E = 0.43429448190325182765;
 
  static constexpr double LN2 = 0.69314718055994530942;
 
  static constexpr double LN10 = 2.30258509299404568402;
 
  static constexpr double PI = 3.14159265358979323846;
 
  static constexpr double PI_2 = 1.57079632679489661923;
 
  static constexpr double PI_4 = 0.78539816339744830962;
 
  static constexpr double SQRT2 = 1.41421356237309504880;
 
  static constexpr double SQRT1_2 = 0.70710678118654752440;
 
  bool
  is_finite(const double x);
 
  bool
  is_finite(const std::complex<double> &x);
 
  bool
  is_finite(const std::complex<float> &x);
 
  bool
  is_finite(const std::complex<long double> &x);
 
  template <typename Number1, typename Number2>
  constexpr DEAL_II_HOST_DEVICE bool
  values_are_equal(const Number1 &value_1, const Number2 &value_2);
 
  template <typename Number1, typename Number2>
  bool
  values_are_not_equal(const Number1 &value_1, const Number2 &value_2);
 
  template <typename Number>
  constexpr DEAL_II_HOST_DEVICE bool
  value_is_zero(const Number &value);
 
  template <typename Number1, typename Number2>
  bool
  value_is_less_than(const Number1 &value_1, const Number2 &value_2);
 
  template <typename Number1, typename Number2>
  bool
  value_is_less_than_or_equal_to(const Number1 &value_1,
                                 const Number2 &value_2);
 
 
 
  template <typename Number1, typename Number2>
  bool
  value_is_greater_than(const Number1 &value_1, const Number2 &value_2);
 
  template <typename Number1, typename Number2>
  bool
  value_is_greater_than_or_equal_to(const Number1 &value_1,
                                    const Number2 &value_2);
 
  template <typename number>
  struct NumberTraits
  {
    static constexpr bool is_complex = false;
 
    using real_type = number;
 
    using double_type = double;
 
    static constexpr DEAL_II_HOST_DEVICE const number &
                                               conjugate(const number &x);
 
    static constexpr DEAL_II_HOST_DEVICE real_type
    abs_square(const number &x);
 
    static real_type
    abs(const number &x);
  };
 
 
  template <typename number>
  struct NumberTraits<std::complex<number>>
  {
    static constexpr bool is_complex = true;
 
    using real_type = number;
 
    using double_type = std::complex<double>;
 
    static constexpr std::complex<number>
    conjugate(const std::complex<number> &x);
 
    static constexpr real_type
    abs_square(const std::complex<number> &x);
 
 
    static real_type
    abs(const std::complex<number> &x);
  };
 
  // --------------- inline and template functions ---------------- //
 
  inline bool
  is_nan(const double x)
  {
    return std::isnan(x);
  }
 
 
 
  inline bool
  is_finite(const double x)
  {
    return std::isfinite(x);
  }
 
 
 
  inline bool
  is_finite(const std::complex<double> &x)
  {
    // Check complex numbers for infinity
    // by testing real and imaginary part
    return (is_finite(x.real()) && is_finite(x.imag()));
  }
 
 
 
  inline bool
  is_finite(const std::complex<float> &x)
  {
    // Check complex numbers for infinity
    // by testing real and imaginary part
    return (is_finite(x.real()) && is_finite(x.imag()));
  }
 
 
 
  inline bool
  is_finite(const std::complex<long double> &x)
  {
    // Same for std::complex<long double>
    return (is_finite(x.real()) && is_finite(x.imag()));
  }
 
 
  template <typename number>
  constexpr DEAL_II_HOST_DEVICE const number &
  NumberTraits<number>::conjugate(const number &x)
  {
    return x;
  }
 
 
 
  template <typename number>
  constexpr DEAL_II_HOST_DEVICE typename NumberTraits<number>::real_type
  NumberTraits<number>::abs_square(const number &x)
  {
    return x * x;
  }
 
 
 
  template <typename number>
  typename NumberTraits<number>::real_type
  NumberTraits<number>::abs(const number &x)
  {
    // Make things work with AD types
    using std::abs;
#ifdef DEAL_II_WITH_ADOLC
    // This one is a little tricky - we have our own abs function in ::,
    // prototyped with forward-declared types in this file, but it only exists
    // if we have ADOL-C: hence we only add this using statement in that
    // situation
    using ::abs;
#endif
    return abs(x);
  }
 
 
 
  template <typename number>
  constexpr std::complex<number>
  NumberTraits<std::complex<number>>::conjugate(const std::complex<number> &x)
  {
    return std::conj(x);
  }
 
 
 
  template <typename number>
  typename NumberTraits<std::complex<number>>::real_type
  NumberTraits<std::complex<number>>::abs(const std::complex<number> &x)
  {
    // Make things work with AD types
    using std::abs;
#ifdef DEAL_II_WITH_ADOLC
    // Same comment as the non-complex case holds here
    using ::abs;
#endif
    return abs(x);
  }
 
 
 
  template <typename number>
  constexpr typename NumberTraits<std::complex<number>>::real_type
  NumberTraits<std::complex<number>>::abs_square(const std::complex<number> &x)
  {
    return std::norm(x);
  }
 
} // namespace numbers
 
 
// Forward declarations
namespace Differentiation
{
  namespace AD
  {
    namespace internal
    {
      // Defined in differentiation/ad/ad_number_traits.h
      template <typename T>
      struct NumberType;
    } // namespace internal
 
    // Defined in differentiation/ad/ad_number_traits.h
    template <typename NumberType>
    struct is_ad_number;
  } // namespace AD
} // namespace Differentiation
 
 
namespace internal
{
  template <typename From, typename To>
  struct is_explicitly_convertible
  {
    // Source: https://stackoverflow.com/a/16944130
  private:
    template <typename T>
    static void f(T);
 
    template <typename F, typename T>
    static constexpr auto
    test(int) -> decltype(f(static_cast<T>(std::declval<F>())), true)
    {
      return true;
    }
 
    template <typename F, typename T>
    static constexpr auto
    test(...) -> bool
    {
      return false;
    }
 
  public:
    static bool const value = test<From, To>(0);
  };
 
  /*
   * The structs below are needed to convert between some special number types.
   * Also see tensor.h for another specialization.
   */
  template <typename T>
  struct NumberType
  {
    static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE const T &
                                                             value(const T &t)
    {
      return t;
    }
 
    // Below are generic functions that allows an overload for any
    // type U that is transformable to type T. This is particularly
    // useful when needing to cast exotic number types
    // (e.g. auto-differentiable or symbolic numbers) to a floating
    // point one, such as might  happen when converting between tensor
    // types.
 
    // Type T is constructible from F.
    template <typename F>
    static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T
    value(const F &f,
          std::enable_if_t<!std::is_same<typename std::decay<T>::type,
                                         typename std::decay<F>::type>::value &&
                           std::is_constructible<T, F>::value> * = nullptr)
    {
      return T(f);
    }
 
    // Type T is explicitly convertible (but not constructible) from F.
    template <typename F>
    static constexpr DEAL_II_HOST_DEVICE_ALWAYS_INLINE T
    value(const F &f,
          std::enable_if_t<!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> * =
            nullptr)
    {
      return static_cast<T>(f);
    }
 
    // Sacado doesn't provide any conversion operators, so we have
    // to extract the value and perform further conversions from there.
    // To be safe, we extend this to other possible AD numbers that
    // might fall into the same category.
    template <typename F>
    static T
    value(
      const F &f,
      std::enable_if_t<!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> * = nullptr)
    {
      return Differentiation::AD::internal::NumberType<T>::value(f);
    }
  };
 
  template <typename T>
  struct NumberType<std::complex<T>>
  {
    static constexpr const std::complex<T> &
    value(const std::complex<T> &t)
    {
      return t;
    }
 
    static constexpr std::complex<T>
    value(const T &t)
    {
      return std::complex<T>(t);
    }
 
    // Facilitate cast from complex<double> to complex<float>
    template <typename U>
    static constexpr std::complex<T>
    value(const std::complex<U> &t)
    {
      return std::complex<T>(NumberType<T>::value(t.real()),
                             NumberType<T>::value(t.imag()));
    }
  };
 
#ifdef DEAL_II_WITH_CUDA
  template <>
  struct NumberType<cuComplex>
  {
    static cuComplex
    value(const float t)
    {
      return make_cuComplex(t, 0.f);
    }
  };
 
  template <>
  struct NumberType<cuDoubleComplex>
  {
    static cuDoubleComplex
    value(const double t)
    {
      return make_cuDoubleComplex(t, 0.);
    }
  };
#endif
} // namespace internal
 
namespace numbers
{
#ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
 
  // Defined in differentiation/ad/adolc_number_types.cc
  bool
  values_are_equal(const adouble &value_1, const adouble &value_2);
 
 
  template <typename Number>
  bool
  values_are_equal(const adouble &value_1, const Number &value_2)
  {
    // Use the specialized definition for two ADOL-C taped types
    return values_are_equal(
      value_1, ::internal::NumberType<adouble>::value(value_2));
  }
 
 
  template <typename Number>
  bool
  values_are_equal(const Number &value_1, const adouble &value_2)
  {
    // Use the above definition
    return values_are_equal(value_2, value_1);
  }
 
  // Defined in differentiation/ad/adolc_number_types.cc
  bool
  value_is_less_than(const adouble &value_1, const adouble &value_2);
 
 
  template <typename Number>
  bool
  value_is_less_than(const adouble &value_1, const Number &value_2)
  {
    // Use the specialized definition for two ADOL-C taped types
    return value_is_less_than(
      value_1, ::internal::NumberType<adouble>::value(value_2));
  }
 
 
  template <typename Number>
  bool
  value_is_less_than(const Number &value_1, const adouble &value_2)
  {
    // Use the specialized definition for two ADOL-C taped types
    return value_is_less_than(
      ::internal::NumberType<adouble>::value(value_1), value_2);
  }
 
#endif
 
 
  template <typename Number1, typename Number2>
  constexpr DEAL_II_HOST_DEVICE bool
  values_are_equal(const Number1 &value_1, const Number2 &value_2)
  {
    return (value_1 == ::internal::NumberType<Number1>::value(value_2));
  }
 
 
  template <typename Number1, typename Number2>
  inline bool
  values_are_not_equal(const Number1 &value_1, const Number2 &value_2)
  {
    return !(values_are_equal(value_1, value_2));
  }
 
 
  template <typename Number>
  constexpr DEAL_II_HOST_DEVICE bool
  value_is_zero(const Number &value)
  {
    return values_are_equal(value, 0.0);
  }
 
 
  template <typename Number1, typename Number2>
  inline bool
  value_is_less_than(const Number1 &value_1, const Number2 &value_2)
  {
    return (value_1 < ::internal::NumberType<Number1>::value(value_2));
  }
 
 
  template <typename Number1, typename Number2>
  inline bool
  value_is_less_than_or_equal_to(const Number1 &value_1, const Number2 &value_2)
  {
    return (value_is_less_than(value_1, value_2) ||
            values_are_equal(value_1, value_2));
  }
 
 
  template <typename Number1, typename Number2>
  bool
  value_is_greater_than(const Number1 &value_1, const Number2 &value_2)
  {
    return !(value_is_less_than_or_equal_to(value_1, value_2));
  }
 
 
  template <typename Number1, typename Number2>
  inline bool
  value_is_greater_than_or_equal_to(const Number1 &value_1,
                                    const Number2 &value_2)
  {
    return !(value_is_less_than(value_1, value_2));
  }
} // namespace numbers
 
DEAL_II_NAMESPACE_CLOSE
 
#endif