From 5c7a31cc3edb897ab33d93129a77202425c6fa5e Mon Sep 17 00:00:00 2001
From: Daniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Date: Thu, 25 Apr 2019 18:54:54 +0200
Subject: [PATCH] Deal with empty sizes in CUDA kernels

---
 include/deal.II/base/partitioner.templates.h  | 26 ++++++-------
 .../lac/la_parallel_vector.templates.h        |  9 ++---
 .../deal.II/lac/vector_operations_internal.h  | 38 +++++++++----------
 3 files changed, 36 insertions(+), 37 deletions(-)

diff --git a/include/deal.II/base/partitioner.templates.h b/include/deal.II/base/partitioner.templates.h
index 4d4b55108d..62072f1c02 100644
--- a/include/deal.II/base/partitioner.templates.h
+++ b/include/deal.II/base/partitioner.templates.h
@@ -116,16 +116,16 @@ namespace Utilities
       defined(DEAL_II_WITH_CUDA_AWARE_MPI)
           if (std::is_same<MemorySpaceType, MemorySpace::CUDA>::value)
             {
-              const int n_blocks =
-                1 + (import_indices_plain_dev[i].second - 1) /
-                      (::dealii::CUDAWrappers::chunk_size *
-                       ::dealii::CUDAWrappers::block_size);
+              const auto chunk_size = import_indices_plain_dev[i].second;
+              const int  n_blocks =
+                1 + chunk_size / (::dealii::CUDAWrappers::chunk_size *
+                                  ::dealii::CUDAWrappers::block_size);
               ::dealii::LinearAlgebra::CUDAWrappers::kernel::
                 gather<<<n_blocks, ::dealii::CUDAWrappers::block_size>>>(
                   temp_array_ptr,
                   import_indices_plain_dev[i].first.get(),
                   locally_owned_array.data(),
-                  import_indices_plain_dev[i].second);
+                  chunk_size);
             }
           else
 #    endif
@@ -611,8 +611,8 @@ namespace Utilities
                 {
                   const auto chunk_size = import_indices_plain.second;
                   const int n_blocks =
-                    1 + (chunk_size - 1) / (::dealii::CUDAWrappers::chunk_size *
-                                            ::dealii::CUDAWrappers::block_size);
+                    1 + chunk_size / (::dealii::CUDAWrappers::chunk_size *
+                                      ::dealii::CUDAWrappers::block_size);
                   dealii::LinearAlgebra::CUDAWrappers::kernel::
                     masked_vector_bin_op<Number,
                                          dealii::LinearAlgebra::CUDAWrappers::
@@ -631,8 +631,8 @@ namespace Utilities
                 {
                   const auto chunk_size = import_indices_plain.second;
                   const int n_blocks =
-                    1 + (chunk_size - 1) / (::dealii::CUDAWrappers::chunk_size *
-                                            ::dealii::CUDAWrappers::block_size);
+                    1 + chunk_size / (::dealii::CUDAWrappers::chunk_size *
+                                      ::dealii::CUDAWrappers::block_size);
                   dealii::LinearAlgebra::CUDAWrappers::kernel::
                     masked_vector_bin_op<
                       Number,
@@ -651,8 +651,8 @@ namespace Utilities
                 {
                   const auto chunk_size = import_indices_plain.second;
                   const int n_blocks =
-                    1 + (chunk_size - 1) / (::dealii::CUDAWrappers::chunk_size *
-                                            ::dealii::CUDAWrappers::block_size);
+                    1 + chunk_size / (::dealii::CUDAWrappers::chunk_size *
+                                      ::dealii::CUDAWrappers::block_size);
                   dealii::LinearAlgebra::CUDAWrappers::kernel::
                     masked_vector_bin_op<
                       Number,
@@ -671,8 +671,8 @@ namespace Utilities
                 {
                   const auto chunk_size = import_indices_plain.second;
                   const int n_blocks =
-                    1 + (chunk_size - 1) / (::dealii::CUDAWrappers::chunk_size *
-                                            ::dealii::CUDAWrappers::block_size);
+                    1 + chunk_size / (::dealii::CUDAWrappers::chunk_size *
+                                      ::dealii::CUDAWrappers::block_size);
                   dealii::LinearAlgebra::CUDAWrappers::kernel::
                     set_permutated<<<n_blocks,
                                      dealii::CUDAWrappers::block_size>>>(
diff --git a/include/deal.II/lac/la_parallel_vector.templates.h b/include/deal.II/lac/la_parallel_vector.templates.h
index 22bfea5531..ff2ce72f83 100644
--- a/include/deal.II/lac/la_parallel_vector.templates.h
+++ b/include/deal.II/lac/la_parallel_vector.templates.h
@@ -297,8 +297,8 @@ namespace LinearAlgebra
 
           // Set the values in tmp_vector
           const int n_blocks =
-            1 + (n_elements - 1) / (::dealii::CUDAWrappers::chunk_size *
-                                    ::dealii::CUDAWrappers::block_size);
+            1 + n_elements / (::dealii::CUDAWrappers::chunk_size *
+                              ::dealii::CUDAWrappers::block_size);
           ::dealii::LinearAlgebra::CUDAWrappers::kernel::set_permutated<Number>
             <<<n_blocks, ::dealii::CUDAWrappers::block_size>>>(
               indices_dev, tmp_vector.begin(), V_dev, n_elements);
@@ -351,9 +351,8 @@ namespace LinearAlgebra
           AssertCuda(error_code);
           error_code = cudaMemset(result_device, 0, sizeof(Number));
 
-          const int n_blocks =
-            1 + (size - 1) / (::dealii::CUDAWrappers::chunk_size *
-                              ::dealii::CUDAWrappers::block_size);
+          const int n_blocks = 1 + size / (::dealii::CUDAWrappers::chunk_size *
+                                           ::dealii::CUDAWrappers::block_size);
           ::dealii::LinearAlgebra::CUDAWrappers::kernel::reduction<
             Number,
             ::dealii::LinearAlgebra::CUDAWrappers::kernel::LInfty<Number>>
diff --git a/include/deal.II/lac/vector_operations_internal.h b/include/deal.II/lac/vector_operations_internal.h
index 47bef5ac61..ab1a9b954f 100644
--- a/include/deal.II/lac/vector_operations_internal.h
+++ b/include/deal.II/lac/vector_operations_internal.h
@@ -2026,7 +2026,7 @@ namespace internal
                                                  ::dealii::MemorySpace::CUDA>
             &data)
       {
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::set<Number>
           <<<n_blocks, block_size>>>(data.values_dev.get(), s, size);
 
@@ -2048,7 +2048,7 @@ namespace internal
                                                ::dealii::MemorySpace::CUDA>
           &data)
       {
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::add_aV<Number>
           <<<n_blocks, block_size>>>(data.values_dev.get(),
                                      1.,
@@ -2073,7 +2073,7 @@ namespace internal
                                                ::dealii::MemorySpace::CUDA>
           &data)
       {
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::add_aV<Number>
           <<<n_blocks, block_size>>>(data.values_dev.get(),
                                      -1.,
@@ -2097,7 +2097,7 @@ namespace internal
                                                ::dealii::MemorySpace::CUDA>
           &data)
       {
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::vec_add<Number>
           <<<n_blocks, block_size>>>(data.values_dev.get(), a, size);
 
@@ -2120,7 +2120,7 @@ namespace internal
                                                ::dealii::MemorySpace::CUDA>
           &data)
       {
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::add_aV<Number>
           <<<n_blocks, block_size>>>(data.values_dev.get(),
                                      a,
@@ -2149,7 +2149,7 @@ namespace internal
                                                ::dealii::MemorySpace::CUDA>
           &data)
       {
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::add_aVbW<Number>
           <<<dim3(n_blocks, 1), dim3(block_size)>>>(data.values_dev.get(),
                                                     a,
@@ -2177,7 +2177,7 @@ namespace internal
                                                ::dealii::MemorySpace::CUDA>
           &data)
       {
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::sadd<Number>
           <<<dim3(n_blocks, 1), dim3(block_size)>>>(
             x, data.values_dev.get(), 1., v_data.values_dev.get(), size);
@@ -2202,7 +2202,7 @@ namespace internal
                                                ::dealii::MemorySpace::CUDA>
           &data)
       {
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::sadd<Number>
           <<<dim3(n_blocks, 1), dim3(block_size)>>>(
             x, data.values_dev.get(), a, v_data.values_dev.get(), size);
@@ -2230,7 +2230,7 @@ namespace internal
                                                ::dealii::MemorySpace::CUDA>
           &data)
       {
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::sadd<Number>
           <<<dim3(n_blocks, 1), dim3(block_size)>>>(x,
                                                     data.values_dev.get(),
@@ -2257,7 +2257,7 @@ namespace internal
                                                ::dealii::MemorySpace::CUDA>
           &data)
       {
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::vec_scale<Number>
           <<<n_blocks, block_size>>>(data.values_dev.get(), factor, size);
 
@@ -2279,7 +2279,7 @@ namespace internal
                                                ::dealii::MemorySpace::CUDA>
           &data)
       {
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::scale<Number>
           <<<dim3(n_blocks, 1), dim3(block_size)>>>(data.values_dev.get(),
                                                     v_data.values_dev.get(),
@@ -2304,7 +2304,7 @@ namespace internal
                                                ::dealii::MemorySpace::CUDA>
           &data)
       {
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::equ<Number>
           <<<dim3(n_blocks, 1), dim3(block_size)>>>(data.values_dev.get(),
                                                     a,
@@ -2333,7 +2333,7 @@ namespace internal
                                                ::dealii::MemorySpace::CUDA>
           &data)
       {
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::equ<Number>
           <<<dim3(n_blocks, 1), dim3(block_size)>>>(data.values_dev.get(),
                                                     a,
@@ -2360,12 +2360,12 @@ namespace internal
             &data)
       {
         Number *    result_device;
-        cudaError_t error_code =
-          cudaMalloc(&result_device, size * sizeof(Number));
+        cudaError_t error_code = cudaMalloc(&result_device, sizeof(Number));
         AssertCuda(error_code);
         error_code = cudaMemset(result_device, 0, sizeof(Number));
+        AssertCuda(error_code);
 
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::double_vector_reduction<
           Number,
           ::dealii::LinearAlgebra::CUDAWrappers::kernel::DotProduct<Number>>
@@ -2416,7 +2416,7 @@ namespace internal
         AssertCuda(error_code);
         error_code = cudaMemset(result_device, 0, sizeof(Number));
 
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::reduction<
           Number,
           ::dealii::LinearAlgebra::CUDAWrappers::kernel::ElemSum<Number>>
@@ -2453,7 +2453,7 @@ namespace internal
         AssertCuda(error_code);
         error_code = cudaMemset(result_device, 0, sizeof(Number));
 
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::reduction<
           Number,
           ::dealii::LinearAlgebra::CUDAWrappers::kernel::L1Norm<Number>>
@@ -2504,7 +2504,7 @@ namespace internal
         error_code = cudaMemset(res_d, 0, sizeof(Number));
         AssertCuda(error_code);
 
-        const int n_blocks = 1 + (size - 1) / (chunk_size * block_size);
+        const int n_blocks = 1 + size / (chunk_size * block_size);
         ::dealii::LinearAlgebra::CUDAWrappers::kernel::add_and_dot<Number>
           <<<dim3(n_blocks, 1), dim3(block_size)>>>(res_d,
                                                     data.values_dev.get(),
-- 
2.39.5