simd.h

Go to the documentation of this file.

//
// SPDX-License-Identifier: MIT
// Copyright (C) 2020 - 2023 by the ryujin authors
//
 
#pragma once
 
#include <compile_time_options.h>
 
#include <deal.II/base/tensor.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/vectorization.h>
 
namespace ryujin
{
 
  template <typename T>
  struct get_value_type {
    using type = T;
  };
 
 
  template <typename T, std::size_t width>
  struct get_value_type<dealii::VectorizedArray<T, width>> {
    using type = T;
  };
 
 
  template <typename T>
  unsigned int get_stride_size = 1;
 
  template <typename T, std::size_t width>
  unsigned int get_stride_size<dealii::VectorizedArray<T, width>> = width;
 
 
#ifndef DOXYGEN
  namespace
  {
    template <typename Functor, size_t... Is>
    auto generate_iterators_impl(Functor f, std::index_sequence<Is...>)
        -> std::array<decltype(f(0)), sizeof...(Is)>
    {
      return {{f(Is)...}};
    }
  } /* namespace */
#endif
 
 
  template <unsigned int length, typename Functor>
  DEAL_II_ALWAYS_INLINE inline auto generate_iterators(Functor f)
      -> std::array<decltype(f(0)), length>
  {
    return generate_iterators_impl<>(f, std::make_index_sequence<length>());
  }
 
 
  template <typename T>
  DEAL_II_ALWAYS_INLINE inline void increment_iterators(T &iterators)
  {
    for (auto &it : iterators)
      it++;
  }
 
 
 
  template <typename Number>
  inline DEAL_II_ALWAYS_INLINE Number positive_part(const Number number)
  {
    return Number(0.5) * (std::abs(number) + number);
  }
 
 
  template <typename Number>
  inline DEAL_II_ALWAYS_INLINE Number negative_part(const Number number)
  {
    return Number(0.5) * (std::abs(number) - number);
  }
 
 
  template <int N, typename T>
  inline T fixed_power(const T x)
  {
    return dealii::Utilities::fixed_power<N, T>(x);
  }
 
 
  template <typename T>
  T pow(const T x, const T b);
 
 
  template <typename T, std::size_t width>
  dealii::VectorizedArray<T, width>
  pow(const dealii::VectorizedArray<T, width> x, const T b);
 
 
  template <typename T, std::size_t width>
  dealii::VectorizedArray<T, width>
  pow(const dealii::VectorizedArray<T, width> x,
      const dealii::VectorizedArray<T, width> b);
 
 
  enum class Bias {
    none,
 
    max,
 
    min
  };
 
 
  template <typename T>
  T fast_pow(const T x, const T b, const Bias bias = Bias::none);
 
 
  template <typename T, std::size_t width>
  dealii::VectorizedArray<T, width>
  fast_pow(const dealii::VectorizedArray<T, width> x,
           const T b,
           const Bias bias = Bias::none);
 
 
  template <typename T, std::size_t width>
  dealii::VectorizedArray<T, width>
  fast_pow(const dealii::VectorizedArray<T, width> x,
           const dealii::VectorizedArray<T, width> b,
           const Bias bias = Bias::none);
 
 
 
  template <typename T, typename V>
  DEAL_II_ALWAYS_INLINE inline T load_value(const V &vector, unsigned int i)
  {
    static_assert(std::is_same_v<typename get_value_type<T>::type,
                                 typename V::value_type>,
                  "type mismatch");
    T result;
 
    if constexpr (std::is_same_v<T, typename get_value_type<T>::type>) {
      /* Non-vectorized sequential access. */
      result = vector.local_element(i);
    } else {
      /* Vectorized fast access. index must be divisible by simd_length */
      result.load(vector.get_values() + i);
    }
 
    return result;
  }
 
 
  template <typename T, typename T2>
  DEAL_II_ALWAYS_INLINE inline T load_value(const std::vector<T2> &vector,
                                            unsigned int i)
  {
    if constexpr (std::is_same_v<typename get_value_type<T>::type, T2>) {
      /* Optimized default for source and destination with same type: */
 
      T result;
      if constexpr (std::is_same_v<T, typename get_value_type<T>::type>) {
        /* Non-vectorized sequential access. */
        result = vector[i];
      } else {
        /* Vectorized fast access. index must be divisible by simd_length */
        result.load(vector.data() + i);
      }
      return result;
 
    } else {
      /* Fallback for mismatched types (float vs double): */
      T result;
      if constexpr (std::is_same_v<T, typename get_value_type<T>::type>) {
        result = vector[i];
      } else {
        // FIXME: suboptimal
        for (unsigned int k = 0; k < T::size(); ++k)
          result[k] = vector[i + k];
      }
      return result;
    }
  }
 
 
  template <typename T, typename V>
  DEAL_II_ALWAYS_INLINE inline T load_value(const V &vector,
                                            const unsigned int *js)
  {
    static_assert(std::is_same_v<typename get_value_type<T>::type,
                                 typename V::value_type>,
                  "type mismatch");
    T result;
 
    if constexpr (std::is_same_v<T, typename get_value_type<T>::type>) {
      /* Non-vectorized sequential access. */
      result = vector.local_element(js[0]);
    } else {
      /* Vectorized fast access. index must be divisible by simd_length */
      result.gather(vector.get_values(), js);
    }
 
    return result;
  }
 
 
  template <typename T, typename T2>
  DEAL_II_ALWAYS_INLINE inline T load_value(const std::vector<T2> &vector,
                                            const unsigned int *js)
  {
    static_assert(std::is_same_v<typename get_value_type<T>::type, T2>,
                  "type mismatch");
    T result;
 
    if constexpr (std::is_same_v<T, typename get_value_type<T>::type>) {
      /* Non-vectorized sequential access. */
      result = vector[js[0]];
    } else {
      /* Vectorized fast access. index must be divisible by simd_length */
      result.load(vector.data(), js);
    }
 
    return result;
  }
 
 
  template <typename T, typename V>
  DEAL_II_ALWAYS_INLINE inline void
  store_value(V &vector, const T &values, unsigned int i)
  {
    static_assert(std::is_same_v<typename get_value_type<T>::type,
                                 typename V::value_type>,
                  "type mismatch");
 
    if constexpr (std::is_same_v<T, typename get_value_type<T>::type>) {
      /* Non-vectorized sequential access. */
      vector.local_element(i) = values;
    } else {
      /* Vectorized fast access. index must be divisible by simd_length */
      values.store(vector.get_values() + i);
    }
  }
 
 
  template <typename T, typename T2>
  DEAL_II_ALWAYS_INLINE inline void
  store_value(std::vector<T2> &vector, const T &values, unsigned int i)
  {
    if constexpr (std::is_same_v<typename get_value_type<T>::type, T2>) {
      /* Optimized default for source and destination with same type: */
 
      if constexpr (std::is_same_v<T, typename get_value_type<T>::type>) {
        /* Non-vectorized sequential access. */
        vector[i] = values;
      } else {
        /* Vectorized fast access. index must be divisible by simd_length */
        values.store(vector.data() + i);
      }
 
    } else {
      /* Fallback for mismatched types (float vs double): */
      if constexpr (std::is_same_v<T, typename get_value_type<T>::type>) {
        vector[i] = values;
      } else {
        // FIXME: suboptimal
        for (unsigned int k = 0; k < T::size(); ++k)
          vector[i + k] = values[k];
      }
    }
  }
 
 
  template <int rank, int dim, std::size_t width, typename Number>
  DEAL_II_ALWAYS_INLINE inline dealii::Tensor<rank, dim, Number>
  serialize_tensor(
      const dealii::Tensor<rank, dim, dealii::VectorizedArray<Number, width>>
          &vectorized,
      const unsigned int k)
  {
    Assert(k < width, dealii::ExcMessage("Index past VectorizedArray width"));
    dealii::Tensor<rank, dim, Number> result;
    if constexpr (rank == 1) {
      for (unsigned int d = 0; d < dim; ++d)
        result[d] = vectorized[d][k];
    } else {
      for (unsigned int d = 0; d < dim; ++d)
        result[d] = serialize_tensor(vectorized[d], k);
    }
    return result;
  }
 
 
  template <int rank, int dim, typename Number>
  DEAL_II_ALWAYS_INLINE inline dealii::Tensor<rank, dim, Number>
  serialize_tensor(const dealii::Tensor<rank, dim, Number> &serial,
                   const unsigned int k [[maybe_unused]])
  {
    Assert(k == 0,
           dealii::ExcMessage(
               "The given index k must be zero for a serial tensor"));
    return serial;
  }
 
 
  template <int rank, int dim, std::size_t width, typename Number>
  DEAL_II_ALWAYS_INLINE inline void assign_serial_tensor(
      dealii::Tensor<rank, dim, dealii::VectorizedArray<Number, width>> &result,
      const dealii::Tensor<rank, dim, Number> &serial,
      const unsigned int k)
  {
    Assert(k < width, dealii::ExcMessage("Index past VectorizedArray width"));
    if constexpr (rank == 1) {
      for (unsigned int d = 0; d < dim; ++d)
        result[d][k] = serial[d];
    } else {
      for (unsigned int d = 0; d < dim; ++d)
        assign_serial_tensor(result[d], serial[d], k);
    }
  }
 
 
  template <int rank, int dim, typename Number>
  DEAL_II_ALWAYS_INLINE inline void
  assign_serial_tensor(dealii::Tensor<rank, dim, Number> &result,
                       const dealii::Tensor<rank, dim, Number> &serial,
                       const unsigned int k [[maybe_unused]])
  {
    Assert(k == 0,
           dealii::ExcMessage(
               "The given index k must be zero for a serial tensor"));
 
    result = serial;
  }
 
 
} // namespace ryujin