Improve grainsize of vectors, unify vector loops

diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h
index 8a3bc24aa5d19e17e71151fe7b4a0cc51ff39171..3bd2faca6ed3a2e24307ae63237c3aa601a30993 100644 (file)
--- a/include/deal.II/lac/vector.templates.h
+++ b/include/deal.II/lac/vector.templates.h
@@ -122,10 +122,11 @@ namespace internal
   void vectorized_transform(Functor &functor,
                             size_type vec_size)
   {
-#ifndef DEAL_II_WITH_THREADS
-    functor(0,vec_size);
-#else
-    if (vec_size>internal::Vector::minimum_parallel_grain_size)
+#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)
       {
         tbb::parallel_for (tbb::blocked_range<size_type> (0,
                                                           vec_size,
@@ -134,12 +135,62 @@ namespace internal
                            tbb::auto_partitioner());
       }
     else if (vec_size > 0)
-      functor(0,vec_size);
 #endif
+      functor(0,vec_size);
   }
 
 
-  // Define the functors necessary to use SIMD with TBB.
+  // 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;
+
+#ifdef DEAL_II_WITH_THREADS
+    void operator() (const tbb::blocked_range<size_type> &range) const
+    {
+      operator()(range.begin(),range.end());
+    }
+#endif
+
+    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;
+
+#ifdef DEAL_II_WITH_THREADS
+    void operator() (const tbb::blocked_range<size_type> &range) const
+    {
+      operator()(range.begin(),range.end());
+    }
+#endif
+
+    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
   {
@@ -554,6 +605,26 @@ namespace internal
     }
   };
 
+  // function used in the header file
+  namespace Vector
+  {
+    template <typename T, typename U>
+    void copy_vector (const dealii::Vector<T> &src,
+                      dealii::Vector<U>       &dst)
+    {
+      if (PointerComparison::equal(&src, &dst))
+        return;
+
+      const typename dealii::Vector<T>::size_type vec_size = src.size();
+      const typename dealii::Vector<U>::size_type dst_size = dst.size();
+      if (dst_size != vec_size)
+        dst.reinit (src, true);
+      dealii::internal::Vector_copy<U,T> copier;
+      copier.dst = dst.begin();
+      copier.src = src.begin();
+      internal::vectorized_transform(copier,vec_size);
+    }
+  }
 }
 
 
@@ -761,76 +832,6 @@ Vector<Number>::is_non_negative () const
 
 
 
-namespace internal
-{
-  namespace Vector
-  {
-    template <typename T>
-    void set_subrange (const T            s,
-                       const typename dealii::Vector<T>::size_type begin,
-                       const typename dealii::Vector<T>::size_type end,
-                       dealii::Vector<T> &dst)
-    {
-      if (s == T())
-        std::memset ((dst.begin()+begin),0,(end-begin)*sizeof(T));
-      else
-        std::fill (&*(dst.begin()+begin), &*(dst.begin()+end), s);
-    }
-
-
-    template <typename T>
-    void copy_subrange (const typename dealii::Vector<T>::size_type         begin,
-                        const typename dealii::Vector<T>::size_type         end,
-                        const dealii::Vector<T> &src,
-                        dealii::Vector<T>       &dst)
-    {
-      memcpy(&*(dst.begin()+begin), &*(src.begin()+begin),
-             (end-begin)*sizeof(T));
-    }
-
-
-    template <typename T, typename U>
-    void copy_subrange (const typename dealii::Vector<T>::size_type         begin,
-                        const typename dealii::Vector<T>::size_type         end,
-                        const dealii::Vector<T> &src,
-                        dealii::Vector<U>       &dst)
-    {
-      const T *src_ptr = src.begin();
-      U *dst_ptr = dst.begin();
-      DEAL_II_OPENMP_SIMD_PRAGMA
-      for (typename dealii::Vector<T>::size_type i=begin; i<end; ++i)
-        dst_ptr[i] = src_ptr[i];
-    }
-
-
-    template <typename T, typename U>
-    void copy_vector (const dealii::Vector<T> &src,
-                      dealii::Vector<U>       &dst)
-    {
-      if (PointerComparison::equal(&src, &dst))
-        return;
-
-      const typename dealii::Vector<T>::size_type vec_size = src.size();
-      const typename dealii::Vector<U>::size_type dst_size = dst.size();
-      if (dst_size != vec_size)
-        dst.reinit (vec_size, true);
-      if (vec_size>internal::Vector::minimum_parallel_grain_size)
-        parallel::apply_to_subranges (0U, vec_size,
-                                      std_cxx11::bind(&internal::Vector::template
-                                                      copy_subrange <T,U>,
-                                                      std_cxx11::_1,
-                                                      std_cxx11::_2,
-                                                      std_cxx11::cref(src),
-                                                      std_cxx11::ref(dst)),
-                                      internal::Vector::minimum_parallel_grain_size);
-      else if (vec_size > 0)
-        copy_subrange (0U, vec_size, src, dst);
-    }
-  }
-}
-
-
-
 template <typename Number>
 Vector<Number> &
 Vector<Number>::operator= (const Number s)
@@ -838,14 +839,12 @@ Vector<Number>::operator= (const Number s)
   AssertIsFinite(s);
   if (s != Number())
     Assert (vec_size!=0, ExcEmptyObject());
-  if (vec_size>internal::Vector::minimum_parallel_grain_size)
-    parallel::apply_to_subranges (0U, vec_size,
-                                  std_cxx11::bind(&internal::Vector::template
-                                                  set_subrange<Number>,
-                                                  s, std_cxx11::_1, std_cxx11::_2, std_cxx11::ref(*this)),
-                                  internal::Vector::minimum_parallel_grain_size);
-  else if (vec_size > 0)
-    internal::Vector::set_subrange<Number>(s, 0U, vec_size, *this);
+
+  internal::Vector_set<Number> setter;
+  setter.dst = val;
+  setter.value = s;
+
+  internal::vectorized_transform(setter,vec_size);
 
   return *this;
 }

diff --git a/source/base/parallel.cc b/source/base/parallel.cc
index 0ca6d65756a38ebdde77fbffdc185cd5a7c05161..09a175aa9a82d94f9a2028c08f5f5be931dd2f88 100644 (file)
--- a/source/base/parallel.cc
+++ b/source/base/parallel.cc
@@ -23,33 +23,28 @@ namespace internal
 {
   namespace Vector
   {
-    // set minimum grain size. this value is
-    // roughly in accordance with the curve
-    // in the TBB book (fig 3.2) that shows
-    // run time as a function of grain size
-    // -- there, values from 200 upward are
-    // so that the scheduling overhead
-    // amortizes well (for very large values
-    // in that example, the grain size is too
-    // large to split the work load into
-    // enough chunks and the problem becomes
-    // badly balanced)
-    unsigned int minimum_parallel_grain_size = 1024;
+    // set minimum grain size. this value has been determined by experiments
+    // with the actual dealii::Vector implementation and takes vectorization
+    // that is done inside most functions of dealii::Vector into account
+    // (without vectorization, we would land at around 1000). for smaller
+    // values than this, the scheduling overhead will be significant. if the
+    // value becomes too large, on the other hand, the work load cannot be
+    // split into enough chunks and the problem becomes badly balanced. the
+    // tests are documented at https://github.com/dealii/dealii/issues/2496
+    unsigned int minimum_parallel_grain_size = 4096;
   }
 
 
   namespace SparseMatrix
   {
-    // set this value to 1/4 of the value of
-    // the minimum grain size of
-    // vectors. this rests on the fact that
-    // we have to do a lot more work per row
-    // of a matrix than per element of a
-    // vector. it could possibly be reduced
-    // even further but that doesn't appear
-    // worth it any more for anything but
-    // very small matrices that we don't care
-    // that much about anyway.
+    // set this value to 1/16 of the value of the minimum grain size of
+    // vectors (a factor of 4 because we assume that sparse matrix-vector
+    // products do not vectorize and the other factor of 4 because we expect
+    // at least four entries per row). this rests on the fact that we have to
+    // do a lot more work per row of a matrix than per element of a vector. it
+    // could possibly be reduced even further but that doesn't appear worth it
+    // any more for anything but very small matrices that we don't care that
+    // much about anyway.
     unsigned int minimum_parallel_grain_size = 256;
   }
 }