MG two-level/matrix-free infrastructure and std::complex data types.

[bangerth]

I tried to instantiate the matrix-free stuff with std::complex vectors, and utterly failed :-) Some of the minor edits necessary are captured in #16717, but the major blocker is that in mg_transfer_global_coarsening.h we have this kind of thing:

class MGTwoLevelTransfer<dim, LinearAlgebra::distributed::Vector<Number>>
  : public MGTwoLevelTransferBase<LinearAlgebra::distributed::Vector<Number>>
{
  using VectorizedArrayType = VectorizedArray<Number>;

  [...]

  struct MGTransferScheme
  {
    [...]

    /**
     * Prolongation matrix for non-tensor-product elements.
     */
    AlignedVector<VectorizedArrayType> prolongation_matrix;

    /**
     * 1d prolongation matrix for tensor-product elements.
     */
    AlignedVector<VectorizedArrayType> prolongation_matrix_1d;

    /**
     * Restriction matrix for non-tensor-product elements.
     */
    AlignedVector<VectorizedArrayType> restriction_matrix;

    /**
     * 1d restriction matrix for tensor-product elements.
     */
    AlignedVector<VectorizedArrayType> 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;
  };

Note how here we store prolongation/restriction matrices (essentially shape functions evaluated at specific points) with data type VectorizedArrayType, which depends on the type of the vectors we want to transfer. This data type is then also passed down to the matrix-free infrastructure that the MGTwoLevelTransfer class uses, specifically to CellProlongator and CellRestrictor in constraint_info.h.

This does not make conceptual sense: The shape functions (and consequently the restriction and prolongation matrices) have nothing to do with the vector type in question, and should be stored as either double or float. In particular, if one had complex-valued vectors, it does not make sense to store shape function values as std::complex. So I tried this kind of patch:

diff --git a/include/deal.II/multigrid/mg_transfer_global_coarsening.h b/include/deal.II/multigrid/mg_transfer_global_coarsening.h
index 49a0f8ccef..04c11c6378 100644
--- a/include/deal.II/multigrid/mg_transfer_global_coarsening.h
+++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.h
@@ -651,28 +651,28 @@ private:
     /**
      * Prolongation matrix for non-tensor-product elements.
      */
-    AlignedVector<VectorizedArrayType> prolongation_matrix;
+    AlignedVector<VectorizedArray<double>> prolongation_matrix;
 
     /**
      * 1d prolongation matrix for tensor-product elements.
      */
-    AlignedVector<VectorizedArrayType> prolongation_matrix_1d;
+    AlignedVector<VectorizedArray<double>> prolongation_matrix_1d;
 
     /**
      * Restriction matrix for non-tensor-product elements.
      */
-    AlignedVector<VectorizedArrayType> restriction_matrix;
+    AlignedVector<VectorizedArray<double>> restriction_matrix;
 
     /**
      * 1d restriction matrix for tensor-product elements.
      */
-    AlignedVector<VectorizedArrayType> restriction_matrix_1d;
+    AlignedVector<VectorizedArray<double>> restriction_matrix_1d;
 
     /**
      * ShapeInfo description of the coarse cell. Needed during the
      * fast application of hanging-node constraints.
      */
-    internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>
+    internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<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..4800bf94cd 100644
--- a/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
+++ b/include/deal.II/multigrid/mg_transfer_global_coarsening.templates.h
@@ -119,10 +119,11 @@ namespace
   class CellProlongator
   {
   public:
-    CellProlongator(const AlignedVector<Number> &prolongation_matrix,
-                    const AlignedVector<Number> &prolongation_matrix_1d,
-                    const Number                *evaluation_data_coarse,
-                    Number                      *evaluation_data_fine)
+    CellProlongator(
+      const AlignedVector<VectorizedArray<double>> &prolongation_matrix,
+      const AlignedVector<VectorizedArray<double>> &prolongation_matrix_1d,
+      const Number                                 *evaluation_data_coarse,
+      Number                                       *evaluation_data_fine)
       : prolongation_matrix(prolongation_matrix)
       , prolongation_matrix_1d(prolongation_matrix_1d)
       , evaluation_data_coarse(evaluation_data_coarse)
@@ -168,10 +169,10 @@ namespace
     }
 
   private:
-    const AlignedVector<Number> &prolongation_matrix;
-    const AlignedVector<Number> &prolongation_matrix_1d;
-    const Number                *evaluation_data_coarse;
-    Number                      *evaluation_data_fine;
+    const AlignedVector<VectorizedArray<double>> &prolongation_matrix;
+    const AlignedVector<VectorizedArray<double>> &prolongation_matrix_1d;
+    const Number                                 *evaluation_data_coarse;
+    Number                                       *evaluation_data_fine;
   };
 
   /**
@@ -181,10 +182,11 @@ namespace
   class CellRestrictor
   {
   public:
-    CellRestrictor(const AlignedVector<Number> &prolongation_matrix,
-                   const AlignedVector<Number> &prolongation_matrix_1d,
-                   Number                      *evaluation_data_fine,
-                   Number                      *evaluation_data_coarse)
+    CellRestrictor(
+      const AlignedVector<VectorizedArray<double>> &prolongation_matrix,
+      const AlignedVector<VectorizedArray<double>> &prolongation_matrix_1d,
+      Number                                       *evaluation_data_fine,
+      Number                                       *evaluation_data_coarse)
       : prolongation_matrix(prolongation_matrix)
       , prolongation_matrix_1d(prolongation_matrix_1d)
       , evaluation_data_fine(evaluation_data_fine)
@@ -233,10 +235,10 @@ namespace
     }
 
   private:
-    const AlignedVector<Number> &prolongation_matrix;
-    const AlignedVector<Number> &prolongation_matrix_1d;
-    Number                      *evaluation_data_fine;
-    Number                      *evaluation_data_coarse;
+    const AlignedVector<VectorizedArray<double>> &prolongation_matrix;
+    const AlignedVector<VectorizedArray<double>> &prolongation_matrix_1d;
+    Number                                       *evaluation_data_fine;
+    Number                                       *evaluation_data_coarse;
   };
 
   class CellProlongatorTest

I think this is conceptually the right approach, but it also helped me understand where the current design comes from: We want to do arithmetic between these matrices and vectorized arrays of DoF values, and we want to do that between vectorized arrays of the same length (VectorizedArray<double> may not have the same size of VectorizedArray<float>). Because the matrix-free framework does that, the patch above does not go very far, with errors of the following kind:

/home/bangerth/p/deal.II/1/dealii/include/deal.II/matrix_free/fe_point_evaluation.h:2378:29: error: cannot convert ‘const std::complex<double>’ to ‘dealii::Tensor<1, 1>::value_type’ {aka ‘double’} in assignment
 2378 |           unit_points[v][d] =
      |           ~~~~~~~~~~~~~~~~~~^
 2379 |             this->unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal];

and

/home/bangerth/p/deal.II/1/dealii/include/deal.II/matrix_free/tensor_product_kernels.h:155:51: note:   ‘const dealii::VectorizedArray<double, 4>’ i
s not derived from ‘const dealii::SymmetricTensor<rank, dim, T>’
  155 |               res0 += matrix[col * n_columns + i] * x[i];
      |                       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~^~~~

I gave up at this point, but I still think that it would be nice if we could make these kinds of things work.

[kronbichler]

I agree that we should be able to make this work. I have some ideas that could help addressing this but limited time this week. I think the first step is to switch the type AlignedVector<VectorizedArray<double>> prolongation_matrix; for the various arrays to AlignedVector<double> prolongation_matrix;. I will post that PR in a second. It will not resolve the other problems, but I think that is a better starting point for further investigations.

[kronbichler]

@bangerth can you try with #16721 and see if that brings you closer to the goal? I anticipate further problems, but I think taking one problem at a time will allow us to eventually get to the desired support.

[bangerth]

I will try tonight (I have meetings 9am - 5pm today and still need to read part of a thesis for a defense today :-( ). I should say that this is not super urgent -- I believe that for what I'm really after (using MGTwoLevelTransfer for interpolating up and down between meshes), everything works for real-valued vectors, and that's all my student needs.