]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add TensorProductMatrixSymmetricSumCache 14289/head
authorPeter Munch <peterrmuench@gmail.com>
Mon, 19 Sep 2022 14:36:45 +0000 (16:36 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Tue, 20 Sep 2022 13:47:29 +0000 (15:47 +0200)
include/deal.II/lac/tensor_product_matrix.h

index 9386deac7c2b4cef34319911a0567567393e7961..a54d900d1673d0536be01efb6ba0c5f021960f39 100644 (file)
@@ -272,6 +272,210 @@ private:
 };
 
 
+
+namespace internal
+{
+  namespace TensorProductMatrixSymmetricSum
+  {
+    template <typename Number>
+    struct MatrixPairComparator
+    {
+      using MatrixPairType = std::pair<Table<2, Number>, Table<2, Number>>;
+
+      bool
+      operator()(const MatrixPairType &left, const MatrixPairType &right) const
+      {
+        const auto &M_0 = left.first;
+        const auto &K_0 = left.second;
+        const auto &M_1 = right.first;
+        const auto &K_1 = right.second;
+
+        const unsigned int n_lanes = Number::size();
+
+        std::array<bool, n_lanes> mask;
+
+        using NumberScalar = typename Number::value_type;
+
+        for (unsigned int v = 0; v < n_lanes; ++v)
+          {
+            NumberScalar a = 0.0;
+            NumberScalar b = 0.0;
+
+            for (unsigned int i = 0; i < M_0.size(0); ++i)
+              for (unsigned int j = 0; j < M_0.size(1); ++j)
+                {
+                  a += std::abs(M_0[i][j][v]);
+                  b += std::abs(M_1[i][j][v]);
+                }
+
+            mask[v] = (a != 0.0) && (b != 0.0);
+          }
+
+        const auto eps = std::pow<NumberScalar>(
+          static_cast<NumberScalar>(10.0),
+          static_cast<NumberScalar>(
+            std::log10(std::numeric_limits<NumberScalar>::epsilon()) / 2.0));
+
+        const auto less = [eps](const auto a, const auto b) -> bool {
+          return (b - a) > std::max(eps, std::abs(a + b) * eps);
+        };
+
+        const auto greater = [eps](const auto a, const auto b) -> bool {
+          return (a - b) > std::max(eps, std::abs(a + b) * eps);
+        };
+
+        for (unsigned int v = 0; v < n_lanes; ++v)
+          if (mask[v])
+            for (unsigned int i = 0; i < M_0.size(0); ++i)
+              for (unsigned int j = 0; j < M_0.size(1); ++j)
+                {
+                  if (less(M_0[i][j][v], M_1[i][j][v]))
+                    return true;
+                  else if (greater(M_0[i][j][v], M_1[i][j][v]))
+                    return false;
+                }
+
+        for (unsigned int v = 0; v < n_lanes; ++v)
+          if (mask[v])
+            for (unsigned int i = 0; i < K_0.size(0); ++i)
+              for (unsigned int j = 0; j < K_0.size(1); ++j)
+                {
+                  if (less(K_0[i][j][v], K_1[i][j][v]))
+                    return true;
+                  else if (greater(K_0[i][j][v], K_1[i][j][v]))
+                    return false;
+                }
+
+        return false;
+      }
+    };
+  } // namespace TensorProductMatrixSymmetricSum
+} // namespace internal
+
+
+
+/**
+ * A class similar to TensorProductMatrixSymmetricSum.
+ *
+ * The class TensorProductMatrixSymmetricSum stores a
+ * 1D mass matrix, 1D stiffness matrix, eigenvalues and eigenvectors
+ * for each direction. If one uses one TensorProductMatrixSymmetricSum
+ * instance for, e.g., each cell, these quantities are stored
+ * for each cell. There is no possibility to reuse quantities between
+ * TensorProductMatrixSymmetricSum instances even if the values of the
+ * internal data structures might be the same. This class targets the case
+ * of many TensorProductMatrixSymmetricSum instances, where some of them might
+ * possibly share the underlying 1D matrices and hence re-use the same data.
+ *
+ * This class is flexible and allows to interpret the parameter
+ * @p index arbitrarily. In the case of an element-centric patch
+ * smoother, the index might correspond to the cell index and,
+ * in the case of a vertex patch smoother, the index might
+ * correspond to the vertex index defining the vertex patch. If
+ * @p n_rows_1d is set to -1, the sizes of the mass matrices and
+ * the stiffness matrices can differ between cells, which might be
+ * useful if the class is used in the context of hp-refinement to
+ * construct a smoother.
+ */
+template <int dim, typename Number, int n_rows_1d = -1>
+class TensorProductMatrixSymmetricSumCollection
+{
+  using MatrixPairType = std::pair<Table<2, Number>, Table<2, Number>>;
+
+public:
+  /**
+   * Allocate memory. The parameter @p specifies the maximum value
+   * of the index used in invert() and apply_inverse().
+   */
+  void
+  reserve(const unsigned int size);
+
+  /**
+   * For a given @p index, attach the mass matrices @p Ms and
+   * stiffness matrices @p Ks to the stored data, looking out for
+   * compression possibilities.
+   */
+  void
+  insert(const unsigned int                       index,
+         const std::array<Table<2, Number>, dim> &Ms,
+         const std::array<Table<2, Number>, dim> &Ks);
+
+  /**
+   * Finalize setup. This function computes, e.g., the
+   * eigenvalues and the eigenvectors.
+   */
+  void
+  finalize();
+
+  /**
+   * Apply the inverse matrix for a given @p index.
+   */
+  void
+  apply_inverse(const unsigned int             index,
+                const ArrayView<Number> &      dst_in,
+                const ArrayView<const Number> &src_in,
+                AlignedVector<Number> &        tmp_array) const;
+
+  /**
+   * Return the memory consumption of this class in bytes.
+   */
+  std::size_t
+  memory_consumption() const;
+
+  /**
+   * Return the number of 1D matrices of each type stored internally.
+   * In the case that no compression could be performed, its value
+   * is the parameter passed to the function reserve() times the
+   * number of dimension. If compression could be performed, the
+   * value returned is less. In the optimal case (uniform Cartesian
+   * mesh), the value is one.
+   */
+  std::size_t
+  storage_size() const;
+
+private:
+  /**
+   * Container used during setup to determine the unique 1D
+   * matrices. The memory is freed during finalize().
+   */
+  std::map<
+    MatrixPairType,
+    unsigned int,
+    internal::TensorProductMatrixSymmetricSum::MatrixPairComparator<Number>>
+    cache;
+
+  /**
+   * Map from index to the storage position within vector_mass_matrix,
+   * vector_derivative_matrix, vector_eigenvectors, and
+   * vector_eigenvalues. If compression was not successful, this
+   * vector is empty, since the vectors can be access directly
+   * with the given index.
+   */
+  std::vector<unsigned int> indices;
+
+  /**
+   * Vector of 1D mass matrices.
+   */
+  std::vector<Table<2, Number>> vector_mass_matrix;
+
+  /**
+   * Vector of 1D derivative matrices.
+   */
+  std::vector<Table<2, Number>> vector_derivative_matrix;
+
+  /**
+   * Vector of eigenvectors.
+   */
+  std::vector<Table<2, Number>> vector_eigenvectors;
+
+  /**
+   * Vector of eigenvalues.
+   */
+  std::vector<AlignedVector<Number>> vector_eigenvalues;
+};
+
+
+
 /*----------------------- Inline functions ----------------------------------*/
 
 #ifndef DOXYGEN
@@ -838,6 +1042,154 @@ TensorProductMatrixSymmetricSum<dim, Number, n_rows_1d>::reinit(
 
 
 
+template <int dim, typename Number, int n_rows_1d>
+void
+TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::reserve(
+  const unsigned int size)
+{
+  indices.assign(size * dim, numbers::invalid_unsigned_int);
+}
+
+
+
+template <int dim, typename Number, int n_rows_1d>
+void
+TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::insert(
+  const unsigned int                       index,
+  const std::array<Table<2, Number>, dim> &Ms,
+  const std::array<Table<2, Number>, dim> &Ks)
+{
+  for (unsigned int d = 0; d < dim; ++d)
+    {
+      const MatrixPairType matrix(Ms[d], Ks[d]);
+
+      const auto ptr = cache.find(matrix);
+
+      if (ptr != cache.end())
+        indices[index * dim + d] = ptr->second;
+      else
+        {
+          indices[index * dim + d] = cache.size();
+          cache[matrix]            = cache.size();
+        }
+    }
+}
+
+
+
+template <int dim, typename Number, int n_rows_1d>
+void
+TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::finalize()
+{
+  vector_mass_matrix.resize(cache.size());
+  vector_derivative_matrix.resize(cache.size());
+  vector_eigenvectors.resize(cache.size());
+  vector_eigenvalues.resize(cache.size());
+
+  const auto store = [&](const unsigned int    index,
+                         const MatrixPairType &M_and_K) {
+    std::array<Table<2, Number>, 1> mass_matrices;
+    mass_matrices[0] = M_and_K.first;
+
+    std::array<Table<2, Number>, 1> derivative_matrices;
+    derivative_matrices[0] = M_and_K.second;
+
+    std::array<Table<2, Number>, 1>      eigenvectors;
+    std::array<AlignedVector<Number>, 1> eigenvalues;
+
+    internal::TensorProductMatrixSymmetricSum::setup(mass_matrices,
+                                                     derivative_matrices,
+                                                     eigenvectors,
+                                                     eigenvalues);
+
+    vector_mass_matrix[index]       = M_and_K.first;
+    vector_derivative_matrix[index] = M_and_K.second;
+    vector_eigenvectors[index]      = eigenvectors[0];
+    vector_eigenvalues[index]       = eigenvalues[0];
+  };
+
+  if (cache.size() == indices.size())
+    {
+      std::map<unsigned int, MatrixPairType> inverted_cache;
+
+      for (const auto &i : cache)
+        inverted_cache[i.second] = i.first;
+
+      for (unsigned int i = 0; i < indices.size(); ++i)
+        store(i, inverted_cache[indices[i]]);
+
+      indices.clear();
+    }
+  else
+    {
+      for (const auto &i : cache)
+        store(i.second, i.first);
+    }
+
+  cache.clear();
+}
+
+
+
+template <int dim, typename Number, int n_rows_1d>
+void
+TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::
+  apply_inverse(const unsigned int             index,
+                const ArrayView<Number> &      dst_in,
+                const ArrayView<const Number> &src_in,
+                AlignedVector<Number> &        tmp_array) const
+{
+  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)
+    {
+      const unsigned int translated_index =
+        (indices.size() > 0) ? indices[dim * index + d] : (dim * index + d);
+
+      eigenvectors[d]         = &vector_eigenvectors[translated_index](0, 0);
+      eigenvalues[d]          = vector_eigenvalues[translated_index].data();
+      n_rows_1d_non_templated = vector_eigenvalues[translated_index].size();
+    }
+
+  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);
+}
+
+
+
+template <int dim, typename Number, int n_rows_1d>
+std::size_t
+TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::
+  memory_consumption() const
+{
+  return MemoryConsumption::memory_consumption(indices) +
+         MemoryConsumption::memory_consumption(vector_mass_matrix) +
+         MemoryConsumption::memory_consumption(vector_derivative_matrix) +
+         MemoryConsumption::memory_consumption(vector_eigenvectors) +
+         MemoryConsumption::memory_consumption(vector_eigenvalues);
+}
+
+
+
+template <int dim, typename Number, int n_rows_1d>
+std::size_t
+TensorProductMatrixSymmetricSumCollection<dim, Number, n_rows_1d>::
+  storage_size() const
+{
+  return vector_mass_matrix.size();
+}
+
+
+
 #endif
 
 DEAL_II_NAMESPACE_CLOSE

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.