Global coarsening: compress weights

diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h
index ae7c1c81e3851a83c6aa80fe09460bed7dfc110c..5948e7c355704903c27b88ed3d039d1587552d1d 100644 (file)
--- a/include/deal.II/matrix_free/tensor_product_kernels.h
+++ b/include/deal.II/matrix_free/tensor_product_kernels.h
@@ -2586,6 +2586,40 @@ namespace internal
   }
 
 
+  template <int dim, int loop_length_template, typename Number>
+  inline void
+  weight_fe_q_dofs_by_entity(const VectorizedArray<Number> *weights,
+                             const unsigned int             n_components,
+                             const int                loop_length_non_template,
+                             VectorizedArray<Number> *data)
+  {
+    const int loop_length = loop_length_template != -1 ?
+                              loop_length_template :
+                              loop_length_non_template;
+
+    Assert(loop_length > 0, ExcNotImplemented());
+    Assert(loop_length < 100, ExcNotImplemented());
+    unsigned int degree_to_3[100];
+    degree_to_3[0] = 0;
+    for (int i = 1; i < loop_length - 1; ++i)
+      degree_to_3[i] = 1;
+    degree_to_3[loop_length - 1] = 2;
+    for (unsigned int c = 0; c < n_components; ++c)
+      for (int k = 0; k < (dim > 2 ? loop_length : 1); ++k)
+        for (int j = 0; j < (dim > 1 ? loop_length : 1); ++j)
+          {
+            const unsigned int shift = 9 * degree_to_3[k] + 3 * degree_to_3[j];
+            data[0] *= weights[shift];
+            // loop bound as int avoids compiler warnings in case loop_length
+            // == 1 (polynomial degree 0)
+            for (int i = 1; i < loop_length - 1; ++i)
+              data[i] *= weights[shift + 1];
+            data[loop_length - 1] *= weights[shift + 2];
+            data += loop_length;
+          }
+  }
+
+
 } // end of namespace internal
 
 

diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.h
index e1a42451cb4c808ac47beaa2a6cae30711546866..35767cde1a52d264293c067d2e6a479d24502dd1 100644 (file)
--- a/include/deal.II/multigrid/mg_transfer_global_coarsening.h
+++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.h
@@ -419,6 +419,13 @@ private:
    */
   std::vector<Number> weights;
 
+  /**
+   * Weights for continuous elements, compressed into 3^dim doubles per
+   * cell if possible.
+   */
+  std::vector<std::array<VectorizedArray<Number>, Utilities::pow(3, dim)>>
+    weights_compressed;
+
   /**
    * DoF indices of the coarse cells, expressed in indices local to the MPI
    * rank.

diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
index ed233633cd56f8423ceb936448f56359a28cb1c3..b3dfc48ecce7ea6fdc4fd14371df84a42a436fad 100644 (file)
--- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
+++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
@@ -40,6 +40,7 @@
 #include <deal.II/hp/dof_handler.h>
 
 #include <deal.II/matrix_free/evaluation_kernels.h>
+#include <deal.II/matrix_free/tensor_product_kernels.h>
 
 #include <deal.II/multigrid/mg_tools.h>
 #include <deal.II/multigrid/mg_transfer_global_coarsening.h>
@@ -1079,6 +1080,121 @@ namespace internal
 
 
 
+    /**
+     * Try to compress weights to Utilities::pow(3, dim) doubles.
+     * Weights of a cell can be compressed if all components have the
+     * same weights.
+     */
+    template <int dim, typename Number>
+    static void
+    compress_weights(
+      MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>
+        &transfer)
+    {
+      const Number *weights = transfer.weights.data();
+
+      // Helper function that extracts the weight. Before doing so,
+      //  it checks the values already set.
+      const auto set = [](auto &mask, const auto &weight) {
+        if (mask == -1.0 || mask == weight)
+          {
+            mask = weight;
+            return true;
+          }
+
+        return false;
+      };
+
+      std::vector<std::array<Number, Utilities::pow(3, dim)>> masks;
+
+      // loop over all schemes
+      for (const auto &scheme : transfer.schemes)
+        {
+          const int    loop_length = scheme.degree_fine + 1;
+          unsigned int degree_to_3[100];
+          degree_to_3[0] = 0;
+          for (int i = 1; i < loop_length - 1; ++i)
+            degree_to_3[i] = 1;
+          degree_to_3[loop_length - 1] = 2;
+
+          // loop over all cells
+          for (unsigned int cell = 0; cell < scheme.n_coarse_cells; ++cell)
+            {
+              std::vector<std::array<Number, Utilities::pow(3, dim)>> mask(
+                transfer.n_components);
+
+              // loop over all components
+              for (unsigned int c = 0; c < transfer.n_components; ++c)
+                mask[c].fill(-1.0);
+
+              for (unsigned int c = 0; c < transfer.n_components; ++c)
+                for (int k = 0; k < (dim > 2 ? loop_length : 1); ++k)
+                  for (int j = 0; j < (dim > 1 ? loop_length : 1); ++j)
+                    {
+                      const unsigned int shift =
+                        9 * degree_to_3[k] + 3 * degree_to_3[j];
+
+                      if (!set(mask[c][shift], weights[0]))
+                        return; // early return, since the entry already
+                                // had a different weight
+
+                      for (int i = 1; i < loop_length - 1; ++i)
+                        if (!set(mask[c][shift + 1], weights[i]))
+                          return;
+
+                      if (!set(mask[c][shift + 2], weights[loop_length - 1]))
+                        return;
+
+                      weights += loop_length;
+                    }
+
+              if (std::find_if(mask.begin(), mask.end(), [&](const auto &a) {
+                    return a != mask[0];
+                  }) != mask.end())
+                return;
+
+              masks.push_back(mask[0]);
+            }
+        }
+
+      // vectorize weights
+      std::vector<std::array<VectorizedArray<Number>, Utilities::pow(3, dim)>>
+        masks_vectorized;
+
+      const unsigned int n_lanes = VectorizedArray<Number>::size();
+
+      unsigned int c = 0;
+
+      for (const auto &scheme : transfer.schemes)
+        {
+          for (unsigned int cell = 0; cell < scheme.n_coarse_cells;
+               cell += n_lanes)
+            {
+              const unsigned int n_lanes_filled =
+                (cell + n_lanes > scheme.n_coarse_cells) ?
+                  (scheme.n_coarse_cells - cell) :
+                  n_lanes;
+
+              std::array<VectorizedArray<Number>, Utilities::pow(3, dim)> mask;
+
+              for (unsigned int v = 0; v < n_lanes_filled; ++v, ++c)
+                {
+                  for (unsigned int i = 0;
+                       i < Utilities::pow<unsigned int>(3, dim);
+                       ++i)
+                    mask[i][v] = masks[c][i];
+                }
+
+              masks_vectorized.push_back(mask);
+            }
+        }
+
+      // copy result
+      transfer.weights_compressed = masks_vectorized;
+    }
+
+
+
   public:
     template <int dim, typename Number>
     static void
@@ -1592,6 +1708,9 @@ namespace internal
                   weights_1 += transfer.schemes[1].n_dofs_per_cell_fine;
                 }
             });
+
+          if (is_feq)
+            compress_weights(transfer);
         }
     }
 
@@ -1669,6 +1788,16 @@ namespace internal
              ExcNotImplemented());
 #endif
 
+      const bool is_feq =
+        std::all_of(dof_handler_fine.get_fe_collection().begin(),
+                    dof_handler_fine.get_fe_collection().end(),
+                    [](const auto &fe) {
+                      return fe.n_base_elements() == 1 &&
+                             (dynamic_cast<const FE_Q<dim> *>(
+                                &fe.base_element(0)) != nullptr);
+                      ;
+                    });
+
       const auto process_cells = [&](const auto &fu) {
         loop_over_active_or_level_cells(
           dof_handler_coarse, mg_level_coarse, [&](const auto &cell_coarse) {
@@ -2096,6 +2225,9 @@ namespace internal
             weights_[fe_pair_no] +=
               transfer.schemes[fe_pair_no].n_dofs_per_cell_fine;
           });
+
+          if (is_feq)
+            compress_weights(transfer);
         }
     }
   };
@@ -2360,9 +2492,14 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
   const unsigned int *indices_coarse = level_dof_indices_coarse.data();
   const unsigned int *indices_fine   = level_dof_indices_fine.data();
   const Number *      weights        = nullptr;
+  const std::array<VectorizedArray<Number>, Utilities::pow(3, dim)>
+    *weights_compressed = nullptr;
 
   if (fine_element_is_continuous)
-    weights = this->weights.data();
+    {
+      weights            = this->weights.data();
+      weights_compressed = this->weights_compressed.data();
+    }
 
   for (const auto &scheme : schemes)
     {
@@ -2370,29 +2507,44 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
       if (scheme.prolongation_matrix.size() == 0 &&
           scheme.prolongation_matrix_1d.size() == 0)
         {
-          for (unsigned int cell = 0; cell < scheme.n_coarse_cells; ++cell)
+          for (unsigned int cell = 0; cell < scheme.n_coarse_cells;
+               cell += n_lanes)
             {
-              if ((scheme.n_dofs_per_cell_fine != 0) &&
-                  (scheme.n_dofs_per_cell_coarse != 0))
+              const unsigned int n_lanes_filled =
+                (cell + n_lanes > scheme.n_coarse_cells) ?
+                  (scheme.n_coarse_cells - cell) :
+                  n_lanes;
+
+              // read from source vector
+              for (unsigned int v = 0; v < n_lanes_filled; ++v)
                 {
+                  if ((scheme.n_dofs_per_cell_fine != 0) &&
+                      (scheme.n_dofs_per_cell_coarse != 0))
+                    {
+                      if (fine_element_is_continuous)
+                        for (unsigned int i = 0;
+                             i < scheme.n_dofs_per_cell_fine;
+                             ++i)
+                          vec_fine_ptr->local_element(indices_fine[i]) +=
+                            read_dof_values(indices_coarse[i], vec_coarse) *
+                            weights[i];
+                      else
+                        for (unsigned int i = 0;
+                             i < scheme.n_dofs_per_cell_fine;
+                             ++i)
+                          vec_fine_ptr->local_element(indices_fine[i]) +=
+                            read_dof_values(indices_coarse[i], vec_coarse);
+                    }
+
+                  indices_fine += scheme.n_dofs_per_cell_fine;
+                  indices_coarse += scheme.n_dofs_per_cell_coarse;
+
                   if (fine_element_is_continuous)
-                    for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine;
-                         ++i)
-                      vec_fine_ptr->local_element(indices_fine[i]) +=
-                        read_dof_values(indices_coarse[i], vec_coarse) *
-                        weights[i];
-                  else
-                    for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine;
-                         ++i)
-                      vec_fine_ptr->local_element(indices_fine[i]) +=
-                        read_dof_values(indices_coarse[i], vec_coarse);
+                    weights += scheme.n_dofs_per_cell_fine;
                 }
 
-              indices_fine += scheme.n_dofs_per_cell_fine;
-              indices_coarse += scheme.n_dofs_per_cell_coarse;
-
               if (fine_element_is_continuous)
-                weights += scheme.n_dofs_per_cell_fine;
+                weights_compressed += 1;
             }
 
           continue;
@@ -2444,9 +2596,20 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
           // ------------------------------ fine -------------------------------
 
           // weight and write into dst vector
+          if (fine_element_is_continuous && this->weights_compressed.size() > 0)
+            {
+              internal::weight_fe_q_dofs_by_entity<dim, -1, Number>(
+                weights_compressed->data(),
+                n_components,
+                scheme.degree_fine + 1,
+                evaluation_data_fine.begin());
+              weights_compressed += 1;
+            }
+
           for (unsigned int v = 0; v < n_lanes_filled; ++v)
             {
-              if (fine_element_is_continuous)
+              if (fine_element_is_continuous &&
+                  this->weights_compressed.size() == 0)
                 for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; ++i)
                   vec_fine_ptr->local_element(indices_fine[i]) +=
                     evaluation_data_fine[i][v] * weights[i];
@@ -2519,9 +2682,14 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
   const unsigned int *indices_coarse = level_dof_indices_coarse.data();
   const unsigned int *indices_fine   = level_dof_indices_fine.data();
   const Number *      weights        = nullptr;
+  const std::array<VectorizedArray<Number>, Utilities::pow(3, dim)>
+    *weights_compressed = nullptr;
 
   if (fine_element_is_continuous)
-    weights = this->weights.data();
+    {
+      weights            = this->weights.data();
+      weights_compressed = this->weights_compressed.data();
+    }
 
   for (const auto &scheme : schemes)
     {
@@ -2529,33 +2697,47 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
       if (scheme.prolongation_matrix.size() == 0 &&
           scheme.prolongation_matrix_1d.size() == 0)
         {
-          for (unsigned int cell = 0; cell < scheme.n_coarse_cells; ++cell)
+          for (unsigned int cell = 0; cell < scheme.n_coarse_cells;
+               cell += n_lanes)
             {
-              if ((scheme.n_dofs_per_cell_fine != 0) &&
-                  (scheme.n_dofs_per_cell_coarse != 0))
+              const unsigned int n_lanes_filled =
+                (cell + n_lanes > scheme.n_coarse_cells) ?
+                  (scheme.n_coarse_cells - cell) :
+                  n_lanes;
+
+              for (unsigned int v = 0; v < n_lanes_filled; ++v)
                 {
+                  if ((scheme.n_dofs_per_cell_fine != 0) &&
+                      (scheme.n_dofs_per_cell_coarse != 0))
+                    {
+                      if (fine_element_is_continuous)
+                        for (unsigned int i = 0;
+                             i < scheme.n_dofs_per_cell_fine;
+                             ++i)
+                          distribute_local_to_global(
+                            indices_coarse[i],
+                            vec_fine_ptr->local_element(indices_fine[i]) *
+                              weights[i],
+                            this->vec_coarse);
+                      else
+                        for (unsigned int i = 0;
+                             i < scheme.n_dofs_per_cell_fine;
+                             ++i)
+                          distribute_local_to_global(
+                            indices_coarse[i],
+                            vec_fine_ptr->local_element(indices_fine[i]),
+                            this->vec_coarse);
+                    }
+
+                  indices_fine += scheme.n_dofs_per_cell_fine;
+                  indices_coarse += scheme.n_dofs_per_cell_coarse;
+
                   if (fine_element_is_continuous)
-                    for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine;
-                         ++i)
-                      distribute_local_to_global(indices_coarse[i],
-                                                 vec_fine_ptr->local_element(
-                                                   indices_fine[i]) *
-                                                   weights[i],
-                                                 this->vec_coarse);
-                  else
-                    for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine;
-                         ++i)
-                      distribute_local_to_global(indices_coarse[i],
-                                                 vec_fine_ptr->local_element(
-                                                   indices_fine[i]),
-                                                 this->vec_coarse);
+                    weights += scheme.n_dofs_per_cell_fine;
                 }
 
-              indices_fine += scheme.n_dofs_per_cell_fine;
-              indices_coarse += scheme.n_dofs_per_cell_coarse;
-
               if (fine_element_is_continuous)
-                weights += scheme.n_dofs_per_cell_fine;
+                weights_compressed += 1;
             }
 
           continue;
@@ -2583,7 +2765,8 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
           // read from source vector and weight
           for (unsigned int v = 0; v < n_lanes_filled; ++v)
             {
-              if (fine_element_is_continuous)
+              if (fine_element_is_continuous &&
+                  this->weights_compressed.size() == 0)
                 for (unsigned int i = 0; i < scheme.n_dofs_per_cell_fine; ++i)
                   evaluation_data_fine[i][v] =
                     vec_fine_ptr->local_element(indices_fine[i]) * weights[i];
@@ -2598,6 +2781,16 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
                 weights += scheme.n_dofs_per_cell_fine;
             }
 
+          if (fine_element_is_continuous && this->weights_compressed.size() > 0)
+            {
+              internal::weight_fe_q_dofs_by_entity<dim, -1, Number>(
+                weights_compressed->data(),
+                n_components,
+                scheme.degree_fine + 1,
+                evaluation_data_fine.begin());
+              weights_compressed += 1;
+            }
+
           // ------------------------------ fine -------------------------------
           for (int c = n_components - 1; c >= 0; --c)
             {

diff --git a/source/multigrid/mg_transfer_matrix_free.cc b/source/multigrid/mg_transfer_matrix_free.cc
index 9762e2eea782d73b8cab87a1f037824805dc48bc..8010b0fa6f028ccf4dd4bc7ff079b80aa2f92c81 100644 (file)
--- a/source/multigrid/mg_transfer_matrix_free.cc
+++ b/source/multigrid/mg_transfer_matrix_free.cc
@@ -29,6 +29,7 @@
 #include <deal.II/lac/la_parallel_vector.h>
 
 #include <deal.II/matrix_free/evaluation_kernels.h>
+#include <deal.II/matrix_free/tensor_product_kernels.h>
 
 #include <deal.II/multigrid/mg_tools.h>
 #include <deal.II/multigrid/mg_transfer_internal.h>
@@ -355,41 +356,6 @@ MGTransferMatrixFree<dim, Number>::restrict_and_add(
 
 
 
-namespace
-{
-  template <int dim, int degree, typename Number>
-  void
-  weight_dofs_on_child(const VectorizedArray<Number> *weights,
-                       const unsigned int             n_components,
-                       const unsigned int             fe_degree,
-                       VectorizedArray<Number> *      data)
-  {
-    Assert(fe_degree > 0, ExcNotImplemented());
-    Assert(fe_degree < 100, ExcNotImplemented());
-    const int loop_length = degree != -1 ? 2 * degree + 1 : 2 * fe_degree + 1;
-    unsigned int degree_to_3[100];
-    degree_to_3[0] = 0;
-    for (int i = 1; i < loop_length - 1; ++i)
-      degree_to_3[i] = 1;
-    degree_to_3[loop_length - 1] = 2;
-    for (unsigned int c = 0; c < n_components; ++c)
-      for (int k = 0; k < (dim > 2 ? loop_length : 1); ++k)
-        for (int j = 0; j < (dim > 1 ? loop_length : 1); ++j)
-          {
-            const unsigned int shift = 9 * degree_to_3[k] + 3 * degree_to_3[j];
-            data[0] *= weights[shift];
-            // loop bound as int avoids compiler warnings in case loop_length
-            // == 1 (polynomial degree 0)
-            for (int i = 1; i < loop_length - 1; ++i)
-              data[i] *= weights[shift + 1];
-            data[loop_length - 1] *= weights[shift + 2];
-            data += loop_length;
-          }
-  }
-} // namespace
-
-
-
 template <int dim, typename Number>
 template <int degree>
 void
@@ -494,10 +460,13 @@ MGTransferMatrixFree<dim, Number>::do_prolongate_add(
                                                      c * n_scalar_cell_dofs,
                                                    fe_degree + 1,
                                                    2 * fe_degree + 1);
-          weight_dofs_on_child<dim, degree, Number>(
+          internal::weight_fe_q_dofs_by_entity<dim,
+                                               degree != -1 ? 2 * degree + 1 :
+                                                              -1,
+                                               Number>(
             &weights_on_refined[to_level - 1][(cell / vec_size) * three_to_dim],
             n_components,
-            fe_degree,
+            2 * fe_degree + 1,
             evaluation_data.begin());
         }
       else
@@ -575,12 +544,14 @@ MGTransferMatrixFree<dim, Number>::do_restrict_add(
       // perform tensorized operation
       if (element_is_continuous)
         {
-          weight_dofs_on_child<dim, degree, Number>(
-            &weights_on_refined[from_level - 1]
-                               [(cell / vec_size) * three_to_dim],
-            n_components,
-            fe_degree,
-            evaluation_data.data());
+          internal::weight_fe_q_dofs_by_entity<
+            dim,
+            degree != -1 ? 2 * degree + 1 : -1,
+            Number>(&weights_on_refined[from_level - 1]
+                                       [(cell / vec_size) * three_to_dim],
+                    n_components,
+                    2 * fe_degree + 1,
+                    evaluation_data.data());
           for (unsigned int c = 0; c < n_components; ++c)
             internal::FEEvaluationImplBasisChange<
               internal::evaluate_general,