]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Modularized MatrixFreeTools::compute_diagonal() 11802/head
authorPeter Munch <peterrmuench@gmail.com>
Wed, 24 Feb 2021 20:54:28 +0000 (21:54 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Fri, 26 Feb 2021 06:34:30 +0000 (07:34 +0100)
include/deal.II/matrix_free/tools.h

index 55376a701ba4e1032b21deb8afbc75f07818d1cd..d3498fba298e215e6d512ecb2d87dbec813b93d8 100644 (file)
@@ -247,110 +247,130 @@ namespace MatrixFreeTools
       std::vector<unsigned int> col;
       std::vector<Number>       val;
     };
-  } // namespace internal
 
-  template <int dim,
-            int fe_degree,
-            int n_q_points_1d,
-            int n_components,
-            typename Number,
-            typename VectorizedArrayType>
-  void
-  compute_diagonal(
-    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
-    LinearAlgebra::distributed::Vector<Number> &        diagonal_global,
-    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)
-  {
-    using VectorType = LinearAlgebra::distributed::Vector<Number>;
+    template <int dim,
+              int fe_degree,
+              int n_q_points_1d,
+              int n_components,
+              typename Number,
+              typename VectorizedArrayType>
+    class ComputeDiagonalHelper
+    {
+    public:
+      static const unsigned int n_lanes = VectorizedArrayType::size();
+
+      ComputeDiagonalHelper(FEEvaluation<dim,
+                                         fe_degree,
+                                         n_q_points_1d,
+                                         n_components,
+                                         Number,
+                                         VectorizedArrayType> &phi)
+        : phi(phi)
+      {}
 
-    // initialize vector
-    matrix_free.initialize_dof_vector(diagonal_global);
+      void
+      reinit(const unsigned int cell)
+      {
+        this->phi.reinit(cell);
+        // STEP 1: get relevant information from FEEvaluation
+        const unsigned int first_selected_component =
+          phi.get_first_selected_component();
+        const auto &       dof_info        = phi.get_dof_info();
+        const unsigned int n_fe_components = dof_info.start_components.back();
+        const unsigned int dofs_per_component = phi.dofs_per_component;
+        const auto &       matrix_free        = phi.get_matrix_free();
+
+        const unsigned int n_lanes_filled =
+          matrix_free.n_active_entries_per_cell_batch(cell);
+
+        std::array<const unsigned int *, n_lanes> dof_indices{};
+        {
+          for (unsigned int v = 0; v < n_lanes_filled; ++v)
+            dof_indices[v] =
+              dof_info.dof_indices.data() +
+              dof_info
+                .row_starts[(cell * n_lanes + v) * n_fe_components +
+                            first_selected_component]
+                .first;
+        }
 
-    int dummy;
+        // STEP 2: setup CSR storage of transposed locally-relevant
+        //   constraint matrix
+        c_pools = std::array<internal::LocalCSR<Number>, n_lanes>();
 
-    matrix_free.template cell_loop<VectorType, int>(
-      [&](const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
-          LinearAlgebra::distributed::Vector<Number> &        diagonal_global,
-          const int &,
-          const std::pair<unsigned int, unsigned int> &range) mutable {
-        FEEvaluation<dim,
-                     fe_degree,
-                     n_q_points_1d,
-                     n_components,
-                     Number,
-                     VectorizedArrayType>
-          phi(matrix_free, range, dof_no, quad_no, first_selected_component);
+        for (unsigned int v = 0; v < n_lanes_filled; ++v)
+          {
+            unsigned int index_indicators, next_index_indicators;
+
+            index_indicators =
+              dof_info
+                .row_starts[(cell * n_lanes + v) * n_fe_components +
+                            first_selected_component]
+                .second;
+            next_index_indicators =
+              dof_info
+                .row_starts[(cell * n_lanes + v) * n_fe_components +
+                            first_selected_component + 1]
+                .second;
+
+            // STEP 2a: setup locally-relevant constraint matrix in a
+            //   coordinate list (COO)
+            std::vector<std::tuple<unsigned int, unsigned int, Number>>
+              locally_relevant_constrains; // (constrained local index,
+                                           // global index of dof which
+                                           // constrains, weight)
+
+            if (n_components == 1 || n_fe_components == 1)
+              {
+                AssertDimension(n_components,
+                                1); // TODO: currently no block vector supported
 
-        const unsigned int n_lanes = VectorizedArrayType::size();
+                unsigned int ind_local = 0;
+                for (; index_indicators != next_index_indicators;
+                     ++index_indicators, ++ind_local)
+                  {
+                    const std::pair<unsigned short, unsigned short> indicator =
+                      dof_info.constraint_indicator[index_indicators];
 
-        // local storage: buffer so that we access the global vector once
-        // note: may be larger then dofs_per_cell in the presence of
-        // constraints!
-        std::array<std::vector<Number>, n_lanes> diagonals_local_constrained;
+                    for (unsigned int j = 0; j < indicator.first;
+                         ++j, ++ind_local)
+                      locally_relevant_constrains.emplace_back(
+                        ind_local, dof_indices[v][j], 1.0);
 
-        for (unsigned int cell = range.first; cell < range.second; ++cell)
-          {
-            // STEP 1: get relevant information from FEEvaluation
-            const auto &       dof_info = matrix_free.get_dof_info(dof_no);
-            const unsigned int n_fe_components =
-              dof_info.start_components.back();
-            const unsigned int dofs_per_component = phi.dofs_per_component;
+                    dof_indices[v] += indicator.first;
 
-            const unsigned int n_lanes_filled =
-              matrix_free.n_active_entries_per_cell_batch(cell);
+                    const Number *data_val =
+                      matrix_free.constraint_pool_begin(indicator.second);
+                    const Number *end_pool =
+                      matrix_free.constraint_pool_end(indicator.second);
 
-            std::array<const unsigned int *, n_lanes> dof_indices{};
-            {
-              for (unsigned int v = 0; v < n_lanes_filled; ++v)
-                dof_indices[v] =
-                  dof_info.dof_indices.data() +
-                  dof_info
-                    .row_starts[(cell * n_lanes + v) * n_fe_components +
-                                first_selected_component]
-                    .first;
-            }
+                    for (; data_val != end_pool; ++data_val, ++dof_indices[v])
+                      locally_relevant_constrains.emplace_back(ind_local,
+                                                               *dof_indices[v],
+                                                               *data_val);
+                  }
 
-            // STEP 2: setup CSR storage of transposed locally-relevant
-            //   constraint matrix
-            std::array<internal::LocalCSR<Number>, n_lanes> c_pools;
+                AssertIndexRange(ind_local, dofs_per_component + 1);
 
-            for (unsigned int v = 0; v < n_lanes_filled; ++v)
+                for (; ind_local < dofs_per_component;
+                     ++dof_indices[v], ++ind_local)
+                  locally_relevant_constrains.emplace_back(ind_local,
+                                                           *dof_indices[v],
+                                                           1.0);
+              }
+            else
               {
-                unsigned int index_indicators, next_index_indicators;
-
-                index_indicators =
-                  dof_info
-                    .row_starts[(cell * n_lanes + v) * n_fe_components +
-                                first_selected_component]
-                    .second;
-                next_index_indicators =
-                  dof_info
-                    .row_starts[(cell * n_lanes + v) * n_fe_components +
-                                first_selected_component + 1]
-                    .second;
-
-                // STEP 2a: setup locally-relevant constraint matrix in a
-                //   coordinate list (COO)
-                std::vector<std::tuple<unsigned int, unsigned int, Number>>
-                  locally_relevant_constrains; // (constrained local index,
-                                               // global index of dof which
-                                               // constrains, weight)
-
-                if (n_components == 1 || n_fe_components == 1)
+                // case with vector-valued finite elements where all
+                // components are included in one single vector. Assumption:
+                // first come all entries to the first component, then all
+                // entries to the second one, and so on. This is ensured by
+                // the way MatrixFree reads out the indices.
+                for (unsigned int comp = 0; comp < n_components; ++comp)
                   {
-                    AssertDimension(
-                      n_components,
-                      1); // TODO: currently no block vector supported
-
                     unsigned int ind_local = 0;
+
+                    // check whether there is any constraint on the current
+                    // cell
                     for (; index_indicators != next_index_indicators;
                          ++index_indicators, ++ind_local)
                       {
@@ -358,11 +378,13 @@ namespace MatrixFreeTools
                           indicator =
                             dof_info.constraint_indicator[index_indicators];
 
+                        // run through values up to next constraint
                         for (unsigned int j = 0; j < indicator.first;
                              ++j, ++ind_local)
                           locally_relevant_constrains.emplace_back(
-                            ind_local, dof_indices[v][j], 1.0);
-
+                            comp * dofs_per_component + ind_local,
+                            dof_indices[v][j],
+                            1.0);
                         dof_indices[v] += indicator.first;
 
                         const Number *data_val =
@@ -373,218 +395,257 @@ namespace MatrixFreeTools
                         for (; data_val != end_pool;
                              ++data_val, ++dof_indices[v])
                           locally_relevant_constrains.emplace_back(
-                            ind_local, *dof_indices[v], *data_val);
+                            comp * dofs_per_component + ind_local,
+                            *dof_indices[v],
+                            *data_val);
                       }
 
                     AssertIndexRange(ind_local, dofs_per_component + 1);
 
+                    // get the dof values past the last constraint
                     for (; ind_local < dofs_per_component;
                          ++dof_indices[v], ++ind_local)
-                      locally_relevant_constrains.emplace_back(ind_local,
-                                                               *dof_indices[v],
-                                                               1.0);
-                  }
-                else
-                  {
-                    // case with vector-valued finite elements where all
-                    // components are included in one single vector. Assumption:
-                    // first come all entries to the first component, then all
-                    // entries to the second one, and so on. This is ensured by
-                    // the way MatrixFree reads out the indices.
-                    for (unsigned int comp = 0; comp < n_components; ++comp)
-                      {
-                        unsigned int ind_local = 0;
-
-                        // check whether there is any constraint on the current
-                        // cell
-                        for (; index_indicators != next_index_indicators;
-                             ++index_indicators, ++ind_local)
-                          {
-                            const std::pair<unsigned short, unsigned short>
-                              indicator =
-                                dof_info.constraint_indicator[index_indicators];
-
-                            // run through values up to next constraint
-                            for (unsigned int j = 0; j < indicator.first;
-                                 ++j, ++ind_local)
-                              locally_relevant_constrains.emplace_back(
-                                comp * dofs_per_component + ind_local,
-                                dof_indices[v][j],
-                                1.0);
-                            dof_indices[v] += indicator.first;
-
-                            const Number *data_val =
-                              matrix_free.constraint_pool_begin(
-                                indicator.second);
-                            const Number *end_pool =
-                              matrix_free.constraint_pool_end(indicator.second);
-
-                            for (; data_val != end_pool;
-                                 ++data_val, ++dof_indices[v])
-                              locally_relevant_constrains.emplace_back(
-                                comp * dofs_per_component + ind_local,
-                                *dof_indices[v],
-                                *data_val);
-                          }
-
-                        AssertIndexRange(ind_local, dofs_per_component + 1);
-
-                        // get the dof values past the last constraint
-                        for (; ind_local < dofs_per_component;
-                             ++dof_indices[v], ++ind_local)
-                          locally_relevant_constrains.emplace_back(
-                            comp * dofs_per_component + ind_local,
-                            *dof_indices[v],
-                            1.0);
+                      locally_relevant_constrains.emplace_back(
+                        comp * dofs_per_component + ind_local,
+                        *dof_indices[v],
+                        1.0);
 
-                        if (comp + 1 < n_components)
-                          {
-                            next_index_indicators =
-                              dof_info
-                                .row_starts[(cell * n_lanes + v) *
-                                              n_fe_components +
-                                            first_selected_component + comp + 2]
-                                .second;
-                          }
+                    if (comp + 1 < n_components)
+                      {
+                        next_index_indicators =
+                          dof_info
+                            .row_starts[(cell * n_lanes + v) * n_fe_components +
+                                        first_selected_component + comp + 2]
+                            .second;
                       }
                   }
+              }
 
-                // STEP 2b: transpose COO
-
-                // presort vector for transposed access
-                std::sort(locally_relevant_constrains.begin(),
-                          locally_relevant_constrains.end(),
-                          [](const auto &a, const auto &b) {
-                            if (std::get<1>(a) < std::get<1>(b))
-                              return true;
-                            return (std::get<1>(a) == std::get<1>(b)) &&
-                                   (std::get<0>(a) < std::get<0>(b));
-                          });
-
-                // make sure that all entries are unique
-                locally_relevant_constrains.erase(
-                  unique(locally_relevant_constrains.begin(),
-                         locally_relevant_constrains.end(),
-                         [](const auto &a, const auto &b) {
-                           return (std::get<1>(a) == std::get<1>(b)) &&
-                                  (std::get<0>(a) == std::get<0>(b));
-                         }),
-                  locally_relevant_constrains.end());
-
-                // STEP 2c: translate COO to CRS
-                auto &c_pool = c_pools[v];
+            // STEP 2b: transpose COO
+
+            // presort vector for transposed access
+            std::sort(locally_relevant_constrains.begin(),
+                      locally_relevant_constrains.end(),
+                      [](const auto &a, const auto &b) {
+                        if (std::get<1>(a) < std::get<1>(b))
+                          return true;
+                        return (std::get<1>(a) == std::get<1>(b)) &&
+                               (std::get<0>(a) < std::get<0>(b));
+                      });
+
+            // make sure that all entries are unique
+            locally_relevant_constrains.erase(
+              unique(locally_relevant_constrains.begin(),
+                     locally_relevant_constrains.end(),
+                     [](const auto &a, const auto &b) {
+                       return (std::get<1>(a) == std::get<1>(b)) &&
+                              (std::get<0>(a) == std::get<0>(b));
+                     }),
+              locally_relevant_constrains.end());
+
+            // STEP 2c: translate COO to CRS
+            auto &c_pool = c_pools[v];
+            {
+              if (locally_relevant_constrains.size() > 0)
+                c_pool.row_lid_to_gid.emplace_back(
+                  std::get<1>(locally_relevant_constrains.front()));
+              for (const auto &j : locally_relevant_constrains)
                 {
-                  if (locally_relevant_constrains.size() > 0)
-                    c_pool.row_lid_to_gid.emplace_back(
-                      std::get<1>(locally_relevant_constrains.front()));
-                  for (const auto &j : locally_relevant_constrains)
+                  if (c_pool.row_lid_to_gid.back() != std::get<1>(j))
                     {
-                      if (c_pool.row_lid_to_gid.back() != std::get<1>(j))
-                        {
-                          c_pool.row_lid_to_gid.push_back(std::get<1>(j));
-                          c_pool.row.push_back(c_pool.val.size());
-                        }
-
-                      c_pool.col.emplace_back(std::get<0>(j));
-                      c_pool.val.emplace_back(std::get<2>(j));
+                      c_pool.row_lid_to_gid.push_back(std::get<1>(j));
+                      c_pool.row.push_back(c_pool.val.size());
                     }
 
-                  if (c_pool.val.size() > 0)
-                    c_pool.row.push_back(c_pool.val.size());
+                  c_pool.col.emplace_back(std::get<0>(j));
+                  c_pool.val.emplace_back(std::get<2>(j));
                 }
-              }
 
-            // STEP 3: compute element matrix A_e, apply
-            //   locally-relevant constraints C_e^T * A_e * C_e, and get the
-            //   the diagonal entry
-            //     (C_e^T * A_e * C_e)(i,i)
-            //   or
-            //     C_e^T(i,:) * A_e * C_e(:,i).
-            //
-            //   Since, we compute the element matrix column-by-column and as a
-            //   result never actually have the full element matrix, we actually
-            //   perform following steps:
-            //    1) loop over all columns of the element matrix
-            //     a) compute column i
-            //     b) compute for each j (rows of C_e^T):
-            //          (C_e^T(j,:) * A_e(:,i)) * C_e(i,j)
-            //       or
-            //          (C_e^T(j,:) * A_e(:,i)) * C_e^T(j,i)
-            //       This gives a contribution the the j-th entry of the
-            //       locally-relevant diagonal and comprises the multiplication
-            //       by the locally-relevant constraint matrix from the left and
-            //       the right. There is no contribution to the j-th vector
-            //       entry if the j-th row of C_e^T is empty or C_e^T(j,i) is
-            //       zero.
-
-            // set size locally-relevant diagonal
-            for (unsigned int v = 0; v < n_lanes_filled; ++v)
-              diagonals_local_constrained[v].assign(
-                c_pools[v].row_lid_to_gid.size(), Number(0.0));
-
-            phi.reinit(cell);
-
-            // loop over all columns of element stiffness matrix
-            for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
+              if (c_pool.val.size() > 0)
+                c_pool.row.push_back(c_pool.val.size());
+            }
+          }
+        // STEP 3: compute element matrix A_e, apply
+        //   locally-relevant constraints C_e^T * A_e * C_e, and get the
+        //   the diagonal entry
+        //     (C_e^T * A_e * C_e)(i,i)
+        //   or
+        //     C_e^T(i,:) * A_e * C_e(:,i).
+        //
+        //   Since, we compute the element matrix column-by-column and as a
+        //   result never actually have the full element matrix, we actually
+        //   perform following steps:
+        //    1) loop over all columns of the element matrix
+        //     a) compute column i
+        //     b) compute for each j (rows of C_e^T):
+        //          (C_e^T(j,:) * A_e(:,i)) * C_e(i,j)
+        //       or
+        //          (C_e^T(j,:) * A_e(:,i)) * C_e^T(j,i)
+        //       This gives a contribution the the j-th entry of the
+        //       locally-relevant diagonal and comprises the multiplication
+        //       by the locally-relevant constraint matrix from the left and
+        //       the right. There is no contribution to the j-th vector
+        //       entry if the j-th row of C_e^T is empty or C_e^T(j,i) is
+        //       zero.
+
+        // set size locally-relevant diagonal
+        for (unsigned int v = 0; v < n_lanes_filled; ++v)
+          diagonals_local_constrained[v].assign(
+            c_pools[v].row_lid_to_gid.size(), Number(0.0));
+      }
+
+      void
+      prepare_basis_vector(const unsigned int i)
+      {
+        this->i = i;
+
+        // compute i-th column of element stiffness matrix:
+        // this could be simply performed as done at the moment with
+        // matrix-free operator evaluation applied to a ith-basis vector
+        for (unsigned int j = 0; j < phi.dofs_per_cell; ++j)
+          phi.begin_dof_values()[j] = static_cast<Number>(i == j);
+      }
+
+      void
+      submit()
+      {
+        const auto ith_column = phi.begin_dof_values();
+
+        // apply local constraint matrix from left and from right:
+        // loop over all rows of transposed constrained matrix
+        for (unsigned int v = 0;
+             v < phi.get_matrix_free().n_active_entries_per_cell_batch(
+                   phi.get_current_cell_index());
+             ++v)
+          {
+            const auto &c_pool = c_pools[v];
+
+            for (unsigned int j = 0; j < c_pool.row.size() - 1; ++j)
               {
-                // compute i-th column of element stiffness matrix:
-                // this could be simply performed as done at the moment with
-                // matrix-free operator evaluation applied to a ith-basis vector
+                // check if the result will be zero, so that we can skip
+                // the following computations -> binary search
+                const auto scale_iterator =
+                  std::lower_bound(c_pool.col.begin() + c_pool.row[j],
+                                   c_pool.col.begin() + c_pool.row[j + 1],
+                                   i);
+
+                // explanation: j-th row of C_e^T is empty (see above)
+                if (scale_iterator == c_pool.col.begin() + c_pool.row[j + 1])
+                  continue;
+
+                // explanation: C_e^T(j,i) is zero (see above)
+                if (*scale_iterator != i)
+                  continue;
+
+                // apply constraint matrix from the left
+                Number temp = 0.0;
+                for (unsigned int k = c_pool.row[j]; k < c_pool.row[j + 1]; ++k)
+                  temp += c_pool.val[k] * ith_column[c_pool.col[k]][v];
+
+                // apply constraint matrix from the right
+                diagonals_local_constrained[v][j] +=
+                  temp *
+                  c_pool.val[std::distance(c_pool.col.begin(), scale_iterator)];
+              }
+          }
+      }
 
-                for (unsigned int j = 0; j < phi.dofs_per_cell; ++j)
-                  phi.begin_dof_values()[j] = static_cast<Number>(i == j);
+      void
+      distribute_local_to_global(
+        LinearAlgebra::distributed::Vector<Number> &diagonal_global)
+      {
+        // STEP 4: assembly results: add into global vector
+        for (unsigned int v = 0;
+             v < phi.get_matrix_free().n_active_entries_per_cell_batch(
+                   phi.get_current_cell_index());
+             ++v)
+          for (unsigned int j = 0; j < c_pools[v].row.size() - 1; ++j)
+            ::dealii::internal::vector_access_add(
+              diagonal_global,
+              c_pools[v].row_lid_to_gid[j],
+              diagonals_local_constrained[v][j]);
+      }
 
-                local_vmult(phi);
+    private:
+      FEEvaluation<dim,
+                   fe_degree,
+                   n_q_points_1d,
+                   n_components,
+                   Number,
+                   VectorizedArrayType> &phi;
 
-                const auto ith_column = phi.begin_dof_values();
+      unsigned int i;
 
-                // apply local constraint matrix from left and from right:
-                // loop over all rows of transposed constrained matrix
-                for (unsigned int v = 0; v < n_lanes_filled; ++v)
-                  {
-                    const auto &c_pool = c_pools[v];
+      std::array<internal::LocalCSR<Number>, n_lanes> c_pools;
 
-                    for (unsigned int j = 0; j < c_pool.row.size() - 1; ++j)
-                      {
-                        // check if the result will be zero, so that we can skip
-                        // the following computations -> binary search
-                        const auto scale_iterator =
-                          std::lower_bound(c_pool.col.begin() + c_pool.row[j],
-                                           c_pool.col.begin() +
-                                             c_pool.row[j + 1],
-                                           i);
-
-                        // explanation: j-th row of C_e^T is empty (see above)
-                        if (scale_iterator ==
-                            c_pool.col.begin() + c_pool.row[j + 1])
-                          continue;
-
-                        // explanation: C_e^T(j,i) is zero (see above)
-                        if (*scale_iterator != i)
-                          continue;
-
-                        // apply constraint matrix from the left
-                        Number temp = 0.0;
-                        for (unsigned int k = c_pool.row[j];
-                             k < c_pool.row[j + 1];
-                             ++k)
-                          temp += c_pool.val[k] * ith_column[c_pool.col[k]][v];
-
-                        // apply constraint matrix from the right
-                        diagonals_local_constrained[v][j] +=
-                          temp * c_pool.val[std::distance(c_pool.col.begin(),
-                                                          scale_iterator)];
-                      }
-                  }
+      // local storage: buffer so that we access the global vector once
+      // note: may be larger then dofs_per_cell in the presence of
+      // constraints!
+      std::array<std::vector<Number>, n_lanes> diagonals_local_constrained;
+    };
+
+  } // namespace internal
+
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename Number,
+            typename VectorizedArrayType>
+  void
+  compute_diagonal(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    LinearAlgebra::distributed::Vector<Number> &        diagonal_global,
+    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)
+  {
+    using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+    // initialize vector
+    matrix_free.initialize_dof_vector(diagonal_global, dof_no);
+
+    int dummy;
+
+    matrix_free.template cell_loop<VectorType, int>(
+      [&](const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+          LinearAlgebra::distributed::Vector<Number> &        diagonal_global,
+          const int &,
+          const std::pair<unsigned int, unsigned int> &range) mutable {
+        FEEvaluation<dim,
+                     fe_degree,
+                     n_q_points_1d,
+                     n_components,
+                     Number,
+                     VectorizedArrayType>
+          phi(matrix_free, range, dof_no, quad_no, first_selected_component);
+
+        internal::ComputeDiagonalHelper<dim,
+                                        fe_degree,
+                                        n_q_points_1d,
+                                        n_components,
+                                        Number,
+                                        VectorizedArrayType>
+          helper(phi);
+
+        for (unsigned int cell = range.first; cell < range.second; ++cell)
+          {
+            helper.reinit(cell);
+
+            for (unsigned int i = 0; i < phi.dofs_per_cell; ++i)
+              {
+                helper.prepare_basis_vector(i);
+                local_vmult(phi);
+                helper.submit();
               }
 
-            // STEP 4: assembly results: add into global vector
-            for (unsigned int v = 0; v < n_lanes_filled; ++v)
-              for (unsigned int j = 0; j < c_pools[v].row.size() - 1; ++j)
-                ::dealii::internal::vector_access_add(
-                  diagonal_global,
-                  c_pools[v].row_lid_to_gid[j],
-                  diagonals_local_constrained[v][j]);
+            helper.distribute_local_to_global(diagonal_global);
           }
       },
       diagonal_global,

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.