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

.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

Docs/sphinx_documentation/source/Basics.rst

@@ -3066,6 +3066,76 @@ of signatures (including the :cpp:`CompileTimeOptions` variants) and
 is identical to :cpp:`ParallelFor` on GPU, but does not add the SIMD
 pragma on CPU.
 
+.. _sec:basics:simdmath:
+
+SIMD Math Functions
+===================
+
+When AMReX is built with SIMD support (CMake option ``AMReX_SIMD=ON``, see
+:ref:`sec:build:cmake`), the math functions in ``AMReX_Math.H`` accept SIMD
+variables in addition to :cpp:`float` and :cpp:`double`. A kernel can therefore
+be written once and instantiated for both, which is the same single-source style
+used by :cpp:`amrex::ParallelForSIMD` and :cpp:`amrex::ParticleReduceSIMD`:
+
+.. highlight:: c++
+
+::
+
+    #include <AMReX_Math.H>
+
+    // T_Real is amrex::ParticleReal in a scalar build and a SIMD type in a
+    // vectorized one
+    template <typename T_Real>
+    AMREX_FORCE_INLINE
+    void focus (T_Real & AMREX_RESTRICT y, T_Real & AMREX_RESTRICT py,
+                T_Real const & AMREX_RESTRICT omega, amrex::Real ds)
+    {
+        T_Real const ch = amrex::Math::cosh(omega * T_Real(ds));
+        T_Real const sh = amrex::Math::sinh(omega * T_Real(ds));
+
+        T_Real const y0 = y;
+        y  = ch * y0 + sh / omega * py;
+        py = omega * sh * y0 + ch * py;
+    }
+
+The following functions have SIMD overloads: :cpp:`sin`, :cpp:`cos`,
+:cpp:`tan`, :cpp:`asin`, :cpp:`acos`, :cpp:`atan`, :cpp:`atan2`, :cpp:`sinh`,
+:cpp:`cosh`, :cpp:`tanh`, :cpp:`asinh`, :cpp:`acosh`, :cpp:`atanh`,
+:cpp:`exp`, :cpp:`exp2`, :cpp:`expm1`, :cpp:`log`, :cpp:`log2`, :cpp:`log10`,
+:cpp:`log1p`, :cpp:`pow`, :cpp:`sqrt`, :cpp:`cbrt`, :cpp:`hypot`, :cpp:`erf`,
+:cpp:`erfc`, :cpp:`abs`, :cpp:`sincos` and :cpp:`sincospi`.
+
+.. note::
+
+   Call these functions **fully qualified**, as :cpp:`amrex::Math::sinh(x)`. An
+   unqualified :cpp:`sinh(x)` on a SIMD argument resolves to the SIMD library's
+   own overload through argument-dependent lookup, which evaluates the function
+   one lane at a time. A :cpp:`using amrex::Math::sinh;` declaration does not
+   change that, because the library overload is the more specialized candidate.
+
+Vectorizing the transcendentals
+-------------------------------
+
+SIMD hardware has instructions for :cpp:`sqrt` and :cpp:`abs`, but not for the
+transcendental functions. Those have to be evaluated by a vector math library,
+such as glibc's ``libmvec``, which computes a whole SIMD register worth of
+results per call. AMReX does not implement that itself: the SIMD overloads
+forward to whatever the SIMD provider offers, which is a vector math library
+call where the provider has one, and one scalar call per lane otherwise.
+
+Nothing has to be configured for this, and no compiler flags are involved.
+Whether the fast path is available depends on the provider and the platform,
+and both cases give correct results.
+
+.. warning::
+
+   Vector math libraries trade accuracy for speed: glibc's ``libmvec``
+   documents a maximum error of 4 ULP, where its scalar routines stay below
+   1 ULP. Where they are used, results differ slightly from a scalar build and
+   are not bit-wise reproducible against one. Consult the SIMD provider's
+   documentation if your application needs the accuracy of the scalar routines
+   or checks ``errno`` after math calls.
+
 Ghost Cells
 ===========
 

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