floating_point_comparator.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------
 
#ifndef dealii_base_floating_point_copmerator_h
#define dealii_base_floating_point_copmerator_h
 
#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
 
template <typename Number>
struct FloatingPointComparator
{
  using ScalarNumber =
    typename ::internal::VectorizedArrayTrait<Number>::value_type;
  static constexpr std::size_t width =
    ::internal::VectorizedArrayTrait<Number>::width();
 
  enum class ComparisonResult
  {
    less,
    equal,
    greater
  };
 
  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;
 
  template <typename T>
  bool
  operator()(const T &object1, const T &object2) const;
 
  template <typename T>
  ComparisonResult
  compare(const std::vector<T> &v1, const std::vector<T> &v2) const;
 
  template <std::size_t dim, typename T>
  ComparisonResult
  compare(const std::array<T, dim> &t1, const std::array<T, dim> &t2) const;
 
  template <int rank, int dim, typename T>
  ComparisonResult
  compare(const Tensor<rank, dim, T> &t1, const Tensor<rank, dim, T> &t2) const;
 
  template <typename T>
  ComparisonResult
  compare(const Table<2, T> &t1, const Table<2, T> &t2) const;
 
  ComparisonResult
  compare(const ScalarNumber s1, const ScalarNumber s2) const;
 
  ComparisonResult
  compare(const VectorizedArray<ScalarNumber, width> &v1,
          const VectorizedArray<ScalarNumber, width> &v2) const;
 
private:
  const ScalarNumber       tolerance;
  const bool               use_absolute_tolerance;
  const std::bitset<width> mask;
};
 
 
/* ------------------------------------------------------------------ */
 
 
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)
{
  Assert(mask.count() > 0, ExcMessage("No component selected"));
}
 
 
 
template <typename Number>
template <typename T>
bool
FloatingPointComparator<Number>::operator()(const T &object1,
                                            const T &object2) const
{
  return compare(object1, object2) == ComparisonResult::less;
}
 
 
 
template <typename Number>
template <typename T>
typename FloatingPointComparator<Number>::ComparisonResult
FloatingPointComparator<Number>::compare(const std::vector<T> &v1,
                                         const std::vector<T> &v2) const
{
  const unsigned int s1 = v1.size(), s2 = v2.size();
  if (s1 < s2)
    return ComparisonResult::less;
  else if (s1 > s2)
    return ComparisonResult::greater;
  else
    for (unsigned int i = 0; i < s1; ++i)
      {
        const ComparisonResult result = compare(v1[i], v2[i]);
        if (result != ComparisonResult::equal)
          return result;
      }
  return ComparisonResult::equal;
}
 
 
 
template <typename Number>
template <std::size_t dim, typename T>
typename FloatingPointComparator<Number>::ComparisonResult
FloatingPointComparator<Number>::compare(const std::array<T, dim> &t1,
                                         const std::array<T, dim> &t2) const
{
  for (unsigned int i = 0; i < t1.size(); ++i)
    {
      const ComparisonResult result = compare(t1[i], t2[i]);
      if (result != ComparisonResult::equal)
        return result;
    }
  return ComparisonResult::equal;
}
 
 
 
template <typename Number>
template <int rank, int dim, typename T>
typename FloatingPointComparator<Number>::ComparisonResult
FloatingPointComparator<Number>::compare(const Tensor<rank, dim, T> &t1,
                                         const Tensor<rank, dim, T> &t2) const
{
  for (unsigned int i = 0; i < dim; ++i)
    {
      const ComparisonResult result = compare(t1[i], t2[i]);
      if (result != ComparisonResult::equal)
        return result;
    }
  return ComparisonResult::equal;
}
 
 
 
template <typename Number>
template <typename T>
typename FloatingPointComparator<Number>::ComparisonResult
FloatingPointComparator<Number>::compare(const Table<2, T> &t1,
                                         const Table<2, T> &t2) const
{
  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)
      {
        const ComparisonResult result = compare(t1[i][j], t2[i][j]);
        if (result != ComparisonResult::equal)
          return result;
      }
  return ComparisonResult::equal;
}
 
 
 
template <typename Number>
typename FloatingPointComparator<Number>::ComparisonResult
FloatingPointComparator<Number>::compare(const ScalarNumber s1,
                                         const ScalarNumber s2) const
{
  if (width == 1 || mask[0])
    {
      const ScalarNumber tolerance =
        use_absolute_tolerance ?
          this->tolerance :
          ((std::abs(s1) + std::abs(s2)) * this->tolerance);
 
      if (s1 < s2 - tolerance)
        return ComparisonResult::less;
      else if (s1 > s2 + tolerance)
        return ComparisonResult::greater;
    }
 
  return ComparisonResult::equal;
}
 
template <typename Number>
typename FloatingPointComparator<Number>::ComparisonResult
FloatingPointComparator<Number>::compare(
  const VectorizedArray<ScalarNumber, width> &v1,
  const VectorizedArray<ScalarNumber, width> &v2) const
{
  for (unsigned int i = 0; i < width; ++i)
    if (mask[i])
      {
        const ComparisonResult result = compare(v1[i], v2[i]);
        if (result != ComparisonResult::equal)
          return result;
      }
 
  return ComparisonResult::equal;
}
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif