From 157a9c1d1f31a515e9be459b62686b69ee81dbbe Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <martin.kronbichler@rub.de>
Date: Wed, 6 Mar 2024 10:19:03 +0100
Subject: [PATCH] Matrix-free MG transfer: Switch interpolation matrices to
 double type

---
 .../multigrid/mg_transfer_global_coarsening.h | 11 ++--
 .../mg_transfer_global_coarsening.templates.h | 56 ++++++++++---------
 2 files changed, 36 insertions(+), 31 deletions(-)

diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.h
index 8820af597b..6c369fb284 100644
--- a/include/deal.II/multigrid/mg_transfer_global_coarsening.h
+++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.h
@@ -651,29 +651,28 @@ private:
     /**
      * Prolongation matrix for non-tensor-product elements.
      */
-    AlignedVector<VectorizedArrayType> prolongation_matrix;
+    AlignedVector<double> prolongation_matrix;
 
     /**
      * 1d prolongation matrix for tensor-product elements.
      */
-    AlignedVector<VectorizedArrayType> prolongation_matrix_1d;
+    AlignedVector<double> prolongation_matrix_1d;
 
     /**
      * Restriction matrix for non-tensor-product elements.
      */
-    AlignedVector<VectorizedArrayType> restriction_matrix;
+    AlignedVector<double> restriction_matrix;
 
     /**
      * 1d restriction matrix for tensor-product elements.
      */
-    AlignedVector<VectorizedArrayType> restriction_matrix_1d;
+    AlignedVector<double> restriction_matrix_1d;
 
     /**
      * ShapeInfo description of the coarse cell. Needed during the
      * fast application of hanging-node constraints.
      */
-    internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>
-      shape_info_coarse;
+    internal::MatrixFreeFunctions::ShapeInfo<double> shape_info_coarse;
   };
 
   /**
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 f035321886..ec20cb9d82 100644
--- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
+++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
@@ -115,14 +115,14 @@ namespace
   /**
    * Helper class containing the cell-wise prolongation operation.
    */
-  template <int dim, typename Number>
+  template <int dim, typename Number, typename Number2>
   class CellProlongator
   {
   public:
     CellProlongator(const AlignedVector<Number> &prolongation_matrix,
                     const AlignedVector<Number> &prolongation_matrix_1d,
-                    const Number                *evaluation_data_coarse,
-                    Number                      *evaluation_data_fine)
+                    const Number2               *evaluation_data_coarse,
+                    Number2                     *evaluation_data_fine)
       : prolongation_matrix(prolongation_matrix)
       , prolongation_matrix_1d(prolongation_matrix_1d)
       , evaluation_data_coarse(evaluation_data_coarse)
@@ -170,21 +170,21 @@ namespace
   private:
     const AlignedVector<Number> &prolongation_matrix;
     const AlignedVector<Number> &prolongation_matrix_1d;
-    const Number                *evaluation_data_coarse;
-    Number                      *evaluation_data_fine;
+    const Number2               *evaluation_data_coarse;
+    Number2                     *evaluation_data_fine;
   };
 
   /**
    * Helper class containing the cell-wise restriction operation.
    */
-  template <int dim, typename Number>
+  template <int dim, typename Number, typename Number2>
   class CellRestrictor
   {
   public:
     CellRestrictor(const AlignedVector<Number> &prolongation_matrix,
                    const AlignedVector<Number> &prolongation_matrix_1d,
-                   Number                      *evaluation_data_fine,
-                   Number                      *evaluation_data_coarse)
+                   Number2                     *evaluation_data_fine,
+                   Number2                     *evaluation_data_coarse)
       : prolongation_matrix(prolongation_matrix)
       , prolongation_matrix_1d(prolongation_matrix_1d)
       , evaluation_data_fine(evaluation_data_fine)
@@ -235,8 +235,8 @@ namespace
   private:
     const AlignedVector<Number> &prolongation_matrix;
     const AlignedVector<Number> &prolongation_matrix_1d;
-    Number                      *evaluation_data_fine;
-    Number                      *evaluation_data_coarse;
+    Number2                     *evaluation_data_fine;
+    Number2                     *evaluation_data_coarse;
   };
 
   class CellProlongatorTest
@@ -2975,11 +2975,13 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
           if (needs_interpolation)
             for (int c = n_components - 1; c >= 0; --c)
               {
-                CellProlongator<dim, VectorizedArrayType> cell_prolongator(
-                  scheme.prolongation_matrix,
-                  scheme.prolongation_matrix_1d,
-                  evaluation_data_coarse.begin() + c * n_scalar_dofs_coarse,
-                  evaluation_data_fine.begin() + c * n_scalar_dofs_fine);
+                CellProlongator<dim, double, VectorizedArrayType>
+                  cell_prolongator(scheme.prolongation_matrix,
+                                   scheme.prolongation_matrix_1d,
+                                   evaluation_data_coarse.begin() +
+                                     c * n_scalar_dofs_coarse,
+                                   evaluation_data_fine.begin() +
+                                     c * n_scalar_dofs_fine);
 
                 if (scheme.prolongation_matrix_1d.size() > 0)
                   cell_transfer.run(cell_prolongator);
@@ -3173,11 +3175,13 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
           if (needs_interpolation)
             for (int c = n_components - 1; c >= 0; --c)
               {
-                CellRestrictor<dim, VectorizedArrayType> cell_restrictor(
-                  scheme.prolongation_matrix,
-                  scheme.prolongation_matrix_1d,
-                  evaluation_data_fine.begin() + c * n_scalar_dofs_fine,
-                  evaluation_data_coarse.begin() + c * n_scalar_dofs_coarse);
+                CellRestrictor<dim, double, VectorizedArrayType>
+                  cell_restrictor(scheme.prolongation_matrix,
+                                  scheme.prolongation_matrix_1d,
+                                  evaluation_data_fine.begin() +
+                                    c * n_scalar_dofs_fine,
+                                  evaluation_data_coarse.begin() +
+                                    c * n_scalar_dofs_coarse);
 
                 if (scheme.prolongation_matrix_1d.size() > 0)
                   cell_transfer.run(cell_restrictor);
@@ -3296,11 +3300,13 @@ MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>::
           if (needs_interpolation)
             for (int c = n_components - 1; c >= 0; --c)
               {
-                CellRestrictor<dim, VectorizedArrayType> cell_restrictor(
-                  scheme.restriction_matrix,
-                  scheme.restriction_matrix_1d,
-                  evaluation_data_fine.begin() + c * n_scalar_dofs_fine,
-                  evaluation_data_coarse.begin() + c * n_scalar_dofs_coarse);
+                CellRestrictor<dim, double, VectorizedArrayType>
+                  cell_restrictor(scheme.restriction_matrix,
+                                  scheme.restriction_matrix_1d,
+                                  evaluation_data_fine.begin() +
+                                    c * n_scalar_dofs_fine,
+                                  evaluation_data_coarse.begin() +
+                                    c * n_scalar_dofs_coarse);
 
                 if (scheme.restriction_matrix_1d.size() > 0)
                   cell_transfer.run(cell_restrictor);
-- 
2.39.5