Allow to use atomic operations instead of graph coloring in

CUDAWrappers::MatrixFree

diff --git a/include/deal.II/matrix_free/cuda_fe_evaluation.h b/include/deal.II/matrix_free/cuda_fe_evaluation.h
index 755b6fdcda6692e48190d1420fe4f6dc224a9288..2bbc299a5cbec693cb06e5868888f25fed8c54cc 100644 (file)
--- a/include/deal.II/matrix_free/cuda_fe_evaluation.h
+++ b/include/deal.II/matrix_free/cuda_fe_evaluation.h
@@ -23,6 +23,7 @@
 #  include <deal.II/base/tensor.h>
 #  include <deal.II/base/utilities.h>
 
+#  include <deal.II/lac/cuda_atomic.h>
 #  include <deal.II/lac/cuda_vector.h>
 
 #  include <deal.II/matrix_free/cuda_hanging_nodes_internal.h>
@@ -199,6 +200,8 @@ namespace CUDAWrappers
 
     const unsigned int constraint_mask;
 
+    const bool use_coloring;
+
     Number *inv_jac;
     Number *JxW;
 
@@ -222,6 +225,7 @@ namespace CUDAWrappers
     : n_cells(data->n_cells)
     , padding_length(data->padding_length)
     , constraint_mask(data->constraint_mask[cell_id])
+    , use_coloring(data->use_coloring)
     , values(shdata->values)
   {
     local_to_global = data->local_to_global + padding_length * cell_id;
@@ -282,7 +286,11 @@ namespace CUDAWrappers
       (dim > 2 ? threadIdx.z : 0) * n_q_points_1d * n_q_points_1d;
     const types::global_dof_index destination_idx = local_to_global[idx];
 
-    dst[destination_idx] += values[idx];
+    if (use_coloring)
+      dst[destination_idx] += values[idx];
+    else
+      LinearAlgebra::CUDAWrappers::atomicAdd_wrapper(&dst[destination_idx],
+                                                     values[idx]);
   }
 
 

diff --git a/include/deal.II/matrix_free/cuda_matrix_free.h b/include/deal.II/matrix_free/cuda_matrix_free.h
index 2015252c16a51208f45622b14c9b1a47e99e05a1..d1e07eeaa080d6475415b32b39373504ea11d6d7 100644 (file)
--- a/include/deal.II/matrix_free/cuda_matrix_free.h
+++ b/include/deal.II/matrix_free/cuda_matrix_free.h
@@ -91,14 +91,19 @@ namespace CUDAWrappers
       parallel_over_elem
     };
 
+    /**
+     * Standardized data struct to pipe additional data to MatrixFree.
+     */
     struct AdditionalData
     {
       AdditionalData(
         const ParallelizationScheme parallelization_scheme = parallel_in_elem,
         const UpdateFlags           mapping_update_flags   = update_gradients |
-                                                 update_JxW_values)
+                                                 update_JxW_values,
+        const bool use_coloring = false)
         : parallelization_scheme(parallelization_scheme)
         , mapping_update_flags(mapping_update_flags)
+        , use_coloring(use_coloring)
       {}
 
       /**
@@ -120,6 +125,13 @@ namespace CUDAWrappers
        * must be specified by this field.
        */
       UpdateFlags mapping_update_flags;
+
+      /**
+       * If true, use graph coloring. Otherwise, use atomic operations. Graph
+       * coloring ensures bitwise reproducibility but is slower on Pascal and
+       * newer architectures.
+       */
+      bool use_coloring;
     };
 
     /**
@@ -136,6 +148,7 @@ namespace CUDAWrappers
       unsigned int             padding_length;
       unsigned int             row_start;
       unsigned int *           constraint_mask;
+      bool                     use_coloring;
     };
 
     /**
@@ -377,6 +390,13 @@ namespace CUDAWrappers
      */
     ParallelizationScheme parallelization_scheme;
 
+    /**
+     * If true, use graph coloring. Otherwise, use atomic operations. Graph
+     * coloring ensures bitwise reproducibility but is slower on Pascal and
+     * newer architectures.
+     */
+    bool use_coloring;
+
     /**
      * Total number of degrees of freedom.
      */

diff --git a/include/deal.II/matrix_free/cuda_matrix_free.templates.h b/include/deal.II/matrix_free/cuda_matrix_free.templates.h
index 2be980af56920322e3a04351cef9e617e07b577c..561f5b18619a7eab7b4f9a69b738afe69d1c01f3 100644 (file)
--- a/include/deal.II/matrix_free/cuda_matrix_free.templates.h
+++ b/include/deal.II/matrix_free/cuda_matrix_free.templates.h
@@ -588,6 +588,7 @@ namespace CUDAWrappers
     data_copy.n_cells         = n_cells[color];
     data_copy.padding_length  = padding_length;
     data_copy.row_start       = row_start[color];
+    data_copy.use_coloring    = use_coloring;
 
     return data_copy;
   }
@@ -803,6 +804,7 @@ namespace CUDAWrappers
       AssertThrow(false, ExcMessage("Invalid parallelization scheme."));
 
     this->parallelization_scheme = additional_data.parallelization_scheme;
+    this->use_coloring           = additional_data.use_coloring;
 
     // TODO: only free if we actually need arrays of different length
     free();
@@ -861,13 +863,27 @@ namespace CUDAWrappers
     CellFilter begin(IteratorFilters::LocallyOwnedCell(),
                      dof_handler.begin_active());
     CellFilter end(IteratorFilters::LocallyOwnedCell(), dof_handler.end());
-    const auto fun = [&](const CellFilter &filter) {
-      return internal::get_conflict_indices<dim, Number>(filter, constraints);
-    };
-
     std::vector<std::vector<CellFilter>> graph;
+
     if (begin != end)
-      graph = GraphColoring::make_graph_coloring(begin, end, fun);
+      {
+        if (additional_data.use_coloring)
+          {
+            const auto fun = [&](const CellFilter &filter) {
+              return internal::get_conflict_indices<dim, Number>(filter,
+                                                                 constraints);
+            };
+            graph = GraphColoring::make_graph_coloring(begin, end, fun);
+          }
+        else
+          {
+            // If we are not using coloring, all the cells belong to the same
+            // color.
+            graph.resize(1, std::vector<CellFilter>());
+            for (auto cell = begin; cell != end; ++cell)
+              graph[0].emplace_back(cell);
+          }
+      }
     n_colors = graph.size();
 
     helper.setup_color_arrays(n_colors);