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 a32a1bc6fdbbd2d1f745873264f46420ad4a833a..2a7f64ca9f7a4248404c016b3eb76eec20ae7778 100644 (file)
--- a/include/deal.II/base/floating_point_comparator.h
+++ b/include/deal.II/base/floating_point_comparator.h
@@ -18,9 +18,11 @@
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/base/table.h>
 #include <deal.II/base/tensor.h>
 #include <deal.II/base/vectorization.h>
 
+#include <bitset>
 #include <vector>
 
 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 <typename VectorizedArrayType>
+template <typename Number>
 struct FloatingPointComparator
 {
-  using Number = typename dealii::internal::VectorizedArrayTrait<
-    VectorizedArrayType>::value_type;
+  using ScalarNumber =
+    typename dealii::internal::VectorizedArrayTrait<Number>::value_type;
   static constexpr std::size_t width =
-    dealii::internal::VectorizedArrayTrait<VectorizedArrayType>::width;
+    dealii::internal::VectorizedArrayTrait<Number>::width;
 
-  FloatingPointComparator(const Number scaling);
+  /**
+   * Constructor.
+   */
+  FloatingPointComparator(
+    const ScalarNumber        tolerance,
+    const bool                use_absolute_tolerance = true,
+    const std::bitset<width> &mask = std::bitset<width>().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 <typename T>
+  bool
+  operator()(const std::vector<T> &v1, const std::vector<T> &v2) const;
+
+  /**
+   * Compare two arrays.
+   */
+  template <std::size_t dim, typename T>
   bool
-  operator()(const std::vector<Number> &v1,
-             const std::vector<Number> &v2) const;
+  operator()(const std::array<T, dim> &t1, const std::array<T, dim> &t2) const;
 
   /**
-   * Compare two vectorized arrays (stored as tensors to avoid alignment
-   * issues).
+   * Compare two tensors.
    */
+  template <int rank, int dim, typename T>
   bool
-  operator()(const Tensor<1, width, Number> &t1,
-             const Tensor<1, width, Number> &t2) const;
+  operator()(const Tensor<rank, dim, T> &t1,
+             const Tensor<rank, dim, T> &t2) const;
 
   /**
-   * Compare two rank-1 tensors of vectorized arrays (stored as tensors to
-   * avoid alignment issues).
+   * Compare two tables.
    */
-  template <int dim>
+  template <typename T>
   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 <int dim>
   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 <int dim>
   bool
-  operator()(const std::array<Tensor<2, dim, Number>, dim + 1> &t1,
-             const std::array<Tensor<2, dim, Number>, dim + 1> &t2) const;
+  operator()(const VectorizedArray<ScalarNumber, width> &v1,
+             const VectorizedArray<ScalarNumber, width> &v2) const;
 
-  Number tolerance;
+private:
+  const ScalarNumber       tolerance;
+  const bool               use_absolute_tolerance;
+  const std::bitset<width> mask;
 };
 
 
 /* ------------------------------------------------------------------ */
 
 
-template <typename VectorizedArrayType>
-FloatingPointComparator<VectorizedArrayType>::FloatingPointComparator(
-  const Number scaling)
-  : tolerance(scaling * std::numeric_limits<double>::epsilon() * 1024.)
+template <typename Number>
+FloatingPointComparator<Number>::FloatingPointComparator(
+  const ScalarNumber        tolerance,
+  const bool                use_absolute_tolerance,
+  const std::bitset<width> &mask)
+  : tolerance(tolerance)
+  , use_absolute_tolerance(use_absolute_tolerance)
+  , mask(mask)
 {}
 
 
 
-template <typename VectorizedArrayType>
+template <typename Number>
+template <typename T>
 bool
-FloatingPointComparator<VectorizedArrayType>::operator()(
-  const std::vector<Number> &v1,
-  const std::vector<Number> &v2) const
+FloatingPointComparator<Number>::operator()(const std::vector<T> &v1,
+                                            const std::vector<T> &v2) const
 {
   const unsigned int s1 = v1.size(), s2 = v2.size();
   if (s1 < s2)
@@ -119,82 +139,105 @@ FloatingPointComparator<VectorizedArrayType>::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 <typename VectorizedArrayType>
+template <typename Number>
+template <std::size_t dim, typename T>
 bool
-FloatingPointComparator<VectorizedArrayType>::operator()(
-  const Tensor<1, width, Number> &t1,
-  const Tensor<1, width, Number> &t2) const
+FloatingPointComparator<Number>::operator()(const std::array<T, dim> &t1,
+                                            const std::array<T, dim> &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 <typename VectorizedArrayType>
-template <int dim>
+template <typename Number>
+template <int rank, int dim, typename T>
 bool
-FloatingPointComparator<VectorizedArrayType>::operator()(
-  const Tensor<1, dim, Tensor<1, width, Number>> &t1,
-  const Tensor<1, dim, Tensor<1, width, Number>> &t2) const
+FloatingPointComparator<Number>::operator()(
+  const Tensor<rank, dim, T> &t1,
+  const Tensor<rank, dim, T> &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 <typename VectorizedArrayType>
-template <int dim>
+template <typename Number>
+template <typename T>
 bool
-FloatingPointComparator<VectorizedArrayType>::operator()(
-  const Tensor<2, dim, Tensor<1, width, Number>> &t1,
-  const Tensor<2, dim, Tensor<1, width, Number>> &t2) const
+FloatingPointComparator<Number>::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 <typename Number>
+bool
+FloatingPointComparator<Number>::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 <typename VectorizedArrayType>
-template <int dim>
+template <typename Number>
 bool
-FloatingPointComparator<VectorizedArrayType>::operator()(
-  const std::array<Tensor<2, dim, Number>, dim + 1> &t1,
-  const std::array<Tensor<2, dim, Number>, dim + 1> &t2) const
+FloatingPointComparator<Number>::operator()(
+  const VectorizedArray<ScalarNumber, width> &v1,
+  const VectorizedArray<ScalarNumber, width> &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 165789bf18b0414a4ad813e3654e6fbdad6076ea..0c167ea0c90368638f128504204d6769b8ab94da 100644 (file)
--- a/include/deal.II/lac/tensor_product_matrix.h
+++ b/include/deal.II/lac/tensor_product_matrix.h
@@ -20,6 +20,7 @@
 #include <deal.II/base/config.h>
 
 #include <deal.II/base/array_view.h>
+#include <deal.II/base/floating_point_comparator.h>
 #include <deal.II/base/mutex.h>
 #include <deal.II/base/vectorization.h>
 
@@ -281,6 +282,14 @@ namespace internal
     struct MatrixPairComparator
     {
       using MatrixPairType = std::pair<Table<2, Number>, Table<2, Number>>;
+      using VectorizedArrayTrait =
+        dealii::internal::VectorizedArrayTrait<Number>;
+      using ScalarNumber = typename VectorizedArrayTrait::value_type;
+      static constexpr std::size_t width = VectorizedArrayTrait::width;
+
+      MatrixPairComparator()
+        : eps(std::sqrt(std::numeric_limits<ScalarNumber>::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<bool, n_lanes> mask;
+        std::bitset<width> 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<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;
+        const FloatingPointComparator<Number> 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<dim, Number, n_rows_1d>::
          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<dim, Number, n_rows_1d>::
   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 47e964a92ed88acefc10226291ba15a919cfb599..7acb48b8b9e61dd591da89e1e06269875a3356d7 100644 (file)
--- 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 <typename Number>
     ConstraintValues<Number>::ConstraintValues()
-      : constraints(FloatingPointComparator<VectorizedArray<Number>>(1.))
+      : constraints(FloatingPointComparator<VectorizedArray<Number>>(
+          1. * std::numeric_limits<double>::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 75faac421c014603442c9c69034dbf8ab819ea23..bb584128e9855cd5f467c9d3a885403ec479fc5b 100644 (file)
--- 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<VectorizedArrayType>(expected_size))
+          : data(FloatingPointComparator<VectorizedArrayType>(
+              expected_size * std::numeric_limits<double>::epsilon() * 1024.))
         {}
 
         std::map<Tensor<2, dim, Tensor<1, VectorizedArrayType::size(), Number>>,
@@ -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<VectorizedArray<double>> comparator(
-          1e4 * jacobian_size);
         std::map<std::array<Tensor<2, dim>, dim + 1>,
                  unsigned int,
                  FloatingPointComparator<VectorizedArray<double>>>
-          compressed_jacobians(comparator);
+          compressed_jacobians(FloatingPointComparator<VectorizedArray<double>>(
+            1e4 * jacobian_size * std::numeric_limits<double>::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<VectorizedArrayType>(512. /
-                                                              jacobian_size))
+          : data(FloatingPointComparator<VectorizedArrayType>(
+              512. / jacobian_size * std::numeric_limits<double>::epsilon() *
+              1024.))
           , jacobian_size(jacobian_size)
         {}