From: Bruno Turcksin <bruno.turcksin@gmail.com>
Date: Fri, 8 Nov 2019 13:34:29 +0000 (-0500)
Subject: Do not hardcode the warp size
X-Git-Tag: v9.2.0-rc1~891^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bc5308e777870cd0a5e1b365b0c5f16208edab0d;p=dealii.git

Do not hardcode the warp size
---

diff --git a/include/deal.II/base/cuda_size.h b/include/deal.II/base/cuda_size.h
index 4080012775..09e12eacbe 100644
--- a/include/deal.II/base/cuda_size.h
+++ b/include/deal.II/base/cuda_size.h
@@ -33,6 +33,11 @@ namespace CUDAWrappers
    * changed depending on the architecture the code is running on.
    */
   constexpr int chunk_size = 1;
+
+  /**
+   * Define the number of threads in a warp.
+   */
+  constexpr int warp_size = 32;
 } // namespace CUDAWrappers
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/include/deal.II/lac/cuda_kernels.h b/include/deal.II/lac/cuda_kernels.h
index 61b1f8fcb9..564a461f80 100644
--- a/include/deal.II/lac/cuda_kernels.h
+++ b/include/deal.II/lac/cuda_kernels.h
@@ -41,6 +41,7 @@ namespace LinearAlgebra
     {
       using ::dealii::CUDAWrappers::block_size;
       using ::dealii::CUDAWrappers::chunk_size;
+      using ::dealii::CUDAWrappers::warp_size;
       using size_type = types::global_dof_index;
 
       /**
diff --git a/include/deal.II/lac/cuda_kernels.templates.h b/include/deal.II/lac/cuda_kernels.templates.h
index 69ce2c7823..8c962b0e93 100644
--- a/include/deal.II/lac/cuda_kernels.templates.h
+++ b/include/deal.II/lac/cuda_kernels.templates.h
@@ -189,38 +189,6 @@ namespace LinearAlgebra
 
 
 
-      template <typename Number, typename Operation>
-      __device__ void
-      reduce_within_warp(volatile Number *result_buffer, size_type local_idx)
-      {
-        if (block_size >= 64)
-          result_buffer[local_idx] =
-            Operation::reduction_op(result_buffer[local_idx],
-                                    result_buffer[local_idx + 32]);
-        if (block_size >= 32)
-          result_buffer[local_idx] =
-            Operation::reduction_op(result_buffer[local_idx],
-                                    result_buffer[local_idx + 16]);
-        if (block_size >= 16)
-          result_buffer[local_idx] =
-            Operation::reduction_op(result_buffer[local_idx],
-                                    result_buffer[local_idx + 8]);
-        if (block_size >= 8)
-          result_buffer[local_idx] =
-            Operation::reduction_op(result_buffer[local_idx],
-                                    result_buffer[local_idx + 4]);
-        if (block_size >= 4)
-          result_buffer[local_idx] =
-            Operation::reduction_op(result_buffer[local_idx],
-                                    result_buffer[local_idx + 2]);
-        if (block_size >= 2)
-          result_buffer[local_idx] =
-            Operation::reduction_op(result_buffer[local_idx],
-                                    result_buffer[local_idx + 1]);
-      }
-
-
-
       template <typename Number, typename Operation>
       __device__ void
       reduce(Number *        result,
@@ -229,7 +197,7 @@ namespace LinearAlgebra
              const size_type global_idx,
              const size_type N)
       {
-        for (size_type s = block_size / 2; s > 32; s = s >> 1)
+        for (size_type s = block_size / 2; s > warp_size; s = s >> 1)
           {
             if (local_idx < s)
               result_buffer[local_idx] =
@@ -238,8 +206,15 @@ namespace LinearAlgebra
             __syncthreads();
           }
 
-        if (local_idx < 32)
-          reduce_within_warp<Number, Operation>(result_buffer, local_idx);
+        if (local_idx < warp_size)
+          {
+            for (size_type s = warp_size; s > 0; s = s >> 1)
+              {
+                result_buffer[local_idx] =
+                  Operation::reduction_op(result_buffer[local_idx],
+                                          result_buffer[local_idx + s]);
+              }
+          }
 
         if (local_idx == 0)
           Operation::atomic_op(result, result_buffer[0]);
diff --git a/include/deal.II/matrix_free/cuda_matrix_free.h b/include/deal.II/matrix_free/cuda_matrix_free.h
index e749c854e2..7b93e43eeb 100644
--- a/include/deal.II/matrix_free/cuda_matrix_free.h
+++ b/include/deal.II/matrix_free/cuda_matrix_free.h
@@ -21,6 +21,7 @@
 
 #ifdef DEAL_II_COMPILER_CUDA_AWARE
 
+#  include <deal.II/base/cuda_size.h>
 #  include <deal.II/base/mpi.h>
 #  include <deal.II/base/quadrature.h>
 #  include <deal.II/base/tensor.h>
@@ -602,13 +603,17 @@ namespace CUDAWrappers
            cells_per_block_shmem(int dim, int fe_degree)
   {
     /* clang-format off */
-    return dim==2 ? (fe_degree==1 ? 32 :
-                     fe_degree==2 ? 8 :
-                     fe_degree==3 ? 4 :
-                     fe_degree==4 ? 4 :
+    // We are limiting the number of threads according to the
+    // following formulas:
+    //  - in 2D: `threads = cells * (k+1)^d <= 4*CUDAWrappers::warp_size`
+    //  - in 3D: `threads = cells * (k+1)^d <= 2*CUDAWrappers::warp_size`
+    return dim==2 ? (fe_degree==1 ? CUDAWrappers::warp_size :    // 128
+                     fe_degree==2 ? CUDAWrappers::warp_size/4 :  //  72
+                     fe_degree==3 ? CUDAWrappers::warp_size/8 :  //  64
+                     fe_degree==4 ? CUDAWrappers::warp_size/8 :  // 100
                      1) :
-           dim==3 ? (fe_degree==1 ? 8 :
-                     fe_degree==2 ? 2 :
+           dim==3 ? (fe_degree==1 ? CUDAWrappers::warp_size/4 :  //  64
+                     fe_degree==2 ? CUDAWrappers::warp_size/16 : //  54
                      1) : 1;
     /* clang-format on */
   }