vectorization.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2012 - 2025 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
 
#ifndef dealii_vectorization_h
#define dealii_vectorization_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/numbers.h>
#include <deal.II/base/template_constraints.h>
 
#include <algorithm>
#include <array>
#include <cmath>
 
// Note:
// The flag DEAL_II_VECTORIZATION_WIDTH_IN_BITS is essentially constructed
// according to the following scheme (on x86-based architectures)
// #ifdef __AVX512F__
// #define DEAL_II_VECTORIZATION_WIDTH_IN_BITS 512
// #elif defined (__AVX__)
// #define DEAL_II_VECTORIZATION_WIDTH_IN_BITS 256
// #elif defined (__SSE2__)
// #define DEAL_II_VECTORIZATION_WIDTH_IN_BITS 128
// #else
// #define DEAL_II_VECTORIZATION_WIDTH_IN_BITS 0
// #endif
// In addition to checking the flags __AVX512F__, __AVX__ and __SSE2__, a CMake
// test, 'check_01_cpu_features.cmake', ensures that these feature are not only
// present in the compilation unit but also working properly.
 
#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS > 0
 
// These error messages try to detect the case that deal.II was compiled with
// a wider instruction set extension as the current compilation unit, for
// example because deal.II was compiled with AVX, but a user project does not
// add -march=native or similar flags, making it fall to SSE2. This leads to
// very strange errors as the size of data structures differs between the
// compiled deal.II code sitting in libdeal_II.so and the user code if not
// detected.
#  if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && !defined(__AVX__)
#    error \
      "Mismatch in vectorization capabilities: AVX was detected during configuration of deal.II and switched on, but it is apparently not available for the file you are trying to compile at the moment. Check compilation flags controlling the instruction set, such as -march=native."
#  endif
#  if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && !defined(__AVX512F__)
#    error \
      "Mismatch in vectorization capabilities: AVX-512F was detected during configuration of deal.II and switched on, but it is apparently not available for the file you are trying to compile at the moment. Check compilation flags controlling the instruction set, such as -march=native."
#  endif
 
#  ifdef _MSC_VER
#    include <intrin.h>
#  elif defined(__ALTIVEC__)
#    include <altivec.h>
 
// altivec.h defines vector, pixel, bool, but we do not use them, so undefine
// them before they make trouble
#    undef vector
#    undef pixel
#    undef bool
#  elif defined(__ARM_NEON)
#    include <arm_neon.h>
#  elif defined(__x86_64__)
#    include <x86intrin.h>
#  endif
 
#endif
 
 
DEAL_II_NAMESPACE_OPEN
 
 
// Enable the EnableIfScalar type trait for VectorizedArray<Number> such
// that it can be used as a Number type in Tensor<rank,dim,Number>, etc.
 
template <typename Number, std::size_t width>
struct EnableIfScalar<VectorizedArray<Number, width>>
{
  using type = VectorizedArray<typename EnableIfScalar<Number>::type, width>;
};
 
 

template <typename T>
class VectorizedArrayIterator
{
public:
  constexpr VectorizedArrayIterator(T &data, const std::size_t lane)
    : data(&data)
    , lane(lane)
  {}

  constexpr bool
  operator==(const VectorizedArrayIterator<T> &other) const
  {
    Assert(this->data == other.data,
           ExcMessage(
             "You are trying to compare iterators into different arrays."));
    return this->lane == other.lane;
  }

  constexpr bool
  operator!=(const VectorizedArrayIterator<T> &other) const
  {
    Assert(this->data == other.data,
           ExcMessage(
             "You are trying to compare iterators into different arrays."));
    return this->lane != other.lane;
  }

  constexpr const typename T::value_type &
  operator*() const
  {
    AssertIndexRange(lane, T::size());
    return (*data)[lane];
  }
 

  template <typename U = T>
  constexpr std::enable_if_t<!std::is_same_v<U, const U>,
                             typename T::value_type> &
  operator*()
  {
    AssertIndexRange(lane, T::size());
    return (*data)[lane];
  }

  constexpr VectorizedArrayIterator<T> &
  operator++()
  {
    AssertIndexRange(lane + 1, T::size() + 1);
    ++lane;
    return *this;
  }

  constexpr VectorizedArrayIterator<T> &
  operator+=(const std::size_t offset)
  {
    AssertIndexRange(lane + offset, T::size() + 1);
    lane += offset;
    return *this;
  }

  constexpr VectorizedArrayIterator<T> &
  operator--()
  {
    Assert(
      lane > 0,
      ExcMessage(
        "You can't decrement an iterator that is already at the beginning of the range."));
    --lane;
    return *this;
  }

  constexpr VectorizedArrayIterator<T>
  operator+(const std::size_t &offset) const
  {
    AssertIndexRange(lane + offset, T::size() + 1);
    return VectorizedArrayIterator<T>(*data, lane + offset);
  }

  constexpr std::ptrdiff_t
  operator-(const VectorizedArrayIterator<T> &other) const
  {
    return static_cast<std::ptrdiff_t>(lane) -
           static_cast<std::ptrdiff_t>(other.lane);
  }
 
private:
  T *data;

  std::size_t lane;
};
 
 

template <typename VectorizedArrayType, std::size_t width>
class VectorizedArrayBase
{
public:
  constexpr VectorizedArrayBase() = default;

  template <typename U>
  constexpr VectorizedArrayBase(const std::initializer_list<U> &list)
  {
    const unsigned int n_initializers = list.size();
    Assert(n_initializers <= size(),
           ExcMessage("The initializer list must have at most "
                      "as many elements as the vector length."));
 
    // Copy what's in the list.
    std::copy_n(list.begin(), n_initializers, this->begin());
 
    // Then add zero padding where necessary.
    if (n_initializers < size())
      std::fill(this->begin() + n_initializers, this->end(), 0.0);
  }

  static constexpr std::size_t
  size()
  {
    return width;
  }

  constexpr VectorizedArrayIterator<VectorizedArrayType>
  begin()
  {
    return VectorizedArrayIterator<VectorizedArrayType>(
      static_cast<VectorizedArrayType &>(*this), 0);
  }

  constexpr VectorizedArrayIterator<const VectorizedArrayType>
  begin() const
  {
    return VectorizedArrayIterator<const VectorizedArrayType>(
      static_cast<const VectorizedArrayType &>(*this), 0);
  }

  constexpr VectorizedArrayIterator<VectorizedArrayType>
  end()
  {
    return VectorizedArrayIterator<VectorizedArrayType>(
      static_cast<VectorizedArrayType &>(*this), width);
  }

  constexpr VectorizedArrayIterator<const VectorizedArrayType>
  end() const
  {
    return VectorizedArrayIterator<const VectorizedArrayType>(
      static_cast<const VectorizedArrayType &>(*this), width);
  }

  auto
  dot_product(const VectorizedArrayType &v) const
  {
    VectorizedArrayType p = static_cast<const VectorizedArrayType &>(*this);
    p *= v;
    return p.sum();
  }
};
 
 

template <typename Number, std::size_t width>
class VectorizedArray
  : public VectorizedArrayBase<VectorizedArray<Number, width>, 1>
{
public:
  using value_type = Number;

  static constexpr bool is_implemented = (width == 1);

  VectorizedArray() = default;

  VectorizedArray(const Number scalar)
  {
    static_assert(width == 1,
                  "You specified an illegal width that is not supported.");
 
    this->operator=(scalar);
  }

  template <typename U>
  VectorizedArray(const std::initializer_list<U> &list)
    : VectorizedArrayBase<VectorizedArray<Number, width>, 1>(list)
  {
    static_assert(width == 1,
                  "You specified an illegal width that is not supported.");
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator=(const Number scalar) &
  {
    data = scalar;
    return *this;
  }

  VectorizedArray &
  operator=(const Number scalar) && = delete;

  DEAL_II_ALWAYS_INLINE
  Number &
  operator[](const unsigned int comp)
  {
    (void)comp;
    AssertIndexRange(comp, 1);
    return data;
  }

  DEAL_II_ALWAYS_INLINE
  const Number &
  operator[](const unsigned int comp) const
  {
    (void)comp;
    AssertIndexRange(comp, 1);
    return data;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator+=(const VectorizedArray &vec)
  {
    data += vec.data;
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator-=(const VectorizedArray &vec)
  {
    data -= vec.data;
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator*=(const VectorizedArray &vec)
  {
    data *= vec.data;
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator/=(const VectorizedArray &vec)
  {
    data /= vec.data;
    return *this;
  }

  template <typename OtherNumber>
  DEAL_II_ALWAYS_INLINE void
  load(const OtherNumber *ptr)
  {
    data = *ptr;
  }

  template <typename OtherNumber>
  DEAL_II_ALWAYS_INLINE void
  store(OtherNumber *ptr) const
  {
    *ptr = data;
  }

  DEAL_II_ALWAYS_INLINE
  void
  streaming_store(Number *ptr) const
  {
    *ptr = data;
  }

  DEAL_II_ALWAYS_INLINE
  void
  gather(const Number *base_ptr, const unsigned int *offsets)
  {
    this->operator=(Number(0));
    if (offsets[0] != numbers::invalid_unsigned_int)
      data = base_ptr[offsets[0]];
  }

  DEAL_II_ALWAYS_INLINE
  void
  scatter(const unsigned int *offsets, Number *base_ptr) const
  {
    if (offsets[0] != numbers::invalid_unsigned_int)
      base_ptr[offsets[0]] = data;
  }

  DEAL_II_ALWAYS_INLINE
  Number
  sum() const
  {
    return data;
  }

  VectorizedArray
  get_floor() const
  {
    VectorizedArray res;
    res.data = std::floor(data);
    return res;
  }

  Number data;
 
private:
  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_sqrt() const
  {
    VectorizedArray res;
    res.data = std::sqrt(data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_abs() const
  {
    VectorizedArray res;
    res.data = std::fabs(data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_max(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = std::max(data, other.data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_min(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = std::min(data, other.data);
    return res;
  }
 
  // Make a few functions friends.
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::sqrt(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::abs(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::max(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::min(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
};
 
 


template <typename Number,
          std::size_t width =
            internal::VectorizedArrayWidthSpecifier<Number>::max_width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number, width>
                             make_vectorized_array(const Number &u)
{
  VectorizedArray<Number, width> result = u;
  return result;
}
 
 

template <typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
make_vectorized_array(const typename VectorizedArrayType::value_type &u)
{
  static_assert(
    std::is_same_v<VectorizedArrayType,
                   VectorizedArray<typename VectorizedArrayType::value_type,
                                   VectorizedArrayType::size()>>,
    "VectorizedArrayType is not a VectorizedArray.");
 
  VectorizedArrayType result = u;
  return result;
}
 
 

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE void
gather(VectorizedArray<Number, width>          &out,
       const std::array<const Number *, width> &ptrs,
       const unsigned int                       offset)
{
  for (unsigned int v = 0; v < width; ++v)
    out.data[v] = ptrs[v][offset];
}
 
 

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE void
vectorized_load_and_transpose(const unsigned int              n_entries,
                              const Number                   *in,
                              const unsigned int             *offsets,
                              VectorizedArray<Number, width> *out)
{
  for (unsigned int i = 0; i < n_entries; ++i)
    for (unsigned int v = 0; v < VectorizedArray<Number, width>::size(); ++v)
      out[i][v] = in[offsets[v] + i];
}
 

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE void
vectorized_load_and_transpose(const unsigned int n_entries,
                              const std::array<const Number *, width> &in,
                              VectorizedArray<Number, width>          *out)
{
  for (unsigned int i = 0; i < n_entries; ++i)
    for (unsigned int v = 0; v < VectorizedArray<Number, width>::size(); ++v)
      out[i][v] = in[v][i];
}
 
 

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE void
vectorized_transpose_and_store(const bool                            add_into,
                               const unsigned int                    n_entries,
                               const VectorizedArray<Number, width> *in,
                               const unsigned int                   *offsets,
                               Number                               *out)
{
  if (add_into)
    for (unsigned int i = 0; i < n_entries; ++i)
      for (unsigned int v = 0; v < VectorizedArray<Number, width>::size(); ++v)
        out[offsets[v] + i] += in[i][v];
  else
    for (unsigned int i = 0; i < n_entries; ++i)
      for (unsigned int v = 0; v < VectorizedArray<Number, width>::size(); ++v)
        out[offsets[v] + i] = in[i][v];
}
 

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE void
vectorized_transpose_and_store(const bool                            add_into,
                               const unsigned int                    n_entries,
                               const VectorizedArray<Number, width> *in,
                               const std::array<Number *, width>    &out)
{
  if (add_into)
    for (unsigned int i = 0; i < n_entries; ++i)
      for (unsigned int v = 0; v < VectorizedArray<Number, width>::size(); ++v)
        out[v][i] += in[i][v];
  else
    for (unsigned int i = 0; i < n_entries; ++i)
      for (unsigned int v = 0; v < VectorizedArray<Number, width>::size(); ++v)
        out[v][i] = in[i][v];
}
 

 
#ifndef DOXYGEN
 
#  if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ARM_NEON)

template <>
class VectorizedArray<double, 2>
  : public VectorizedArrayBase<VectorizedArray<double, 2>, 2>
{
public:
  using value_type = double;

  static constexpr bool is_implemented = true;

  VectorizedArray() = default;

  VectorizedArray(const double scalar)
  {
    this->operator=(scalar);
  }

  template <typename U>
  VectorizedArray(const std::initializer_list<U> &list)
    : VectorizedArrayBase<VectorizedArray<double, 2>, 2>(list)
  {}

  VectorizedArray &
  operator=(const double x) &
  {
    data = vdupq_n_f64(x);
    return *this;
  }

  VectorizedArray &
  operator=(const double scalar) && = delete;

  double &
  operator[](const unsigned int comp)
  {
    return *(reinterpret_cast<double *>(&data) + comp);
  }

  const double &
  operator[](const unsigned int comp) const
  {
    return *(reinterpret_cast<const double *>(&data) + comp);
  }

  VectorizedArray &
  operator+=(const VectorizedArray &vec)
  {
    data = vaddq_f64(data, vec.data);
    return *this;
  }

  VectorizedArray &
  operator-=(const VectorizedArray &vec)
  {
    data = vsubq_f64(data, vec.data);
    return *this;
  }

  VectorizedArray &
  operator*=(const VectorizedArray &vec)
  {
    data = vmulq_f64(data, vec.data);
    return *this;
  }

  VectorizedArray &
  operator/=(const VectorizedArray &vec)
  {
    data = vdivq_f64(data, vec.data);
    return *this;
  }

  void
  load(const double *ptr)
  {
    data = vld1q_f64(ptr);
  }
 
  DEAL_II_ALWAYS_INLINE
  void
  load(const float *ptr)
  {
    DEAL_II_OPENMP_SIMD_PRAGMA
    for (unsigned int i = 0; i < 2; ++i)
      data[i] = ptr[i];
  }

  void
  store(double *ptr) const
  {
    vst1q_f64(ptr, data);
  }
 
  DEAL_II_ALWAYS_INLINE
  void
  store(float *ptr) const
  {
    DEAL_II_OPENMP_SIMD_PRAGMA
    for (unsigned int i = 0; i < 2; ++i)
      ptr[i] = data[i];
  }

  DEAL_II_ALWAYS_INLINE
  void
  streaming_store(double *ptr) const
  {
    Assert(reinterpret_cast<std::size_t>(ptr) % 16 == 0,
           ExcMessage("Memory not aligned"));
    vst1q_f64(ptr, data);
  }

  void
  gather(const double *base_ptr, const unsigned int *offsets)
  {
    this->operator=(0.);
    for (unsigned int i = 0; i < 2; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        *(reinterpret_cast<double *>(&data) + i) = base_ptr[offsets[i]];
  }

  void
  scatter(const unsigned int *offsets, double *base_ptr) const
  {
    for (unsigned int i = 0; i < 2; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        base_ptr[offsets[i]] = *(reinterpret_cast<const double *>(&data) + i);
  }

  double
  sum() const
  {
    return vaddvq_f64(data);
  }

  VectorizedArray
  get_floor() const
  {
    VectorizedArray res;
    res.data = vrndmq_f64(data);
    return res;
  }

  mutable float64x2_t data;
 
private:
  VectorizedArray
  get_sqrt() const
  {
    VectorizedArray res;
    res.data = vsqrtq_f64(data);
    return res;
  }

  VectorizedArray
  get_abs() const
  {
    VectorizedArray res;
    res.data = vabsq_f64(data);
    return res;
  }

  VectorizedArray
  get_max(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = vmaxq_f64(data, other.data);
    return res;
  }

  VectorizedArray
  get_min(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = vminq_f64(data, other.data);
    return res;
  }
 
  // Make a few functions friends.
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::sqrt(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::abs(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::max(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::min(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
};

template <>
class VectorizedArray<float, 4>
  : public VectorizedArrayBase<VectorizedArray<float, 4>, 4>
{
public:
  using value_type = float;

  static constexpr bool is_implemented = true;

  VectorizedArray() = default;

  VectorizedArray(const float scalar)
  {
    this->operator=(scalar);
  }

  template <typename U>
  VectorizedArray(const std::initializer_list<U> &list)
    : VectorizedArrayBase<VectorizedArray<float, 4>, 4>(list)
  {}

  VectorizedArray &
  operator=(const float x) &
  {
    data = vdupq_n_f32(x);
    return *this;
  }

  VectorizedArray &
  operator=(const float scalar) && = delete;

  value_type &
  operator[](const unsigned int comp)
  {
    return *(reinterpret_cast<float *>(&data) + comp);
  }

  const value_type &
  operator[](const unsigned int comp) const
  {
    return *(reinterpret_cast<const float *>(&data) + comp);
  }

  VectorizedArray &
  operator+=(const VectorizedArray &vec)
  {
    data = vaddq_f32(data, vec.data);
    return *this;
  }

  VectorizedArray &
  operator-=(const VectorizedArray &vec)
  {
    data = vsubq_f32(data, vec.data);
    return *this;
  }

  VectorizedArray &
  operator*=(const VectorizedArray &vec)
  {
    data = vmulq_f32(data, vec.data);
    return *this;
  }

  VectorizedArray &
  operator/=(const VectorizedArray &vec)
  {
    data = vdivq_f32(data, vec.data);
    return *this;
  }

  void
  load(const float *ptr)
  {
    data = vld1q_f32(ptr);
  }

  void
  store(float *ptr) const
  {
    vst1q_f32(ptr, data);
  }

  DEAL_II_ALWAYS_INLINE
  void
  streaming_store(float *ptr) const
  {
    Assert(reinterpret_cast<std::size_t>(ptr) % 16 == 0,
           ExcMessage("Memory not aligned"));
    vst1q_f32(ptr, data);
  }

  void
  gather(const float *base_ptr, const unsigned int *offsets)
  {
    this->operator=(0.f);
    for (unsigned int i = 0; i < 4; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        *(reinterpret_cast<float *>(&data) + i) = base_ptr[offsets[i]];
  }

  void
  scatter(const unsigned int *offsets, float *base_ptr) const
  {
    for (unsigned int i = 0; i < 4; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        base_ptr[offsets[i]] = *(reinterpret_cast<const float *>(&data) + i);
  }

  float
  sum() const
  {
    return vaddvq_f32(data);
  }
 

  VectorizedArray
  get_floor() const
  {
    VectorizedArray res;
    res.data = vrndmq_f32(data);
    return res;
  }

  mutable float32x4_t data;
 
private:
  VectorizedArray
  get_sqrt() const
  {
    VectorizedArray res;
    res.data = vsqrtq_f32(data);
    return res;
  }

  VectorizedArray
  get_abs() const
  {
    VectorizedArray res;
    res.data = vabsq_f32(data);
    return res;
  }

  VectorizedArray
  get_max(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = vmaxq_f32(data, other.data);
    return res;
  }

  VectorizedArray
  get_min(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = vminq_f32(data, other.data);
    return res;
  }
 
  // Make a few functions friends.
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::sqrt(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::abs(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::max(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::min(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
};
 
 
#  endif
 
#  if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__)

template <>
class VectorizedArray<double, 2>
  : public VectorizedArrayBase<VectorizedArray<double, 2>, 2>
{
public:
  using value_type = double;

  static constexpr bool is_implemented = true;

  VectorizedArray() = default;

  VectorizedArray(const double scalar)
  {
    this->operator=(scalar);
  }

  template <typename U>
  VectorizedArray(const std::initializer_list<U> &list)
    : VectorizedArrayBase<VectorizedArray<double, 2>, 2>(list)
  {}

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator=(const double x) &
  {
    data = _mm_set1_pd(x);
    return *this;
  }

  VectorizedArray &
  operator=(const double scalar) && = delete;

  DEAL_II_ALWAYS_INLINE
  double &
  operator[](const unsigned int comp)
  {
    AssertIndexRange(comp, 2);
    return *(reinterpret_cast<double *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  const double &
  operator[](const unsigned int comp) const
  {
    AssertIndexRange(comp, 2);
    return *(reinterpret_cast<const double *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator+=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data += vec.data;
#    else
    data = _mm_add_pd(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator-=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data -= vec.data;
#    else
    data = _mm_sub_pd(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator*=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data *= vec.data;
#    else
    data = _mm_mul_pd(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator/=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data /= vec.data;
#    else
    data = _mm_div_pd(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  void
  load(const double *ptr)
  {
    data = _mm_loadu_pd(ptr);
  }
 
  DEAL_II_ALWAYS_INLINE
  void
  load(const float *ptr)
  {
    DEAL_II_OPENMP_SIMD_PRAGMA
    for (unsigned int i = 0; i < 2; ++i)
      data[i] = ptr[i];
  }

  DEAL_II_ALWAYS_INLINE
  void
  store(double *ptr) const
  {
    _mm_storeu_pd(ptr, data);
  }
 
  DEAL_II_ALWAYS_INLINE
  void
  store(float *ptr) const
  {
    DEAL_II_OPENMP_SIMD_PRAGMA
    for (unsigned int i = 0; i < 2; ++i)
      ptr[i] = data[i];
  }

  DEAL_II_ALWAYS_INLINE
  void
  streaming_store(double *ptr) const
  {
    Assert(reinterpret_cast<std::size_t>(ptr) % 16 == 0,
           ExcMessage("Memory not aligned"));
    _mm_stream_pd(ptr, data);
  }

  DEAL_II_ALWAYS_INLINE
  void
  gather(const double *base_ptr, const unsigned int *offsets)
  {
    this->operator=(0.);
    for (unsigned int i = 0; i < 2; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        *(reinterpret_cast<double *>(&data) + i) = base_ptr[offsets[i]];
  }

  DEAL_II_ALWAYS_INLINE
  void
  scatter(const unsigned int *offsets, double *base_ptr) const
  {
    for (unsigned int i = 0; i < 2; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        base_ptr[offsets[i]] = *(reinterpret_cast<const double *>(&data) + i);
  }

  double
  sum() const
  {
    __m128d t1 = _mm_unpackhi_pd(data, data);
    __m128d t2 = _mm_add_pd(data, t1);
    return _mm_cvtsd_f64(t2);
  }

  VectorizedArray
  get_floor() const
  {
    VectorizedArray res;
    for (std::size_t i = 0; i < 2; ++i)
      res.data[i] = std::floor(data[i]);
    return res;
  }

  __m128d data;
 
private:
  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_sqrt() const
  {
    VectorizedArray res;
    res.data = _mm_sqrt_pd(data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_abs() const
  {
    // to compute the absolute value, perform
    // bitwise andnot with -0. This will leave all
    // value and exponent bits unchanged but force
    // the sign value to +.
    __m128d         mask = _mm_set1_pd(-0.);
    VectorizedArray res;
    res.data = _mm_andnot_pd(mask, data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_max(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = _mm_max_pd(data, other.data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_min(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = _mm_min_pd(data, other.data);
    return res;
  }
 
  // Make a few functions friends.
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::sqrt(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::abs(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::max(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::min(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
};
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_load_and_transpose(const unsigned int                   n_entries,
                              const std::array<const double *, 2> &in,
                              VectorizedArray<double, 2>          *out)
{
  const unsigned int n_chunks = n_entries / 2;
  for (unsigned int i = 0; i < n_chunks; ++i)
    {
      __m128d u0          = _mm_loadu_pd(in[0] + 2 * i);
      __m128d u1          = _mm_loadu_pd(in[1] + 2 * i);
      out[2 * i + 0].data = _mm_unpacklo_pd(u0, u1);
      out[2 * i + 1].data = _mm_unpackhi_pd(u0, u1);
    }
 
  // remainder loop of work that does not divide by 2
  for (unsigned int i = 2 * n_chunks; i < n_entries; ++i)
    for (unsigned int v = 0; v < 2; ++v)
      out[i][v] = in[v][i];
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_load_and_transpose(const unsigned int          n_entries,
                              const double               *in,
                              const unsigned int         *offsets,
                              VectorizedArray<double, 2> *out)
{
  vectorized_load_and_transpose(n_entries,
                                {{in + offsets[0], in + offsets[1]}},
                                out);
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_transpose_and_store(const bool                        add_into,
                               const unsigned int                n_entries,
                               const VectorizedArray<double, 2> *in,
                               const std::array<double *, 2>    &out)
{
  const unsigned int n_chunks = n_entries / 2;
  if (add_into)
    {
      for (unsigned int i = 0; i < n_chunks; ++i)
        {
          __m128d u0   = in[2 * i + 0].data;
          __m128d u1   = in[2 * i + 1].data;
          __m128d res0 = _mm_unpacklo_pd(u0, u1);
          __m128d res1 = _mm_unpackhi_pd(u0, u1);
          _mm_storeu_pd(out[0] + 2 * i,
                        _mm_add_pd(_mm_loadu_pd(out[0] + 2 * i), res0));
          _mm_storeu_pd(out[1] + 2 * i,
                        _mm_add_pd(_mm_loadu_pd(out[1] + 2 * i), res1));
        }
 
      // remainder loop for lengths that do not divide by 2
      for (unsigned int i = 2 * n_chunks; i < n_entries; ++i)
        for (unsigned int v = 0; v < 2; ++v)
          out[v][i] += in[i][v];
    }
  else
    {
      for (unsigned int i = 0; i < n_chunks; ++i)
        {
          __m128d u0   = in[2 * i + 0].data;
          __m128d u1   = in[2 * i + 1].data;
          __m128d res0 = _mm_unpacklo_pd(u0, u1);
          __m128d res1 = _mm_unpackhi_pd(u0, u1);
          _mm_storeu_pd(out[0] + 2 * i, res0);
          _mm_storeu_pd(out[1] + 2 * i, res1);
        }
 
      // remainder loop for lengths that do not divide by 2
      for (unsigned int i = 2 * n_chunks; i < n_entries; ++i)
        for (unsigned int v = 0; v < 2; ++v)
          out[v][i] = in[i][v];
    }
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_transpose_and_store(const bool                        add_into,
                               const unsigned int                n_entries,
                               const VectorizedArray<double, 2> *in,
                               const unsigned int               *offsets,
                               double                           *out)
{
  vectorized_transpose_and_store(add_into,
                                 n_entries,
                                 in,
                                 {{out + offsets[0], out + offsets[1]}});
}
 
 

template <>
class VectorizedArray<float, 4>
  : public VectorizedArrayBase<VectorizedArray<float, 4>, 4>
{
public:
  using value_type = float;

  static constexpr bool is_implemented = true;

  VectorizedArray() = default;

  VectorizedArray(const float scalar)
  {
    this->operator=(scalar);
  }

  template <typename U>
  VectorizedArray(const std::initializer_list<U> &list)
    : VectorizedArrayBase<VectorizedArray<float, 4>, 4>(list)
  {}

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator=(const float x) &
  {
    data = _mm_set1_ps(x);
    return *this;
  }

  VectorizedArray &
  operator=(const float scalar) && = delete;

  DEAL_II_ALWAYS_INLINE
  float &
  operator[](const unsigned int comp)
  {
    AssertIndexRange(comp, 4);
    return *(reinterpret_cast<float *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  const float &
  operator[](const unsigned int comp) const
  {
    AssertIndexRange(comp, 4);
    return *(reinterpret_cast<const float *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator+=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data += vec.data;
#    else
    data = _mm_add_ps(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator-=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data -= vec.data;
#    else
    data = _mm_sub_ps(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator*=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data *= vec.data;
#    else
    data = _mm_mul_ps(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator/=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data /= vec.data;
#    else
    data = _mm_div_ps(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  void
  load(const float *ptr)
  {
    data = _mm_loadu_ps(ptr);
  }

  DEAL_II_ALWAYS_INLINE
  void
  store(float *ptr) const
  {
    _mm_storeu_ps(ptr, data);
  }

  DEAL_II_ALWAYS_INLINE
  void
  streaming_store(float *ptr) const
  {
    Assert(reinterpret_cast<std::size_t>(ptr) % 16 == 0,
           ExcMessage("Memory not aligned"));
    _mm_stream_ps(ptr, data);
  }

  DEAL_II_ALWAYS_INLINE
  void
  gather(const float *base_ptr, const unsigned int *offsets)
  {
    this->operator=(0.f);
    for (unsigned int i = 0; i < 4; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        *(reinterpret_cast<float *>(&data) + i) = base_ptr[offsets[i]];
  }

  DEAL_II_ALWAYS_INLINE
  void
  scatter(const unsigned int *offsets, float *base_ptr) const
  {
    for (unsigned int i = 0; i < 4; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        base_ptr[offsets[i]] = *(reinterpret_cast<const float *>(&data) + i);
  }

  float
  sum() const
  {
    __m128 t1 = _mm_movehl_ps(data, data);
    __m128 t2 = _mm_add_ps(data, t1);
    __m128 t3 = _mm_shuffle_ps(t2, t2, 1);
    __m128 t4 = _mm_add_ss(t2, t3);
    return _mm_cvtss_f32(t4);
  }

  VectorizedArray
  get_floor() const
  {
    VectorizedArray res;
    for (std::size_t i = 0; i < 4; ++i)
      res.data[i] = std::floor(data[i]);
    return res;
  }

  __m128 data;
 
private:
  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_sqrt() const
  {
    VectorizedArray res;
    res.data = _mm_sqrt_ps(data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_abs() const
  {
    // to compute the absolute value, perform bitwise andnot with -0. This
    // will leave all value and exponent bits unchanged but force the sign
    // value to +.
    __m128          mask = _mm_set1_ps(-0.f);
    VectorizedArray res;
    res.data = _mm_andnot_ps(mask, data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_max(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = _mm_max_ps(data, other.data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_min(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = _mm_min_ps(data, other.data);
    return res;
  }
 
  // Make a few functions friends.
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::sqrt(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::abs(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::max(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::min(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
};
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_load_and_transpose(const unsigned int                  n_entries,
                              const std::array<const float *, 4> &in,
                              VectorizedArray<float, 4>          *out)
{
  const unsigned int n_chunks = n_entries / 4;
  for (unsigned int i = 0; i < n_chunks; ++i)
    {
      __m128 u0           = _mm_loadu_ps(in[0] + 4 * i);
      __m128 u1           = _mm_loadu_ps(in[1] + 4 * i);
      __m128 u2           = _mm_loadu_ps(in[2] + 4 * i);
      __m128 u3           = _mm_loadu_ps(in[3] + 4 * i);
      __m128 v0           = _mm_shuffle_ps(u0, u1, 0x44);
      __m128 v1           = _mm_shuffle_ps(u0, u1, 0xee);
      __m128 v2           = _mm_shuffle_ps(u2, u3, 0x44);
      __m128 v3           = _mm_shuffle_ps(u2, u3, 0xee);
      out[4 * i + 0].data = _mm_shuffle_ps(v0, v2, 0x88);
      out[4 * i + 1].data = _mm_shuffle_ps(v0, v2, 0xdd);
      out[4 * i + 2].data = _mm_shuffle_ps(v1, v3, 0x88);
      out[4 * i + 3].data = _mm_shuffle_ps(v1, v3, 0xdd);
    }
 
  // remainder loop for lengths not divisible by 4
  for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
    for (unsigned int v = 0; v < 4; ++v)
      out[i][v] = in[v][i];
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_load_and_transpose(const unsigned int         n_entries,
                              const float               *in,
                              const unsigned int        *offsets,
                              VectorizedArray<float, 4> *out)
{
  vectorized_load_and_transpose(
    n_entries,
    {{in + offsets[0], in + offsets[1], in + offsets[2], in + offsets[3]}},
    out);
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_transpose_and_store(const bool                       add_into,
                               const unsigned int               n_entries,
                               const VectorizedArray<float, 4> *in,
                               const std::array<float *, 4>    &out)
{
  const unsigned int n_chunks = n_entries / 4;
  for (unsigned int i = 0; i < n_chunks; ++i)
    {
      __m128 u0 = in[4 * i + 0].data;
      __m128 u1 = in[4 * i + 1].data;
      __m128 u2 = in[4 * i + 2].data;
      __m128 u3 = in[4 * i + 3].data;
      __m128 t0 = _mm_shuffle_ps(u0, u1, 0x44);
      __m128 t1 = _mm_shuffle_ps(u0, u1, 0xee);
      __m128 t2 = _mm_shuffle_ps(u2, u3, 0x44);
      __m128 t3 = _mm_shuffle_ps(u2, u3, 0xee);
      u0        = _mm_shuffle_ps(t0, t2, 0x88);
      u1        = _mm_shuffle_ps(t0, t2, 0xdd);
      u2        = _mm_shuffle_ps(t1, t3, 0x88);
      u3        = _mm_shuffle_ps(t1, t3, 0xdd);
 
      // Cannot use the same store instructions in both paths of the 'if'
      // because the compiler cannot know that there is no aliasing between
      // pointers
      if (add_into)
        {
          u0 = _mm_add_ps(_mm_loadu_ps(out[0] + 4 * i), u0);
          _mm_storeu_ps(out[0] + 4 * i, u0);
          u1 = _mm_add_ps(_mm_loadu_ps(out[1] + 4 * i), u1);
          _mm_storeu_ps(out[1] + 4 * i, u1);
          u2 = _mm_add_ps(_mm_loadu_ps(out[2] + 4 * i), u2);
          _mm_storeu_ps(out[2] + 4 * i, u2);
          u3 = _mm_add_ps(_mm_loadu_ps(out[3] + 4 * i), u3);
          _mm_storeu_ps(out[3] + 4 * i, u3);
        }
      else
        {
          _mm_storeu_ps(out[0] + 4 * i, u0);
          _mm_storeu_ps(out[1] + 4 * i, u1);
          _mm_storeu_ps(out[2] + 4 * i, u2);
          _mm_storeu_ps(out[3] + 4 * i, u3);
        }
    }
 
  // remainder loop of work that does not divide by 4
  if (add_into)
    for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
      for (unsigned int v = 0; v < 4; ++v)
        out[v][i] += in[i][v];
  else
    for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
      for (unsigned int v = 0; v < 4; ++v)
        out[v][i] = in[i][v];
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_transpose_and_store(const bool                       add_into,
                               const unsigned int               n_entries,
                               const VectorizedArray<float, 4> *in,
                               const unsigned int              *offsets,
                               float                           *out)
{
  vectorized_transpose_and_store(
    add_into,
    n_entries,
    in,
    {{out + offsets[0], out + offsets[1], out + offsets[2], out + offsets[3]}});
}
 
 
 
#  endif // if DEAL_II_VECTORIZATION_WIDTH_IN_BITS > 0 && defined(__SSE2__)
 
#  if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)

template <>
class VectorizedArray<double, 4>
  : public VectorizedArrayBase<VectorizedArray<double, 4>, 4>
{
public:
  using value_type = double;

  static constexpr bool is_implemented = true;

  VectorizedArray() = default;

  VectorizedArray(const double scalar)
  {
    this->operator=(scalar);
  }

  template <typename U>
  VectorizedArray(const std::initializer_list<U> &list)
    : VectorizedArrayBase<VectorizedArray<double, 4>, 4>(list)
  {}

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator=(const double x) &
  {
    data = _mm256_set1_pd(x);
    return *this;
  }

  VectorizedArray &
  operator=(const double scalar) && = delete;

  DEAL_II_ALWAYS_INLINE
  double &
  operator[](const unsigned int comp)
  {
    AssertIndexRange(comp, 4);
    return *(reinterpret_cast<double *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  const double &
  operator[](const unsigned int comp) const
  {
    AssertIndexRange(comp, 4);
    return *(reinterpret_cast<const double *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator+=(const VectorizedArray &vec)
  {
    // if the compiler supports vector arithmetic, we can simply use +=
    // operator on the given data type. this allows the compiler to combine
    // additions with multiplication (fused multiply-add) if those
    // instructions are available. Otherwise, we need to use the built-in
    // intrinsic command for __m256d
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data += vec.data;
#    else
    data = _mm256_add_pd(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator-=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data -= vec.data;
#    else
    data = _mm256_sub_pd(data, vec.data);
#    endif
    return *this;
  }
  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator*=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data *= vec.data;
#    else
    data = _mm256_mul_pd(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator/=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data /= vec.data;
#    else
    data = _mm256_div_pd(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  void
  load(const double *ptr)
  {
    data = _mm256_loadu_pd(ptr);
  }
 
  DEAL_II_ALWAYS_INLINE
  void
  load(const float *ptr)
  {
    data = _mm256_cvtps_pd(_mm_loadu_ps(ptr));
  }

  DEAL_II_ALWAYS_INLINE
  void
  store(double *ptr) const
  {
    _mm256_storeu_pd(ptr, data);
  }
 
  DEAL_II_ALWAYS_INLINE
  void
  store(float *ptr) const
  {
    _mm_storeu_ps(ptr, _mm256_cvtpd_ps(data));
  }

  DEAL_II_ALWAYS_INLINE
  void
  streaming_store(double *ptr) const
  {
    Assert(reinterpret_cast<std::size_t>(ptr) % 32 == 0,
           ExcMessage("Memory not aligned"));
    _mm256_stream_pd(ptr, data);
  }

  DEAL_II_ALWAYS_INLINE
  void
  gather(const double *base_ptr, const unsigned int *offsets)
  {
#    if defined(__AVX2__) && defined(DEAL_II_USE_VECTORIZATION_GATHER)
    // unfortunately, there does not appear to be a 128 bit integer load, so
    // do it by some reinterpret casts here. this is allowed because the Intel
    // API allows aliasing between different vector types.
    const __m128 index_val =
      _mm_loadu_ps(reinterpret_cast<const float *>(offsets));
    const __m128i index = *reinterpret_cast<const __m128i *>(&index_val);
 
    // work around a warning with gcc-12 about an uninitialized initial state
    // for gather by starting with a zero guess, even though all lanes will be
    // overwritten
    __m256d zero  = _mm256_setzero_pd();
    __m128i neq32 = _mm_andnot_si128(_mm_cmpeq_epi32(index, _mm_set1_epi32(-1)),
                                     _mm_set1_epi32(-1));
 
    __m256d mask = _mm256_castsi256_pd(_mm256_cvtepi32_epi64(neq32));
 
    data = _mm256_mask_i32gather_pd(zero, base_ptr, index, mask, 8);
#    else
    this->operator=(0.);
    for (unsigned int i = 0; i < 4; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        *(reinterpret_cast<double *>(&data) + i) = base_ptr[offsets[i]];
#    endif
  }

  DEAL_II_ALWAYS_INLINE
  void
  scatter(const unsigned int *offsets, double *base_ptr) const
  {
    // no scatter operation in AVX/AVX2
    for (unsigned int i = 0; i < 4; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        base_ptr[offsets[i]] = *(reinterpret_cast<const double *>(&data) + i);
  }

  double
  sum() const
  {
    VectorizedArray<double, 2> t1;
    t1.data = _mm_add_pd(this->get_lower(), this->get_upper());
    return t1.sum();
  }

  VectorizedArray
  get_floor() const
  {
    VectorizedArray res;
    res.data = _mm256_floor_pd(data);
    return res;
  }

  __m256d data;
 
private:
  DEAL_II_ALWAYS_INLINE
  __m128d
  get_lower() const
  {
    return _mm256_castpd256_pd128(data);
  }

  DEAL_II_ALWAYS_INLINE
  __m128d
  get_upper() const
  {
    return _mm256_extractf128_pd(data, 1);
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_sqrt() const
  {
    VectorizedArray res;
    res.data = _mm256_sqrt_pd(data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_abs() const
  {
    // to compute the absolute value, perform bitwise andnot with -0. This
    // will leave all value and exponent bits unchanged but force the sign
    // value to +.
    __m256d         mask = _mm256_set1_pd(-0.);
    VectorizedArray res;
    res.data = _mm256_andnot_pd(mask, data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_max(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = _mm256_max_pd(data, other.data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_min(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = _mm256_min_pd(data, other.data);
    return res;
  }
 
  // Make a few functions friends.
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::sqrt(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::abs(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::max(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::min(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
};
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_load_and_transpose(const unsigned int                   n_entries,
                              const std::array<const double *, 4> &in,
                              VectorizedArray<double, 4>          *out)
{
  const unsigned int n_chunks = n_entries / 4;
  const double      *in0      = in[0];
  const double      *in1      = in[1];
  const double      *in2      = in[2];
  const double      *in3      = in[3];
 
  for (unsigned int i = 0; i < n_chunks; ++i)
    {
      __m256d u0          = _mm256_loadu_pd(in0 + 4 * i);
      __m256d u1          = _mm256_loadu_pd(in1 + 4 * i);
      __m256d u2          = _mm256_loadu_pd(in2 + 4 * i);
      __m256d u3          = _mm256_loadu_pd(in3 + 4 * i);
      __m256d t0          = _mm256_permute2f128_pd(u0, u2, 0x20);
      __m256d t1          = _mm256_permute2f128_pd(u1, u3, 0x20);
      __m256d t2          = _mm256_permute2f128_pd(u0, u2, 0x31);
      __m256d t3          = _mm256_permute2f128_pd(u1, u3, 0x31);
      out[4 * i + 0].data = _mm256_unpacklo_pd(t0, t1);
      out[4 * i + 1].data = _mm256_unpackhi_pd(t0, t1);
      out[4 * i + 2].data = _mm256_unpacklo_pd(t2, t3);
      out[4 * i + 3].data = _mm256_unpackhi_pd(t2, t3);
    }
 
  // Specialized remainder code
  const unsigned int i         = 4 * n_chunks;
  const unsigned int remainder = n_entries - i;
  if (remainder == 3)
    {
      // Use masked load of 3 doubles (this avoids accessing in[x][3], which
      // might be an invalid address)
      __m256i mask{-1, -1, -1, 0};
      __m256d u0      = _mm256_maskload_pd(in[0] + i, mask);
      __m256d u1      = _mm256_maskload_pd(in[1] + i, mask);
      __m256d u2      = _mm256_maskload_pd(in[2] + i, mask);
      __m256d u3      = _mm256_maskload_pd(in[3] + i, mask);
      __m256d t0      = _mm256_permute2f128_pd(u0, u2, 0x20);
      __m256d t1      = _mm256_permute2f128_pd(u1, u3, 0x20);
      __m256d t2      = _mm256_permute2f128_pd(u0, u2, 0x31);
      __m256d t3      = _mm256_permute2f128_pd(u1, u3, 0x31);
      out[i + 0].data = _mm256_unpacklo_pd(t0, t1);
      out[i + 1].data = _mm256_unpackhi_pd(t0, t1);
      out[i + 2].data = _mm256_unpacklo_pd(t2, t3);
    }
  else if (remainder == 2)
    {
      // Use 128 bit loads
      __m256d t0, t1 = {};
      t0 = _mm256_insertf128_pd(t1, _mm_loadu_pd(in[0] + i), 0);
      t0 = _mm256_insertf128_pd(t0, _mm_loadu_pd(in[2] + i), 1);
      t1 = _mm256_insertf128_pd(t1, _mm_loadu_pd(in[1] + i), 0);
      t1 = _mm256_insertf128_pd(t1, _mm_loadu_pd(in[3] + i), 1);
 
      out[i + 0].data = _mm256_unpacklo_pd(t0, t1);
      out[i + 1].data = _mm256_unpackhi_pd(t0, t1);
    }
  else if (remainder == 1)
    for (unsigned int v = 0; v < 4; ++v)
      out[i][v] = in[v][i];
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_load_and_transpose(const unsigned int          n_entries,
                              const double               *in,
                              const unsigned int         *offsets,
                              VectorizedArray<double, 4> *out)
{
  vectorized_load_and_transpose(
    n_entries,
    {{in + offsets[0], in + offsets[1], in + offsets[2], in + offsets[3]}},
    out);
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_transpose_and_store(const bool                        add_into,
                               const unsigned int                n_entries,
                               const VectorizedArray<double, 4> *in,
                               const std::array<double *, 4>    &out)
{
  const unsigned int n_chunks = n_entries / 4;
  double            *out0     = out[0];
  double            *out1     = out[1];
  double            *out2     = out[2];
  double            *out3     = out[3];
  for (unsigned int i = 0; i < n_chunks; ++i)
    {
      __m256d u0   = in[4 * i + 0].data;
      __m256d u1   = in[4 * i + 1].data;
      __m256d u2   = in[4 * i + 2].data;
      __m256d u3   = in[4 * i + 3].data;
      __m256d t0   = _mm256_permute2f128_pd(u0, u2, 0x20);
      __m256d t1   = _mm256_permute2f128_pd(u1, u3, 0x20);
      __m256d t2   = _mm256_permute2f128_pd(u0, u2, 0x31);
      __m256d t3   = _mm256_permute2f128_pd(u1, u3, 0x31);
      __m256d res0 = _mm256_unpacklo_pd(t0, t1);
      __m256d res1 = _mm256_unpackhi_pd(t0, t1);
      __m256d res2 = _mm256_unpacklo_pd(t2, t3);
      __m256d res3 = _mm256_unpackhi_pd(t2, t3);
 
      // Cannot use the same store instructions in both paths of the 'if'
      // because the compiler cannot know that there is no aliasing between
      // pointers
      if (add_into)
        {
          res0 = _mm256_add_pd(_mm256_loadu_pd(out0 + 4 * i), res0);
          _mm256_storeu_pd(out0 + 4 * i, res0);
          res1 = _mm256_add_pd(_mm256_loadu_pd(out1 + 4 * i), res1);
          _mm256_storeu_pd(out1 + 4 * i, res1);
          res2 = _mm256_add_pd(_mm256_loadu_pd(out2 + 4 * i), res2);
          _mm256_storeu_pd(out2 + 4 * i, res2);
          res3 = _mm256_add_pd(_mm256_loadu_pd(out3 + 4 * i), res3);
          _mm256_storeu_pd(out3 + 4 * i, res3);
        }
      else
        {
          _mm256_storeu_pd(out0 + 4 * i, res0);
          _mm256_storeu_pd(out1 + 4 * i, res1);
          _mm256_storeu_pd(out2 + 4 * i, res2);
          _mm256_storeu_pd(out3 + 4 * i, res3);
        }
    }
 
  // remainder loop of work that does not divide by 4
  if (add_into)
    for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
      for (unsigned int v = 0; v < 4; ++v)
        out[v][i] += in[i][v];
  else
    for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
      for (unsigned int v = 0; v < 4; ++v)
        out[v][i] = in[i][v];
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_transpose_and_store(const bool                        add_into,
                               const unsigned int                n_entries,
                               const VectorizedArray<double, 4> *in,
                               const unsigned int               *offsets,
                               double                           *out)
{
  vectorized_transpose_and_store(
    add_into,
    n_entries,
    in,
    {{out + offsets[0], out + offsets[1], out + offsets[2], out + offsets[3]}});
}
 
 

template <>
class VectorizedArray<float, 8>
  : public VectorizedArrayBase<VectorizedArray<float, 8>, 8>
{
public:
  using value_type = float;

  static constexpr bool is_implemented = true;

  VectorizedArray() = default;

  VectorizedArray(const float scalar)
  {
    this->operator=(scalar);
  }

  template <typename U>
  VectorizedArray(const std::initializer_list<U> &list)
    : VectorizedArrayBase<VectorizedArray<float, 8>, 8>(list)
  {}

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator=(const float x) &
  {
    data = _mm256_set1_ps(x);
    return *this;
  }

  VectorizedArray &
  operator=(const float scalar) && = delete;

  DEAL_II_ALWAYS_INLINE
  float &
  operator[](const unsigned int comp)
  {
    AssertIndexRange(comp, 8);
    return *(reinterpret_cast<float *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  const float &
  operator[](const unsigned int comp) const
  {
    AssertIndexRange(comp, 8);
    return *(reinterpret_cast<const float *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator+=(const VectorizedArray &vec)
  {
    // if the compiler supports vector arithmetic, we can simply use +=
    // operator on the given data type. this allows the compiler to combine
    // additions with multiplication (fused multiply-add) if those
    // instructions are available. Otherwise, we need to use the built-in
    // intrinsic command for __m256d
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data += vec.data;
#    else
    data = _mm256_add_ps(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator-=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data -= vec.data;
#    else
    data = _mm256_sub_ps(data, vec.data);
#    endif
    return *this;
  }
  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator*=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data *= vec.data;
#    else
    data = _mm256_mul_ps(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator/=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data /= vec.data;
#    else
    data = _mm256_div_ps(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  void
  load(const float *ptr)
  {
    data = _mm256_loadu_ps(ptr);
  }

  DEAL_II_ALWAYS_INLINE
  void
  store(float *ptr) const
  {
    _mm256_storeu_ps(ptr, data);
  }

  DEAL_II_ALWAYS_INLINE
  void
  streaming_store(float *ptr) const
  {
    Assert(reinterpret_cast<std::size_t>(ptr) % 32 == 0,
           ExcMessage("Memory not aligned"));
    _mm256_stream_ps(ptr, data);
  }

  DEAL_II_ALWAYS_INLINE
  void
  gather(const float *base_ptr, const unsigned int *offsets)
  {
#    if defined(__AVX2__) && defined(DEAL_II_USE_VECTORIZATION_GATHER)
    // unfortunately, there does not appear to be a 256 bit integer load, so
    // do it by some reinterpret casts here. this is allowed because the Intel
    // API allows aliasing between different vector types.
    const __m256 index_val =
      _mm256_loadu_ps(reinterpret_cast<const float *>(offsets));
    const __m256i index = *reinterpret_cast<const __m256i *>(&index_val);
 
    // work around a warning with gcc-12 about an uninitialized initial state
    // for gather by starting with a zero guess, even though all lanes will be
    // overwritten
    __m256  zero         = _mm256_setzero_ps();
    __m256i invalid      = _mm256_set1_epi32(numbers::invalid_unsigned_int);
    __m256i inverse_mask = _mm256_cmpeq_epi32(index, invalid);
    __m256i mask         = _mm256_xor_si256(invalid, inverse_mask);
 
    data = _mm256_mask_i32gather_ps(
      zero, base_ptr, index, *reinterpret_cast<__m256 *>(&mask), 4);
#    else
    this->operator=(0.f);
    for (unsigned int i = 0; i < 8; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        *(reinterpret_cast<float *>(&data) + i) = base_ptr[offsets[i]];
#    endif
  }

  DEAL_II_ALWAYS_INLINE
  void
  scatter(const unsigned int *offsets, float *base_ptr) const
  {
    // no scatter operation in AVX/AVX2
    for (unsigned int i = 0; i < 8; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        base_ptr[offsets[i]] = *(reinterpret_cast<const float *>(&data) + i);
  }

  float
  sum() const
  {
    VectorizedArray<float, 4> t1;
    t1.data = _mm_add_ps(this->get_lower(), this->get_upper());
    return t1.sum();
  }

  VectorizedArray
  get_floor() const
  {
    VectorizedArray res;
    res.data = _mm256_floor_ps(data);
    return res;
  }

  __m256 data;
 
private:
  DEAL_II_ALWAYS_INLINE
  __m128
  get_lower() const
  {
    return _mm256_castps256_ps128(data);
  }

  DEAL_II_ALWAYS_INLINE
  __m128
  get_upper() const
  {
    return _mm256_extractf128_ps(data, 1);
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_sqrt() const
  {
    VectorizedArray res;
    res.data = _mm256_sqrt_ps(data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_abs() const
  {
    // to compute the absolute value, perform bitwise andnot with -0. This
    // will leave all value and exponent bits unchanged but force the sign
    // value to +.
    __m256          mask = _mm256_set1_ps(-0.f);
    VectorizedArray res;
    res.data = _mm256_andnot_ps(mask, data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_max(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = _mm256_max_ps(data, other.data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_min(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = _mm256_min_ps(data, other.data);
    return res;
  }
 
  // Make a few functions friends.
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::sqrt(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::abs(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::max(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::min(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
};
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_load_and_transpose(const unsigned int                  n_entries,
                              const std::array<const float *, 8> &in,
                              VectorizedArray<float, 8>          *out)
{
  const unsigned int n_chunks = n_entries / 4;
  for (unsigned int i = 0; i < n_chunks; ++i)
    {
      // To avoid warnings about uninitialized variables, need to initialize
      // one variable with zero before using it.
      __m256 t0, t1, t2, t3 = {};
      t0 = _mm256_insertf128_ps(t3, _mm_loadu_ps(in[0] + 4 * i), 0);
      t0 = _mm256_insertf128_ps(t0, _mm_loadu_ps(in[4] + 4 * i), 1);
      t1 = _mm256_insertf128_ps(t3, _mm_loadu_ps(in[1] + 4 * i), 0);
      t1 = _mm256_insertf128_ps(t1, _mm_loadu_ps(in[5] + 4 * i), 1);
      t2 = _mm256_insertf128_ps(t3, _mm_loadu_ps(in[2] + 4 * i), 0);
      t2 = _mm256_insertf128_ps(t2, _mm_loadu_ps(in[6] + 4 * i), 1);
      t3 = _mm256_insertf128_ps(t3, _mm_loadu_ps(in[3] + 4 * i), 0);
      t3 = _mm256_insertf128_ps(t3, _mm_loadu_ps(in[7] + 4 * i), 1);
 
      __m256 v0           = _mm256_shuffle_ps(t0, t1, 0x44);
      __m256 v1           = _mm256_shuffle_ps(t0, t1, 0xee);
      __m256 v2           = _mm256_shuffle_ps(t2, t3, 0x44);
      __m256 v3           = _mm256_shuffle_ps(t2, t3, 0xee);
      out[4 * i + 0].data = _mm256_shuffle_ps(v0, v2, 0x88);
      out[4 * i + 1].data = _mm256_shuffle_ps(v0, v2, 0xdd);
      out[4 * i + 2].data = _mm256_shuffle_ps(v1, v3, 0x88);
      out[4 * i + 3].data = _mm256_shuffle_ps(v1, v3, 0xdd);
    }
 
  // Specialized remainder code
  const unsigned int i         = 4 * n_chunks;
  const unsigned int remainder = n_entries - i;
  if (remainder == 3)
    {
      // Use masked load of 3 floats (this avoids accessing in[x][3], which
      // might be an invalid address)
      __m128i mask = _mm_setr_epi32(-1, -1, -1, 0);
      __m256  t0, t1, t2, t3 = {};
      t0        = _mm256_insertf128_ps(t3, _mm_maskload_ps(in[0] + i, mask), 0);
      t0        = _mm256_insertf128_ps(t0, _mm_maskload_ps(in[4] + i, mask), 1);
      t1        = _mm256_insertf128_ps(t3, _mm_maskload_ps(in[1] + i, mask), 0);
      t1        = _mm256_insertf128_ps(t1, _mm_maskload_ps(in[5] + i, mask), 1);
      t2        = _mm256_insertf128_ps(t3, _mm_maskload_ps(in[2] + i, mask), 0);
      t2        = _mm256_insertf128_ps(t2, _mm_maskload_ps(in[6] + i, mask), 1);
      t3        = _mm256_insertf128_ps(t3, _mm_maskload_ps(in[3] + i, mask), 0);
      t3        = _mm256_insertf128_ps(t3, _mm_maskload_ps(in[7] + i, mask), 1);
      __m256 v0 = _mm256_shuffle_ps(t0, t1, 0x44);
      __m256 v1 = _mm256_shuffle_ps(t0, t1, 0xee);
      __m256 v2 = _mm256_shuffle_ps(t2, t3, 0x44);
      __m256 v3 = _mm256_shuffle_ps(t2, t3, 0xee);
      out[i + 0].data = _mm256_shuffle_ps(v0, v2, 0x88);
      out[i + 1].data = _mm256_shuffle_ps(v0, v2, 0xdd);
      out[i + 2].data = _mm256_shuffle_ps(v1, v3, 0x88);
    }
  else if (remainder == 2)
    {
      // Use 64-bit integer loads
      __m128i t0, t1;
      t0 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[0] + i));
      t1 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[1] + i));
      __m256i v0 = _mm256_castsi128_si256(_mm_unpacklo_epi64(t0, t1));
      t0 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[4] + i));
      t1 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[5] + i));
      v0 = _mm256_inserti128_si256(v0, _mm_unpacklo_epi64(t0, t1), 1);
      t0 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[2] + i));
      t1 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[3] + i));
      __m256i v1 = _mm256_castsi128_si256(_mm_unpacklo_epi64(t0, t1));
      t0 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[6] + i));
      t1 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[7] + i));
      v1 = _mm256_inserti128_si256(v1, _mm_unpacklo_epi64(t0, t1), 1);
      out[i + 0].data = _mm256_shuffle_ps(_mm256_castsi256_ps(v0),
                                          _mm256_castsi256_ps(v1),
                                          0x88);
      out[i + 1].data = _mm256_shuffle_ps(_mm256_castsi256_ps(v0),
                                          _mm256_castsi256_ps(v1),
                                          0xdd);
    }
  else if (remainder == 1)
    for (unsigned int v = 0; v < 8; ++v)
      out[i][v] = in[v][i];
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_load_and_transpose(const unsigned int         n_entries,
                              const float               *in,
                              const unsigned int        *offsets,
                              VectorizedArray<float, 8> *out)
{
  vectorized_load_and_transpose(n_entries,
                                {{in + offsets[0],
                                  in + offsets[1],
                                  in + offsets[2],
                                  in + offsets[3],
                                  in + offsets[4],
                                  in + offsets[5],
                                  in + offsets[6],
                                  in + offsets[7]}},
                                out);
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_transpose_and_store(const bool                       add_into,
                               const unsigned int               n_entries,
                               const VectorizedArray<float, 8> *in,
                               const std::array<float *, 8>    &out)
{
  const unsigned int n_chunks = n_entries / 4;
  for (unsigned int i = 0; i < n_chunks; ++i)
    {
      __m256 u0   = in[4 * i + 0].data;
      __m256 u1   = in[4 * i + 1].data;
      __m256 u2   = in[4 * i + 2].data;
      __m256 u3   = in[4 * i + 3].data;
      __m256 t0   = _mm256_shuffle_ps(u0, u1, 0x44);
      __m256 t1   = _mm256_shuffle_ps(u0, u1, 0xee);
      __m256 t2   = _mm256_shuffle_ps(u2, u3, 0x44);
      __m256 t3   = _mm256_shuffle_ps(u2, u3, 0xee);
      u0          = _mm256_shuffle_ps(t0, t2, 0x88);
      u1          = _mm256_shuffle_ps(t0, t2, 0xdd);
      u2          = _mm256_shuffle_ps(t1, t3, 0x88);
      u3          = _mm256_shuffle_ps(t1, t3, 0xdd);
      __m128 res0 = _mm256_extractf128_ps(u0, 0);
      __m128 res4 = _mm256_extractf128_ps(u0, 1);
      __m128 res1 = _mm256_extractf128_ps(u1, 0);
      __m128 res5 = _mm256_extractf128_ps(u1, 1);
      __m128 res2 = _mm256_extractf128_ps(u2, 0);
      __m128 res6 = _mm256_extractf128_ps(u2, 1);
      __m128 res3 = _mm256_extractf128_ps(u3, 0);
      __m128 res7 = _mm256_extractf128_ps(u3, 1);
 
      // Cannot use the same store instructions in both paths of the 'if'
      // because the compiler cannot know that there is no aliasing between
      // pointers
      if (add_into)
        {
          res0 = _mm_add_ps(_mm_loadu_ps(out[0] + 4 * i), res0);
          _mm_storeu_ps(out[0] + 4 * i, res0);
          res1 = _mm_add_ps(_mm_loadu_ps(out[1] + 4 * i), res1);
          _mm_storeu_ps(out[1] + 4 * i, res1);
          res2 = _mm_add_ps(_mm_loadu_ps(out[2] + 4 * i), res2);
          _mm_storeu_ps(out[2] + 4 * i, res2);
          res3 = _mm_add_ps(_mm_loadu_ps(out[3] + 4 * i), res3);
          _mm_storeu_ps(out[3] + 4 * i, res3);
          res4 = _mm_add_ps(_mm_loadu_ps(out[4] + 4 * i), res4);
          _mm_storeu_ps(out[4] + 4 * i, res4);
          res5 = _mm_add_ps(_mm_loadu_ps(out[5] + 4 * i), res5);
          _mm_storeu_ps(out[5] + 4 * i, res5);
          res6 = _mm_add_ps(_mm_loadu_ps(out[6] + 4 * i), res6);
          _mm_storeu_ps(out[6] + 4 * i, res6);
          res7 = _mm_add_ps(_mm_loadu_ps(out[7] + 4 * i), res7);
          _mm_storeu_ps(out[7] + 4 * i, res7);
        }
      else
        {
          _mm_storeu_ps(out[0] + 4 * i, res0);
          _mm_storeu_ps(out[1] + 4 * i, res1);
          _mm_storeu_ps(out[2] + 4 * i, res2);
          _mm_storeu_ps(out[3] + 4 * i, res3);
          _mm_storeu_ps(out[4] + 4 * i, res4);
          _mm_storeu_ps(out[5] + 4 * i, res5);
          _mm_storeu_ps(out[6] + 4 * i, res6);
          _mm_storeu_ps(out[7] + 4 * i, res7);
        }
    }
 
  if (add_into)
    for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
      for (unsigned int v = 0; v < 8; ++v)
        out[v][i] += in[i][v];
  else
    for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
      for (unsigned int v = 0; v < 8; ++v)
        out[v][i] = in[i][v];
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_transpose_and_store(const bool                       add_into,
                               const unsigned int               n_entries,
                               const VectorizedArray<float, 8> *in,
                               const unsigned int              *offsets,
                               float                           *out)
{
  vectorized_transpose_and_store(add_into,
                                 n_entries,
                                 in,
                                 {{out + offsets[0],
                                   out + offsets[1],
                                   out + offsets[2],
                                   out + offsets[3],
                                   out + offsets[4],
                                   out + offsets[5],
                                   out + offsets[6],
                                   out + offsets[7]}});
}
 
#  endif
 
// for safety, also check that __AVX512F__ is defined in case the user manually
// set some conflicting compile flags which prevent compilation
 
#  if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && defined(__AVX512F__)

template <>
class VectorizedArray<double, 8>
  : public VectorizedArrayBase<VectorizedArray<double, 8>, 8>
{
public:
  using value_type = double;

  static constexpr bool is_implemented = true;

  VectorizedArray() = default;

  VectorizedArray(const double scalar)
  {
    this->operator=(scalar);
  }

  template <typename U>
  VectorizedArray(const std::initializer_list<U> &list)
    : VectorizedArrayBase<VectorizedArray<double, 8>, 8>(list)
  {}

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator=(const double x) &
  {
    data = _mm512_set1_pd(x);
    return *this;
  }
 

  VectorizedArray &
  operator=(const double scalar) && = delete;

  DEAL_II_ALWAYS_INLINE
  double &
  operator[](const unsigned int comp)
  {
    AssertIndexRange(comp, 8);
    return *(reinterpret_cast<double *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  const double &
  operator[](const unsigned int comp) const
  {
    AssertIndexRange(comp, 8);
    return *(reinterpret_cast<const double *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator+=(const VectorizedArray &vec)
  {
    // if the compiler supports vector arithmetic, we can simply use +=
    // operator on the given data type. this allows the compiler to combine
    // additions with multiplication (fused multiply-add) if those
    // instructions are available. Otherwise, we need to use the built-in
    // intrinsic command for __m512d
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data += vec.data;
#    else
    data = _mm512_add_pd(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator-=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data -= vec.data;
#    else
    data = _mm512_sub_pd(data, vec.data);
#    endif
    return *this;
  }
  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator*=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data *= vec.data;
#    else
    data = _mm512_mul_pd(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator/=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data /= vec.data;
#    else
    data = _mm512_div_pd(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  void
  load(const double *ptr)
  {
    data = _mm512_loadu_pd(ptr);
  }
 
  DEAL_II_ALWAYS_INLINE
  void
  load(const float *ptr)
  {
    data = _mm512_cvtps_pd(_mm256_loadu_ps(ptr));
  }

  DEAL_II_ALWAYS_INLINE
  void
  store(double *ptr) const
  {
    _mm512_storeu_pd(ptr, data);
  }
 
  DEAL_II_ALWAYS_INLINE
  void
  store(float *ptr) const
  {
    _mm256_storeu_ps(ptr, _mm512_cvtpd_ps(data));
  }

  DEAL_II_ALWAYS_INLINE
  void
  streaming_store(double *ptr) const
  {
    Assert(reinterpret_cast<std::size_t>(ptr) % 64 == 0,
           ExcMessage("Memory not aligned"));
    _mm512_stream_pd(ptr, data);
  }

  DEAL_II_ALWAYS_INLINE
  void
  gather(const double *base_ptr, const unsigned int *offsets)
  {
#    ifdef DEAL_II_USE_VECTORIZATION_GATHER
    // unfortunately, there does not appear to be a 256 bit integer load, so
    // do it by some reinterpret casts here. this is allowed because the Intel
    // API allows aliasing between different vector types.
    const __m256 index_val =
      _mm256_loadu_ps(reinterpret_cast<const float *>(offsets));
    const __m256i index = *reinterpret_cast<const __m256i *>(&index_val);
 
    // work around a warning with gcc-12 about an uninitialized initial state
    // for gather by starting with a zero guess, even though all lanes will be
    // overwritten
    __m512d       zero    = {};
    const __m256i invalid = _mm256_set1_epi32(numbers::invalid_unsigned_int);
    const __m256i equal   = _mm256_cmpeq_epi32(invalid, index);
    __mmask8      mask =
      static_cast<__mmask8>(~(_mm256_movemask_ps(_mm256_castsi256_ps(equal))));
 
    data = _mm512_mask_i32gather_pd(zero, mask, index, base_ptr, 8);
#    else
    this->operator=(0.);
    for (unsigned int i = 0; i < 8; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        *(reinterpret_cast<double *>(&data) + i) = base_ptr[offsets[i]];
#    endif
  }

  DEAL_II_ALWAYS_INLINE
  void
  scatter(const unsigned int *offsets, double *base_ptr) const
  {
#    ifdef DEAL_II_USE_VECTORIZATION_GATHER
    for (unsigned int i = 0; i < 8; ++i)
      for (unsigned int j = i + 1; j < 8; ++j)
        Assert(offsets[i] == numbers::invalid_unsigned_int ||
                 offsets[i] != offsets[j],
               ExcMessage("Result of scatter undefined if two offset elements"
                          " point to the same position"));
 
    // unfortunately, there does not appear to be a 256 bit integer load, so
    // do it by some reinterpret casts here. this is allowed because the Intel
    // API allows aliasing between different vector types.
    const __m256 index_val =
      _mm256_loadu_ps(reinterpret_cast<const float *>(offsets));
    const __m256i index   = *reinterpret_cast<const __m256i *>(&index_val);
    const __m256i invalid = _mm256_set1_epi32(numbers::invalid_unsigned_int);
    const __m256i equal   = _mm256_cmpeq_epi32(invalid, index);
    __mmask8      mask =
      static_cast<__mmask8>(~(_mm256_movemask_ps(_mm256_castsi256_ps(equal))));
    _mm512_mask_i32scatter_pd(base_ptr, mask, index, data, 8);
 
#    else
    for (unsigned int i = 0; i < 8; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        base_ptr[offsets[i]] = *(reinterpret_cast<const double *>(&data) + i);
#    endif
  }

  double
  sum() const
  {
    VectorizedArray<double, 4> t1;
    t1.data = _mm256_add_pd(this->get_lower(), this->get_upper());
    return t1.sum();
  }

  VectorizedArray
  get_floor() const
  {
    VectorizedArray res;
    res.data = _mm512_roundscale_pd(data, _MM_FROUND_TO_NEG_INF);
    return res;
  }

  __m512d data;
 
private:
  DEAL_II_ALWAYS_INLINE
  __m256d
  get_lower() const
  {
    return _mm512_castpd512_pd256(data);
  }

  DEAL_II_ALWAYS_INLINE
  __m256d
  get_upper() const
  {
    return _mm512_extractf64x4_pd(data, 1);
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_sqrt() const
  {
    VectorizedArray res;
    res.data = _mm512_sqrt_pd(data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_abs() const
  {
    // to compute the absolute value, perform bitwise andnot with -0. This
    // will leave all value and exponent bits unchanged but force the sign
    // value to +. Since there is no andnot for AVX512, we interpret the data
    // as 64 bit integers and do the andnot on those types (note that andnot
    // is a bitwise operation so the data type does not matter)
    __m512d         mask = _mm512_set1_pd(-0.);
    VectorizedArray res;
    res.data = reinterpret_cast<__m512d>(
      _mm512_andnot_epi64(reinterpret_cast<__m512i>(mask),
                          reinterpret_cast<__m512i>(data)));
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_max(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = _mm512_max_pd(data, other.data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_min(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = _mm512_min_pd(data, other.data);
    return res;
  }
 
  // Make a few functions friends.
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::sqrt(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::abs(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::max(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::min(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
};
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_load_and_transpose(const unsigned int                   n_entries,
                              const std::array<const double *, 8> &in,
                              VectorizedArray<double, 8>          *out)
{
  // do not do full transpose because the code is long and will most
  // likely not pay off because many processors have two load units
  // (for the top 8 instructions) but only 1 permute unit (for the 8
  // shuffle/unpack instructions). rather start the transposition on the
  // vectorized array of half the size with 256 bits
  const unsigned int n_chunks = n_entries / 4;
  for (unsigned int i = 0; i < n_chunks; ++i)
    {
      __m512d t0, t1, t2, t3 = {};
 
      t0 = _mm512_insertf64x4(t3, _mm256_loadu_pd(in[0] + 4 * i), 0);
      t0 = _mm512_insertf64x4(t0, _mm256_loadu_pd(in[2] + 4 * i), 1);
      t1 = _mm512_insertf64x4(t3, _mm256_loadu_pd(in[1] + 4 * i), 0);
      t1 = _mm512_insertf64x4(t1, _mm256_loadu_pd(in[3] + 4 * i), 1);
      t2 = _mm512_insertf64x4(t3, _mm256_loadu_pd(in[4] + 4 * i), 0);
      t2 = _mm512_insertf64x4(t2, _mm256_loadu_pd(in[6] + 4 * i), 1);
      t3 = _mm512_insertf64x4(t3, _mm256_loadu_pd(in[5] + 4 * i), 0);
      t3 = _mm512_insertf64x4(t3, _mm256_loadu_pd(in[7] + 4 * i), 1);
 
      __m512d v0          = _mm512_shuffle_f64x2(t0, t2, 0x88);
      __m512d v1          = _mm512_shuffle_f64x2(t0, t2, 0xdd);
      __m512d v2          = _mm512_shuffle_f64x2(t1, t3, 0x88);
      __m512d v3          = _mm512_shuffle_f64x2(t1, t3, 0xdd);
      out[4 * i + 0].data = _mm512_unpacklo_pd(v0, v2);
      out[4 * i + 1].data = _mm512_unpackhi_pd(v0, v2);
      out[4 * i + 2].data = _mm512_unpacklo_pd(v1, v3);
      out[4 * i + 3].data = _mm512_unpackhi_pd(v1, v3);
    }
 
  // Specialized remainder code
  const unsigned int i         = 4 * n_chunks;
  const unsigned int remainder = n_entries - i;
  if (remainder == 3)
    {
      // Use masked loads to access the first 3 doubles along each pointer
      // without touching the possibly invalid addressed at in[x][3].
      __m256i mask{-1, -1, -1, 0};
      __m512d t0, t1, t2, t3 = {};
 
      t0 = _mm512_insertf64x4(t3, _mm256_maskload_pd(in[0] + i, mask), 0);
      t0 = _mm512_insertf64x4(t0, _mm256_maskload_pd(in[2] + i, mask), 1);
      t1 = _mm512_insertf64x4(t3, _mm256_maskload_pd(in[1] + i, mask), 0);
      t1 = _mm512_insertf64x4(t1, _mm256_maskload_pd(in[3] + i, mask), 1);
      t2 = _mm512_insertf64x4(t3, _mm256_maskload_pd(in[4] + i, mask), 0);
      t2 = _mm512_insertf64x4(t2, _mm256_maskload_pd(in[6] + i, mask), 1);
      t3 = _mm512_insertf64x4(t3, _mm256_maskload_pd(in[5] + i, mask), 0);
      t3 = _mm512_insertf64x4(t3, _mm256_maskload_pd(in[7] + i, mask), 1);
 
      __m512d v0      = _mm512_shuffle_f64x2(t0, t2, 0x88);
      __m512d v1      = _mm512_shuffle_f64x2(t0, t2, 0xdd);
      __m512d v2      = _mm512_shuffle_f64x2(t1, t3, 0x88);
      __m512d v3      = _mm512_shuffle_f64x2(t1, t3, 0xdd);
      out[i + 0].data = _mm512_unpacklo_pd(v0, v2);
      out[i + 1].data = _mm512_unpackhi_pd(v0, v2);
      out[i + 2].data = _mm512_unpacklo_pd(v1, v3);
    }
  else if (remainder == 2)
    {
      // Use 128-bit loads of 2 doubles
      __m256d t0, t1, t2, t3 = {};
      t0 = _mm256_insertf128_pd(t3, _mm_loadu_pd(in[0] + i), 0);
      t0 = _mm256_insertf128_pd(t0, _mm_loadu_pd(in[2] + i), 1);
      t1 = _mm256_insertf128_pd(t3, _mm_loadu_pd(in[1] + i), 0);
      t1 = _mm256_insertf128_pd(t1, _mm_loadu_pd(in[3] + i), 1);
      t2 = _mm256_insertf128_pd(t3, _mm_loadu_pd(in[4] + i), 0);
      t2 = _mm256_insertf128_pd(t2, _mm_loadu_pd(in[6] + i), 1);
      t3 = _mm256_insertf128_pd(t3, _mm_loadu_pd(in[5] + i), 0);
      t3 = _mm256_insertf128_pd(t3, _mm_loadu_pd(in[7] + i), 1);
 
      __m512d v0, v1 = {};
      v0              = _mm512_insertf64x4(v1, t0, 0);
      v0              = _mm512_insertf64x4(v0, t2, 1);
      v1              = _mm512_insertf64x4(v1, t1, 0);
      v1              = _mm512_insertf64x4(v1, t3, 1);
      out[i + 0].data = _mm512_unpacklo_pd(v0, v1);
      out[i + 1].data = _mm512_unpackhi_pd(v0, v1);
    }
  else if (remainder == 1)
    for (unsigned int v = 0; v < 8; ++v)
      out[i][v] = in[v][i];
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_load_and_transpose(const unsigned int          n_entries,
                              const double               *in,
                              const unsigned int         *offsets,
                              VectorizedArray<double, 8> *out)
{
  vectorized_load_and_transpose(n_entries,
                                {{in + offsets[0],
                                  in + offsets[1],
                                  in + offsets[2],
                                  in + offsets[3],
                                  in + offsets[4],
                                  in + offsets[5],
                                  in + offsets[6],
                                  in + offsets[7]}},
                                out);
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_transpose_and_store(const bool                        add_into,
                               const unsigned int                n_entries,
                               const VectorizedArray<double, 8> *in,
                               const std::array<double *, 8>    &out)
{
  // as for the load, we split the store operations into 256 bit units to
  // better balance between code size, shuffle instructions, and stores
  const unsigned int n_chunks = n_entries / 4;
  __m512i mask1 = _mm512_set_epi64(0xd, 0xc, 0x5, 0x4, 0x9, 0x8, 0x1, 0x0);
  __m512i mask2 = _mm512_set_epi64(0xf, 0xe, 0x7, 0x6, 0xb, 0xa, 0x3, 0x2);
  for (unsigned int i = 0; i < n_chunks; ++i)
    {
      __m512d t0   = _mm512_unpacklo_pd(in[i * 4].data, in[i * 4 + 1].data);
      __m512d t1   = _mm512_unpackhi_pd(in[i * 4].data, in[i * 4 + 1].data);
      __m512d t2   = _mm512_unpacklo_pd(in[i * 4 + 2].data, in[i * 4 + 3].data);
      __m512d t3   = _mm512_unpackhi_pd(in[i * 4 + 2].data, in[i * 4 + 3].data);
      __m512d v0   = _mm512_permutex2var_pd(t0, mask1, t2);
      __m512d v1   = _mm512_permutex2var_pd(t0, mask2, t2);
      __m512d v2   = _mm512_permutex2var_pd(t1, mask1, t3);
      __m512d v3   = _mm512_permutex2var_pd(t1, mask2, t3);
      __m256d res0 = _mm512_extractf64x4_pd(v0, 0);
      __m256d res4 = _mm512_extractf64x4_pd(v0, 1);
      __m256d res1 = _mm512_extractf64x4_pd(v2, 0);
      __m256d res5 = _mm512_extractf64x4_pd(v2, 1);
      __m256d res2 = _mm512_extractf64x4_pd(v1, 0);
      __m256d res6 = _mm512_extractf64x4_pd(v1, 1);
      __m256d res3 = _mm512_extractf64x4_pd(v3, 0);
      __m256d res7 = _mm512_extractf64x4_pd(v3, 1);
 
      // Cannot use the same store instructions in both paths of the 'if'
      // because the compiler cannot know that there is no aliasing
      // between pointers
      if (add_into)
        {
          res0 = _mm256_add_pd(_mm256_loadu_pd(out[0] + 4 * i), res0);
          _mm256_storeu_pd(out[0] + 4 * i, res0);
          res1 = _mm256_add_pd(_mm256_loadu_pd(out[1] + 4 * i), res1);
          _mm256_storeu_pd(out[1] + 4 * i, res1);
          res2 = _mm256_add_pd(_mm256_loadu_pd(out[2] + 4 * i), res2);
          _mm256_storeu_pd(out[2] + 4 * i, res2);
          res3 = _mm256_add_pd(_mm256_loadu_pd(out[3] + 4 * i), res3);
          _mm256_storeu_pd(out[3] + 4 * i, res3);
          res4 = _mm256_add_pd(_mm256_loadu_pd(out[4] + 4 * i), res4);
          _mm256_storeu_pd(out[4] + 4 * i, res4);
          res5 = _mm256_add_pd(_mm256_loadu_pd(out[5] + 4 * i), res5);
          _mm256_storeu_pd(out[5] + 4 * i, res5);
          res6 = _mm256_add_pd(_mm256_loadu_pd(out[6] + 4 * i), res6);
          _mm256_storeu_pd(out[6] + 4 * i, res6);
          res7 = _mm256_add_pd(_mm256_loadu_pd(out[7] + 4 * i), res7);
          _mm256_storeu_pd(out[7] + 4 * i, res7);
        }
      else
        {
          _mm256_storeu_pd(out[0] + 4 * i, res0);
          _mm256_storeu_pd(out[1] + 4 * i, res1);
          _mm256_storeu_pd(out[2] + 4 * i, res2);
          _mm256_storeu_pd(out[3] + 4 * i, res3);
          _mm256_storeu_pd(out[4] + 4 * i, res4);
          _mm256_storeu_pd(out[5] + 4 * i, res5);
          _mm256_storeu_pd(out[6] + 4 * i, res6);
          _mm256_storeu_pd(out[7] + 4 * i, res7);
        }
    }
 
  if (add_into)
    for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
      for (unsigned int v = 0; v < 8; ++v)
        out[v][i] += in[i][v];
  else
    for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
      for (unsigned int v = 0; v < 8; ++v)
        out[v][i] = in[i][v];
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_transpose_and_store(const bool                        add_into,
                               const unsigned int                n_entries,
                               const VectorizedArray<double, 8> *in,
                               const unsigned int               *offsets,
                               double                           *out)
{
  vectorized_transpose_and_store(add_into,
                                 n_entries,
                                 in,
                                 {{out + offsets[0],
                                   out + offsets[1],
                                   out + offsets[2],
                                   out + offsets[3],
                                   out + offsets[4],
                                   out + offsets[5],
                                   out + offsets[6],
                                   out + offsets[7]}});
}
 
 

template <>
class VectorizedArray<float, 16>
  : public VectorizedArrayBase<VectorizedArray<float, 16>, 16>
{
public:
  using value_type = float;

  static constexpr bool is_implemented = true;

  VectorizedArray() = default;

  VectorizedArray(const float scalar)
  {
    this->operator=(scalar);
  }

  template <typename U>
  VectorizedArray(const std::initializer_list<U> &list)
    : VectorizedArrayBase<VectorizedArray<float, 16>, 16>(list)
  {}

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator=(const float x) &
  {
    data = _mm512_set1_ps(x);
    return *this;
  }

  VectorizedArray &
  operator=(const float scalar) && = delete;

  DEAL_II_ALWAYS_INLINE
  float &
  operator[](const unsigned int comp)
  {
    AssertIndexRange(comp, 16);
    return *(reinterpret_cast<float *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  const float &
  operator[](const unsigned int comp) const
  {
    AssertIndexRange(comp, 16);
    return *(reinterpret_cast<const float *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator+=(const VectorizedArray &vec)
  {
    // if the compiler supports vector arithmetic, we can simply use +=
    // operator on the given data type. this allows the compiler to combine
    // additions with multiplication (fused multiply-add) if those
    // instructions are available. Otherwise, we need to use the built-in
    // intrinsic command for __m512d
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data += vec.data;
#    else
    data = _mm512_add_ps(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator-=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data -= vec.data;
#    else
    data = _mm512_sub_ps(data, vec.data);
#    endif
    return *this;
  }
  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator*=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data *= vec.data;
#    else
    data = _mm512_mul_ps(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator/=(const VectorizedArray &vec)
  {
#    ifdef DEAL_II_COMPILER_USE_VECTOR_ARITHMETICS
    data /= vec.data;
#    else
    data = _mm512_div_ps(data, vec.data);
#    endif
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  void
  load(const float *ptr)
  {
    data = _mm512_loadu_ps(ptr);
  }

  DEAL_II_ALWAYS_INLINE
  void
  store(float *ptr) const
  {
    _mm512_storeu_ps(ptr, data);
  }

  DEAL_II_ALWAYS_INLINE
  void
  streaming_store(float *ptr) const
  {
    Assert(reinterpret_cast<std::size_t>(ptr) % 64 == 0,
           ExcMessage("Memory not aligned"));
    _mm512_stream_ps(ptr, data);
  }

  DEAL_II_ALWAYS_INLINE
  void
  gather(const float *base_ptr, const unsigned int *offsets)
  {
#    ifdef DEAL_II_USE_VECTORIZATION_GATHER
    // unfortunately, there does not appear to be a 512 bit integer load, so
    // do it by some reinterpret casts here. this is allowed because the Intel
    // API allows aliasing between different vector types.
    const __m512 index_val =
      _mm512_loadu_ps(reinterpret_cast<const float *>(offsets));
    const __m512i index = *reinterpret_cast<const __m512i *>(&index_val);
 
    // work around a warning with gcc-12 about an uninitialized initial state
    // for gather by starting with a zero guess, even though all lanes will be
    // overwritten
    __m512        zero    = {};
    const __m512i invalid = _mm512_set1_epi32(numbers::invalid_unsigned_int);
    __mmask16     mask    = _mm512_cmpneq_epu32_mask(invalid, index);
 
    data = _mm512_mask_i32gather_ps(zero, mask, index, base_ptr, 4);
#    else
    this->operator=(0.f);
    for (unsigned int i = 0; i < 16; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        *(reinterpret_cast<float *>(&data) + i) = base_ptr[offsets[i]];
#    endif
  }

  DEAL_II_ALWAYS_INLINE
  void
  scatter(const unsigned int *offsets, float *base_ptr) const
  {
#    ifdef DEAL_II_USE_VECTORIZATION_GATHER
    for (unsigned int i = 0; i < 16; ++i)
      for (unsigned int j = i + 1; j < 16; ++j)
        Assert(offsets[i] == numbers::invalid_unsigned_int ||
                 offsets[i] != offsets[j],
               ExcMessage("Result of scatter undefined if two offset elements"
                          " point to the same position"));
 
    // unfortunately, there does not appear to be a 512 bit integer load, so
    // do it by some reinterpret casts here. this is allowed because the Intel
    // API allows aliasing between different vector types.
    const __m512 index_val =
      _mm512_loadu_ps(reinterpret_cast<const float *>(offsets));
    const __m512i index   = *reinterpret_cast<const __m512i *>(&index_val);
    const __m512i invalid = _mm512_set1_epi32(numbers::invalid_unsigned_int);
    __mmask16     mask    = _mm512_cmpneq_epu32_mask(invalid, index);
    _mm512_mask_i32scatter_ps(base_ptr, mask, index, data, 4);
#    else
    for (unsigned int i = 0; i < 16; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        base_ptr[offsets[i]] = *(reinterpret_cast<const float *>(&data) + i);
#    endif
  }

  float
  sum() const
  {
    VectorizedArray<float, 8> t1;
    t1.data = _mm256_add_ps(this->get_lower(), this->get_upper());
    return t1.sum();
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_floor() const
  {
    VectorizedArray res;
    res.data = _mm512_roundscale_ps(data, _MM_FROUND_TO_NEG_INF);
    return res;
  }

  __m512 data;
 
private:
  DEAL_II_ALWAYS_INLINE
  __m256
  get_lower() const
  {
    return _mm512_castps512_ps256(data);
  }

  DEAL_II_ALWAYS_INLINE
  __m256
  get_upper() const
  {
    return _mm256_castpd_ps(_mm512_extractf64x4_pd(_mm512_castps_pd(data), 1));
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_sqrt() const
  {
    VectorizedArray res;
    res.data = _mm512_sqrt_ps(data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_abs() const
  {
    // to compute the absolute value, perform bitwise andnot with -0. This
    // will leave all value and exponent bits unchanged but force the sign
    // value to +. Since there is no andnot for AVX512, we interpret the data
    // as 32 bit integers and do the andnot on those types (note that andnot
    // is a bitwise operation so the data type does not matter)
    __m512          mask = _mm512_set1_ps(-0.f);
    VectorizedArray res;
    res.data = reinterpret_cast<__m512>(
      _mm512_andnot_epi32(reinterpret_cast<__m512i>(mask),
                          reinterpret_cast<__m512i>(data)));
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_max(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = _mm512_max_ps(data, other.data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_min(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = _mm512_min_ps(data, other.data);
    return res;
  }
 
  // Make a few functions friends.
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::sqrt(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::abs(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::max(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::min(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
};
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_load_and_transpose(const unsigned int                   n_entries,
                              const std::array<const float *, 16> &in,
                              VectorizedArray<float, 16>          *out)
{
  // Similar to the double case, we perform the work on smaller entities. In
  // this case, we start from 128 bit arrays and insert them into a full 512
  // bit index. This reduces the code size and register pressure because we do
  // shuffles on 4 numbers rather than 16.
  const unsigned int n_chunks = n_entries / 4;
 
  for (unsigned int i = 0; i < n_chunks; ++i)
    {
      __m512 t0, t1, t2, t3 = {};
 
      t0 = _mm512_insertf32x4(t3, _mm_loadu_ps(in[0] + 4 * i), 0);
      t0 = _mm512_insertf32x4(t0, _mm_loadu_ps(in[4] + 4 * i), 1);
      t0 = _mm512_insertf32x4(t0, _mm_loadu_ps(in[8] + 4 * i), 2);
      t0 = _mm512_insertf32x4(t0, _mm_loadu_ps(in[12] + 4 * i), 3);
      t1 = _mm512_insertf32x4(t3, _mm_loadu_ps(in[1] + 4 * i), 0);
      t1 = _mm512_insertf32x4(t1, _mm_loadu_ps(in[5] + 4 * i), 1);
      t1 = _mm512_insertf32x4(t1, _mm_loadu_ps(in[9] + 4 * i), 2);
      t1 = _mm512_insertf32x4(t1, _mm_loadu_ps(in[13] + 4 * i), 3);
      t2 = _mm512_insertf32x4(t3, _mm_loadu_ps(in[2] + 4 * i), 0);
      t2 = _mm512_insertf32x4(t2, _mm_loadu_ps(in[6] + 4 * i), 1);
      t2 = _mm512_insertf32x4(t2, _mm_loadu_ps(in[10] + 4 * i), 2);
      t2 = _mm512_insertf32x4(t2, _mm_loadu_ps(in[14] + 4 * i), 3);
      t3 = _mm512_insertf32x4(t3, _mm_loadu_ps(in[3] + 4 * i), 0);
      t3 = _mm512_insertf32x4(t3, _mm_loadu_ps(in[7] + 4 * i), 1);
      t3 = _mm512_insertf32x4(t3, _mm_loadu_ps(in[11] + 4 * i), 2);
      t3 = _mm512_insertf32x4(t3, _mm_loadu_ps(in[15] + 4 * i), 3);
 
      __m512 v0 = _mm512_shuffle_ps(t0, t1, 0x44);
      __m512 v1 = _mm512_shuffle_ps(t0, t1, 0xee);
      __m512 v2 = _mm512_shuffle_ps(t2, t3, 0x44);
      __m512 v3 = _mm512_shuffle_ps(t2, t3, 0xee);
 
      out[4 * i + 0].data = _mm512_shuffle_ps(v0, v2, 0x88);
      out[4 * i + 1].data = _mm512_shuffle_ps(v0, v2, 0xdd);
      out[4 * i + 2].data = _mm512_shuffle_ps(v1, v3, 0x88);
      out[4 * i + 3].data = _mm512_shuffle_ps(v1, v3, 0xdd);
    }
 
  // Specialized remainder code
  const unsigned int i         = 4 * n_chunks;
  const unsigned int remainder = n_entries - i;
  if (remainder == 3)
    {
      __m512 t0, t1, t2, t3 = {};
 
      // Use masked load of 3 floats (this avoids accessing in[x][3], which
      // might be an invalid address)
      __m128i mask = _mm_setr_epi32(-1, -1, -1, 0);
      t0        = _mm512_insertf32x4(t3, _mm_maskload_ps(in[0] + i, mask), 0);
      t0        = _mm512_insertf32x4(t0, _mm_maskload_ps(in[4] + i, mask), 1);
      t0        = _mm512_insertf32x4(t0, _mm_maskload_ps(in[8] + i, mask), 2);
      t0        = _mm512_insertf32x4(t0, _mm_maskload_ps(in[12] + i, mask), 3);
      t1        = _mm512_insertf32x4(t3, _mm_maskload_ps(in[1] + i, mask), 0);
      t1        = _mm512_insertf32x4(t1, _mm_maskload_ps(in[5] + i, mask), 1);
      t1        = _mm512_insertf32x4(t1, _mm_maskload_ps(in[9] + i, mask), 2);
      t1        = _mm512_insertf32x4(t1, _mm_maskload_ps(in[13] + i, mask), 3);
      t2        = _mm512_insertf32x4(t3, _mm_maskload_ps(in[2] + i, mask), 0);
      t2        = _mm512_insertf32x4(t2, _mm_maskload_ps(in[6] + i, mask), 1);
      t2        = _mm512_insertf32x4(t2, _mm_maskload_ps(in[10] + i, mask), 2);
      t2        = _mm512_insertf32x4(t2, _mm_maskload_ps(in[14] + i, mask), 3);
      t3        = _mm512_insertf32x4(t3, _mm_maskload_ps(in[3] + i, mask), 0);
      t3        = _mm512_insertf32x4(t3, _mm_maskload_ps(in[7] + i, mask), 1);
      t3        = _mm512_insertf32x4(t3, _mm_maskload_ps(in[11] + i, mask), 2);
      t3        = _mm512_insertf32x4(t3, _mm_maskload_ps(in[15] + i, mask), 3);
      __m512 v0 = _mm512_shuffle_ps(t0, t1, 0x44);
      __m512 v1 = _mm512_shuffle_ps(t0, t1, 0xee);
      __m512 v2 = _mm512_shuffle_ps(t2, t3, 0x44);
      __m512 v3 = _mm512_shuffle_ps(t2, t3, 0xee);
      out[i + 0].data = _mm512_shuffle_ps(v0, v2, 0x88);
      out[i + 1].data = _mm512_shuffle_ps(v0, v2, 0xdd);
      out[i + 2].data = _mm512_shuffle_ps(v1, v3, 0x88);
    }
  else if (remainder == 2)
    {
      __m512 v0, v1 = {};
 
      // Use 64-bit integer loads
      __m128i t0, t1;
      t0 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[0] + i));
      t1 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[1] + i));
      v0 =
        _mm512_insertf32x4(v1, _mm_castsi128_ps(_mm_unpacklo_epi64(t0, t1)), 0);
      t0 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[4] + i));
      t1 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[5] + i));
      v0 =
        _mm512_insertf32x4(v0, _mm_castsi128_ps(_mm_unpacklo_epi64(t0, t1)), 1);
      t0 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[8] + i));
      t1 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[9] + i));
      v0 =
        _mm512_insertf32x4(v0, _mm_castsi128_ps(_mm_unpacklo_epi64(t0, t1)), 2);
      t0 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[12] + i));
      t1 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[13] + i));
      v0 =
        _mm512_insertf32x4(v0, _mm_castsi128_ps(_mm_unpacklo_epi64(t0, t1)), 3);
      t0 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[2] + i));
      t1 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[3] + i));
      v1 =
        _mm512_insertf32x4(v1, _mm_castsi128_ps(_mm_unpacklo_epi64(t0, t1)), 0);
      t0 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[6] + i));
      t1 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[7] + i));
      v1 =
        _mm512_insertf32x4(v1, _mm_castsi128_ps(_mm_unpacklo_epi64(t0, t1)), 1);
      t0 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[10] + i));
      t1 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[11] + i));
      v1 =
        _mm512_insertf32x4(v1, _mm_castsi128_ps(_mm_unpacklo_epi64(t0, t1)), 2);
      t0 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[14] + i));
      t1 = _mm_loadl_epi64(reinterpret_cast<const __m128i *>(in[15] + i));
      v1 =
        _mm512_insertf32x4(v1, _mm_castsi128_ps(_mm_unpacklo_epi64(t0, t1)), 3);
      out[i + 0].data = _mm512_shuffle_ps(v0, v1, 0x88);
      out[i + 1].data = _mm512_shuffle_ps(v0, v1, 0xdd);
    }
  else if (remainder == 1)
    for (unsigned int v = 0; v < 16; ++v)
      out[i][v] = in[v][i];
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_load_and_transpose(const unsigned int          n_entries,
                              const float                *in,
                              const unsigned int         *offsets,
                              VectorizedArray<float, 16> *out)
{
  vectorized_load_and_transpose(n_entries,
                                {{in + offsets[0],
                                  in + offsets[1],
                                  in + offsets[2],
                                  in + offsets[3],
                                  in + offsets[4],
                                  in + offsets[5],
                                  in + offsets[6],
                                  in + offsets[7],
                                  in + offsets[8],
                                  in + offsets[9],
                                  in + offsets[10],
                                  in + offsets[11],
                                  in + offsets[12],
                                  in + offsets[13],
                                  in + offsets[14],
                                  in + offsets[15]}},
                                out);
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_transpose_and_store(const bool                        add_into,
                               const unsigned int                n_entries,
                               const VectorizedArray<float, 16> *in,
                               const std::array<float *, 16>    &out)
{
  const unsigned int n_chunks = n_entries / 4;
  for (unsigned int i = 0; i < n_chunks; ++i)
    {
      __m512 t0 = _mm512_shuffle_ps(in[4 * i].data, in[1 + 4 * i].data, 0x44);
      __m512 t1 = _mm512_shuffle_ps(in[4 * i].data, in[1 + 4 * i].data, 0xee);
      __m512 t2 =
        _mm512_shuffle_ps(in[2 + 4 * i].data, in[3 + 4 * i].data, 0x44);
      __m512 t3 =
        _mm512_shuffle_ps(in[2 + 4 * i].data, in[3 + 4 * i].data, 0xee);
      __m512 u0 = _mm512_shuffle_ps(t0, t2, 0x88);
      __m512 u1 = _mm512_shuffle_ps(t0, t2, 0xdd);
      __m512 u2 = _mm512_shuffle_ps(t1, t3, 0x88);
      __m512 u3 = _mm512_shuffle_ps(t1, t3, 0xdd);
 
      __m128 res0  = _mm512_extractf32x4_ps(u0, 0);
      __m128 res4  = _mm512_extractf32x4_ps(u0, 1);
      __m128 res8  = _mm512_extractf32x4_ps(u0, 2);
      __m128 res12 = _mm512_extractf32x4_ps(u0, 3);
      __m128 res1  = _mm512_extractf32x4_ps(u1, 0);
      __m128 res5  = _mm512_extractf32x4_ps(u1, 1);
      __m128 res9  = _mm512_extractf32x4_ps(u1, 2);
      __m128 res13 = _mm512_extractf32x4_ps(u1, 3);
      __m128 res2  = _mm512_extractf32x4_ps(u2, 0);
      __m128 res6  = _mm512_extractf32x4_ps(u2, 1);
      __m128 res10 = _mm512_extractf32x4_ps(u2, 2);
      __m128 res14 = _mm512_extractf32x4_ps(u2, 3);
      __m128 res3  = _mm512_extractf32x4_ps(u3, 0);
      __m128 res7  = _mm512_extractf32x4_ps(u3, 1);
      __m128 res11 = _mm512_extractf32x4_ps(u3, 2);
      __m128 res15 = _mm512_extractf32x4_ps(u3, 3);
 
      // Cannot use the same store instructions in both paths of the 'if'
      // because the compiler cannot know that there is no aliasing between
      // pointers
      if (add_into)
        {
          res0 = _mm_add_ps(_mm_loadu_ps(out[0] + 4 * i), res0);
          _mm_storeu_ps(out[0] + 4 * i, res0);
          res1 = _mm_add_ps(_mm_loadu_ps(out[1] + 4 * i), res1);
          _mm_storeu_ps(out[1] + 4 * i, res1);
          res2 = _mm_add_ps(_mm_loadu_ps(out[2] + 4 * i), res2);
          _mm_storeu_ps(out[2] + 4 * i, res2);
          res3 = _mm_add_ps(_mm_loadu_ps(out[3] + 4 * i), res3);
          _mm_storeu_ps(out[3] + 4 * i, res3);
          res4 = _mm_add_ps(_mm_loadu_ps(out[4] + 4 * i), res4);
          _mm_storeu_ps(out[4] + 4 * i, res4);
          res5 = _mm_add_ps(_mm_loadu_ps(out[5] + 4 * i), res5);
          _mm_storeu_ps(out[5] + 4 * i, res5);
          res6 = _mm_add_ps(_mm_loadu_ps(out[6] + 4 * i), res6);
          _mm_storeu_ps(out[6] + 4 * i, res6);
          res7 = _mm_add_ps(_mm_loadu_ps(out[7] + 4 * i), res7);
          _mm_storeu_ps(out[7] + 4 * i, res7);
          res8 = _mm_add_ps(_mm_loadu_ps(out[8] + 4 * i), res8);
          _mm_storeu_ps(out[8] + 4 * i, res8);
          res9 = _mm_add_ps(_mm_loadu_ps(out[9] + 4 * i), res9);
          _mm_storeu_ps(out[9] + 4 * i, res9);
          res10 = _mm_add_ps(_mm_loadu_ps(out[10] + 4 * i), res10);
          _mm_storeu_ps(out[10] + 4 * i, res10);
          res11 = _mm_add_ps(_mm_loadu_ps(out[11] + 4 * i), res11);
          _mm_storeu_ps(out[11] + 4 * i, res11);
          res12 = _mm_add_ps(_mm_loadu_ps(out[12] + 4 * i), res12);
          _mm_storeu_ps(out[12] + 4 * i, res12);
          res13 = _mm_add_ps(_mm_loadu_ps(out[13] + 4 * i), res13);
          _mm_storeu_ps(out[13] + 4 * i, res13);
          res14 = _mm_add_ps(_mm_loadu_ps(out[14] + 4 * i), res14);
          _mm_storeu_ps(out[14] + 4 * i, res14);
          res15 = _mm_add_ps(_mm_loadu_ps(out[15] + 4 * i), res15);
          _mm_storeu_ps(out[15] + 4 * i, res15);
        }
      else
        {
          _mm_storeu_ps(out[0] + 4 * i, res0);
          _mm_storeu_ps(out[1] + 4 * i, res1);
          _mm_storeu_ps(out[2] + 4 * i, res2);
          _mm_storeu_ps(out[3] + 4 * i, res3);
          _mm_storeu_ps(out[4] + 4 * i, res4);
          _mm_storeu_ps(out[5] + 4 * i, res5);
          _mm_storeu_ps(out[6] + 4 * i, res6);
          _mm_storeu_ps(out[7] + 4 * i, res7);
          _mm_storeu_ps(out[8] + 4 * i, res8);
          _mm_storeu_ps(out[9] + 4 * i, res9);
          _mm_storeu_ps(out[10] + 4 * i, res10);
          _mm_storeu_ps(out[11] + 4 * i, res11);
          _mm_storeu_ps(out[12] + 4 * i, res12);
          _mm_storeu_ps(out[13] + 4 * i, res13);
          _mm_storeu_ps(out[14] + 4 * i, res14);
          _mm_storeu_ps(out[15] + 4 * i, res15);
        }
    }
 
  if (add_into)
    for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
      for (unsigned int v = 0; v < 16; ++v)
        out[v][i] += in[i][v];
  else
    for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
      for (unsigned int v = 0; v < 16; ++v)
        out[v][i] = in[i][v];
}
 
 

template <>
inline DEAL_II_ALWAYS_INLINE void
vectorized_transpose_and_store(const bool                        add_into,
                               const unsigned int                n_entries,
                               const VectorizedArray<float, 16> *in,
                               const unsigned int               *offsets,
                               float                            *out)
{
  vectorized_transpose_and_store(add_into,
                                 n_entries,
                                 in,
                                 {{out + offsets[0],
                                   out + offsets[1],
                                   out + offsets[2],
                                   out + offsets[3],
                                   out + offsets[4],
                                   out + offsets[5],
                                   out + offsets[6],
                                   out + offsets[7],
                                   out + offsets[8],
                                   out + offsets[9],
                                   out + offsets[10],
                                   out + offsets[11],
                                   out + offsets[12],
                                   out + offsets[13],
                                   out + offsets[14],
                                   out + offsets[15]}});
}
 
#  endif
 
#  if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ALTIVEC__) && \
    defined(__VSX__)
 
template <>
class VectorizedArray<double, 2>
  : public VectorizedArrayBase<VectorizedArray<double, 2>, 2>
{
public:
  using value_type = double;

  static constexpr bool is_implemented = true;

  VectorizedArray() = default;

  VectorizedArray(const double scalar)
  {
    this->operator=(scalar);
  }

  template <typename U>
  VectorizedArray(const std::initializer_list<U> &list)
    : VectorizedArrayBase<VectorizedArray<double, 2>, 2>(list)
  {}

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator=(const double x) &
  {
    data = vec_splats(x);
 
    // Some compilers believe that vec_splats sets 'x', but that's not true.
    // They then warn about setting a variable and not using it. Suppress the
    // warning by "using" the variable:
    (void)x;
    return *this;
  }

  VectorizedArray &
  operator=(const double scalar) && = delete;

  DEAL_II_ALWAYS_INLINE
  double &
  operator[](const unsigned int comp)
  {
    AssertIndexRange(comp, 2);
    return *(reinterpret_cast<double *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  const double &
  operator[](const unsigned int comp) const
  {
    AssertIndexRange(comp, 2);
    return *(reinterpret_cast<const double *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator+=(const VectorizedArray &vec)
  {
    data = vec_add(data, vec.data);
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator-=(const VectorizedArray &vec)
  {
    data = vec_sub(data, vec.data);
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator*=(const VectorizedArray &vec)
  {
    data = vec_mul(data, vec.data);
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator/=(const VectorizedArray &vec)
  {
    data = vec_div(data, vec.data);
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  void
  load(const double *ptr)
  {
    data = vec_vsx_ld(0, ptr);
  }

  DEAL_II_ALWAYS_INLINE
  void
  store(double *ptr) const
  {
    vec_vsx_st(data, 0, ptr);
  }

  DEAL_II_ALWAYS_INLINE
  void
  streaming_store(double *ptr) const
  {
    store(ptr);
  }

  DEAL_II_ALWAYS_INLINE
  void
  gather(const double *base_ptr, const unsigned int *offsets)
  {
    this->operator=(0.);
    for (unsigned int i = 0; i < 2; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        *(reinterpret_cast<double *>(&data) + i) = base_ptr[offsets[i]];
  }

  DEAL_II_ALWAYS_INLINE
  void
  scatter(const unsigned int *offsets, double *base_ptr) const
  {
    for (unsigned int i = 0; i < 2; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        base_ptr[offsets[i]] = *(reinterpret_cast<const double *>(&data) + i);
  }

  __vector double data;
 
private:
  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_sqrt() const
  {
    VectorizedArray res;
    res.data = vec_sqrt(data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_abs() const
  {
    VectorizedArray res;
    res.data = vec_abs(data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_max(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = vec_max(data, other.data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_min(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = vec_min(data, other.data);
    return res;
  }
 
  // Make a few functions friends.
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::sqrt(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::abs(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::max(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::min(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
};
 
 
 
template <>
class VectorizedArray<float, 4>
  : public VectorizedArrayBase<VectorizedArray<float, 4>, 4>
{
public:
  using value_type = float;

  static constexpr bool is_implemented = true;

  VectorizedArray() = default;

  VectorizedArray(const float scalar)
  {
    this->operator=(scalar);
  }

  template <typename U>
  VectorizedArray(const std::initializer_list<U> &list)
    : VectorizedArrayBase<VectorizedArray<float, 4>, 4>(list)
  {}

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator=(const float x) &
  {
    data = vec_splats(x);
 
    // Some compilers believe that vec_splats sets 'x', but that's not true.
    // They then warn about setting a variable and not using it. Suppress the
    // warning by "using" the variable:
    (void)x;
    return *this;
  }

  VectorizedArray &
  operator=(const float scalar) && = delete;

  DEAL_II_ALWAYS_INLINE
  float &
  operator[](const unsigned int comp)
  {
    AssertIndexRange(comp, 4);
    return *(reinterpret_cast<float *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  const float &
  operator[](const unsigned int comp) const
  {
    AssertIndexRange(comp, 4);
    return *(reinterpret_cast<const float *>(&data) + comp);
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator+=(const VectorizedArray &vec)
  {
    data = vec_add(data, vec.data);
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator-=(const VectorizedArray &vec)
  {
    data = vec_sub(data, vec.data);
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator*=(const VectorizedArray &vec)
  {
    data = vec_mul(data, vec.data);
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray &
  operator/=(const VectorizedArray &vec)
  {
    data = vec_div(data, vec.data);
    return *this;
  }

  DEAL_II_ALWAYS_INLINE
  void
  load(const float *ptr)
  {
    data = vec_vsx_ld(0, ptr);
  }

  DEAL_II_ALWAYS_INLINE
  void
  store(float *ptr) const
  {
    vec_vsx_st(data, 0, ptr);
  }

  DEAL_II_ALWAYS_INLINE
  void
  streaming_store(float *ptr) const
  {
    store(ptr);
  }

  DEAL_II_ALWAYS_INLINE
  void
  gather(const float *base_ptr, const unsigned int *offsets)
  {
    this->operator=(0.f);
    for (unsigned int i = 0; i < 4; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        *(reinterpret_cast<float *>(&data) + i) = base_ptr[offsets[i]];
  }

  DEAL_II_ALWAYS_INLINE
  void
  scatter(const unsigned int *offsets, float *base_ptr) const
  {
    for (unsigned int i = 0; i < 4; ++i)
      if (offsets[i] != numbers::invalid_unsigned_int)
        base_ptr[offsets[i]] = *(reinterpret_cast<const float *>(&data) + i);
  }

  __vector float data;
 
private:
  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_sqrt() const
  {
    VectorizedArray res;
    res.data = vec_sqrt(data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_abs() const
  {
    VectorizedArray res;
    res.data = vec_abs(data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_max(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = vec_max(data, other.data);
    return res;
  }

  DEAL_II_ALWAYS_INLINE
  VectorizedArray
  get_min(const VectorizedArray &other) const
  {
    VectorizedArray res;
    res.data = vec_min(data, other.data);
    return res;
  }
 
  // Make a few functions friends.
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::sqrt(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::abs(const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::max(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
  template <typename Number2, std::size_t width2>
  friend VectorizedArray<Number2, width2>
  std::min(const VectorizedArray<Number2, width2> &,
           const VectorizedArray<Number2, width2> &);
};
 
#  endif // if DEAL_II_VECTORIZATION_LEVEL >=1 && defined(__ALTIVEC__) &&
         // defined(__VSX__)
 
 
#endif // DOXYGEN
 
 


template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE bool
operator==(const VectorizedArray<Number, width> &lhs,
           const VectorizedArray<Number, width> &rhs)
{
  for (unsigned int i = 0; i < VectorizedArray<Number, width>::size(); ++i)
    if (lhs[i] != rhs[i])
      return false;
 
  return true;
}
 

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number, width>
                             operator+(const VectorizedArray<Number, width> &u,
          const VectorizedArray<Number, width> &v)
{
  VectorizedArray<Number, width> tmp = u;
  return tmp += v;
}

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number, width>
                             operator-(const VectorizedArray<Number, width> &u,
          const VectorizedArray<Number, width> &v)
{
  VectorizedArray<Number, width> tmp = u;
  return tmp -= v;
}

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number, width>
                             operator*(const VectorizedArray<Number, width> &u,
          const VectorizedArray<Number, width> &v)
{
  VectorizedArray<Number, width> tmp = u;
  return tmp *= v;
}

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number, width>
                             operator/(const VectorizedArray<Number, width> &u,
          const VectorizedArray<Number, width> &v)
{
  VectorizedArray<Number, width> tmp = u;
  return tmp /= v;
}

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number, width>
operator+(const Number &u, const VectorizedArray<Number, width> &v)
{
  VectorizedArray<Number, width> tmp = u;
  return tmp += v;
}

template <std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<float, width>
operator+(const double u, const VectorizedArray<float, width> &v)
{
  VectorizedArray<float, width> tmp = u;
  return tmp += v;
}

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number, width>
operator+(const VectorizedArray<Number, width> &v, const Number &u)
{
  return u + v;
}

template <std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<float, width>
operator+(const VectorizedArray<float, width> &v, const double u)
{
  return u + v;
}

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number, width>
operator-(const Number &u, const VectorizedArray<Number, width> &v)
{
  VectorizedArray<Number, width> tmp = u;
  return tmp -= v;
}

template <std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<float, width>
operator-(const double u, const VectorizedArray<float, width> &v)
{
  VectorizedArray<float, width> tmp = static_cast<float>(u);
  return tmp -= v;
}

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number, width>
operator-(const VectorizedArray<Number, width> &v, const Number &u)
{
  VectorizedArray<Number, width> tmp = u;
  return v - tmp;
}

template <std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<float, width>
operator-(const VectorizedArray<float, width> &v, const double u)
{
  VectorizedArray<float, width> tmp = static_cast<float>(u);
  return v - tmp;
}

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number, width>
operator*(const Number &u, const VectorizedArray<Number, width> &v)
{
  VectorizedArray<Number, width> tmp = u;
  return tmp *= v;
}

template <std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<float, width>
operator*(const double u, const VectorizedArray<float, width> &v)
{
  VectorizedArray<float, width> tmp = static_cast<float>(u);
  return tmp *= v;
}

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number, width>
operator*(const VectorizedArray<Number, width> &v, const Number &u)
{
  return u * v;
}

template <std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<float, width>
operator*(const VectorizedArray<float, width> &v, const double u)
{
  return u * v;
}

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number, width>
operator/(const Number &u, const VectorizedArray<Number, width> &v)
{
  VectorizedArray<Number, width> tmp = u;
  return tmp /= v;
}

template <std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<float, width>
operator/(const double u, const VectorizedArray<float, width> &v)
{
  VectorizedArray<float, width> tmp = static_cast<float>(u);
  return tmp /= v;
}

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number, width>
operator/(const VectorizedArray<Number, width> &v, const Number &u)
{
  VectorizedArray<Number, width> tmp = u;
  return v / tmp;
}

template <std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<float, width>
operator/(const VectorizedArray<float, width> &v, const double u)
{
  VectorizedArray<float, width> tmp = static_cast<float>(u);
  return v / tmp;
}

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number, width>
                             operator+(const VectorizedArray<Number, width> &u)
{
  return u;
}

template <typename Number, std::size_t width>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number, width>
                             operator-(const VectorizedArray<Number, width> &u)
{
  // to get a negative sign, subtract the input from zero (could also
  // multiply by -1, but this one is slightly simpler)
  return VectorizedArray<Number, width>() - u;
}

template <typename Number, std::size_t width>
inline std::ostream &
operator<<(std::ostream &out, const VectorizedArray<Number, width> &p)
{
  constexpr unsigned int n = VectorizedArray<Number, width>::size();
  for (unsigned int i = 0; i < n - 1; ++i)
    out << p[i] << ' ';
  out << p[n - 1];
 
  return out;
}



enum class SIMDComparison : int
{
#if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)
  equal                 = _CMP_EQ_OQ,
  not_equal             = _CMP_NEQ_OQ,
  less_than             = _CMP_LT_OQ,
  less_than_or_equal    = _CMP_LE_OQ,
  greater_than          = _CMP_GT_OQ,
  greater_than_or_equal = _CMP_GE_OQ
#else
  equal,
  not_equal,
  less_than,
  less_than_or_equal,
  greater_than,
  greater_than_or_equal
#endif
};
 

template <SIMDComparison predicate, typename Number>
DEAL_II_ALWAYS_INLINE inline Number
compare_and_apply_mask(const Number &left,
                       const Number &right,
                       const Number &true_value,
                       const Number &false_value)
{
  bool mask;
  switch (predicate)
    {
      case SIMDComparison::equal:
        mask = (left == right);
        break;
      case SIMDComparison::not_equal:
        mask = (left != right);
        break;
      case SIMDComparison::less_than:
        mask = (left < right);
        break;
      case SIMDComparison::less_than_or_equal:
        mask = (left <= right);
        break;
      case SIMDComparison::greater_than:
        mask = (left > right);
        break;
      case SIMDComparison::greater_than_or_equal:
        mask = (left >= right);
        break;
    }
 
  return mask ? true_value : false_value;
}
 

template <SIMDComparison predicate, typename Number>
DEAL_II_ALWAYS_INLINE inline VectorizedArray<Number, 1>
compare_and_apply_mask(const VectorizedArray<Number, 1> &left,
                       const VectorizedArray<Number, 1> &right,
                       const VectorizedArray<Number, 1> &true_value,
                       const VectorizedArray<Number, 1> &false_value)
{
  VectorizedArray<Number, 1> result;
  result.data = compare_and_apply_mask<predicate, Number>(left.data,
                                                          right.data,
                                                          true_value.data,
                                                          false_value.data);
  return result;
}

 
#ifndef DOXYGEN
#  if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 512 && defined(__AVX512F__)
 
template <SIMDComparison predicate>
DEAL_II_ALWAYS_INLINE inline VectorizedArray<float, 16>
compare_and_apply_mask(const VectorizedArray<float, 16> &left,
                       const VectorizedArray<float, 16> &right,
                       const VectorizedArray<float, 16> &true_values,
                       const VectorizedArray<float, 16> &false_values)
{
  const __mmask16 mask =
    _mm512_cmp_ps_mask(left.data, right.data, static_cast<int>(predicate));
  VectorizedArray<float, 16> result;
  result.data = _mm512_mask_mov_ps(false_values.data, mask, true_values.data);
  return result;
}
 
 
 
template <SIMDComparison predicate>
DEAL_II_ALWAYS_INLINE inline VectorizedArray<double, 8>
compare_and_apply_mask(const VectorizedArray<double, 8> &left,
                       const VectorizedArray<double, 8> &right,
                       const VectorizedArray<double, 8> &true_values,
                       const VectorizedArray<double, 8> &false_values)
{
  const __mmask16 mask =
    _mm512_cmp_pd_mask(left.data, right.data, static_cast<int>(predicate));
  VectorizedArray<double, 8> result;
  result.data = _mm512_mask_mov_pd(false_values.data, mask, true_values.data);
  return result;
}
 
#  endif
 
#  if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 256 && defined(__AVX__)
 
template <SIMDComparison predicate>
DEAL_II_ALWAYS_INLINE inline VectorizedArray<float, 8>
compare_and_apply_mask(const VectorizedArray<float, 8> &left,
                       const VectorizedArray<float, 8> &right,
                       const VectorizedArray<float, 8> &true_values,
                       const VectorizedArray<float, 8> &false_values)
{
  const auto mask =
    _mm256_cmp_ps(left.data, right.data, static_cast<int>(predicate));
 
  VectorizedArray<float, 8> result;
  result.data = _mm256_blendv_ps(false_values.data, true_values.data, mask);
  return result;
}
 
 
template <SIMDComparison predicate>
DEAL_II_ALWAYS_INLINE inline VectorizedArray<double, 4>
compare_and_apply_mask(const VectorizedArray<double, 4> &left,
                       const VectorizedArray<double, 4> &right,
                       const VectorizedArray<double, 4> &true_values,
                       const VectorizedArray<double, 4> &false_values)
{
  const auto mask =
    _mm256_cmp_pd(left.data, right.data, static_cast<int>(predicate));
 
  VectorizedArray<double, 4> result;
  result.data = _mm256_blendv_pd(false_values.data, true_values.data, mask);
  return result;
}
 
#  endif
 
#  if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__SSE2__)
 
template <SIMDComparison predicate>
DEAL_II_ALWAYS_INLINE inline VectorizedArray<float, 4>
compare_and_apply_mask(const VectorizedArray<float, 4> &left,
                       const VectorizedArray<float, 4> &right,
                       const VectorizedArray<float, 4> &true_values,
                       const VectorizedArray<float, 4> &false_values)
{
  __m128 mask;
  switch (predicate)
    {
      case SIMDComparison::equal:
        mask = _mm_cmpeq_ps(left.data, right.data);
        break;
      case SIMDComparison::not_equal:
        mask = _mm_cmpneq_ps(left.data, right.data);
        break;
      case SIMDComparison::less_than:
        mask = _mm_cmplt_ps(left.data, right.data);
        break;
      case SIMDComparison::less_than_or_equal:
        mask = _mm_cmple_ps(left.data, right.data);
        break;
      case SIMDComparison::greater_than:
        mask = _mm_cmpgt_ps(left.data, right.data);
        break;
      case SIMDComparison::greater_than_or_equal:
        mask = _mm_cmpge_ps(left.data, right.data);
        break;
    }
 
  VectorizedArray<float, 4> result;
  result.data = _mm_or_ps(_mm_and_ps(mask, true_values.data),
                          _mm_andnot_ps(mask, false_values.data));
 
  return result;
}
 
 
template <SIMDComparison predicate>
DEAL_II_ALWAYS_INLINE inline VectorizedArray<double, 2>
compare_and_apply_mask(const VectorizedArray<double, 2> &left,
                       const VectorizedArray<double, 2> &right,
                       const VectorizedArray<double, 2> &true_values,
                       const VectorizedArray<double, 2> &false_values)
{
  __m128d mask;
  switch (predicate)
    {
      case SIMDComparison::equal:
        mask = _mm_cmpeq_pd(left.data, right.data);
        break;
      case SIMDComparison::not_equal:
        mask = _mm_cmpneq_pd(left.data, right.data);
        break;
      case SIMDComparison::less_than:
        mask = _mm_cmplt_pd(left.data, right.data);
        break;
      case SIMDComparison::less_than_or_equal:
        mask = _mm_cmple_pd(left.data, right.data);
        break;
      case SIMDComparison::greater_than:
        mask = _mm_cmpgt_pd(left.data, right.data);
        break;
      case SIMDComparison::greater_than_or_equal:
        mask = _mm_cmpge_pd(left.data, right.data);
        break;
    }
 
  VectorizedArray<double, 2> result;
  result.data = _mm_or_pd(_mm_and_pd(mask, true_values.data),
                          _mm_andnot_pd(mask, false_values.data));
 
  return result;
}
 
#  endif
 
#  if DEAL_II_VECTORIZATION_WIDTH_IN_BITS >= 128 && defined(__ARM_NEON)
 
template <SIMDComparison predicate>
DEAL_II_ALWAYS_INLINE inline VectorizedArray<float, 4>
compare_and_apply_mask(const VectorizedArray<float, 4> &left,
                       const VectorizedArray<float, 4> &right,
                       const VectorizedArray<float, 4> &true_values,
                       const VectorizedArray<float, 4> &false_values)
{
  uint32x4_t mask;
  switch (predicate)
    {
      case SIMDComparison::equal:
        mask = vceqq_f32(left.data, right.data);
        break;
      case SIMDComparison::not_equal:
        mask = vmvnq_u32(vceqq_f32(left.data, right.data));
        break;
      case SIMDComparison::less_than:
        mask = vcltq_f32(left.data, right.data);
        break;
      case SIMDComparison::less_than_or_equal:
        mask = vcleq_f32(left.data, right.data);
        break;
      case SIMDComparison::greater_than:
        mask = vcgtq_f32(left.data, right.data);
        break;
      case SIMDComparison::greater_than_or_equal:
        mask = vcgeq_f32(left.data, right.data);
        break;
    }
 
  VectorizedArray<float, 4> result;
  result.data = vreinterpretq_f32_u32(vorrq_u32(
    vandq_u32(mask, vreinterpretq_u32_f32(true_values.data)),
    vandq_u32(vmvnq_u32(mask), vreinterpretq_u32_f32(false_values.data))));
 
  return result;
}
 
 
template <SIMDComparison predicate>
DEAL_II_ALWAYS_INLINE inline VectorizedArray<double, 2>
compare_and_apply_mask(const VectorizedArray<double, 2> &left,
                       const VectorizedArray<double, 2> &right,
                       const VectorizedArray<double, 2> &true_values,
                       const VectorizedArray<double, 2> &false_values)
{
  uint64x2_t mask;
  switch (predicate)
    {
      case SIMDComparison::equal:
        mask = vceqq_f64(left.data, right.data);
        break;
      case SIMDComparison::not_equal:
        mask = vreinterpretq_u64_u32(
          vmvnq_u32(vreinterpretq_u32_u64(vceqq_f64(left.data, right.data))));
        break;
      case SIMDComparison::less_than:
        mask = vcltq_f64(left.data, right.data);
        break;
      case SIMDComparison::less_than_or_equal:
        mask = vcleq_f64(left.data, right.data);
        break;
      case SIMDComparison::greater_than:
        mask = vcgtq_f64(left.data, right.data);
        break;
      case SIMDComparison::greater_than_or_equal:
        mask = vcgeq_f64(left.data, right.data);
        break;
    }
 
  VectorizedArray<double, 2> result;
  result.data = vreinterpretq_f64_u64(vorrq_u64(
    vandq_u64(mask, vreinterpretq_u64_f64(true_values.data)),
    vandq_u64(vreinterpretq_u64_u32(vmvnq_u32(vreinterpretq_u32_u64(mask))),
              vreinterpretq_u64_f64(false_values.data))));
 
  return result;
}
 
#  endif
#endif // DOXYGEN
 
 
namespace internal
{
  template <typename T>
  struct VectorizedArrayTrait
  {
    using value_type = T;

    static constexpr std::size_t
    width()
    {
      return 1;
    }

    using vectorized_value_type = VectorizedArray<T>;

    static constexpr std::size_t
    stride()
    {
      return vectorized_value_type::size();
    }

    static value_type &
    get(value_type &value, unsigned int c)
    {
      AssertIndexRange(c, 1);
      (void)c;
 
      return value;
    }

    static const value_type &
    get(const value_type &value, unsigned int c)
    {
      AssertIndexRange(c, 1);
      (void)c;
 
      return value;
    }

    static value_type &
    get_from_vectorized(vectorized_value_type &values, unsigned int c)
    {
      AssertIndexRange(c, stride());
 
      return values[c];
    }

    static const value_type &
    get_from_vectorized(const vectorized_value_type &values, unsigned int c)
    {
      AssertIndexRange(c, stride());
 
      return values[c];
    }
  };
 
  template <typename T, std::size_t width_>
  struct VectorizedArrayTrait<VectorizedArray<T, width_>>
  {
    using value_type = T;

    static constexpr std::size_t
    width()
    {
      return width_;
    }

    using vectorized_value_type = VectorizedArray<T, width_>;

    static constexpr std::size_t
    stride()
    {
      return 1;
    }

    static value_type &
    get(vectorized_value_type &values, unsigned int c)
    {
      AssertIndexRange(c, width_);
 
      return values[c];
    }

    static const value_type &
    get(const vectorized_value_type &values, unsigned int c)
    {
      AssertIndexRange(c, width_);
 
      return values[c];
    }

    static vectorized_value_type &
    get_from_vectorized(vectorized_value_type &values, unsigned int c)
    {
      (void)c;
      AssertIndexRange(c, stride());
 
      return values;
    }

    static const vectorized_value_type &
    get_from_vectorized(const vectorized_value_type &values, unsigned int c)
    {
      (void)c;
      AssertIndexRange(c, stride());
 
      return values;
    }
  };
} // namespace internal
 
 
DEAL_II_NAMESPACE_CLOSE

namespace std
{
  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  sin(const ::VectorizedArray<Number, width> &x)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::sin(x[i]);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  cos(const ::VectorizedArray<Number, width> &x)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::cos(x[i]);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  tan(const ::VectorizedArray<Number, width> &x)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::tan(x[i]);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  acos(const ::VectorizedArray<Number, width> &x)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::acos(x[i]);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  asin(const ::VectorizedArray<Number, width> &x)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::asin(x[i]);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  atan(const ::VectorizedArray<Number, width> &x)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::atan(x[i]);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  cosh(const ::VectorizedArray<Number, width> &x)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::cosh(x[i]);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  sinh(const ::VectorizedArray<Number, width> &x)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::sinh(x[i]);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  tanh(const ::VectorizedArray<Number, width> &x)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::tanh(x[i]);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  acosh(const ::VectorizedArray<Number, width> &x)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::acosh(x[i]);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  asinh(const ::VectorizedArray<Number, width> &x)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::asinh(x[i]);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  atanh(const ::VectorizedArray<Number, width> &x)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::atanh(x[i]);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  exp(const ::VectorizedArray<Number, width> &x)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::exp(x[i]);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  log(const ::VectorizedArray<Number, width> &x)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::log(x[i]);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  sqrt(const ::VectorizedArray<Number, width> &x)
  {
    return x.get_sqrt();
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  pow(const ::VectorizedArray<Number, width> &x, const Number p)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::pow(x[i], p);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  pow(const ::VectorizedArray<Number, width> &x,
      const ::VectorizedArray<Number, width> &p)
  {
    ::VectorizedArray<Number, width> out;
    for (unsigned int i = 0; i < ::VectorizedArray<Number, width>::size();
         ++i)
      out[i] = std::pow(x[i], p[i]);
    return out;
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  abs(const ::VectorizedArray<Number, width> &x)
  {
    return x.get_abs();
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  max(const ::VectorizedArray<Number, width> &x,
      const ::VectorizedArray<Number, width> &y)
  {
    return x.get_max(y);
  }
 
 

  template <typename Number, std::size_t width>
  inline ::VectorizedArray<Number, width>
  min(const ::VectorizedArray<Number, width> &x,
      const ::VectorizedArray<Number, width> &y)
  {
    return x.get_min(y);
  }
 
 

  template <class T>
  struct iterator_traits<::VectorizedArrayIterator<T>>
  {
#ifdef DEAL_II_HAVE_CXX20
    using iterator_category = contiguous_iterator_tag;
#else
    using iterator_category = random_access_iterator_tag;
#endif
    using value_type      = typename T::value_type;
    using difference_type = std::ptrdiff_t;
  };
 
} // namespace std
 
#endif