]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement MatrixFreeTools::compute_matrix() 11347/head
authorPeter Munch <peterrmuench@gmail.com>
Wed, 9 Dec 2020 12:49:21 +0000 (13:49 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Fri, 11 Dec 2020 06:51:14 +0000 (07:51 +0100)
include/deal.II/matrix_free/tools.h
tests/matrix_free/compute_diagonal_01.cc
tests/matrix_free/compute_diagonal_02.cc
tests/matrix_free/compute_diagonal_03.cc
tests/matrix_free/compute_diagonal_util.h

index 6ff02f395023e57815490ccf14443ec8ef2e1f38..4d46c108c7d2dfe371f8cc3f49ef5a2f25fd3f6c 100644 (file)
@@ -72,6 +72,64 @@ namespace MatrixFreeTools
     const unsigned int first_selected_component = 0);
 
 
+  /**
+   * Compute the matrix representation of a linear operator (@p matrix), given
+   * @p matrix_free and the local cell integral operation @p local_vmult. The
+   * vector is initialized to the right size in the function.
+   *
+   * The parameters @p dof_no, @p quad_no, and @p first_selected_component are
+   * passed to the constructor of the FEEvaluation that is internally set up.
+   */
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename Number,
+            typename VectorizedArrayType,
+            typename MatrixType>
+  void
+  compute_matrix(
+    const MatrixFree<dim, Number, VectorizedArrayType> &            matrix_free,
+    const AffineConstraints<Number> &                               constraints,
+    MatrixType &                                                    matrix,
+    const std::function<void(FEEvaluation<dim,
+                                          fe_degree,
+                                          n_q_points_1d,
+                                          n_components,
+                                          Number,
+                                          VectorizedArrayType> &)> &local_vmult,
+    const unsigned int                                              dof_no  = 0,
+    const unsigned int                                              quad_no = 0,
+    const unsigned int first_selected_component = 0);
+
+  /**
+   * Same as above but with a class and a function pointer.
+   */
+  template <typename CLASS,
+            int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename Number,
+            typename VectorizedArrayType,
+            typename MatrixType>
+  void
+  compute_matrix(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const AffineConstraints<Number> &                   constraints,
+    MatrixType &                                        matrix,
+    void (CLASS::*cell_operation)(FEEvaluation<dim,
+                                               fe_degree,
+                                               n_q_points_1d,
+                                               n_components,
+                                               Number,
+                                               VectorizedArrayType> &) const,
+    const CLASS *      owning_class,
+    const unsigned int dof_no                   = 0,
+    const unsigned int quad_no                  = 0,
+    const unsigned int first_selected_component = 0);
+
+
   // implementations
 
 
@@ -509,6 +567,156 @@ namespace MatrixFreeTools
       false);
   }
 
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename Number,
+            typename VectorizedArrayType,
+            typename MatrixType>
+  void
+  compute_matrix(
+    const MatrixFree<dim, Number, VectorizedArrayType> &            matrix_free,
+    const AffineConstraints<Number> &                               constraints,
+    MatrixType &                                                    matrix,
+    const std::function<void(FEEvaluation<dim,
+                                          fe_degree,
+                                          n_q_points_1d,
+                                          n_components,
+                                          Number,
+                                          VectorizedArrayType> &)> &local_vmult,
+    const unsigned int                                              dof_no,
+    const unsigned int                                              quad_no,
+    const unsigned int first_selected_component)
+  {
+    matrix_free.template cell_loop<MatrixType, MatrixType>(
+      [&](const auto &, auto &dst, const auto &, const auto range) {
+        FEEvaluation<dim,
+                     fe_degree,
+                     n_q_points_1d,
+                     n_components,
+                     Number,
+                     VectorizedArrayType>
+          integrator(
+            matrix_free, range, dof_no, quad_no, first_selected_component);
+
+        unsigned int const dofs_per_cell = integrator.dofs_per_cell;
+
+        std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
+        std::vector<types::global_dof_index> dof_indices_mf(dofs_per_cell);
+
+        std::array<FullMatrix<typename MatrixType::value_type>,
+                   VectorizedArrayType::size()>
+          matrices;
+
+        std::fill_n(matrices.begin(),
+                    VectorizedArrayType::size(),
+                    FullMatrix<typename MatrixType::value_type>(dofs_per_cell,
+                                                                dofs_per_cell));
+
+        const auto lexicographic_numbering =
+          matrix_free
+            .get_shape_info(dof_no,
+                            quad_no,
+                            first_selected_component,
+                            integrator.get_active_fe_index(),
+                            integrator.get_active_quadrature_index())
+            .lexicographic_numbering;
+
+        for (auto cell = range.first; cell < range.second; ++cell)
+          {
+            integrator.reinit(cell);
+
+            unsigned int const n_filled_lanes =
+              matrix_free.n_active_entries_per_cell_batch(cell);
+
+            for (unsigned int v = 0; v < n_filled_lanes; ++v)
+              matrices[v] = 0.0;
+
+            for (unsigned int j = 0; j < dofs_per_cell; ++j)
+              {
+                for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                  integrator.begin_dof_values()[i] =
+                    static_cast<Number>(i == j);
+
+                local_vmult(integrator);
+
+                for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                  for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                    matrices[v](i, j) = integrator.begin_dof_values()[i][v];
+              }
+
+            for (unsigned int v = 0; v < n_filled_lanes; ++v)
+              {
+                const auto cell_v = matrix_free.get_cell_iterator(cell, v);
+
+                if (matrix_free.get_mg_level() != numbers::invalid_unsigned_int)
+                  cell_v->get_mg_dof_indices(dof_indices);
+                else
+                  cell_v->get_dof_indices(dof_indices);
+
+                for (unsigned int j = 0; j < dof_indices.size(); ++j)
+                  dof_indices_mf[j] = dof_indices[lexicographic_numbering[j]];
+
+                constraints.distribute_local_to_global(matrices[v],
+                                                       dof_indices_mf,
+                                                       dof_indices_mf,
+                                                       dst);
+              }
+          }
+      },
+      matrix,
+      matrix);
+
+    matrix.compress(VectorOperation::add);
+
+    for (unsigned int i = 0; i < matrix.m(); ++i)
+      if (matrix.get_row_length(i) > 0 && matrix(i, i) == 0.0 &&
+          constraints.is_constrained(i))
+        matrix.add(i, i, 1);
+  }
+
+  template <typename CLASS,
+            int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename Number,
+            typename VectorizedArrayType,
+            typename MatrixType>
+  void
+  compute_matrix(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const AffineConstraints<Number> &                   constraints,
+    MatrixType &                                        matrix,
+    void (CLASS::*cell_operation)(FEEvaluation<dim,
+                                               fe_degree,
+                                               n_q_points_1d,
+                                               n_components,
+                                               Number,
+                                               VectorizedArrayType> &) const,
+    const CLASS *      owning_class,
+    const unsigned int dof_no,
+    const unsigned int quad_no,
+    const unsigned int first_selected_component)
+  {
+    compute_matrix<dim,
+                   fe_degree,
+                   n_q_points_1d,
+                   n_components,
+                   Number,
+                   VectorizedArrayType,
+                   MatrixType>(matrix_free,
+                               constraints,
+                               matrix,
+                               [&](auto &feeval) {
+                                 (owning_class->*cell_operation)(feeval);
+                               },
+                               dof_no,
+                               quad_no,
+                               first_selected_component);
+  }
+
 } // namespace MatrixFreeTools
 
 DEAL_II_NAMESPACE_CLOSE
index 5013d77e72510bdc6d419b34ae28649791bee0cf..1a9d46ba9c1ef26140a3602e08754e41831aaa91 100644 (file)
@@ -65,6 +65,7 @@ test()
 
   Test<dim, fe_degree, n_points, n_components, Number, VectorizedArrayType>
     test(matrix_free,
+         constraint,
          [](FEEvaluation<dim,
                          fe_degree,
                          n_points,
index f7949e0fc7625d693df0ca8368c3df4f0689e186..25e5985b4bad7cbf2a4cf06466bb37a0011fe730 100644 (file)
@@ -104,6 +104,7 @@ test()
 
   Test<dim, fe_degree, n_points, n_components, Number, VectorizedArrayType>
     test(matrix_free,
+         constraints,
          [](FEEvaluation<dim,
                          fe_degree,
                          n_points,
index 0fd9c6e9b3176b298183868c929f729713b2b1b2..2426b6137bd82e8793c6745cd511e72dc0d4ce43 100644 (file)
@@ -63,6 +63,7 @@ test()
 
   Test<dim, fe_degree, n_points, n_components, Number, VectorizedArrayType>
     test(matrix_free,
+         constraint,
          [](FEEvaluation<dim,
                          fe_degree,
                          n_points,
index c0095bd82a2d32448773b64f9f2f0d2ea563ec19..75e9ba2a7fff31ce645a0b761cb75ef9b3535999 100644 (file)
@@ -20,6 +20,7 @@
 #include <deal.II/distributed/tria.h>
 
 #include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
 
 #include <deal.II/fe/fe_dgq.h>
 #include <deal.II/fe/fe_q.h>
@@ -30,6 +31,9 @@
 #include <deal.II/grid/grid_generator.h>
 #include <deal.II/grid/grid_out.h>
 
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+
 #include <deal.II/matrix_free/fe_evaluation.h>
 #include <deal.II/matrix_free/matrix_free.h>
 #include <deal.II/matrix_free/tools.h>
@@ -55,7 +59,8 @@ class Test
   using VectorType = LinearAlgebra::distributed::Vector<Number>;
 
 public:
-  Test(const MatrixFree<dim, Number, VectorizedArrayType> matrix_free,
+  Test(const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+       const AffineConstraints<Number> &                   constraints,
        const std::function<void(FEEvaluation<dim,
                                              fe_degree,
                                              n_points,
@@ -64,6 +69,7 @@ public:
                                              VectorizedArrayType> &)>
          &cell_operation)
     : matrix_free(matrix_free)
+    , constraints(constraints)
     , cell_operation(cell_operation)
   {}
 
@@ -71,8 +77,27 @@ public:
   do_test()
   {
     // compute diagonal with the new function
+    VectorType diagonal_global, diagonal_global_reference;
+
+    SparseMatrix<Number> A1, A2, A_ref;
+    SparsityPattern      sparsity_pattern;
+
+    const bool test_matrix = Utilities::MPI::n_mpi_processes(
+                               matrix_free.get_task_info().communicator) == 1;
+
+    if (test_matrix)
+      {
+        DynamicSparsityPattern dsp(matrix_free.get_dof_handler().n_dofs());
+        DoFTools::make_sparsity_pattern(matrix_free.get_dof_handler(),
+                                        dsp,
+                                        constraints);
+        sparsity_pattern.copy_from(dsp);
+        A1.reinit(sparsity_pattern);
+        A2.reinit(sparsity_pattern);
+        A_ref.reinit(sparsity_pattern);
+      }
+
     {
-      VectorType diagonal_global;
       MatrixFreeTools::compute_diagonal<dim,
                                         fe_degree,
                                         n_points,
@@ -87,12 +112,32 @@ public:
       deallog << diagonal_global.l2_norm() << std::endl;
     }
 
+    if (test_matrix)
+      {
+        MatrixFreeTools::compute_matrix<dim,
+                                        fe_degree,
+                                        n_points,
+                                        n_components,
+                                        Number,
+                                        VectorizedArrayType,
+                                        SparseMatrix<Number>>(
+          matrix_free, constraints, A1, [&](auto &phi) {
+            this->cell_operation(phi);
+          });
+      }
+
+    if (test_matrix)
+      {
+        MatrixFreeTools::compute_matrix(
+          matrix_free, constraints, A2, &Test::cell_function, this);
+      }
+
     // compute diagonal globally
     {
-      VectorType src, dst, temp;
+      VectorType src, temp;
 
       matrix_free.initialize_dof_vector(src);
-      matrix_free.initialize_dof_vector(dst);
+      matrix_free.initialize_dof_vector(diagonal_global_reference);
       matrix_free.initialize_dof_vector(temp);
 
       for (unsigned int i = 0; i < src.size(); i++)
@@ -104,15 +149,40 @@ public:
 
           if (src.get_partitioner()->in_local_range(i))
             {
-              dst[i] = temp[i];
-              src[i] = 0.0;
+              diagonal_global_reference[i] = temp[i];
+              src[i]                       = 0.0;
             }
+
+          if (test_matrix)
+            {
+              for (unsigned int j = 0; j < src.size(); j++)
+                if (temp[j] != 0.0)
+                  A_ref(j, i) = temp[j];
+                else if (i == j)
+                  A_ref(j, i) = 1.0;
+            }
+
           temp = 0.0;
         }
 
-      dst.print(deallog.get_file_stream());
-      deallog << dst.l2_norm() << std::endl;
+      diagonal_global_reference.print(deallog.get_file_stream());
+      deallog << diagonal_global_reference.l2_norm() << std::endl;
     }
+
+    if (test_matrix)
+      {
+        auto a1    = A1.begin();
+        auto a2    = A2.begin();
+        auto a_ref = A_ref.begin();
+
+        for (; a1 != A1.end(); ++a1, ++a2, ++a_ref)
+          {
+            Assert(std::abs(a1->value() - a_ref->value()) < 1e-6,
+                   ExcNotImplemented());
+            Assert(std::abs(a2->value() - a_ref->value()) < 1e-6,
+                   ExcNotImplemented());
+          }
+      }
   }
 
   void
@@ -137,7 +207,19 @@ public:
       }
   }
 
-  const MatrixFree<dim, Number, VectorizedArrayType> matrix_free;
+  void
+  cell_function(FEEvaluation<dim,
+                             fe_degree,
+                             n_points,
+                             n_components,
+                             Number,
+                             VectorizedArrayType> &phi) const
+  {
+    this->cell_operation(phi);
+  }
+
+  const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free;
+  const AffineConstraints<Number> &                   constraints;
   const std::function<void(FEEvaluation<dim,
                                         fe_degree,
                                         n_points,

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.