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 40800127755e5c250097139b95ba4285d2136a98..09e12eacbeff2a83d7909416c514c245df4e3504 100644 (file)
--- 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 61b1f8fcb98fa0e022cc8327220ba06c82d24dce..564a461f800fc14a90129dec0d6810922bb052f4 100644 (file)
--- 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 69ce2c7823f80e46996ed4198c1941fd5ea32bd2..8c962b0e93b84f5a715ddf99ecb73659996ec15e 100644 (file)
--- 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 e749c854e2e79a6d52d77068ca3eabd98cb66571..7b93e43eebfcbc08152415b4df3411db10493aa3 100644 (file)
--- 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 */
   }