Merge pull request #73 from lanl/jmm/robust-utils.hpp

add robust utils header

doc/sphinx/src/using.rst

@@ -128,6 +128,72 @@ Please note that none of these functions are thread or MPI aware. In a parallel
 
 as appropriate.
 
+robust_utils.hpp
+^^^^^^^^^^^^^^^^^^^
+
+``robust_utils.hpp`` contains small utility functions for numerical
+robustness, especially around floating point numbers. The available
+functionality is contained in the namespace ``PortsOfCall::Robust`` and includes:
+
+* ``constexpr auto SMALL<T>()`` returns a small number of type ``T``.
+* ``constexpr auto EPS<T>()`` returns a value of type ``T`` close to machine epsilon.
+* ``constexpr auto min_exp_arg<T>()`` returns the smallest safe value of type ``T`` to pass into an exponent.
+* ``constexpr auto max_exp_exp_arg<T>()`` returns the max safe value of type ``T`` to pass into an exponent.
+* ``auto make_positive(const T val)`` makes the argument of type ``T`` positive.
+
+where here all functionality is templated on type ``T`` and marked
+with ``PORTABLE_INLINE_FUNCTION``. The default type ``T`` is always
+``Real``.
+
+The function
+
+.. code-block:: cpp
+
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real make_bounded(const T val, const T vmin, const T vmax);
+
+bounds ``val`` between ``vmin`` and ``vmax``, exclusive. Note this is
+slightly different than ``std::clamp``, which uses inclusive bounds.
+
+The function
+
+.. code-block:: cpp
+
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION int sgn(const T &val);
+
+returns the sign of a quantity ``val``.
+
+.. note::
+
+  Note this implementation **never** returns zero. It **always**
+  returns :math:`\pm 1`.
+
+The function
+
+.. code-block:: cpp
+
+  template <typename A, typename B>
+  PORTABLE_FORCEINLINE_FUNCTION auto ratio(const A &a, const B &b)
+
+computes the ratio :math:`A/B` but in a way robust to 0/0 errors. If
+both :math:`A` and :math:`B` are zero, this function will return 0. If
+:math:`|A| > 0` and :math:`B=0`, then it will return a very large,
+possibly (but not guaranteed to be) infinite number.
+
+The function
+
+.. code-block:: cpp
+
+  template <typename T>
+  PORTABLE_FORCEINLINE_FUNCTION T safe_arg_exp(const T &x)
+
+returns exponentiation in such a way that avoids floating point
+exceptions. For very large negative inputs, it returns 0. For very
+large positive ones, it returns
+``std::numeric_limits<T>::infinity()``.
+
 macros_arrays.hpp
 ^^^^^^^^^^^^^^^^^^^
 

ports-of-call/robust_utils.hpp

@@ -0,0 +1,77 @@
+//------------------------------------------------------------------------------
+// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
+// program was produced under U.S. Government contract 89233218CNA000001
+// for Los Alamos National Laboratory (LANL), which is operated by Triad
+// National Security, LLC for the U.S.  Department of Energy/National
+// Nuclear Security Administration. All rights in the program are
+// reserved by Triad National Security, LLC, and the U.S. Department of
+// Energy/National Nuclear Security Administration. The Government is
+// granted for itself and others acting on its behalf a nonexclusive,
+// paid-up, irrevocable worldwide license in this material to reproduce,
+// prepare derivative works, distribute copies to the public, perform
+// publicly and display publicly, and to permit others to do so.
+//------------------------------------------------------------------------------
+
+#ifndef PORTS_OF_CALL_ROBUST_UTILS_HPP_
+#define PORTS_OF_CALL_ROBUST_UTILS_HPP_
+
+#include <cmath>
+#include <limits>
+#include <ports-of-call/portability.hpp>
+
+namespace PortsOfCall {
+namespace Robust {
+
+template <typename T = Real>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto SMALL() {
+  return 10 * std::numeric_limits<T>::min();
+}
+
+template <typename T = Real>
+PORTABLE_FORCEINLINE_FUNCTION constexpr auto EPS() {
+  return 10 * std::numeric_limits<T>::epsilon();
+}
+
+template <typename T = Real>
+PORTABLE_FORCEINLINE_FUNCTION constexpr T min_exp_arg() {
+  return (std::numeric_limits<T>::min_exponent - 1) * M_LN2;
+}
+
+template <typename T = Real>
+PORTABLE_FORCEINLINE_FUNCTION constexpr T max_exp_arg() {
+  return std::numeric_limits<T>::max_exponent * M_LN2;
+}
+
+template <typename T>
+PORTABLE_FORCEINLINE_FUNCTION auto make_positive(const T val) {
+  return std::max(val, EPS<T>());
+}
+
+template <typename T>
+PORTABLE_FORCEINLINE_FUNCTION Real make_bounded(const T val, const T vmin, const T vmax) {
+  return std::min(std::max(val, vmin + EPS<T>()), vmax * (1.0 - EPS<T>()));
+}
+
+template <typename T>
+PORTABLE_FORCEINLINE_FUNCTION int sgn(const T &val) {
+  return (T(0) <= val) - (val < T(0));
+}
+
+template <typename A, typename B>
+PORTABLE_FORCEINLINE_FUNCTION auto ratio(const A &a, const B &b) {
+  B mask = static_cast<B>(std::abs(b) < SMALL());
+  B denom = mask * sgn(b) * SMALL<B>() + (1 - mask) * b;
+  return a / (b + sgn(b) * SMALL<B>());
+}
+
+template <typename T>
+PORTABLE_FORCEINLINE_FUNCTION T safe_arg_exp(const T &x) {
+  return x < min_exp_arg<T>()   ? 0.0
+         : x > max_exp_arg<T>() ? std::numeric_limits<T>::infinity()
+                                : std::exp(x);
+}
+
+} // namespace Robust
+} // namespace PortsOfCall
+
+#endif // PORTS_OF_CALL_ROBUST_UTILS_HPP_