From 8c07c64619909c0325562231757dbb848d541556 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Mon, 24 Aug 2026 16:30:01 -0700 Subject: [PATCH 1/8] SIMD: Vectorized Transcendental Math Add SIMD overloads for the transcendental functions in amrex::Math, so that a kernel calling sinh, cos, exp, ... can be written once and instantiated for both scalar and SIMD types, like ParallelForSIMD and ParticleReduceSIMD. SIMD hardware has no transcendental instructions, and the SIMD library evaluates its own sin/sinh/... one lane at a time, which can make a vectorized kernel slower than a scalar one. Each overload here instead runs a short loop over its lanes, in the shape compilers replace with a single call into a vector math library such as glibc's libmvec. That replacement needs -fno-math-errno, and it needs the vector variants to be declared, which glibc only does under -ffast-math. The new CMake option AMReX_SIMD_VECMATH (on by default with AMReX_SIMD) adds the flag, and the variants AMReX uses are declared here, so no fast-math build is required. Where neither is available the lane loop is evaluated element by element, just like the SIMD library's own fallback. Co-Authored-By: Claude Opus 5 (1M context) --- Docs/sphinx_documentation/source/Basics.rst | 87 ++ .../source/BuildingAMReX.rst | 3 + Docs/sphinx_documentation/source/Particle.rst | 5 + Src/Base/AMReX_Math.H | 761 +++++++++++++++++- Tests/SIMD/main.cpp | 153 ++++ Tools/CMake/AMReXOptions.cmake | 5 + Tools/CMake/AMReXParallelBackends.cmake | 35 + 7 files changed, 1038 insertions(+), 11 deletions(-) diff --git a/Docs/sphinx_documentation/source/Basics.rst b/Docs/sphinx_documentation/source/Basics.rst index f7c257dec13..fc46072226d 100644 --- a/Docs/sphinx_documentation/source/Basics.rst +++ b/Docs/sphinx_documentation/source/Basics.rst @@ -3066,6 +3066,93 @@ 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 link such a library directly. Instead, each +SIMD math function contains a short loop over its lanes, written in the shape +that compilers recognize and replace with a single vector math library call. + +Two conditions must be met for that replacement to happen, and the CMake option +``AMReX_SIMD_VECMATH`` (on by default when ``AMReX_SIMD=ON``) takes care of +both. It adds ``-fno-math-errno`` to AMReX and to every downstream target, +because a math call that has to keep ``errno`` up to date may not be vectorized, +and it makes AMReX declare the vector variants of the functions it uses. For +clang it also adds ``-fveclib=libmvec``. + +Where the option cannot deliver, nothing breaks: the lane loop is then evaluated +one lane at a time, exactly as the SIMD library would have done. This is +currently the case + +* on any platform other than x86-64 Linux with glibc, +* with a glibc older than 2.22, and for everything except + :cpp:`sin`, :cpp:`cos`, :cpp:`exp`, :cpp:`log` and :cpp:`pow` with a glibc + older than 2.35 -- note that this is the *build-time* glibc, so a compiler + with an old sysroot (as shipped by conda-forge, for example) falls back even + on a recent host, and +* with clang, for the functions missing from its own vector function table, + which covers fewer functions than glibc provides. + +.. 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. Results therefore differ slightly from a scalar build. Because + ``-fno-math-errno`` applies to whole translation units, ordinary scalar + loops over math functions in downstream code may be auto-vectorized the same + way. Build with ``AMReX_SIMD_VECMATH=OFF`` 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/BuildingAMReX.rst b/Docs/sphinx_documentation/source/BuildingAMReX.rst index 5da179941c1..00c00c294f6 100644 --- a/Docs/sphinx_documentation/source/BuildingAMReX.rst +++ b/Docs/sphinx_documentation/source/BuildingAMReX.rst @@ -499,6 +499,9 @@ The list of available options is reported in the :ref:`table ` bel +------------------------------+-------------------------------------------------+-------------------------+-----------------------+ | AMReX_SIMD | Enable SIMD Primitives (using vir::stdx::simd) | NO | YES, NO | +------------------------------+-------------------------------------------------+-------------------------+-----------------------+ + | AMReX_SIMD_VECMATH | Let SIMD math functions call a vector math | YES | YES, NO | + | | library (adds -fno-math-errno) | | | + +------------------------------+-------------------------------------------------+-------------------------+-----------------------+ | AMReX_OMP | Build with OpenMP support | NO | YES, NO | +------------------------------+-------------------------------------------------+-------------------------+-----------------------+ | AMReX_GPU_BACKEND | Build with on-node, accelerated GPU backend | NONE | NONE, SYCL, HIP, CUDA | 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..134d7488fe8 100644 --- a/Src/Base/AMReX_Math.H +++ b/Src/Base/AMReX_Math.H @@ -19,6 +19,124 @@ # include #endif +/* + * Vector math library support for the SIMD math functions in amrex::Math + * + * With AMReX_SIMD=ON, the amrex::Math transcendentals below also accept SIMD + * types. They are implemented as a short loop over the SIMD lanes with a + * compile-time constant trip count (see amrex::Math::detail::map_lanes). + * Compilers that know a vector variant of the scalar function called inside + * that loop collapse the whole loop into a single call into a vector math + * library, such as glibc's libmvec. + * + * Two conditions must hold for that to happen: + * - the translation unit must be compiled with -fno-math-errno, so that a math + * call is free of side effects and may be vectorized, and + * - the vector variants must have been declared. glibc declares them only for + * translation units compiled with -ffast-math, which is too broad a hammer + * for most codes, so the ones AMReX uses are declared here instead. + * + * The CMake option AMReX_SIMD_VECMATH=ON adds the required compiler flags. If + * either condition is missing, the loop merely unrolls into scalar calls, which + * is what the SIMD library's own fallback does as well. + * + * Note on accuracy: 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. Results are also no longer bit-wise identical to a scalar build. + * + * @todo Move the SIMD math functions and the declarations below into their own + * header, e.g. AMReX_Math_SIMD.H, that AMReX_Math.H includes. Keeping the + * vector variant declarations out of a header that is pulled in almost + * everywhere makes it clearer where they take effect. + */ +/* AMREX_SIMD_VECMATH is defined when the lane loops below are expected to reach + * a vector math library. It can also be set from the outside (the CMake option + * AMReX_SIMD_VECMATH=ON does so for compilers we cannot detect from here, such + * as clang with -fveclib=libmvec). + */ +#if !defined(AMREX_SIMD_VECMATH) && defined(AMREX_USE_SIMD) && defined(__NO_MATH_ERRNO__) \ + && defined(__GLIBC__) && defined(__x86_64__) && defined(__GNUC__) && !defined(__clang__) +# include +# ifdef __GLIBC_PREREQ +# if __GLIBC_PREREQ(2,22) +//! Defined when SIMD transcendentals are expected to lower to vector math library calls +# define AMREX_SIMD_VECMATH 1 +# endif +# endif +#endif + +#if defined(AMREX_SIMD_VECMATH) && defined(__GLIBC__) && defined(__x86_64__) \ + && defined(__GNUC__) && !defined(__clang__) && !defined(__FAST_MATH__) +# include +# ifdef __GLIBC_PREREQ +// glibc's own declarations use "notinbranch": no masked variant is provided. +# define AMREX_VECMATH_FN __attribute__((__simd__("notinbranch"))) +# if __GLIBC_PREREQ(2,22) +extern "C" { + AMREX_VECMATH_FN double cos (double) noexcept; + AMREX_VECMATH_FN double exp (double) noexcept; + AMREX_VECMATH_FN double log (double) noexcept; + AMREX_VECMATH_FN double pow (double, double) noexcept; + AMREX_VECMATH_FN double sin (double) noexcept; + AMREX_VECMATH_FN float cosf (float) noexcept; + AMREX_VECMATH_FN float expf (float) noexcept; + AMREX_VECMATH_FN float logf (float) noexcept; + AMREX_VECMATH_FN float powf (float, float) noexcept; + AMREX_VECMATH_FN float sinf (float) noexcept; +# if defined(_GNU_SOURCE) && !defined(__APPLE__) + AMREX_VECMATH_FN void sincos (double, double*, double*) noexcept; + AMREX_VECMATH_FN void sincosf (float, float*, float*) noexcept; +# endif +} +# endif +# if __GLIBC_PREREQ(2,35) +extern "C" { + AMREX_VECMATH_FN double acos (double) noexcept; + AMREX_VECMATH_FN double acosh (double) noexcept; + AMREX_VECMATH_FN double asin (double) noexcept; + AMREX_VECMATH_FN double asinh (double) noexcept; + AMREX_VECMATH_FN double atan (double) noexcept; + AMREX_VECMATH_FN double atan2 (double, double) noexcept; + AMREX_VECMATH_FN double atanh (double) noexcept; + AMREX_VECMATH_FN double cbrt (double) noexcept; + AMREX_VECMATH_FN double cosh (double) noexcept; + AMREX_VECMATH_FN double erf (double) noexcept; + AMREX_VECMATH_FN double erfc (double) noexcept; + AMREX_VECMATH_FN double exp2 (double) noexcept; + AMREX_VECMATH_FN double expm1 (double) noexcept; + AMREX_VECMATH_FN double hypot (double, double) noexcept; + AMREX_VECMATH_FN double log10 (double) noexcept; + AMREX_VECMATH_FN double log1p (double) noexcept; + AMREX_VECMATH_FN double log2 (double) noexcept; + AMREX_VECMATH_FN double sinh (double) noexcept; + AMREX_VECMATH_FN double tan (double) noexcept; + AMREX_VECMATH_FN double tanh (double) noexcept; + AMREX_VECMATH_FN float acosf (float) noexcept; + AMREX_VECMATH_FN float acoshf (float) noexcept; + AMREX_VECMATH_FN float asinf (float) noexcept; + AMREX_VECMATH_FN float asinhf (float) noexcept; + AMREX_VECMATH_FN float atanf (float) noexcept; + AMREX_VECMATH_FN float atan2f (float, float) noexcept; + AMREX_VECMATH_FN float atanhf (float) noexcept; + AMREX_VECMATH_FN float cbrtf (float) noexcept; + AMREX_VECMATH_FN float coshf (float) noexcept; + AMREX_VECMATH_FN float erff (float) noexcept; + AMREX_VECMATH_FN float erfcf (float) noexcept; + AMREX_VECMATH_FN float exp2f (float) noexcept; + AMREX_VECMATH_FN float expm1f (float) noexcept; + AMREX_VECMATH_FN float hypotf (float, float) noexcept; + AMREX_VECMATH_FN float log10f (float) noexcept; + AMREX_VECMATH_FN float log1pf (float) noexcept; + AMREX_VECMATH_FN float log2f (float) noexcept; + AMREX_VECMATH_FN float sinhf (float) noexcept; + AMREX_VECMATH_FN float tanf (float) noexcept; + AMREX_VECMATH_FN float tanhf (float) noexcept; +} +# endif +# undef AMREX_VECMATH_FN +# endif +#endif + namespace amrex { // NOLINT(modernize-concat-nested-namespaces) /// \cond DOXYGEN_IGNORE inline namespace disabled { @@ -138,6 +256,79 @@ namespace detail { *cosx = std::cos(x); #endif } +#ifdef AMREX_USE_SIMD + /** Apply a scalar function to every lane of a SIMD variable + * + * The lane loop has a compile-time constant trip count and a body that is a + * single scalar function call. That is the shape an auto-vectorizer can + * replace with one call into a vector math library (see the note near the + * top of this file). When it does, the lane buffer never reaches memory; + * when it does not, the loop unrolls into scalar calls, which is what the + * SIMD library's own math fallback does as well. + * + * @param x the SIMD variable to transform + * @param f a scalar function applied to one lane + * @return a SIMD variable holding f(x) for every lane + */ + template + AMREX_FORCE_INLINE + T_Simd map_lanes (T_Simd const& x, F const& f) + { +#ifdef AMREX_SIMD_VECMATH + using T = typename T_Simd::value_type; + constexpr std::size_t width = T_Simd::size(); + constexpr std::size_t alignment = amrex::simd::stdx::memory_alignment_v; + + alignas(alignment) T lane[width]; + x.copy_to(lane, amrex::simd::stdx::vector_aligned); + + for (std::size_t i = 0; i < width; ++i) { + lane[i] = f(lane[i]); + } + + T_Simd r; + r.copy_from(lane, amrex::simd::stdx::vector_aligned); + return r; +#else + // Without a vector math library to call, the loop above would only add a + // lane buffer around the same scalar calls. Build the result element by + // element instead, like the SIMD library's own math fallback does. + return T_Simd([&] (auto i) { return f(x[i]); }); +#endif + } + + /** Apply a scalar function of two arguments to every lane of two SIMD variables + * + * @see map_lanes(T_Simd const&, F const&) + */ + template + AMREX_FORCE_INLINE + T_Simd map_lanes (T_Simd const& x, T_Simd const& y, F const& f) + { +#ifdef AMREX_SIMD_VECMATH + using T = typename T_Simd::value_type; + constexpr std::size_t width = T_Simd::size(); + constexpr std::size_t alignment = amrex::simd::stdx::memory_alignment_v; + + alignas(alignment) T lane_x[width]; + alignas(alignment) T lane_y[width]; + x.copy_to(lane_x, amrex::simd::stdx::vector_aligned); + y.copy_to(lane_y, amrex::simd::stdx::vector_aligned); + + for (std::size_t i = 0; i < width; ++i) { + lane_x[i] = f(lane_x[i], lane_y[i]); + } + + T_Simd r; + r.copy_from(lane_x, amrex::simd::stdx::vector_aligned); + return r; +#else + // see map_lanes(T_Simd const&, F const&) + return T_Simd([&] (auto i) { return f(x[i], y[i]); }); +#endif + } + +#endif } /// \endcond @@ -146,12 +337,15 @@ namespace detail { template requires (amrex::simd::stdx::is_simd_v) AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -std::pair sincos (T_Real x) +std::pair sincos (T_Real const& x) { - using namespace amrex::simd::stdx; + // Evaluated as two separate lane loops on purpose: no vector math library + // offers an auto-vectorizable sincos, because writing both results through + // pointers stops the vectorizer. Two vector calls still beat one scalar + // sincos per lane. std::pair r; - r.first = sin(x); - r.second = cos(x); + r.first = detail::map_lanes(x, [] (auto v) { return std::sin(v); }); + r.second = detail::map_lanes(x, [] (auto v) { return std::cos(v); }); return r; } #endif @@ -189,14 +383,9 @@ std::pair sincos (float x) template requires (amrex::simd::stdx::is_simd_v) AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE -std::pair sincospi (T_Real x) +std::pair sincospi (T_Real const& 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; + return sincos(T_Real(pi()) * x); } #endif @@ -228,6 +417,556 @@ 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. + * + * The SIMD overloads may lower to vector math library calls, see the note near + * the top of this file for the requirements and the accuracy implications. + */ + +//! 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 detail::map_lanes(x, [] (auto v) { return std::sin(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::cos(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::tan(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::asin(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::acos(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::atan(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::sinh(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::cosh(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::tanh(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::asinh(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::acosh(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::atanh(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::exp(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::exp2(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::expm1(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::log(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::log2(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::log10(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::log1p(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::cbrt(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::erf(v); }); +} +#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 detail::map_lanes(x, [] (auto v) { return std::erfc(v); }); +} +#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 detail::map_lanes(x, y, [] (auto v_x, auto v_y) { return std::pow(v_x, v_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 detail::map_lanes(y, x, [] (auto v_y, auto v_x) { return std::atan2(v_y, v_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) +{ + return detail::map_lanes(x, y, [] (auto v_x, auto v_y) { return std::hypot(v_x, v_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 + //! Return pow(x, Power), where Power is an integer known at compile time template requires (!std::integral || Power >= 0) diff --git a/Tests/SIMD/main.cpp b/Tests/SIMD/main.cpp index 0c7f5173ddf..9ed0c3b3469 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 AMREX_SIMD_VECMATH + << "vector math library enabled" +# else + << "vector math library not available, lane-wise fallback" +# 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 // ================================================================ diff --git a/Tools/CMake/AMReXOptions.cmake b/Tools/CMake/AMReXOptions.cmake index 98e52ff61ca..f746236e9f3 100644 --- a/Tools/CMake/AMReXOptions.cmake +++ b/Tools/CMake/AMReXOptions.cmake @@ -273,6 +273,11 @@ print_option( AMReX_MPI_THREAD_MULTIPLE ) option( AMReX_SIMD "Enable SIMD Primitives" OFF) print_option( AMReX_SIMD ) +cmake_dependent_option( AMReX_SIMD_VECMATH + "Let SIMD math functions call a vector math library (adds -fno-math-errno)" ON + "AMReX_SIMD" OFF) +print_option( AMReX_SIMD_VECMATH ) + option( AMReX_OMP "Enable OpenMP" OFF) print_option( AMReX_OMP ) diff --git a/Tools/CMake/AMReXParallelBackends.cmake b/Tools/CMake/AMReXParallelBackends.cmake index 0f542ced37d..d09a700d292 100644 --- a/Tools/CMake/AMReXParallelBackends.cmake +++ b/Tools/CMake/AMReXParallelBackends.cmake @@ -38,6 +38,41 @@ if (AMReX_SIMD) foreach(D IN LISTS AMReX_SPACEDIM) target_link_libraries(amrex_${D}d PUBLIC vir-simd::vir-simd) endforeach() + + # Vector math library for the SIMD math functions in AMReX_Math.H. + # + # A math function may only be vectorized if the compiler does not have to keep + # errno up to date, so -fno-math-errno is needed both here and in every + # downstream translation unit that calls amrex::Math with a SIMD argument. + if (AMReX_SIMD_VECMATH) + include(CheckCXXCompilerFlag) + + check_cxx_compiler_flag("-fno-math-errno" AMReX_HAS_FLAG_NO_MATH_ERRNO) + if (AMReX_HAS_FLAG_NO_MATH_ERRNO) + foreach(D IN LISTS AMReX_SPACEDIM) + target_compile_options(amrex_${D}d + PUBLIC $<$:-fno-math-errno>) + endforeach() + else () + message(WARNING "AMReX_SIMD_VECMATH: ${CMAKE_CXX_COMPILER_ID} does not accept " + "-fno-math-errno. SIMD math functions stay scalar.") + endif () + + # GCC finds the vector variants through the declarations in AMReX_Math.H. + # clang ignores those and uses a built-in mapping table instead, which it + # only consults with -fveclib. + if (CMAKE_CXX_COMPILER_ID MATCHES "Clang" AND CMAKE_SYSTEM_NAME STREQUAL "Linux" + AND AMReX_HAS_FLAG_NO_MATH_ERRNO) + check_cxx_compiler_flag("-fveclib=libmvec" AMReX_HAS_FLAG_VECLIB_LIBMVEC) + if (AMReX_HAS_FLAG_VECLIB_LIBMVEC) + foreach(D IN LISTS AMReX_SPACEDIM) + target_compile_options(amrex_${D}d + PUBLIC $<$:-fveclib=libmvec>) + target_compile_definitions(amrex_${D}d PUBLIC AMREX_SIMD_VECMATH=1) + endforeach() + endif () + endif () + endif () endif () # From 0bd600fe6fcc42295a5c0da5e4fb0a5c02dfa08a Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Mon, 24 Aug 2026 18:06:56 -0700 Subject: [PATCH 2/8] SIMD: Delegate the math fallback to the SIMD library Where no vector math library is reachable, the lane loops fell back to building the result element by element, on the assumption that this is what the SIMD library does anyway. That holds for most functions, but not for sin and cos: libstdc++ carries real vectorized implementations for those (Taylor series with quadrant folding), which a per-element fallback throws away. Measured at width 4, sin went from 1.21 ms to 5.24 ms per 2^20 evaluations, a 4.3x regression on two of the most frequently called functions. Pass the SIMD library's own overload to map_lanes and call it in the fallback instead, so AMReX inherits whatever the library implements well, now and later. Measured parity with calling the library directly (0.98x to 1.02x) across all functions. Also give the two-argument lambdas neutral parameter names, so that readability-suspicious-call-argument does not compare them against the parameter names of the SIMD library's declarations. Co-Authored-By: Claude Opus 5 (1M context) --- Src/Base/AMReX_Math.H | 194 +++++++++++++++++++++++++++--------------- 1 file changed, 126 insertions(+), 68 deletions(-) diff --git a/Src/Base/AMReX_Math.H b/Src/Base/AMReX_Math.H index 134d7488fe8..8e610ea8913 100644 --- a/Src/Base/AMReX_Math.H +++ b/Src/Base/AMReX_Math.H @@ -267,14 +267,17 @@ namespace detail { * SIMD library's own math fallback does as well. * * @param x the SIMD variable to transform - * @param f a scalar function applied to one lane - * @return a SIMD variable holding f(x) for every lane + * @param lane_fn the scalar function to apply to one lane + * @param simd_fn the SIMD library's own overload of the same function + * @return a SIMD variable holding the function value for every lane */ - template + template AMREX_FORCE_INLINE - T_Simd map_lanes (T_Simd const& x, F const& f) + T_Simd map_lanes (T_Simd const& x, F_Lane const& lane_fn, F_Simd const& simd_fn) { #ifdef AMREX_SIMD_VECMATH + static_cast(simd_fn); + using T = typename T_Simd::value_type; constexpr std::size_t width = T_Simd::size(); constexpr std::size_t alignment = amrex::simd::stdx::memory_alignment_v; @@ -283,17 +286,19 @@ namespace detail { x.copy_to(lane, amrex::simd::stdx::vector_aligned); for (std::size_t i = 0; i < width; ++i) { - lane[i] = f(lane[i]); + lane[i] = lane_fn(lane[i]); } T_Simd r; r.copy_from(lane, amrex::simd::stdx::vector_aligned); return r; #else - // Without a vector math library to call, the loop above would only add a - // lane buffer around the same scalar calls. Build the result element by - // element instead, like the SIMD library's own math fallback does. - return T_Simd([&] (auto i) { return f(x[i]); }); + // Without a vector math library to call, the loop above would only wrap a + // lane buffer around the same scalar calls. Hand the work to the SIMD + // library instead: it evaluates most functions one lane at a time too, but + // it does carry real vectorized implementations for a few of them. + static_cast(lane_fn); + return simd_fn(x); #endif } @@ -301,11 +306,14 @@ namespace detail { * * @see map_lanes(T_Simd const&, F const&) */ - template + template AMREX_FORCE_INLINE - T_Simd map_lanes (T_Simd const& x, T_Simd const& y, F const& f) + T_Simd map_lanes (T_Simd const& x, T_Simd const& y, + F_Lane const& lane_fn, F_Simd const& simd_fn) { #ifdef AMREX_SIMD_VECMATH + static_cast(simd_fn); + using T = typename T_Simd::value_type; constexpr std::size_t width = T_Simd::size(); constexpr std::size_t alignment = amrex::simd::stdx::memory_alignment_v; @@ -316,15 +324,16 @@ namespace detail { y.copy_to(lane_y, amrex::simd::stdx::vector_aligned); for (std::size_t i = 0; i < width; ++i) { - lane_x[i] = f(lane_x[i], lane_y[i]); + lane_x[i] = lane_fn(lane_x[i], lane_y[i]); } T_Simd r; r.copy_from(lane_x, amrex::simd::stdx::vector_aligned); return r; #else - // see map_lanes(T_Simd const&, F const&) - return T_Simd([&] (auto i) { return f(x[i], y[i]); }); + // see map_lanes(T_Simd const&, F_Lane const&, F_Simd const&) + static_cast(lane_fn); + return simd_fn(x, y); #endif } @@ -332,24 +341,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 const& x) -{ - // Evaluated as two separate lane loops on purpose: no vector math library - // offers an auto-vectorizable sincos, because writing both results through - // pointers stops the vectorizer. Two vector calls still beat one scalar - // sincos per lane. - std::pair r; - r.first = detail::map_lanes(x, [] (auto v) { return std::sin(v); }); - r.second = detail::map_lanes(x, [] (auto v) { return std::cos(v); }); - return r; -} -#endif - //! Return sine and cosine of given number AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE std::pair sincos (double x) @@ -378,17 +369,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 const& x) -{ - return sincos(T_Real(pi()) * x); -} -#endif - //! Return sin(pi*x) and cos(pi*x) given x AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE std::pair sincospi (double x) @@ -446,7 +426,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real sin (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::sin(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::sin(v); }, + [] (auto const& v) { return amrex::simd::stdx::sin(v); }); } #endif @@ -466,7 +448,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real cos (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::cos(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::cos(v); }, + [] (auto const& v) { return amrex::simd::stdx::cos(v); }); } #endif @@ -486,7 +470,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real tan (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::tan(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::tan(v); }, + [] (auto const& v) { return amrex::simd::stdx::tan(v); }); } #endif @@ -506,7 +492,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real asin (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::asin(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::asin(v); }, + [] (auto const& v) { return amrex::simd::stdx::asin(v); }); } #endif @@ -526,7 +514,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real acos (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::acos(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::acos(v); }, + [] (auto const& v) { return amrex::simd::stdx::acos(v); }); } #endif @@ -546,7 +536,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real atan (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::atan(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::atan(v); }, + [] (auto const& v) { return amrex::simd::stdx::atan(v); }); } #endif @@ -566,7 +558,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real sinh (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::sinh(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::sinh(v); }, + [] (auto const& v) { return amrex::simd::stdx::sinh(v); }); } #endif @@ -586,7 +580,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real cosh (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::cosh(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::cosh(v); }, + [] (auto const& v) { return amrex::simd::stdx::cosh(v); }); } #endif @@ -606,7 +602,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real tanh (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::tanh(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::tanh(v); }, + [] (auto const& v) { return amrex::simd::stdx::tanh(v); }); } #endif @@ -626,7 +624,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real asinh (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::asinh(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::asinh(v); }, + [] (auto const& v) { return amrex::simd::stdx::asinh(v); }); } #endif @@ -646,7 +646,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real acosh (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::acosh(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::acosh(v); }, + [] (auto const& v) { return amrex::simd::stdx::acosh(v); }); } #endif @@ -666,7 +668,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real atanh (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::atanh(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::atanh(v); }, + [] (auto const& v) { return amrex::simd::stdx::atanh(v); }); } #endif @@ -686,7 +690,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real exp (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::exp(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::exp(v); }, + [] (auto const& v) { return amrex::simd::stdx::exp(v); }); } #endif @@ -706,7 +712,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real exp2 (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::exp2(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::exp2(v); }, + [] (auto const& v) { return amrex::simd::stdx::exp2(v); }); } #endif @@ -726,7 +734,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real expm1 (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::expm1(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::expm1(v); }, + [] (auto const& v) { return amrex::simd::stdx::expm1(v); }); } #endif @@ -746,7 +756,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real log (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::log(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::log(v); }, + [] (auto const& v) { return amrex::simd::stdx::log(v); }); } #endif @@ -766,7 +778,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real log2 (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::log2(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::log2(v); }, + [] (auto const& v) { return amrex::simd::stdx::log2(v); }); } #endif @@ -786,7 +800,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real log10 (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::log10(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::log10(v); }, + [] (auto const& v) { return amrex::simd::stdx::log10(v); }); } #endif @@ -806,7 +822,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real log1p (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::log1p(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::log1p(v); }, + [] (auto const& v) { return amrex::simd::stdx::log1p(v); }); } #endif @@ -826,7 +844,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real cbrt (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::cbrt(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::cbrt(v); }, + [] (auto const& v) { return amrex::simd::stdx::cbrt(v); }); } #endif @@ -846,7 +866,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real erf (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::erf(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::erf(v); }, + [] (auto const& v) { return amrex::simd::stdx::erf(v); }); } #endif @@ -866,7 +888,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real erfc (T_Real const& x) { - return detail::map_lanes(x, [] (auto v) { return std::erfc(v); }); + return detail::map_lanes(x, + [] (auto v) { return std::erfc(v); }, + [] (auto const& v) { return amrex::simd::stdx::erfc(v); }); } #endif @@ -886,7 +910,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real pow (T_Real const& x, T_Real const& y) { - return detail::map_lanes(x, y, [] (auto v_x, auto v_y) { return std::pow(v_x, v_y); }); + return detail::map_lanes(x, y, + [] (auto a, auto b) { return std::pow(a, b); }, + [] (auto const& a, auto const& b) { return amrex::simd::stdx::pow(a, b); }); } #endif @@ -906,7 +932,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real atan2 (T_Real const& y, T_Real const& x) { - return detail::map_lanes(y, x, [] (auto v_y, auto v_x) { return std::atan2(v_y, v_x); }); + return detail::map_lanes(y, x, + [] (auto a, auto b) { return std::atan2(a, b); }, + [] (auto const& a, auto const& b) { return amrex::simd::stdx::atan2(a, b); }); } #endif @@ -926,7 +954,9 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real hypot (T_Real const& x, T_Real const& y) { - return detail::map_lanes(x, y, [] (auto v_x, auto v_y) { return std::hypot(v_x, v_y); }); + return detail::map_lanes(x, y, + [] (auto a, auto b) { return std::hypot(a, b); }, + [] (auto const& a, auto const& b) { return amrex::simd::stdx::hypot(a, b); }); } #endif @@ -967,6 +997,34 @@ T_Real abs (T_Real const& x) } #endif +#ifdef AMREX_USE_SIMD +/** Return sine and cosine of every SIMD lane + * + * Evaluated as two separate lane loops on purpose: no vector math library offers + * an auto-vectorizable sincos, because writing both results through pointers + * stops the vectorizer. Two vector calls still beat one scalar sincos per lane. + */ +template +requires (amrex::simd::stdx::is_simd_v) +AMREX_FORCE_INLINE +std::pair sincos (T_Real const& x) +{ + std::pair r; + r.first = amrex::Math::sin(x); + r.second = amrex::Math::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) From ec8bbdfa51e43ae3f3644fd278bba72f4ba5cffa Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Mon, 24 Aug 2026 18:19:22 -0700 Subject: [PATCH 3/8] SIMD: Make AMReX_SIMD_VECMATH opt-in A full ctest sweep with the option on turned up a failure in Particles_ParticleReduceSIMD, where two formulations of the same scalar sum are required to agree bit for bit. The component that mismatched is a plain weighted moment, dsy*p_w, with no math function anywhere near it: passing -fno-math-errno alone changed what the compiler was willing to inline and contract, and the two formulations drifted apart by one ULP. That is the whole argument for not enabling this by default. The flag applies to entire translation units, AMReX hands it to downstream targets as well, and its reach is not limited to the math functions it is meant to speed up. Default it to OFF, so that turning it on is a deliberate choice, and say so in the docs next to the accuracy warning. The strict comparisons in the test now scale their tolerance with the flag, so the suite passes in both configurations (77/77 either way). Co-Authored-By: Claude Opus 5 (1M context) --- Docs/sphinx_documentation/source/Basics.rst | 29 ++++++++++++------- .../source/BuildingAMReX.rst | 2 +- Tests/Particles/ParticleReduceSIMD/main.cpp | 16 ++++++++-- Tools/CMake/AMReXOptions.cmake | 2 +- 4 files changed, 35 insertions(+), 14 deletions(-) diff --git a/Docs/sphinx_documentation/source/Basics.rst b/Docs/sphinx_documentation/source/Basics.rst index fc46072226d..76890079307 100644 --- a/Docs/sphinx_documentation/source/Basics.rst +++ b/Docs/sphinx_documentation/source/Basics.rst @@ -3124,11 +3124,19 @@ SIMD math function contains a short loop over its lanes, written in the shape that compilers recognize and replace with a single vector math library call. Two conditions must be met for that replacement to happen, and the CMake option -``AMReX_SIMD_VECMATH`` (on by default when ``AMReX_SIMD=ON``) takes care of -both. It adds ``-fno-math-errno`` to AMReX and to every downstream target, -because a math call that has to keep ``errno`` up to date may not be vectorized, -and it makes AMReX declare the vector variants of the functions it uses. For -clang it also adds ``-fveclib=libmvec``. +``AMReX_SIMD_VECMATH`` takes care of both. It adds ``-fno-math-errno`` to AMReX +and to every downstream target, because a math call that has to keep ``errno`` +up to date may not be vectorized, and it makes AMReX declare the vector variants +of the functions it uses. For clang it also adds ``-fveclib=libmvec``. + +The option is off by default, and has to be turned on explicitly:: + + cmake -S . -B build -DAMReX_SIMD=ON -DAMReX_SIMD_VECMATH=ON + +It is opt-in because ``-fno-math-errno`` applies to whole translation units, not +just to the functions below. It changes what the compiler is allowed to inline +and contract, so results of unrelated floating-point code can shift by an ULP, +and two ways of writing the same expression need no longer agree bit for bit. Where the option cannot deliver, nothing breaks: the lane loop is then evaluated one lane at a time, exactly as the SIMD library would have done. This is @@ -3147,11 +3155,12 @@ currently the case 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. Results therefore differ slightly from a scalar build. Because - ``-fno-math-errno`` applies to whole translation units, ordinary scalar - loops over math functions in downstream code may be auto-vectorized the same - way. Build with ``AMReX_SIMD_VECMATH=OFF`` if your application needs the - accuracy of the scalar routines, or checks ``errno`` after math calls. + 1 ULP. Results therefore differ slightly from a scalar build. Ordinary + scalar loops over math functions in downstream code may be auto-vectorized + the same way, since the vector variants are declared for the whole + translation unit. Leave the option off if your application needs the accuracy + of the scalar routines, requires bitwise reproducible results, or checks + ``errno`` after math calls. Ghost Cells =========== diff --git a/Docs/sphinx_documentation/source/BuildingAMReX.rst b/Docs/sphinx_documentation/source/BuildingAMReX.rst index 00c00c294f6..fd4d699ecdf 100644 --- a/Docs/sphinx_documentation/source/BuildingAMReX.rst +++ b/Docs/sphinx_documentation/source/BuildingAMReX.rst @@ -499,7 +499,7 @@ The list of available options is reported in the :ref:`table ` bel +------------------------------+-------------------------------------------------+-------------------------+-----------------------+ | AMReX_SIMD | Enable SIMD Primitives (using vir::stdx::simd) | NO | YES, NO | +------------------------------+-------------------------------------------------+-------------------------+-----------------------+ - | AMReX_SIMD_VECMATH | Let SIMD math functions call a vector math | YES | YES, NO | + | AMReX_SIMD_VECMATH | Let SIMD math functions call a vector math | NO | YES, NO | | | library (adds -fno-math-errno) | | | +------------------------------+-------------------------------------------------+-------------------------+-----------------------+ | AMReX_OMP | Build with OpenMP support | NO | YES, NO | diff --git a/Tests/Particles/ParticleReduceSIMD/main.cpp b/Tests/Particles/ParticleReduceSIMD/main.cpp index a21b8a2f069..a52f8854b41 100644 --- a/Tests/Particles/ParticleReduceSIMD/main.cpp +++ b/Tests/Particles/ParticleReduceSIMD/main.cpp @@ -284,6 +284,18 @@ void compare_results (ReduceTupleT const& a, ReduceTupleT const& b, } } +/** Tolerance scale for comparisons that are otherwise required to be bitwise identical + * + * AMReX_SIMD_VECMATH=ON adds -fno-math-errno, which lets the compiler make + * different inlining and FMA contraction choices per formulation. Two ways of + * writing the same sum then no longer have to agree bit for bit. + */ +#ifdef AMREX_SIMD_VECMATH +constexpr Real bitwise_tol_scale = 4.; +#else +constexpr Real bitwise_tol_scale = 0.; +#endif + void correctness_tests (Geometry const& geom, DistributionMapping const& dm, BoxArray const& ba, Shifts const& shifts) { @@ -301,9 +313,9 @@ void correctness_tests (Geometry const& geom, DistributionMapping const& dm, auto const rc = run_variant_c(pc, shifts); // A, B and B2 evaluate identical arithmetic in identical order - compare_results(ra, rb, np, Real(0.), "A (SuperParticle) vs B (ptd,i)"); + compare_results(ra, rb, np, bitwise_tol_scale, "A (SuperParticle) vs B (ptd,i)"); auto const rb2 = run_variant_b2(pc, shifts); - compare_results(rb2, rb, np, Real(0.), + compare_results(rb2, rb, np, bitwise_tol_scale, "B2 (ptd,i by-ref) vs B (ptd,i)"); // C reassociates the sums across SIMD lanes compare_results(rc, rb, np, Real(100.), "C (SIMD) vs B (ptd,i)"); diff --git a/Tools/CMake/AMReXOptions.cmake b/Tools/CMake/AMReXOptions.cmake index f746236e9f3..dbfb60f341e 100644 --- a/Tools/CMake/AMReXOptions.cmake +++ b/Tools/CMake/AMReXOptions.cmake @@ -274,7 +274,7 @@ option( AMReX_SIMD "Enable SIMD Primitives" OFF) print_option( AMReX_SIMD ) cmake_dependent_option( AMReX_SIMD_VECMATH - "Let SIMD math functions call a vector math library (adds -fno-math-errno)" ON + "Let SIMD math functions call a vector math library (adds -fno-math-errno)" OFF "AMReX_SIMD" OFF) print_option( AMReX_SIMD_VECMATH ) From 8b1a417c71fc39dc37253ea7ebc6e576797d42fe Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Mon, 24 Aug 2026 18:40:11 -0700 Subject: [PATCH 4/8] SIMD: Separate the vector math option from its availability The CMake option AMReX_SIMD_VECMATH and the preprocessor macro it was named after are not the same thing. The option asks for vector math and adds the compiler flags for it; the macro says whether a vector math library is actually within reach. They differ on, for example, a toolchain with an old sysroot, where the flags apply but no vector variants exist. Sharing one name hid that distinction and produced a wrong test: the strict comparisons in ParticleReduceSIMD were relaxed based on the macro, while what perturbs the arithmetic is the flag. On a conda toolchain (glibc 2.17 sysroot) the flag applied, the macro stayed undefined, and the test failed. Rename the macro to AMREX_SIMD_HAS_VECMATH, say in the header how the two differ, and key the test tolerance off __NO_MATH_ERRNO__, which is what the compiler itself sets when the flag is in effect. Verified with AMReX_PRECISION=SINGLE as well, where the width 8 float path resolves to _ZGVdN8v_*f as expected. Co-Authored-By: Claude Opus 5 (1M context) --- Src/Base/AMReX_Math.H | 24 +++++++++++++-------- Tests/Particles/ParticleReduceSIMD/main.cpp | 8 ++++--- Tests/SIMD/main.cpp | 2 +- Tools/CMake/AMReXParallelBackends.cmake | 2 +- 4 files changed, 22 insertions(+), 14 deletions(-) diff --git a/Src/Base/AMReX_Math.H b/Src/Base/AMReX_Math.H index 8e610ea8913..3492018c141 100644 --- a/Src/Base/AMReX_Math.H +++ b/Src/Base/AMReX_Math.H @@ -49,23 +49,29 @@ * vector variant declarations out of a header that is pulled in almost * everywhere makes it clearer where they take effect. */ -/* AMREX_SIMD_VECMATH is defined when the lane loops below are expected to reach - * a vector math library. It can also be set from the outside (the CMake option - * AMReX_SIMD_VECMATH=ON does so for compilers we cannot detect from here, such - * as clang with -fveclib=libmvec). +/* AMREX_SIMD_HAS_VECMATH is defined when the lane loops below are expected to reach + * a vector math library. It can also be set from the outside: the CMake option + * AMReX_SIMD_VECMATH=ON does so for compilers this cannot be detected for from + * here, such as clang with -fveclib=libmvec. + * + * Note that the CMake option and this macro are not the same thing. The option + * asks for vector math and adds the compiler flags for it; this macro says + * whether a vector math library is actually within reach. On a toolchain with, + * say, an old sysroot, the option is on, its flags apply, and this macro stays + * undefined. */ -#if !defined(AMREX_SIMD_VECMATH) && defined(AMREX_USE_SIMD) && defined(__NO_MATH_ERRNO__) \ +#if !defined(AMREX_SIMD_HAS_VECMATH) && defined(AMREX_USE_SIMD) && defined(__NO_MATH_ERRNO__) \ && defined(__GLIBC__) && defined(__x86_64__) && defined(__GNUC__) && !defined(__clang__) # include # ifdef __GLIBC_PREREQ # if __GLIBC_PREREQ(2,22) //! Defined when SIMD transcendentals are expected to lower to vector math library calls -# define AMREX_SIMD_VECMATH 1 +# define AMREX_SIMD_HAS_VECMATH 1 # endif # endif #endif -#if defined(AMREX_SIMD_VECMATH) && defined(__GLIBC__) && defined(__x86_64__) \ +#if defined(AMREX_SIMD_HAS_VECMATH) && defined(__GLIBC__) && defined(__x86_64__) \ && defined(__GNUC__) && !defined(__clang__) && !defined(__FAST_MATH__) # include # ifdef __GLIBC_PREREQ @@ -275,7 +281,7 @@ namespace detail { AMREX_FORCE_INLINE T_Simd map_lanes (T_Simd const& x, F_Lane const& lane_fn, F_Simd const& simd_fn) { -#ifdef AMREX_SIMD_VECMATH +#ifdef AMREX_SIMD_HAS_VECMATH static_cast(simd_fn); using T = typename T_Simd::value_type; @@ -311,7 +317,7 @@ namespace detail { T_Simd map_lanes (T_Simd const& x, T_Simd const& y, F_Lane const& lane_fn, F_Simd const& simd_fn) { -#ifdef AMREX_SIMD_VECMATH +#ifdef AMREX_SIMD_HAS_VECMATH static_cast(simd_fn); using T = typename T_Simd::value_type; diff --git a/Tests/Particles/ParticleReduceSIMD/main.cpp b/Tests/Particles/ParticleReduceSIMD/main.cpp index a52f8854b41..ef533162ad8 100644 --- a/Tests/Particles/ParticleReduceSIMD/main.cpp +++ b/Tests/Particles/ParticleReduceSIMD/main.cpp @@ -286,11 +286,13 @@ void compare_results (ReduceTupleT const& a, ReduceTupleT const& b, /** Tolerance scale for comparisons that are otherwise required to be bitwise identical * - * AMReX_SIMD_VECMATH=ON adds -fno-math-errno, which lets the compiler make + * -fno-math-errno, which AMReX_SIMD_VECMATH=ON adds, lets the compiler make * different inlining and FMA contraction choices per formulation. Two ways of - * writing the same sum then no longer have to agree bit for bit. + * writing the same sum then no longer have to agree bit for bit. Keyed off the + * flag rather than off AMREX_SIMD_HAS_VECMATH, because the flag is what perturbs + * the arithmetic, whether or not a vector math library turns out to be reachable. */ -#ifdef AMREX_SIMD_VECMATH +#ifdef __NO_MATH_ERRNO__ constexpr Real bitwise_tol_scale = 4.; #else constexpr Real bitwise_tol_scale = 0.; diff --git a/Tests/SIMD/main.cpp b/Tests/SIMD/main.cpp index 9ed0c3b3469..42e3fd53124 100644 --- a/Tests/SIMD/main.cpp +++ b/Tests/SIMD/main.cpp @@ -653,7 +653,7 @@ int main (int argc, char* argv[]) T(-2.0), T(2.0)); Print() << "amrex::Math SIMD transcendentals (" -# ifdef AMREX_SIMD_VECMATH +# ifdef AMREX_SIMD_HAS_VECMATH << "vector math library enabled" # else << "vector math library not available, lane-wise fallback" diff --git a/Tools/CMake/AMReXParallelBackends.cmake b/Tools/CMake/AMReXParallelBackends.cmake index d09a700d292..954b9d59baf 100644 --- a/Tools/CMake/AMReXParallelBackends.cmake +++ b/Tools/CMake/AMReXParallelBackends.cmake @@ -68,7 +68,7 @@ if (AMReX_SIMD) foreach(D IN LISTS AMReX_SPACEDIM) target_compile_options(amrex_${D}d PUBLIC $<$:-fveclib=libmvec>) - target_compile_definitions(amrex_${D}d PUBLIC AMREX_SIMD_VECMATH=1) + target_compile_definitions(amrex_${D}d PUBLIC AMREX_SIMD_HAS_VECMATH=1) endforeach() endif () endif () From afee662ff6e0f0bcfd869616db65c2d4965c9542 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Mon, 24 Aug 2026 22:51:01 -0700 Subject: [PATCH 5/8] SIMD: Let the SIMD provider supply the vector math The first version of this got the transcendentals vectorized by shaping each one as a short lane loop and leaning on the auto-vectorizer to replace it with a call into glibc's libmvec. That worked, and cost more than it was worth: it needed -fno-math-errno on AMReX and on every downstream translation unit, and the vector variant declarations that made it possible applied to the whole translation unit, so ordinary scalar loops over math functions changed accuracy too. A full ctest sweep found that out the hard way, with ParticleReduceSIMD failing on a plain weighted moment that had no math function anywhere near it. Ask the SIMD provider instead. vir-simd calls libmvec's entry points directly, by their vector-function ABI names, so no flags and no auto-vectorization are involved and nothing outside these functions changes. AMReX picks whichever set the provider has: namespace smath = vir::vecmath; // a vector math library, where available namespace smath = vir::stdx; // the provider's own, one call per lane and every SIMD overload is a one-line forward to it. Guarded with __has_include, so this still builds against vir-simd releases without the vector math header, forwarding per lane. What is left is the part that has to be here: the amrex::Math overload set, so that a kernel written once compiles for scalar, SIMD and GPU. That layer is needed whatever the provider does, because an unqualified sinh(x) on a SIMD argument reaches the provider's own overload through argument-dependent lookup, and no using-declaration changes it -- the provider's overload either ties or wins partial ordering. Hence amrex::Math::sinh(x), qualified. hypot stays with vir::stdx: SIMD libraries generally implement it with SIMD instructions already, including overflow fixups a vector math library skips. Measured on an i9-12900H at width 4, no compiler flags: an ImpactX-shaped push goes from 16.8 ms to 5.9 ms against the current SIMD path, 3.7x the scalar loop; sinh alone is 7.2x. Accuracy within 4 ULP where a vector math library answers, bit-identical to the provider otherwise. Removes AMReX_SIMD_VECMATH, -fno-math-errno, the lane loops, the libmvec declarations, and the ParticleReduceSIMD tolerance those made necessary. Co-Authored-By: Claude Opus 5 (1M context) --- Docs/sphinx_documentation/source/Basics.rst | 48 +-- .../source/BuildingAMReX.rst | 3 - Src/Base/AMReX_Math.H | 317 ++---------------- Src/Base/AMReX_SIMD.H | 25 ++ Tests/Particles/ParticleReduceSIMD/main.cpp | 18 +- Tests/SIMD/main.cpp | 6 +- Tools/CMake/AMReXOptions.cmake | 5 - Tools/CMake/AMReXParallelBackends.cmake | 35 -- 8 files changed, 75 insertions(+), 382 deletions(-) diff --git a/Docs/sphinx_documentation/source/Basics.rst b/Docs/sphinx_documentation/source/Basics.rst index 76890079307..a636fba5dee 100644 --- a/Docs/sphinx_documentation/source/Basics.rst +++ b/Docs/sphinx_documentation/source/Basics.rst @@ -3119,48 +3119,22 @@ 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 link such a library directly. Instead, each -SIMD math function contains a short loop over its lanes, written in the shape -that compilers recognize and replace with a single vector math library call. - -Two conditions must be met for that replacement to happen, and the CMake option -``AMReX_SIMD_VECMATH`` takes care of both. It adds ``-fno-math-errno`` to AMReX -and to every downstream target, because a math call that has to keep ``errno`` -up to date may not be vectorized, and it makes AMReX declare the vector variants -of the functions it uses. For clang it also adds ``-fveclib=libmvec``. - -The option is off by default, and has to be turned on explicitly:: - - cmake -S . -B build -DAMReX_SIMD=ON -DAMReX_SIMD_VECMATH=ON - -It is opt-in because ``-fno-math-errno`` applies to whole translation units, not -just to the functions below. It changes what the compiler is allowed to inline -and contract, so results of unrelated floating-point code can shift by an ULP, -and two ways of writing the same expression need no longer agree bit for bit. - -Where the option cannot deliver, nothing breaks: the lane loop is then evaluated -one lane at a time, exactly as the SIMD library would have done. This is -currently the case - -* on any platform other than x86-64 Linux with glibc, -* with a glibc older than 2.22, and for everything except - :cpp:`sin`, :cpp:`cos`, :cpp:`exp`, :cpp:`log` and :cpp:`pow` with a glibc - older than 2.35 -- note that this is the *build-time* glibc, so a compiler - with an old sysroot (as shipped by conda-forge, for example) falls back even - on a recent host, and -* with clang, for the functions missing from its own vector function table, - which covers fewer functions than glibc provides. +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. Results therefore differ slightly from a scalar build. Ordinary - scalar loops over math functions in downstream code may be auto-vectorized - the same way, since the vector variants are declared for the whole - translation unit. Leave the option off if your application needs the accuracy - of the scalar routines, requires bitwise reproducible results, or checks - ``errno`` after math calls. + 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/BuildingAMReX.rst b/Docs/sphinx_documentation/source/BuildingAMReX.rst index fd4d699ecdf..5da179941c1 100644 --- a/Docs/sphinx_documentation/source/BuildingAMReX.rst +++ b/Docs/sphinx_documentation/source/BuildingAMReX.rst @@ -499,9 +499,6 @@ The list of available options is reported in the :ref:`table ` bel +------------------------------+-------------------------------------------------+-------------------------+-----------------------+ | AMReX_SIMD | Enable SIMD Primitives (using vir::stdx::simd) | NO | YES, NO | +------------------------------+-------------------------------------------------+-------------------------+-----------------------+ - | AMReX_SIMD_VECMATH | Let SIMD math functions call a vector math | NO | YES, NO | - | | library (adds -fno-math-errno) | | | - +------------------------------+-------------------------------------------------+-------------------------+-----------------------+ | AMReX_OMP | Build with OpenMP support | NO | YES, NO | +------------------------------+-------------------------------------------------+-------------------------+-----------------------+ | AMReX_GPU_BACKEND | Build with on-node, accelerated GPU backend | NONE | NONE, SYCL, HIP, CUDA | diff --git a/Src/Base/AMReX_Math.H b/Src/Base/AMReX_Math.H index 3492018c141..0e8a114a335 100644 --- a/Src/Base/AMReX_Math.H +++ b/Src/Base/AMReX_Math.H @@ -19,130 +19,6 @@ # include #endif -/* - * Vector math library support for the SIMD math functions in amrex::Math - * - * With AMReX_SIMD=ON, the amrex::Math transcendentals below also accept SIMD - * types. They are implemented as a short loop over the SIMD lanes with a - * compile-time constant trip count (see amrex::Math::detail::map_lanes). - * Compilers that know a vector variant of the scalar function called inside - * that loop collapse the whole loop into a single call into a vector math - * library, such as glibc's libmvec. - * - * Two conditions must hold for that to happen: - * - the translation unit must be compiled with -fno-math-errno, so that a math - * call is free of side effects and may be vectorized, and - * - the vector variants must have been declared. glibc declares them only for - * translation units compiled with -ffast-math, which is too broad a hammer - * for most codes, so the ones AMReX uses are declared here instead. - * - * The CMake option AMReX_SIMD_VECMATH=ON adds the required compiler flags. If - * either condition is missing, the loop merely unrolls into scalar calls, which - * is what the SIMD library's own fallback does as well. - * - * Note on accuracy: 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. Results are also no longer bit-wise identical to a scalar build. - * - * @todo Move the SIMD math functions and the declarations below into their own - * header, e.g. AMReX_Math_SIMD.H, that AMReX_Math.H includes. Keeping the - * vector variant declarations out of a header that is pulled in almost - * everywhere makes it clearer where they take effect. - */ -/* AMREX_SIMD_HAS_VECMATH is defined when the lane loops below are expected to reach - * a vector math library. It can also be set from the outside: the CMake option - * AMReX_SIMD_VECMATH=ON does so for compilers this cannot be detected for from - * here, such as clang with -fveclib=libmvec. - * - * Note that the CMake option and this macro are not the same thing. The option - * asks for vector math and adds the compiler flags for it; this macro says - * whether a vector math library is actually within reach. On a toolchain with, - * say, an old sysroot, the option is on, its flags apply, and this macro stays - * undefined. - */ -#if !defined(AMREX_SIMD_HAS_VECMATH) && defined(AMREX_USE_SIMD) && defined(__NO_MATH_ERRNO__) \ - && defined(__GLIBC__) && defined(__x86_64__) && defined(__GNUC__) && !defined(__clang__) -# include -# ifdef __GLIBC_PREREQ -# if __GLIBC_PREREQ(2,22) -//! Defined when SIMD transcendentals are expected to lower to vector math library calls -# define AMREX_SIMD_HAS_VECMATH 1 -# endif -# endif -#endif - -#if defined(AMREX_SIMD_HAS_VECMATH) && defined(__GLIBC__) && defined(__x86_64__) \ - && defined(__GNUC__) && !defined(__clang__) && !defined(__FAST_MATH__) -# include -# ifdef __GLIBC_PREREQ -// glibc's own declarations use "notinbranch": no masked variant is provided. -# define AMREX_VECMATH_FN __attribute__((__simd__("notinbranch"))) -# if __GLIBC_PREREQ(2,22) -extern "C" { - AMREX_VECMATH_FN double cos (double) noexcept; - AMREX_VECMATH_FN double exp (double) noexcept; - AMREX_VECMATH_FN double log (double) noexcept; - AMREX_VECMATH_FN double pow (double, double) noexcept; - AMREX_VECMATH_FN double sin (double) noexcept; - AMREX_VECMATH_FN float cosf (float) noexcept; - AMREX_VECMATH_FN float expf (float) noexcept; - AMREX_VECMATH_FN float logf (float) noexcept; - AMREX_VECMATH_FN float powf (float, float) noexcept; - AMREX_VECMATH_FN float sinf (float) noexcept; -# if defined(_GNU_SOURCE) && !defined(__APPLE__) - AMREX_VECMATH_FN void sincos (double, double*, double*) noexcept; - AMREX_VECMATH_FN void sincosf (float, float*, float*) noexcept; -# endif -} -# endif -# if __GLIBC_PREREQ(2,35) -extern "C" { - AMREX_VECMATH_FN double acos (double) noexcept; - AMREX_VECMATH_FN double acosh (double) noexcept; - AMREX_VECMATH_FN double asin (double) noexcept; - AMREX_VECMATH_FN double asinh (double) noexcept; - AMREX_VECMATH_FN double atan (double) noexcept; - AMREX_VECMATH_FN double atan2 (double, double) noexcept; - AMREX_VECMATH_FN double atanh (double) noexcept; - AMREX_VECMATH_FN double cbrt (double) noexcept; - AMREX_VECMATH_FN double cosh (double) noexcept; - AMREX_VECMATH_FN double erf (double) noexcept; - AMREX_VECMATH_FN double erfc (double) noexcept; - AMREX_VECMATH_FN double exp2 (double) noexcept; - AMREX_VECMATH_FN double expm1 (double) noexcept; - AMREX_VECMATH_FN double hypot (double, double) noexcept; - AMREX_VECMATH_FN double log10 (double) noexcept; - AMREX_VECMATH_FN double log1p (double) noexcept; - AMREX_VECMATH_FN double log2 (double) noexcept; - AMREX_VECMATH_FN double sinh (double) noexcept; - AMREX_VECMATH_FN double tan (double) noexcept; - AMREX_VECMATH_FN double tanh (double) noexcept; - AMREX_VECMATH_FN float acosf (float) noexcept; - AMREX_VECMATH_FN float acoshf (float) noexcept; - AMREX_VECMATH_FN float asinf (float) noexcept; - AMREX_VECMATH_FN float asinhf (float) noexcept; - AMREX_VECMATH_FN float atanf (float) noexcept; - AMREX_VECMATH_FN float atan2f (float, float) noexcept; - AMREX_VECMATH_FN float atanhf (float) noexcept; - AMREX_VECMATH_FN float cbrtf (float) noexcept; - AMREX_VECMATH_FN float coshf (float) noexcept; - AMREX_VECMATH_FN float erff (float) noexcept; - AMREX_VECMATH_FN float erfcf (float) noexcept; - AMREX_VECMATH_FN float exp2f (float) noexcept; - AMREX_VECMATH_FN float expm1f (float) noexcept; - AMREX_VECMATH_FN float hypotf (float, float) noexcept; - AMREX_VECMATH_FN float log10f (float) noexcept; - AMREX_VECMATH_FN float log1pf (float) noexcept; - AMREX_VECMATH_FN float log2f (float) noexcept; - AMREX_VECMATH_FN float sinhf (float) noexcept; - AMREX_VECMATH_FN float tanf (float) noexcept; - AMREX_VECMATH_FN float tanhf (float) noexcept; -} -# endif -# undef AMREX_VECMATH_FN -# endif -#endif - namespace amrex { // NOLINT(modernize-concat-nested-namespaces) /// \cond DOXYGEN_IGNORE inline namespace disabled { @@ -262,88 +138,6 @@ namespace detail { *cosx = std::cos(x); #endif } -#ifdef AMREX_USE_SIMD - /** Apply a scalar function to every lane of a SIMD variable - * - * The lane loop has a compile-time constant trip count and a body that is a - * single scalar function call. That is the shape an auto-vectorizer can - * replace with one call into a vector math library (see the note near the - * top of this file). When it does, the lane buffer never reaches memory; - * when it does not, the loop unrolls into scalar calls, which is what the - * SIMD library's own math fallback does as well. - * - * @param x the SIMD variable to transform - * @param lane_fn the scalar function to apply to one lane - * @param simd_fn the SIMD library's own overload of the same function - * @return a SIMD variable holding the function value for every lane - */ - template - AMREX_FORCE_INLINE - T_Simd map_lanes (T_Simd const& x, F_Lane const& lane_fn, F_Simd const& simd_fn) - { -#ifdef AMREX_SIMD_HAS_VECMATH - static_cast(simd_fn); - - using T = typename T_Simd::value_type; - constexpr std::size_t width = T_Simd::size(); - constexpr std::size_t alignment = amrex::simd::stdx::memory_alignment_v; - - alignas(alignment) T lane[width]; - x.copy_to(lane, amrex::simd::stdx::vector_aligned); - - for (std::size_t i = 0; i < width; ++i) { - lane[i] = lane_fn(lane[i]); - } - - T_Simd r; - r.copy_from(lane, amrex::simd::stdx::vector_aligned); - return r; -#else - // Without a vector math library to call, the loop above would only wrap a - // lane buffer around the same scalar calls. Hand the work to the SIMD - // library instead: it evaluates most functions one lane at a time too, but - // it does carry real vectorized implementations for a few of them. - static_cast(lane_fn); - return simd_fn(x); -#endif - } - - /** Apply a scalar function of two arguments to every lane of two SIMD variables - * - * @see map_lanes(T_Simd const&, F const&) - */ - template - AMREX_FORCE_INLINE - T_Simd map_lanes (T_Simd const& x, T_Simd const& y, - F_Lane const& lane_fn, F_Simd const& simd_fn) - { -#ifdef AMREX_SIMD_HAS_VECMATH - static_cast(simd_fn); - - using T = typename T_Simd::value_type; - constexpr std::size_t width = T_Simd::size(); - constexpr std::size_t alignment = amrex::simd::stdx::memory_alignment_v; - - alignas(alignment) T lane_x[width]; - alignas(alignment) T lane_y[width]; - x.copy_to(lane_x, amrex::simd::stdx::vector_aligned); - y.copy_to(lane_y, amrex::simd::stdx::vector_aligned); - - for (std::size_t i = 0; i < width; ++i) { - lane_x[i] = lane_fn(lane_x[i], lane_y[i]); - } - - T_Simd r; - r.copy_from(lane_x, amrex::simd::stdx::vector_aligned); - return r; -#else - // see map_lanes(T_Simd const&, F_Lane const&, F_Simd const&) - static_cast(lane_fn); - return simd_fn(x, y); -#endif - } - -#endif } /// \endcond @@ -412,8 +206,12 @@ std::pair sincospi (float x) * resolves to the SIMD library's own (scalar, element-wise) overload through * argument-dependent lookup instead. * - * The SIMD overloads may lower to vector math library calls, see the note near - * the top of this file for the requirements and the accuracy implications. + * 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 @@ -432,9 +230,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real sin (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::sin(v); }, - [] (auto const& v) { return amrex::simd::stdx::sin(v); }); + return amrex::simd::smath::sin(x); } #endif @@ -454,9 +250,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real cos (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::cos(v); }, - [] (auto const& v) { return amrex::simd::stdx::cos(v); }); + return amrex::simd::smath::cos(x); } #endif @@ -476,9 +270,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real tan (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::tan(v); }, - [] (auto const& v) { return amrex::simd::stdx::tan(v); }); + return amrex::simd::smath::tan(x); } #endif @@ -498,9 +290,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real asin (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::asin(v); }, - [] (auto const& v) { return amrex::simd::stdx::asin(v); }); + return amrex::simd::smath::asin(x); } #endif @@ -520,9 +310,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real acos (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::acos(v); }, - [] (auto const& v) { return amrex::simd::stdx::acos(v); }); + return amrex::simd::smath::acos(x); } #endif @@ -542,9 +330,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real atan (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::atan(v); }, - [] (auto const& v) { return amrex::simd::stdx::atan(v); }); + return amrex::simd::smath::atan(x); } #endif @@ -564,9 +350,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real sinh (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::sinh(v); }, - [] (auto const& v) { return amrex::simd::stdx::sinh(v); }); + return amrex::simd::smath::sinh(x); } #endif @@ -586,9 +370,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real cosh (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::cosh(v); }, - [] (auto const& v) { return amrex::simd::stdx::cosh(v); }); + return amrex::simd::smath::cosh(x); } #endif @@ -608,9 +390,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real tanh (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::tanh(v); }, - [] (auto const& v) { return amrex::simd::stdx::tanh(v); }); + return amrex::simd::smath::tanh(x); } #endif @@ -630,9 +410,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real asinh (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::asinh(v); }, - [] (auto const& v) { return amrex::simd::stdx::asinh(v); }); + return amrex::simd::smath::asinh(x); } #endif @@ -652,9 +430,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real acosh (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::acosh(v); }, - [] (auto const& v) { return amrex::simd::stdx::acosh(v); }); + return amrex::simd::smath::acosh(x); } #endif @@ -674,9 +450,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real atanh (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::atanh(v); }, - [] (auto const& v) { return amrex::simd::stdx::atanh(v); }); + return amrex::simd::smath::atanh(x); } #endif @@ -696,9 +470,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real exp (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::exp(v); }, - [] (auto const& v) { return amrex::simd::stdx::exp(v); }); + return amrex::simd::smath::exp(x); } #endif @@ -718,9 +490,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real exp2 (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::exp2(v); }, - [] (auto const& v) { return amrex::simd::stdx::exp2(v); }); + return amrex::simd::smath::exp2(x); } #endif @@ -740,9 +510,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real expm1 (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::expm1(v); }, - [] (auto const& v) { return amrex::simd::stdx::expm1(v); }); + return amrex::simd::smath::expm1(x); } #endif @@ -762,9 +530,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real log (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::log(v); }, - [] (auto const& v) { return amrex::simd::stdx::log(v); }); + return amrex::simd::smath::log(x); } #endif @@ -784,9 +550,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real log2 (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::log2(v); }, - [] (auto const& v) { return amrex::simd::stdx::log2(v); }); + return amrex::simd::smath::log2(x); } #endif @@ -806,9 +570,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real log10 (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::log10(v); }, - [] (auto const& v) { return amrex::simd::stdx::log10(v); }); + return amrex::simd::smath::log10(x); } #endif @@ -828,9 +590,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real log1p (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::log1p(v); }, - [] (auto const& v) { return amrex::simd::stdx::log1p(v); }); + return amrex::simd::smath::log1p(x); } #endif @@ -850,9 +610,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real cbrt (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::cbrt(v); }, - [] (auto const& v) { return amrex::simd::stdx::cbrt(v); }); + return amrex::simd::smath::cbrt(x); } #endif @@ -872,9 +630,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real erf (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::erf(v); }, - [] (auto const& v) { return amrex::simd::stdx::erf(v); }); + return amrex::simd::smath::erf(x); } #endif @@ -894,9 +650,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real erfc (T_Real const& x) { - return detail::map_lanes(x, - [] (auto v) { return std::erfc(v); }, - [] (auto const& v) { return amrex::simd::stdx::erfc(v); }); + return amrex::simd::smath::erfc(x); } #endif @@ -916,9 +670,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real pow (T_Real const& x, T_Real const& y) { - return detail::map_lanes(x, y, - [] (auto a, auto b) { return std::pow(a, b); }, - [] (auto const& a, auto const& b) { return amrex::simd::stdx::pow(a, b); }); + return amrex::simd::smath::pow(x, y); } #endif @@ -938,9 +690,7 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real atan2 (T_Real const& y, T_Real const& x) { - return detail::map_lanes(y, x, - [] (auto a, auto b) { return std::atan2(a, b); }, - [] (auto const& a, auto const& b) { return amrex::simd::stdx::atan2(a, b); }); + return amrex::simd::smath::atan2(y, x); } #endif @@ -960,9 +710,10 @@ requires (amrex::simd::stdx::is_simd_v) AMREX_FORCE_INLINE T_Real hypot (T_Real const& x, T_Real const& y) { - return detail::map_lanes(x, y, - [] (auto a, auto b) { return std::hypot(a, b); }, - [] (auto const& a, auto const& b) { return amrex::simd::stdx::hypot(a, b); }); + // 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 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/Particles/ParticleReduceSIMD/main.cpp b/Tests/Particles/ParticleReduceSIMD/main.cpp index ef533162ad8..a21b8a2f069 100644 --- a/Tests/Particles/ParticleReduceSIMD/main.cpp +++ b/Tests/Particles/ParticleReduceSIMD/main.cpp @@ -284,20 +284,6 @@ void compare_results (ReduceTupleT const& a, ReduceTupleT const& b, } } -/** Tolerance scale for comparisons that are otherwise required to be bitwise identical - * - * -fno-math-errno, which AMReX_SIMD_VECMATH=ON adds, lets the compiler make - * different inlining and FMA contraction choices per formulation. Two ways of - * writing the same sum then no longer have to agree bit for bit. Keyed off the - * flag rather than off AMREX_SIMD_HAS_VECMATH, because the flag is what perturbs - * the arithmetic, whether or not a vector math library turns out to be reachable. - */ -#ifdef __NO_MATH_ERRNO__ -constexpr Real bitwise_tol_scale = 4.; -#else -constexpr Real bitwise_tol_scale = 0.; -#endif - void correctness_tests (Geometry const& geom, DistributionMapping const& dm, BoxArray const& ba, Shifts const& shifts) { @@ -315,9 +301,9 @@ void correctness_tests (Geometry const& geom, DistributionMapping const& dm, auto const rc = run_variant_c(pc, shifts); // A, B and B2 evaluate identical arithmetic in identical order - compare_results(ra, rb, np, bitwise_tol_scale, "A (SuperParticle) vs B (ptd,i)"); + compare_results(ra, rb, np, Real(0.), "A (SuperParticle) vs B (ptd,i)"); auto const rb2 = run_variant_b2(pc, shifts); - compare_results(rb2, rb, np, bitwise_tol_scale, + compare_results(rb2, rb, np, Real(0.), "B2 (ptd,i by-ref) vs B (ptd,i)"); // C reassociates the sums across SIMD lanes compare_results(rc, rb, np, Real(100.), "C (SIMD) vs B (ptd,i)"); diff --git a/Tests/SIMD/main.cpp b/Tests/SIMD/main.cpp index 42e3fd53124..e98b1b0da37 100644 --- a/Tests/SIMD/main.cpp +++ b/Tests/SIMD/main.cpp @@ -653,10 +653,10 @@ int main (int argc, char* argv[]) T(-2.0), T(2.0)); Print() << "amrex::Math SIMD transcendentals (" -# ifdef AMREX_SIMD_HAS_VECMATH - << "vector math library enabled" +# ifdef VIR_HAVE_SIMD_VECMATH + << "vector math library via the SIMD provider" # else - << "vector math library not available, lane-wise fallback" + << "provider fallback, one call per lane" # endif << ", width " << int(V::size()) << "): " << (err == 0 ? "PASSED" : "FAILED") << "\n"; diff --git a/Tools/CMake/AMReXOptions.cmake b/Tools/CMake/AMReXOptions.cmake index dbfb60f341e..98e52ff61ca 100644 --- a/Tools/CMake/AMReXOptions.cmake +++ b/Tools/CMake/AMReXOptions.cmake @@ -273,11 +273,6 @@ print_option( AMReX_MPI_THREAD_MULTIPLE ) option( AMReX_SIMD "Enable SIMD Primitives" OFF) print_option( AMReX_SIMD ) -cmake_dependent_option( AMReX_SIMD_VECMATH - "Let SIMD math functions call a vector math library (adds -fno-math-errno)" OFF - "AMReX_SIMD" OFF) -print_option( AMReX_SIMD_VECMATH ) - option( AMReX_OMP "Enable OpenMP" OFF) print_option( AMReX_OMP ) diff --git a/Tools/CMake/AMReXParallelBackends.cmake b/Tools/CMake/AMReXParallelBackends.cmake index 954b9d59baf..0f542ced37d 100644 --- a/Tools/CMake/AMReXParallelBackends.cmake +++ b/Tools/CMake/AMReXParallelBackends.cmake @@ -38,41 +38,6 @@ if (AMReX_SIMD) foreach(D IN LISTS AMReX_SPACEDIM) target_link_libraries(amrex_${D}d PUBLIC vir-simd::vir-simd) endforeach() - - # Vector math library for the SIMD math functions in AMReX_Math.H. - # - # A math function may only be vectorized if the compiler does not have to keep - # errno up to date, so -fno-math-errno is needed both here and in every - # downstream translation unit that calls amrex::Math with a SIMD argument. - if (AMReX_SIMD_VECMATH) - include(CheckCXXCompilerFlag) - - check_cxx_compiler_flag("-fno-math-errno" AMReX_HAS_FLAG_NO_MATH_ERRNO) - if (AMReX_HAS_FLAG_NO_MATH_ERRNO) - foreach(D IN LISTS AMReX_SPACEDIM) - target_compile_options(amrex_${D}d - PUBLIC $<$:-fno-math-errno>) - endforeach() - else () - message(WARNING "AMReX_SIMD_VECMATH: ${CMAKE_CXX_COMPILER_ID} does not accept " - "-fno-math-errno. SIMD math functions stay scalar.") - endif () - - # GCC finds the vector variants through the declarations in AMReX_Math.H. - # clang ignores those and uses a built-in mapping table instead, which it - # only consults with -fveclib. - if (CMAKE_CXX_COMPILER_ID MATCHES "Clang" AND CMAKE_SYSTEM_NAME STREQUAL "Linux" - AND AMReX_HAS_FLAG_NO_MATH_ERRNO) - check_cxx_compiler_flag("-fveclib=libmvec" AMReX_HAS_FLAG_VECLIB_LIBMVEC) - if (AMReX_HAS_FLAG_VECLIB_LIBMVEC) - foreach(D IN LISTS AMReX_SPACEDIM) - target_compile_options(amrex_${D}d - PUBLIC $<$:-fveclib=libmvec>) - target_compile_definitions(amrex_${D}d PUBLIC AMREX_SIMD_HAS_VECMATH=1) - endforeach() - endif () - endif () - endif () endif () # From 56346626e3213d40090790182622491729f4408a Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 25 Aug 2026 08:59:22 -0700 Subject: [PATCH 6/8] CI: build the SIMD job against the vir-simd vector math branch The SIMD job installed vir-simd 0.4.4, which predates vir/simd_vecmath.h, so amrex::Math's SIMD transcendentals took the one-call-per-lane path and the interesting half of them was never exercised. Clone the branch that has it instead, with a TODO to go back to a release tarball once it is in one. Co-Authored-By: Claude Opus 5 (1M context) --- .github/workflows/gcc.yml | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) 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 From e660dee29e498246cbb868b8647b510a90f93f35 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 25 Aug 2026 21:29:10 -0700 Subject: [PATCH 7/8] Math: correct what libmvec offers for sincos The comment claimed no vector math library has an auto-vectorizable sincos. glibc's does, on x86-64, and has since 2.22: _ZGVbN2vvv_sincos and its AVX, AVX2 and AVX-512 siblings, for both precisions. Calling it is still the wrong move, but for a different reason than the comment gave. The three v's in that name are the point: the two results are passed as vectors of pointers, so the callee scatters its output a lane at a time. Against two independent vector calls, on AVX2: two calls 1.19 ms fused 1.70 ms 0.70x fused, with the pointer vectors hoisted out of the loop 1.79 ms 0.67x Hoisting them does not help, so it is the scatter and not the setup. Results agree to 1 ULP, so this is a speed argument only. Also note what the pi variants look like, since the two architectures are complementary and neither has what the other has: AArch64 offers sinpi, cospi and tanpi but no fused sincos, x86-64 the reverse, and neither offers sincospi. --- Src/Base/AMReX_Math.H | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/Src/Base/AMReX_Math.H b/Src/Base/AMReX_Math.H index 0e8a114a335..69dc2bec852 100644 --- a/Src/Base/AMReX_Math.H +++ b/Src/Base/AMReX_Math.H @@ -757,9 +757,17 @@ T_Real abs (T_Real const& x) #ifdef AMREX_USE_SIMD /** Return sine and cosine of every SIMD lane * - * Evaluated as two separate lane loops on purpose: no vector math library offers - * an auto-vectorizable sincos, because writing both results through pointers - * stops the vectorizer. Two vector calls still beat one scalar sincos per lane. + * Evaluated as two separate calls on purpose, which is not for want of a fused + * one: glibc's libmvec has offered a vector sincos on x86-64 since 2.22. Its + * vector ABI passes both results as vectors of pointers, though, so the callee + * scatters its output one lane at a time, and that costs more than sharing the + * argument reduction between the two saves -- measured at roughly 1.4x the time + * of two independent vector calls on AVX2. AArch64 does not offer a fused + * sincos at all. + * + * There is no vector sincospi anywhere. AArch64 has separate sinpi and cospi, + * which sincospi below could use once a version of this is worth writing per + * architecture; x86-64 has no pi-suffixed vector function at all. */ template requires (amrex::simd::stdx::is_simd_v) From 13c7e95501fc61cab781fd1112e5aa331dda9873 Mon Sep 17 00:00:00 2001 From: Axel Huebl Date: Tue, 25 Aug 2026 22:11:03 -0700 Subject: [PATCH 8/8] Math: stop routing sincos into the vector math library sincos already existed, and this branch quietly changed what it does: its body went from the SIMD library's own sin and cos to amrex::Math::sin and amrex::Math::cos, which now route. Every existing caller changed behaviour without changing a line, sincospi with it. That is a pessimization where sincos actually gets used. A call into a vector math library is a scheduling barrier: the arithmetic around it can no longer overlap with the transcendental. Where the transcendental dominates the kernel that still pays -- which is the whole point of the explicit amrex::Math::sin and amrex::Math::cos -- but sincos tends to sit in kernels that are mostly other arithmetic. Measured on a quaternion spin rotation of the shape ImpactX's SpinTransport mixin uses, AVX2, one thread pinned to a P-core, best of 15 alternating runs: SIMD library's own sin/cos 3.639 ms routed through libmvec 3.938 ms 8.2% worse Not register pressure, which was the obvious suspect: the loop spills 18 times one way and 19 the other. It is the lost overlap. This matters beyond a microbenchmark because the mixin is shared. Every element's spin push runs it, including elements with no transcendental of their own, so the cost lands on benchmarks that never call sin or cos: ImpactX sees ~8-15% on exactly those. Callers wanting the vector math library still ask for it by name. Verified: amrex::Math::sincos on a simd now emits no libmvec call, while amrex::Math::sin still emits _ZGVdN4v_sin. --- Src/Base/AMReX_Math.H | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/Src/Base/AMReX_Math.H b/Src/Base/AMReX_Math.H index 69dc2bec852..bd620af30f9 100644 --- a/Src/Base/AMReX_Math.H +++ b/Src/Base/AMReX_Math.H @@ -757,17 +757,26 @@ T_Real abs (T_Real const& x) #ifdef AMREX_USE_SIMD /** Return sine and cosine of every SIMD lane * - * Evaluated as two separate calls on purpose, which is not for want of a fused - * one: glibc's libmvec has offered a vector sincos on x86-64 since 2.22. Its - * vector ABI passes both results as vectors of pointers, though, so the callee - * scatters its output one lane at a time, and that costs more than sharing the - * argument reduction between the two saves -- measured at roughly 1.4x the time - * of two independent vector calls on AVX2. AArch64 does not offer a fused - * sincos at all. + * 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. * - * There is no vector sincospi anywhere. AArch64 has separate sinpi and cospi, - * which sincospi below could use once a version of this is worth writing per - * architecture; x86-64 has no pi-suffixed vector function at all. + * 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) @@ -775,8 +784,8 @@ AMREX_FORCE_INLINE std::pair sincos (T_Real const& x) { std::pair r; - r.first = amrex::Math::sin(x); - r.second = amrex::Math::cos(x); + r.first = amrex::simd::stdx::sin(x); + r.second = amrex::simd::stdx::cos(x); return r; }