]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Group cells by parents in matrix-free setup
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 2 May 2020 13:52:12 +0000 (15:52 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 6 May 2020 17:37:45 +0000 (19:37 +0200)
include/deal.II/matrix_free/matrix_free.templates.h
include/deal.II/matrix_free/task_info.h
source/matrix_free/task_info.cc

index ca57050fc0cd8faad3d2c342d4191c86ac6bea39..f26f18b63b9f11332dd0c67d5b547f048e34bdd6 100644 (file)
@@ -460,11 +460,11 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
       initialize_dof_handlers(dof_handler, additional_data);
       std::vector<unsigned int>  dummy;
       std::vector<unsigned char> dummy2;
-      task_info.collect_boundary_cells(cell_level_index.size(),
-                                       cell_level_index.size(),
-                                       VectorizedArrayType::size(),
-                                       dummy);
-      task_info.create_blocks_serial(dummy, 1, dummy, false, dummy, dummy2);
+      task_info.vectorization_length = VectorizedArrayType::size();
+      task_info.n_active_cells       = cell_level_index.size();
+      task_info.create_blocks_serial(
+        dummy, 1, dummy, false, dummy, dummy, dummy2);
+
       for (unsigned int i = 0; i < dof_info.size(); ++i)
         {
           Assert(dof_handler[i]->get_fe_collection().size() == 1,
@@ -970,11 +970,10 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
         subdomain_boundary_cells.push_back(counter);
     }
 
-  const unsigned int n_lanes = VectorizedArrayType::size();
-  task_info.collect_boundary_cells(cell_level_index_end_local,
-                                   n_active_cells,
-                                   n_lanes,
-                                   subdomain_boundary_cells);
+  const unsigned int n_lanes     = VectorizedArrayType::size();
+  task_info.n_active_cells       = cell_level_index_end_local;
+  task_info.n_ghost_cells        = n_active_cells - cell_level_index_end_local;
+  task_info.vectorization_length = n_lanes;
 
   // Finalize the creation of the ghost indices
   {
@@ -1022,21 +1021,53 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
   std::vector<unsigned char> irregular_cells;
   if (task_info.scheme == internal::MatrixFreeFunctions::TaskInfo::none)
     {
-      const bool strict_categories =
+      bool strict_categories =
         additional_data.cell_vectorization_categories_strict ||
         dof_handlers.active_dof_handler == DoFHandlers::hp;
       unsigned int dofs_per_cell = 0;
       for (const auto &info : dof_info)
         dofs_per_cell = std::max(dofs_per_cell, info.dofs_per_cell[0]);
+
+      // Detect cells with the same parent to make sure they get scheduled
+      // together in the loop, which increases data locality.
+      std::vector<unsigned int> parent_relation(task_info.n_active_cells +
+                                                  task_info.n_ghost_cells,
+                                                numbers::invalid_unsigned_int);
+      std::map<std::pair<int, int>, std::vector<unsigned int>> cell_parents;
+      for (unsigned int c = 0; c < cell_level_index_end_local; ++c)
+        if (cell_level_index[c].first > 0)
+          {
+            typename Triangulation<dim>::cell_iterator cell(
+              dof_handlers.active_dof_handler == DoFHandlers::usual ?
+                &dof_handlers.dof_handler[0]->get_triangulation() :
+                &dof_handlers.hp_dof_handler[0]->get_triangulation(),
+              cell_level_index[c].first,
+              cell_level_index[c].second);
+            Assert(cell->level() > 0, ExcInternalError());
+            cell_parents[std::make_pair(cell->parent()->level(),
+                                        cell->parent()->index())]
+              .push_back(c);
+          }
+      unsigned int position = 0;
+      for (const auto &it : cell_parents)
+        if (it.second.size() == GeometryInfo<dim>::max_children_per_cell)
+          {
+            for (auto i : it.second)
+              parent_relation[i] = position;
+            ++position;
+          }
       task_info.create_blocks_serial(subdomain_boundary_cells,
                                      dofs_per_cell,
                                      dof_info[0].cell_active_fe_index,
                                      strict_categories,
+                                     parent_relation,
                                      renumbering,
                                      irregular_cells);
     }
   else
     {
+      task_info.make_boundary_cells_divisible(subdomain_boundary_cells);
+
       // For strategy with blocking before partitioning: reorganize the indices
       // in order to overlap communication in MPI with computations: Place all
       // cells with ghost indices into one chunk. Also reorder cells so that we
index cfc972ccd99d63338e504308dddad3ffed7340c0..ed9a128b0270c27b40279d6cef12f2057a74538d 100644 (file)
@@ -135,14 +135,11 @@ namespace internal
       loop(MFWorkerInterface &worker) const;
 
       /**
-       * Determines the position of cells with ghosts for distributed-memory
-       * calculations.
+       * Make the number of cells which can only be treated in the
+       * communication overlap divisible by the vectorization length.
        */
       void
-      collect_boundary_cells(const unsigned int n_active_cells,
-                             const unsigned int n_active_and_ghost_cells,
-                             const unsigned int vectorization_length,
-                             std::vector<unsigned int> &boundary_cells);
+      make_boundary_cells_divisible(std::vector<unsigned int> &boundary_cells);
 
       /**
        * Sets up the blocks for running the cell loop based on the options
@@ -166,6 +163,10 @@ namespace internal
        * strictly or whether it is allowed to insert lower categories into the
        * next high one(s).
        *
+       * @param parent_relation This data field is used to specify which cells
+       * have the same parent cell. Cells with the same ancestor are grouped
+       * together into the same batch(es) with vectorization across cells.
+       *
        * @param renumbering When leaving this function, the vector contains a
        * new numbering of the cells that aligns with the grouping stored in
        * this class.
@@ -182,6 +183,7 @@ namespace internal
         const unsigned int               dofs_per_cell,
         const std::vector<unsigned int> &cell_vectorization_categories,
         const bool                       cell_vectorization_categories_strict,
+        const std::vector<unsigned int> &parent_relation,
         std::vector<unsigned int> &      renumbering,
         std::vector<unsigned char> &     incompletely_filled_vectorization);
 
index cbae4ea508e60c3fa95312de07e14c87c8efc42f..877395765d5faa13c929b3064b5af6dff68e9db2 100644 (file)
@@ -695,16 +695,9 @@ namespace internal
 
 
     void
-    TaskInfo::collect_boundary_cells(
-      const unsigned int         n_active_cells_in,
-      const unsigned int         n_active_and_ghost_cells,
-      const unsigned int         vectorization_length_in,
+    TaskInfo::make_boundary_cells_divisible(
       std::vector<unsigned int> &boundary_cells)
     {
-      vectorization_length = vectorization_length_in;
-      n_active_cells       = n_active_cells_in;
-      n_ghost_cells        = n_active_and_ghost_cells - n_active_cells;
-
       // try to make the number of boundary cells divisible by the number of
       // vectors in vectorization
       unsigned int fillup_needed =
@@ -779,6 +772,7 @@ namespace internal
       const unsigned int               dofs_per_cell,
       const std::vector<unsigned int> &cell_vectorization_categories,
       const bool                       cell_vectorization_categories_strict,
+      const std::vector<unsigned int> &parent_relation,
       std::vector<unsigned int> &      renumbering,
       std::vector<unsigned char> &     incompletely_filled_vectorization)
     {
@@ -786,7 +780,6 @@ namespace internal
         (n_active_cells + vectorization_length - 1) / vectorization_length;
       const unsigned int n_ghost_slots =
         (n_ghost_cells + vectorization_length - 1) / vectorization_length;
-      const unsigned int n_boundary_cells = boundary_cells.size();
 
       incompletely_filled_vectorization.resize(n_macro_cells + n_ghost_slots);
       renumbering.resize(n_active_cells + n_ghost_cells,
@@ -799,28 +792,115 @@ namespace internal
       else
         partition_row_index.resize(5);
 
-      // Initially mark the cells according to the MPI ranking
+      int max_parent_index = -1;
+      for (unsigned int i : parent_relation)
+        if (i != numbers::invalid_unsigned_int)
+          max_parent_index = std::max(static_cast<int>(i), max_parent_index);
+      unsigned int expected_group_size =
+        max_parent_index != -1 ?
+          std::count(parent_relation.begin(), parent_relation.end(), 0) :
+          1;
+
       std::vector<unsigned char> cell_marked(n_active_cells + n_ghost_cells, 0);
       if (n_procs > 1)
         {
-          for (unsigned int i = 0; i < n_boundary_cells; ++i)
-            cell_marked[boundary_cells[i]] = 2;
+          // This lambda is used to mark the siblings (belong to the same
+          // parent) if a particular cell was touched in a pass as well
+          const auto mark_siblings =
+            [&](const unsigned int       mark,
+                const std::vector<bool> &relevant_parents) {
+              for (unsigned int i = 0; i < n_active_cells; ++i)
+                if (cell_marked[i] == 0 &&
+                    parent_relation[i] != numbers::invalid_unsigned_int &&
+                    relevant_parents[parent_relation[i]])
+                  cell_marked[i] = mark;
+            };
+
+          // This lambda makes the cells at the processor boundary divisible
+          // by the vectorization length; we start the fillup with the more
+          // unstructured cells without a parent to increase chances that the
+          // cells sharing the parent get placed together
+          const auto fill_up_vectorization = [&](const unsigned int mark) {
+            unsigned int n_marked_cells =
+              std::count(cell_marked.begin(),
+                         cell_marked.begin() + n_active_cells,
+                         mark);
+            unsigned int n_available_cells =
+              std::count(cell_marked.begin(),
+                         cell_marked.begin() + n_active_cells,
+                         0);
+            if (n_marked_cells % vectorization_length > 0 &&
+                n_available_cells > 0)
+              {
+                unsigned int n_missing =
+                  vectorization_length -
+                  (n_marked_cells % vectorization_length);
+                for (unsigned int i = 0; i < n_active_cells; ++i)
+                  if (cell_marked[i] == 0 &&
+                      parent_relation[i] == numbers::invalid_unsigned_int)
+                    {
+                      cell_marked[i] = mark;
+                      --n_missing;
+                      --n_available_cells;
+                      ++n_marked_cells;
+                      if (n_missing == 0)
+                        break;
+                    }
+                for (unsigned int i = 0; i < n_active_cells; ++i)
+                  if (cell_marked[n_active_cells - 1 - i] == 0)
+                    {
+                      cell_marked[n_active_cells - 1 - i] = mark;
+                      --n_missing;
+                      --n_available_cells;
+                      ++n_marked_cells;
+                      if (n_missing == 0)
+                        break;
+                    }
+              }
 
-          Assert(boundary_cells.size() % vectorization_length == 0 ||
-                   boundary_cells.size() == n_active_cells,
-                 ExcInternalError());
+            Assert(n_marked_cells % vectorization_length == 0 ||
+                     n_available_cells == 0,
+                   ExcInternalError("error " + std::to_string(n_marked_cells) +
+                                    " " + std::to_string(n_available_cells)));
+            return n_marked_cells;
+          };
 
+          // Mark all cells needing data exchange as well as those belonging
+          // to the same parent with a special number
+          {
+            std::vector<bool> parent_at_boundary(max_parent_index + 1);
+            for (const unsigned int cell : boundary_cells)
+              {
+                cell_marked[cell] = 2;
+                if (parent_relation[cell] != numbers::invalid_unsigned_int)
+                  parent_at_boundary[parent_relation[cell]] = true;
+              }
+            mark_siblings(2, parent_at_boundary);
+          }
+          const unsigned int n_boundary_cells = fill_up_vectorization(2);
+
+          // Mark the cells that get placed before the cells at processor
+          // boundaries
           const unsigned int n_second_slot =
             ((n_active_cells - n_boundary_cells) / 2 / vectorization_length) *
             vectorization_length;
-          unsigned int count = 0;
-          unsigned int c     = 0;
-          for (; c < n_active_cells && count < n_second_slot; ++c)
-            if (cell_marked[c] == 0)
-              {
-                cell_marked[c] = 1;
-                ++count;
-              }
+          unsigned int c = 0;
+          {
+            unsigned int      count = 0;
+            std::vector<bool> parent_marked(max_parent_index + 1, false);
+            for (; c < n_active_cells && count < n_second_slot; ++c)
+              if (cell_marked[c] == 0)
+                {
+                  if (parent_relation[c] != numbers::invalid_unsigned_int)
+                    parent_marked[parent_relation[c]] = true;
+                  cell_marked[c] = 1;
+                  ++count;
+                }
+            mark_siblings(1, parent_marked);
+            fill_up_vectorization(1);
+          }
+
+          // Finally, mark the remaining cells
           for (; c < n_active_cells; ++c)
             if (cell_marked[c] == 0)
               cell_marked[c] = 3;
@@ -912,8 +992,78 @@ namespace internal
           // step 3: append cells according to categories
           for (unsigned int j = 0; j < n_categories; ++j)
             {
-              for (const unsigned int cell : renumbering_category[j])
-                renumbering[counter++] = cell;
+              {
+                // Among the current category, we need to distinguish cells
+                // which we want to have grouped together and other cells. We
+                // first start by setting up these two categories
+                std::vector<std::pair<unsigned int, unsigned int>>
+                                          grouped_cells_tmp;
+                std::vector<unsigned int> other_cells;
+                for (const unsigned int cell : renumbering_category[j])
+                  if (parent_relation[cell] == numbers::invalid_unsigned_int)
+                    other_cells.push_back(cell);
+                  else
+                    grouped_cells_tmp.emplace_back(parent_relation[cell], cell);
+
+                // Create a CRS data structure to identify each of the chunks
+                std::sort(grouped_cells_tmp.begin(), grouped_cells_tmp.end());
+                std::vector<unsigned int> crs_group(1);
+                for (unsigned int i = 1; i < grouped_cells_tmp.size(); ++i)
+                  if (grouped_cells_tmp[i].first !=
+                      grouped_cells_tmp[i - 1].first)
+                    crs_group.push_back(i);
+                crs_group.push_back(grouped_cells_tmp.size());
+
+                // Move groups that do not have the complete size (due to
+                // categories) to the 'other_cells'
+                std::vector<unsigned int> grouped_cells;
+                for (unsigned int i = 0; i < crs_group.size() - 1; ++i)
+                  if (crs_group[i + 1] - crs_group[i] < expected_group_size)
+                    for (unsigned int j = crs_group[i]; j < crs_group[i + 1];
+                         ++j)
+                      other_cells.push_back(grouped_cells_tmp[j].second);
+                  else
+                    for (unsigned int j = crs_group[i]; j < crs_group[i + 1];
+                         ++j)
+                      grouped_cells.push_back(grouped_cells_tmp[j].second);
+
+                // Sort the remaining cells
+                std::sort(other_cells.begin(), other_cells.end());
+
+                // Now fill in the cells from the two slots, the one with
+                // groups and the one without
+                auto regular = grouped_cells.begin();
+                auto fillup  = other_cells.begin();
+                while (regular != grouped_cells.end() ||
+                       fillup != other_cells.end())
+                  {
+                    // Case 1: Fill up until the next expected group size
+                    while (counter % expected_group_size &&
+                           fillup != other_cells.end())
+                      renumbering[counter++] = *fillup++;
+
+                    // Case 2: If the start of the next group has a larger
+                    // index than all indices we have queued from the
+                    // irregular
+                    if (fillup + expected_group_size <= other_cells.end() &&
+                        (regular == grouped_cells.end() ||
+                         *(fillup + expected_group_size - 1) < *regular))
+                      for (unsigned int j = 0; j < expected_group_size; ++j)
+                        renumbering[counter++] = *fillup++;
+
+                    // Case 3: Add a group at once
+                    if (regular != grouped_cells.end())
+                      for (unsigned int j = 0; j < expected_group_size; ++j)
+                        renumbering[counter++] = *regular++;
+
+                    // Case 4: The groups are empty, so fill up from the other
+                    // chunk
+                    else
+                      while (fillup != other_cells.end())
+                        renumbering[counter++] = *fillup++;
+                  }
+              }
+
               unsigned int remainder =
                 renumbering_category[j].size() % vectorization_length;
               if (remainder)

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.