Add the first vectorized add and add flags to use OpenMP SIMD pragmas.

diff --git a/cmake/setup_compiler_flags_gnu.cmake b/cmake/setup_compiler_flags_gnu.cmake
index 6fb9d81ab9698764e83aee691d9d91c223382897..f5017fd2cfa3fe4fa97241527762185f5776f697 100644 (file)
--- a/cmake/setup_compiler_flags_gnu.cmake
+++ b/cmake/setup_compiler_flags_gnu.cmake
@@ -85,6 +85,13 @@ ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-Wno-long-long")
 ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-Wno-deprecated")
 ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-Wno-deprecated-declarations")
 
+#
+# OpenMP 4.0 can be used for vectorization. Only the vectorization instructions
+# are allowed, the threading must done through TBB. Disable unknown-pragmas
+# warning, if OpenMP 4.0 is not supported.
+#
+ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-fopenmp-simd")
+ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-Wno-unkown-pragmas")
 
 IF(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
   #

diff --git a/cmake/setup_compiler_flags_intel.cmake b/cmake/setup_compiler_flags_intel.cmake
index 33f19fc8fc195a01290fe733a95a944a1bb95fd7..4b8080cd5e6b606e5fbb5cc49e4b7c1f17d5a56e 100644 (file)
--- a/cmake/setup_compiler_flags_intel.cmake
+++ b/cmake/setup_compiler_flags_intel.cmake
@@ -119,6 +119,7 @@ ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-wd21")
 ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-wd2536")
 
 
+#
 # Also disable the following warnings that we frequently
 # trigger writing dimension independent code:
 #   -w111 statement is unreachable
@@ -130,12 +131,22 @@ ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-wd2536")
 #         that is executed only for one specific dimension
 #   -w280 selector expression is constant
 #         When writing 'switch(dim)'
+#
 ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-wd111")
 ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-wd128")
 ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-wd185")
 ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-wd280")
 
 
+#
+# OpenMP 4.0 can be used for vectorization. Only the vectorization instructions
+# are allowed, the threading must done through TBB. Disable unknown-pragmas
+# warning, if OpenMP 4.0 is not supported.
+#
+ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-openmp-simd")
+ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-Wno-unkown-pragmas")
+
+
 
 IF(DEAL_II_STATIC_EXECUTABLE)
   #

diff --git a/cmake/setup_compiler_flags_msvc.cmake b/cmake/setup_compiler_flags_msvc.cmake
index 7529844b4bd13671adccdd337b30b2526024fcb5..db7fb77ef861f4db32a5ff7c76a4da0fc8716322 100644 (file)
--- a/cmake/setup_compiler_flags_msvc.cmake
+++ b/cmake/setup_compiler_flags_msvc.cmake
@@ -45,6 +45,9 @@ ADD_FLAGS(DEAL_II_CXX_FLAGS "/W3")
 # Globally disable some legacy min and max macros that cause problems:
 ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "/NOMINMAX")
 
+# Disable warning about unknown pragmas
+ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "/wd4068")
+
 #############################
 #                           #
 #    For Release target:    #

diff --git a/include/deal.II/lac/vector.h b/include/deal.II/lac/vector.h
index 6d41c07df919988ea1c58512690651b1425fa175..59205f0dc252c87d55d094aea6cf16973af66c9a 100644 (file)
--- a/include/deal.II/lac/vector.h
+++ b/include/deal.II/lac/vector.h
@@ -586,7 +586,7 @@ public:
   Vector<Number> &operator -= (const Vector<Number> &V);
 
   /**
-   * A collective add operation: This funnction adds a whole set of values
+   * A collective add operation: This function adds a whole set of values
    * stored in @p values to the vector components specified by @p indices.
    */
   template <typename OtherNumber>
@@ -627,11 +627,13 @@ public:
   void add (const Vector<Number> &V);
 
   /**
-   * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
+   * Simple vectorized  (SIMD) vector addition, equal to the <tt>operator +=</tt>. 
+   * This function should not be used to add the vector to itself, i.e to
+   * multiply the current vector by two.
    *
    * @dealiiOperationIsMultithreaded
    */
-  void add (const Number a, const Vector<Number> &V);
+  void vectorized_add (const Vector<Number> &V);
 
   /**
    * Multiple addition of scaled vectors, i.e. <tt>*this += a*V+b*W</tt>.
@@ -641,6 +643,13 @@ public:
   void add (const Number a, const Vector<Number> &V,
             const Number b, const Vector<Number> &W);
 
+  /**
+   * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
+   *
+   * @dealiiOperationIsMultithreaded
+   */
+  void add (const Number a, const Vector<Number> &V);
+
   /**
    * Scaling and simple vector addition, i.e.  <tt>*this = s*(*this)+V</tt>.
    *

diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h
index 881dd58c5c0e1b8f19d9d0960e4312bd45f1225b..0289f655b73e0082bcfa6783d4bff0030790e209 100644 (file)
--- a/include/deal.II/lac/vector.templates.h
+++ b/include/deal.II/lac/vector.templates.h
@@ -47,6 +47,8 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace internal
 {
+  typedef types::global_dof_index size_type;
+
   template <typename T>
   bool is_non_negative (const T &t)
   {
@@ -115,6 +117,22 @@ namespace internal
     Assert (false, ExcMessage ("Can't convert a vector of complex numbers "
                                "into a vector of reals/doubles"));
   }
+
+#ifdef DEAL_II_WITH_THREADS
+  template <typename Number>
+  struct Vectorization_add_1
+  {
+    Number* val;
+    Number* v_val;
+    void operator() (const tbb::blocked_range<size_type> &range) const
+    {
+#pragma omp simd
+      for (size_type i=range.begin(); i!=range.end(); ++i)
+        val[i] += v_val[i];
+    }
+  };
+#endif
+   
 }
 
 
@@ -992,6 +1010,36 @@ void Vector<Number>::add (const Number a, const Vector<Number> &v,
 }
 
 
+template <typename Number>
+void Vector<Number>::vectorized_add (const Vector<Number> &v)
+{
+  Assert (vec_size!=0, ExcEmptyObject());
+  Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size));
+
+#ifndef DEAL_II_WITH_THREADS
+#pragma omp simd
+  for (size_type i=0; i<vec_size; ++i)
+    val[i] += v.val[i];
+#else
+  if (vec_size>internal::Vector::minimum_parallel_grain_size)
+  {
+    internal::Vectorization_add_1<Number> vector_add;
+    vector_add.val = val;
+    vector_add.v_val = v.val;
+    tbb::parallel_for (tbb::blocked_range<size_type> (0, 
+                                                      vec_size, 
+                                                      internal::Vector::minimum_parallel_grain_size),
+                       vector_add,
+                       tbb::auto_partitioner());
+  }
+  else if (vec_size > 0)
+#pragma omp simd
+    for (size_type i=0; i<vec_size; ++i)
+      val[i] += v.val[i];
+#endif
+}
+
+
 template <typename Number>
 void Vector<Number>::sadd (const Number x,
                            const Vector<Number> &v)