diff --git a/.github/workflows/gcc.yml b/.github/workflows/gcc.yml index 25b216e9bfb..cbe895c1c80 100644 --- a/.github/workflows/gcc.yml +++ b/.github/workflows/gcc.yml @@ -254,11 +254,13 @@ jobs: .github/workflows/dependencies/dependencies_clang-tidy-apt-llvm.sh 21 .github/workflows/dependencies/dependencies_ccache.sh - name: install vir-simd + # TODO: back to the release tarball once the vector math header is in one. + # Until then this branch is what provides vir/simd_vecmath.h, which + # amrex::Math's SIMD transcendentals forward to. Without it AMReX + # still builds, and evaluates them one lane at a time. run: | - wget https://github.com/mattkretz/vir-simd/archive/refs/tags/v0.4.4.tar.gz - tar -xvf v0.4.4.tar.gz - rm -rf v0.4.4.tar.gz - cmake -S vir-simd-0.4.4 -B vir-simd-build + git clone --depth 1 --branch topic-vecmath https://github.com/ax3l/vir-simd.git vir-simd-src + cmake -S vir-simd-src -B vir-simd-build sudo cmake --build vir-simd-build --target install - name: Set Up Cache uses: actions/cache@v6 diff --git a/Docs/sphinx_documentation/source/Basics.rst b/Docs/sphinx_documentation/source/Basics.rst index f7c257dec13..a636fba5dee 100644 --- a/Docs/sphinx_documentation/source/Basics.rst +++ b/Docs/sphinx_documentation/source/Basics.rst @@ -3066,6 +3066,76 @@ of signatures (including the :cpp:`CompileTimeOptions` variants) and is identical to :cpp:`ParallelFor` on GPU, but does not add the SIMD pragma on CPU. +.. _sec:basics:simdmath: + +SIMD Math Functions +=================== + +When AMReX is built with SIMD support (CMake option ``AMReX_SIMD=ON``, see +:ref:`sec:build:cmake`), the math functions in ``AMReX_Math.H`` accept SIMD +variables in addition to :cpp:`float` and :cpp:`double`. A kernel can therefore +be written once and instantiated for both, which is the same single-source style +used by :cpp:`amrex::ParallelForSIMD` and :cpp:`amrex::ParticleReduceSIMD`: + +.. highlight:: c++ + +:: + + #include + + // T_Real is amrex::ParticleReal in a scalar build and a SIMD type in a + // vectorized one + template + AMREX_FORCE_INLINE + void focus (T_Real & AMREX_RESTRICT y, T_Real & AMREX_RESTRICT py, + T_Real const & AMREX_RESTRICT omega, amrex::Real ds) + { + T_Real const ch = amrex::Math::cosh(omega * T_Real(ds)); + T_Real const sh = amrex::Math::sinh(omega * T_Real(ds)); + + T_Real const y0 = y; + y = ch * y0 + sh / omega * py; + py = omega * sh * y0 + ch * py; + } + +The following functions have SIMD overloads: :cpp:`sin`, :cpp:`cos`, +:cpp:`tan`, :cpp:`asin`, :cpp:`acos`, :cpp:`atan`, :cpp:`atan2`, :cpp:`sinh`, +:cpp:`cosh`, :cpp:`tanh`, :cpp:`asinh`, :cpp:`acosh`, :cpp:`atanh`, +:cpp:`exp`, :cpp:`exp2`, :cpp:`expm1`, :cpp:`log`, :cpp:`log2`, :cpp:`log10`, +:cpp:`log1p`, :cpp:`pow`, :cpp:`sqrt`, :cpp:`cbrt`, :cpp:`hypot`, :cpp:`erf`, +:cpp:`erfc`, :cpp:`abs`, :cpp:`sincos` and :cpp:`sincospi`. + +.. note:: + + Call these functions **fully qualified**, as :cpp:`amrex::Math::sinh(x)`. An + unqualified :cpp:`sinh(x)` on a SIMD argument resolves to the SIMD library's + own overload through argument-dependent lookup, which evaluates the function + one lane at a time. A :cpp:`using amrex::Math::sinh;` declaration does not + change that, because the library overload is the more specialized candidate. + +Vectorizing the transcendentals +------------------------------- + +SIMD hardware has instructions for :cpp:`sqrt` and :cpp:`abs`, but not for the +transcendental functions. Those have to be evaluated by a vector math library, +such as glibc's ``libmvec``, which computes a whole SIMD register worth of +results per call. AMReX does not implement that itself: the SIMD overloads +forward to whatever the SIMD provider offers, which is a vector math library +call where the provider has one, and one scalar call per lane otherwise. + +Nothing has to be configured for this, and no compiler flags are involved. +Whether the fast path is available depends on the provider and the platform, +and both cases give correct results. + +.. warning:: + + Vector math libraries trade accuracy for speed: glibc's ``libmvec`` + documents a maximum error of 4 ULP, where its scalar routines stay below + 1 ULP. Where they are used, results differ slightly from a scalar build and + are not bit-wise reproducible against one. Consult the SIMD provider's + documentation if your application needs the accuracy of the scalar routines + or checks ``errno`` after math calls. + Ghost Cells =========== diff --git a/Docs/sphinx_documentation/source/Particle.rst b/Docs/sphinx_documentation/source/Particle.rst index 19950c500f8..2e48e53a017 100644 --- a/Docs/sphinx_documentation/source/Particle.rst +++ b/Docs/sphinx_documentation/source/Particle.rst @@ -622,6 +622,11 @@ with an :cpp:`amrex::SIMDindex`, provided by ``AMReX_ReduceSIMD.H``). See For a complete example, including a benchmark against the scalar entry points, see ``Tests/Particles/ParticleReduceSIMD``. +Kernels that call transcendental functions on SIMD variables should use the +math functions in ``AMReX_Math.H``, see :ref:`sec:basics:simdmath`. The SIMD +library evaluates its own :cpp:`sin`, :cpp:`sinh` and friends one lane at a +time, which can make a vectorized kernel slower than the scalar one. + .. _sec:Particles:Interacting: Interacting with Mesh Data diff --git a/Src/Base/AMReX_Math.H b/Src/Base/AMReX_Math.H index 8b1addf8271..bd620af30f9 100644 --- a/Src/Base/AMReX_Math.H +++ b/Src/Base/AMReX_Math.H @@ -141,21 +141,6 @@ namespace detail { } /// \endcond -#ifdef AMREX_USE_SIMD -//! Return sine and cosine of given number -template -requires (amrex::simd::stdx::is_simd_v) -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -std::pair sincos (T_Real x) -{ - using namespace amrex::simd::stdx; - std::pair r; - r.first = sin(x); - r.second = cos(x); - return r; -} -#endif - //! Return sine and cosine of given number AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE std::pair sincos (double x) @@ -184,22 +169,6 @@ std::pair sincos (float x) return r; } -#ifdef AMREX_USE_SIMD -//! Return sin(pi*x) and cos(pi*x) given x -template -requires (amrex::simd::stdx::is_simd_v) -AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -std::pair sincospi (T_Real x) -{ - using namespace amrex::simd::stdx; - T_Real const px = pi() * x; - std::pair r; - r.first = sin(px); - r.second = cos(px); - return r; -} -#endif - //! Return sin(pi*x) and cos(pi*x) given x AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE std::pair sincospi (double x) @@ -228,6 +197,608 @@ std::pair sincospi (float x) return r; } +/* + * Transcendental functions + * + * Each function comes in a scalar and, with AMReX_SIMD=ON, a SIMD overload, so + * that a compute kernel written once can be instantiated for both. Call them + * qualified as amrex::Math::sinh(x): an unqualified call on a SIMD argument + * resolves to the SIMD library's own (scalar, element-wise) overload through + * argument-dependent lookup instead. + * + * SIMD hardware has instructions for sqrt and abs but not for the + * transcendental functions, so the SIMD overloads are only as fast as what the + * SIMD provider offers for them: a call into a vector math library where it has + * one, otherwise one scalar call per lane. See amrex::simd::smath. Where a + * vector math library answers, results are accurate to a few ULP rather than + * correctly rounded, and are not bit-wise identical to a scalar build. + */ + +//! Return the sine of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real sin (T_Real x) +{ + return std::sin(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the sine of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real sin (T_Real const& x) +{ + return amrex::simd::smath::sin(x); +} +#endif + +//! Return the cosine of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real cos (T_Real x) +{ + return std::cos(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the cosine of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real cos (T_Real const& x) +{ + return amrex::simd::smath::cos(x); +} +#endif + +//! Return the tangent of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real tan (T_Real x) +{ + return std::tan(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the tangent of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real tan (T_Real const& x) +{ + return amrex::simd::smath::tan(x); +} +#endif + +//! Return the arc sine of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real asin (T_Real x) +{ + return std::asin(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the arc sine of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real asin (T_Real const& x) +{ + return amrex::simd::smath::asin(x); +} +#endif + +//! Return the arc cosine of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real acos (T_Real x) +{ + return std::acos(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the arc cosine of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real acos (T_Real const& x) +{ + return amrex::simd::smath::acos(x); +} +#endif + +//! Return the arc tangent of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real atan (T_Real x) +{ + return std::atan(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the arc tangent of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real atan (T_Real const& x) +{ + return amrex::simd::smath::atan(x); +} +#endif + +//! Return the hyperbolic sine of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real sinh (T_Real x) +{ + return std::sinh(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the hyperbolic sine of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real sinh (T_Real const& x) +{ + return amrex::simd::smath::sinh(x); +} +#endif + +//! Return the hyperbolic cosine of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real cosh (T_Real x) +{ + return std::cosh(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the hyperbolic cosine of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real cosh (T_Real const& x) +{ + return amrex::simd::smath::cosh(x); +} +#endif + +//! Return the hyperbolic tangent of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real tanh (T_Real x) +{ + return std::tanh(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the hyperbolic tangent of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real tanh (T_Real const& x) +{ + return amrex::simd::smath::tanh(x); +} +#endif + +//! Return the inverse hyperbolic sine of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real asinh (T_Real x) +{ + return std::asinh(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the inverse hyperbolic sine of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real asinh (T_Real const& x) +{ + return amrex::simd::smath::asinh(x); +} +#endif + +//! Return the inverse hyperbolic cosine of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real acosh (T_Real x) +{ + return std::acosh(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the inverse hyperbolic cosine of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real acosh (T_Real const& x) +{ + return amrex::simd::smath::acosh(x); +} +#endif + +//! Return the inverse hyperbolic tangent of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real atanh (T_Real x) +{ + return std::atanh(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the inverse hyperbolic tangent of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real atanh (T_Real const& x) +{ + return amrex::simd::smath::atanh(x); +} +#endif + +//! Return the base-e exponential of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real exp (T_Real x) +{ + return std::exp(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the base-e exponential of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real exp (T_Real const& x) +{ + return amrex::simd::smath::exp(x); +} +#endif + +//! Return the base-2 exponential of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real exp2 (T_Real x) +{ + return std::exp2(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the base-2 exponential of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real exp2 (T_Real const& x) +{ + return amrex::simd::smath::exp2(x); +} +#endif + +//! Return the base-e exponential minus one of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real expm1 (T_Real x) +{ + return std::expm1(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the base-e exponential minus one of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real expm1 (T_Real const& x) +{ + return amrex::simd::smath::expm1(x); +} +#endif + +//! Return the natural logarithm of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real log (T_Real x) +{ + return std::log(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the natural logarithm of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real log (T_Real const& x) +{ + return amrex::simd::smath::log(x); +} +#endif + +//! Return the base-2 logarithm of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real log2 (T_Real x) +{ + return std::log2(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the base-2 logarithm of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real log2 (T_Real const& x) +{ + return amrex::simd::smath::log2(x); +} +#endif + +//! Return the base-10 logarithm of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real log10 (T_Real x) +{ + return std::log10(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the base-10 logarithm of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real log10 (T_Real const& x) +{ + return amrex::simd::smath::log10(x); +} +#endif + +//! Return the natural logarithm of one plus the argument of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real log1p (T_Real x) +{ + return std::log1p(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the natural logarithm of one plus the argument of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real log1p (T_Real const& x) +{ + return amrex::simd::smath::log1p(x); +} +#endif + +//! Return the cube root of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real cbrt (T_Real x) +{ + return std::cbrt(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the cube root of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real cbrt (T_Real const& x) +{ + return amrex::simd::smath::cbrt(x); +} +#endif + +//! Return the error function of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real erf (T_Real x) +{ + return std::erf(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the error function of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real erf (T_Real const& x) +{ + return amrex::simd::smath::erf(x); +} +#endif + +//! Return the complementary error function of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real erfc (T_Real x) +{ + return std::erfc(x); +} + +#ifdef AMREX_USE_SIMD +//! Return the complementary error function of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real erfc (T_Real const& x) +{ + return amrex::simd::smath::erfc(x); +} +#endif + +//! Return x raised to the power y +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real pow (T_Real x, T_Real y) +{ + return std::pow(x, y); +} + +#ifdef AMREX_USE_SIMD +//! Return x raised to the power y, for every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real pow (T_Real const& x, T_Real const& y) +{ + return amrex::simd::smath::pow(x, y); +} +#endif + +//! Return arc tangent of y/x, using their signs +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real atan2 (T_Real y, T_Real x) +{ + return std::atan2(y, x); +} + +#ifdef AMREX_USE_SIMD +//! Return arc tangent of y/x, using their signs, for every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real atan2 (T_Real const& y, T_Real const& x) +{ + return amrex::simd::smath::atan2(y, x); +} +#endif + +//! Return square root of x*x + y*y +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real hypot (T_Real x, T_Real y) +{ + return std::hypot(x, y); +} + +#ifdef AMREX_USE_SIMD +//! Return square root of x*x + y*y, for every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real hypot (T_Real const& x, T_Real const& y) +{ + // not amrex::simd::smath: hypot is one of the functions a SIMD library + // usually implements itself, with SIMD instructions and the overflow + // fixups a vector math library's version would skip + return amrex::simd::stdx::hypot(x, y); +} +#endif + +//! Return the square root of the given number +template +requires (std::is_floating_point_v) +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +T_Real sqrt (T_Real x) +{ + return std::sqrt(x); +} + +#ifdef AMREX_USE_SIMD +/** Return the square root of every SIMD lane + * + * Unlike the transcendentals above, this maps onto a hardware instruction and + * is handed to the SIMD library directly. No vector math library is involved + * and the result is correctly rounded. + */ +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real sqrt (T_Real const& x) +{ + return amrex::simd::stdx::sqrt(x); +} + +/** Return the absolute value of every SIMD lane + * + * @see sqrt for why no vector math library is involved + */ +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +T_Real abs (T_Real const& x) +{ + return amrex::simd::stdx::abs(x); +} +#endif + +#ifdef AMREX_USE_SIMD +/** Return sine and cosine of every SIMD lane + * + * Evaluated with the SIMD library's own sin and cos, which inline into the + * caller. This is the one place that does not hand the work to a vector math + * library: amrex::Math::sin and amrex::Math::cos do. + * + * The difference is what each tends to be used for. A call into a vector math + * library is a scheduling barrier, so the surrounding arithmetic cannot overlap + * with the transcendental. Where the transcendental dominates a kernel that is + * worth paying for, which is why the single-result overloads route. sincos + * usually sits in kernels that are mostly other arithmetic instead -- a + * rotation, a coordinate transform -- and there the barrier costs more than the + * faster transcendental saves: on a quaternion spin rotation, AVX2, one thread + * pinned to a core, 3.64 ms this way against 3.94 ms routed. + * + * For a kernel the transcendental does dominate, call amrex::Math::sin and + * amrex::Math::cos separately. + * + * No fused vector sincos exists to call in any case. glibc's libmvec has one on + * x86-64, but its ABI returns both results through vectors of pointers, so the + * callee scatters its output a lane at a time -- about 1.4x the time of two + * independent vector calls. AArch64 has none, and neither has sincospi. + */ +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +std::pair sincos (T_Real const& x) +{ + std::pair r; + r.first = amrex::simd::stdx::sin(x); + r.second = amrex::simd::stdx::cos(x); + return r; +} + +//! Return sin(pi*x) and cos(pi*x) of every SIMD lane +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +std::pair sincospi (T_Real const& x) +{ + return amrex::Math::sincos(T_Real(pi()) * x); +} +#endif + //! Return pow(x, Power), where Power is an integer known at compile time template requires (!std::integral || Power >= 0) diff --git a/Src/Base/AMReX_SIMD.H b/Src/Base/AMReX_SIMD.H index 18daf09d612..8f82ca66ea6 100644 --- a/Src/Base/AMReX_SIMD.H +++ b/Src/Base/AMReX_SIMD.H @@ -12,6 +12,11 @@ # if __cplusplus >= 202002L # include # endif +// Transcendentals evaluated by a vector math library, where the provider has +// them. Older providers do not, hence the __has_include. +# if __has_include() +# include +# endif #endif #include @@ -193,6 +198,26 @@ namespace amrex::simd #endif } + /** Where the SIMD math functions in amrex::Math come from + * + * SIMD hardware has instructions for sqrt and abs but not for the + * transcendental functions, so a SIMD library has to evaluate those one + * lane at a time unless it can hand a whole register to a vector math + * library such as glibc's libmvec. Providers that can offer them + * separately, and this alias names whichever set is available: the fast + * one where the provider has it, its own otherwise. Either way the results + * are correct; only the speed and the last few ULP differ. + * + * @see amrex::Math::sinh and the other transcendentals + */ +#ifdef AMREX_USE_SIMD +# ifdef VIR_HAVE_SIMD_VECMATH + namespace smath = vir::vecmath; +# else + namespace smath = vir::stdx; +# endif +#endif + // TODO: move to AMReX_REAL.H? #ifdef AMREX_USE_SIMD diff --git a/Tests/SIMD/main.cpp b/Tests/SIMD/main.cpp index 0c7f5173ddf..e98b1b0da37 100644 --- a/Tests/SIMD/main.cpp +++ b/Tests/SIMD/main.cpp @@ -11,8 +11,12 @@ #include #include +#include +#include + #include #include +#include #include #include @@ -55,6 +59,51 @@ void func_mc (ParticleReal& x, ParticleReal const& y) { x += y; } void func_cc (ParticleReal const& /*x*/, ParticleReal const& /*y*/) {} void func_mm (ParticleReal& x, ParticleReal& y) { x += y; y *= ParticleReal(2); } +#ifdef AMREX_USE_SIMD +// Compare a SIMD math overload against its scalar counterpart over a range. +// +// A vector math library is allowed to be less accurate than scalar libm; glibc's +// libmvec documents a maximum error of 4 ULP, so this uses a tolerance of 8 ULP. +template +int check_simd_math (char const* name, F_Simd const& f_simd, F_Scalar const& f_scalar, + typename T_Simd::value_type lo, typename T_Simd::value_type hi) +{ + using T = typename T_Simd::value_type; + constexpr std::size_t width = T_Simd::size(); + constexpr int nchunk = 16; + constexpr int npoint = nchunk * int(width); + constexpr T max_ulp = T(8); + + int err = 0; + for (int c = 0; c < nchunk; ++c) { + T in[width]; + for (std::size_t i = 0; i < width; ++i) { + in[i] = lo + (hi - lo) * T(c * int(width) + int(i)) / T(npoint - 1); + } + T_Simd x; + x.copy_from(in, simd::stdx::element_aligned); + + T_Simd const y = f_simd(x); + + for (std::size_t i = 0; i < width; ++i) { + T const ref = f_scalar(in[i]); + T const got = y[i]; + T const tol = max_ulp * std::numeric_limits::epsilon() + * amrex::max(Math::abs(ref), T(1)); + if (!(Math::abs(got - ref) <= tol)) { + ++err; + if (err <= 2) { + Print() << " " << name << " mismatch at x=" << double(in[i]) + << ": got " << double(got) << ", expected " << double(ref) << "\n"; + } + } + } + } + if (err != 0) { Print() << " " << name << ": FAILED (" << err << " lanes)\n"; } + return err; +} +#endif + // --------------------------------------------------------------------------- int main (int argc, char* argv[]) { @@ -514,6 +563,110 @@ int main (int argc, char* argv[]) } #endif // !AMREX_USE_GPU + // ================================================================ + // Test: amrex::Math transcendentals, scalar and SIMD + // ================================================================ + { + int err = 0; + + // The scalar overloads must exist for every build, so that a kernel + // written once compiles with AMReX_SIMD both ON and OFF. + { + constexpr Real x = Real(0.5); + if (!amrex::almostEqual(Math::sinh(x), std::sinh(x))) { ++err; } + if (!amrex::almostEqual(Math::cosh(x), std::cosh(x))) { ++err; } + if (!amrex::almostEqual(Math::exp(x), std::exp(x))) { ++err; } + if (!amrex::almostEqual(Math::sqrt(x), std::sqrt(x))) { ++err; } + if (!amrex::almostEqual(Math::pow(x, Real(3)), std::pow(x, Real(3)))) { ++err; } + auto const [s, c] = Math::sincos(x); + if (!amrex::almostEqual(s, std::sin(x))) { ++err; } + if (!amrex::almostEqual(c, std::cos(x))) { ++err; } + } + +#ifdef AMREX_USE_SIMD + using V = simd::SIMDReal<>; + using T = Real; + +# define AMREX_CHECK_SIMD_MATH(FUNC, LO, HI) \ + err += check_simd_math(#FUNC, \ + [] (V const& v) { return Math::FUNC(v); }, \ + [] (T const v) { return std::FUNC(v); }, \ + T(LO), T(HI)) + + AMREX_CHECK_SIMD_MATH(sin, -6.0, 6.0); + AMREX_CHECK_SIMD_MATH(cos, -6.0, 6.0); + AMREX_CHECK_SIMD_MATH(tan, -1.5, 1.5); + AMREX_CHECK_SIMD_MATH(asin, -1.0, 1.0); + AMREX_CHECK_SIMD_MATH(acos, -1.0, 1.0); + AMREX_CHECK_SIMD_MATH(atan, -10.0, 10.0); + AMREX_CHECK_SIMD_MATH(sinh, -5.0, 5.0); + AMREX_CHECK_SIMD_MATH(cosh, -5.0, 5.0); + AMREX_CHECK_SIMD_MATH(tanh, -5.0, 5.0); + AMREX_CHECK_SIMD_MATH(asinh, -5.0, 5.0); + AMREX_CHECK_SIMD_MATH(acosh, 1.0, 10.0); + AMREX_CHECK_SIMD_MATH(atanh, -0.9, 0.9); + AMREX_CHECK_SIMD_MATH(exp, -5.0, 5.0); + AMREX_CHECK_SIMD_MATH(exp2, -5.0, 5.0); + AMREX_CHECK_SIMD_MATH(expm1, -1.0, 1.0); + AMREX_CHECK_SIMD_MATH(log, 0.1, 20.0); + AMREX_CHECK_SIMD_MATH(log2, 0.1, 20.0); + AMREX_CHECK_SIMD_MATH(log10, 0.1, 20.0); + AMREX_CHECK_SIMD_MATH(log1p, -0.9, 9.0); + AMREX_CHECK_SIMD_MATH(cbrt, -20.0, 20.0); + AMREX_CHECK_SIMD_MATH(erf, -3.0, 3.0); + AMREX_CHECK_SIMD_MATH(erfc, -3.0, 3.0); + AMREX_CHECK_SIMD_MATH(sqrt, 0.0, 20.0); + AMREX_CHECK_SIMD_MATH(abs, -20.0, 20.0); + +# undef AMREX_CHECK_SIMD_MATH + + // two-argument functions + err += check_simd_math("pow", + [] (V const& v) { return Math::pow(v, V(T(2.5))); }, + [] (T const v) { return std::pow(v, T(2.5)); }, + T(0.1), T(10.0)); + err += check_simd_math("atan2", + [] (V const& v) { return Math::atan2(v, V(T(2.0))); }, + [] (T const v) { return std::atan2(v, T(2.0)); }, + T(-10.0), T(10.0)); + err += check_simd_math("hypot", + [] (V const& v) { return Math::hypot(v, V(T(3.0))); }, + [] (T const v) { return std::hypot(v, T(3.0)); }, + T(-10.0), T(10.0)); + + // sincos and sincospi return both results at once + err += check_simd_math("sincos (sin)", + [] (V const& v) { return Math::sincos(v).first; }, + [] (T const v) { return std::sin(v); }, + T(-6.0), T(6.0)); + err += check_simd_math("sincos (cos)", + [] (V const& v) { return Math::sincos(v).second; }, + [] (T const v) { return std::cos(v); }, + T(-6.0), T(6.0)); + err += check_simd_math("sincospi (sin)", + [] (V const& v) { return Math::sincospi(v).first; }, + [] (T const v) { return std::sin(Math::pi() * v); }, + T(-2.0), T(2.0)); + err += check_simd_math("sincospi (cos)", + [] (V const& v) { return Math::sincospi(v).second; }, + [] (T const v) { return std::cos(Math::pi() * v); }, + T(-2.0), T(2.0)); + + Print() << "amrex::Math SIMD transcendentals (" +# ifdef VIR_HAVE_SIMD_VECMATH + << "vector math library via the SIMD provider" +# else + << "provider fallback, one call per lane" +# endif + << ", width " << int(V::size()) << "): " + << (err == 0 ? "PASSED" : "FAILED") << "\n"; +#else + Print() << "amrex::Math scalar transcendentals: " + << (err == 0 ? "PASSED" : "FAILED") << "\n"; +#endif + nerrors += err; + } + // ================================================================ // Final report // ================================================================