From 98bfdd0a10c38c5a0069674c9812e5809b957303 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Mon, 17 Nov 2014 10:35:29 +0100
Subject: [PATCH] Enable vectorization of vector updates in Chebyshev
 preconditioner

Use DEAL_II_OPENMP_SIMD_PRAGMA annotation. To avoid a bug in gcc, use an additional
boolean in the 'if' statement that avoids checking factor1 == Number().
---
 include/deal.II/lac/precondition.h | 109 ++++++++++++++++++-----------
 1 file changed, 67 insertions(+), 42 deletions(-)

diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h
index b1a322f817..a2daa67096 100644
--- a/include/deal.II/lac/precondition.h
+++ b/include/deal.II/lac/precondition.h
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 1999 - 2013 by the deal.II authors
+// Copyright (C) 1999 - 2014 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -1581,55 +1581,50 @@ namespace internal
         }
     }
 
-    // worker loop for deal.II vectors
+    // worker routine for deal.II vectors. Because of vectorization, we need
+    // to put the loop into an extra structure because the virtual function of
+    // VectorUpdatesRange prevents the compiler from applying vectorization.
     template <typename Number>
-    struct VectorUpdatesRange : public parallel::ParallelForInteger
+    struct VectorUpdater
     {
-      typedef types::global_dof_index size_type;
-
-      VectorUpdatesRange (const size_t  size,
-                          const Number *src,
-                          const Number *matrix_diagonal_inverse,
-                          const bool    start_zero,
-                          const Number  factor1,
-                          const Number  factor2,
-                          Number       *update1,
-                          Number       *update2,
-                          Number       *dst)
+      VectorUpdater (const Number *src,
+                     const Number *matrix_diagonal_inverse,
+                     const bool    start_zero,
+                     const Number  factor1,
+                     const Number  factor2,
+                     Number       *update1,
+                     Number       *update2,
+                     Number       *dst)
         :
         src (src),
         matrix_diagonal_inverse (matrix_diagonal_inverse),
+        do_startup (factor1 == Number()),
         start_zero (start_zero),
         factor1 (factor1),
         factor2 (factor2),
         update1 (update1),
         update2 (update2),
         dst (dst)
-      {
-        if (size < internal::Vector::minimum_parallel_grain_size)
-          apply_to_subrange (0, size);
-        else
-          apply_parallel (0, size,
-                          internal::Vector::minimum_parallel_grain_size);
-      }
-
-      ~VectorUpdatesRange()
       {}
 
-      virtual void
-      apply_to_subrange (const size_t begin,
-                         const size_t end) const
+      void
+      apply_to_subrange (const std::size_t begin,
+                         const std::size_t end) const
       {
-        if (factor1 == Number())
+        const Number factor1 = this->factor1;
+        const Number factor2 = this->factor2;
+        if (do_startup)
           {
             if (start_zero)
-              for (size_type i=begin; i<end; ++i)
+              DEAL_II_OPENMP_SIMD_PRAGMA
+              for (std::size_t i=begin; i<end; ++i)
                 {
                   dst[i] = factor2 * src[i] * matrix_diagonal_inverse[i];
                   update1[i] = -dst[i];
                 }
             else
-              for (size_type i=begin; i<end; ++i)
+              DEAL_II_OPENMP_SIMD_PRAGMA
+              for (std::size_t i=begin; i<end; ++i)
                 {
                   update1[i] = ((update2[i]-src[i]) *
                                 factor2*matrix_diagonal_inverse[i]);
@@ -1637,17 +1632,20 @@ namespace internal
                 }
           }
         else
-          for (size_type i=begin; i<end; ++i)
+          DEAL_II_OPENMP_SIMD_PRAGMA
+          for (std::size_t i=begin; i<end; ++i)
             {
-              const Number update2i = ((update2[i] - src[i]) *
-                                       matrix_diagonal_inverse[i]);
-              update1[i] = factor1 * update1[i] + factor2 * update2i;
-              dst[i] -= update1[i];
+              const Number update =
+                factor1 * update1[i] + factor2 *
+                ((update2[i] - src[i]) * matrix_diagonal_inverse[i]);
+              update1[i] = update;
+              dst[i] -= update;
             }
       }
 
       const Number *src;
       const Number *matrix_diagonal_inverse;
+      const bool do_startup;
       const bool start_zero;
       const Number factor1;
       const Number factor2;
@@ -1656,6 +1654,33 @@ namespace internal
       mutable Number *dst;
     };
 
+    template<typename Number>
+    struct VectorUpdatesRange : public parallel::ParallelForInteger
+    {
+      VectorUpdatesRange(const VectorUpdater<Number> &updater,
+                         const std::size_t size)
+        :
+        updater (updater)
+      {
+        if (size < internal::Vector::minimum_parallel_grain_size)
+          apply_to_subrange (0, size);
+        else
+          apply_parallel (0, size,
+                          internal::Vector::minimum_parallel_grain_size);
+      }
+
+      ~VectorUpdatesRange() {}
+
+      virtual void
+      apply_to_subrange (const std::size_t begin,
+                         const std::size_t end) const
+      {
+        updater.apply_to_subrange(begin, end);
+      }
+
+      const VectorUpdater<Number> &updater;
+    };
+
     // selection for deal.II vector
     template <typename Number>
     inline
@@ -1669,10 +1694,10 @@ namespace internal
                     ::dealii::Vector<Number> &update2,
                     ::dealii::Vector<Number> &dst)
     {
-      VectorUpdatesRange<Number>(src.size(), src.begin(),
-                                 matrix_diagonal_inverse.begin(),
-                                 start_zero, factor1, factor2,
-                                 update1.begin(), update2.begin(), dst.begin());
+      VectorUpdater<Number> upd(src.begin(), matrix_diagonal_inverse.begin(),
+                                start_zero, factor1, factor2,
+                                update1.begin(), update2.begin(), dst.begin());
+      VectorUpdatesRange<Number>(upd, src.size());
     }
 
     // selection for parallel deal.II vector
@@ -1688,10 +1713,10 @@ namespace internal
                     parallel::distributed::Vector<Number> &update2,
                     parallel::distributed::Vector<Number> &dst)
     {
-      VectorUpdatesRange<Number>(src.local_size(), src.begin(),
-                                 matrix_diagonal_inverse.begin(),
-                                 start_zero, factor1, factor2,
-                                 update1.begin(), update2.begin(), dst.begin());
+      VectorUpdater<Number> upd(src.begin(), matrix_diagonal_inverse.begin(),
+                                start_zero, factor1, factor2,
+                                update1.begin(), update2.begin(), dst.begin());
+      VectorUpdatesRange<Number>(upd, src.local_size());
     }
 
     template <typename VECTOR>
-- 
2.39.5