From: Peter Munch Date: Wed, 21 Sep 2022 18:57:43 +0000 (+0200) Subject: Use FloatingPointComparator in TensorProductMatrixSymmetricSumCollection X-Git-Tag: v9.5.0-rc1~928^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d41588156f388bd18db8e16d962d390a2ccb897e;p=dealii.git Use FloatingPointComparator in TensorProductMatrixSymmetricSumCollection --- diff --git a/include/deal.II/base/floating_point_comparator.h b/include/deal.II/base/floating_point_comparator.h index a32a1bc6fd..2a7f64ca9f 100644 --- a/include/deal.II/base/floating_point_comparator.h +++ b/include/deal.II/base/floating_point_comparator.h @@ -18,9 +18,11 @@ #include +#include #include #include +#include #include DEAL_II_NAMESPACE_OPEN @@ -38,79 +40,97 @@ DEAL_II_NAMESPACE_OPEN * satisfied). This is not a problem in the use cases for this class, but be * careful when using it in other contexts. */ -template +template struct FloatingPointComparator { - using Number = typename dealii::internal::VectorizedArrayTrait< - VectorizedArrayType>::value_type; + using ScalarNumber = + typename dealii::internal::VectorizedArrayTrait::value_type; static constexpr std::size_t width = - dealii::internal::VectorizedArrayTrait::width; + dealii::internal::VectorizedArrayTrait::width; - FloatingPointComparator(const Number scaling); + /** + * Constructor. + */ + FloatingPointComparator( + const ScalarNumber tolerance, + const bool use_absolute_tolerance = true, + const std::bitset &mask = std::bitset().flip()); + + FloatingPointComparator(const FloatingPointComparator &rhs) = default; + + FloatingPointComparator(FloatingPointComparator &&rhs) noexcept = default; /** * Compare two vectors of numbers (not necessarily of the same length), * where vectors of different lengths are first sorted by their length and * then by the entries. */ + template + bool + operator()(const std::vector &v1, const std::vector &v2) const; + + /** + * Compare two arrays. + */ + template bool - operator()(const std::vector &v1, - const std::vector &v2) const; + operator()(const std::array &t1, const std::array &t2) const; /** - * Compare two vectorized arrays (stored as tensors to avoid alignment - * issues). + * Compare two tensors. */ + template bool - operator()(const Tensor<1, width, Number> &t1, - const Tensor<1, width, Number> &t2) const; + operator()(const Tensor &t1, + const Tensor &t2) const; /** - * Compare two rank-1 tensors of vectorized arrays (stored as tensors to - * avoid alignment issues). + * Compare two tables. */ - template + template bool - operator()(const Tensor<1, dim, Tensor<1, width, Number>> &t1, - const Tensor<1, dim, Tensor<1, width, Number>> &t2) const; + operator()(const Table<2, T> &t1, const Table<2, T> &t2) const; /** - * Compare two rank-2 tensors of vectorized arrays (stored as tensors to - * avoid alignment issues). + * Compare two scalar numbers. */ - template bool - operator()(const Tensor<2, dim, Tensor<1, width, Number>> &t1, - const Tensor<2, dim, Tensor<1, width, Number>> &t2) const; + operator()(const ScalarNumber &s1, const ScalarNumber &s2) const; /** - * Compare two arrays of tensors. + * Compare two VectorizedArray instances. */ - template bool - operator()(const std::array, dim + 1> &t1, - const std::array, dim + 1> &t2) const; + operator()(const VectorizedArray &v1, + const VectorizedArray &v2) const; - Number tolerance; +private: + const ScalarNumber tolerance; + const bool use_absolute_tolerance; + const std::bitset mask; }; /* ------------------------------------------------------------------ */ -template -FloatingPointComparator::FloatingPointComparator( - const Number scaling) - : tolerance(scaling * std::numeric_limits::epsilon() * 1024.) +template +FloatingPointComparator::FloatingPointComparator( + const ScalarNumber tolerance, + const bool use_absolute_tolerance, + const std::bitset &mask) + : tolerance(tolerance) + , use_absolute_tolerance(use_absolute_tolerance) + , mask(mask) {} -template +template +template bool -FloatingPointComparator::operator()( - const std::vector &v1, - const std::vector &v2) const +FloatingPointComparator::operator()(const std::vector &v1, + const std::vector &v2) const { const unsigned int s1 = v1.size(), s2 = v2.size(); if (s1 < s2) @@ -119,82 +139,105 @@ FloatingPointComparator::operator()( return false; else for (unsigned int i = 0; i < s1; ++i) - if (v1[i] < v2[i] - tolerance) + if (this->operator()(v1[i], v2[i])) return true; - else if (v1[i] > v2[i] + tolerance) + else if (this->operator()(v2[i], v1[i])) return false; return false; } -template +template +template bool -FloatingPointComparator::operator()( - const Tensor<1, width, Number> &t1, - const Tensor<1, width, Number> &t2) const +FloatingPointComparator::operator()(const std::array &t1, + const std::array &t2) const { - for (unsigned int k = 0; k < width; ++k) - if (t1[k] < t2[k] - tolerance) + for (unsigned int i = 0; i < t1.size(); ++i) + if (this->operator()(t1[i], t2[i])) return true; - else if (t1[k] > t2[k] + tolerance) + else if (this->operator()(t2[i], t1[i])) return false; return false; } -template -template +template +template bool -FloatingPointComparator::operator()( - const Tensor<1, dim, Tensor<1, width, Number>> &t1, - const Tensor<1, dim, Tensor<1, width, Number>> &t2) const +FloatingPointComparator::operator()( + const Tensor &t1, + const Tensor &t2) const { - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int k = 0; k < width; ++k) - if (t1[d][k] < t2[d][k] - tolerance) - return true; - else if (t1[d][k] > t2[d][k] + tolerance) - return false; + for (unsigned int k = 0; k < dim; ++k) + if (this->operator()(t1[k], t2[k])) + return true; + else if (this->operator()(t2[k], t1[k])) + return false; return false; } -template -template +template +template bool -FloatingPointComparator::operator()( - const Tensor<2, dim, Tensor<1, width, Number>> &t1, - const Tensor<2, dim, Tensor<1, width, Number>> &t2) const +FloatingPointComparator::operator()(const Table<2, T> &t1, + const Table<2, T> &t2) const { - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int e = 0; e < dim; ++e) - for (unsigned int k = 0; k < width; ++k) - if (t1[d][e][k] < t2[d][e][k] - tolerance) - return true; - else if (t1[d][e][k] > t2[d][e][k] + tolerance) - return false; + AssertDimension(t1.size(0), t2.size(0)); + AssertDimension(t1.size(1), t2.size(1)); + + for (unsigned int i = 0; i < t1.size(0); ++i) + for (unsigned int j = 0; j < t1.size(1); ++j) + if (this->operator()(t1[i][j], t2[i][j])) + return true; + else if (this->operator()(t2[i][j], t1[i][j])) + return false; return false; } +template +bool +FloatingPointComparator::operator()(const ScalarNumber &s1, + const ScalarNumber &s2) const +{ + if (mask[0]) + { + const ScalarNumber tolerance = use_absolute_tolerance ? + this->tolerance : + (std::abs(s1 + s2) * this->tolerance); + if ((s1 < s2 - tolerance)) + return true; + else + return false; + } + + return false; +} -template -template +template bool -FloatingPointComparator::operator()( - const std::array, dim + 1> &t1, - const std::array, dim + 1> &t2) const +FloatingPointComparator::operator()( + const VectorizedArray &v1, + const VectorizedArray &v2) const { - for (unsigned int i = 0; i < t1.size(); ++i) - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int e = 0; e < dim; ++e) - if (t1[i][d][e] < t2[i][d][e] - tolerance) + for (unsigned int v = 0; v < width; ++v) + if (mask[v]) + { + const ScalarNumber tolerance = + use_absolute_tolerance ? this->tolerance : + (std::abs(v1[v] + v2[v]) * this->tolerance); + + if (v1[v] < v2[v] - tolerance) return true; - else if (t1[i][d][e] > t2[i][d][e] + tolerance) + if (v1[v] > v2[v] + tolerance) return false; + } + return false; } diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h index 165789bf18..0c167ea0c9 100644 --- a/include/deal.II/lac/tensor_product_matrix.h +++ b/include/deal.II/lac/tensor_product_matrix.h @@ -20,6 +20,7 @@ #include #include +#include #include #include @@ -281,6 +282,14 @@ namespace internal struct MatrixPairComparator { using MatrixPairType = std::pair, Table<2, Number>>; + using VectorizedArrayTrait = + dealii::internal::VectorizedArrayTrait; + using ScalarNumber = typename VectorizedArrayTrait::value_type; + static constexpr std::size_t width = VectorizedArrayTrait::width; + + MatrixPairComparator() + : eps(std::sqrt(std::numeric_limits::epsilon())) + {} bool operator()(const MatrixPairType &left, const MatrixPairType &right) const @@ -290,64 +299,42 @@ namespace internal const auto &M_1 = right.first; const auto &K_1 = right.second; - const unsigned int n_lanes = Number::size(); - - std::array mask; + std::bitset mask; - using NumberScalar = typename Number::value_type; + using ScalarNumber = typename Number::value_type; - for (unsigned int v = 0; v < n_lanes; ++v) + for (unsigned int v = 0; v < width; ++v) { - NumberScalar a = 0.0; - NumberScalar b = 0.0; + ScalarNumber a = 0.0; + ScalarNumber 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]); + a += std::abs(VectorizedArrayTrait::get(M_0[i][j], v)); + a += std::abs(VectorizedArrayTrait::get(K_0[i][j], v)); + b += std::abs(VectorizedArrayTrait::get(M_1[i][j], v)); + b += std::abs(VectorizedArrayTrait::get(K_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; + const FloatingPointComparator comparator( + eps, false /*use relativ torlerance*/, mask); + + if (comparator(M_0, M_1)) + return true; + else if (comparator(M_1, M_0)) + return false; + else if (comparator(K_0, K_1)) + return true; + else + return false; } + + private: + const ScalarNumber eps; }; } // namespace TensorProductMatrixSymmetricSum } // namespace internal @@ -1231,7 +1218,9 @@ TensorProductMatrixSymmetricSumCollection:: MemoryConsumption::memory_consumption(mass_matrices) + MemoryConsumption::memory_consumption(derivative_matrices) + MemoryConsumption::memory_consumption(eigenvectors) + - MemoryConsumption::memory_consumption(eigenvalues); + MemoryConsumption::memory_consumption(eigenvalues) + + MemoryConsumption::memory_consumption(matrix_ptr) + + MemoryConsumption::memory_consumption(vector_ptr); } @@ -1241,7 +1230,10 @@ std::size_t TensorProductMatrixSymmetricSumCollection:: storage_size() const { - return mass_matrices.size(); + if (matrix_ptr.size() == 0) + return 0; // if not initialized + + return matrix_ptr.size() - 1; } diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h index 47e964a92e..7acb48b8b9 100644 --- a/include/deal.II/matrix_free/dof_info.templates.h +++ b/include/deal.II/matrix_free/dof_info.templates.h @@ -41,7 +41,8 @@ namespace internal { template ConstraintValues::ConstraintValues() - : constraints(FloatingPointComparator>(1.)) + : constraints(FloatingPointComparator>( + 1. * std::numeric_limits::epsilon() * 1024.)) {} diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index 75faac421c..bb584128e9 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -348,7 +348,8 @@ namespace internal struct CompressedCellData { CompressedCellData(const double expected_size) - : data(FloatingPointComparator(expected_size)) + : data(FloatingPointComparator( + expected_size * std::numeric_limits::epsilon() * 1024.)) {} std::map>, @@ -1078,12 +1079,12 @@ namespace internal // first quadrature point on the cell - we use a relatively coarse // tolerance to account for some inaccuracies in the manifold // evaluation - const FloatingPointComparator> comparator( - 1e4 * jacobian_size); std::map, dim + 1>, unsigned int, FloatingPointComparator>> - compressed_jacobians(comparator); + compressed_jacobians(FloatingPointComparator>( + 1e4 * jacobian_size * std::numeric_limits::epsilon() * + 1024.)); unsigned int n_data_buckets = 0; for (unsigned int cell = 0; cell < jacobians_on_stencil.size(); ++cell) @@ -1505,8 +1506,9 @@ namespace internal // CompressedCellData) and add another factor of 512 to account for // some roundoff effects. CompressedFaceData(const Number jacobian_size) - : data(FloatingPointComparator(512. / - jacobian_size)) + : data(FloatingPointComparator( + 512. / jacobian_size * std::numeric_limits::epsilon() * + 1024.)) , jacobian_size(jacobian_size) {}