vectorization.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------


#ifndef dealii_vectorization_h
#define dealii_vectorization_h

#include <deal.II/base/config.h>

#include <deal.II/base/exceptions.h>
#include <deal.II/base/template_constraints.h>

#include <cmath>

// Note:
// The flag DEAL_II_COMPILER_VECTORIZATION_LEVEL is essentially constructed
// according to the following scheme
// #ifdef __AVX512F__
// #define DEAL_II_COMPILER_VECTORIZATION_LEVEL 3
// #elif defined (__AVX__)
// #define DEAL_II_COMPILER_VECTORIZATION_LEVEL 2
// #elif defined (__SSE2__)
// #define DEAL_II_COMPILER_VECTORIZATION_LEVEL 1
// #else
// #define DEAL_II_COMPILER_VECTORIZATION_LEVEL 0
// #endif
// In addition to checking the flags __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_COMPILER_VECTORIZATION_LEVEL >= 2 && !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_COMPILER_VECTORIZATION_LEVEL >= 3 && !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

#if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 2 // AVX, AVX-512
#  include <immintrin.h>
#elif DEAL_II_COMPILER_VECTORIZATION_LEVEL == 1 // SSE2
#  include <emmintrin.h>
#endif


DEAL_II_NAMESPACE_OPEN


namespace internal
{
  template <typename T>
  struct NumberType<VectorizedArray<T>>
  {
    static const VectorizedArray<T> &
    value(const VectorizedArray<T> &t)
    {
      return t;
    }

    static VectorizedArray<T>
    value(const T &t)
    {
      VectorizedArray<T> tmp;
      tmp = t;
      return tmp;
    }
  };
} // namespace internal


// 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>
struct EnableIfScalar<VectorizedArray<Number>>
{
  using type = VectorizedArray<typename EnableIfScalar<Number>::type>;
};



template <typename Number>
class VectorizedArray
{
public:
  static const unsigned int n_array_elements = 1;

  // POD means that there should be no user-defined constructors, destructors
  // and copy functions (the standard is somewhat relaxed in C++2011, though).

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

  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<Number> &vec)
  {
    data += vec.data;
    return *this;
  }

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

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

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

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

  DEAL_II_ALWAYS_INLINE
  void
  store(Number *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)
  {
    data = base_ptr[offsets[0]];
  }

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

  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;
  }

  template <typename Number2>
  friend VectorizedArray<Number2>
  std::sqrt(const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::abs(const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::max(const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::min(const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
};



template <typename Number>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
                             make_vectorized_array(const Number &u)
{
  VectorizedArray<Number> result;
  result = u;
  return result;
}



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



template <typename Number>
inline void
vectorized_transpose_and_store(const bool                     add_into,
                               const unsigned int             n_entries,
                               const VectorizedArray<Number> *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>::n_array_elements;
           ++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>::n_array_elements;
           ++v)
        out[offsets[v] + i] = in[i][v];
}



// for safety, also check that __AVX512F__ is defined in case the user manually
// set some conflicting compile flags which prevent compilation

#if DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 3 && defined(__AVX512F__)

template <>
class VectorizedArray<double>
{
public:
  static const unsigned int n_array_elements = 8;

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

  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 arithmetics, 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
  store(double *ptr) const
  {
    _mm512_storeu_pd(ptr, 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)
  {
    // 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((const float *)offsets);
    const __m256i index     = *((__m256i *)(&index_val));
    data                    = _mm512_i32gather_pd(index, base_ptr, 8);
  }

  DEAL_II_ALWAYS_INLINE
  void
  scatter(const unsigned int *offsets, double *base_ptr) const
  {
    for (unsigned int i = 0; i < 8; ++i)
      for (unsigned int j = i + 1; j < 8; ++j)
        Assert(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((const float *)offsets);
    const __m256i index     = *((__m256i *)(&index_val));
    _mm512_i32scatter_pd(base_ptr, index, data, 8);
  }

  __m512d data;

private:
  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 = (__m512d)_mm512_andnot_epi64((__m512i)mask, (__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;
  }

  template <typename Number2>
  friend VectorizedArray<Number2>
  std::sqrt(const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::abs(const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::max(const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::min(const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
};



template <>
inline void
vectorized_load_and_transpose(const unsigned int       n_entries,
                              const double *           in,
                              const unsigned int *     offsets,
                              VectorizedArray<double> *out)
{
  const unsigned int n_chunks = n_entries / 4;
  for (unsigned int outer = 0; outer < 8; outer += 4)
    {
      const double *in0 = in + offsets[0 + outer];
      const double *in1 = in + offsets[1 + outer];
      const double *in2 = in + offsets[2 + outer];
      const double *in3 = in + offsets[3 + outer];

      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);
          *(__m256d *)((double *)(&out[4 * i + 0].data) + outer) =
            _mm256_unpacklo_pd(t0, t1);
          *(__m256d *)((double *)(&out[4 * i + 1].data) + outer) =
            _mm256_unpackhi_pd(t0, t1);
          *(__m256d *)((double *)(&out[4 * i + 2].data) + outer) =
            _mm256_unpacklo_pd(t2, t3);
          *(__m256d *)((double *)(&out[4 * i + 3].data) + outer) =
            _mm256_unpackhi_pd(t2, t3);
        }
      for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
        for (unsigned int v = 0; v < 4; ++v)
          out[i][outer + v] = in[offsets[v + outer] + i];
    }
}



template <>
inline void
vectorized_transpose_and_store(const bool                     add_into,
                               const unsigned int             n_entries,
                               const VectorizedArray<double> *in,
                               const unsigned int *           offsets,
                               double *                       out)
{
  const unsigned int n_chunks = n_entries / 4;
  // do not do full transpose because the code is too long and will most
  // likely not pay off. rather do the transposition on the vectorized array
  // on size smaller, mm256d
  for (unsigned int outer = 0; outer < 8; outer += 4)
    {
      double *out0 = out + offsets[0 + outer];
      double *out1 = out + offsets[1 + outer];
      double *out2 = out + offsets[2 + outer];
      double *out3 = out + offsets[3 + outer];
      for (unsigned int i = 0; i < n_chunks; ++i)
        {
          __m256d u0 =
            *(const __m256d *)((const double *)(&in[4 * i + 0].data) + outer);
          __m256d u1 =
            *(const __m256d *)((const double *)(&in[4 * i + 1].data) + outer);
          __m256d u2 =
            *(const __m256d *)((const double *)(&in[4 * i + 2].data) + outer);
          __m256d u3 =
            *(const __m256d *)((const double *)(&in[4 * i + 3].data) + outer);
          __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);
            }
        }
      if (add_into)
        for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
          for (unsigned int v = 0; v < 4; ++v)
            out[offsets[v + outer] + i] += in[i][v + outer];
      else
        for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
          for (unsigned int v = 0; v < 4; ++v)
            out[offsets[v + outer] + i] = in[i][v + outer];
    }
}



template <>
class VectorizedArray<float>
{
public:
  static const unsigned int n_array_elements = 16;

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

  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 arithmetics, 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)
  {
    // 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((const float *)offsets);
    const __m512i index     = *((__m512i *)(&index_val));
    data                    = _mm512_i32gather_ps(index, base_ptr, 4);
  }

  DEAL_II_ALWAYS_INLINE
  void
  scatter(const unsigned int *offsets, float *base_ptr) const
  {
    for (unsigned int i = 0; i < 16; ++i)
      for (unsigned int j = i + 1; j < 16; ++j)
        Assert(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((const float *)offsets);
    const __m512i index     = *((__m512i *)(&index_val));
    _mm512_i32scatter_ps(base_ptr, index, data, 4);
  }

  __m512 data;

private:
  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 = (__m512)_mm512_andnot_epi32((__m512i)mask, (__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;
  }

  template <typename Number2>
  friend VectorizedArray<Number2>
  std::sqrt(const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::abs(const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::max(const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::min(const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
};



template <>
inline void
vectorized_load_and_transpose(const unsigned int      n_entries,
                              const float *           in,
                              const unsigned int *    offsets,
                              VectorizedArray<float> *out)
{
  const unsigned int n_chunks = n_entries / 4;
  for (unsigned int outer = 0; outer < 16; outer += 8)
    {
      for (unsigned int i = 0; i < n_chunks; ++i)
        {
          __m128 u0 = _mm_loadu_ps(in + 4 * i + offsets[0 + outer]);
          __m128 u1 = _mm_loadu_ps(in + 4 * i + offsets[1 + outer]);
          __m128 u2 = _mm_loadu_ps(in + 4 * i + offsets[2 + outer]);
          __m128 u3 = _mm_loadu_ps(in + 4 * i + offsets[3 + outer]);
          __m128 u4 = _mm_loadu_ps(in + 4 * i + offsets[4 + outer]);
          __m128 u5 = _mm_loadu_ps(in + 4 * i + offsets[5 + outer]);
          __m128 u6 = _mm_loadu_ps(in + 4 * i + offsets[6 + outer]);
          __m128 u7 = _mm_loadu_ps(in + 4 * i + offsets[7 + outer]);
          // To avoid warnings about uninitialized variables, need to initialize
          // one variable with zero before using it.
          __m256 t0, t1, t2, t3 = _mm256_set1_ps(0.F);
          t0        = _mm256_insertf128_ps(t3, u0, 0);
          t0        = _mm256_insertf128_ps(t0, u4, 1);
          t1        = _mm256_insertf128_ps(t3, u1, 0);
          t1        = _mm256_insertf128_ps(t1, u5, 1);
          t2        = _mm256_insertf128_ps(t3, u2, 0);
          t2        = _mm256_insertf128_ps(t2, u6, 1);
          t3        = _mm256_insertf128_ps(t3, u3, 0);
          t3        = _mm256_insertf128_ps(t3, u7, 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);
          *(__m256 *)((float *)(&out[4 * i + 0].data) + outer) =
            _mm256_shuffle_ps(v0, v2, 0x88);
          *(__m256 *)((float *)(&out[4 * i + 1].data) + outer) =
            _mm256_shuffle_ps(v0, v2, 0xdd);
          *(__m256 *)((float *)(&out[4 * i + 2].data) + outer) =
            _mm256_shuffle_ps(v1, v3, 0x88);
          *(__m256 *)((float *)(&out[4 * i + 3].data) + outer) =
            _mm256_shuffle_ps(v1, v3, 0xdd);
        }
      for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
        for (unsigned int v = 0; v < 8; ++v)
          out[i][v + outer] = in[offsets[v + outer] + i];
    }
}



template <>
inline void
vectorized_transpose_and_store(const bool                    add_into,
                               const unsigned int            n_entries,
                               const VectorizedArray<float> *in,
                               const unsigned int *          offsets,
                               float *                       out)
{
  const unsigned int n_chunks = n_entries / 4;
  for (unsigned int outer = 0; outer < 16; outer += 8)
    {
      for (unsigned int i = 0; i < n_chunks; ++i)
        {
          __m256 u0 =
            *(const __m256 *)((const float *)(&in[4 * i + 0].data) + outer);
          __m256 u1 =
            *(const __m256 *)((const float *)(&in[4 * i + 1].data) + outer);
          __m256 u2 =
            *(const __m256 *)((const float *)(&in[4 * i + 2].data) + outer);
          __m256 u3 =
            *(const __m256 *)((const float *)(&in[4 * i + 3].data) + outer);
          __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 + 4 * i + offsets[0 + outer]),
                                res0);
              _mm_storeu_ps(out + 4 * i + offsets[0 + outer], res0);
              res1 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[1 + outer]),
                                res1);
              _mm_storeu_ps(out + 4 * i + offsets[1 + outer], res1);
              res2 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[2 + outer]),
                                res2);
              _mm_storeu_ps(out + 4 * i + offsets[2 + outer], res2);
              res3 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[3 + outer]),
                                res3);
              _mm_storeu_ps(out + 4 * i + offsets[3 + outer], res3);
              res4 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[4 + outer]),
                                res4);
              _mm_storeu_ps(out + 4 * i + offsets[4 + outer], res4);
              res5 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[5 + outer]),
                                res5);
              _mm_storeu_ps(out + 4 * i + offsets[5 + outer], res5);
              res6 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[6 + outer]),
                                res6);
              _mm_storeu_ps(out + 4 * i + offsets[6 + outer], res6);
              res7 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[7 + outer]),
                                res7);
              _mm_storeu_ps(out + 4 * i + offsets[7 + outer], res7);
            }
          else
            {
              _mm_storeu_ps(out + 4 * i + offsets[0 + outer], res0);
              _mm_storeu_ps(out + 4 * i + offsets[1 + outer], res1);
              _mm_storeu_ps(out + 4 * i + offsets[2 + outer], res2);
              _mm_storeu_ps(out + 4 * i + offsets[3 + outer], res3);
              _mm_storeu_ps(out + 4 * i + offsets[4 + outer], res4);
              _mm_storeu_ps(out + 4 * i + offsets[5 + outer], res5);
              _mm_storeu_ps(out + 4 * i + offsets[6 + outer], res6);
              _mm_storeu_ps(out + 4 * i + offsets[7 + outer], res7);
            }
        }
      if (add_into)
        for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
          for (unsigned int v = 0; v < 8; ++v)
            out[offsets[v + outer] + i] += in[i][v + outer];
      else
        for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
          for (unsigned int v = 0; v < 8; ++v)
            out[offsets[v + outer] + i] = in[i][v + outer];
    }
}



#elif DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 2 && defined(__AVX__)

template <>
class VectorizedArray<double>
{
public:
  static const unsigned int n_array_elements = 4;

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

  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 arithmetics, 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
  store(double *ptr) const
  {
    _mm256_storeu_pd(ptr, 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)
  {
#  ifdef __AVX2__
    // 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((const float *)offsets);
    const __m128i index     = *((__m128i *)(&index_val));
    data                    = _mm256_i32gather_pd(base_ptr, index, 8);
#  else
    for (unsigned int i = 0; i < 4; ++i)
      *(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)
      base_ptr[offsets[i]] = *(reinterpret_cast<const double *>(&data) + i);
  }

  __m256d data;

private:
  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;
  }

  template <typename Number2>
  friend VectorizedArray<Number2>
  std::sqrt(const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::abs(const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::max(const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::min(const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
};



template <>
inline void
vectorized_load_and_transpose(const unsigned int       n_entries,
                              const double *           in,
                              const unsigned int *     offsets,
                              VectorizedArray<double> *out)
{
  const unsigned int n_chunks = n_entries / 4;
  const double *     in0      = in + offsets[0];
  const double *     in1      = in + offsets[1];
  const double *     in2      = in + offsets[2];
  const double *     in3      = in + offsets[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);
    }
  for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
    for (unsigned int v = 0; v < 4; ++v)
      out[i][v] = in[offsets[v] + i];
}



template <>
inline void
vectorized_transpose_and_store(const bool                     add_into,
                               const unsigned int             n_entries,
                               const VectorizedArray<double> *in,
                               const unsigned int *           offsets,
                               double *                       out)
{
  const unsigned int n_chunks = n_entries / 4;
  double *           out0     = out + offsets[0];
  double *           out1     = out + offsets[1];
  double *           out2     = out + offsets[2];
  double *           out3     = out + offsets[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);
        }
    }
  if (add_into)
    for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
      for (unsigned int v = 0; v < 4; ++v)
        out[offsets[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[offsets[v] + i] = in[i][v];
}



template <>
class VectorizedArray<float>
{
public:
  static const unsigned int n_array_elements = 8;

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

  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 arithmetics, 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)
  {
#  ifdef __AVX2__
    // 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((const float *)offsets);
    const __m256i index     = *((__m256i *)(&index_val));
    data                    = _mm256_i32gather_ps(base_ptr, index, 4);
#  else
    for (unsigned int i = 0; i < 8; ++i)
      *(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)
      base_ptr[offsets[i]] = *(reinterpret_cast<const float *>(&data) + i);
  }

  __m256 data;

private:
  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;
  }

  template <typename Number2>
  friend VectorizedArray<Number2>
  std::sqrt(const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::abs(const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::max(const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::min(const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
};



template <>
inline void
vectorized_load_and_transpose(const unsigned int      n_entries,
                              const float *           in,
                              const unsigned int *    offsets,
                              VectorizedArray<float> *out)
{
  const unsigned int n_chunks = n_entries / 4;
  for (unsigned int i = 0; i < n_chunks; ++i)
    {
      __m128 u0 = _mm_loadu_ps(in + 4 * i + offsets[0]);
      __m128 u1 = _mm_loadu_ps(in + 4 * i + offsets[1]);
      __m128 u2 = _mm_loadu_ps(in + 4 * i + offsets[2]);
      __m128 u3 = _mm_loadu_ps(in + 4 * i + offsets[3]);
      __m128 u4 = _mm_loadu_ps(in + 4 * i + offsets[4]);
      __m128 u5 = _mm_loadu_ps(in + 4 * i + offsets[5]);
      __m128 u6 = _mm_loadu_ps(in + 4 * i + offsets[6]);
      __m128 u7 = _mm_loadu_ps(in + 4 * i + offsets[7]);
      // To avoid warnings about uninitialized variables, need to initialize
      // one variable with zero before using it.
      __m256 t0, t1, t2, t3 = _mm256_set1_ps(0.F);
      t0                  = _mm256_insertf128_ps(t3, u0, 0);
      t0                  = _mm256_insertf128_ps(t0, u4, 1);
      t1                  = _mm256_insertf128_ps(t3, u1, 0);
      t1                  = _mm256_insertf128_ps(t1, u5, 1);
      t2                  = _mm256_insertf128_ps(t3, u2, 0);
      t2                  = _mm256_insertf128_ps(t2, u6, 1);
      t3                  = _mm256_insertf128_ps(t3, u3, 0);
      t3                  = _mm256_insertf128_ps(t3, u7, 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);
    }
  for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
    for (unsigned int v = 0; v < 8; ++v)
      out[i][v] = in[offsets[v] + i];
}



template <>
inline void
vectorized_transpose_and_store(const bool                    add_into,
                               const unsigned int            n_entries,
                               const VectorizedArray<float> *in,
                               const unsigned int *          offsets,
                               float *                       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 + 4 * i + offsets[0]), res0);
          _mm_storeu_ps(out + 4 * i + offsets[0], res0);
          res1 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[1]), res1);
          _mm_storeu_ps(out + 4 * i + offsets[1], res1);
          res2 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[2]), res2);
          _mm_storeu_ps(out + 4 * i + offsets[2], res2);
          res3 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[3]), res3);
          _mm_storeu_ps(out + 4 * i + offsets[3], res3);
          res4 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[4]), res4);
          _mm_storeu_ps(out + 4 * i + offsets[4], res4);
          res5 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[5]), res5);
          _mm_storeu_ps(out + 4 * i + offsets[5], res5);
          res6 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[6]), res6);
          _mm_storeu_ps(out + 4 * i + offsets[6], res6);
          res7 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[7]), res7);
          _mm_storeu_ps(out + 4 * i + offsets[7], res7);
        }
      else
        {
          _mm_storeu_ps(out + 4 * i + offsets[0], res0);
          _mm_storeu_ps(out + 4 * i + offsets[1], res1);
          _mm_storeu_ps(out + 4 * i + offsets[2], res2);
          _mm_storeu_ps(out + 4 * i + offsets[3], res3);
          _mm_storeu_ps(out + 4 * i + offsets[4], res4);
          _mm_storeu_ps(out + 4 * i + offsets[5], res5);
          _mm_storeu_ps(out + 4 * i + offsets[6], res6);
          _mm_storeu_ps(out + 4 * i + offsets[7], res7);
        }
    }
  if (add_into)
    for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
      for (unsigned int v = 0; v < 8; ++v)
        out[offsets[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[offsets[v] + i] = in[i][v];
}



// for safety, also check that __SSE2__ is defined in case the user manually
// set some conflicting compile flags which prevent compilation

#elif DEAL_II_COMPILER_VECTORIZATION_LEVEL >= 1

template <>
class VectorizedArray<double>
{
public:
  static const unsigned int n_array_elements = 2;

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

  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
  store(double *ptr) const
  {
    _mm_storeu_pd(ptr, data);
  }

  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)
  {
    for (unsigned int i = 0; i < 2; ++i)
      *(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)
      base_ptr[offsets[i]] = *(reinterpret_cast<const double *>(&data) + i);
  }

  __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;
  }

  template <typename Number2>
  friend VectorizedArray<Number2>
  std::sqrt(const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::abs(const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::max(const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::min(const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
};



template <>
inline void
vectorized_load_and_transpose(const unsigned int       n_entries,
                              const double *           in,
                              const unsigned int *     offsets,
                              VectorizedArray<double> *out)
{
  const unsigned int n_chunks = n_entries / 2;
  for (unsigned int i = 0; i < n_chunks; ++i)
    {
      __m128d u0          = _mm_loadu_pd(in + 2 * i + offsets[0]);
      __m128d u1          = _mm_loadu_pd(in + 2 * i + offsets[1]);
      out[2 * i + 0].data = _mm_unpacklo_pd(u0, u1);
      out[2 * i + 1].data = _mm_unpackhi_pd(u0, u1);
    }
  for (unsigned int i = 2 * n_chunks; i < n_entries; ++i)
    for (unsigned int v = 0; v < 2; ++v)
      out[i][v] = in[offsets[v] + i];
}



template <>
inline void
vectorized_transpose_and_store(const bool                     add_into,
                               const unsigned int             n_entries,
                               const VectorizedArray<double> *in,
                               const unsigned int *           offsets,
                               double *                       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 + 2 * i + offsets[0],
                        _mm_add_pd(_mm_loadu_pd(out + 2 * i + offsets[0]),
                                   res0));
          _mm_storeu_pd(out + 2 * i + offsets[1],
                        _mm_add_pd(_mm_loadu_pd(out + 2 * i + offsets[1]),
                                   res1));
        }
      for (unsigned int i = 2 * n_chunks; i < n_entries; ++i)
        for (unsigned int v = 0; v < 2; ++v)
          out[offsets[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 + 2 * i + offsets[0], res0);
          _mm_storeu_pd(out + 2 * i + offsets[1], res1);
        }
      for (unsigned int i = 2 * n_chunks; i < n_entries; ++i)
        for (unsigned int v = 0; v < 2; ++v)
          out[offsets[v] + i] = in[i][v];
    }
}



template <>
class VectorizedArray<float>
{
public:
  static const unsigned int n_array_elements = 4;

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

  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)
  {
    for (unsigned int i = 0; i < 4; ++i)
      *(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)
      base_ptr[offsets[i]] = *(reinterpret_cast<const float *>(&data) + i);
  }

  __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;
  }

  template <typename Number2>
  friend VectorizedArray<Number2>
  std::sqrt(const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::abs(const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::max(const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
  template <typename Number2>
  friend VectorizedArray<Number2>
  std::min(const VectorizedArray<Number2> &, const VectorizedArray<Number2> &);
};



template <>
inline void
vectorized_load_and_transpose(const unsigned int      n_entries,
                              const float *           in,
                              const unsigned int *    offsets,
                              VectorizedArray<float> *out)
{
  const unsigned int n_chunks = n_entries / 4;
  for (unsigned int i = 0; i < n_chunks; ++i)
    {
      __m128 u0           = _mm_loadu_ps(in + 4 * i + offsets[0]);
      __m128 u1           = _mm_loadu_ps(in + 4 * i + offsets[1]);
      __m128 u2           = _mm_loadu_ps(in + 4 * i + offsets[2]);
      __m128 u3           = _mm_loadu_ps(in + 4 * i + offsets[3]);
      __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);
    }
  for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
    for (unsigned int v = 0; v < 4; ++v)
      out[i][v] = in[offsets[v] + i];
}



template <>
inline void
vectorized_transpose_and_store(const bool                    add_into,
                               const unsigned int            n_entries,
                               const VectorizedArray<float> *in,
                               const unsigned int *          offsets,
                               float *                       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 + 4 * i + offsets[0]), u0);
          _mm_storeu_ps(out + 4 * i + offsets[0], u0);
          u1 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[1]), u1);
          _mm_storeu_ps(out + 4 * i + offsets[1], u1);
          u2 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[2]), u2);
          _mm_storeu_ps(out + 4 * i + offsets[2], u2);
          u3 = _mm_add_ps(_mm_loadu_ps(out + 4 * i + offsets[3]), u3);
          _mm_storeu_ps(out + 4 * i + offsets[3], u3);
        }
      else
        {
          _mm_storeu_ps(out + 4 * i + offsets[0], u0);
          _mm_storeu_ps(out + 4 * i + offsets[1], u1);
          _mm_storeu_ps(out + 4 * i + offsets[2], u2);
          _mm_storeu_ps(out + 4 * i + offsets[3], u3);
        }
    }
  if (add_into)
    for (unsigned int i = 4 * n_chunks; i < n_entries; ++i)
      for (unsigned int v = 0; v < 4; ++v)
        out[offsets[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[offsets[v] + i] = in[i][v];
}



#endif // if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0


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

  return true;
}


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

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

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

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

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

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

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

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

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

inline DEAL_II_ALWAYS_INLINE VectorizedArray<float>
                             operator-(const double &u, const VectorizedArray<float> &v)
{
  VectorizedArray<float> tmp;
  tmp = float(u);
  return tmp -= v;
}

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

inline DEAL_II_ALWAYS_INLINE VectorizedArray<float>
                             operator-(const VectorizedArray<float> &v, const double &u)
{
  VectorizedArray<float> tmp;
  tmp = float(u);
  return v - tmp;
}

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

inline DEAL_II_ALWAYS_INLINE VectorizedArray<float>
                             operator*(const double &u, const VectorizedArray<float> &v)
{
  VectorizedArray<float> tmp;
  tmp = float(u);
  return tmp *= v;
}

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

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

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

inline DEAL_II_ALWAYS_INLINE VectorizedArray<float>
                             operator/(const double &u, const VectorizedArray<float> &v)
{
  VectorizedArray<float> tmp;
  tmp = float(u);
  return tmp /= v;
}

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

inline DEAL_II_ALWAYS_INLINE VectorizedArray<float>
                             operator/(const VectorizedArray<float> &v, const double &u)
{
  VectorizedArray<float> tmp;
  tmp = float(u);
  return v / tmp;
}

template <typename Number>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
                             operator+(const VectorizedArray<Number> &u)
{
  return u;
}

template <typename Number>
inline DEAL_II_ALWAYS_INLINE VectorizedArray<Number>
                             operator-(const VectorizedArray<Number> &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>() - u;
}


DEAL_II_NAMESPACE_CLOSE


namespace std
{
  template <typename Number>
  inline ::VectorizedArray<Number>
  sin(const ::VectorizedArray<Number> &x)
  {
    // put values in an array and later read in that array with an unaligned
    // read. This should save some instructions as compared to directly
    // setting the individual elements and also circumvents a compiler
    // optimization bug in gcc-4.6 with SSE2 (see also deal.II developers list
    // from April 2014, topic "matrix_free/step-48 Test").
    Number values[::VectorizedArray<Number>::n_array_elements];
    for (unsigned int i = 0;
         i < ::VectorizedArray<Number>::n_array_elements;
         ++i)
      values[i] = std::sin(x[i]);
    ::VectorizedArray<Number> out;
    out.load(&values[0]);
    return out;
  }



  template <typename Number>
  inline ::VectorizedArray<Number>
  cos(const ::VectorizedArray<Number> &x)
  {
    Number values[::VectorizedArray<Number>::n_array_elements];
    for (unsigned int i = 0;
         i < ::VectorizedArray<Number>::n_array_elements;
         ++i)
      values[i] = std::cos(x[i]);
    ::VectorizedArray<Number> out;
    out.load(&values[0]);
    return out;
  }



  template <typename Number>
  inline ::VectorizedArray<Number>
  tan(const ::VectorizedArray<Number> &x)
  {
    Number values[::VectorizedArray<Number>::n_array_elements];
    for (unsigned int i = 0;
         i < ::VectorizedArray<Number>::n_array_elements;
         ++i)
      values[i] = std::tan(x[i]);
    ::VectorizedArray<Number> out;
    out.load(&values[0]);
    return out;
  }



  template <typename Number>
  inline ::VectorizedArray<Number>
  exp(const ::VectorizedArray<Number> &x)
  {
    Number values[::VectorizedArray<Number>::n_array_elements];
    for (unsigned int i = 0;
         i < ::VectorizedArray<Number>::n_array_elements;
         ++i)
      values[i] = std::exp(x[i]);
    ::VectorizedArray<Number> out;
    out.load(&values[0]);
    return out;
  }



  template <typename Number>
  inline ::VectorizedArray<Number>
  log(const ::VectorizedArray<Number> &x)
  {
    Number values[::VectorizedArray<Number>::n_array_elements];
    for (unsigned int i = 0;
         i < ::VectorizedArray<Number>::n_array_elements;
         ++i)
      values[i] = std::log(x[i]);
    ::VectorizedArray<Number> out;
    out.load(&values[0]);
    return out;
  }



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



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



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



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



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

} // namespace std

#endif