numbers.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2006 - 2020 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_COMPILER_CUDA_AWARE
#  include <cuComplex.h>
#endif
 
#include <cmath>
#include <complex>
#include <cstddef>
#include <type_traits>
 
#ifdef DEAL_II_COMPILER_CUDA_AWARE
#  define DEAL_II_CUDA_HOST_DEV __host__ __device__
#else
#  define DEAL_II_CUDA_HOST_DEV
#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
 
DEAL_II_NAMESPACE_CLOSE
 
// Declare / Import auto-differentiable math functions in(to) standard
// namespace before numbers::NumberTraits is defined
#ifdef DEAL_II_WITH_ADOLC
#  include <deal.II/differentiation/ad/adolc_math.h>
 
#  include <adolc/adouble.h> // Taped double
#endif
// Ideally we'd like to #include <deal.II/differentiation/ad/sacado_math.h>
// but header indirectly references numbers.h. We therefore simply
// import the whole Sacado header at this point to get the math
// functions imported into the standard namespace.
#ifdef DEAL_II_TRILINOS_WITH_SACADO
#  include <Sacado.hpp>
#endif
 
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;
 
  template <typename Number, typename = void>
  struct is_cuda_compatible : std::true_type
  {};
 
  template <typename Number>
  struct is_cuda_compatible<std::complex<Number>, void> : std::false_type
  {};
 
  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 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 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;
 
    static constexpr DEAL_II_CUDA_HOST_DEV const number &
                                                 conjugate(const number &x);
 
    template <typename Dummy = number>
    static constexpr DEAL_II_CUDA_HOST_DEV
      typename std::enable_if<std::is_same<Dummy, number>::value &&
                                is_cuda_compatible<Dummy>::value,
                              real_type>::type
      abs_square(const number &x);
 
    template <typename Dummy = number>
    static constexpr
      typename std::enable_if<std::is_same<Dummy, number>::value &&
                                !is_cuda_compatible<Dummy>::value,
                              real_type>::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;
 
    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_CUDA_HOST_DEV const number &
                                        NumberTraits<number>::conjugate(const number &x)
  {
    return x;
  }
 
 
 
  template <typename number>
  template <typename Dummy>
  constexpr DEAL_II_CUDA_HOST_DEV
    typename std::enable_if<std::is_same<Dummy, number>::value &&
                              is_cuda_compatible<Dummy>::value,
                            typename NumberTraits<number>::real_type>::type
    NumberTraits<number>::abs_square(const number &x)
  {
    return x * x;
  }
 
 
 
  template <typename number>
  template <typename Dummy>
  constexpr
    typename std::enable_if<std::is_same<Dummy, number>::value &&
                              !is_cuda_compatible<Dummy>::value,
                            typename NumberTraits<number>::real_type>::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)
  {
    return std::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)
  {
    return std::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_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV 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_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV 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)
    {
      return T(f);
    }
 
    // Type T is explicitly convertible (but not constructible) from F.
    template <typename F>
    static constexpr DEAL_II_ALWAYS_INLINE 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)
    {
      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,
          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)
    {
      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_COMPILER_CUDA_AWARE
  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 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 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