From: Peter Munch Date: Mon, 19 Sep 2022 14:36:45 +0000 (+0200) Subject: Add TensorProductMatrixSymmetricSumCache X-Git-Tag: v9.5.0-rc1~946^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a4cb107f8dd8aaafb228a5d0fd77d72587dff3ca;p=dealii.git Add TensorProductMatrixSymmetricSumCache --- diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h index 9386deac7c..a54d900d16 100644 --- a/include/deal.II/lac/tensor_product_matrix.h +++ b/include/deal.II/lac/tensor_product_matrix.h @@ -272,6 +272,210 @@ private: }; + +namespace internal +{ + namespace TensorProductMatrixSymmetricSum + { + template + struct MatrixPairComparator + { + using MatrixPairType = std::pair, 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 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( + static_cast(10.0), + static_cast( + std::log10(std::numeric_limits::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 +class TensorProductMatrixSymmetricSumCollection +{ + using MatrixPairType = std::pair, 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, dim> &Ms, + const std::array, 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 & dst_in, + const ArrayView &src_in, + AlignedVector & 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> + 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 indices; + + /** + * Vector of 1D mass matrices. + */ + std::vector> vector_mass_matrix; + + /** + * Vector of 1D derivative matrices. + */ + std::vector> vector_derivative_matrix; + + /** + * Vector of eigenvectors. + */ + std::vector> vector_eigenvectors; + + /** + * Vector of eigenvalues. + */ + std::vector> vector_eigenvalues; +}; + + + /*----------------------- Inline functions ----------------------------------*/ #ifndef DOXYGEN @@ -838,6 +1042,154 @@ TensorProductMatrixSymmetricSum::reinit( +template +void +TensorProductMatrixSymmetricSumCollection::reserve( + const unsigned int size) +{ + indices.assign(size * dim, numbers::invalid_unsigned_int); +} + + + +template +void +TensorProductMatrixSymmetricSumCollection::insert( + const unsigned int index, + const std::array, dim> &Ms, + const std::array, 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 +void +TensorProductMatrixSymmetricSumCollection::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, 1> mass_matrices; + mass_matrices[0] = M_and_K.first; + + std::array, 1> derivative_matrices; + derivative_matrices[0] = M_and_K.second; + + std::array, 1> eigenvectors; + std::array, 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 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 +void +TensorProductMatrixSymmetricSumCollection:: + apply_inverse(const unsigned int index, + const ArrayView & dst_in, + const ArrayView &src_in, + AlignedVector & tmp_array) const +{ + Number * dst = dst_in.begin(); + const Number *src = src_in.begin(); + + std::array 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 +std::size_t +TensorProductMatrixSymmetricSumCollection:: + 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 +std::size_t +TensorProductMatrixSymmetricSumCollection:: + storage_size() const +{ + return vector_mass_matrix.size(); +} + + + #endif DEAL_II_NAMESPACE_CLOSE