Precompute inverse diagonal (FDM)

diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h
index 47e6a1aa776a0ebe070a4519ec1549b96ead91b2..dcb542fffbee4707dac0fd87952fc58f8a2f1f8c 100644 (file)
--- a/include/deal.II/lac/tensor_product_matrix.h
+++ b/include/deal.II/lac/tensor_product_matrix.h
@@ -376,12 +376,18 @@ public:
     /**
      * Constructor.
      */
-    AdditionalData(const bool compress_matrices = true);
+    AdditionalData(const bool compress_matrices           = true,
+                   const bool precompute_inverse_diagonal = true);
 
     /**
      * Try to compress internal matrices. Default: true.
      */
     bool compress_matrices;
+
+    /**
+     * Precompute inverse diagonal.
+     */
+    bool precompute_inverse_diagonal;
   };
 
   /**
@@ -445,6 +451,11 @@ private:
    */
   const bool compress_matrices;
 
+  /**
+   * Precompute inverse diagonal.
+   */
+  const bool precompute_inverse_diagonal;
+
   /**
    * Container used to collect 1d matrices if no compression is
    * requested. The memory is freed during finalize().
@@ -490,6 +501,11 @@ private:
    */
   AlignedVector<Number> eigenvalues;
 
+  /**
+   * Vector of inverted eigenvalues.
+   */
+  AlignedVector<Number> inverted_eigenvalues;
+
   /**
    * Pointer into mass_matrices, derivative_matrices, and eigenvalues.
    */
@@ -499,6 +515,11 @@ private:
    * Pointer into mass_matrices, derivative_matrices, and eigenvalues.
    */
   std::vector<unsigned int> matrix_ptr;
+
+  /**
+   * Number of rows in 1 of each cell.
+   */
+  std::vector<unsigned int> vector_n_rows_1d;
 };
 
 
@@ -810,7 +831,8 @@ namespace internal
                   AlignedVector<Number> &tmp,
                   const unsigned int     n_rows_1d_non_templated,
                   const std::array<const Number *, dim> &eigenvectors,
-                  const std::array<const Number *, dim> &eigenvalues)
+                  const std::array<const Number *, dim> &eigenvalues,
+                  const Number *inverted_eigenvalues = nullptr)
     {
       const unsigned int n_rows_1d = n_rows_1d_templated == 0 ?
                                        n_rows_1d_non_templated :
@@ -836,8 +858,13 @@ namespace internal
         {
           const Number *S = eigenvectors[0];
           eval.template apply<0, true, false>(S, src, t);
+
           for (unsigned int i = 0; i < n_rows_1d; ++i)
-            t[i] /= eigenvalues[0][i];
+            if (inverted_eigenvalues)
+              t[i] *= inverted_eigenvalues[i];
+            else
+              t[i] /= eigenvalues[0][i];
+
           eval.template apply<0, false, false>(S, t, dst);
         }
 
@@ -847,9 +874,14 @@ namespace internal
           const Number *S1 = eigenvectors[1];
           eval.template apply<0, true, false>(S0, src, t);
           eval.template apply<1, true, false>(S1, t, dst);
+
           for (unsigned int i1 = 0, c = 0; i1 < n_rows_1d; ++i1)
             for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
-              dst[c] /= (eigenvalues[1][i1] + eigenvalues[0][i0]);
+              if (inverted_eigenvalues)
+                dst[c] *= inverted_eigenvalues[c];
+              else
+                dst[c] /= (eigenvalues[1][i1] + eigenvalues[0][i0]);
+
           eval.template apply<0, false, false>(S0, dst, t);
           eval.template apply<1, false, false>(S1, t, dst);
         }
@@ -862,11 +894,16 @@ namespace internal
           eval.template apply<0, true, false>(S0, src, t);
           eval.template apply<1, true, false>(S1, t, dst);
           eval.template apply<2, true, false>(S2, dst, t);
+
           for (unsigned int i2 = 0, c = 0; i2 < n_rows_1d; ++i2)
             for (unsigned int i1 = 0; i1 < n_rows_1d; ++i1)
               for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
-                t[c] /= (eigenvalues[2][i2] + eigenvalues[1][i1] +
-                         eigenvalues[0][i0]);
+                if (inverted_eigenvalues)
+                  t[c] *= inverted_eigenvalues[c];
+                else
+                  t[c] /= (eigenvalues[2][i2] + eigenvalues[1][i1] +
+                           eigenvalues[0][i0]);
+
           eval.template apply<0, false, false>(S0, t, dst);
           eval.template apply<1, false, false>(S1, dst, t);
           eval.template apply<2, false, false>(S2, t, dst);
@@ -896,7 +933,8 @@ namespace internal
                          AlignedVector<Number> &                tmp,
                          const unsigned int                     n_rows_1d,
                          const std::array<const Number *, dim> &eigenvectors,
-                         const std::array<const Number *, dim> &eigenvalues);
+                         const std::array<const Number *, dim> &eigenvalues,
+                         const Number *inverted_eigenvalues = nullptr);
   } // namespace TensorProductMatrixSymmetricSum
 } // namespace internal
 
@@ -1073,8 +1111,10 @@ TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::reinit(
 
 template <int dim, typename Number, int n_rows_1d>
 TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::
-  AdditionalData::AdditionalData(const bool compress_matrices)
+  AdditionalData::AdditionalData(const bool compress_matrices,
+                                 const bool precompute_inverse_diagonal)
   : compress_matrices(compress_matrices)
+  , precompute_inverse_diagonal(precompute_inverse_diagonal)
 {}
 
 
@@ -1084,6 +1124,7 @@ TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::
   TensorProductMatrixSymmetricSumCollection(
     const AdditionalData &additional_data)
   : compress_matrices(additional_data.compress_matrices)
+  , precompute_inverse_diagonal(additional_data.precompute_inverse_diagonal)
 {}
 
 
@@ -1248,7 +1289,7 @@ TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::finalize()
     }
   else
     {
-      // case 2) compression requested but none possible
+      // case 3) compress
 
       this->vector_ptr.resize(cache.size() + 1);
       this->matrix_ptr.resize(cache.size() + 1);
@@ -1277,6 +1318,137 @@ TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::finalize()
 
       cache.clear();
     }
+
+  if (precompute_inverse_diagonal)
+    {
+      if (dim == 1)
+        {
+          // 1D case: simply invert 1D eigenvalues
+          for (unsigned int i = 0; i < this->eigenvalues.size(); ++i)
+            this->eigenvalues[i] = Number(1.0) / this->eigenvalues[i];
+          std::swap(this->inverted_eigenvalues, eigenvalues);
+        }
+      else
+        {
+          // 2D and 3D case: we have 2 or 3 1d eigenvalues so that we
+          // need to combine these
+
+          // step 1) if eigenvalues/eigenvectors are compressed, we
+          // need to compress the diagonal (the combination of ev
+          // indices) as well. This is an optional step.
+          std::vector<unsigned int> indices_ev;
+
+          if (indices.size() > 0)
+            {
+              // 1a) create cache (ev indics -> diag index)
+              const unsigned int n_cells = indices.size() / dim;
+              std::map<std::array<unsigned int, dim>, unsigned int> cache_ev;
+              std::vector<unsigned int> cache_ev_idx(n_cells);
+
+              for (unsigned int i = 0, c = 0; i < n_cells; ++i)
+                {
+                  std::array<unsigned int, dim> id;
+
+                  for (unsigned int d = 0; d < dim; ++d, ++c)
+                    id[d] = indices[c];
+
+                  const auto id_ptr = cache_ev.find(id);
+
+                  if (id_ptr == cache_ev.end())
+                    {
+                      const auto size = cache_ev.size();
+                      cache_ev_idx[i] = size;
+                      cache_ev[id]    = size;
+                    }
+                  else
+                    {
+                      cache_ev_idx[i] = id_ptr->second;
+                    }
+                }
+
+              // 1b) store diagonal indices for each cell
+              std::vector<unsigned int> new_indices;
+              new_indices.reserve(indices.size() / dim * (dim + 1));
+
+              for (unsigned int i = 0, c = 0; i < n_cells; ++i)
+                {
+                  for (unsigned int d = 0; d < dim; ++d, ++c)
+                    new_indices.push_back(indices[c]);
+                  new_indices.push_back(cache_ev_idx[i]);
+                }
+
+              // 1c) transpose cache (diag index -> ev indices)
+              indices_ev.resize(cache_ev.size() * dim);
+              for (const auto &entry : cache_ev)
+                for (unsigned int d = 0; d < dim; ++d)
+                  indices_ev[entry.second * dim + d] = entry.first[d];
+
+              std::swap(this->indices, new_indices);
+            }
+
+          // step 2) allocate memory and set pointers
+          const unsigned int n_diag =
+            ((indices_ev.size() > 0) ? indices_ev.size() :
+                                       (matrix_ptr.size() - 1)) /
+            dim;
+
+          std::vector<unsigned int> new_vector_ptr(n_diag + 1, 0);
+          std::vector<unsigned int> new_vector_n_rows_1d(n_diag, 0);
+
+          for (unsigned int i = 0; i < n_diag; ++i)
+            {
+              const unsigned int c = (indices_ev.size() > 0) ?
+                                       indices_ev[dim * i + 0] :
+                                       (dim * i + 0);
+
+              const unsigned int n_rows = vector_ptr[c + 1] - vector_ptr[c];
+
+              new_vector_n_rows_1d[i] = n_rows;
+              new_vector_ptr[i + 1]   = Utilities::pow(n_rows, dim);
+            }
+
+          for (unsigned int i = 0; i < n_diag; ++i)
+            new_vector_ptr[i + 1] += new_vector_ptr[i];
+
+          this->inverted_eigenvalues.resize(new_vector_ptr.back());
+
+          // step 3) loop over all unique diagonal entries and invert
+          for (unsigned int i = 0; i < n_diag; ++i)
+            {
+              std::array<Number *, dim> evs;
+
+              for (unsigned int d = 0; d < dim; ++d)
+                evs[d] =
+                  &this
+                     ->eigenvalues[this->vector_ptr[(indices_ev.size() > 0) ?
+                                                      indices_ev[dim * i + d] :
+                                                      (dim * i + d)]];
+
+              const unsigned int mm = new_vector_n_rows_1d[i];
+              if (dim == 2)
+                {
+                  for (unsigned int i1 = 0, c = 0; i1 < mm; ++i1)
+                    for (unsigned int i0 = 0; i0 < mm; ++i0, ++c)
+                      this->inverted_eigenvalues[new_vector_ptr[i] + c] =
+                        Number(1.0) / (evs[1][i1] + evs[0][i0]);
+                }
+              else
+                {
+                  for (unsigned int i2 = 0, c = 0; i2 < mm; ++i2)
+                    for (unsigned int i1 = 0; i1 < mm; ++i1)
+                      for (unsigned int i0 = 0; i0 < mm; ++i0, ++c)
+                        this->inverted_eigenvalues[new_vector_ptr[i] + c] =
+                          Number(1.0) / (evs[2][i2] + evs[1][i1] + evs[0][i0]);
+                }
+            }
+
+          // step 4) clean up
+          std::swap(this->vector_ptr, new_vector_ptr);
+          std::swap(this->vector_n_rows_1d, new_vector_n_rows_1d);
+        }
+
+      this->eigenvalues.clear();
+    }
 }
 
 
@@ -1292,28 +1464,92 @@ TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::
   Number *      dst = dst_in.begin();
   const Number *src = src_in.begin();
 
-  std::array<const Number *, dim> eigenvectors, eigenvalues;
-  unsigned int                    n_rows_1d_non_templated = 0;
-
-  for (unsigned int d = 0; d < dim; ++d)
+  if (this->eigenvalues.empty() == false)
     {
-      const unsigned int translated_index =
-        (indices.size() > 0) ? indices[dim * index + d] : (dim * index + d);
-
-      eigenvectors[d] =
-        this->eigenvectors.data() + matrix_ptr[translated_index];
-      eigenvalues[d] = this->eigenvalues.data() + vector_ptr[translated_index];
-      n_rows_1d_non_templated =
-        vector_ptr[translated_index + 1] - vector_ptr[translated_index];
-    }
+      std::array<const Number *, dim> eigenvectors;
+      std::array<const Number *, dim> eigenvalues;
+      unsigned int                    n_rows_1d_non_templated = 0;
 
-  if (n_rows_1d != -1)
-    internal::TensorProductMatrixSymmetricSum::apply_inverse<
-      n_rows_1d == -1 ? 0 : n_rows_1d>(
-      dst, src, tmp_array, n_rows_1d_non_templated, eigenvectors, eigenvalues);
+      for (unsigned int d = 0; d < dim; ++d)
+        {
+          const unsigned int translated_index =
+            (indices.size() > 0) ? indices[dim * index + d] : (dim * index + d);
+
+          eigenvectors[d] =
+            this->eigenvectors.data() + matrix_ptr[translated_index];
+          eigenvalues[d] =
+            this->eigenvalues.data() + vector_ptr[translated_index];
+          n_rows_1d_non_templated =
+            vector_ptr[translated_index + 1] - vector_ptr[translated_index];
+        }
+
+      if (n_rows_1d != -1)
+        internal::TensorProductMatrixSymmetricSum::apply_inverse<
+          n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
+                                           src,
+                                           tmp_array,
+                                           n_rows_1d_non_templated,
+                                           eigenvectors,
+                                           eigenvalues);
+      else
+        internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
+          dst,
+          src,
+          tmp_array,
+          n_rows_1d_non_templated,
+          eigenvectors,
+          eigenvalues);
+    }
   else
-    internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
-      dst, src, tmp_array, n_rows_1d_non_templated, eigenvectors, eigenvalues);
+    {
+      std::array<const Number *, dim> eigenvectors;
+      const Number *                  inverted_eigenvalues    = nullptr;
+      unsigned int                    n_rows_1d_non_templated = 0;
+
+      for (unsigned int d = 0; d < dim; ++d)
+        {
+          const unsigned int translated_index =
+            (indices.size() > 0) ?
+              indices[((dim == 1) ? 1 : (dim + 1)) * index + d] :
+              (dim * index + d);
+
+          eigenvectors[d] =
+            this->eigenvectors.data() + matrix_ptr[translated_index];
+        }
+
+      {
+        const unsigned int translated_index =
+          ((indices.size() > 0) && (dim != 1)) ?
+            indices[(dim + 1) * index + dim] :
+            index;
+
+        inverted_eigenvalues =
+          this->inverted_eigenvalues.data() + vector_ptr[translated_index];
+        n_rows_1d_non_templated =
+          (dim == 1) ?
+            (vector_ptr[translated_index + 1] - vector_ptr[translated_index]) :
+            vector_n_rows_1d[translated_index];
+      }
+
+      if (n_rows_1d != -1)
+        internal::TensorProductMatrixSymmetricSum::apply_inverse<
+          n_rows_1d == -1 ? 0 : n_rows_1d>(dst,
+                                           src,
+                                           tmp_array,
+                                           n_rows_1d_non_templated,
+                                           eigenvectors,
+                                           {},
+                                           inverted_eigenvalues);
+      else
+        internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
+          dst,
+          src,
+          tmp_array,
+          n_rows_1d_non_templated,
+          eigenvectors,
+          {},
+          inverted_eigenvalues);
+    }
 }
 
 

diff --git a/include/deal.II/lac/tensor_product_matrix.templates.h b/include/deal.II/lac/tensor_product_matrix.templates.h
index c7b026fea049a736ee870c2eebc0996f7ff91706..90a16f2e325cfcb843c2f5cc3f9d578420226c1e 100644 (file)
--- a/include/deal.II/lac/tensor_product_matrix.templates.h
+++ b/include/deal.II/lac/tensor_product_matrix.templates.h
@@ -69,16 +69,34 @@ namespace internal
                          AlignedVector<Number> &                tmp,
                          const unsigned int                     n_rows_1d,
                          const std::array<const Number *, dim> &eigenvectors,
-                         const std::array<const Number *, dim> &eigenvalues)
+                         const std::array<const Number *, dim> &eigenvalues,
+                         const Number *inverted_eigenvalues)
     {
       if (n_rows_1d_templated == n_rows_1d)
-        apply_inverse<n_rows_1d_templated>(
-          dst, src, tmp, n_rows_1d, eigenvectors, eigenvalues);
+        apply_inverse<n_rows_1d_templated>(dst,
+                                           src,
+                                           tmp,
+                                           n_rows_1d,
+                                           eigenvectors,
+                                           eigenvalues,
+                                           inverted_eigenvalues);
       else if (n_rows_1d_templated < FDM_N_ROWS_MAX)
         select_apply_inverse<std::min(n_rows_1d_templated + 1, FDM_N_ROWS_MAX)>(
-          dst, src, tmp, n_rows_1d, eigenvectors, eigenvalues);
+          dst,
+          src,
+          tmp,
+          n_rows_1d,
+          eigenvectors,
+          eigenvalues,
+          inverted_eigenvalues);
       else
-        apply_inverse<0>(dst, src, tmp, n_rows_1d, eigenvectors, eigenvalues);
+        apply_inverse<0>(dst,
+                         src,
+                         tmp,
+                         n_rows_1d,
+                         eigenvectors,
+                         eigenvalues,
+                         inverted_eigenvalues);
     }
   } // namespace TensorProductMatrixSymmetricSum
 } // namespace internal

diff --git a/source/lac/tensor_product_matrix.inst.in b/source/lac/tensor_product_matrix.inst.in
index 6607de8aca7c560481edd542dbd80b96f14ccf29..ec5c1c41cadb2cfceafd69871a41060b4c5a46bc 100644 (file)
--- a/source/lac/tensor_product_matrix.inst.in
+++ b/source/lac/tensor_product_matrix.inst.in
@@ -35,7 +35,8 @@ for (deal_II_dimension : DIMENSIONS;
       const std::array<const deal_II_scalar_vectorized *, deal_II_dimension>
         &eigenvectors,
       const std::array<const deal_II_scalar_vectorized *, deal_II_dimension>
-        &eigenvalues);
+        &                              eigenvalues,
+      const deal_II_scalar_vectorized *inverted_eigenvalues);
   }
 
 for (deal_II_dimension : DIMENSIONS; deal_II_scalar : REAL_SCALARS)
@@ -55,5 +56,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_scalar : REAL_SCALARS)
       AlignedVector<deal_II_scalar> &                              tmp,
       const unsigned int                                           n_rows,
       const std::array<const deal_II_scalar *, deal_II_dimension> &eigenvectors,
-      const std::array<const deal_II_scalar *, deal_II_dimension> &eigenvalues);
+      const std::array<const deal_II_scalar *, deal_II_dimension> &eigenvalues,
+      const deal_II_scalar *inverted_eigenvalues);
   }

diff --git a/tests/lac/tensor_product_matrix_08.cc b/tests/lac/tensor_product_matrix_08.cc
index 508869c3b9a190336211f07422452bafa8305572..bbe162923af4254bd7ebc05b01e7877dcee039cc 100644 (file)
--- a/tests/lac/tensor_product_matrix_08.cc
+++ b/tests/lac/tensor_product_matrix_08.cc
@@ -53,11 +53,15 @@ do_test_mesh(const Mapping<dim> &mapping, const Triangulation<dim> &tria)
   const auto harmonic_patch_extent =
     GridTools::compute_harmonic_patch_extent(mapping, tria, quadrature_face);
 
-  FDM collection_0(typename FDM::AdditionalData(true));
-  FDM collection_1(typename FDM::AdditionalData(false));
+  FDM collection_0(typename FDM::AdditionalData(true, false));
+  FDM collection_1(typename FDM::AdditionalData(false, false));
+  FDM collection_2(typename FDM::AdditionalData(true, true));
+  FDM collection_3(typename FDM::AdditionalData(false, true));
 
   collection_0.reserve(tria.n_active_cells());
   collection_1.reserve(tria.n_active_cells());
+  collection_2.reserve(tria.n_active_cells());
+  collection_3.reserve(tria.n_active_cells());
 
   for (const auto &cell : tria.active_cell_iterators())
     {
@@ -80,10 +84,18 @@ do_test_mesh(const Mapping<dim> &mapping, const Triangulation<dim> &tria)
       collection_1.insert(cell->active_cell_index(),
                           M_and_K.first,
                           M_and_K.second);
+      collection_2.insert(cell->active_cell_index(),
+                          M_and_K.first,
+                          M_and_K.second);
+      collection_3.insert(cell->active_cell_index(),
+                          M_and_K.first,
+                          M_and_K.second);
     }
 
   collection_0.finalize();
   collection_1.finalize();
+  collection_2.finalize();
+  collection_3.finalize();
 
   deallog << "Storage sizes: " << collection_0.storage_size() << " "
           << collection_1.storage_size() << std::endl;
@@ -93,6 +105,8 @@ do_test_mesh(const Mapping<dim> &mapping, const Triangulation<dim> &tria)
   AlignedVector<Number> tmp;
   FullMatrix<Number>    matrix_0(fe.n_dofs_per_cell(), fe.n_dofs_per_cell());
   FullMatrix<Number>    matrix_1(fe.n_dofs_per_cell(), fe.n_dofs_per_cell());
+  FullMatrix<Number>    matrix_2(fe.n_dofs_per_cell(), fe.n_dofs_per_cell());
+  FullMatrix<Number>    matrix_3(fe.n_dofs_per_cell(), fe.n_dofs_per_cell());
 
   for (unsigned int cell = 0; cell < tria.n_active_cells(); ++cell)
     {
@@ -114,14 +128,34 @@ do_test_mesh(const Mapping<dim> &mapping, const Triangulation<dim> &tria)
                                      tmp);
           for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
             matrix_1[j][i] = dst[j];
+
+          collection_2.apply_inverse(cell,
+                                     make_array_view(dst),
+                                     make_array_view(src),
+                                     tmp);
+          for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
+            matrix_2[j][i] = dst[j];
+
+          collection_3.apply_inverse(cell,
+                                     make_array_view(dst),
+                                     make_array_view(src),
+                                     tmp);
+          for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
+            matrix_3[j][i] = dst[j];
         }
 
 
-      FloatingPointComparator<Number> comp(1e-5, false);
+      FloatingPointComparator<Number> comp(1e-5, true);
 
       Assert((comp.compare(matrix_0, matrix_1) ==
               FloatingPointComparator<Number>::ComparisonResult::equal),
              ExcInternalError());
+      Assert((comp.compare(matrix_0, matrix_2) ==
+              FloatingPointComparator<Number>::ComparisonResult::equal),
+             ExcInternalError());
+      Assert((comp.compare(matrix_0, matrix_3) ==
+              FloatingPointComparator<Number>::ComparisonResult::equal),
+             ExcInternalError());
     }
 
   deallog << "OK!" << std::endl;