vectorization.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2011 - 2017 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------


#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 // 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 VectorizedArray<T> value (const T &t)
    {
      VectorizedArray<T> tmp;
      tmp=t;
      return tmp;
    }
  };
}


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



#ifndef DEAL_II_WITH_CXX11
// Specify the types for the implemented multiplications explicitly

template <typename Number>
struct ProductType<Number, VectorizedArray<Number> >
{
  typedef VectorizedArray<Number> type;
};

template <typename Number>
struct ProductType<VectorizedArray<Number>, Number>
{
  typedef VectorizedArray<Number> type;
};

// In contrast to scalar types for which the product of a float and a double
// variable would be a double variable, the implemented type here really is
// VectorizedArray<float>. Since VectorizedArray<double> is only half as
// wide as VectorizedArray<float>, we would have to throw away half of the
// vector otherwise.
template<>
struct ProductType<double, VectorizedArray<float> >
{
  typedef VectorizedArray<float> type;
};

template<>
struct ProductType<VectorizedArray<float>, double>
{
  typedef VectorizedArray<float> type;
};
#endif



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 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 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 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 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
  {
    for (unsigned int i=0; i<4; ++i)
      for (unsigned int j=i+1; j<4; ++j)
        Assert(offsets[i] != offsets[j],
               ExcMessage("Result of scatter undefined if two offset elements"
                          " point to the same position"));

    // 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 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
  {
    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"));

    // 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 && defined(__SSE2__)

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 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
  {
    Assert(offsets[0] != offsets[1],
           ExcMessage("Result of scatter undefined if two offset elements"
                      " point to the same position"));

    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 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)
      for (unsigned int j=i+1; j<4; ++j)
        Assert(offsets[i] != offsets[j],
               ExcMessage("Result of scatter undefined if two offset elements"
                          " point to the same position"));

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

}

#endif