From b19fd2c957d68d9186a7a5596f1aa2fbceeefc24 Mon Sep 17 00:00:00 2001
From: Bruno Turcksin <bruno.turcksin@gmail.com>
Date: Tue, 3 May 2016 17:02:00 -0400
Subject: [PATCH] Move the functions of the internal namespace of
 vector.templates to their own file.

---
 include/deal.II/lac/vector.templates.h | 1048 +----------------------
 include/deal.II/lac/vector_internal.h  | 1068 ++++++++++++++++++++++++
 2 files changed, 1069 insertions(+), 1047 deletions(-)
 create mode 100644 include/deal.II/lac/vector_internal.h

diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h
index 053dda0c78..78f89d2c8d 100644
--- a/include/deal.II/lac/vector.templates.h
+++ b/include/deal.II/lac/vector.templates.h
@@ -25,6 +25,7 @@
 #include <deal.II/base/vectorization.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/vector_internal.h>
 
 #ifdef DEAL_II_WITH_PETSC
 #  include <deal.II/lac/petsc_vector.h>
@@ -44,1053 +45,6 @@
 
 DEAL_II_NAMESPACE_OPEN
 
-
-namespace internal
-{
-  typedef types::global_dof_index size_type;
-
-  template <typename T>
-  bool is_non_negative (const T &t)
-  {
-    return t >= 0;
-  }
-
-
-  template <typename T>
-  bool is_non_negative (const std::complex<T> &)
-  {
-    Assert (false,
-            ExcMessage ("Complex numbers do not have an ordering."));
-
-    return false;
-  }
-
-
-  template <typename T>
-  void print (const T    &t,
-              const char *format)
-  {
-    if (format != 0)
-      std::printf (format, t);
-    else
-      std::printf (" %5.2f", double(t));
-  }
-
-
-
-  template <typename T>
-  void print (const std::complex<T> &t,
-              const char            *format)
-  {
-    if (format != 0)
-      std::printf (format, t.real(), t.imag());
-    else
-      std::printf (" %5.2f+%5.2fi",
-                   double(t.real()), double(t.imag()));
-  }
-
-  // call std::copy, except for in
-  // the case where we want to copy
-  // from std::complex to a
-  // non-complex type
-  template <typename T, typename U>
-  void copy (const T *begin,
-             const T *end,
-             U       *dest)
-  {
-    std::copy (begin, end, dest);
-  }
-
-  template <typename T, typename U>
-  void copy (const std::complex<T> *begin,
-             const std::complex<T> *end,
-             std::complex<U>       *dest)
-  {
-    std::copy (begin, end, dest);
-  }
-
-  template <typename T, typename U>
-  void copy (const std::complex<T> *,
-             const std::complex<T> *,
-             U *)
-  {
-    Assert (false, ExcMessage ("Can't convert a vector of complex numbers "
-                               "into a vector of reals/doubles"));
-  }
-
-
-
-#ifdef DEAL_II_WITH_THREADS
-  /**
-   * This struct takes the loop range from the tbb parallel for loop and
-   * translates it to the actual ranges of the for loop within the vector. It
-   * encodes the grain size but might choose larger values of chunks than the
-   * minimum grain size. The minimum grain size given to tbb is then simple
-   * 1. For affinity reasons, the layout in this loop must be kept in sync
-   * with the respective class for reductions further down.
-   */
-  template <typename Functor>
-  struct TBBForFunctor
-  {
-    TBBForFunctor(Functor &functor,
-                  const size_type vec_size)
-      :
-      functor(functor),
-      vec_size(vec_size)
-    {
-      // set chunk size for sub-tasks
-      const unsigned int gs = internal::Vector::minimum_parallel_grain_size;
-      n_chunks = std::min(static_cast<size_type>(4*MultithreadInfo::n_threads()),
-                          vec_size / gs);
-      chunk_size = vec_size / n_chunks;
-
-      // round to next multiple of 512 (or minimum grain size if that happens
-      // to be smaller). this is advantageous because our accumulation
-      // algorithms favor lengths of a power of 2 due to pairwise summation ->
-      // at most one 'oddly' sized chunk
-      if (chunk_size > 512)
-        chunk_size = ((chunk_size + 511)/512)*512;
-      n_chunks = (vec_size + chunk_size - 1) / chunk_size;
-      AssertIndexRange((n_chunks-1)*chunk_size, vec_size);
-      AssertIndexRange(vec_size, n_chunks*chunk_size+1);
-    };
-
-    void operator() (const tbb::blocked_range<size_type> &range) const
-    {
-      const size_type begin = range.begin()*chunk_size;
-      const size_type end = std::min(range.end()*chunk_size, vec_size);
-      functor(begin, end);
-    }
-
-    Functor &functor;
-    const size_type vec_size;
-    unsigned int n_chunks;
-    size_type chunk_size;
-  };
-#endif
-
-  template <typename Functor>
-  void parallel_for(Functor &functor,
-                    size_type vec_size,
-                    std_cxx11::shared_ptr<parallel::internal::TBBPartitioner> &partitioner)
-  {
-#ifdef DEAL_II_WITH_THREADS
-    // only go to the parallel function in case there are at least 4 parallel
-    // items, otherwise the overhead is too large
-    if (vec_size >= 4*internal::Vector::minimum_parallel_grain_size &&
-        MultithreadInfo::n_threads() > 1)
-      {
-        Assert(partitioner.get() != NULL,
-               ExcInternalError("Unexpected initialization of Vector that does "
-                                "not set the TBB partitioner to a usable state."));
-        std_cxx11::shared_ptr<tbb::affinity_partitioner> tbb_partitioner =
-          partitioner->acquire_one_partitioner();
-
-        TBBForFunctor<Functor> generic_functor(functor, vec_size);
-        tbb::parallel_for (tbb::blocked_range<size_type> (0,
-                                                          generic_functor.n_chunks,
-                                                          1),
-                           generic_functor,
-                           *tbb_partitioner);
-        partitioner->release_one_partitioner(tbb_partitioner);
-      }
-    else if (vec_size > 0)
-      functor(0,vec_size);
-#else
-    functor(0,vec_size);
-    (void)partitioner;
-#endif
-  }
-
-
-  // Define the functors necessary to use SIMD with TBB. we also include the
-  // simple copy and set operations
-
-  template <typename Number>
-  struct Vector_set
-  {
-    Number *dst;
-    Number value;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (value == Number())
-        std::memset (dst+begin,0,(end-begin)*sizeof(Number));
-      else
-        std::fill (dst+begin, dst+end, value);
-    }
-  };
-
-  template <typename Number, typename OtherNumber>
-  struct Vector_copy
-  {
-    const OtherNumber *src;
-    Number *dst;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (types_are_equal<Number,OtherNumber>::value)
-        std::memcpy(dst+begin, src+begin, (end-begin)*sizeof(Number));
-      else
-        {
-          DEAL_II_OPENMP_SIMD_PRAGMA
-          for (typename dealii::Vector<Number>::size_type i=begin; i<end; ++i)
-            dst[i] = src[i];
-        }
-    }
-  };
-
-  template <typename Number>
-  struct Vectorization_multiply_factor
-  {
-    Number *val;
-    Number factor;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
-        {
-          DEAL_II_OPENMP_SIMD_PRAGMA
-          for (size_type i=begin; i<end; ++i)
-            val[i] *= factor;
-        }
-      else
-        {
-          for (size_type i=begin; i<end; ++i)
-            val[i] *= factor;
-        }
-    }
-  };
-
-  template <typename Number>
-  struct Vectorization_add_av
-  {
-    Number *val;
-    Number *v_val;
-    Number factor;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
-        {
-          DEAL_II_OPENMP_SIMD_PRAGMA
-          for (size_type i=begin; i<end; ++i)
-            val[i] += factor*v_val[i];
-        }
-      else
-        {
-          for (size_type i=begin; i<end; ++i)
-            val[i] += factor*v_val[i];
-        }
-    }
-  };
-
-  template <typename Number>
-  struct Vectorization_sadd_xav
-  {
-    Number *val;
-    Number *v_val;
-    Number a;
-    Number x;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
-        {
-          DEAL_II_OPENMP_SIMD_PRAGMA
-          for (size_type i=begin; i<end; ++i)
-            val[i] = x*val[i] + a*v_val[i];
-        }
-      else
-        {
-          for (size_type i=begin; i<end; ++i)
-            val[i] = x*val[i] + a*v_val[i];
-        }
-    }
-  };
-
-  template <typename Number>
-  struct Vectorization_subtract_v
-  {
-    Number *val;
-    Number *v_val;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
-        {
-          DEAL_II_OPENMP_SIMD_PRAGMA
-          for (size_type i=begin; i<end; ++i)
-            val[i] -= v_val[i];
-        }
-      else
-        {
-          for (size_type i=begin; i<end; ++i)
-            val[i] -= v_val[i];
-        }
-    }
-  };
-
-  template <typename Number>
-  struct Vectorization_add_factor
-  {
-    Number *val;
-    Number factor;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
-        {
-          DEAL_II_OPENMP_SIMD_PRAGMA
-          for (size_type i=begin; i<end; ++i)
-            val[i] += factor;
-        }
-      else
-        {
-          for (size_type i=begin; i<end; ++i)
-            val[i] += factor;
-        }
-    }
-  };
-
-  template <typename Number>
-  struct Vectorization_add_v
-  {
-    Number *val;
-    Number *v_val;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
-        {
-          DEAL_II_OPENMP_SIMD_PRAGMA
-          for (size_type i=begin; i<end; ++i)
-            val[i] += v_val[i];
-        }
-      else
-        {
-          for (size_type i=begin; i<end; ++i)
-            val[i] += v_val[i];
-        }
-    }
-  };
-
-  template <typename Number>
-  struct Vectorization_add_avpbw
-  {
-    Number *val;
-    Number *v_val;
-    Number *w_val;
-    Number a;
-    Number b;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
-        {
-          DEAL_II_OPENMP_SIMD_PRAGMA
-          for (size_type i=begin; i<end; ++i)
-            val[i] = val[i] + a*v_val[i] + b*w_val[i];
-        }
-      else
-        {
-          for (size_type i=begin; i<end; ++i)
-            val[i] = val[i] + a*v_val[i] + b*w_val[i];
-        }
-    }
-  };
-
-  template <typename Number>
-  struct Vectorization_sadd_xv
-  {
-    Number *val;
-    Number *v_val;
-    Number x;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
-        {
-          DEAL_II_OPENMP_SIMD_PRAGMA
-          for (size_type i=begin; i<end; ++i)
-            val[i] = x*val[i] + v_val[i];
-        }
-      else
-        {
-          for (size_type i=begin; i<end; ++i)
-            val[i] = x*val[i] + v_val[i];
-        }
-    }
-  };
-
-  template <typename Number>
-  struct Vectorization_sadd_xavbw
-  {
-    Number *val;
-    Number *v_val;
-    Number *w_val;
-    Number x;
-    Number a;
-    Number b;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
-        {
-          DEAL_II_OPENMP_SIMD_PRAGMA
-          for (size_type i=begin; i<end; ++i)
-            val[i] = x*val[i] + a*v_val[i] + b*w_val[i];
-        }
-      else
-        {
-          for (size_type i=begin; i<end; ++i)
-            val[i] = x*val[i] + a*v_val[i] + b*w_val[i];
-        }
-    }
-  };
-
-  template <typename Number>
-  struct Vectorization_scale
-  {
-    Number *val;
-    Number *v_val;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
-        {
-          DEAL_II_OPENMP_SIMD_PRAGMA
-          for (size_type i=begin; i<end; ++i)
-            val[i] *= v_val[i];
-        }
-      else
-        {
-          for (size_type i=begin; i<end; ++i)
-            val[i] *= v_val[i];
-        }
-    }
-  };
-
-  template <typename Number>
-  struct Vectorization_equ_au
-  {
-    Number *val;
-    Number *u_val;
-    Number a;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
-        {
-          DEAL_II_OPENMP_SIMD_PRAGMA
-          for (size_type i=begin; i<end; ++i)
-            val[i] = a*u_val[i];
-        }
-      else
-        {
-          for (size_type i=begin; i<end; ++i)
-            val[i] = a*u_val[i];
-        }
-    }
-  };
-
-  template <typename Number>
-  struct Vectorization_equ_aubv
-  {
-    Number *val;
-    Number *u_val;
-    Number *v_val;
-    Number a;
-    Number b;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
-        {
-          DEAL_II_OPENMP_SIMD_PRAGMA
-          for (size_type i=begin; i<end; ++i)
-            val[i] = a*u_val[i] + b*v_val[i];
-        }
-      else
-        {
-          for (size_type i=begin; i<end; ++i)
-            val[i] = a*u_val[i] + b*v_val[i];
-        }
-    }
-  };
-
-  template <typename Number>
-  struct Vectorization_equ_aubvcw
-  {
-    Number *val;
-    Number *u_val;
-    Number *v_val;
-    Number *w_val;
-    Number a;
-    Number b;
-    Number c;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
-        {
-          DEAL_II_OPENMP_SIMD_PRAGMA
-          for (size_type i=begin; i<end; ++i)
-            val[i] = a*u_val[i] + b*v_val[i] + c*w_val[i];
-        }
-      else
-        {
-          for (size_type i=begin; i<end; ++i)
-            val[i] = a*u_val[i] + b*v_val[i] + c*w_val[i];
-        }
-    }
-  };
-
-  template <typename Number>
-  struct Vectorization_ratio
-  {
-    Number *val;
-    Number *a_val;
-    Number *b_val;
-
-    void operator() (const size_type begin, const size_type end) const
-    {
-      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
-        {
-          DEAL_II_OPENMP_SIMD_PRAGMA
-          for (size_type i=begin; i<end; ++i)
-            val[i] = a_val[i]/b_val[i];
-        }
-      else
-        {
-          for (size_type i=begin; i<end; ++i)
-            val[i] = a_val[i]/b_val[i];
-        }
-    }
-  };
-
-
-
-  // All sums over all the vector entries (l2-norm, inner product, etc.) are
-  // performed with the same code, using a templated operation defined
-  // here. There are always two versions defined, a standard one that covers
-  // most cases and a vectorized one which is only for equal types and float
-  // and double.
-  template <typename Number, typename Number2>
-  struct Dot
-  {
-    static const bool vectorizes = types_are_equal<Number,Number2>::value &&
-                                   (VectorizedArray<Number>::n_array_elements > 1);
-
-    Number
-    operator() (const size_type i) const
-    {
-      return X[i] * Number(numbers::NumberTraits<Number2>::conjugate(Y[i]));
-    }
-
-    VectorizedArray<Number>
-    do_vectorized(const size_type i) const
-    {
-      VectorizedArray<Number> x, y;
-      x.load(X+i);
-      y.load(Y+i);
-      return x * y;
-    }
-
-    const Number  *X;
-    const Number2 *Y;
-  };
-
-  template <typename Number, typename RealType>
-  struct Norm2
-  {
-    static const bool vectorizes = VectorizedArray<Number>::n_array_elements > 1;
-
-    RealType
-    operator() (const size_type i) const
-    {
-      return numbers::NumberTraits<Number>::abs_square(X[i]);
-    }
-
-    VectorizedArray<Number>
-    do_vectorized(const size_type i) const
-    {
-      VectorizedArray<Number> x;
-      x.load(X+i);
-      return x * x;
-    }
-
-    const Number *X;
-  };
-
-  template <typename Number, typename RealType>
-  struct Norm1
-  {
-    static const bool vectorizes = VectorizedArray<Number>::n_array_elements > 1;
-
-    RealType
-    operator() (const size_type i) const
-    {
-      return numbers::NumberTraits<Number>::abs(X[i]);
-    }
-
-    VectorizedArray<Number>
-    do_vectorized(const size_type i) const
-    {
-      VectorizedArray<Number> x;
-      x.load(X+i);
-      return std::abs(x);
-    }
-
-    const Number *X;
-  };
-
-  template <typename Number, typename RealType>
-  struct NormP
-  {
-    static const bool vectorizes = VectorizedArray<Number>::n_array_elements > 1;
-
-    RealType
-    operator() (const size_type i) const
-    {
-      return std::pow(numbers::NumberTraits<Number>::abs(X[i]), p);
-    }
-
-    VectorizedArray<Number>
-    do_vectorized(const size_type i) const
-    {
-      VectorizedArray<Number> x;
-      x.load(X+i);
-      return std::pow(std::abs(x),p);
-    }
-
-    const Number *X;
-    RealType p;
-  };
-
-  template <typename Number>
-  struct MeanValue
-  {
-    static const bool vectorizes = VectorizedArray<Number>::n_array_elements > 1;
-
-    Number
-    operator() (const size_type i) const
-    {
-      return X[i];
-    }
-
-    VectorizedArray<Number>
-    do_vectorized(const size_type i) const
-    {
-      VectorizedArray<Number> x;
-      x.load(X+i);
-      return x;
-    }
-
-    const Number *X;
-  };
-
-  template <typename Number>
-  struct AddAndDot
-  {
-    static const bool vectorizes = VectorizedArray<Number>::n_array_elements > 1;
-
-    Number
-    operator() (const size_type i) const
-    {
-      X[i] += a * V[i];
-      return X[i] * Number(numbers::NumberTraits<Number>::conjugate(W[i]));
-    }
-
-    VectorizedArray<Number>
-    do_vectorized(const size_type i) const
-    {
-      VectorizedArray<Number> x, w, v;
-      x.load(X+i);
-      v.load(V+i);
-      x += a * v;
-      x.store(X+i);
-      // may only load from W after storing in X because the pointers might
-      // point to the same memory
-      w.load(W+i);
-      return x * w;
-    }
-
-    Number *X;
-    const Number *V, *W;
-    Number a;
-  };
-
-
-
-  // this is the main working loop for all vector sums using the templated
-  // operation above. it accumulates the sums using a block-wise summation
-  // algorithm with post-update. this blocked algorithm has been proposed in
-  // a similar form by Castaldo, Whaley and Chronopoulos (SIAM
-  // J. Sci. Comput. 31, 1156-1174, 2008) and we use the smallest possible
-  // block size, 2. Sometimes it is referred to as pairwise summation. The
-  // worst case error made by this algorithm is on the order O(eps *
-  // log2(vec_size)), whereas a naive summation is O(eps * vec_size). Even
-  // though the Kahan summation is even more accurate with an error O(eps)
-  // by carrying along remainders not captured by the main sum, that involves
-  // additional costs which are not worthwhile. See the Wikipedia article on
-  // the Kahan summation algorithm.
-
-  // The algorithm implemented here has the additional benefit that it is
-  // easily parallelized without changing the order of how the elements are
-  // added (floating point addition is not associative). For the same vector
-  // size and minimum_parallel_grainsize, the blocks are always the
-  // same and added pairwise.
-
-  // The depth of recursion is controlled by the 'magic' parameter
-  // vector_accumulation_recursion_threshold: If the length is below
-  // vector_accumulation_recursion_threshold * 32 (32 is the part of code we
-  // unroll), a straight loop instead of recursion will be used.  At the
-  // innermost level, eight values are added consecutively in order to better
-  // balance multiplications and additions.
-
-  // The code returns the result as the last argument in order to make
-  // spawning tasks simpler and use automatic template deduction.
-
-  const unsigned int vector_accumulation_recursion_threshold = 128;
-
-  template <typename Operation, typename ResultType>
-  void accumulate_recursive (const Operation   &op,
-                             const size_type    first,
-                             const size_type    last,
-                             ResultType        &result)
-  {
-    const size_type vec_size = last - first;
-    if (vec_size <= vector_accumulation_recursion_threshold * 32)
-      {
-        // the vector is short enough so we perform the summation. first
-        // work on the regular part. The innermost 32 values are expanded in
-        // order to obtain known loop bounds for most of the work.
-        size_type index = first;
-        ResultType outer_results [vector_accumulation_recursion_threshold];
-        size_type n_chunks = vec_size / 32;
-        const size_type remainder = vec_size % 32;
-        Assert (remainder == 0 || n_chunks < vector_accumulation_recursion_threshold,
-                ExcInternalError());
-
-        // Select between the regular version and vectorized version based
-        // on the number types we are given. To choose the vectorized
-        // version often enough, we need to have all tasks but the last one
-        // to be divisible by the vectorization length
-        accumulate_regular(op, n_chunks, index, outer_results,
-                           internal::bool2type<Operation::vectorizes>());
-
-        // now work on the remainder, i.e., the last up to 32 values. Use
-        // switch statement with fall-through to work on these values.
-        if (remainder > 0)
-          {
-            AssertIndexRange(n_chunks, vector_accumulation_recursion_threshold+1);
-            const size_type inner_chunks = remainder / 8;
-            Assert (inner_chunks <= 3, ExcInternalError());
-            const size_type remainder_inner = remainder % 8;
-            ResultType r0 = ResultType(), r1 = ResultType(),
-                       r2 = ResultType();
-            switch (inner_chunks)
-              {
-              case 3:
-                r2 = op(index++);
-                for (size_type j=1; j<8; ++j)
-                  r2 += op(index++);
-              // no break
-              case 2:
-                r1 = op(index++);
-                for (size_type j=1; j<8; ++j)
-                  r1 += op(index++);
-                r1 += r2;
-              // no break
-              case 1:
-                r2 = op(index++);
-                for (size_type j=1; j<8; ++j)
-                  r2 += op(index++);
-              // no break
-              default:
-                for (size_type j=0; j<remainder_inner; ++j)
-                  r0 += op(index++);
-                r0 += r2;
-                r0 += r1;
-                if (n_chunks == vector_accumulation_recursion_threshold)
-                  outer_results[vector_accumulation_recursion_threshold-1] += r0;
-                else
-                  {
-                    outer_results[n_chunks] = r0;
-                    n_chunks++;
-                  }
-                break;
-              }
-          }
-        AssertDimension(index, last);
-
-        // now sum the results from the chunks
-        // recursively
-        while (n_chunks > 1)
-          {
-            if (n_chunks % 2 == 1)
-              outer_results[n_chunks++] = ResultType();
-            for (size_type i=0; i<n_chunks; i+=2)
-              outer_results[i/2] = outer_results[i] + outer_results[i+1];
-            n_chunks /= 2;
-          }
-        result = outer_results[0];
-      }
-    else
-      {
-        // split vector into four pieces and work on the pieces
-        // recursively. Make pieces (except last) divisible by one fourth the
-        // recursion threshold.
-        const size_type new_size =
-          (vec_size / (vector_accumulation_recursion_threshold * 32)) *
-          vector_accumulation_recursion_threshold * 8;
-        ResultType r0, r1, r2, r3;
-        accumulate_recursive (op, first, first+new_size, r0);
-        accumulate_recursive (op, first+new_size, first+2*new_size, r1);
-        accumulate_recursive (op, first+2*new_size, first+3*new_size, r2);
-        accumulate_recursive (op, first+3*new_size, last, r3);
-        r0 += r1;
-        r2 += r3;
-        result = r0 + r2;
-      }
-  }
-
-
-  // this is the inner working routine for the accumulation loops
-  // below. This is the standard case where the loop bounds are known. We
-  // pulled this function out of the regular accumulate routine because we
-  // might do this thing vectorized (see specialized function below)
-  template <typename Operation, typename ResultType>
-  void
-  accumulate_regular(const Operation &op,
-                     size_type       &n_chunks,
-                     size_type       &index,
-                     ResultType (&outer_results)[vector_accumulation_recursion_threshold],
-                     internal::bool2type<false>)
-  {
-    for (size_type i=0; i<n_chunks; ++i)
-      {
-        ResultType r0 = op(index);
-        ResultType r1 = op(index+1);
-        ResultType r2 = op(index+2);
-        ResultType r3 = op(index+3);
-        index += 4;
-        for (size_type j=1; j<8; ++j, index += 4)
-          {
-            r0 += op(index);
-            r1 += op(index+1);
-            r2 += op(index+2);
-            r3 += op(index+3);
-          }
-        r0 += r1;
-        r2 += r3;
-        outer_results[i] = r0 + r2;
-      }
-  }
-
-
-
-  // this is the inner working routine for the accumulation loops
-  // below. This is the specialized case where the loop bounds are known and
-  // where we can vectorize. In that case, we request the 'do_vectorized'
-  // routine of the operation instead of the regular one which does several
-  // operations at once.
-  template <typename Operation, typename Number>
-  void
-  accumulate_regular(const Operation &op,
-                     size_type       &n_chunks,
-                     size_type       &index,
-                     Number (&outer_results)[vector_accumulation_recursion_threshold],
-                     internal::bool2type<true>)
-  {
-    const unsigned int nvecs = VectorizedArray<Number>::n_array_elements;
-    const size_type regular_chunks = n_chunks/nvecs;
-    for (size_type i=0; i<regular_chunks; ++i)
-      {
-        VectorizedArray<Number> r0 = op.do_vectorized(index);
-        VectorizedArray<Number> r1 = op.do_vectorized(index+nvecs);
-        VectorizedArray<Number> r2 = op.do_vectorized(index+2*nvecs);
-        VectorizedArray<Number> r3 = op.do_vectorized(index+3*nvecs);
-        index += nvecs*4;
-        for (size_type j=1; j<8; ++j, index += nvecs*4)
-          {
-            r0 += op.do_vectorized(index);
-            r1 += op.do_vectorized(index+nvecs);
-            r2 += op.do_vectorized(index+2*nvecs);
-            r3 += op.do_vectorized(index+3*nvecs);
-          }
-        r0 += r1;
-        r2 += r3;
-        r0 += r2;
-        r0.store(&outer_results[i*VectorizedArray<Number>::n_array_elements]);
-      }
-
-    // If we are treating a case where the vector length is not divisible by
-    // the vectorization length, need a cleanup loop
-    AssertIndexRange(VectorizedArray<Number>::n_array_elements,
-                     17);
-    if (n_chunks % VectorizedArray<Number>::n_array_elements != 0)
-      {
-        VectorizedArray<Number> r0 = VectorizedArray<Number>(),
-                                r1 = VectorizedArray<Number>();
-        const size_type start_irreg = regular_chunks * nvecs;
-        for (size_type c=start_irreg; c<n_chunks; ++c)
-          for (size_type j=0; j<32; j+=2*nvecs, index+=2*nvecs)
-            {
-              r0 += op.do_vectorized(index);
-              r1 += op.do_vectorized(index+nvecs);
-            }
-        r0 += r1;
-        r0.store(&outer_results[start_irreg]);
-        n_chunks = start_irreg + VectorizedArray<Number>::n_array_elements;
-      }
-  }
-
-
-
-#ifdef DEAL_II_WITH_THREADS
-  /**
-   * This struct takes the loop range from the tbb parallel for loop and
-   * translates it to the actual ranges of the reduction loop inside the
-   * vector. It encodes the grain size but might choose larger values of
-   * chunks than the minimum grain size. The minimum grain size given to tbb
-   * is 1. For affinity reasons, the layout in this loop must be kept in sync
-   * with the respective class for plain for loops further up.
-   *
-   * Due to this construction, TBB usually only sees a loop of length
-   * 4*num_threads with grain size 1. The actual ranges inside the vector are
-   * computed outside of TBB because otherwise TBB would split the ranges in
-   * some unpredictable position which destroys exact bitwise
-   * reproducibility. An important part of this is that inside
-   * TBBReduceFunctor::operator() the recursive calls to accumulate are done
-   * sequentially on one item a time (even though we could directly run it on
-   * the whole range given through the tbb::blocked_range times the chunk size
-   * - but that would be unpredictable). Thus, the values we cannot control
-   * are the positions in the array that gets filled - but up to that point
-   * the algorithm TBB sees is just a parallel for and nothing unpredictable
-   * can happen.
-   *
-   * To sum up: Once the number of threads and the vector size are fixed, we
-   * have an exact layout of how the calls into the recursive function will
-   * happen. Inside the recursive function, we again only depend on the
-   * length. Finally, the concurrent threads write into different positions in
-   * a result vector in a thread-safe way and the addition in the short array
-   * is again serial.
-   */
-  template <typename Operation, typename ResultType>
-  struct TBBReduceFunctor
-  {
-    static const unsigned int threshold_array_allocate = 512;
-
-    TBBReduceFunctor(const Operation   &op,
-                     const size_type    vec_size)
-      :
-      op(op),
-      vec_size(vec_size)
-    {
-      // set chunk size for sub-tasks
-      const unsigned int gs = internal::Vector::minimum_parallel_grain_size;
-      n_chunks = std::min(static_cast<size_type>(4*MultithreadInfo::n_threads()),
-                          vec_size / gs);
-      chunk_size = vec_size / n_chunks;
-
-      // round to next multiple of 512 (or leave it at the minimum grain size
-      // if that happens to be smaller). this is advantageous because our
-      // algorithm favors lengths of a power of 2 due to pairwise summation ->
-      // at most one 'oddly' sized chunk
-      if (chunk_size > 512)
-        chunk_size = ((chunk_size + 511)/512)*512;
-      n_chunks = (vec_size + chunk_size - 1) / chunk_size;
-      AssertIndexRange((n_chunks-1)*chunk_size, vec_size);
-      AssertIndexRange(vec_size, n_chunks*chunk_size+1);
-
-      if (n_chunks > threshold_array_allocate)
-        {
-          large_array.resize(n_chunks);
-          array_ptr = &large_array[0];
-        }
-      else
-        array_ptr = &small_array[0];
-    };
-
-    void operator() (const tbb::blocked_range<size_type> &range) const
-    {
-      for (size_type i = range.begin(); i < range.end(); ++i)
-        accumulate_recursive(op, i*chunk_size, std::min((i+1)*chunk_size, vec_size),
-                             array_ptr[i]);
-    }
-
-    ResultType do_sum() const
-    {
-      while (n_chunks > 1)
-        {
-          if (n_chunks % 2 == 1)
-            array_ptr[n_chunks++] = ResultType();
-          for (size_type i=0; i<n_chunks; i+=2)
-            array_ptr[i/2] = array_ptr[i] + array_ptr[i+1];
-          n_chunks /= 2;
-        }
-      return array_ptr[0];
-    }
-
-    const Operation &op;
-    const size_type vec_size;
-
-    mutable unsigned int n_chunks;
-    unsigned int chunk_size;
-    ResultType small_array [threshold_array_allocate];
-    std::vector<ResultType> large_array;
-    // this variable either points to small_array or large_array depending on
-    // the number of threads we want to feed
-    mutable ResultType *array_ptr;
-  };
-#endif
-
-
-
-  /**
-   * This is the general caller for parallel reduction operations that work in
-   * parallel.
-   */
-  template <typename Operation, typename ResultType>
-  void parallel_reduce (const Operation   &op,
-                        const size_type    vec_size,
-                        ResultType        &result,
-                        std_cxx11::shared_ptr<parallel::internal::TBBPartitioner> &partitioner)
-  {
-#ifdef DEAL_II_WITH_THREADS
-    // only go to the parallel function in case there are at least 4 parallel
-    // items, otherwise the overhead is too large
-    if (vec_size >= 4*internal::Vector::minimum_parallel_grain_size &&
-        MultithreadInfo::n_threads() > 1)
-      {
-        Assert(partitioner.get() != NULL,
-               ExcInternalError("Unexpected initialization of Vector that does "
-                                "not set the TBB partitioner to a usable state."));
-        std_cxx11::shared_ptr<tbb::affinity_partitioner> tbb_partitioner =
-          partitioner->acquire_one_partitioner();
-
-        TBBReduceFunctor<Operation,ResultType> generic_functor(op, vec_size);
-        tbb::parallel_for (tbb::blocked_range<size_type> (0,
-                                                          generic_functor.n_chunks,
-                                                          1),
-                           generic_functor,
-                           *tbb_partitioner);
-        partitioner->release_one_partitioner(tbb_partitioner);
-        result = generic_functor.do_sum();
-      }
-    else if (vec_size > 0)
-      accumulate_recursive(op,0,vec_size,result);
-#else
-    accumulate_recursive(op,0,vec_size,result);
-    (void)partitioner;
-#endif
-  }
-}
-
-
-
 template <typename Number>
 Vector<Number>::Vector (const Vector<Number> &v)
   :
diff --git a/include/deal.II/lac/vector_internal.h b/include/deal.II/lac/vector_internal.h
new file mode 100644
index 0000000000..eeedf23be4
--- /dev/null
+++ b/include/deal.II/lac/vector_internal.h
@@ -0,0 +1,1068 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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__vector_internal_h
+#define dealii__vector_internal_h
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace internal
+{
+  typedef types::global_dof_index size_type;
+
+  template <typename T>
+  bool is_non_negative (const T &t)
+  {
+    return t >= 0;
+  }
+
+
+  template <typename T>
+  bool is_non_negative (const std::complex<T> &)
+  {
+    Assert (false,
+            ExcMessage ("Complex numbers do not have an ordering."));
+
+    return false;
+  }
+
+
+  template <typename T>
+  void print (const T    &t,
+              const char *format)
+  {
+    if (format != 0)
+      std::printf (format, t);
+    else
+      std::printf (" %5.2f", double(t));
+  }
+
+
+
+  template <typename T>
+  void print (const std::complex<T> &t,
+              const char            *format)
+  {
+    if (format != 0)
+      std::printf (format, t.real(), t.imag());
+    else
+      std::printf (" %5.2f+%5.2fi",
+                   double(t.real()), double(t.imag()));
+  }
+
+  // call std::copy, except for in
+  // the case where we want to copy
+  // from std::complex to a
+  // non-complex type
+  template <typename T, typename U>
+  void copy (const T *begin,
+             const T *end,
+             U       *dest)
+  {
+    std::copy (begin, end, dest);
+  }
+
+  template <typename T, typename U>
+  void copy (const std::complex<T> *begin,
+             const std::complex<T> *end,
+             std::complex<U>       *dest)
+  {
+    std::copy (begin, end, dest);
+  }
+
+  template <typename T, typename U>
+  void copy (const std::complex<T> *,
+             const std::complex<T> *,
+             U *)
+  {
+    Assert (false, ExcMessage ("Can't convert a vector of complex numbers "
+                               "into a vector of reals/doubles"));
+  }
+
+
+
+#ifdef DEAL_II_WITH_THREADS
+  /**
+   * This struct takes the loop range from the tbb parallel for loop and
+   * translates it to the actual ranges of the for loop within the vector. It
+   * encodes the grain size but might choose larger values of chunks than the
+   * minimum grain size. The minimum grain size given to tbb is then simple
+   * 1. For affinity reasons, the layout in this loop must be kept in sync
+   * with the respective class for reductions further down.
+   */
+  template <typename Functor>
+  struct TBBForFunctor
+  {
+    TBBForFunctor(Functor &functor,
+                  const size_type vec_size)
+      :
+      functor(functor),
+      vec_size(vec_size)
+    {
+      // set chunk size for sub-tasks
+      const unsigned int gs = internal::Vector::minimum_parallel_grain_size;
+      n_chunks = std::min(static_cast<size_type>(4*MultithreadInfo::n_threads()),
+                          vec_size / gs);
+      chunk_size = vec_size / n_chunks;
+
+      // round to next multiple of 512 (or minimum grain size if that happens
+      // to be smaller). this is advantageous because our accumulation
+      // algorithms favor lengths of a power of 2 due to pairwise summation ->
+      // at most one 'oddly' sized chunk
+      if (chunk_size > 512)
+        chunk_size = ((chunk_size + 511)/512)*512;
+      n_chunks = (vec_size + chunk_size - 1) / chunk_size;
+      AssertIndexRange((n_chunks-1)*chunk_size, vec_size);
+      AssertIndexRange(vec_size, n_chunks*chunk_size+1);
+    };
+
+    void operator() (const tbb::blocked_range<size_type> &range) const
+    {
+      const size_type begin = range.begin()*chunk_size;
+      const size_type end = std::min(range.end()*chunk_size, vec_size);
+      functor(begin, end);
+    }
+
+    Functor &functor;
+    const size_type vec_size;
+    unsigned int n_chunks;
+    size_type chunk_size;
+  };
+#endif
+
+  template <typename Functor>
+  void parallel_for(Functor &functor,
+                    size_type vec_size,
+                    std_cxx11::shared_ptr<parallel::internal::TBBPartitioner> &partitioner)
+  {
+#ifdef DEAL_II_WITH_THREADS
+    // only go to the parallel function in case there are at least 4 parallel
+    // items, otherwise the overhead is too large
+    if (vec_size >= 4*internal::Vector::minimum_parallel_grain_size &&
+        MultithreadInfo::n_threads() > 1)
+      {
+        Assert(partitioner.get() != NULL,
+               ExcInternalError("Unexpected initialization of Vector that does "
+                                "not set the TBB partitioner to a usable state."));
+        std_cxx11::shared_ptr<tbb::affinity_partitioner> tbb_partitioner =
+          partitioner->acquire_one_partitioner();
+
+        TBBForFunctor<Functor> generic_functor(functor, vec_size);
+        tbb::parallel_for (tbb::blocked_range<size_type> (0,
+                                                          generic_functor.n_chunks,
+                                                          1),
+                           generic_functor,
+                           *tbb_partitioner);
+        partitioner->release_one_partitioner(tbb_partitioner);
+      }
+    else if (vec_size > 0)
+      functor(0,vec_size);
+#else
+    functor(0,vec_size);
+    (void)partitioner;
+#endif
+  }
+
+
+  // Define the functors necessary to use SIMD with TBB. we also include the
+  // simple copy and set operations
+
+  template <typename Number>
+  struct Vector_set
+  {
+    Number *dst;
+    Number value;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (value == Number())
+        std::memset (dst+begin,0,(end-begin)*sizeof(Number));
+      else
+        std::fill (dst+begin, dst+end, value);
+    }
+  };
+
+  template <typename Number, typename OtherNumber>
+  struct Vector_copy
+  {
+    const OtherNumber *src;
+    Number *dst;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (types_are_equal<Number,OtherNumber>::value)
+        std::memcpy(dst+begin, src+begin, (end-begin)*sizeof(Number));
+      else
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (typename dealii::Vector<Number>::size_type i=begin; i<end; ++i)
+            dst[i] = src[i];
+        }
+    }
+  };
+
+  template <typename Number>
+  struct Vectorization_multiply_factor
+  {
+    Number *val;
+    Number factor;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (size_type i=begin; i<end; ++i)
+            val[i] *= factor;
+        }
+      else
+        {
+          for (size_type i=begin; i<end; ++i)
+            val[i] *= factor;
+        }
+    }
+  };
+
+  template <typename Number>
+  struct Vectorization_add_av
+  {
+    Number *val;
+    Number *v_val;
+    Number factor;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (size_type i=begin; i<end; ++i)
+            val[i] += factor*v_val[i];
+        }
+      else
+        {
+          for (size_type i=begin; i<end; ++i)
+            val[i] += factor*v_val[i];
+        }
+    }
+  };
+
+  template <typename Number>
+  struct Vectorization_sadd_xav
+  {
+    Number *val;
+    Number *v_val;
+    Number a;
+    Number x;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (size_type i=begin; i<end; ++i)
+            val[i] = x*val[i] + a*v_val[i];
+        }
+      else
+        {
+          for (size_type i=begin; i<end; ++i)
+            val[i] = x*val[i] + a*v_val[i];
+        }
+    }
+  };
+
+  template <typename Number>
+  struct Vectorization_subtract_v
+  {
+    Number *val;
+    Number *v_val;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (size_type i=begin; i<end; ++i)
+            val[i] -= v_val[i];
+        }
+      else
+        {
+          for (size_type i=begin; i<end; ++i)
+            val[i] -= v_val[i];
+        }
+    }
+  };
+
+  template <typename Number>
+  struct Vectorization_add_factor
+  {
+    Number *val;
+    Number factor;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (size_type i=begin; i<end; ++i)
+            val[i] += factor;
+        }
+      else
+        {
+          for (size_type i=begin; i<end; ++i)
+            val[i] += factor;
+        }
+    }
+  };
+
+  template <typename Number>
+  struct Vectorization_add_v
+  {
+    Number *val;
+    Number *v_val;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (size_type i=begin; i<end; ++i)
+            val[i] += v_val[i];
+        }
+      else
+        {
+          for (size_type i=begin; i<end; ++i)
+            val[i] += v_val[i];
+        }
+    }
+  };
+
+  template <typename Number>
+  struct Vectorization_add_avpbw
+  {
+    Number *val;
+    Number *v_val;
+    Number *w_val;
+    Number a;
+    Number b;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (size_type i=begin; i<end; ++i)
+            val[i] = val[i] + a*v_val[i] + b*w_val[i];
+        }
+      else
+        {
+          for (size_type i=begin; i<end; ++i)
+            val[i] = val[i] + a*v_val[i] + b*w_val[i];
+        }
+    }
+  };
+
+  template <typename Number>
+  struct Vectorization_sadd_xv
+  {
+    Number *val;
+    Number *v_val;
+    Number x;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (size_type i=begin; i<end; ++i)
+            val[i] = x*val[i] + v_val[i];
+        }
+      else
+        {
+          for (size_type i=begin; i<end; ++i)
+            val[i] = x*val[i] + v_val[i];
+        }
+    }
+  };
+
+  template <typename Number>
+  struct Vectorization_sadd_xavbw
+  {
+    Number *val;
+    Number *v_val;
+    Number *w_val;
+    Number x;
+    Number a;
+    Number b;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (size_type i=begin; i<end; ++i)
+            val[i] = x*val[i] + a*v_val[i] + b*w_val[i];
+        }
+      else
+        {
+          for (size_type i=begin; i<end; ++i)
+            val[i] = x*val[i] + a*v_val[i] + b*w_val[i];
+        }
+    }
+  };
+
+  template <typename Number>
+  struct Vectorization_scale
+  {
+    Number *val;
+    Number *v_val;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (size_type i=begin; i<end; ++i)
+            val[i] *= v_val[i];
+        }
+      else
+        {
+          for (size_type i=begin; i<end; ++i)
+            val[i] *= v_val[i];
+        }
+    }
+  };
+
+  template <typename Number>
+  struct Vectorization_equ_au
+  {
+    Number *val;
+    Number *u_val;
+    Number a;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (size_type i=begin; i<end; ++i)
+            val[i] = a*u_val[i];
+        }
+      else
+        {
+          for (size_type i=begin; i<end; ++i)
+            val[i] = a*u_val[i];
+        }
+    }
+  };
+
+  template <typename Number>
+  struct Vectorization_equ_aubv
+  {
+    Number *val;
+    Number *u_val;
+    Number *v_val;
+    Number a;
+    Number b;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (size_type i=begin; i<end; ++i)
+            val[i] = a*u_val[i] + b*v_val[i];
+        }
+      else
+        {
+          for (size_type i=begin; i<end; ++i)
+            val[i] = a*u_val[i] + b*v_val[i];
+        }
+    }
+  };
+
+  template <typename Number>
+  struct Vectorization_equ_aubvcw
+  {
+    Number *val;
+    Number *u_val;
+    Number *v_val;
+    Number *w_val;
+    Number a;
+    Number b;
+    Number c;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (size_type i=begin; i<end; ++i)
+            val[i] = a*u_val[i] + b*v_val[i] + c*w_val[i];
+        }
+      else
+        {
+          for (size_type i=begin; i<end; ++i)
+            val[i] = a*u_val[i] + b*v_val[i] + c*w_val[i];
+        }
+    }
+  };
+
+  template <typename Number>
+  struct Vectorization_ratio
+  {
+    Number *val;
+    Number *a_val;
+    Number *b_val;
+
+    void operator() (const size_type begin, const size_type end) const
+    {
+      if (parallel::internal::EnableOpenMPSimdFor<Number>::value)
+        {
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (size_type i=begin; i<end; ++i)
+            val[i] = a_val[i]/b_val[i];
+        }
+      else
+        {
+          for (size_type i=begin; i<end; ++i)
+            val[i] = a_val[i]/b_val[i];
+        }
+    }
+  };
+
+
+
+  // All sums over all the vector entries (l2-norm, inner product, etc.) are
+  // performed with the same code, using a templated operation defined
+  // here. There are always two versions defined, a standard one that covers
+  // most cases and a vectorized one which is only for equal types and float
+  // and double.
+  template <typename Number, typename Number2>
+  struct Dot
+  {
+    static const bool vectorizes = types_are_equal<Number,Number2>::value &&
+                                   (VectorizedArray<Number>::n_array_elements > 1);
+
+    Number
+    operator() (const size_type i) const
+    {
+      return X[i] * Number(numbers::NumberTraits<Number2>::conjugate(Y[i]));
+    }
+
+    VectorizedArray<Number>
+    do_vectorized(const size_type i) const
+    {
+      VectorizedArray<Number> x, y;
+      x.load(X+i);
+      y.load(Y+i);
+      return x * y;
+    }
+
+    const Number  *X;
+    const Number2 *Y;
+  };
+
+  template <typename Number, typename RealType>
+  struct Norm2
+  {
+    static const bool vectorizes = VectorizedArray<Number>::n_array_elements > 1;
+
+    RealType
+    operator() (const size_type i) const
+    {
+      return numbers::NumberTraits<Number>::abs_square(X[i]);
+    }
+
+    VectorizedArray<Number>
+    do_vectorized(const size_type i) const
+    {
+      VectorizedArray<Number> x;
+      x.load(X+i);
+      return x * x;
+    }
+
+    const Number *X;
+  };
+
+  template <typename Number, typename RealType>
+  struct Norm1
+  {
+    static const bool vectorizes = VectorizedArray<Number>::n_array_elements > 1;
+
+    RealType
+    operator() (const size_type i) const
+    {
+      return numbers::NumberTraits<Number>::abs(X[i]);
+    }
+
+    VectorizedArray<Number>
+    do_vectorized(const size_type i) const
+    {
+      VectorizedArray<Number> x;
+      x.load(X+i);
+      return std::abs(x);
+    }
+
+    const Number *X;
+  };
+
+  template <typename Number, typename RealType>
+  struct NormP
+  {
+    static const bool vectorizes = VectorizedArray<Number>::n_array_elements > 1;
+
+    RealType
+    operator() (const size_type i) const
+    {
+      return std::pow(numbers::NumberTraits<Number>::abs(X[i]), p);
+    }
+
+    VectorizedArray<Number>
+    do_vectorized(const size_type i) const
+    {
+      VectorizedArray<Number> x;
+      x.load(X+i);
+      return std::pow(std::abs(x),p);
+    }
+
+    const Number *X;
+    RealType p;
+  };
+
+  template <typename Number>
+  struct MeanValue
+  {
+    static const bool vectorizes = VectorizedArray<Number>::n_array_elements > 1;
+
+    Number
+    operator() (const size_type i) const
+    {
+      return X[i];
+    }
+
+    VectorizedArray<Number>
+    do_vectorized(const size_type i) const
+    {
+      VectorizedArray<Number> x;
+      x.load(X+i);
+      return x;
+    }
+
+    const Number *X;
+  };
+
+  template <typename Number>
+  struct AddAndDot
+  {
+    static const bool vectorizes = VectorizedArray<Number>::n_array_elements > 1;
+
+    Number
+    operator() (const size_type i) const
+    {
+      X[i] += a * V[i];
+      return X[i] * Number(numbers::NumberTraits<Number>::conjugate(W[i]));
+    }
+
+    VectorizedArray<Number>
+    do_vectorized(const size_type i) const
+    {
+      VectorizedArray<Number> x, w, v;
+      x.load(X+i);
+      v.load(V+i);
+      x += a * v;
+      x.store(X+i);
+      // may only load from W after storing in X because the pointers might
+      // point to the same memory
+      w.load(W+i);
+      return x * w;
+    }
+
+    Number *X;
+    const Number *V, *W;
+    Number a;
+  };
+
+
+
+  // this is the main working loop for all vector sums using the templated
+  // operation above. it accumulates the sums using a block-wise summation
+  // algorithm with post-update. this blocked algorithm has been proposed in
+  // a similar form by Castaldo, Whaley and Chronopoulos (SIAM
+  // J. Sci. Comput. 31, 1156-1174, 2008) and we use the smallest possible
+  // block size, 2. Sometimes it is referred to as pairwise summation. The
+  // worst case error made by this algorithm is on the order O(eps *
+  // log2(vec_size)), whereas a naive summation is O(eps * vec_size). Even
+  // though the Kahan summation is even more accurate with an error O(eps)
+  // by carrying along remainders not captured by the main sum, that involves
+  // additional costs which are not worthwhile. See the Wikipedia article on
+  // the Kahan summation algorithm.
+
+  // The algorithm implemented here has the additional benefit that it is
+  // easily parallelized without changing the order of how the elements are
+  // added (floating point addition is not associative). For the same vector
+  // size and minimum_parallel_grainsize, the blocks are always the
+  // same and added pairwise.
+
+  // The depth of recursion is controlled by the 'magic' parameter
+  // vector_accumulation_recursion_threshold: If the length is below
+  // vector_accumulation_recursion_threshold * 32 (32 is the part of code we
+  // unroll), a straight loop instead of recursion will be used.  At the
+  // innermost level, eight values are added consecutively in order to better
+  // balance multiplications and additions.
+
+  // The code returns the result as the last argument in order to make
+  // spawning tasks simpler and use automatic template deduction.
+
+  const unsigned int vector_accumulation_recursion_threshold = 128;
+
+  template <typename Operation, typename ResultType>
+  void accumulate_recursive (const Operation   &op,
+                             const size_type    first,
+                             const size_type    last,
+                             ResultType        &result)
+  {
+    const size_type vec_size = last - first;
+    if (vec_size <= vector_accumulation_recursion_threshold * 32)
+      {
+        // the vector is short enough so we perform the summation. first
+        // work on the regular part. The innermost 32 values are expanded in
+        // order to obtain known loop bounds for most of the work.
+        size_type index = first;
+        ResultType outer_results [vector_accumulation_recursion_threshold];
+        size_type n_chunks = vec_size / 32;
+        const size_type remainder = vec_size % 32;
+        Assert (remainder == 0 || n_chunks < vector_accumulation_recursion_threshold,
+                ExcInternalError());
+
+        // Select between the regular version and vectorized version based
+        // on the number types we are given. To choose the vectorized
+        // version often enough, we need to have all tasks but the last one
+        // to be divisible by the vectorization length
+        accumulate_regular(op, n_chunks, index, outer_results,
+                           internal::bool2type<Operation::vectorizes>());
+
+        // now work on the remainder, i.e., the last up to 32 values. Use
+        // switch statement with fall-through to work on these values.
+        if (remainder > 0)
+          {
+            AssertIndexRange(n_chunks, vector_accumulation_recursion_threshold+1);
+            const size_type inner_chunks = remainder / 8;
+            Assert (inner_chunks <= 3, ExcInternalError());
+            const size_type remainder_inner = remainder % 8;
+            ResultType r0 = ResultType(), r1 = ResultType(),
+                       r2 = ResultType();
+            switch (inner_chunks)
+              {
+              case 3:
+                r2 = op(index++);
+                for (size_type j=1; j<8; ++j)
+                  r2 += op(index++);
+              // no break
+              case 2:
+                r1 = op(index++);
+                for (size_type j=1; j<8; ++j)
+                  r1 += op(index++);
+                r1 += r2;
+              // no break
+              case 1:
+                r2 = op(index++);
+                for (size_type j=1; j<8; ++j)
+                  r2 += op(index++);
+              // no break
+              default:
+                for (size_type j=0; j<remainder_inner; ++j)
+                  r0 += op(index++);
+                r0 += r2;
+                r0 += r1;
+                if (n_chunks == vector_accumulation_recursion_threshold)
+                  outer_results[vector_accumulation_recursion_threshold-1] += r0;
+                else
+                  {
+                    outer_results[n_chunks] = r0;
+                    n_chunks++;
+                  }
+                break;
+              }
+          }
+        AssertDimension(index, last);
+
+        // now sum the results from the chunks
+        // recursively
+        while (n_chunks > 1)
+          {
+            if (n_chunks % 2 == 1)
+              outer_results[n_chunks++] = ResultType();
+            for (size_type i=0; i<n_chunks; i+=2)
+              outer_results[i/2] = outer_results[i] + outer_results[i+1];
+            n_chunks /= 2;
+          }
+        result = outer_results[0];
+      }
+    else
+      {
+        // split vector into four pieces and work on the pieces
+        // recursively. Make pieces (except last) divisible by one fourth the
+        // recursion threshold.
+        const size_type new_size =
+          (vec_size / (vector_accumulation_recursion_threshold * 32)) *
+          vector_accumulation_recursion_threshold * 8;
+        ResultType r0, r1, r2, r3;
+        accumulate_recursive (op, first, first+new_size, r0);
+        accumulate_recursive (op, first+new_size, first+2*new_size, r1);
+        accumulate_recursive (op, first+2*new_size, first+3*new_size, r2);
+        accumulate_recursive (op, first+3*new_size, last, r3);
+        r0 += r1;
+        r2 += r3;
+        result = r0 + r2;
+      }
+  }
+
+
+  // this is the inner working routine for the accumulation loops
+  // below. This is the standard case where the loop bounds are known. We
+  // pulled this function out of the regular accumulate routine because we
+  // might do this thing vectorized (see specialized function below)
+  template <typename Operation, typename ResultType>
+  void
+  accumulate_regular(const Operation &op,
+                     size_type       &n_chunks,
+                     size_type       &index,
+                     ResultType (&outer_results)[vector_accumulation_recursion_threshold],
+                     internal::bool2type<false>)
+  {
+    for (size_type i=0; i<n_chunks; ++i)
+      {
+        ResultType r0 = op(index);
+        ResultType r1 = op(index+1);
+        ResultType r2 = op(index+2);
+        ResultType r3 = op(index+3);
+        index += 4;
+        for (size_type j=1; j<8; ++j, index += 4)
+          {
+            r0 += op(index);
+            r1 += op(index+1);
+            r2 += op(index+2);
+            r3 += op(index+3);
+          }
+        r0 += r1;
+        r2 += r3;
+        outer_results[i] = r0 + r2;
+      }
+  }
+
+
+
+  // this is the inner working routine for the accumulation loops
+  // below. This is the specialized case where the loop bounds are known and
+  // where we can vectorize. In that case, we request the 'do_vectorized'
+  // routine of the operation instead of the regular one which does several
+  // operations at once.
+  template <typename Operation, typename Number>
+  void
+  accumulate_regular(const Operation &op,
+                     size_type       &n_chunks,
+                     size_type       &index,
+                     Number (&outer_results)[vector_accumulation_recursion_threshold],
+                     internal::bool2type<true>)
+  {
+    const unsigned int nvecs = VectorizedArray<Number>::n_array_elements;
+    const size_type regular_chunks = n_chunks/nvecs;
+    for (size_type i=0; i<regular_chunks; ++i)
+      {
+        VectorizedArray<Number> r0 = op.do_vectorized(index);
+        VectorizedArray<Number> r1 = op.do_vectorized(index+nvecs);
+        VectorizedArray<Number> r2 = op.do_vectorized(index+2*nvecs);
+        VectorizedArray<Number> r3 = op.do_vectorized(index+3*nvecs);
+        index += nvecs*4;
+        for (size_type j=1; j<8; ++j, index += nvecs*4)
+          {
+            r0 += op.do_vectorized(index);
+            r1 += op.do_vectorized(index+nvecs);
+            r2 += op.do_vectorized(index+2*nvecs);
+            r3 += op.do_vectorized(index+3*nvecs);
+          }
+        r0 += r1;
+        r2 += r3;
+        r0 += r2;
+        r0.store(&outer_results[i*VectorizedArray<Number>::n_array_elements]);
+      }
+
+    // If we are treating a case where the vector length is not divisible by
+    // the vectorization length, need a cleanup loop
+    AssertIndexRange(VectorizedArray<Number>::n_array_elements,
+                     17);
+    if (n_chunks % VectorizedArray<Number>::n_array_elements != 0)
+      {
+        VectorizedArray<Number> r0 = VectorizedArray<Number>(),
+                                r1 = VectorizedArray<Number>();
+        const size_type start_irreg = regular_chunks * nvecs;
+        for (size_type c=start_irreg; c<n_chunks; ++c)
+          for (size_type j=0; j<32; j+=2*nvecs, index+=2*nvecs)
+            {
+              r0 += op.do_vectorized(index);
+              r1 += op.do_vectorized(index+nvecs);
+            }
+        r0 += r1;
+        r0.store(&outer_results[start_irreg]);
+        n_chunks = start_irreg + VectorizedArray<Number>::n_array_elements;
+      }
+  }
+
+
+
+#ifdef DEAL_II_WITH_THREADS
+  /**
+   * This struct takes the loop range from the tbb parallel for loop and
+   * translates it to the actual ranges of the reduction loop inside the
+   * vector. It encodes the grain size but might choose larger values of
+   * chunks than the minimum grain size. The minimum grain size given to tbb
+   * is 1. For affinity reasons, the layout in this loop must be kept in sync
+   * with the respective class for plain for loops further up.
+   *
+   * Due to this construction, TBB usually only sees a loop of length
+   * 4*num_threads with grain size 1. The actual ranges inside the vector are
+   * computed outside of TBB because otherwise TBB would split the ranges in
+   * some unpredictable position which destroys exact bitwise
+   * reproducibility. An important part of this is that inside
+   * TBBReduceFunctor::operator() the recursive calls to accumulate are done
+   * sequentially on one item a time (even though we could directly run it on
+   * the whole range given through the tbb::blocked_range times the chunk size
+   * - but that would be unpredictable). Thus, the values we cannot control
+   * are the positions in the array that gets filled - but up to that point
+   * the algorithm TBB sees is just a parallel for and nothing unpredictable
+   * can happen.
+   *
+   * To sum up: Once the number of threads and the vector size are fixed, we
+   * have an exact layout of how the calls into the recursive function will
+   * happen. Inside the recursive function, we again only depend on the
+   * length. Finally, the concurrent threads write into different positions in
+   * a result vector in a thread-safe way and the addition in the short array
+   * is again serial.
+   */
+  template <typename Operation, typename ResultType>
+  struct TBBReduceFunctor
+  {
+    static const unsigned int threshold_array_allocate = 512;
+
+    TBBReduceFunctor(const Operation   &op,
+                     const size_type    vec_size)
+      :
+      op(op),
+      vec_size(vec_size)
+    {
+      // set chunk size for sub-tasks
+      const unsigned int gs = internal::Vector::minimum_parallel_grain_size;
+      n_chunks = std::min(static_cast<size_type>(4*MultithreadInfo::n_threads()),
+                          vec_size / gs);
+      chunk_size = vec_size / n_chunks;
+
+      // round to next multiple of 512 (or leave it at the minimum grain size
+      // if that happens to be smaller). this is advantageous because our
+      // algorithm favors lengths of a power of 2 due to pairwise summation ->
+      // at most one 'oddly' sized chunk
+      if (chunk_size > 512)
+        chunk_size = ((chunk_size + 511)/512)*512;
+      n_chunks = (vec_size + chunk_size - 1) / chunk_size;
+      AssertIndexRange((n_chunks-1)*chunk_size, vec_size);
+      AssertIndexRange(vec_size, n_chunks*chunk_size+1);
+
+      if (n_chunks > threshold_array_allocate)
+        {
+          large_array.resize(n_chunks);
+          array_ptr = &large_array[0];
+        }
+      else
+        array_ptr = &small_array[0];
+    };
+
+    void operator() (const tbb::blocked_range<size_type> &range) const
+    {
+      for (size_type i = range.begin(); i < range.end(); ++i)
+        accumulate_recursive(op, i*chunk_size, std::min((i+1)*chunk_size, vec_size),
+                             array_ptr[i]);
+    }
+
+    ResultType do_sum() const
+    {
+      while (n_chunks > 1)
+        {
+          if (n_chunks % 2 == 1)
+            array_ptr[n_chunks++] = ResultType();
+          for (size_type i=0; i<n_chunks; i+=2)
+            array_ptr[i/2] = array_ptr[i] + array_ptr[i+1];
+          n_chunks /= 2;
+        }
+      return array_ptr[0];
+    }
+
+    const Operation &op;
+    const size_type vec_size;
+
+    mutable unsigned int n_chunks;
+    unsigned int chunk_size;
+    ResultType small_array [threshold_array_allocate];
+    std::vector<ResultType> large_array;
+    // this variable either points to small_array or large_array depending on
+    // the number of threads we want to feed
+    mutable ResultType *array_ptr;
+  };
+#endif
+
+
+
+  /**
+   * This is the general caller for parallel reduction operations that work in
+   * parallel.
+   */
+  template <typename Operation, typename ResultType>
+  void parallel_reduce (const Operation   &op,
+                        const size_type    vec_size,
+                        ResultType        &result,
+                        std_cxx11::shared_ptr<parallel::internal::TBBPartitioner> &partitioner)
+  {
+#ifdef DEAL_II_WITH_THREADS
+    // only go to the parallel function in case there are at least 4 parallel
+    // items, otherwise the overhead is too large
+    if (vec_size >= 4*internal::Vector::minimum_parallel_grain_size &&
+        MultithreadInfo::n_threads() > 1)
+      {
+        Assert(partitioner.get() != NULL,
+               ExcInternalError("Unexpected initialization of Vector that does "
+                                "not set the TBB partitioner to a usable state."));
+        std_cxx11::shared_ptr<tbb::affinity_partitioner> tbb_partitioner =
+          partitioner->acquire_one_partitioner();
+
+        TBBReduceFunctor<Operation,ResultType> generic_functor(op, vec_size);
+        tbb::parallel_for (tbb::blocked_range<size_type> (0,
+                                                          generic_functor.n_chunks,
+                                                          1),
+                           generic_functor,
+                           *tbb_partitioner);
+        partitioner->release_one_partitioner(tbb_partitioner);
+        result = generic_functor.do_sum();
+      }
+    else if (vec_size > 0)
+      accumulate_recursive(op,0,vec_size,result);
+#else
+    accumulate_recursive(op,0,vec_size,result);
+    (void)partitioner;
+#endif
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
-- 
2.39.5