[WIP] SIMD: Transcendental Math for amrex::Math (sinh, cos, exp, ...) by ax3l · Pull Request #5644 · AMReX-Codes/amrex

ax3l · 2026-08-24T23:31:06Z

Summary

SIMD hardware has instructions for sqrt and abs, but not for the transcendental functions. std::experimental::simd (and so vir::stdx) therefore evaluates sin, sinh, exp and friends one lane at a time, which makes a vectorized kernel that calls them slower than the scalar one it replaced.

Uses mattkretz/vir-simd#53

See

https://eel.is/c++draft/simd#math
https://www.open-std.org/jtc1/sc22/wg21/docs/papers/2024/p1928r15.pdf
and for fun:
https://github.com/NoNaeAbC/std_simd and
for a hit piece https://lucisqr.substack.com/p/c26-shipped-a-simd-library-nobody (obviously vibe-blogged and wrong in some places, e.g., C++26 has vectorized math).

Measured here at width 4, per 2²⁰ evaluations: stdx::sinh takes 8.8 ms where a plain scalar loop takes 8.0 ms. Using SIMD makes it worse.

This PR adds SIMD overloads for the transcendentals in amrex::Math, so that a kernel can be written once and instantiated for scalar, SIMD and GPU, in the same single-source style as ParallelForSIMD and ParticleReduceSIMD:

// T_Real is amrex::ParticleReal in a scalar or GPU build, a SIMD type in a vectorized one
template <typename T_Real>
AMREX_GPU_HOST_DEVICE 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));
    ...
}

Covered: sin cos tan asin acos atan atan2 sinh cosh tanh asinh acosh atanh exp exp2 expm1 log log2 log10 log1p pow sqrt cbrt hypot erf erfc abs sincos sincospi.

Where the speed comes from

AMReX does not implement vector math itself. The SIMD overloads forward to the SIMD provider, through one alias in AMReX_SIMD.H:

#ifdef VIR_HAVE_SIMD_VECMATH
    namespace smath = vir::vecmath;   // a vector math library, e.g. glibc's libmvec
#else
    namespace smath = vir::stdx;      // the provider's own, one call per lane
#endif

so each overload is a one-line forward. The vector math support itself lives in vir-simd (https://github.com/mattkretz/vir-simd), where it belongs: it calls libmvec's entry points directly by their x86-64 vector-function-ABI names. No compiler flags are involved, nothing depends on the auto-vectorizer, and the header is guarded with __has_include(<vir/simd_vecmath.h>), so AMReX still builds against vir-simd releases that do not have it — the overloads then simply forward to the per-lane implementation.

hypot deliberately stays with vir::stdx: SIMD libraries generally implement it with SIMD instructions already, including overflow fixups a vector math library's version skips. sqrt and abs map onto hardware instructions and are likewise handed straight to the provider.

sincos and sincospi stay with vir::stdx too, for a different reason. A call into a vector math library is a scheduling barrier: the surrounding arithmetic can no longer overlap with the transcendental. Where the transcendental dominates a kernel that is worth paying for — which is exactly why the single-result overloads route — but sincos tends to sit in kernels that are mostly other arithmetic, a rotation or a coordinate transform, and there the barrier costs more than the faster transcendental saves. 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):

`amrex::Math::sincos` evaluated with	ms
the SIMD provider's own `sin`/`cos` (this PR)	3.64
routed to the vector math library	3.94

Not register pressure, which is the obvious suspect: the loop spills 18 times one way and 19 the other. It is the lost overlap. This matters more than a microbenchmark suggests, because sincos is usually reached from shared code — ImpactX's spin-transport mixin runs for every element, so routing it moved benchmarks containing no transcendental of their own by 8–15%. A caller whose kernel really is dominated by the transcendental should call amrex::Math::sin and amrex::Math::cos separately.

There is also no fused vector sincos worth calling. glibc's libmvec has offered one on x86-64 since 2.22 (_ZGVdN4vvv_sincos), but its vector ABI returns both results through vectors of pointers, so the callee scatters its output one lane at a time — 1.70 ms against 1.19 ms for two independent vector calls (AVX2, 2²⁰ evaluations), and hoisting the pointer vectors out of the loop does not help (1.79 ms), which pins the cost on the scatter rather than the setup. AArch64 offers no fused sincos at all, and neither architecture offers a vector sincospi.

Note on calling convention

Call these fully qualified, as amrex::Math::sinh(x). An unqualified sinh(x) on a SIMD argument resolves to the SIMD library's own overload through argument-dependent lookup, and neither a using amrex::Math::sinh; declaration nor a using-directive changes that: the library's overload is either equally or more specialized, so it ties (ambiguous) or wins (silently slow). Both outcomes were verified. This is a property of ADL, not of any particular implementation, and is why the amrex::Math layer exists rather than users calling stdx:: directly.

Measurements

i9-12900H (AVX2, width 4), GCC 13.3, -O3 -march=native, one thread pinned to a P-core, N = 2²⁰, best of 5 batches × 15 passes. No compiler flags beyond the usual.

Per function, ms per pass — "stdx" is what an unqualified call resolves to today, "Math" is amrex::Math:

fn	scalar loop	stdx	`amrex::Math`	vs stdx
`sinh`	8.0	8.8	1.2	7.2×
`erf`	7.5	7.8	1.2	6.4×
`tan`	4.4	4.8	0.9	5.7×
`exp`	2.9	3.9	0.7	5.7×
`cbrt`	8.3	8.6	1.6	5.5×
`tanh`	7.1	8.4	1.9	4.5×
`log`	2.6	3.8	1.0	3.8×
`cosh`	3.8	4.9	1.4	3.6×
`atan`	3.9	4.3	1.3	3.4×
`sin`	4.9	1.2	0.7	1.7×
`cos`	5.1	1.2	0.7	1.7×
`sqrt`	1.3	0.65	0.65	1.0×

sin/cos start out closer because libstdc++ is the one SIMD library that carries its own vectorized implementations of exactly those two; everything else it evaluates per lane. sqrt is a hardware instruction and identical everywhere, as expected.

ImpactX-shaped kernel (per particle: two sqrt, sinh, cosh, plus transport arithmetic):

code path	ms
scalar loop	22.3
SIMD, unqualified → `stdx` (status quo)	16.8
SIMD, `amrex::Math` (this PR)	5.9

2.8× over the current SIMD path, 3.7× over the scalar loop.

Accuracy: within 4 ULP of scalar libm where a vector math library answers (glibc documents 4 for libmvec; ≤3 measured across 2¹⁸ points per function), and bit-identical to the provider otherwise. Results are therefore not bit-wise reproducible against a scalar build — noted with a warning in the User's Guide.

Testing

Tests/SIMD compares every SIMD overload against its scalar counterpart over a range, with an 8 ULP tolerance, plus scalar-overload checks that run in every build:

cmake -S . -B build -DCMAKE_BUILD_TYPE=Release -DAMReX_SPACEDIM=3 \
      -DAMReX_MPI=OFF -DAMReX_OMP=OFF -DAMReX_FORTRAN=OFF \
      -DAMReX_SIMD=ON -DAMReX_ENABLE_TESTS=ON -DAMReX_TEST_TYPE=All
cmake --build build -j6 && ctest --test-dir build

Passing with AMReX_SIMD=ON (both with and without a vir-simd that has the vector math header), with AMReX_SIMD=OFF, in single precision (AMReX_PRECISION=SINGLE, where the width 8 float path resolves to _ZGVdN8v_*f), and with a conda toolchain whose sysroot is too old for libmvec, where it forwards per lane. Full ctest green. Clean under the SIMD CI job's warning set and clang-tidy-21.

Dependency

The fast path needs a vir-simd that provides vir/simd_vecmath.h, which is not in a release yet — proposed upstream separately. Until then this PR is a no-op in performance terms and a pure API addition: amrex::Math gains the SIMD overloads, and they forward per lane.

WIP / to do

Wait on the vir-simd side landing upstream, then bump the minimum vir-simd version and drop the __has_include guard.
Move the SIMD math functions into their own header, e.g. AMReX_Math_SIMD.H, included by AMReX_Math.H. A @todo marks the spot.
~~Decide whether sincos should stay two calls.~~ Decided: two calls, and evaluated with the SIMD provider's own sin/cos rather than routed — see "Where the speed comes from". (An earlier revision of this PR routed sincos as well, which silently changed every existing caller and cost 8% where it is used.)

🤖 Generated with Claude Code

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) <noreply@anthropic.com>

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) <noreply@anthropic.com>

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) <noreply@anthropic.com>

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) <noreply@anthropic.com>

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) <noreply@anthropic.com>

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) <noreply@anthropic.com>

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.

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.

## Summary `clang-tidy`'s `modernize-pass-by-value` fires on `detail::gpu_tuple_element`'s constructor and asks for the argument to be taken by value and moved. Following that advice would be a pessimization here. This adds a `NOLINTNEXTLINE` with the reasoning written down, because it is not obvious from the code. ## Why by value is wrong here The forwarding constructor immediately below already takes everything that can be deduced — lvalues and rvalues of `T` alike. The const-reference overload is reached only by arguments it *cannot* deduce, a braced initializer list above all: ```cpp struct P { int a, b; }; P p{1,2}; E<P> a(p); // forwarding constructor E<P> b(P{3,4}); // forwarding constructor E<P> c({5,6}); // this one -- U cannot be deduced from a braced-init-list ``` (verified by instrumenting the two constructors and running the three cases) So taking it by value would add a move for exactly the callers this overload exists to serve, and save a copy for nobody, since rvalues never reach it. Removing the overload instead is not an option either — the braced-init case would stop compiling. ## Why it is showing up now Nothing changed in `AMReX_Tuple.H`. The SIMD CI job runs `clang-tidy` over ccache *misses*, so which headers get linted depends on what was recompiled. A change to a widely included header — #5644 touches `AMReX_Math.H` — pulls `AMReX_Tuple.H` into the linted set and the finding appears. It is pre-existing on `development` and unrelated to that work, which is why it is split out here: this is a two-line comment change that can be reviewed on its own, and #5644 can rebase on it afterwards. ## Testing ``` clang-tidy-21 --config-file=.clang-tidy --header-filter='.*AMReX.*' ``` over `Tests/SIMD/main.cpp` and `Tests/Particles/ParticleReduceSIMD/main.cpp`: `modernize-pass-by-value` findings across all AMReX headers go from 1 to 0, and this was the only one. `AMReX_Tuple.H` still compiles and behaves the same (`amrex::makeTuple` with scalar and aggregate members). 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-authored-by: Claude Opus 5 (1M context) <noreply@anthropic.com>

ax3l added the performance label Aug 24, 2026

ax3l mentioned this pull request Aug 24, 2026

How to Vectorize C++20 Transcendental Std Math Functions? mattkretz/vir-simd#52

Open

ax3l and others added 4 commits August 24, 2026 18:06

ax3l changed the title ~~[WIP] SIMD: Vectorized Transcendental Math (sinh, cos, exp, ...)~~ [WIP] SIMD: Transcendental Math for amrex::Math (sinh, cos, exp, ...) Aug 25, 2026

ax3l mentioned this pull request Aug 25, 2026

[WIP] SIMD: Vectorized Transcendentals in Element Pushes BLAST-ImpactX/impactx#1632

Open

3 tasks

ax3l force-pushed the topic-simd-vecmath branch from d653ddb to 5634662 Compare August 25, 2026 20:48

ax3l mentioned this pull request Aug 25, 2026

Tuple: Silence modernize-pass-by-value #5646

Merged

ax3l added 2 commits August 25, 2026 21:29

Navigation Menu

Uh oh!

[WIP] SIMD: Transcendental Math for amrex::Math (sinh, cos, exp, ...) - #5644

[WIP] SIMD: Transcendental Math for amrex::Math (sinh, cos, exp, ...)#5644
ax3l wants to merge 8 commits into
AMReX-Codes:developmentfrom
ax3l:topic-simd-vecmath

ax3l commented Aug 24, 2026 •

edited

Loading

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Uh oh!

Conversation

ax3l commented Aug 24, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

Where the speed comes from

Note on calling convention

Measurements

Testing

Dependency

WIP / to do

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

ax3l commented Aug 24, 2026 •

edited

Loading