From c449fe0cb71a07ec3c50ac8fa75e3ad3884ec5a1 Mon Sep 17 00:00:00 2001 From: John Travers Date: Sat, 24 Oct 2026 13:12:22 +0100 Subject: [PATCH 1/4] Enable NEON FFTW and tune persistent Raman plan selection --- docs/guide/performance/building-native.md | 40 +++++++++++++-- docs/guide/performance/tuning.md | 11 ++++- native/ci/build-fftw.sh | 5 +- native/include/solprop_native/raman_fft.hpp | 30 ++++++----- native/src/bindings.cpp | 4 +- native/tests/native/test_fft.cpp | 55 +++++++++++++++++++-- native/tests/test_raman_tail.py | 43 ++++++++++++---- 7 files changed, 154 insertions(+), 34 deletions(-) diff --git a/docs/guide/performance/building-native.md b/docs/guide/performance/building-native.md index 961277fe6..ab99918b5 100644 --- a/docs/guide/performance/building-native.md +++ b/docs/guide/performance/building-native.md @@ -77,6 +77,38 @@ export SDKROOT=/Library/Developer/CommandLineTools/SDKs/MacOSX.sdk export CMAKE_ARGS="-DCMAKE_OSX_SYSROOT=$SDKROOT" ``` +For an Apple silicon performance build, the pinned dependency recipe enables FFTW +3.3.11's NEON codelets in **both precisions**. Compiling scalar code with `-mcpu` +does not enable those codelets. The recipe also retains the Arm cycle counter needed +by measured planning, the Apple M1 distribution baseline and the macOS 14 deployment +target. Its FFTW threading libraries and the native extension share one source-built +LLVM OpenMP runtime; Accelerate remains the BLAS provider. + +From the repository root, use a fresh, owned dependency prefix: + +```sh +native_deps="$PWD/native/build/neon-dependencies" +MACOSX_DEPLOYMENT_TARGET=14.0 native/ci/build-deps.sh "$native_deps" +export CMAKE_ARGS="-DFFTW_ROOT=$native_deps -DSOLPROP_NATIVE_LIBOMP_LIBRARY=$native_deps/lib/libomp.dylib -DSOLPROP_NATIVE_LIBOMP_INCLUDE_DIR=$native_deps/include -DSOLPROP_NATIVE_BUNDLED_LICENCES=$native_deps/licences" +uv sync --extra native --reinstall-package solprop-native +uv run --no-sync python -m solprop_native +``` + +Check the reported `fft_library` and `openmp_library` before benchmarking. The +[pinned upstream configuration](https://github.com/FFTW/fftw3/blob/fftw-3.3.11/configure.ac) +supports `--enable-neon` for double and single precision on AArch64. The installed +libraries identify their codelet family in `fftw_version` and `fftwf_version`, which +both end in `-neon` for this recipe. Keep these prefixes separate from a comparison's +baseline libraries. A pyFFTW build used in the same process must also share this +runtime; see [OpenMP runtimes](#openmp-runtimes-and-pyfftw) below. + +For a distributable wheel, `native/ci/build-wheel-macos.sh OUT_DIR DEPS_PREFIX` +builds the distribution profile, repairs its library references with delocate, and +includes the FFTW and OpenMP licences. Use a fresh dependency prefix when changing +the recipe: that helper intentionally reuses an existing prefix. Validate the repaired +wheel in a clean environment with `python -I -m solprop_native`, its linked libraries, +the native tests and the import-order/coexistence checks before distributing it. + **Linux**, on Debian or Ubuntu. Ubuntu 24.04's default compiler is GCC 13, which the build refuses; install GCC 14 from the distribution's archive: @@ -310,7 +342,7 @@ The options are read by `native/CMakeLists.txt` and `native/cmake/HostLibraries. |---|---|---| | `SOLPROP_NATIVE_FFT_PROVIDER` | `FFTW` | `FFTW`, or `MKL_FFTW` (MKL's FFTW3 interface, Linux). Host builds only | | `SOLPROP_NATIVE_FFTW_THREADING` | `OPENMP` | `OPENMP` (`libfftw3_omp`) or `PTHREADS` (`libfftw3_threads`); IntelLLVM with FFTW requires `PTHREADS` | -| `SOLPROP_NATIVE_RAMAN_FFT_PROVIDER` | `ACCELERATE` on macOS CPU builds, else `KOKKOSFFT` | the Raman convolution transforms through Accelerate's vDSP (contiguous power-of-two lengths; short serial transforms use FFTW) or through Kokkos-FFT. `ACCELERATE` is macOS CPU only | +| `SOLPROP_NATIVE_RAMAN_FFT_PROVIDER` | `ACCELERATE` on macOS CPU builds, else `KOKKOSFFT` | construction-time Raman selection: with NEON FFTW, vDSP handles short parallel power-of-two rows (up to 16,384 samples); serial and longer rows use Kokkos-FFT. Without NEON, only short serial rows use Kokkos-FFT. `KOKKOSFFT` forces that provider; `ACCELERATE` is macOS CPU only | | `FFTW_ROOT` | the Homebrew prefix on macOS | a non-default FFTW installation | **Development and diagnostics** @@ -471,7 +503,8 @@ A process can safely load only one OpenMP runtime. Two different runtimes (LLVM' `libomp` and GNU's `libgomp`, or two copies of one) each start their own thread pools, oversubscribe the cores, and can crash when the second initialises. A CPU build therefore links every native part — Kokkos, FFTW's threaded library, an OpenMP BLAS — -against one runtime: on macOS, Homebrew's `libomp`, which Homebrew's FFTW also uses. +against one runtime: on macOS, either Homebrew's `libomp` with its FFTW, or the one +runtime in the owned dependency prefix built above. solprop's `numba_cpu` backend uses [pyFFTW](https://github.com/pyFFTW/pyFFTW), whose published wheels bundle their own OpenMP runtimes, so running both backends in one @@ -506,7 +539,8 @@ CC=gcc CXX=g++ PYFFTW_USE_PTHREADS=1 PYFFTW_FFTW_PREFIX=/usr \ ``` For an FFTW installed elsewhere, set `PYFFTW_FFTW_PREFIX` to its prefix; pyFFTW needs the -single-precision FFTW libraries as well as the double-precision ones. With AOCL's FFTW, +single-precision FFTW libraries as well as the double-precision ones. On macOS also +change both library paths in `LDFLAGS` to that prefix's OpenMP library directory. With AOCL's FFTW, build pyFFTW with AOCC and without `PYFFTW_USE_PTHREADS`, so that it uses AOCL's OpenMP interface; do not combine a standard FFTW pthreads library with AOCL's FFTW core. diff --git a/docs/guide/performance/tuning.md b/docs/guide/performance/tuning.md index 3007453a3..d6c82d671 100644 --- a/docs/guide/performance/tuning.md +++ b/docs/guide/performance/tuning.md @@ -34,7 +34,7 @@ The backend decides which pools exist at all (see [Choosing a backend](backends. |---|---|---| | `numpy_ref` | NumPy's BLAS only | the reference; single-threaded apart from matrix products | | `numba_cpu` | Numba kernels, FFTW, BLAS (and FFTW for the transverse transforms of free-space runs) | sizes each pool from the problem | -| native CPU (`native`) | OpenMP kernels, FFTW (modal and spatial transforms), the Raman convolution transforms (Accelerate's vDSP on macOS, Kokkos-FFT over FFTW elsewhere), a vendor BLAS | uses every logical CPU unless told otherwise | +| native CPU (`native`) | OpenMP kernels, FFTW (modal and spatial transforms), the Raman convolution transforms (Kokkos-FFT over FFTW, with selected macOS shapes on Accelerate's vDSP), a vendor BLAS | uses every logical CPU unless told otherwise | | native GPU (`native`, CUDA or HIP build) | the device | CPU pools are irrelevant; one device per process | `backend="auto"` chooses the native backend when `solprop-native` is installed and @@ -240,6 +240,15 @@ Free-space transforms are planned at `SOLPROP_FFTW_PLAN`'s effort short runs and tests. Answers at different thread counts may differ by rounding; solprop holds them to its rounding bounds, not to identical bits ([SEP 25](../../dev/seps/archive/sep-0025-rounding-level-reproducibility.md)). +Contiguous modal and Raman FFTW plans default to measured planning at one thread and +estimate planning at wider widths; the explicit effort overrides both, including +long Raman rows and partial chunks. Measured planning can add substantial construction +time even when the warmed propagation is faster. Use `SOLPROP_FFTW_PLAN=estimate` +when construction dominates a short run. On macOS with NEON FFTW, serial and long +Raman rows use Kokkos-FFT, while parallel power-of-two rows of at most 16,384 samples +retain vDSP. The [build guide](building-native.md) describes the coherent dependency +recipe and the fallback for a library without NEON codelets. + `OMP_NUM_THREADS` and `KOKKOS_NUM_THREADS` do **not** size the Kokkos pool — the backend initialises Kokkos with an explicit thread count, which takes precedence; set `SOLPROP_NATIVE_THREADS` instead. (Kokkos prints a warning at start-up if diff --git a/native/ci/build-fftw.sh b/native/ci/build-fftw.sh index ea7f7a7a7..3e152e7a7 100755 --- a/native/ci/build-fftw.sh +++ b/native/ci/build-fftw.sh @@ -14,7 +14,8 @@ # x86-64 FFTW's SSE2 and AVX codelets are chosen at run time (AVX2 and AVX-512 codelets were not a # uniform gain when measured, so they are left out). On Arm the cycle counter is # enabled, so that measured planning (FFTW_MEASURE) works; without one FFTW plans by -# estimate only. +# estimate only. Apple silicon also enables FFTW's explicit NEON codelets in both +# precisions; compiler auto-vectorisation does not enable these implementations. # # OpenMP: on Linux the compiler's own runtime (GCC's libgomp in the manylinux image); # on macOS the libomp of LIBOMP_PREFIX (native/ci/build-libomp.sh), the runtime the @@ -49,7 +50,7 @@ case "$(uname -s)/$(uname -m)" in : "${MACOSX_DEPLOYMENT_TARGET:?set MACOSX_DEPLOYMENT_TARGET}" : "${LIBOMP_PREFIX:?set LIBOMP_PREFIX to the libomp built by build-libomp.sh}" cflags="-O3 -mcpu=apple-m1 -mmacosx-version-min=$MACOSX_DEPLOYMENT_TARGET" - simd=(--enable-armv8-cntvct-el0) + simd=(--enable-armv8-cntvct-el0 --enable-neon) # The extension's compiler, Homebrew's Clang (native/pyproject.toml), which takes # -fopenmp and links the libomp found first on the library path: this one. # Dead-stripping keeps libomp out of the libraries that do not use it. diff --git a/native/include/solprop_native/raman_fft.hpp b/native/include/solprop_native/raman_fft.hpp index af7baba3d..f9970c3c2 100644 --- a/native/include/solprop_native/raman_fft.hpp +++ b/native/include/solprop_native/raman_fft.hpp @@ -13,6 +13,7 @@ #include #include #include +#include #include namespace solprop_native { @@ -44,17 +45,22 @@ template class RamanFftPlan { RamanFftPlan(const Execution& execution, const Input& input, const Output& output, const KokkosFFT::Direction direction, const int axis) : execution_(execution) { - constexpr bool forward = - std::is_same_v; - const auto samples = forward ? input.extent(1) : output.extent(1); +#if !defined(KOKKOS_ENABLE_CUDA) && !defined(KOKKOS_ENABLE_HIP) const bool contiguous = (axis == -1 || axis == 1) && input.stride(1) == 1 && output.stride(1) == 1; - // At one thread the short Raman rows favour measured FFTW; vDSP retains - // the large-field and parallel workloads for which it was selected. - const bool short_serial = - contiguous && execution.concurrency() == 1 && samples <= 16384; #if defined(SOLPROP_NATIVE_RAMAN_ACCELERATE) - if (!short_serial && (axis == -1 || axis == 1) && input.stride(1) == 1 && + constexpr bool forward = + std::is_same_v; + const auto samples = forward ? input.extent(1) : output.extent(1); + // NEON FFTW wins the long rows, but short threaded plans can still be slower + // than vDSP. Preserve the previous policy for libraries without NEON codelets; + // a local build may still resolve such an FFTW even on an Arm CPU. + const bool neon = + std::string_view(fftw_version).find("-neon") != std::string_view::npos; + const bool prefer_fftw = + (execution.concurrency() == 1 && samples <= 16384) || + (neon && (execution.concurrency() == 1 || samples > 16384)); + if (!prefer_fftw && (axis == -1 || axis == 1) && input.stride(1) == 1 && output.stride(1) == 1 && samples >= 8 && std::has_single_bit(samples)) { const auto frequencies = forward ? output.extent(1) : input.extent(1); if (input.extent(0) != output.extent(0) || frequencies != samples / 2 + 1) @@ -73,11 +79,11 @@ template class RamanFftPlan { return; } #endif -#if !defined(KOKKOS_ENABLE_CUDA) && !defined(KOKKOS_ENABLE_HIP) - if (short_serial) + // Every FFTW fallback must honour the planner effort, including the long + // doubled rows used when the Accelerate provider is disabled. The shared + // helper retains serial-only measured planning unless explicitly overridden. + if (contiguous) seed_row_fft_wisdom(execution, input, output, direction); -#else - (void)short_serial; #endif fallback_ = std::make_unique(execution, input, output, direction, axis); } diff --git a/native/src/bindings.cpp b/native/src/bindings.cpp index 3dab730cc..d1b890387 100644 --- a/native/src/bindings.cpp +++ b/native/src/bindings.cpp @@ -1319,8 +1319,8 @@ nb::dict build_record() { #else result["fft_planner"] = std::string("Kokkos-FFT persistent Plan / ") + SOLPROP_NATIVE_HOST_FFT_PLANNER + - "; general modal defaults: serial measure, threaded " - "estimate; short serial Raman: measure; " + "; contiguous modal/Raman defaults: serial measure, " + "threaded estimate; " "SOLPROP_FFTW_PLAN overrides"; #endif result["layout"] = "LayoutRight"; diff --git a/native/tests/native/test_fft.cpp b/native/tests/native/test_fft.cpp index 41f0a5ccc..c8fc9e508 100644 --- a/native/tests/native/test_fft.cpp +++ b/native/tests/native/test_fft.cpp @@ -10,6 +10,7 @@ #include #include +#include #include #include #include @@ -175,6 +176,42 @@ TEST(FftProvider, RowWisdomUsesSerialDefaultAndHonorsExplicitEffort) { EXPECT_EQ(static_cast(witness), measured); } } + +TEST(RamanFftProvider, LongFallbackRowsHonorPlannerEffort) { + using Wisdom = std::unique_ptr; + const Wisdom supported(fftw_export_wisdom_to_string(), fftw_free); + if (!supported) + GTEST_SKIP() << "The linked FFT provider does not export FFTW wisdom"; + using Plan = + std::unique_ptr, decltype(&fftw_destroy_plan)>; + constexpr int rows = 1; + const int samples = 32768; + const int frequencies = samples / 2 + 1; + const RealRows input("raman_fft.wisdom_input", rows, samples); + const ComplexRows output("raman_fft.wisdom_output", rows, frequencies); + const char* requested = std::getenv("SOLPROP_FFTW_PLAN"); + for (const int threads : {1, native::ExecutionSpace().concurrency()}) { + SCOPED_TRACE(threads); + const auto execution = native::execution_instance(threads); + fftw_forget_wisdom(); + native::RamanFftPlan forward(execution, input, output, + KokkosFFT::Direction::forward, -1); + if (std::string_view(forward.provider()) != "KokkosFFT") + GTEST_SKIP() << "This library uses the vDSP provider for long Raman rows"; + fftw_plan_with_nthreads(execution.concurrency()); + // ESTIMATE alone cannot create MEASURE wisdom. This witnesses the Raman + // wrapper's seeding for a length beyond the former short-row cutoff. + const Plan witness(fftw_plan_many_dft_r2c( + 1, &samples, rows, input.data(), &samples, 1, samples, + reinterpret_cast(output.data()), &frequencies, + 1, frequencies, FFTW_WISDOM_ONLY | FFTW_MEASURE), + fftw_destroy_plan); + const bool measured = requested == nullptr + ? execution.concurrency() == 1 + : std::string_view(requested) != "estimate"; + EXPECT_EQ(static_cast(witness), measured); + } +} #endif TEST(FftProvider, PartialChunkPlanLeavesUnusedRowsUntouched) { @@ -205,10 +242,15 @@ TEST(FftProvider, PartialChunkPlanLeavesUnusedRowsUntouched) { namespace { TEST(RamanFftProvider, PreservesRawNormalizationInputAndPartialRows) { - constexpr int rows = 7, tail = 3; + // Include the long full/tail geometries whose provider choice differs from + // the short serial rows. Padding remains an exact storage contract. + constexpr std::array shapes{std::array{15, 7, 3}, std::array{32, 7, 3}, + std::array{16384, 7, 3}, std::array{32768, 7, 3}, + std::array{65536, 21, 12}, std::array{131072, 3, 1}}; for (const int threads : {1, native::ExecutionSpace().concurrency()}) { - for (const int length : {15, 32, 16384, 32768}) { + for (const auto& [length, rows, tail] : shapes) { SCOPED_TRACE(length); + SCOPED_TRACE(threads); const RealRows input("raman_fft.input", rows, length); const RealRows restored("raman_fft.restored", rows, length); const ComplexRows output("raman_fft.output", rows, length / 2 + 1); @@ -242,7 +284,7 @@ TEST(RamanFftProvider, PreservesRawNormalizationInputAndPartialRows) { for (int sample = 0; sample < length; ++sample) EXPECT_NEAR(roundtrip(row, sample), row < tail ? length * input_host(row, sample) : -5.0, - fft_bound(length, length)); + row < tail ? fft_bound(length, length) : 0.0); for (int sample = 0; sample <= length / 2; ++sample) { const auto expected = row < tail @@ -253,11 +295,14 @@ TEST(RamanFftProvider, PreservesRawNormalizationInputAndPartialRows) { 0) : Complex(-7.0, 11.0); EXPECT_LE(Kokkos::abs(spectrum(row, sample) - expected), - fft_bound(length, length)); + row < tail ? fft_bound(length, length) : 0.0); } } #if defined(SOLPROP_NATIVE_RAMAN_ACCELERATE) - const bool accelerated = length != 15 && (threads > 1 || length > 16384); + const bool neon = + std::string_view(fftw_version).find("-neon") != std::string_view::npos; + const bool accelerated = length != 15 && (neon ? threads > 1 && length <= 16384 + : threads > 1 || length > 16384); EXPECT_STREQ(forward.provider(), accelerated ? "Accelerate_vDSP" : "KokkosFFT"); EXPECT_EQ(forward.tasks(execution.concurrency()), accelerated ? std::min(tail, execution.concurrency()) diff --git a/native/tests/test_raman_tail.py b/native/tests/test_raman_tail.py index 1a96756ff..eb60943ee 100644 --- a/native/tests/test_raman_tail.py +++ b/native/tests/test_raman_tail.py @@ -6,6 +6,7 @@ from __future__ import annotations +import ctypes import dataclasses import math from typing import cast @@ -27,6 +28,20 @@ from solprop_native import NativeBackend, build_info +def _accelerated(samples: int, threads: int) -> bool: + """Check the actual FFTW codelet family used by the host provider policy.""" + info = build_info() + if info.raman_fft_option != "ACCELERATE": + return False + library = ctypes.CDLL(info.fft_library) + version = ctypes.string_at( + ctypes.addressof(ctypes.c_char.in_dll(library, "fftw_version")) + ) + if b"-neon" in version: + return threads > 1 and samples <= 16384 + return threads > 1 or samples > 16384 + + def _bound(samples: int) -> float: """Bound eight fixed steps against another run of them, relative to the peak. @@ -107,15 +122,9 @@ def test_raman_partial_chunk_matches_reference_and_repeats( tail = sizing[ f"raman_tail_fft.raman.{(7 % chunk) * 2 * grid.n_to}" ]["threads"] - accelerated = info.raman_fft_option == "ACCELERATE" and ( - full > 1 or 2 * grid.n_to > 16384 - ) - tail_accelerated = info.raman_fft_option == "ACCELERATE" and ( - tail > 1 or 2 * grid.n_to > 16384 - ) - kernel_accelerated = info.raman_fft_option == "ACCELERATE" and ( - kernel > 1 or 2 * grid.n_to > 16384 - ) + accelerated = _accelerated(2 * grid.n_to, full) + tail_accelerated = _accelerated(2 * grid.n_to, tail) + kernel_accelerated = _accelerated(2 * grid.n_to, kernel) assert metadata["raman_fft_provider"] == ( "Accelerate_vDSP" if accelerated else "KokkosFFT" ) @@ -128,6 +137,22 @@ def test_raman_partial_chunk_matches_reference_and_repeats( } ledger = backend.memory_ledger() assert ledger["node_raman_drive"] == chunk * 2 * grid.n_to * 8 + # Planning another full/tail shape must leave the first state's + # live plans usable, even when both consume shared FFTW wisdom. + for peer_chunk in (1, 4): + peer_problem = dataclasses.replace( + problem, + transform=dataclasses.replace( + problem.transform, node_chunk=peer_chunk + ), + ) + peer = NativeBackend(peer_problem, grid=grid) + try: + peer.set_initial(problem.initial.field0, None) + peer.propagate_to(grid.length / 8) + assert np.isfinite(peer.observe().require_field()).all() + finally: + peer.close() # A repeat of the same prescribed march: rounding only, within the # equivalence bound below (solprop SEP 25). backend.set_initial(problem.initial.field0, None) From 9e644bb9cfe3bc86a129239c744720cf9758bfa0 Mon Sep 17 00:00:00 2001 From: John Travers Date: Sat, 24 Oct 2026 13:12:22 +0100 Subject: [PATCH 2/4] Defer embedded modal quadrature work until it is needed --- docs/dev/backends.md | 18 ++++- native/src/general.cpp | 122 ++++++++++++++++++++++++------- native/tests/test_band_plasma.py | 77 +++++++++++++++++++ 3 files changed, 189 insertions(+), 28 deletions(-) diff --git a/docs/dev/backends.md b/docs/dev/backends.md index 0c6d07145..3ec759eb6 100644 --- a/docs/dev/backends.md +++ b/docs/dev/backends.md @@ -1310,13 +1310,27 @@ at every width; explicit `measure` and experimental `patient` enable timed plann at any width. Estimate can still consume matching wisdom already present in the process; the backend does not clear shared FFTW wisdom. Planning is destructive and precedes state initialization. -Short serial Raman transforms (at most 16384 doubled-axis samples) likewise use -measured FFTW; larger or parallel macOS Raman transforms retain Accelerate. +Contiguous Raman transforms use the same planning policy, including long doubled +rows and partial chunks. On macOS, NEON FFTW handles serial and long Raman rows; +short parallel power-of-two rows (at most 16384 samples) retain Accelerate. Selection +uses the linked FFTW's codelet family and the plan's dimensions and thread width at +construction. A local FFTW without NEON retains the previous short-serial FFTW rule. The full-chunk Raman provider and task counts, planning seconds, and plan heap growth remain in execution metadata. General-modal tail views add no field buffers. These choices preserve explicit transform scaling and the device paths. See the multimode CPU measurements in [benchmarks.md](benchmarks.md#multimode-cpu-step-costs). +For ordinary CPU DOPRI5 modal trials with an embedded transverse rule, only the +final source projects the coarse rule. Earlier stages use fine-only products, +including the fine rows of existing stacked projection scratch. The final coarse +time buffers remain live until accepted-step scoring performs their FFTs and, for +held field-plasma streams, their spectral assembly. Private readiness flags are +invalidated by a new trial, a standalone source evaluation, reset and closure; +accepted diagnostic values and peaks have separate storage and survive rejection. +The validated seven-stage FSAL method data determine eligibility. Adams, low-storage, +general-modal, dispersive and envelope-plasma paths keep their existing scheduling. +No seam fields, persistent field buffers or arithmetic protections change. + Scalar and declared tensor Raman are native too. Persistent doubled-grid batched r2c/c2r plans transform one or three real histories per node and a separate r2c plan transforms one or two effective kernels. The kernel diff --git a/native/src/general.cpp b/native/src/general.cpp index 93100f724..31dd23534 100644 --- a/native/src/general.cpp +++ b/native/src/general.cpp @@ -1820,7 +1820,20 @@ class GeneralWorkspace { } void held_response(std::size_t stage, std::size_t first_node, std::size_t node_count, const TransverseChunk* chunk = nullptr); - void held_transform(); + void held_transform(bool coarse_only = false); + void finish_quadrature(); + bool deferred_quadrature_eligible() const noexcept { +#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP) + return false; +#else + // The validated seven-stage FSAL tableau identifies the ordinary DOPRI5 + // route without another seam field. Other consumers retain eager spectra. + return problem_.first_same_as_last && problem_.n_stages == 7 && + !problem_.low_storage && !problem_.separable_geometry() && !general_ && + !flow_ && !dispersive_ && !problem_.envelope_plasma && !held_radial() && + !problem_.projection_coarse_matrix.empty(); +#endif + } void held_radial_projection(std::size_t first_node, std::size_t node_count); void envelope_cycle(std::size_t channel, std::size_t node_count); void envelope_population(std::size_t channel, std::size_t stage, @@ -1904,6 +1917,12 @@ class GeneralWorkspace { bool has_recycled_{}; std::size_t rk_source_evaluations_{}; bool pending_fsal_{}; + // Only the final trial source supplies the accepted quadrature score. Its + // coarse time buffers survive until scoring, or until a new source replaces + // them. The stored accepted score and peak are independent of these flags. + bool defer_quadrature_{}; + bool project_coarse_{true}; + bool coarse_pending_{}; // Set by a trial for a stage whose finish may run in the next stage's entry kernel; // `finish_pending_` names the stage's finish still to run there. bool defer_finish_{}; @@ -3755,6 +3774,7 @@ void GeneralWorkspace::prepare_trial(const double dz, const TrialScalars& scalar if (!initialized_) throw std::logic_error("no initial field has been loaded"); require_trial_scalars(scalars); + coarse_pending_ = false; upload_trial(scalars); if (flow_) flow_->set_trial(scalars.flow_section, scalars.flow_weights); @@ -3769,23 +3789,30 @@ void GeneralWorkspace::prepare_trial(const double dz, const TrialScalars& scalar // 6.7) unless a coupled flow carries the stage back first, or the stage is the last. finish_pending_ = false; const bool deferrable = !flow_ && !problem_.separable_geometry(); + defer_quadrature_ = deferred_quadrature_eligible(); for (std::size_t stage = first; stage < problem_.n_stages; ++stage) { - if (problem_.raman) - prepare_raman_stage( - stage, scalars.raman_rates.data() + stage * problem_.n_raman_terms, - problem_.n_raman_components > 1 - ? scalars.raman_densities.data() + stage * problem_.n_raman_components - : nullptr); - defer_finish_ = deferrable && stage + 1 < problem_.n_stages; try { + if (problem_.raman) + prepare_raman_stage( + stage, scalars.raman_rates.data() + stage * problem_.n_raman_terms, + problem_.n_raman_components > 1 + ? scalars.raman_densities.data() + stage * problem_.n_raman_components + : nullptr); + defer_finish_ = deferrable && stage + 1 < problem_.n_stages; + project_coarse_ = !defer_quadrature_ || stage + 1 == problem_.n_stages; evaluate_stage(stage, dz); } catch (...) { defer_finish_ = false; finish_pending_ = false; + defer_quadrature_ = false; + project_coarse_ = true; + coarse_pending_ = false; throw; } } defer_finish_ = false; + defer_quadrature_ = false; + project_coarse_ = true; if (finish_pending_) throw std::logic_error("a deferred stage finish was never run"); stages_ready_ = true; @@ -7019,6 +7046,8 @@ void GeneralWorkspace::evaluate_stage(const std::size_t stage, const double dz) } void GeneralWorkspace::evaluate_stage_source(const std::size_t stage, const double dz) { + if (!defer_quadrature_) + coarse_pending_ = false; if (problem_.separable_geometry()) { evaluate_separable_stage(stage, dz); return; @@ -7109,13 +7138,13 @@ bool GeneralWorkspace::modal_core(const std::size_t stage) { if (problem_.family == FieldFamily::envelope) { KokkosFFT::execute(*forward_complex_plan_, projected_complex_, projected_spectrum_, KokkosFFT::Normalization::none); - if (forward_coarse_complex_plan_) + if (forward_coarse_complex_plan_ && !defer_quadrature_) KokkosFFT::execute(*forward_coarse_complex_plan_, projected_coarse_complex_, projected_coarse_spectrum_, KokkosFFT::Normalization::none); } else { KokkosFFT::execute(*forward_real_plan_, projected_real_, projected_spectrum_, KokkosFFT::Normalization::none); - if (forward_coarse_real_plan_) + if (forward_coarse_real_plan_ && !defer_quadrature_) KokkosFFT::execute(*forward_coarse_real_plan_, projected_coarse_real_, projected_coarse_spectrum_, KokkosFFT::Normalization::none); } @@ -7136,6 +7165,8 @@ bool GeneralWorkspace::modal_core(const std::size_t stage) { envelope_transform(); if (held_plasma() && !held_radial()) held_transform(); + if (defer_quadrature_ && project_coarse_) + coarse_pending_ = true; return true; } @@ -9205,12 +9236,13 @@ void GeneralWorkspace::held_response(const std::size_t stage, } } -void GeneralWorkspace::held_transform() { +void GeneralWorkspace::held_transform(const bool coarse_only) { for (std::size_t index = 0; index < 2; ++index) { auto& stream = held_streams_[index]; - KokkosFFT::execute(*stream.forward, stream.projected, stream.spectrum, - KokkosFFT::Normalization::none); - if (stream.coarse_forward) + if (!coarse_only) + KokkosFFT::execute(*stream.forward, stream.projected, stream.spectrum, + KokkosFFT::Normalization::none); + if (stream.coarse_forward && (coarse_only || !defer_quadrature_)) KokkosFFT::execute(*stream.coarse_forward, stream.coarse, stream.coarse_spectrum, KokkosFFT::Normalization::none); } @@ -9223,6 +9255,7 @@ void GeneralWorkspace::held_transform() { const auto inverse = held_inverse_; const auto active = active_mask_; const auto crop = crop_index_; + const bool include_coarse = !coarse.empty() && (coarse_only || !defer_quadrature_); cartesian_for( "solprop_native.held.spectral_assembly", execution_space_, projected.extent(0), problem_.n_omega, KOKKOS_LAMBDA(const std::size_t row, const std::size_t omega) { @@ -9230,14 +9263,16 @@ void GeneralWorkspace::held_transform() { if (active(omega)) { const double value = inverse(omega); const Complex factor(0.0, -value); - projected(row, bin) += - factor * product(row, bin) + value * value * edge(row, bin); - if (!coarse.empty()) + if (!coarse_only) + projected(row, bin) += + factor * product(row, bin) + value * value * edge(row, bin); + if (include_coarse) coarse(row, bin) += factor * coarse_product(row, bin) + value * value * coarse_edge(row, bin); } else { - projected(row, bin) = Complex{}; - if (!coarse.empty()) + if (!coarse_only) + projected(row, bin) = Complex{}; + if (include_coarse) coarse(row, bin) = Complex{}; } }); @@ -9368,7 +9403,7 @@ void GeneralWorkspace::project(const std::size_t first_node, const Interleaved fine(reinterpret_cast(projected_target.data()), problem_.n_modes, 2 * problem_.n_to); engine_gemm(transform_space_for(node_count), scaled_projection_, nodes, fine, 0.0); - if (!problem_.projection_coarse_matrix.empty()) { + if (project_coarse_ && !problem_.projection_coarse_matrix.empty()) { const Interleaved coarse(reinterpret_cast(coarse_target.data()), problem_.n_modes, 2 * problem_.n_to); engine_gemm(transform_space_for(node_count), scaled_projection_coarse_, nodes, @@ -9379,7 +9414,7 @@ void GeneralWorkspace::project(const std::size_t first_node, } if (whole_transform && use_blas_for_projection()) { Kokkos::Profiling::pushRegion("solprop_native.transform.projection.kokkosblas"); - if (problem_.family == FieldFamily::envelope && + if (problem_.family == FieldFamily::envelope && project_coarse_ && !scaled_projection_stacked_complex_.empty()) { selected_gemm(transform_space_for(node_count), scaled_projection_stacked_complex_, node_source, stacked_target, Complex{}, true); @@ -9397,7 +9432,7 @@ void GeneralWorkspace::project(const std::size_t first_node, } else if (problem_.family == FieldFamily::envelope) { selected_gemm(transform_space_for(node_count), scaled_projection_complex_, node_source, projected_target, Complex{}, true); - } else if (!scaled_projection_stacked_.empty()) { + } else if (project_coarse_ && !scaled_projection_stacked_.empty()) { selected_gemm(transform_space_for(node_count), scaled_projection_stacked_, real_source, real_stacked_target, 0.0, true); const auto& stacked = real_stacked_target; @@ -9425,7 +9460,7 @@ void GeneralWorkspace::project(const std::size_t first_node, const auto first_row = first_node * problem_.n_pol_transform; const bool tail = node_count != problem_.node_chunk; const bool last_chunk = first_node + node_count == problem_.n_nodes; - const bool stacked = !problem_.projection_coarse_matrix.empty(); + const bool stacked = project_coarse_ && !problem_.projection_coarse_matrix.empty(); if (problem_.family == FieldFamily::envelope) { const auto polarisation = Kokkos::subview( node_source, std::make_pair(std::size_t{0}, rows), Kokkos::ALL()); @@ -9457,7 +9492,16 @@ void GeneralWorkspace::project(const std::size_t first_node, }); } } else { - const auto packed = tail ? tail_projection_complex_ : chunk_projection_complex_; + // The first half of the existing stacked scratch is a contiguous fine + // matrix. Alias it when this trial stage does not need the coarse rows. + const auto scratch = + tail ? tail_projection_stacked_complex_ : chunk_projection_stacked_complex_; + const ComplexView2 packed = + scratch.empty() + ? (tail ? tail_projection_complex_ : chunk_projection_complex_) + : ComplexView2(Kokkos::subview( + scratch, std::pair{std::size_t{0}, problem_.n_modes}, + Kokkos::ALL())); const auto source = scaled_projection_complex_; Kokkos::parallel_for( "solprop_native.transform.pack_projection_complex", @@ -9498,7 +9542,14 @@ void GeneralWorkspace::project(const std::size_t first_node, }); } } else { - const auto packed = tail ? tail_projection_ : chunk_projection_; + const auto scratch = + tail ? tail_projection_stacked_ : chunk_projection_stacked_; + const RealView2 packed = + scratch.empty() + ? (tail ? tail_projection_ : chunk_projection_) + : RealView2(Kokkos::subview(scratch, + std::pair{std::size_t{0}, problem_.n_modes}, + Kokkos::ALL())); const auto source = scaled_projection_; Kokkos::parallel_for( "solprop_native.transform.pack_projection", @@ -9526,7 +9577,7 @@ void GeneralWorkspace::project(const std::size_t first_node, const auto n_pol_transform = problem_.n_pol_transform; const auto n_to = problem_.n_to; const bool envelope = problem_.family == FieldFamily::envelope; - const bool coarse = !projection_coarse.empty(); + const bool coarse = project_coarse_ && !projection_coarse.empty(); Kokkos::parallel_for( "solprop_native.transform.projection.custom", Range2(transform_space_for(node_count), {0, 0}, {problem_.n_columns(), n_to}), @@ -11056,12 +11107,27 @@ void GeneralWorkspace::finish(const double dz) { execution_space_.fence("solprop_native accepted-step boundary"); } +void GeneralWorkspace::finish_quadrature() { + if (!coarse_pending_) + return; + if (problem_.family == FieldFamily::envelope) + KokkosFFT::execute(*forward_coarse_complex_plan_, projected_coarse_complex_, + projected_coarse_spectrum_, KokkosFFT::Normalization::none); + else + KokkosFFT::execute(*forward_coarse_real_plan_, projected_coarse_real_, + projected_coarse_spectrum_, KokkosFFT::Normalization::none); + if (held_plasma()) + held_transform(true); + coarse_pending_ = false; +} + void GeneralWorkspace::score_quadrature() { quadrature_error_ = Kokkos::Experimental::quiet_NaN_v; quadrature_stale_ = false; if (quadrature_partials_.empty() || (!problem_.nonlinear && !problem_.plasma && !dispersive_)) return; + finish_quadrature(); // SEP 6's dispersive products join the embedded rule's spectrum once, when scored. if (dispersive_) dispersive_->finish_coarse(); @@ -11208,6 +11274,7 @@ void GeneralWorkspace::promote_fsal() { void GeneralWorkspace::clear_fsal() { has_recycled_ = false; pending_fsal_ = false; + coarse_pending_ = false; if (!stash_current_.empty()) Kokkos::deep_copy(execution_space_, stash_current_, 0.0); stash_is_current_ = false; @@ -11388,6 +11455,8 @@ std::vector GeneralWorkspace::projected_to_host(const bool coarse) { throw std::logic_error("this quadrature has no embedded rule"); if (coarse && dispersive_) dispersive_->finish_coarse(); + if (coarse) + finish_quadrature(); execution_space_.fence("solprop_native projected_to_host"); const auto host = Kokkos::create_mirror_view_and_copy( Kokkos::HostSpace{}, coarse ? projected_coarse_spectrum_ : projected_spectrum_); @@ -11938,6 +12007,7 @@ void GeneralWorkspace::close() noexcept { table_ready_ = false; has_recycled_ = false; pending_fsal_ = false; + coarse_pending_ = false; stash_is_current_ = false; raman_kernel_ready_ = false; } diff --git a/native/tests/test_band_plasma.py b/native/tests/test_band_plasma.py index d73809bee..e824835e2 100644 --- a/native/tests/test_band_plasma.py +++ b/native/tests/test_band_plasma.py @@ -23,6 +23,7 @@ import dataclasses import json +import math import os import subprocess import sys @@ -467,6 +468,82 @@ def test_every_stage_of_a_trial_matches_the_reference( assert_regression_guard(got, want, f"{kind} {integration} stage {stage}") +@pytest.mark.parametrize("kind", ["multimode", "analytic", "mixture"]) +@pytest.mark.parametrize("chunk", [1, 4]) +def test_quadrature_keeps_accepted_scores_through_rejection_and_reset( + kind: str, chunk: int +) -> None: + """Score held field currents once, preserving accepted observations on rejection. + + The norm difference divided by the fine norm takes at most eight source + rounding budgets when both projection norms are comparable. Accumulate the + existing seven-stage, two-step plasma budget before either backend runs. + """ + grid, original = modal(kind) + problem = dataclasses.replace( + original, + transform=dataclasses.replace(original.transform, node_chunk=chunk), + stepper=dataclasses.replace(original.stepper, integrator=None), + ) + bound = march_bound(problem, grid, STEPS) + score_bound = 8 * bound + assert bound > 0.0 + specification = StatsSpec(names=("quadrature_error", "electron_density")) + reference = NumpyRefBackend(problem, grid=grid) + backend = NativeBackend(problem, grid=grid) + state = backend._native_state # pyright: ignore[reportPrivateUsage] - score lifetime + dz = float(problem.stepper.save_planes[-1]) / STEPS + try: + for engine in (reference, backend): + engine.set_initial(problem.initial.field0, None) + engine.propagate_to(dz) + accepted = backend.observe().require_field().copy() + score = state.quadrature_error + assert math.isfinite(score) + assert state.quadrature_peak() == score # One accepted step since the reset. + assert math.isnan(state.quadrature_peak()) + expected = reference.stats(specification)["quadrature_error"] + np.testing.assert_allclose(score, expected, rtol=0.0, atol=score_bound) + + # Leave an unaccepted trial in scratch, then observe the accepted field and + # plasma statistics before the next trial replaces that scratch again. + backend._trial_step(dz, 0.5 * dz) # pyright: ignore[reportPrivateUsage] - rejected trial + observed = backend.stats(specification) + assert state.quadrature_error == score # Stored accepted value, no arithmetic. + assert math.isnan(state.quadrature_peak()) + np.testing.assert_array_equal(backend.observe().require_field(), accepted) + np.testing.assert_allclose( + observed["quadrature_error"], expected, rtol=0.0, atol=score_bound + ) + for engine in (reference, backend): + engine.propagate_to(2 * dz) + want = reference.observe().require_field() + got = backend.observe().require_field() + assert np.isfinite(got).all() + np.testing.assert_allclose( + got, want, rtol=0.0, atol=bound * float(np.max(np.abs(want))) + ) + np.testing.assert_allclose( + state.quadrature_error, + reference.stats(specification)["quadrature_error"], + rtol=0.0, + atol=score_bound, + ) + backend.set_initial(problem.initial.field0, None) + assert math.isnan(state.quadrature_error) + assert math.isnan(state.quadrature_peak()) + backend.propagate_to(dz) + np.testing.assert_allclose( + backend.observe().require_field(), + accepted, + rtol=0.0, + atol=bound * float(np.max(np.abs(accepted))), + ) + finally: + backend.close() + reference.close() + + @pytest.mark.gpu def test_a_trapezoid_problem_reports_the_run_diagnostics() -> None: """The trapezoid answers the diagnostics from the band scan with its one tap.""" From b3bdd0f2959910d0bf7f88aed23ea32a5aea1eb1 Mon Sep 17 00:00:00 2001 From: John Travers Date: Sat, 24 Oct 2026 13:12:22 +0100 Subject: [PATCH 3/4] Exercise band-filter boundaries and high dynamic range --- native/tests/native/test_core.cpp | 48 +++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/native/tests/native/test_core.cpp b/native/tests/native/test_core.cpp index 4e5c73e7d..c2c6e3505 100644 --- a/native/tests/native/test_core.cpp +++ b/native/tests/native/test_core.cpp @@ -2217,6 +2217,54 @@ TEST(GeneralState, BandWindowSkipsOnlyExactZeros) { } } +TEST(PlasmaBandScan, FilteringPreservesShortRowsBlockTailsAndDynamicRange) { + const std::vector taps{1.1, -0.1, 0.0625, -0.025, 0.0125}; + const int half = static_cast(taps.size()) - 1; + using TeamPolicy = Kokkos::TeamPolicy; + const TeamPolicy probe(native::ExecutionSpace{}, 1, 1); + const int widest = probe.team_size_max(EmptyTeamBody{}, Kokkos::ParallelForTag()); + for (const int samples : {5, 17, 255, 257, 519}) { + SCOPED_TRACE(samples); + std::vector row(static_cast(samples)); + for (int sample = 0; sample < samples; ++sample) + row[static_cast(sample)] = sample % 3 == 0 ? 1e100 + : sample % 3 == 1 ? -1e100 + : 1e-100; + for (const int team_size : {1, 2, 4}) { + if (team_size > widest) + continue; + const auto scan = windowed_band_scan(row, taps, 0.125, team_size); + const auto reflected = [&](const int index) -> long double { + if (index < 0) + return 2.0L * row.front() - row[static_cast(-index)]; + if (index >= samples) + return 2.0L * row.back() - + row[static_cast(2 * (samples - 1) - index)]; + return row[static_cast(index)]; + }; + for (int sample = 0; sample < samples; ++sample) { + long double expected = taps[0] * reflected(sample); + long double magnitude = std::abs(expected); + for (int tap = 1; tap <= half; ++tap) { + const auto before = reflected(sample - tap); + const auto after = reflected(sample + tap); + expected += taps[static_cast(tap)] * (before + after); + magnitude += std::abs(taps[static_cast(tap)]) * + (std::abs(before) + std::abs(after)); + } + // Centre + ordered tap pairs + reflections, in both computations. This + // also covers platforms whose long double has no extra mantissa bits. + const double bound = (8.0 * half + 16.0) * kEpsilon * + std::max(1.0, static_cast(magnitude)); + ASSERT_GT(bound, 0.0); + const auto actual = scan.filtered[static_cast(sample)]; + ASSERT_TRUE(std::isfinite(actual)); + EXPECT_NEAR(actual, static_cast(expected), bound); + } + } + } +} + TEST(State, RepeatedTrialsOverwriteFullCropAndStageZeroAccumulators) { auto problem = base_problem(native::FieldFamily::envelope); problem.nonlinear = true; From 75d4fe18ce0f9f5f34b677fa2dea33480d368a41 Mon Sep 17 00:00:00 2001 From: John Travers Date: Sat, 24 Oct 2026 14:54:38 +0100 Subject: [PATCH 4/4] Document native CPU measurements and retained policies --- CHANGELOG.md | 6 + docs/dev/benchmarks.md | 228 ++++++++++++++++++++++ docs/dev/cpu-performance.md | 33 +++- docs/guide/performance/building-native.md | 21 +- 4 files changed, 284 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 90e57fa85..32492ea7a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -73,6 +73,12 @@ The first public release. It contains: - Each backend is tested against the reference implementation within derived error bounds. - **The native backend.** + - Native CPU optimisation: a coherent NEON FFTW build for Apple silicon, long-row + Raman planning with a compact FFTW/vDSP selection rule, and deferred embedded + transverse diagnostics for ordinary DOPRI5 modal trials. Accelerate BLAS, kernel + schema 35, frontend signatures, protected arithmetic and state layouts are + unchanged. See the [CPU measurements](https://solprop.readthedocs.io/en/latest/dev/benchmarks.html#native-cpu-fft-and-diagnostic-scheduling) + for the lightweight validation scope and rejected traversal and sizing leads. - `NativeBackend` runs one trial step per native call, with states, views and FFT plans resident in the execution space; `NativeBatchBackend` runs fixed-schedule batches with fused accepted updates. diff --git a/docs/dev/benchmarks.md b/docs/dev/benchmarks.md index 5fa1b5f18..82c757db8 100644 --- a/docs/dev/benchmarks.md +++ b/docs/dev/benchmarks.md @@ -430,3 +430,231 @@ ledgers, the input field adds one `F`; Kokkos keeps another `F` of host state, a export passes through two more `F` at the Python boundary. Host Adams trials form about `4F` of temporary arrays. Capacities and transients do not all become resident at once, so their sum overstates the measured peak. + +## Native CPU FFT and diagnostic scheduling + +This campaign uses the frozen `benchmarks ab` harness and unchanged registry recipes, +under an approved exception to standard-run and whole-suite performance reporting. +These are **lightweight measurements**, not a standard-run adoption verdict. Every +reported propagation window contains ten prescribed steps: `fixed-rk` for modal and +radial cases, `fixed-rk4ip` for Cartesian Kerr. Construction, first-window cost and +observations are recorded separately; no window is extrapolated to a full run. + +The baseline is `9521a88`; the retained implementation is `b3bdd0f`. Each arm has a +separate source/environment/native wheel. On the Apple M5 Max (18 cores, 128 GiB, +macOS 27.0.1), Homebrew Clang 23.1.3 uses Accelerate BLAS and one LLVM OpenMP runtime. +The source-built FFTW 3.3.11 control and NEON library have identical M1 flags, +deployment target 14.0 and cycle-counter support, with both precisions and coherent +threading libraries. The two-socket Xeon Gold 6240 uses exclusive 72-logical-CPU +allocations and one hardware thread per selected physical core. Its independent +bundles are GCC 16.2/OpenBLAS 0.3.34/FFTW 3.3.11 with GNU OpenMP, and IntelLLVM +2026.1.1/MKL 2026.1 (including MKL FFT) with Intel OpenMP. Timed arms run sequentially, +without concurrent builds. Python is 3.14.5 on the Mac and 3.14.7 on the Xeon; +NumPy is 2.5.2, Kokkos 5.2.1, Kokkos Kernels 5.2.0 and Kokkos-FFT 2.0.0. + +### NEON codelets and Raman plans + +A same-source control with matched compiler flags changes only FFTW's explicit NEON +codelets. Its eight-case screen at budgets 1 and 18 has geometric mean ratio 0.933; +this one-round/two-repetition screen isolates the library effect but is not a +confirmation. Full and tail microchecks include H₂'s 65536×21/12 and CtC's 131072×3/1 +rows at serial and parallel widths, raw FFT normalization and exact padding zeros. +ESTIMATE and MEASURE both retain short parallel shapes where vDSP is faster, so +complete vDSP removal is rejected. + +With NEON in both arms, replacing the previous Raman choice with tuned Kokkos-FFT +has the following three-round/five-repetition warmed-window ratios. The rest of the +source is held fixed for this provider/planning comparison. + +| Case | 1 thread | 18 threads | +|---|---:|---:| +| CtC | 0.698 [0.694, 0.699] | 0.966 [0.961, 0.971] | +| H₂ Raman | 0.480 [0.459, 0.485] | 0.788 [0.767, 0.795] | +| GNLSE Raman shock | 0.828 [0.817, 0.830] | 0.839 [0.839, 0.846] | + +Measured planning has a real setup cost. In this ordered Raman comparison, median +serial CtC construction changes from 3.00 to 22.83 s, while a warmed ten-step window +changes from 6.00 to 4.18 s. H₂ construction changes from 1.09 to 1.81 s and its +window from 1.28 to 0.62 s. Wisdom is process-local and previous constructions affect +these costs; they are not universal constants. Explicit `SOLPROP_FFTW_PLAN=estimate` +trades this measured setup for cheaper construction. Parallel default planning stays +ESTIMATE. No full-propagation payback is inferred from these ten-step windows. + +### Retained Mac combination + +The fresh current-library baseline and retained NEON/provider/scheduling candidate +use three alternating rounds and five repetitions at each budget. Ratios are the +median of paired round medians, with the range across rounds in brackets. + +| Case | 1 thread | 18 threads | +|---|---:|---:| +| RDW | 0.888 [0.877, 0.901] | 0.831 [0.810, 0.835] | +| CtC | 0.595 [0.577, 0.604] | 0.859 [0.858, 0.878] | +| Radial 256 | 0.893 [0.853, 0.912] | 0.962 [0.934, 0.967] | +| Resolved argon | 0.954 [0.950, 0.961] | 0.886 [0.879, 0.898] | +| H₂ Raman | 0.472 [0.457, 0.474] | 0.718 [0.712, 0.718] | +| Averaged helium RDW | 0.891 [0.880, 0.901] | 0.902 [0.891, 0.916] | +| Averaged argon | 0.918 [0.909, 0.922] | 0.915 [0.905, 0.920] | +| GNLSE soliton | 0.889 [0.882, 0.893] | 0.870 [0.867, 0.871] | + +The geometric mean is 0.792 at one thread and 0.865 at the production budget of 18 +(0.828 across both budgets). Using each fresh process's warmed median as one sample, +the conservative aggregate `(candidate median + 2 MAD) / (baseline median − 2 MAD)` +is 0.805 and 0.878 respectively. All denominators are positive. These spread checks +support the retained change within the lightweight scope. All fixed windows retain +ten accepted steps, no rejected steps and 61 source evaluations; field differences +are at most 2.1e-16 in relative L2 norm. A one-round/two-repetition guard screen at budget 18 gives radial +512 1.012, Cartesian plasma 128² 0.894, line 0.976, general modal 0.892, scalar wall +0.936 and vector wall 0.990. The 1.2% radial variation does not exceed the A/A spread. + +The retained library also passes the Cartesian Kerr confirmation (three rounds, +two repetitions, ten `fixed-rk4ip` steps). No Cartesian traversal change is included: + +| Case | 1 thread | 18 threads | +|---|---:|---:| +| Kerr 192² | 0.843 [0.832, 0.849] | 0.869 [0.718, 0.924] | +| Kerr 384² | 0.837 [0.834, 0.860] | 0.880 [0.880, 0.895] | + +The geometric mean is 0.840 at one thread and 0.874 at 18 (0.857 combined); +the corresponding conservative aggregates are 0.860 and 0.919. Fields differ by at most 1.4e-18 and work +matches. The 192² parallel spread is wider, so its precise gain is less certain; +all three paired rounds still favour the candidate. These results measure the +library change with the original traversal. + +### Diagnostic scheduling + +The affected set is ordinary DOPRI5 modal CtC, RDW, resolved argon and H₂ Raman. +The implementation defers coarse projections to the final source and their FFTs and +held-current assembly to accepted scoring. General-modal, dispersive, low-storage, +Adams and envelope-plasma scheduling remains eager. A serial full-source profile +reduced RDW's forward real FFT calls from 427 to 274 per ten-step window; CtC's +fell from 1464 to 1311. Persistent field-buffer ledgers are unchanged. FFT plan heap +estimates can vary with planning and allocator state. + +These confirmations use three alternating rounds and five repetitions. Each ratio +is the median of the per-round paired medians; brackets give the range over rounds. +A ratio below one favours the candidate. + +| Case | M5 Max, 18 threads | GCC Xeon, 1 thread | GCC Xeon, 18 threads | Intel Xeon, 1 thread | Intel Xeon, 18 threads | +|---|---:|---:|---:|---:|---:| +| RDW | 0.803 [0.755, 0.845] | 0.930 [0.925, 0.946] | 0.823 [0.820, 0.843] | 0.900 [0.887, 0.916] | 0.822 [0.813, 0.828] | +| CtC | 0.906 [0.888, 0.917] | 0.897 [0.894, 0.901] | 0.946 [0.945, 0.949] | 0.904 [0.893, 0.910] | 0.937 [0.933, 0.937] | +| Resolved argon | — | 0.998 [0.984, 1.002] | 0.949 [0.941, 0.953] | 0.966 [0.963, 0.976] | 0.966 [0.956, 0.972] | +| H₂ Raman | 0.933 [0.914, 0.967] | 1.007 [0.986, 1.082] | 0.982 [0.981, 0.988] | 0.985 [0.976, 1.044] | 0.973 [0.957, 0.977] | + +At one thread the GCC affected-set geometric mean is 0.957 and Intel's is 0.938; +at 18 threads they are 0.923 and 0.922. The corresponding conservative aggregates +from fresh-process medians are 0.972/0.931 for GCC and 0.951/0.929 for Intel, with +positive denominators. Both production-budget gains exceed 5% and satisfy the +conservative spread criterion. Mac's isolated 18-thread ratio is 0.879, with a +conservative aggregate of 0.950. Across both budgets, GCC is 0.940 and Intel 0.930. +The matched fixed windows retain ten accepted steps, no rejections and 61 source evaluations; their +field differences are at rounding level (at most 1.9e-16 in these confirmations). +The Mac's first 18-thread screen and one A/A round showed large host timing +variation and remain inconclusive. The runner's busy-at-launch flag does not establish its cause; load averages also include prior +owned work. A subsequent three-round A/A control gave RDW 1.010 [1.008, 1.054] +and H₂ 1.015 [0.995, 1.040]. +Those controls do not justify interpreting a few-percent change as a gain. + +### Retained Xeon combination + +Each compiler bundle also compares the complete retained candidate with its frozen +baseline. Raman cases use three alternating rounds/five repetitions at budgets 1 +and 18; other main cases and geometry guards use one round/two repetitions. The +Raman confirmations are: + +| Case | GCC, 1 thread | GCC, 18 threads | Intel, 1 thread | Intel, 18 threads | +|---|---:|---:|---:|---:| +| CtC | 0.790 [0.789, 0.792] | 0.963 [0.953, 0.966] | 0.912 [0.908, 0.914] | 0.939 [0.934, 0.942] | +| H₂ Raman | 0.824 [0.817, 0.841] | 1.003 [0.988, 1.030] | 0.977 [0.966, 0.986] | 0.976 [0.974, 0.977] | + +The geometric means at budgets 1/18 are 0.807/0.983 for GCC and 0.944/0.957 for +Intel (combined 0.890 and 0.951). GCC's parallel Raman-only conservative aggregate +is 1.004, so that subset does not establish a gain beyond spread; the confirmed +shared scheduling gain is assessed over its complete affected set above. The +serial GCC improvement includes measured FFTW planning; Intel uses MKL FFT, whose +FFTW compatibility layer ignores planning flags. No benefit is inferred for +unmeasured Linux Arm codelets. + +Main-case guard screens, followed by A/A controls for the small averaged and GNLSE +windows, are reported separately from the confirmed affected-set gains: + +| Case | GCC, 1 thread | GCC, 18 threads | Intel, 1 thread | Intel, 18 threads | +|---|---:|---:|---:|---:| +| RDW | 0.928 | 0.834 | 0.929 | 0.827 | +| Radial 256 | 1.002 | 0.962 | 1.012 | 0.991 | +| Resolved argon | 1.000 | 0.947 | 0.988 | 0.965 | +| Averaged helium RDW | 1.014 | 1.051 | 1.011 | 1.001 | +| Averaged argon | 1.004 | 1.025 | 0.994 | 1.017 | +| GNLSE soliton | 1.008 | 1.005 | 1.019 | 1.039 | + +The initial GCC averaged-helium 18-thread ratio of 1.051 and Intel GNLSE ratio of +1.039 triggered separate three-round/five-repetition A/A and A/B controls: + +| Control case | A/A, 1 thread | A/B, 1 thread | A/A, 18 threads | A/B, 18 threads | +|---|---:|---:|---:|---:| +| GCC averaged helium | 0.997 [0.989, 0.998] | 1.014 [1.005, 1.016] | 1.016 [0.998, 1.128] | 1.012 [1.012, 1.014] | +| GCC averaged argon | 0.987 [0.984, 0.994] | 1.014 [1.009, 1.031] | 1.001 [1.000, 1.001] | 1.010 [0.983, 1.017] | +| GCC GNLSE | 0.991 [0.818, 0.997] | 1.005 [1.003, 1.016] | 0.999 [0.995, 0.999] | 1.010 [0.992, 1.017] | +| Intel GNLSE | 0.981 [0.746, 0.993] | 1.005 [0.996, 1.008] | 1.018 [0.985, 1.061] | 0.999 [0.991, 1.018] | + +The initial 4–5% slowdowns do not repeat. GCC's small-window confirmation aggregate +is 1.011, a roughly 1% cost that remains visible; this does not justify additional +specialisation. Large A/A outliers make those controls inconclusive at fine +resolution. No improvement is claimed for these unchanged small-window routes. + +At budget 18, GCC's six geometry guards span 0.977–1.018 (aggregate 1.000), and +Intel's span 0.998–1.017 (aggregate 1.007). These cover radial 512, Cartesian plasma +128², line, general modal and scalar/vector walls. Cartesian Kerr 192²/384² guards +are 0.997/0.999 for GCC and 1.001/1.004 for Intel, with ten steps and matched work. +They support the absence of a large regression; they do not establish percent-level +improvements on routes whose implementation is unchanged. + +### Verification scope + +Focused checks cover raw FFT normalization, full/tail shapes, threaded execution, +provider selection, planning overrides, coherent library coexistence, repeated +construction/closure, source/step parity, held-plasma diagnostics after discarded +trials, FSAL, reset, accepted values and peaks, and band-filter boundaries and high +dynamic range. Native suites passed on both Xeon bundles: GCC has 137 C++ tests and +378 Python tests passing; Intel has 133 and 376. Provider/platform/device skips +remain explicit. Mac final focused Python checks have 118 passes plus 13 lifetime +and import checks. Independent bounded smooth-law adaptive fixtures pass on both +arms and all three builds against accuracy envelopes fixed from tighter NumPy +references; these small fixtures do not certify full production trajectories. + +The macOS distribution helper was exercised with freshly built dependencies and +wheel repair. A clean installation discovers both bundled FFTW precisions as NEON, +one bundled OpenMP runtime and the required third-party licences, with the Apple M1 +baseline and macOS 14 target. No workflow edit is needed: the existing macOS wheel +job calls this dependency recipe. Full gates, slow and coverage lanes are deferred +to CI at the maintainer's request; GPU/device checks were unavailable on these CPU +hosts. Focused formatting, lint, spelling and diff checks passed. Supplemental +CPU/wall, observation, load and memory receipts are collected outside A/B timing; +no additional speed claim is taken from those resource measurements. + +### Rejected leads + +| Lead | Bounded evidence | Decision | +|---|---|---| +| Remove vDSP completely | NEON FFTW still loses on short parallel rows, including 8192×1 and 16384×1, under both ESTIMATE and MEASURE | Keep vDSP for short parallel power-of-two rows; use FFTW for serial and long rows with NEON | +| Tap-first band filter, preserving each sample's tap order | Mac serial RDW 0.973 [0.970, 0.982], argon 1.079 [1.061, 1.090]; GCC/Intel RDW/argon aggregates 0.911/0.900 across 1/18 threads | Reject the shared traversal and avoid an architecture-specific fork; retain the new boundary/dynamic-range tests | +| Cartesian 8×32 tile plus unused-load elimination | 192² Kerr: serial 0.980 [0.939, 0.997], 18 threads 0.980 [0.959, 1.005], three rounds/two repetitions | Reject the approximately 2% gain; keep 32×16 and the existing loads | +| Force one Accelerate BLAS thread | Five-case confirmation at budget 18: geometric mean 1.000, three rounds/five repetitions | Keep the existing provider-managed startup policy | + +The post-FFTW 192² serial profile assigns 4.62 s to the separable fold, 4.72 s to +complex FFTs and 2.11 s to temporal pack/crop within a 12.16 s instrumented window. +A separate non-IPO compiler diagnostic reports control-flow and call barriers in the +crossing and final reduction loops; the measured production binaries retain IPO. +The tile experiment changes locality without changing stored layouts or per-column +summation order, but its complete-window gain does not meet the greater-than-5% +threshold. Instrumented profiles are diagnostic evidence, not timings used for +retaining a change. + +Independent startup-only pool screens at budget 18 varied kernel, ordinary FFT, +Raman full, Raman tail and BLAS widths separately to 1 or 4. Both Xeon bundles +showed aggregate slowdowns for every reduced-width screen (the smallest was 1.095 +for Intel's four-thread Raman tail). Mac responses were mixed, including opposite +RDW and argon directions. No shared dimension threshold or pool policy is changed; +explicit overrides and global caps retain their existing semantics. A 36-core +sizing confirmation is unnecessary because no sizing change is retained. diff --git a/docs/dev/cpu-performance.md b/docs/dev/cpu-performance.md index 42394cd4b..51970d984 100644 --- a/docs/dev/cpu-performance.md +++ b/docs/dev/cpu-performance.md @@ -32,6 +32,11 @@ time divided by the native backend's: above 1, the native backend is faster. to every logical CPU. - **Accelerate on macOS.** It matched or beat both OpenBLAS variants and native GEMM, except at one thread on the TE/TM wall, where it is the fastest option. +- **NEON FFTW on Apple silicon**, in both precisions, with the Arm cycle counter, + M1 distribution baseline and existing deployment target. FFTW and Kokkos share + one OpenMP runtime. This enables FFTW's explicit codelets; compiler vectorisation + alone does not select them. The [build recipe](../guide/performance/building-native.md) + keeps x86 SIMD and Linux Arm settings unchanged. - **x86-64-v3 distribution baseline.** An x86-64-v2 build's median cost against x86-64-v3 was within about 2 % on a 32-core AVX-512 server (six pairs, seven configurations) and 4 % on an 18-core-per-socket AVX-512 server, on the benchmark core. @@ -105,9 +110,19 @@ time divided by the native backend's: above 1, the native backend is faster. 0.09 s per 40 source evaluations ([measurements](benchmarks.md#multimode-cpu-step-costs)). - **General modal: plans for the active tail rows only.** 2–5 %. -- **Short serial Raman rows (≤ 16384 samples) on measured FFTW instead of Accelerate.** - Faster in both directions at every size tested. Larger and parallel Raman transforms - stay on Accelerate. +- **Raman provider selection with NEON FFTW.** Serial and long rows use Kokkos-FFT; + short parallel power-of-two rows (≤ 16384 samples) keep vDSP. A non-NEON local + FFTW retains the previous short-serial FFTW rule. Every contiguous FFTW fallback, + including full and tail batches, honours the requested planning effort. Serial + measured planning improves warmed execution but can substantially increase + construction time; `SOLPROP_FFTW_PLAN=estimate` remains available for short runs. +- **Deferred ordinary DOPRI5 modal diagnostics.** Only the last trial stage projects + the embedded transverse rule, and accepted-step scoring transforms it. Fine-only + stages avoid stacked coarse GEMM work. Held field-plasma streams follow the same + schedule; other methods and routes retain their existing consumers. Accepted + values and peaks survive rejected trials, and existing scratch is reused without + a new field allocation. The [lightweight measurements](benchmarks.md#native-cpu-fft-and-diagnostic-scheduling) + record the affected windows and their limits. - **Resident Adams history: no dense-output volumes when planes are landed on** ([measurements](benchmarks.md#resident-adams-storage)). Storage went from 29 to 18 field volumes. The projected 384² Kerr-beam peak went from 62.7 to @@ -121,6 +136,18 @@ time divided by the native backend's: above 1, the native backend is faster. did not explain the gap. - **A direct FFTW execution adapter instead of Kokkos-FFT's.** Matched probes showed no wrapper overhead. +- **Complete vDSP removal.** Even with NEON codelets, short parallel FFTW plans + regress former vDSP shapes. Keep the compact construction-time rule above. +- **Tap-first CPU band filtering.** It preserves tap order and passed the boundary, + support-window and high-dynamic-range checks. It improved the x86 windows but + repeatedly slowed the Apple-silicon argon window by 6–9%; an architecture-specific + fork was not justified. The original traversal remains, with stronger tests. +- **An 8×32 Cartesian crossing tile and omitted unused register loads.** Three + alternating rounds on the 192² Kerr window reduced the initial apparent gain to + about 2%, below the complexity threshold. The 32×16 traversal remains. +- **A forced one-thread Accelerate BLAS limit.** An initial RDW gain disappeared + over three alternating rounds on RDW, radial 256, two wall cases and general modal: + the geometric mean ratio was 1.000. Existing process-start settings remain. ## Open diff --git a/docs/guide/performance/building-native.md b/docs/guide/performance/building-native.md index ab99918b5..22671d872 100644 --- a/docs/guide/performance/building-native.md +++ b/docs/guide/performance/building-native.md @@ -100,7 +100,22 @@ supports `--enable-neon` for double and single precision on AArch64. The install libraries identify their codelet family in `fftw_version` and `fftwf_version`, which both end in `-neon` for this recipe. Keep these prefixes separate from a comparison's baseline libraries. A pyFFTW build used in the same process must also share this -runtime; see [OpenMP runtimes](#openmp-runtimes-and-pyfftw) below. +runtime. To build the existing pyFFTW dependency against that same owned prefix: + +```sh +pyfftw_version=$(uv run --no-sync python -c 'from importlib.metadata import version; print(version("pyfftw"))') +CC=/opt/homebrew/opt/llvm/bin/clang \ +CXX=/opt/homebrew/opt/llvm/bin/clang++ \ +PYFFTW_FFTW_PREFIX="$native_deps" \ +LDFLAGS="-L$native_deps/lib -Wl,-rpath,$native_deps/lib" \ +uv pip install --python .venv/bin/python --no-deps --no-binary pyfftw \ + --reinstall --no-cache "pyfftw==$pyfftw_version" +``` + +The uncached source build prevents reuse of a wheel linked to another runtime. +Keep the project's locked pyFFTW version when qualifying a comparison. See +[OpenMP runtimes](#openmp-runtimes-and-pyfftw) below for import ordering and the +limits of mixing independently bundled wheels. For a distributable wheel, `native/ci/build-wheel-macos.sh OUT_DIR DEPS_PREFIX` builds the distribution profile, repairs its library references with delocate, and @@ -180,6 +195,10 @@ object suitable for a run log: The library paths (`fft_library`, `blas_library`, `openmp_library`) are the files actually loaded, the reliable way to confirm which FFT and BLAS a process uses. +`raman_fft_option` identifies the compiled Raman wrapper, which can select FFTW for +individual plans even when its value is `ACCELERATE`. The backend's execution +metadata records the selected full-chunk `raman_fft_provider` and full/tail task +counts; selection also accounts for the actual tail width. `native_arch_enabled` reports the package's option only; architecture flags passed by hand through `CXXFLAGS` appear in `compiler_flags`.