]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement better strategy for data locality detection
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 4 May 2020 11:17:48 +0000 (13:17 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Wed, 6 May 2020 17:37:46 +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
tests/matrix_free/thread_correctness.cc

index f26f18b63b9f11332dd0c67d5b547f048e34bdd6..faaac4e085197c1a7fed54853ba31e828ebd96c7 100644 (file)
@@ -463,7 +463,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::internal_reinit(
       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);
+        dummy, 1, false, dummy, false, dummy, dummy, dummy2);
 
       for (unsigned int i = 0; i < dof_info.size(); ++i)
         {
@@ -1058,6 +1058,8 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
           }
       task_info.create_blocks_serial(subdomain_boundary_cells,
                                      dofs_per_cell,
+                                     dof_handlers.active_dof_handler ==
+                                       DoFHandlers::hp,
                                      dof_info[0].cell_active_fe_index,
                                      strict_categories,
                                      parent_relation,
index ed9a128b0270c27b40279d6cef12f2057a74538d..61299b88a394b5fc9a5ab0fa16cd91d85883f967 100644 (file)
@@ -145,14 +145,20 @@ namespace internal
        * Sets up the blocks for running the cell loop based on the options
        * controlled by the input arguments.
        *
-       * @param boundary_cells A list of cells that need to exchange data prior
-       * to performing computations. These will be given a certain id in the
-       * partitioning.
+       * @param cells_with_comm A list of cells that need to exchange data
+       * prior to performing computations. These will be given a certain id in
+       * the partitioning to make sure cell loops that overlap communication
+       * with communication have the ghost data ready.
        *
        * @param dofs_per_cell Gives an expected value for the number of degrees
        * of freedom on a cell, which is used to determine the block size for
        * interleaving cell and face integrals.
        *
+       * @param categories_are_hp Defines whether
+       * `cell_vectorization_categories` is originating from a hp adaptive
+       * computation with variable polynomial degree or a user-defined
+       * variant.
+       *
        * @param cell_vectorization_categories This set of categories defines
        * the cells that should be grouped together inside the lanes of a
        * vectorized array. This can be the polynomial degree in an hp-element
@@ -179,8 +185,9 @@ namespace internal
        */
       void
       create_blocks_serial(
-        const std::vector<unsigned int> &boundary_cells,
+        const std::vector<unsigned int> &cells_with_comm,
         const unsigned int               dofs_per_cell,
+        const bool                       categories_are_hp,
         const std::vector<unsigned int> &cell_vectorization_categories,
         const bool                       cell_vectorization_categories_strict,
         const std::vector<unsigned int> &parent_relation,
index 877395765d5faa13c929b3064b5af6dff68e9db2..953041fc536e3365196092e2595e48b4e85e4f21 100644 (file)
@@ -768,175 +768,42 @@ namespace internal
 
     void
     TaskInfo::create_blocks_serial(
-      const std::vector<unsigned int> &boundary_cells,
+      const std::vector<unsigned int> &cells_with_comm,
       const unsigned int               dofs_per_cell,
+      const bool                       categories_are_hp,
       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)
     {
-      const unsigned int n_macro_cells =
-        (n_active_cells + vectorization_length - 1) / vectorization_length;
-      const unsigned int n_ghost_slots =
-        (n_ghost_cells + vectorization_length - 1) / vectorization_length;
-
-      incompletely_filled_vectorization.resize(n_macro_cells + n_ghost_slots);
-      renumbering.resize(n_active_cells + n_ghost_cells,
-                         numbers::invalid_unsigned_int);
-
-      // Define the outer number of partitions. In the MPI case, we have three
-      // partitions (part before comm, part with comm, part after comm)
-      if (n_procs == 1)
-        partition_row_index.resize(3);
-      else
-        partition_row_index.resize(5);
-
-      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)
-        {
-          // 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(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 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;
-          for (; c < n_active_cells + n_ghost_cells; ++c)
-            if (cell_marked[c] == 0)
-              cell_marked[c] = 4;
-        }
-      else
-        std::fill(cell_marked.begin(), cell_marked.end(), 1);
-
-      for (const unsigned char marker : cell_marked)
-        {
-          (void)marker;
-          Assert(marker != 0, ExcInternalError());
-        }
-
+      // Give the compiler a chance to detect that vectorization_length is a
+      // power of two, which allows it to replace integer divisions by shifts
+      unsigned int vectorization_length_bits = 0;
+      unsigned int my_length                 = vectorization_length;
+      while (my_length >>= 1)
+        ++vectorization_length_bits;
+      const unsigned int n_lanes = 1 << vectorization_length_bits;
+
+      // Step 1: find tight map of categories for not taking exceeding amounts
+      // of memory below. Sort the new categories by the numbers in the
+      // old one to ensure we respect the given rules
       unsigned int              n_categories = 1;
-      std::vector<unsigned int> tight_category_map;
+      std::vector<unsigned int> tight_category_map(n_active_cells, 0);
       if (cell_vectorization_categories.empty() == false)
         {
           AssertDimension(cell_vectorization_categories.size(),
                           n_active_cells + n_ghost_cells);
 
-          // create a tight map of categories for not taking exceeding amounts
-          // of memory below. Sort the new categories by the numbers in the
-          // old one.
-          tight_category_map.resize(n_active_cells + n_ghost_cells);
           std::set<unsigned int> used_categories;
-          for (unsigned int i = 0; i < n_active_cells + n_ghost_cells; ++i)
+          for (unsigned int i = 0; i < n_active_cells; ++i)
             used_categories.insert(cell_vectorization_categories[i]);
           std::vector<unsigned int> used_categories_vector(
             used_categories.size());
           n_categories = 0;
           for (const auto &it : used_categories)
             used_categories_vector[n_categories++] = it;
-          for (unsigned int i = 0; i < n_active_cells + n_ghost_cells; ++i)
+          for (unsigned int i = 0; i < n_active_cells; ++i)
             {
               const unsigned int index =
                 std::lower_bound(used_categories_vector.begin(),
@@ -946,163 +813,222 @@ namespace internal
               AssertIndexRange(index, used_categories_vector.size());
               tight_category_map[i] = index;
             }
-
-          // leave some more space for empty lanes
-          incompletely_filled_vectorization.resize(
-            incompletely_filled_vectorization.size() + 4 * n_categories);
         }
-      else
-        tight_category_map.resize(n_active_cells + n_ghost_cells, 0);
 
-      cell_partition_data.clear();
-      cell_partition_data.resize(1, 0);
-      unsigned int                           counter = 0;
-      unsigned int                           n_cells = 0;
+      // Step 2: Sort the cells by the category. If we want to fill up the
+      // ranges in vectorization, promote some of the cells to a higher
+      // category
       std::vector<std::vector<unsigned int>> renumbering_category(n_categories);
-      for (unsigned int block = 1; block < (n_procs > 1u ? 5u : 3u); ++block)
-        {
-          // step 1: sort by category
-          for (unsigned int i = 0; i < n_active_cells + n_ghost_cells; ++i)
-            if (cell_marked[i] == block)
-              renumbering_category[tight_category_map[i]].push_back(i);
-
-          // step 2: if we want to fill up the ranges in vectorization, promote
-          // some of the cells to a higher category
-          if (cell_vectorization_categories_strict == false && n_categories > 1)
-            for (unsigned int j = n_categories - 1; j > 0; --j)
+      for (unsigned int i = 0; i < n_active_cells; ++i)
+        renumbering_category[tight_category_map[i]].push_back(i);
+
+      if (cell_vectorization_categories_strict == false && n_categories > 1)
+        for (unsigned int j = n_categories - 1; j > 0; --j)
+          {
+            unsigned int lower_index = j - 1;
+            while (renumbering_category[j].size() % n_lanes)
               {
-                unsigned int lower_index = j - 1;
-                while (renumbering_category[j].size() % vectorization_length)
+                while (renumbering_category[j].size() % n_lanes &&
+                       !renumbering_category[lower_index].empty())
                   {
-                    while (renumbering_category[j].size() %
-                             vectorization_length &&
-                           !renumbering_category[lower_index].empty())
-                      {
-                        renumbering_category[j].push_back(
-                          renumbering_category[lower_index].back());
-                        renumbering_category[lower_index].pop_back();
-                      }
-                    if (lower_index == 0)
-                      break;
-                    else
-                      --lower_index;
+                    renumbering_category[j].push_back(
+                      renumbering_category[lower_index].back());
+                    renumbering_category[lower_index].pop_back();
                   }
+                if (lower_index == 0)
+                  break;
+                else
+                  --lower_index;
               }
+          }
 
-          // step 3: append cells according to categories
-          for (unsigned int j = 0; j < n_categories; ++j)
-            {
+      // Step 3: Use the parent relation to find a good grouping of cells. To
+      // do this, we first put cells of each category defined above into two
+      // bins, those which we know can be grouped together by the given parent
+      // relation and those which cannot
+      std::vector<unsigned int> temporary_numbering;
+      temporary_numbering.reserve(n_active_cells +
+                                  (n_lanes - 1) * n_categories);
+      const unsigned int n_cells_per_parent =
+        std::count(parent_relation.begin(), parent_relation.end(), 0);
+      std::vector<unsigned int> category_size;
+      for (unsigned int j = 0; j < n_categories; ++j)
+        {
+          std::vector<std::pair<unsigned int, unsigned int>> grouped_cells;
+          std::vector<unsigned int>                          other_cells;
+          for (const unsigned int cell : renumbering_category[j])
+            if (parent_relation.empty() ||
+                parent_relation[cell] == numbers::invalid_unsigned_int)
+              other_cells.push_back(cell);
+            else
+              grouped_cells.emplace_back(parent_relation[cell], cell);
+
+          // Compute the number of cells per group
+          std::sort(grouped_cells.begin(), grouped_cells.end());
+          std::vector<unsigned int> n_cells_per_group;
+          unsigned int              length = 0;
+          for (unsigned int i = 0; i < grouped_cells.size(); ++i, ++length)
+            if (i > 0 && grouped_cells[i].first != grouped_cells[i - 1].first)
               {
-                // 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);
+                n_cells_per_group.push_back(length);
+                length = 0;
+              }
+          if (length > 0)
+            n_cells_per_group.push_back(length);
+
+          // Move groups that do not have the complete size (due to
+          // categories) to the 'other_cells'. The cells with correct group
+          // size are immediately appended to the temporary cell numbering
+          auto group_it = grouped_cells.begin();
+          for (unsigned int length : n_cells_per_group)
+            if (length < n_cells_per_parent)
+              for (unsigned int j = 0; j < length; ++j)
+                other_cells.push_back((group_it++)->second);
+            else
+              {
+                // we should not have more cells in a group than in the first
+                // check we did above
+                AssertDimension(length, n_cells_per_parent);
+                for (unsigned int j = 0; j < length; ++j)
+                  temporary_numbering.push_back((group_it++)->second);
+              }
 
-                // Sort the remaining cells
-                std::sort(other_cells.begin(), other_cells.end());
+          // Sort the remaining cells and append them as well
+          std::sort(other_cells.begin(), other_cells.end());
+          temporary_numbering.insert(temporary_numbering.end(),
+                                     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++;
-                  }
-              }
+          while (temporary_numbering.size() % n_lanes != 0)
+            temporary_numbering.push_back(numbers::invalid_unsigned_int);
 
-              unsigned int remainder =
-                renumbering_category[j].size() % vectorization_length;
-              if (remainder)
-                incompletely_filled_vectorization
-                  [renumbering_category[j].size() / vectorization_length +
-                   n_cells] = remainder;
-              const unsigned int n_my_macro_cells =
-                (renumbering_category[j].size() + vectorization_length - 1) /
-                vectorization_length;
-              renumbering_category[j].clear();
-
-              // step 4: create blocks for face integrals, make the number of
-              // cells divisible by 4 if possible
-              const unsigned int block_size =
-                std::max((2048U / dofs_per_cell) / 8 * 4, 2U);
-              if (block < 4)
-                for (unsigned int k = 0; k < n_my_macro_cells; k += block_size)
-                  cell_partition_data.push_back(
-                    n_cells + std::min(k + block_size, n_my_macro_cells));
-              else
-                cell_partition_data.back() += n_my_macro_cells;
-              n_cells += n_my_macro_cells;
-            }
-          partition_row_index[block] = cell_partition_data.size() - 1;
-          if (block == 3 || (block == 1 && n_procs == 1))
-            cell_partition_data.push_back(n_cells);
+          category_size.push_back(temporary_numbering.size());
+        }
+
+      // Step 4: Identify the batches with cells marked as "comm"
+      std::vector<bool> batch_with_comm(temporary_numbering.size() / n_lanes,
+                                        false);
+      std::vector<unsigned int> temporary_numbering_inverse(n_active_cells);
+      for (unsigned int i = 0; i < temporary_numbering.size(); ++i)
+        if (temporary_numbering[i] != numbers::invalid_unsigned_int)
+          temporary_numbering_inverse[temporary_numbering[i]] = i;
+      for (const unsigned int cell : cells_with_comm)
+        batch_with_comm[temporary_numbering_inverse[cell] / n_lanes] = true;
+
+      // Step 5: Sort the batches of cells by their last cell index to get
+      // good locality, assuming that the initial cell order is of good
+      // locality. In case we have hp calculations with categories, we need to
+      // sort also by the category.
+      std::vector<std::array<unsigned int, 3>> batch_order;
+      std::vector<std::array<unsigned int, 3>> batch_order_comm;
+      for (unsigned int i = 0; i < temporary_numbering.size(); i += n_lanes)
+        {
+          unsigned int max_index = 0;
+          for (unsigned int j = 0; j < n_lanes; ++j)
+            if (temporary_numbering[i + j] < numbers::invalid_unsigned_int)
+              max_index = std::max(temporary_numbering[i + j], max_index);
+          const unsigned int category_hp =
+            categories_are_hp ?
+              std::upper_bound(category_size.begin(), category_size.end(), i) -
+                category_size.begin() :
+              0;
+          const std::array<unsigned int, 3> next{category_hp, max_index, i};
+          if (batch_with_comm[i / n_lanes])
+            batch_order_comm.emplace_back(next);
+          else
+            batch_order.emplace_back(next);
         }
-      if (cell_vectorization_categories_strict == true)
+
+      std::sort(batch_order.begin(), batch_order.end());
+      std::sort(batch_order_comm.begin(), batch_order_comm.end());
+
+      // Step 6: Put the cells with communication in the middle of the cell
+      // range. For the MPI case, we need three groups to enable overlap for
+      // communication and computation (part before comm, part with comm, part
+      // after comm), whereas we need one for the other case. And in each
+      // case, we allow for a slot of "ghosted" cells.
+      std::vector<unsigned int> blocks;
+      if (n_procs == 1)
         {
-          Assert(n_cells >= n_macro_cells + n_ghost_slots, ExcInternalError());
+          if (batch_order.empty())
+            std::swap(batch_order_comm, batch_order);
+          Assert(batch_order_comm.empty(), ExcInternalError());
+          partition_row_index.resize(3);
+          blocks = {0, static_cast<unsigned int>(batch_order.size())};
         }
       else
         {
-          AssertDimension(n_cells, n_macro_cells + n_ghost_slots);
+          partition_row_index.resize(5);
+          const unsigned int comm_begin = batch_order.size() / 2;
+          batch_order.insert(batch_order.begin() + comm_begin,
+                             batch_order_comm.begin(),
+                             batch_order_comm.end());
+          const unsigned int comm_end = comm_begin + batch_order_comm.size();
+          const unsigned int end      = batch_order.size();
+          blocks                      = {0, comm_begin, comm_end, end};
         }
-      AssertDimension(cell_partition_data.back(), n_cells);
-      AssertDimension(counter, n_active_cells + n_ghost_cells);
 
-      incompletely_filled_vectorization.resize(cell_partition_data.back());
+      // Step 7: Fill in the data by batches for the locally owned cells.
+      const unsigned int n_cell_batches = batch_order.size();
+      const unsigned int n_ghost_batches =
+        (n_ghost_cells + n_lanes - 1) / n_lanes;
+      incompletely_filled_vectorization.resize(n_cell_batches +
+                                               n_ghost_batches);
+
+      cell_partition_data.clear();
+      cell_partition_data.resize(1, 0);
+
+      renumbering.clear();
+      renumbering.resize(n_active_cells + n_ghost_cells,
+                         numbers::invalid_unsigned_int);
+
+      unsigned int counter = 0;
+      for (unsigned int block = 0; block < blocks.size() - 1; ++block)
+        {
+          const unsigned int grain_size =
+            std::max((2048U / dofs_per_cell) / 8 * 4, 2U);
+          for (unsigned int k = blocks[block]; k < blocks[block + 1];
+               k += grain_size)
+            cell_partition_data.push_back(
+              std::min(k + grain_size, blocks[block + 1]));
+          partition_row_index[block + 1] = cell_partition_data.size() - 1;
+
+          // Set the numbering according to the reordered temporary one
+          for (unsigned int k = blocks[block]; k < blocks[block + 1]; ++k)
+            {
+              const unsigned int pos = batch_order[k][2];
+              unsigned int       j   = 0;
+              for (; j < n_lanes && temporary_numbering[pos + j] !=
+                                      numbers::invalid_unsigned_int;
+                   ++j)
+                renumbering[counter++] = temporary_numbering[pos + j];
+              if (j < n_lanes)
+                incompletely_filled_vectorization[k] = j;
+            }
+        }
+      AssertDimension(counter, n_active_cells);
+
+      // Step 8: Treat the ghost cells
+      for (unsigned int cell = n_active_cells;
+           cell < n_active_cells + n_ghost_cells;
+           ++cell)
+        {
+          if (!cell_vectorization_categories.empty())
+            AssertDimension(cell_vectorization_categories[cell],
+                            cell_vectorization_categories[n_active_cells]);
+          renumbering[cell] = cell;
+        }
+      if (n_ghost_cells % n_lanes)
+        incompletely_filled_vectorization.back() = n_ghost_cells % n_lanes;
+      cell_partition_data.push_back(n_cell_batches + n_ghost_batches);
+      partition_row_index.back() = cell_partition_data.size() - 1;
+
+#ifdef DEBUG
+      std::vector<unsigned int> renumber_cpy(renumbering);
+      std::sort(renumber_cpy.begin(), renumber_cpy.end());
+      for (unsigned int i = 0; i < renumber_cpy.size(); ++i)
+        AssertDimension(i, renumber_cpy[i]);
+#endif
     }
 
 
index 4d9f6681d65359e69eedd6026b194dff6e1c52c2..d8519f72e934287f0d672ec3203e86e1fa21bfee 100644 (file)
@@ -128,6 +128,8 @@ sub_test()
 
       // make 10 sweeps in order to get in some
       // variation to the threaded program
+      const double float_factor =
+        std::is_same<number, float>::value ? 0.01 : 1.;
       for (unsigned int sweep = 0; sweep < 10; ++sweep)
         {
           mf_color.vmult(out_color, in_dist);
@@ -135,14 +137,12 @@ sub_test()
 
           out_color -= out_dist;
           double diff_norm = out_color.linfty_norm();
-          deallog << "Sweep " << sweep
-                  << ", error in partition/color:     " << diff_norm
-                  << std::endl;
+          deallog << "Sweep " << sweep << ", error in partition/color:     "
+                  << diff_norm * float_factor << std::endl;
           out_partition -= out_dist;
           diff_norm = out_partition.linfty_norm();
-          deallog << "Sweep " << sweep
-                  << ", error in partition/partition: " << diff_norm
-                  << std::endl;
+          deallog << "Sweep " << sweep << ", error in partition/partition: "
+                  << diff_norm * float_factor << std::endl;
         }
       deallog << std::endl;
     }

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.