Benchmarks

How solprop’s performance is measured and judged, and what the measurements so far indicate. A performance claim names what it measured, the build, the machine and the result; the rules below say how to take it. The numbers in the indicative sections are indicative: each was measured once, on one machine, for one change. Re-measure before relying on any of them. Rows naming numba_cpu record a backend that SEP 43 retired; they can no longer be re-measured, and stay as what was measured.

Where things are

• The benchmarks are benchmarks/ (SEP 38), run with uv run poe bench (python -m benchmarks). track times complete runs of the tracking set (about a dozen representative workloads, with each interface’s default integrator and tolerances) at one thread and the production width, in about an hour per machine; ab compares two builds while a change is made. Running the benchmarks is the manual.

• The history is benchmarks/history/: one CSV per hardware set, named by its CPU model and core count, with a few rows per release. python -m benchmarks history summarises it. Full run bundles stay in the ignored .benchmarks/ directory.

• CPU performance: what was tried records the compiled CPU backend’s performance experiments: what was adopted, measured and rejected, and what remains open. Read it before repeating an experiment.

• backends.md describes each backend’s execution design and the performance policies it adopts.

Rules

• State thread counts and plan flags; never compare across different ones.

• A hot-path change needs a benchmark entry in its pull request or commit message.

• Estimates are labelled as estimates in the docs until a measurement exists.

• Do not tune for the benchmark: the cases are the target simulations at their real sizes.

• A performance change is judged on the aggregate, not case by case: the equal-weight geometric mean of its per-case timing ratios over the cases it can affect, declared before timing, taken separately for the production thread policy and for one thread, with the complete per-case table beside it and a whole-suite aggregate when the work closes. A change is adopted when the production saving is significant: more than 5 %, and the geometric mean of the conservative per-case ratios (candidate median + 2 MAD) / (before median − 2 MAD) below one, a non-positive denominator counting as unresolved. Small regressions on individual cases, or at one thread, are accepted when the aggregate improves; they stay visible in the table and do not automatically veto the change. Numerical bounds are never relaxed to make a change pass. The A/A spread of the same comparison says how large a difference has to be to count (comparing two builds).

What is timed

• Release tracking times complete runs. The tracking set’s cases are complete guides at production sizes, with the library’s defaults, so the history says how long the usual simulations take now. Construction, translation and backend construction are timed separately; an untimed warm-up builds plans first (what is timed).

• A change in progress may time less. ab --length runs part of each guide from the entrance, which is fine for comparing two builds on the same work. Report the fraction run, and do not extrapolate an entrance window to a complete-run cost: the expensive part of a propagation (compression, ionisation, the Raman shift) often comes late.

• A later state of a long propagation is timed by replaying a saved witness into a short run, not by regenerating the propagation.

• Compare two source trees with alternating pairs (A/B/A/B…), and report the A/A spread beside the ratio; a difference inside the spread is not a result.

Reporting a performance claim

The machine record. Every measurement carries the CPU (or device) model, core count, memory, operating system, Python and dependency versions (including the Kokkos build and its profile), thread counts, FFTW plan flags and wisdom state, the source revision and whether the tree was dirty, and the load before and after. A tracking run’s plan and results record these; a hand-run measurement states them. Results are compared only within one machine.

A quiet machine. Absolute numbers are taken on a dedicated workstation or an exclusively allocated node. A busy laptop is fit only for same-process A/B comparisons, labelled as such: read their ratios, not their times. Report the A/A spread of repeated baseline runs beside every ratio; a difference inside it is not a result.

What goes in a pull request. The measured window and what it includes, the build and machine record, the per-case table with medians and spreads, the declared aggregate and its verdict, and any case that regressed. Concise evidence belongs in the change description; raw runs stay local.

Indicative results

Indicative ratios, not absolute times. Each is a short window (see above) on one machine; the rows’ cases, grids and step lengths differ, so rows share no scale. They say which backend or method to expect to win, and by roughly how much.

Backends and builds

Speed-up over numpy_ref, one thread, ten-step windows on an 18-core AVX-512 server (one socket, one NUMA node) with a local Kokkos build (vendor BLAS, native architecture). The capillary cases have 3 or 4 channels on a 513-bin grid; the general-route cases are a two-channel step-index pair and a graded-index guide at 2304 quadrature nodes.

case	numba_cpu	Kokkos CPU
M1 soliton (GNLSE), reduced	2.6×	2.1×
GNLSE soliton with Raman	2.6×	2.6×
M2 (RDW capillary), reduced	4.9×	5.9×
M2 (RDW capillary), full	7.1×	8.6×
capillary, no wall	2.5×	3.9×
capillary, scalar coherent wall	2.1×	3.5×
capillary, TE/TM coherent wall	4.2×	8.8×
capillary taper with wall (coupled flow)	2.4×	2.8×
general route, step-index pair	2.7×	5.5×
general route, graded index	2.9×	4.7×

Threads. numpy_ref does not thread. At four threads Kokkos CPU ran the general route 2 to 6 times faster than numba_cpu and the capillary cases 1.5 to 3 times faster. On larger general-route and coherent-wall workloads (2 to 32 channels, up to 16 384 fine samples) at 18 threads it was 1.2 to 4 times faster than numba_cpu at 18 threads and 6 to 40 times faster than single-threaded numpy_ref. Small capillary grids stop scaling beyond a few threads on every backend.

Kokkos build profiles, same server, seven ten-step capillary and general-route cases. The vendor BLAS and the instruction set are separate gains. Against the reference build (no vendor BLAS, generic x86-64), linking OpenBLAS alone made the straight capillaries 2.6 to 3.8 times faster at one thread (the tapers 1.7) and left the general route unchanged; compiling for x86-64-v3 alone made the general route 1.4 to 1.5 times faster and left the capillaries unchanged. A distribution-style build (x86-64-v3 and the scipy-openblas32 OpenBLAS, ADR-034) ran within about 5 % of the local build (-march=native, the system OpenBLAS) on every case at one and four threads and beat numba_cpu on all seven; x86-64-v4 added 7 to 14 % on the general route only. The distribution wheels now target a lower baseline (CPU performance). At one thread, in milliseconds per ten steps:

case	reference	distribution-style	local	numba_cpu
scalar wall	61.4	21.8	22.8	35.0
TE/TM wall	1547	390	402	480
taper wall	90.3	54.6	56.6	64.7
general pair	188	118	120	241
general GRIN	358	245	244	454

On the same server MKL was 1.4 to 1.8 times faster than OpenBLAS on the TE/TM wall, where the MKL build keeps the custom complex kernel, and level elsewhere.

GPU against a whole CPU socket. On the same seven cases, the Kokkos CUDA build on a data-centre NVIDIA GPU and the HIP build on an AMD Instinct GPU ran 2.7 to 18 times faster than the Kokkos OpenMP build on all 18 cores of the server, the graded-index general route most. These grids are small for a device; the ratios say nothing about large grids.

The tracking set on the release wheel, without Numba

SEP 43’s final tree against the history’s 0.1.0 rows: python -m benchmarks track with CI’s cp314-abi3 wheel at its x86-64-v3 build, three repeats, passive threading, an exclusively allocated AMD EPYC 7763 node at 1 and 16 threads. The rows aren’t appended, because the history holds releases.

case	threads	now, s	0.1.0 native, s	ratio	0.1.0 numba_cpu over now
gnlse_raman_shock	1	43.57	44.16	0.987	1.24
gnlse_soliton	1	0.583	0.578	1.009	1.60
rdw	1	11.71	11.63	1.007	1.80
he_rdw_averaged	1	1.075	1.058	1.016	2.63
argon_blueshift	1	26.50	26.44	1.002	1.86
argon_blueshift_averaged	1	3.941	3.920	1.005	2.79
argon_blueshift_envelope	1	1.250	1.240	1.008	1.98
h2_raman_averaged	1	6.101	6.119	0.997	2.40
kerr_beam_192	1	13.34	12.07	1.105	1.45
gnlse_raman_shock	16	39.05	37.84	1.032	1.37
gnlse_soliton	16	0.587	0.582	1.007	1.62
rdw	16	5.112	5.148	0.993	3.18
he_rdw_averaged	16	0.770	0.785	0.980	3.70
argon_blueshift	16	9.852	9.861	0.999	2.30
argon_blueshift_averaged	16	3.271	3.261	1.003	3.37
argon_blueshift_envelope	16	1.123	1.322	0.850	2.22
h2_raman_averaged	16	5.163	5.250	0.983	2.82
ctc	16	108.4	111.5	0.972	2.44
grin_nonlinear	16	20.50	21.17	0.969	1.79
kerr_beam_384	16	29.10	25.26	1.152	1.24

The geometric mean is 1.002, with the same steps on all twenty, and Numba’s rows are 2.08 times slower. These are runs on different days, so a few per cent is not a result. radial_filament_256 is left out, since its guide changed at revision 2.

• kerr_beam: the rise isn’t the removal. The branch against the main it came from gave 0.9956 and 0.9994 on these rows, with identical fields. Its size matches the transverse edge’s read cost, +10.5 % on the Cartesian case there, though this comparison didn’t bisect it.

• The same work on the Xeon Gold 6240: at one thread, the same steps and ratios of 0.99 to 1.04 in two runs, apart from kerr_beam_192 at 1.18 to 1.20 and rdw at 0.99 to 1.09. At 36 threads the two runs of one wheel differ from each other by 0.84 to 1.57 per configuration, unpinned across two sockets, so the 36-thread column gives no verdict.

Kokkos CPU operation-size threading

The Kokkos CPU backend chooses each operation’s thread width from its dimensions at construction (tuning; the thresholds are under reference tables). It was calibrated on an 18-core arm64 laptop and a two-socket, 36-core AVX-512 server, against one thread on the small cases and against the previous default build on the large windows (best of complete paired rounds):

window	arm64 laptop	two-socket server, all 72 logical CPUs
small cases (eight case/track configurations), geometric mean	2.00× faster	2.06× faster
averaged helium RDW	same within noise	same within noise
resolved blueshift	2.0–2.6× faster	1.8–2.4× faster
reduced comb-to-continuum (revision 1, 1 ps)	2.8–3.1× faster	2.4–3.6× faster
RDW ten-step window	1.81× faster	2.33× faster
reduced comb-to-continuum ten-step window (revision 1, 1 ps)	3.18× faster	3.85× faster
full comb-to-continuum ten-step window	7.9 % faster	2.7 % faster, noisy
general modal	4–11 % faster	same, noisy
Cartesian and radial filaments	same within noise	within 5 %, noisy
scalar / TE-TM walls	44 % / 17 % faster	TE-TM same

Separate thresholds for copies and arithmetic matter: an 8,192-element temporal copy is slower threaded, while arithmetic stages of the same size gain from the full pool. A wide response chunk followed by a full-width FFT needs both at full width, and a serialised coupled-flow matrix function on the scalar wall took about 20 times its threaded kernel’s time; both were corrected before the comparison above. “Noisy” rows did not meet the predeclared stability check; taper, short-line and planar controls were not requalified.

The exponential Adams trial without a coupled flow is two passes around its source evaluation, a prediction and a correction that also reduces the defects’ norms. An element of them costs about ten stage entries, so they take a class of their own, its work the field’s elements times the history’s rows (seven at the default fifth order). Forced to 1-16 threads of a 16-thread process whose idle threads sleep, over single-column capillaries of 401 to 10,241 elements, on 16 cores of a dual AMD EPYC 7643: four threads lost at 4,487 (0.8×) and gained 1.1-2.1× (median of two rounds) from 5,607 to 71,687, while sixteen beat four on argon from 11,207 but lost to them on helium up to 17,927 and lost to one thread at 8,967 (0.7×), a fork and join of sixteen sleeping threads costing 100-200 µs. The first round had the node to itself; in the second, with other jobs on it, the serial times held within 5 % while the threaded ones slowed by up to 2.3×, the sixteen-thread ones most, yet four threads stayed at least as fast as one from 5,607 up. On a loaded Apple M5 Max (indicative) four gained 1.0-1.6× at 4,487-8,967 and sixteen beat four only at 71,687 (4.7× against 2.8×). The class therefore threads at four from 8,192 and takes the pool from 65,536 on both platforms.

A stable kernel team

GNU OpenMP keeps one pool of threads for the main thread’s parallel regions. When a region asks for fewer threads than the last, the surplus threads end; the next wider region starts them again. A four-thread kernel between full-width ones therefore costs a thread creation per surplus thread each time: alternating 16- and 4-thread regions took about 490 us a pair, against 40 us for two 16-thread regions and 23 us for 16 and 1 (one socket of a Xeon Gold 6240, idle threads sleeping; twelve clone3 calls a pair). FFTW’s threaded loops run on the runtime’s whole team and only hand out fewer iterations, so a four-thread FFT plan costs no resize. At 16 threads, h2_raman_averaged (an exponential Adams trial pass and a response at four threads beside a full-width stage) re-created twelve threads every source evaluation, and was slower than at one thread.

A kernel pool that already holds a full-width instance now hands it to medium-width kernels too; a pool without one keeps four threads (averaged helium, whose other kernels are serial, took 14 % longer on a Xeon Gold 6240 with its Adams passes at full width). Workspaces ask for their stage instance first. Interleaved benchmarks ab at 16 threads against the parent, 3 rounds, --length 0.1, passive-bound placement, the same steps and fields to rounding (at most 1e-13 relative):

case	Xeon Gold 6240 (socket 1)	AMD EPYC 7763 (16 cores)
h2_raman_averaged	0.526 [0.521, 0.531]	0.483 [0.483, 0.510]
argon_blueshift	0.717 [0.700, 1.070]	0.806 [0.800, 0.850]
ctc	1.001	0.979
gnlse_raman_shock, gnlse_soliton, rdw, he_rdw_averaged, argon_blueshift_averaged, argon_blueshift_envelope, radial_filament_256, grin_nonlinear, kerr_beam_384	0.984-1.036	0.983-1.106 (node noise)
geometric mean (12 cases)	0.924	0.938

The unchanged cases have no medium-width kernel or none beside a full-width one; their spread is the A/A spread (about 2 % on a Xeon Gold 6240, up to 10 % on the shared EPYC nodes). On a Xeon Gold 6240 h2_raman_averaged now takes 0.79 s at 16 threads against 1.06-1.12 s at one. Every medium width serial instead lost 2.8 times on averaged helium (EPYC geometric mean 1.08); every medium width full, the FFTs included, gained on resolved cases but lost on the averaged ones (EPYC 0.94, Xeon 0.90 with averaged helium 1.22).

The promotion applies everywhere but macOS (stable_kernel_team in core.hpp). On an 18-core Apple M5 Max, whose cores are of two kinds, the 11207-element Adams passes of argon_blueshift at 18 threads instead of four took the window from 3.0 to 3.2 s, and h2_raman_averaged at 18 threads ran 0.86 of main’s time without the promotion against 0.94 with it. The Windows wheel also runs on LLVM’s libomp, which keeps a narrower region’s threads, but there the promotion still gained: on eight cores of a Google Cloud C4D (AMD EPYC 9B45) h2_raman_averaged took 0.906 [0.899, 0.914] of its time with it at eight threads, every other modal case within 1 % (5 rounds, --length 0.2).

Compact CPU teams

A problem of one column on one node (the GNLSE, the mode-averaged capillary) has few operations large enough to thread, and between them serial transforms and products. At 16 threads gnlse_raman_shock (16384 bins, 6289 source evaluations a --length 0.1 window) ran its stage entry on all sixteen threads and its responses on all of them too, its FFTs on four threads within the runtime’s whole team, and gained nothing over one thread (5.86 against 6.02 s on one socket of a Xeon Gold 6240). Its serial kernels each ran 1.2-2.5 times slower than in a one-thread process; with every pool forced serial in the same 16-thread process they ran at their one-thread times, so the idle pool costs nothing and the wide regions do. A per-thread perf stat of the timed window shows both parts: the main thread ran at about 2.7 GHz instead of 3.5 (waking fifteen cores every 0.9 ms evaluation holds the socket in its many-core turbo bin), and its loads found 16 million lines modified in another core’s cache, against 60 thousand.

Such a problem (compact_problem: one column, one node, one transform polarisation, no general modal route, coupled flow or separable geometry, fewer than 65536 fine samples) now runs on a compact team: its kernels at most eight threads wide, medium widths promoted to that team as in a stable kernel team, and its FFT and Raman plans at most eight (see the end of this section for their sizing). FFTW’s threaded loop normally opens the runtime’s whole team (sixteen threads for a four-thread plan); beside an eight-thread kernel team that made GNU OpenMP end and restart threads at every pair, and a whole run took 2.6 times as long. The native backend therefore installs FFTW’s loop callback (FFTW 3.3.9 and later, the OpenMP FFTW only): outside a scope it runs the jobs on the runtime’s team, as FFTW’s own loop does, and inside the scope that every native call of a compact problem opens it runs them on the problem’s team. A build without the callback keeps a compact problem’s transforms serial. An observation now copies the field on the stage’s instance rather than the default one, which would also widen the team.

Sweeps of the whole native pool on the same 16-thread process (SOLPROP_NATIVE_THREADS, 2 rounds, wall s) chose the team: gnlse_raman_shock 5.84 at 16, 4.72 at 4, 4.64 at 6, 4.76 at 8, 5.42 at 12; h2_raman_averaged 0.71 at 16, 0.61 at 4 and at 8; the 33-node cases lose at 4 (argon_blueshift 1.77 against 1.48 s, rdw 0.90 against 0.75 s) and stay with their own sizing. Interleaved benchmarks ab against the parent on a Xeon Gold 6240, 3 rounds, --length 0.1, passive-bound placement, the same steps and fields to the measured-plan floor (at most 6e-14 relative):

case	1 thread	16 threads
gnlse_raman_shock	0.995	0.801 [0.794, 0.831]
h2_raman_averaged	0.992	0.888 [0.883, 0.891]
he_rdw_averaged	1.027 [0.971, 1.055]	0.979
gnlse_soliton, argon_blueshift_averaged	0.983-0.985	0.995-1.017
argon_blueshift_envelope, rdw (not compact)	0.988-1.009	0.999-1.009

gnlse_raman_shock now runs 1.25 times faster at 16 threads than at one (4.76 against 5.94 s), where it ran at the same speed.

On the compact team the transforms kept a 16-thread sizing: the four-thread plans of the 16384-sample transform left half the team idle. A compact problem’s threaded transform now takes its whole team. Three refinements were measured on both machines, interleaved benchmarks ab at 16 threads, 3 rounds, --length 0.1, the same steps (AMD EPYC 7763: 16 cores of a quiet shared node):

change	Xeon Gold 6240	AMD EPYC 7763
threaded transforms at the team width: gnlse_raman_shock	0.943	0.981
the same: h2_raman_averaged	1.003	0.995
Raman transforms from the modal FFT’s start (not kept): gnlse_raman_shock	0.946	1.056
the same: h2_raman_averaged	0.978	1.17 [1.08, 1.21]
one-node products as the one-mode kernel on the team (not kept): gnlse_raman_shock	0.973	1.045

The Raman transforms therefore keep their own start (65536 samples), and the GNLSE’s 1 × 1 synthesis and projection stay with the vendor GEMM. Threading every compact kernel of at least 4096 elements at the team width gained nothing either (h2_raman_averaged 3 % slower on a Xeon Gold 6240), nor did threading compact transforms from 4096 or 8192 samples (averaged argon 1.6 times slower at 4096).

The final rule against the parent of the compact-team commit, 3 rounds, --length 0.1, the same steps (fields within 9e-14), 16 threads; one thread moved by 0.98-1.02:

case	Xeon Gold 6240	AMD EPYC 7763 (16 cores of a quiet shared node)
gnlse_raman_shock	0.774 [0.763, 0.800]	0.965 [0.960, 0.970]
h2_raman_averaged	0.895 [0.872, 0.905]	0.965 [0.957, 0.980]
gnlse_soliton, he_rdw_averaged, argon_blueshift_averaged	0.98-1.05	0.97-1.04

With the package’s own plan mode (threaded rows planned at FFTW_ESTIMATE, where the suite measures every plan) the Intel Xeon Gold 6240 rows were 0.780 and 0.905. On the EPYC the 16 allocated cores lay in three CCXs; its 16-thread gnlse_raman_shock window ran 1.04 times faster than one thread before and 1.08 times after (Xeon: 1.02 and 1.29).

The products of a compact problem run on one OpenBLAS thread. On 16 cores of a Google Cloud C4A (Neoverse V2) the compact team made gnlse_raman_shock slower at 16 threads than main (1.13) and than itself at eight (29.6 against 19.5 s in a window, with the same sizing). Its 1 × 1 synthesis and projection reach an OpenMP OpenBLAS, which on aarch64 sends them to its threaded GEMV at the process width; between the team’s eight-thread regions, GNU OpenMP ended and re-created eight threads at every call: a short window made 100,836 clone3 calls, against 72 for main and 28 at eight threads. Every native call of a compact problem already opens a team scope (CpuTeamScope); for an OpenMP OpenBLAS the scope now also sets the calling thread’s OpenMP default width to one, which OpenBLAS reads before it threads a call. openblas_set_num_threads is not used for this: it frees the per-thread buffers above the new count and allocates them again when the count returns, which cost h2_raman_averaged 3 % at 16 threads on an Intel Xeon Gold 6240. The library is recognised at import whatever width the environment gave it (the suite’s workers set OMP_NUM_THREADS, which skips fit_blas_threads). On the C4A gnlse_raman_shock now takes 17.6 s at 16 threads, against 25.8 s for main; on an Intel Xeon Gold 6240 the change is within 1 % of the A/A spread on every case.

One-mode products on the response team

A one-mode capillary with a field plasma column (argon_blueshift: 33 radial nodes, 6400 fine samples) synthesises its node fields as an outer product (33 × 1 by 1 × 6400) and projects three node drives a source evaluation as a stacked 2 × 33 by 33 × 6400 product. The product sizing keeps them serial on the engine’s side (211200 work, below the dense threshold), so the vendor GEMM threaded each call over its own pool: at 16 threads on a Xeon Gold 6240 the 4212 calls of a --length 0.1 window took 0.46 s, 108 us each, against 0.29 s at one thread; the band column, which writes the node rows, runs on the full-width response team. Plain loops tried earlier ran on the product’s serial instance and lost to the threaded library. On a host build with a vendor BLAS, one real mode at several nodes in one chunk now synthesises and projects through two small kernels on the response team: the projection one contiguous share of the samples per thread, both stacked rows from one read of the node rows, written straight into the fine and coarse rows (no stacked copy); the synthesis one node row per item. Each element sums the nodes in order from zero (within the GEMM rounding bound; a native test holds both kernels to KokkosBlas::gemm). Label timer at 16 threads: projection 40 us a call instead of 112, synthesis 26 instead of about 100; the band column and the held products that read the synthesised rows also ran faster (0.35 to 0.29 s, 0.10 to 0.06 s). Interleaved benchmarks ab against the compact-team commit, Xeon, 3 rounds, --length 0.1, the same steps, fields within 4.4e-15:

case	1 thread	16 threads
argon_blueshift	0.999 [0.993, 1.000]	0.679 [0.672, 0.688]
rdw (four modes, unchanged path)	1.000	1.000

argon_blueshift now runs 4.1 times faster at 16 threads than at one (1.01 against 4.13 s), where it was 2.8 times. On one thread the kernels are left to the vendor GEMM’s serial ones: in two tracking-set runs the one-thread window was 1.8 % slower with them, and without them it took 0.963 [0.948, 0.969] of the time with them (the four-mode rdw, whose code is unchanged, also moved by 0.964 in that A/B: code-layout noise of this size is present at one thread).

The kernels apply with OpenBLAS only. Apple’s Accelerate runs these small calls well, and on an 18-core Apple M5 Max the kernels on the 18-thread response team took the argon_blueshift window from 3.0 to 3.8 s (1.31 of main at 18 threads in the tracking-set A/B; a bisection of the push found this commit).

Stacked projection splits

A real projection that also forms the embedded rule’s coarse rows projects both in one stacked product and then copies the halves into the fine and coarse rows. The copy ran on the temporal copy’s instance, which the copy sizing keeps serial below 65536 samples: on the four-mode capillary (rdw, 4 × 10240 samples) each took about 90 us at 16 threads on a Xeon Gold 6240, 43 ms of a --length 0.1 window, beside full-width kernels; on ctc 159 us. It now runs on the stage’s instance (full width there, serial where the stage is). Interleaved benchmarks ab against the one-mode commit, 3 rounds, --length 0.05, the same steps, fields within 6.1e-15:

case	1 thread	16 threads
ctc	0.990	0.977 [0.976, 0.980]
rdw	0.987	0.992 [0.962, 0.994]
argon_blueshift_envelope	1.004	1.012 (no stacked real split)

NumPy’s BLAS beside a native run

NumPy’s OpenBLAS is a pthreads library with a pool of its own, sized from OPENBLAS_NUM_THREADS or every CPU (the benchmark suite sets the variable to the run’s width). During a native run the host calls it only for small products: the power weights’ expansion at each saved plane is a 3 × 2561 dgemv on averaged helium, 24 of them in a --length 0.5 window. After each call OpenBLAS’s workers spin before they sleep, and the planes come often enough that they never do: at 16 threads on a Xeon Gold 6240 Gold 6240, its 15 workers took 1.09 s of CPU each in a 1.24 s window, and the four-thread Adams passes ran three times as slowly (0.74 s against 0.22 s, label timer). The native backend now holds every BLAS but its own to one thread while it lives, as numba_cpu holds its BLAS pool. Interleaved benchmarks ab against the parent, 3 rounds, --length 0.1, passive-bound placement, the same steps and fields to rounding:

case	Xeon Gold 6240, 1 thread	16 threads	AMD EPYC 7763, 1 thread	16 threads
rdw	0.996	0.462 [0.456, 0.465]	1.012	0.762 [0.749, 0.780]
he_rdw_averaged	1.011	0.693 [0.670, 0.773]	0.997	0.343 [0.326, 0.816]
radial_filament_256	0.996	0.725 [0.718, 0.726]	0.992	0.610 [0.603, 0.791]
the other nine tracking cases	0.984-1.019	0.979-1.010	0.988-1.020	0.972-1.017
geometric mean (12 cases)	1.000	0.886	1.001	0.855

The EPYC rows are 16 cores of a shared node (its averaged-helium spread is the node’s).

The four-mode capillary now runs 3.1 times faster at 16 threads than at one (0.77 s against 2.40 s), where it was 1.4 times.

Smooth grid sizes

native/benchmarks/size_cliffs.py times the four-mode RDW capillary per source evaluation (short and long capped Runge-Kutta windows, their difference) on 2-3-5-smooth and power-of-two grids, and flags a size whose cost per sample exceeds its nearest power of two’s by more than 25 %. Two size boundaries made cliffs: the stage threshold (8192 elements, just above n_t 4000’s 8004-element stage, whose entry kernel then took 708 us serially against 59 us at full width) with the FFT threshold (32,768 samples, above its 32,000-sample batch), and a node budget that left 33 nodes as a 32-node pass and a serial one-node tail at n_t 8640. With the thresholds at 4096 and 16,384 and the nodes spread evenly over the passes, on one 18-core socket of a Xeon Gold 6240 (GCC 16, FFTW 3.3.11, ms per evaluation, two rounds each):

n_t	4000	4096	4500	5000	6000	7200	8000	8192	8640
before	2.45–2.47	1.24–1.36	1.49–1.56	1.54–1.57	1.86–1.87	2.21–2.26	2.38–2.39	2.29–2.30	3.52–3.55
after	1.32	1.27–1.36	1.53–1.54	1.57–1.61	1.83–1.87	2.17–2.26	2.38–2.40	2.28–2.31	2.54–2.58

At a stage threshold of 1024 a 1351-element SSFS stage regressed (1.2 to 1.7 ms); the reduced comb-to-continuum and SSFS cases did not move at 4096. On a dual EPYC 7643 (Zen 3, idle node, pinned) the same change made n_t 4000 24-27 % faster at 12 and 24 threads (2.36-2.42 to 1.74-1.82 ms), and the even chunks made n_t 8640 4-13 % faster; other sizes moved within 5 %. On an 18-core Apple arm64 laptop, in contrast, n_t 4000 went from 0.88 to 1.04 ms (level with 4096), and even 4096’s threaded stage is dearer than 4000’s serial one: the crossover there lies above 8192. macOS therefore keeps the stage and FFT starts at 8192 and 32,768; its own crossover still needs a sweep on a quiet machine. Its node-chunk gain stands (n_t 8640: 2.35 to 1.73 ms).

Radial analytic chunk dispatch

Admitting partial analytic radial node chunks to the interleaved real GEMM, instead of a scalar contraction that took about 94 % of the source time, made the radial filament’s adaptive Adams window (256 and 512 nodes, 12 accepted steps, construction excluded) 2.8 to 8.8 times faster on Kokkos CPU, on an exclusively allocated two-socket AVX-512 server. At one thread Kokkos then took 0.91 (256 nodes) and 0.95 (512) times numba_cpu’s time, at 36 threads 0.10 and 0.12 times.

Multimode CPU step costs

Sharing real cubic monomials across output components, planning general-modal tail transforms over active rows only and seeding serial transforms with measured FFTW wisdom brought Kokkos CPU to 1.00–1.07 times numba_cpu’s best time at one thread on an arm64 laptop, over 1–2 mm entrance windows of the general modal and reduced comb-to-continuum cases with both adaptive methods. In ten-step profiles the general route’s tensor work fell from 0.44 to 0.15 s and its ordinary FFTs from 1.00 to 0.84 s. Measured FFTW planning adds about 18 ms (general modal) and 120 ms (reduced comb-to-continuum) of construction. A direct FFTW execution adapter, a pooled allocator and vector transcendental functions were measured and not needed (CPU performance).

Threaded FFT plan selection

FFTW’s measured planning sometimes chose a slower threaded plan: on a 32-core AVX-512 server, for 204 rows of 2048 complex samples at 32 threads, a 64 × 32 decomposition instead of 32 × 64 made the warm general-modal RK window about 12 % slower, almost all of it in complex transforms. Measured wisdom is now seeded only for serial transforms; threaded ones default to estimate planning. Against the previous build, the general modal case’s window then took about 0.80 times as long at one thread and the same at 32 threads; reduced comb-to-continuum stayed within its repeatability.

Low-storage interaction-picture RK

The low-storage march (four evaluations a step) against the IP-DOPRI5 replay (six, with first-same-as-last) on a Cartesian envelope Kerr beam with no collar, fixed 0.5 µm steps, ten timed steps after one untimed one, same-process comparisons on an 18-core arm64 workstation. Shape A is 192 × 192 × 300 (one state 177 MB), shape B 384 × 384 × 600 (1.42 GB). Ratios are how many times less time the march takes than the replay on the same backend:

backend	shape, threads	per step	per evaluation	field volumes held, march / replay
numba_cpu	A, 1	1.83×	1.22×	4 / 16
numba_cpu	A, all	1.77×	1.18×	4 / 16
numba_cpu	B, all	1.73×	1.15×	4 / 16
Kokkos CPU	A, 1	2.02×	1.35×	4 / 15
Kokkos CPU	A, all	1.81×	1.21×	4 / 15
Kokkos CPU	B, all	1.85×	1.24×	4 / 15

Recomputing the flow factors instead of tabulating them holds three volumes and costs 3–13 % more. Most of a march step is now the source evaluation’s core (transforms, the Kerr column, packing): 69 % on numba_cpu at one thread, the fused finish 27 %. At all threads the Kokkos march was 1.6 (shape A) to 2.1 (shape B) times faster than numba_cpu’s.

Band-rule plasma integration

The band-exact plasma integrals (SEP 18) against the trapezoid, on the four-mode RDW window (33 nodes, 8192 fine samples, ten fixed steps from the launch), for the band column as first built:

backend	threads	band / trapezoid window
numba_cpu (arm64 laptop)	default	1.51–1.58
numba_cpu (arm64 laptop, or the AVX-512 server, where this column ran serially)	1	2.4–2.5
Kokkos CPU (arm64 laptop)	18 / 4 / 1	1.3 / 1.9 / 2.7
Kokkos CUDA (data-centre NVIDIA GPU), direct halo reads	—	1.10
Kokkos HIP (AMD Instinct GPU), halo staged in team scratch	—	1.25

On a CPU each integral is a nine-tap filter inside a serial traversal; on a device the filters are parallel maps over the row, so the band rule costs about a tenth more than the trapezoid there. Staging the halo in team scratch loses on the NVIDIA device (1.21) and wins on the AMD one (1.36 with direct reads), so each build defaults to its winner. numba_cpu’s column has since been rewritten as whole-row passes, 1.7 times faster than the traversal measured here (algorithms.md). A finer grid, the alternative SEP 18 rejects, costs about twice the whole nonlinear term for a quarter of the error.

Dispersive electronic response

The dispersive Kerr response (SEP 6, backends.md §3.11) is judged by its nonlinear stage time against the constant cubic’s on the same route: the symmetric evaluator’s target is 1.3, the polynomial’s 4. Same-process medians of ten-step windows from the launch, on the RDW grid and on Cartesian, line, planar and general-route problems:

evaluator	numba_cpu	Kokkos CPU	GPUs (NVIDIA and AMD)
symmetric, one species	0.92–1.23 (1.31–1.35 on one averaged node, analytic sector)	0.81–1.21 at one thread	0.86–1.30
symmetric, mixture (one filter from the mixture’s Kerr diagonal)	0.90–1.11	0.97–1.21 at one thread	0.86–1.30
polynomial, one species or equal leg bounds	1.16–3.72	1.52–3.65 at one thread; 4.2–5.0 on the general route	1.01–3.53
polynomial, mixture of distinct leg bounds	4.1–7.6 (3.3–3.7 for helium and argon in the default pools)	3.2–8.3 at one thread, at most 2.2 in the default pools	within 4 (air 1.5–2.1)

Ten-step windows, which carry the integrator’s own work, move less than the stage. A rebuild of the mixture filter when the densities change costs about 2 % of a symmetric stage. A polynomial mixture is one group per distinct leg bound, which is why mixtures such as air miss its target on the CPUs.

Far-field statistics pass

The banded far-field pass (six aperture statistics and the on-axis spectrum, 30 slots) on the last state of a Cartesian envelope Kerr run (192 × 192 × 256 and 384 × 384 × 512), same-process comparisons: about 2.3–2.6 % of one Runge–Kutta step on both numba_cpu and Kokkos CPU, paid at save planes and periodic rows only. On numba_cpu it costs about as much as the energy statistic and is 30 to 37 times faster than the host pass it replaced (which cost 77–84 % of a step). On Kokkos CPU it is compute-bound (square roots and saturated tanh weights), 15–16 times energy.

The per-row profile and shape of kernel schema 45 (#446) left a set of circles’ path alone. Measured on the pass alone, NativeBackend.stats of aperture_spectrum and beam_spectrum called seven times per process on a held 192 × 192 × 192-bin state of a Cartesian grid, interleaved old and new builds in four rounds, native CPU at 8 threads on an Intel Xeon Gold 6240, GCC 16.2. The machine was shared and each round’s numbers moved by up to 60 % between runs, so the fastest round of each arm is the figure (ms per call, before → after): three hard holes 31.75 → 32.16 (+1.3 %), three hard holes on one centre, with the radius shared, 32.26 → 32.69 (+1.3 %), three tanh holes 79.42 → 76.74 (−3.4 %), sixteen rcos holes 103.4 → 89.7 (−13 %). The results of both backends on these sets were identical to the bit. A shape other than a circle forms its offsets and takes one more branch per row and element. That comparison is indicative; the far-field pass runs at save planes and periodic rows only. The whole-run check is python -m benchmarks ab --case rdw --case far_field_stats --threads 1,16 --length 0.1 --rounds 5 --repeats 1, an exclusive AMD EPYC 7763 node, GCC 16.2, base f004a662: candidate/baseline far_field_stats 0.994 at one thread (21.48 → 21.38 s) and 0.994 at 16 (4.699 → 4.696 s), rdw 0.997 at both, geometric mean 0.995, the fields agreeing to 6e-15 and the same work on every row. Those holes are circles; the new shapes were not timed there.

On a device (#364): one stats call of aperture_spectrum on the marched field, rcos edges and quadrature weighting, the median of six in one process, against one low-storage RK4IP step of the same process. NVIDIA H200 with CUDA 13.1 and AMD Instinct MI325X with ROCm 10, 2026-11-04.

Grid	Holes	H200 call / ms	of a step	MI325X call / ms	of a step
384² × 600	4 circles	6.2	11.6 %	6.1	9.9 %
384² × 600	16 circles	13.6	25.3 %	14.6	23.9 %
768² × 600	4 circles	15.4	6.5 %	14.4	6.5 %
768² × 600	16 circles	33.3	14.1 %	33.8	15.2 %
768² × 600	4 mixed shapes	16.0	6.7 %	16.7	7.5 %
768² × 600	16 mixed shapes	36.9	15.6 %	44.7	20.0 %

A step takes 53.9 and 236 ms on the H200, 61.0 and 223 ms on the MI325X. So at a few steps per saved plane the pass is under 5 % of the plane, and SEP 23 §6’s registers and groups have little left to win: the pass costs about 1.5 ms per hole on top of the 9 ms beam_spectrum takes, and even sixteen holes cost less than the 54 ms (H200) and 80 ms (MI325X) of energy on the same field.

Raman kernel rebuild

On a density gradient the Raman kernel is rebuilt at every stage, and its cost scales with the response’s rows. numpy_ref’s rebuild at \(2^{14}\) fine samples (median of 35, same-process, a loaded laptop) took 0.31 ms for a two-term N₂ model (rotation and vibration), 0.47 ms for optmat’s H₂ (4 rows), 6.0 ms for optmat’s N₂ (72 rows) and 10.2 ms for dry air (126 rows), against 0.25 ms for one convolution. With tens of rows the rebuild dominates the step; Kokkos’s native rebuild is slower than the host’s at one thread and overtakes it at four. SEP 10’s low-rank representation is what removes it.

General-route symmetric contraction

Summing the general route’s real cubic over its ten symmetric monomials instead of the dense 81 + 3 products (algorithms.md), on a two-channel step-index pair (288 nodes, 256 samples) at one thread: the respond kernel alone 2.5 times faster on numba_cpu; a ten-step window 1.45 times faster on numba_cpu and 1.35 times on Kokkos CPU. The gain shrinks as channels are added.

Adams carrier frame per step

The scalar carrier frame (SEP 35) adds, per trial, a shifted argument and one complex product per element in the moments pass and one complex product per history row read. Ten replayed steps of a 2 mm, 40 bar hydrogen capillary (16 µm, 515 nm, 2 ps window, 4097 bins, one mode, a weak collar), the stepping alone timed, one thread, median ratio of nine same-process interleaved pairs on an Apple M5 Max, a machine in ordinary use (labelled A/B, not an absolute timing):

source	numba_cpu off → scalar	ratio	native CPU off → scalar	ratio
Kerr only	11.9 → 12.4 ms	1.042 [1.019, 1.096]	9.2 → 9.7 ms	1.057 [1.036, 1.066]
Kerr and Raman	121.5 → 121.9 ms	1.007 [0.995, 1.017]	95.4 → 95.8 ms	1.006 [0.997, 1.009]

Forming exp(x - sigma) as exp(x) exp(-sigma) instead of a second exponential took native’s Kerr-only ratio from 1.076 to 1.057. The per-step cost is a few per cent where the source is cheap and negligible where it is not; the policy’s case rests on the calls it saves, which SEP 35’s phase 4 measures at matched error.

Resident Adams storage

At maximum order five the Kokkos resident Adams engine held 29 complex field volumes, eleven of them dense-output slopes and saved fields allocated even when planes are landed on. Allocating those only for interpolate_planes=True reduces it to 18 (backends.md §8). Peak RSS of a reduced-resolution kerr_beam_384 (same box, 600 temporal bins, ten accepted steps, four threads, median of three fresh processes on an 18-core arm64 laptop), in GiB; the 384² row is a projection from the other two (constant plus transverse area), not a measurement:

grid	numba_cpu	Kokkos, 29 volumes	Kokkos, 18 volumes	18 volumes / numba_cpu
128²	4.31	7.53	5.92	1.37
256²	15.7	28.2	21.8	1.38
384², projected	34.8	62.7	48.2	1.39

Field export adds a further host copy (up to 1.43 times numba_cpu’s peak). No field-sized workspace scales with the thread count. The allocations are itemised under reference tables.

Free-space cases at the package defaults

Per-step time of the free-space cases at the package’s defaults, production thread policy, on a 64-core AVX-512 cloud server (seconds per step, medians of two repeats):

case	stepper	Kokkos CPU	numba_cpu
Kerr beam 384²	march	0.348	0.796
Kerr beam 384²	IP-DOPRI5	0.998	2.49
Kerr beam 768²	march	1.51	3.34
Kerr beam 768²	IP-DOPRI5	4.09	10.1
Cartesian filament 128²	default	0.076	0.300
Cartesian filament 256²	default	0.296	1.19

The march scales from 384² to 768² as 4.33 (Kokkos) and 4.20 (numba_cpu) against \(N \log N\)’s 4.30. Against a per-evaluation reference built from the cost of its parts (transforms and response), the march’s evaluations take 0.94–1.37 times the reference and IP-DOPRI5’s 1.21 (Kokkos) and 1.80 (numba_cpu) (SEP 28).

The boundary monitor’s reading

Each backend reduces the boundary monitor’s five numbers where its state lives, instead of copying the field to the host to be transformed there (ADR-057). The benchmark harness runs without a monitor, so the monitor’s cost is measured by its own same-process comparison: each case built by the harness’s builder and run through Simulation.run with the default check and with boundary_check=None, alternated, five times each after a warm-up of each, the shortest of each arm compared.

Monitored over unmonitored run time, native CPU, one thread, on a two-socket AMD EPYC 7763 node (GCC 16.2, FFTW 3.3.11, Kokkos 5.2.1 OpenMP); before is the parent revision, after is this change:

case, window	accepted steps	before	after
Kerr beam 192², the whole 2 µm guide	4	+7.8 %	+4.0 %
radial filament 256, 3 mm	68	+0.9 %	+0.6 %
GNLSE soliton, a fifth of the guide	81	+5.5 %	+1.5 %

The Kerr beam’s guide is shorter than its steps times the monitor’s 256 checks, so every step is checked. A check now costs one inverse transform of the state where it lives and a pass that sums it; it was three host copies and a single-threaded host transform. One check of the launched state, the backend’s reading against that host path, same process, shortest of seven:

case	EPYC 7763, one thread	NVIDIA A40	NVIDIA L4
Kerr beam 192² (36,864 columns × 384)	100 ms / 218 ms	3.4 ms / 214 ms	5.0 ms / 511 ms
radial filament 256 (256 × 864, field)	1.08 ms / 1.22 ms	0.13 ms / 1.03 ms	0.14 ms / 4.3 ms
GNLSE soliton (one column × 3072)	0.027 ms / 0.071 ms	0.082 ms / 0.085 ms	0.072 ms / 0.095 ms

The GPU columns are CUDA 12.8 builds (an A40 on its node’s host, an L4 on an 8-vCPU cloud host, whose host path is the slower). On a field grid under the split scheme the native reading uses the collar’s inverse real transform, half the analytic signal’s work; through the analytic signal the radial filament’s check took 1.83 ms against the host path’s 1.13 ms at one thread. The L4’s column was measured before that change. On mlx (Apple M5 Max, a busy laptop: ratios only, same process), monitored over unmonitored on the spatial cases of 60 fixed steps: radial 256 1.36 → 1.01, Cartesian 64² 1.21 → 1.05, Cartesian 128² 1.28 → 1.06. The unmonitored path is unchanged: python -m benchmarks ab at one thread on the same node gives 1.000 (GNLSE soliton), 1.005 (radial filament) and 1.005 (Kerr beam) for the three cases above.

What the split finds (the guards judged on the power each step’s split finds, halfway through its factor, and the split collar’s crossing limit, #127 and #149). The finish’s removal pass sums four more numbers of the samples it already holds, on every accepted step under the split scheme while a monitor is attached. Only a check reads them, but a saved plane is checked whatever the cadence said of its step, and the native marches decide their own steps; forming them only on the steps that end at a check would save a few multiply-adds per sample already in hand, at the cost of telling every march which steps those are. They are formed on every step. Measured on an Apple M5 Max (Homebrew Clang 23, NEON FFTW 3.3.11, one thread, a shared machine: ratios only), a tenth of each guide, before the change (its parent on main) and after:

case, a tenth of the guide	backend	monitored / unmonitored, before	after	ab, unmonitored
radial filament 256	native	1.028	1.030	1.018
line	native	1.666	1.660	1.001
planar	native	1.469	1.470	1.001

The first two columns are the same-process comparison above (shortest of five after a warm-up, alternated); the last is python -m benchmarks ab --length 0.1 --rounds 3, the median ratio of candidate to baseline over paired rounds. Within the machine’s spread both ways; the line and planar cases’ monitored cost is the transverse edge’s reads (#396). The steps are the same: on these cases the crossing limit does not bind.

The crossing cap read at every check (#440). Under the split scheme each check also sums the state’s energy per walk-off bucket, one more pass over the state where it lives (a histogram of 128 buckets: NumPy’s bincount on the host, a Kokkos ScatterView in native, a scatter-add on the mlx device). Measured as above, the same machine and window, before (the previous head of this branch) and after:

case, a tenth of the guide	backend	monitored / unmonitored, before	after
radial filament 256	native	1.031	1.039
line	native	1.658	1.698
planar	native	1.470	1.484

One to two and a half per cent of the run on the shared machine, about a third of a check’s own cost. The Cartesian filament 128² takes 1.158 monitored over unmonitored before; its after was not timed, at fifteen minutes a pair at one thread. The steps are unchanged on every case: the cap never fell below a step these runs took.

The transverse edge and a requested collar

SEP 40 puts a k⊥ safety ramp and reads of the transverse edge on every free-space run by default, and builds a real-space collar on request. The default’s cost was measured as the monitor’s was above, with each case built by the harness’s builder and run four ways in the same process:

• with neither the ramp nor reads;

• with the ramp alone;

• with reads alone;

• with both.

The four were interleaved, six rounds each. The setup: a Xeon Gold 6240, one socket, 18 threads, GCC 16.

case	steps (no reads → reads)	both, native	both, numba_cpu
line	27 → 28	+14.0 %	+6.4 %
planar	40 → 40	0.0 %	−1.1 %
radial filament 256	127 → 127	+2.2 %	−0.1 %
Cartesian filament 128, half length (its light reaches the edge)	40 → 41	+10.5 %	+2.6 %

The ramp alone costs nothing measurable. Reading costs in two ways. The read ceiling, min(L/16, 4g/v_q), adds a step on the line and Cartesian cases. Each read also has a host cost: about 1 ms on the line case’s 0.2 s run, and about 180 ms on the native Cartesian filament. Both are open in #396; the edge against 0.1.0 has the later measurements.

A requested collar was measured on the suite’s spatial_absorber case, a 0.9 P_cr argon beam over 8 cm with 64 samples across its 3e-4 m aperture. The table gives the field error on the aperture against a bare wall at 8 R:

transverse room	wall	error	translate	run, native	run, numba_cpu
3 R	bare	1.2e-7	—	0.22 s	1.24 s
3.6 R	collar, f = 2.6	2.3e-6	6.4 s	0.74 s	2.27 s
4.5 R	collar, f = 3.5	5.0e-7	8.2 s	1.08 s	2.24 s

Most of a collar’s cost is its design check at translation. That check also raised the process’s peak memory from 0.45 to 0.68 GB. The SEP’s phase 3 report has the full tables.

The transverse edge against 0.1.0

#453 found the tracking set’s kerr_beam rows slower than at 0.1.0, and ctc’s construction three times longer. These runs compare 0.1.0 with main at 00d83454 and with the changes this section describes.

The setup: python -m benchmarks ab on 16 cores of one NUMA domain of an exclusively allocated two-socket AMD EPYC 7763, SOLPROP_FFTW_PLAN=measure, with main’s harness on both arms. A small local change to benchmarks/builders.py, not merged, gave 0.1.0 the same grids, because it calls collar_fraction window_factor. The builds are GCC 16.2 from source for every 0.1.0 and fix comparison; the first sweep used the published 0.1.0 wheel and CI’s wheel, and agrees with them on the Kerr rows. Every row took the same steps in both arms and no sample was taken on a busy machine.

ctc’s construction is not a regression. Both wheels take about 5.5 s without the planning flag and about 15.5 s with SOLPROP_FFTW_PLAN=measure, construction_s including the backend’s construction, where FFTW plans. The 15.01 s #453 compared had the flag, and the 5.21 s row in the history is close to a run without it, so the two differ in planning, not in code. In the first sweep, with the flag on both arms, the ctc window’s ratio is 1.010. Limiting host BLAS during construction (#382) is not the cause: SOLPROP_BLAS_THREADS=16 made no improvement. --append-history now refuses SOLPROP_FFTW_PLAN, so the history’s rows stay comparable.

The Kerr rows’ cost is the transverse edge, with a smaller part from the split collar’s occupancy. With the monitor off, a Kerr run takes the same time in both arms, so the kernels are unchanged. Two first-parent boundaries account for the rise, with threads bound (--threading passive-bound, below): SEP 40 (05ee1f32 to 39ee3335) is ×1.110 on kerr_beam_192 at one thread and ×1.406 on kerr_beam_384 at 16; #439 (6083d63a to 3e13360e) is ×1.060 and ×1.015.

Two changes reduce it, neither moving a field by more than rounding:

• A read leaves out what cannot reach its resolution (READ_RESOLUTION, transverse_monitor.py). Each read synthesised every frequency plane, and a run’s light occupies from a fifth to a third of them at this level on the filaments and Kerr beams. A read now leaves out the planes whose energy together is at most 1e-3 of the warning. By Parseval a plane cannot put more on a face than it holds, so a face’s share reads low by at most 1e-9 of the launch energy at the default budget, against a warning at 1e-6.

• The occupancy takes its host scatter per row. scatter.access() finds its thread’s copy, and it dominated the loop when taken per element.

Complete runs, passive unbound threads as the history uses, three rounds, medians of the per-round pairs with the fix’s range. Each ratio is from its own comparison; the 0.1.0 times shown are from the comparison with main, and the fix’s comparison’s were within 1 % of them.

row	threads	0.1.0 (s)	main (s)	fix (s)	main / 0.1.0	fix / 0.1.0	range
kerr_beam_192	1	12.05	14.24	13.02	1.182	1.080	1.078–1.090
kerr_beam_192	16	3.17	3.75	3.75	1.185	1.184	1.179–1.184
kerr_beam_384	1	101.13	117.06	106.48	1.157	1.067	1.060–1.067
kerr_beam_384	16	25.33	29.91	28.83	1.180	1.136	1.129–1.140

Against main, the free-space set with threads bound, three rounds, same steps, final fields equal to 2e-15 or better except the Cartesian filament’s 8e-13 at 16 threads. Complete runs unless stated:

case	threads	main (s)	fix (s)	fix / main	range
kerr_beam_192	1	14.13	12.99	0.918	0.917–0.919
kerr_beam_192	16	4.56	3.95	0.871	0.868–0.980
kerr_beam_384	1	117.18	106.81	0.911	0.905–0.915
kerr_beam_384	16	35.77	29.73	0.831	0.830–0.832
radial_filament_256	1	6.85	6.78	0.990	0.979–0.992
radial_filament_256	16	1.88	1.89	1.003	0.992–1.007
line	1	0.905	0.889	0.982	0.978–1.000
line	16	0.275	0.267	0.980	0.962–0.987
planar	1	0.450	0.436	0.969	0.960–0.983
planar	16	0.148	0.149	1.012	1.006–1.047
cartesian_filament_128, half length	1	346.2	342.4	0.990	0.976–0.992
cartesian_filament_128, half length	16	97.3	96.3	0.990	0.986–0.995

Bound threads overstate the transverse edge. SciPy’s transforms take a bound main thread’s single core (0.48 s with 16 workers against 0.11 s unbound, on a 384² × 256 array), so main’s kerr_beam_384 at 16 threads is ×1.425 against 0.1.0 when bound and ×1.180 unbound. #453 measured ×1.152 unbound, from runs on different days.

What remains, proposed as an accepted cost. Each of these is a known part of the default’s work, and reducing it further means moving work rather than removing it:

• kerr_beam_192 at 16 threads does not improve. The monitor takes 0.77 s of 3.70 s against 0.20 s of 3.13 s at 0.1.0. Of the 0.57 s more, five reads are 0.31 s, of which the plane scan is 0.21 s, about what the left-out planes save at this width; six more host copies of the state are 0.12 s; the other 0.14 s, which includes the launch’s transform, was not broken down. Reducing the faces and the plane energies in native would remove most of this, which is the seam change in #396.

• These rows read at every step. A four-step run reads five times, its launch and each step, because L/16 is shorter than a step. radial_filament_256 took 19 reads over 127 steps in #396. How far the share falls on a longer Kerr run has not been measured.

• Construction. kerr_beam_192 goes from 1.63 s to 1.95 s and kerr_beam_384 from 11.25 s to 13.52 s, one-off work per translation. At 384² the launch’s occupancy is about 1.0 s, the walk-off buckets 0.6 s and the edge’s design 0.7 s.

• Memory. Peak resident memory rises by about 1 %: 2.51 to 2.54 GiB and 15.60 to 15.81 GiB.

• radial_filament_256, 3 to 5 % against 0.1.0. Five rounds, one and 16 threads:

boundary	t1	t16
0.1.0 to SEP 40’s parent	0.992	0.997
SEP 40	1.016	1.035
SEP 40 to #439’s parent	1.002	1.002
#439	1.028	1.016
#439 to main	1.003	1.006
main to the fix	0.990	1.003
0.1.0 to the fix	1.030 (1.018–1.035)	1.048 (1.019–1.070)

SEP 40 moves the final field by 1.2e-7 through its ramp. #439 takes no step more and leaves the field exact; its accepted steps now form five sums instead of one, so that the guards are judged on the symmetric splitting’s state. That the sums are its whole cost is read from the code, not profiled.

The rest of the tracking set. The 22 configurations, a quarter of each guide from the entrance, three rounds, 0.1.0 against the fix, threads bound: geometric mean 1.021. Apart from the Kerr and radial rows, the largest are argon_blueshift_envelope at 16 threads at 1.025 (0.12 s; 1.006 and 1.010 for a complete run on main) and grin_nonlinear at 1.017, whose range (0.999–1.020) overlaps its A/A control’s (0.979–1.015). The windows say nothing about what a run does past them. rdw at 16 threads, 1.048 in the first such window, was 1.001 at full length.

Host vector maths

The native CPU kernels’ transcendentals in vector passes (ADR-059) against the parent revision, complete runs of the tracking cases from the entrance (the fraction of the guide in brackets where it is not the whole), the default stepper and RungeKutta(), three rounds of interleaved fresh processes, medians with their median absolute deviations. Every arm took the same steps; final fields differ from the parent’s by 4e-15 to 6e-11 relative, the run-to-run spread of FFTW’s measured plans (a candidate with this change’s arithmetic on the libm route differs as much).

Two-socket AMD EPYC 7763 (Zen 3, AVX2; 16 cores of a shared node, glibc 2.34, GCC 16.2, FFTW 3.3.11, OpenBLAS 0.3.34, Kokkos 5.2.1 OpenMP, a local build, so SLEEF’s AVX2 functions), one thread, this change’s final revision, two rounds:

case	stepper	before (s)	vector maths
gnlse_raman_shock (0.3)	RK	18.84 ± 0.02	−19.4 %
argon_blueshift_averaged	RK	2.947 ± 0.015	−17.2 %
argon_blueshift_averaged	default	8.235 ± 0.003	−6.2 %
he_rdw_averaged	RK	1.365 ± 0.003	−12.4 %
he_rdw_averaged	default	2.338 ± 0.009	−10.5 %
h2_raman_averaged (0.5)	RK	2.866 ± 0.002	−10.1 %
rdw (0.3)	RK	7.679 ± 0.012	−14.4 %
rdw (0.3)	default	5.313 ± 0.012	−11.4 %
radial_filament_256 (0.3)	RK	8.915 ± 0.010	−7.8 %
radial_filament_256 (0.3)	default	6.261 ± 0.025	−5.9 %
argon_blueshift_envelope	RK	1.665 ± 0.006	−5.7 %
geometric mean			−11.1 % (conservative 0.896)

An earlier revision, three rounds, with a third arm built with SOLPROP_NATIVE_VECTOR_MATH=SCALAR (this change’s arithmetic on the libm route, an A/A-like control): that arm’s geometric mean was +0.7 % and its cases −1.2 to +2.2 %.

The free-space routes, the same machine and protocol, one thread, with the low-storage flows and the split factor: spatial_absorber (0.5) −12.2 %, radial_filament_256 (0.3) RK −5.3 % and default −1.8 %, planar (0.5) −1.5 %, line (0.5) −1.1 %, kerr_beam_192 (0.5) −0.9 % (inside its spread); geometric mean −3.9 %. Sixteen-thread runs on that node, which other jobs shared, scattered by more than the effect (one case’s arms spread by 30 %), and are not reported.

AMD EPYC 9B45 (Zen 5, AVX-512; a cloud instance of 16 cores without SMT, glibc 2.39, GCC 14.2, Ubuntu’s FFTW 3.3.10 and OpenMP OpenBLAS 0.3.26), local builds at -march=native; the AVX2 rows with -mno-avx512f too, so the AVX2 base is the reference. One thread, three rounds:

case	stepper	before, AVX2 (s)	vector maths, AVX2	AVX-512F (opt-in)
gnlse_raman_shock (0.3)	RK	10.42 ± 0.03	−21.3 %	−29.2 %
argon_blueshift_averaged	RK	1.581 ± 0.003	−10.9 %	−15.8 %
argon_blueshift_averaged	default	4.424 ± 0.004	−6.2 %	−13.1 %
he_rdw_averaged	RK	0.765 ± 0.011	−6.2 %	−8.7 %
he_rdw_averaged	default	1.295 ± 0.014	−10.8 %	−18.1 %
h2_raman_averaged (0.5)	RK	1.546 ± 0.003	−10.3 %	−14.7 %
rdw (0.3)	RK	3.873 ± 0.012	+4.9 %	+2.2 %
rdw (0.3)	default	3.123 ± 0.001	−2.5 %	−6.6 %
radial_filament_256 (0.3)	RK	5.026 ± 0.012	−4.3 %	−6.1 %
kerr_beam_192 (0.5)	default	5.071 ± 0.036	+2.0 %	+0.5 %
geometric mean			−6.8 % (conservative 0.945)	−11.4 % (0.900)

The AVX-512 base was itself 2.0 % faster than the AVX2 base, so AVX-512F’s functions save 9.8 % against it. At sixteen threads, the AVX2 builds, two rounds, six cases: geometric mean −5.3 % (conservative 0.983), argon_blueshift_averaged −10.9 % (RK) and −11.3 % (default), he_rdw_averaged −6.9 % (RK) and −1.5 % (default), rdw −0.9 % (RK) and +0.2 % (default). An Intel Xeon of the Ice Lake generation (a cloud instance of 16 cores without SMT, the same software), one to two rounds and noisier: the AVX2 functions −3.7 %, the AVX-512F functions +2.9 % (he_rdw_averaged with RK +21 %, rdw with RK +12 %). The AVX2 functions are therefore the default, AVX-512F an opt-in.

rdw with RK lost 3 to 5 % on both cloud machines, against a 7 % gain on the EPYC 7763, whose glibc (2.34) is older and slower. The label timer put it in the band plasma column (1.92 s → 2.21 s), whose block passes staged every sample’s magnitude, rate and fraction: most of the pulse’s samples lie outside the rate table and cost about a nanosecond each, as much as the staging. Writing the rows as the classifying loop goes and gathering only the live samples for the vector calls took the Zen 5 machine’s rdw with RK to −4.7 % (from +4.4 %), argon_blueshift_averaged to −14.1 % (RK) and −6.8 % (default) and he_rdw_averaged with RK to −11.4 %, three rounds; the tables above predate that change.

Apple M5 Max, a machine in ordinary use (labelled A/B), one thread: this change’s default (the system libm, each element’s own calls) against the parent, geometric mean +0.4 % over eight of the cases above, every case within its spread; SLEEF’s AdvSIMD functions forced on (SOLPROP_NATIVE_VECTOR_MATH=SLEEF) +11.6 %, gnlse_raman_shock +56 %.

Exponential Adams host cost: the SEP 39 baseline

The host side of an exponential Multistep() trial on the native CPU backend, before SEP 39’s native march, at the head of the vector maths change (a500d2d). The SEP’s §3.3 reads these figures.

Method.

• Complete runs of the tracking cases, with the fraction of the guide in brackets where it is not the whole, timed by the suite’s worker after its warm-up.

• The default stepper and RungeKutta() interleaved over three rounds of fresh processes, the order reversed every round; medians reported.

• Host time per trial: python -m benchmarks.adams_host CASE [--length F] (SEP 39), the march’s wall time less its native adams_* calls, with no other wrappers.

• The per-piece split comes from one run with perf_counter wrappers on about twenty-five host functions, which inflate the total by about half, so it gives the shares, not the sum.

• Steps and calls are counted.

AMD EPYC 7763 (Zen 3, 16 cores of one NUMA domain of a shared node, GCC 16.2, SLEEF’s AVX2 functions, FFTW 3.3.11, OpenBLAS 0.3.34), one thread:

case	Multistep() (s)	RungeKutta() (s)	host / native per trial (µs)	steps (acc/rej), calls: MS · RK
he_rdw_averaged	2.126	1.528	212 / 1 182	1 434/3, 1 438 · 480/1, 2 887
argon_blueshift_averaged	9.17 (7.83–9.47)	2.460	177 / 632	11 166/2, 11 169 · 1 684/0, 10 105
h2_raman_averaged (0.4)	4.820	2.146	173 / 695	5 417/2, 5 420 · 1 308/0, 7 849
rdw (0.25)	4.157	5.440	372 / 9 245	393/1, 395 · 175/0, 1 051

The node’s sixteen-thread runs, which other jobs shared, were slower than one thread on several cases and are not reported.

Apple M5 Max (Homebrew Clang 23.1.3, FFTW 3.3.11 with NEON codelets, the system libm, an otherwise idle machine):

case	threads	Multistep() (s)	RungeKutta() (s)	host / native per trial (µs)
he_rdw_averaged	1	0.815	0.650	51 / 448
argon_blueshift_averaged	1	2.920	1.419	38 / 213
h2_raman_averaged (0.4)	1	1.755	1.311	40 / 263
rdw (0.25)	1	1.711	2.736	129 / 3 730
he_rdw_averaged	16	0.636	0.650	53 / 282
argon_blueshift_averaged	16	2.474	1.381	36 / 175
h2_raman_averaged (0.4)	16	1.697	1.652	44 / 254
rdw (0.25)	16	0.862	1.256	205 / 1 555

The step and call counts are the EPYC’s on both machines.

Host pieces, helium, one thread (µs per trial, instrumented):

piece	EPYC 7763	M5 Max
stage arrays (bundle, packing, sample)	113	21
… of which the bundle’s Clenshaw reads	73	14
the decision, less the native accept	48	7
rate weights (uncached Clenshaw)	26	4
landing and approach	30	5
the tables call (one native helper)	18	4
order choice and stencil	17	3
sample, hoist, record and the rest	22	4

On the four-mode rdw the step history’s record costs 67 µs a trial on the EPYC: its backend hook scores the quadrature alarm natively.

Native time per trial (µs, host-only Kokkos label timer, one thread):

case	machine	predict pass	correct pass	the source and the rest
he_rdw_averaged	EPYC 7763	557	287	≈290
argon_blueshift_averaged	EPYC 7763	192	112	≈200
he_rdw_averaged	M5 Max	169	100	≈155
argon_blueshift_averaged	M5 Max	51	34	≈110

On the EPYC the two passes cost 1.5–2.9 source evaluations a trial, against 0.6–0.9 for the host. Their per-element cost is 2.9–3.8 times the M5’s, against about 1.9 times for the source.

The native Adams march

SEP 39’s native march against solprop’s Python march on the same build (benchmarks.adams_host --host-march), and against the march with the host packing the stage rows, as before the SEP’s phase 1 (--host-stages). Interleaved fresh processes, the arm order reversed every round, medians; every arm took the same steps and calls (the counts above). Host time per trial as in the baseline above.

Exponential formulation (the default Multistep()):

machine, threads	case	host µs/trial: packed → staged → native	wall, native vs staged	wall, native vs packed
M5 Max, 1	he_rdw_averaged	43 → 13 → 0.6	−3.3 %	−7.6 %
	argon_blueshift_averaged	30 → 11 → 0.4	−6.0 %	−13.6 %
	h2_raman_averaged (0.4)	37 → 11 → 0.5	−3.9 %	−11.6 %
	rdw (0.25)	91 → 47 → 2.7	−0.6 %	+0.4 %
M5 Max, 16	he_rdw_averaged	47 → 13 → 0.6	−4.2 %	−13.2 %
	argon_blueshift_averaged	32 → 12 → 0.4	−9.2 %	−17.2 %
	h2_raman_averaged (0.4)	39 → 12 → 0.6	−6.2 %	−16.2 %
	rdw (0.25)	171 → 122 → 1.9	−2.2 %	−4.5 %
EPYC 7763, 1	he_rdw_averaged	208 → 65 → 1.7	−6.1 %	−14.0 %
	argon_blueshift_averaged	148 → 58 → 1.1	−11.5 %	−24.2 %
	h2_raman_averaged (0.4)	168 → 58 → 1.2	−9.5 %	−22.7 %
	rdw (0.25)	168 → 8 (staged → native)	−1.7 %	—

The EPYC’s first three rows are from a shared node (load 80–97), whose one-thread walls were consistent; its rdw row is from a quiet node, nine rounds, where the sixteen-thread median read +2.7 % and the minimum −3.3 % inside a ±6 % spread.

Interaction formulation (Multistep(formulation="interaction"), M5 Max):

threads	case	host µs/trial: packed → staged → native	wall, native vs staged	wall, native vs packed
1	he_rdw_averaged	41 → 12 → 0.5	−4.0 %	−11.8 %
	argon_blueshift_averaged	33 → 11 → 0.5	−5.4 %	−16.7 %
	h2_raman_averaged (0.4)	40 → 11 → 0.5	−5.6 %	−18.3 %
	rdw (0.25)	85 → 43 → 1.5	−1.3 %	−2.9 %
16	he_rdw_averaged	41 → 12 → 0.5	−3.9 %	−12.9 %
	argon_blueshift_averaged	33 → 11 → 0.4	−8.6 %	−18.2 %
	h2_raman_averaged (0.4)	42 → 11 → 0.5	−34.1 %	−15.0 %
	rdw (0.25)	165 → 119 → 1.5	−2.1 %	−4.7 %

The staged Python march’s sixteen-thread h2_raman_averaged (2.24 s, slower than the packed march’s 1.74 s, in every round) is unexplained; see SEP 39’s phase 3 report.

On an NVIDIA A40 (CUDA 12.8.1; exponential formulation, short windows; the same steps in both arms): host time per trial 64 → 2.2 µs (he_rdw_averaged, 0.3), 61 → 1.6 (argon_blueshift_averaged, 0.1) and 83 → 5.0 (rdw, 0.1), walls −6.7 %, −10.6 % and −1.8 %. A trial costs about a millisecond on this FP64-weak device.

On an AMD Instinct MI325X (HIP 7.14 prerelease, AMD_GFX942; the same build, --host-march against the native march, five rounds, the same steps): host time per trial 21.7 → 1.5 µs (he_rdw_averaged, 0.3), 17.3 → 0.8 (argon_blueshift_averaged, 0.1), 40.5 → 3.0 (rdw, 0.1), 14.9 → 0.5 (h2_raman_averaged, 0.4) and 30.2 → 2.9 (radial_filament_256, 0.1); walls (ab --device, three rounds) −2.9 %, −3.4 %, −2.0 %, −3.0 % and −0.8 %, fields bitwise equal. A trial there is 160 to 550 µs of device time in 25 to 70 kernels, most of it launch-bound (device launches on an MI325X). On an NVIDIA H200 (CUDA 13.1, HOPPER90, three rounds, the same steps): 49 → 2.2 µs (he_rdw_averaged), 44 → 1.2 (argon_blueshift_averaged), 70 → 5.2 (rdw), 41 → 0.8 (h2_raman_averaged) and 65 → 4.7 (radial_filament_256, 0.1); walls −15 %, −21 %, −4 %, −23 % and −3 %.

The Adams passes on x86-64 lanes

The exponential Adams trial’s predict and correct passes on a host block of 32 elements held as lanes (include/solprop_native/adams_phi.hpp, adams_lanes_* in src/general.cpp) against the element-by-element form, at 740b24d. Before the change the passes ran each element through latency-bound scalar strict code (the phi-functions’ series, Horner pass and squarings, the compensation and the defects), GCC calling every strict helper out of line from the contracting kernels; per element they cost 2.9–3.8 times the Apple M5 Max’s on the AMD EPYC 7763 (§ “Exponential Adams host cost” above). Now every lane loop vectorises on x86-64-v3, Zen 3 and Cascade Lake (checked with GCC 16.2’s vectoriser report): the selections are bitwise, the comparisons quiet, and the corrector’s endpoint rate is formed beside the step’s phase, in one pass over the operator’s terms. The strict parts keep each element’s operations and their rounding (a native test holds the lanes’ phi-functions to adams_phi’s bits); the predictor’s and corrector’s sums may contract differently, so the first step whose inputs differ differs by about 2e-16 in its norm. Every case below took the same steps and trials in both arms; final fields differ by 2e-15 to 5e-12 relative (he_rdw_averaged at 0.3 of its length: 4.6e-12, the step sizes apart by up to 5e-5 relative after the pulse’s compression amplifies the rounding).

AMD EPYC 7763 (16 cores of a shared node, GCC 16.2, SLEEF’s AVX2 functions, OpenMP OpenBLAS), one thread, three interleaved rounds, medians; the passes in µs per call (host-only Kokkos label timer), the wall of the timed window:

case (fraction)	predict	correct	wall
he_rdw_averaged (0.3)	562 → 273	293 → 109	−32 %
argon_blueshift_averaged (0.1)	190 → 82	112 → 45	−30 %
h2_raman_averaged (0.2)	269 → 115	164 → 65	−37 %
rdw (0.1)	2430 → 1373	1222 → 530	−16 %
radial_filament_256 (0.05)	13 470 → 7 050	9 890 → 4 430	−16 %
argon_blueshift (0.05)	189 → 83	111 → 48	−7 %

Per element the predict pass now costs about 1.6 times the M5 Max’s figure of § 3.3 of SEP 39, and the correct pass about 1.1 times. A stage clock (rdtsc around each lane step, a profiling build) puts the rest of he_rdw_averaged’s predict pass at about a third in SLEEF’s six complex exponentials an element (the flow, the series argument and four spans), a sixth in the 18-term phase and rate sum, which is memory-bound, an eighth in the phi-functions and the rest in the history sums and stores. That node’s sixteen-thread runs scattered by more than the effect and are not reported.

Intel Xeon Gold 6240 (18 physical cores of one socket, the other socket in use, GCC 16.2 -march=cascadelake, SLEEF’s AVX2 functions, OpenMP OpenBLAS), complete windows at 0.1 of each case’s length, python -m benchmarks ab with passive bound threads, three rounds, medians of the paired ratios (range over rounds); every pair did the same work, final fields within 6.2e-15:

case	one thread	sixteen threads
he_rdw_averaged	0.786 (0.779–0.809)	0.879 (0.775–0.885)
argon_blueshift_averaged	0.752 (0.747–0.759)	0.946 (0.920–1.057)
h2_raman_averaged	0.697 (0.696–0.710)	0.936 (0.921–0.952)
rdw	0.827 (0.816–0.829)	0.974 (0.961–0.974)
radial_filament_256	0.839 (0.836–0.849)	0.937 (0.932–0.945)
argon_blueshift	0.936 (0.934–0.963)	0.953 (0.946–0.975)

The geometric mean of the twelve ratios is 0.867. At sixteen threads the single-column passes run on four threads (the size class above), so less of the gain shows.

SLEEF’s AVX-512F functions (SOLPROP_NATIVE_VECTOR_MATH_AVX512, opt-in) took a further 5–8 % off these windows at one thread on a Xeon Gold 6240, and were 6 % slower on he_rdw_averaged at sixteen threads (one round); they stay opt-in.

Tried and not adopted: blocks of 16 elements (5–15 % slower on the EPYC) and of 64 (within 2 % of 32 on single columns, 5 % faster on rdw).

Device transfers and the free-space Kerr stage

On an NVIDIA A40 (CUDA 12.8.1, AMPERE86, binary64 at about a sixty-fourth of the single-precision rate), python -m benchmarks ab --device, three rounds of interleaved fresh processes, half of kerr_beam_384 (two fixed RK4IP steps, 384² × 750):

change	wall, s	ratio (range)	field difference
one copy at a load and an observation, the free-space split on the stored columns	4.775 → 3.878	0.812 (0.812–0.815)	1.0e-15
the split factor per sample, the pinned single upload	4.686 → 3.892 (with the above)	0.830 (0.829–0.832)	1.0e-15
the rank-3 Kerr stage	3.913 → 3.764	0.963 (0.957–0.964)	1.3e-16
all, against the baseline	4.732 → 3.747	0.792 (0.789–0.798)	1.0e-15

Kernel time fell from 827 to 651 ms a step: the split’s transverse pair and transposes (about 100 ms) and the chunks’ pack, crop and pointwise passes (about 75 ms). An observation of the 1.77 GB state took 0.85 s of host time through three fresh buffers. On an Intel Xeon Gold 6240 at one thread kerr_beam_192 (half) took 0.920 of the time (0.917–0.922, fields to 6e-16). radial_filament_256 is unchanged on the A40 (0.3 and 0.5 of its length, within ±2 %): its trial, 3.3 ms of kernels, is bound by binary64 arithmetic (the synthesis and projection GEMMs near the device’s binary64 peak, the φ-functions), so neither the single upload, the GEMM for its node chunks nor one chunk in place of two (48 → 31 kernels a trial, not adopted) moved it. The retained-spectrum synthesis does, by halving the synthesis’s product: on the A40, half the run, three rounds of three repeats, 0.569 → 0.512 s, ratio 0.902 (0.898–0.903), fields to 3.2e-15, kernel time 238 → 200 ms over a 71-trial window. The device prefactor with one reciprocal in place of two divisions took the separable fold from 29.5 to 21.6 ms (kerr_beam_384, kernels 651 → 620 ms a step, wall 0.981, 0.936–0.984, fields to 2.4e-17).

Device launches and the Adams passes on an AMD Instinct MI325X

On an AMD Instinct MI325X (HIP 7.14 prerelease, AMD Clang 23, AMD_GFX942, band staging on), rocprofv3 --kernel-trace over a timed window showed the device busy for a third of an Adams run. A single-mode trial was 30 to 70 kernels of 2 to 9 µs, launched one by one, and between them Kokkos’s constant-memory launches: a functor of 512 B to 32 kB (most of a trial’s) is copied by a blit kernel behind an event wait on the previous such launch, 4 to 10 a trial, 13 to 18 % of the kernel time. Each change below, ab --device, interleaved fresh processes, three rounds, the same steps:

change	case	ratio (range)	field difference
functors under 3840 B passed as the launch argument (a light-weight hint)	he_rdw_averaged (0.3)	0.900 (0.885–0.935)	1.4e-11
	argon_blueshift_averaged (0.1)	0.829 (0.823–0.849)	3.7e-15
	radial_filament_256 (0.5)	0.938 (0.935–0.947)	3.4e-15
HIP range reductions back to 256-thread blocks with the hint	he_rdw_averaged	0.944 (0.931–0.953)	0 (against no hint)
predict pass’s phi-functions in registers (order dispatched, one to eight)	he_rdw_averaged	0.971 (0.949–0.984)	1.0e-12
correct pass’s history read once into registers	he_rdw_averaged	0.979 (0.973–0.995)	1.3e-11
	radial_filament_256	1.000 (0.983–1.002)	3.8e-15
no fence after an accepted step on a device	he_rdw_averaged	0.959 (0.955–1.005)	0
	argon_blueshift_averaged	0.966 (0.900–0.971)	0
the held currents’ transforms in one plan	he_rdw_averaged	0.964 (0.957–0.967)	0
	argon_blueshift (0.1)	0.903 (0.890–0.938)	0

The light-weight hint first took the range reductions to Kokkos’s largest HIP block (1024 threads), which put a single-mode trial’s predict and correct passes on three compute units, 40 → 56 µs; their reduction order moved the norms, hence the 1.4e-11 over 411 steps. The predict pass held its phi-functions in scratch memory (560 B a thread) because their loops ran to a run-time order: 38.9 → 28.8 µs a trial on he_rdw_averaged and 64.6 → 32.7 µs on radial_filament_256. The correct pass’s register copy is 31.9 → 27.5 µs on the former and 54.2 → 57.7 µs on the latter (its registers doubled); it is kept for the small runs. The unrolled device arithmetic contracts differently, so fields move at rounding level, which the adaptive steps carry to the 1e-12 and 1e-11 above with the same steps.

All of them, against the HIP build with the native march alone (five rounds; radial three repeats):

case	wall, s	ratio (range)	field difference
he_rdw_averaged (0.3)	0.196 → 0.145	0.744 (0.697–0.814)	1.3e-11
argon_blueshift_averaged (0.1)	0.346 → 0.237	0.686 (0.683–0.688)	2.9e-15
h2_raman_averaged (0.4)	1.344 → 0.994	0.740 (0.715–0.791)	1.3e-13
rdw (0.1)	0.196 → 0.171	0.875 (0.861–0.901)	8.1e-15
radial_filament_256 (0.5)	0.182 → 0.165	0.901 (0.889–0.907)	3.7e-15

A trial’s native time (benchmarks.adams_host, medians of five) went 323 → 206 µs, 279 → 178, 224 → 160, 446 → 291 and 547 → 399 µs. Against the previous HIP build, which stepped the runs with the Python march: he_rdw_averaged 0.463 (0.405–0.666; that build’s walls are bimodal on this device, 0.20 or 0.37 s), argon_blueshift_averaged 0.629, h2_raman_averaged 0.624, rdw 0.605, radial_filament_256 0.843, kerr_beam_384 (half) 0.650 and kerr_beam_192 (half) 0.702 (both from the free-space work above), and at a tenth of their length (two rounds) gnlse_raman_shock 0.843, argon_blueshift 0.657, ctc 0.780. What remains is mostly the host’s launch rate: a trial’s kernels take 2 to 4 µs each, and rocFFT runs each of these transforms as four kernels.

On an NVIDIA H200 (CUDA 13.1, HOPPER90) the same changes, of which the light-weight hint is a no-op there (grid-constant arguments), against the build without them (three rounds, the same steps): he_rdw_averaged 0.920 (0.904–0.936, fields to 1.2e-11), argon_blueshift_averaged 0.902 (0.880–0.912), h2_raman_averaged 0.918 (0.916–0.920), rdw 0.969 (0.962–1.035), radial_filament_256 0.989 (0.962–1.134). Over a 411-trial window of the first the correct pass went 21.4 → 14.1 µs and the predict pass 24.7 → 22.2 µs, the transforms 4 → 3 a trial (cuFFT’s one kernel each). Against the previous main build: he_rdw_averaged 0.790, argon_blueshift_averaged 0.882, h2_raman_averaged 0.788, rdw 0.841, radial_filament_256 0.868, kerr_beam_384 (half) 0.722, kerr_beam_192 (half) 0.776.

Pinned trial norms and batched forward rows on a device

On a device, an Adams trial’s defect norms are reduced into pinned host memory. The reduction’s last block writes them there directly, so the device-to-host copy is gone (on HIP it was a blit kernel: 2.16 → 1.16 copies a trial), and the trial keeps its one fence. A modal stage’s in-stage absorber, projection and held-current rows are carved from one buffer and transformed with one batched plan after the absorber is formed. This does not apply to dispersive responses, the general route or coarse rows. The radial route’s three held streams also take one plan.

All A/Bs use ab --device, three rounds, against the tree before the change. Every case took the same steps, and the fields were bitwise equal.

device	case	ratio (range)
AMD Instinct MI325X	he_rdw_averaged (0.3)	0.903 (0.856–0.906)
	argon_blueshift_averaged (0.1)	0.858 (0.854–0.929)
	h2_raman_averaged (0.4)	0.888 (0.836–0.893)
	rdw (0.1, coarse rows, not packed)	0.994 (0.959–1.001)
	radial_filament_256 (0.5)	0.969 (0.963–0.977)
	kerr_beam_192 (0.5, unaffected)	1.004 (0.984–1.024)
NVIDIA A40 (batched rows and pinned norms)	he_rdw_averaged	0.784
	argon_blueshift_averaged	0.755
	h2_raman_averaged	0.878

The pinned norms alone were within noise on the MI325X. When the pack is forced on a CPU build, the fields are bitwise equal.

Device graph replay of an Adams trial’s source

On a CUDA or HIP build, an elementwise Adams trial’s source evaluation runs in a graph region (graph::Region). The first trial captures its kernels, device copies and transforms. Each later trial whose operations repeat, in order and with equal keys, launches the instantiated graph instead. In a microbenchmark, thirty dependent tiny kernels cost 3.36 µs each back to back and 2.18 µs as a replayed graph on an AMD Instinct MI325X, and 2.41 and 1.67 µs on an NVIDIA A40. Capturing every trial and updating the instantiated graph cost more than plain launches (4.04 µs on the MI325X).

All A/Bs below use ab --device, three rounds, against the build without the region. Every case took the same steps, and the fields were bitwise equal.

case	MI325X	NVIDIA H200	NVIDIA A40
he_rdw_averaged (0.3)	0.915 (0.902–0.995)	0.991 (0.937–1.003)	0.972 (0.944–0.976)
argon_blueshift_averaged (0.1)	0.893 (0.880–0.897)	0.918 (0.917–0.962)	0.962 (0.959–0.965)
h2_raman_averaged (0.4)	0.885 (0.864–0.889)	0.964 (0.948–0.976)	0.978 (0.976–0.981)
rdw (0.1)	1.029 (0.978–1.031)	0.987 (0.968–1.039)	0.999 (0.997–1.000)
radial_filament_256 (0.5)	0.990 (0.986–0.994)	1.022 (0.992–1.114)	0.999 (0.997–1.004)
kerr_beam_384 (0.5)	1.001 (0.981–1.003)	1.001 (0.993–1.023)	0.987 (0.986–1.005)

The three mode-averaged runs replay all but their first trial: 410, 1048 and 5418 replays to one capture. The sources of rdw and radial_filament_256 hold a vendor GEMM or a reduction. Their region stops at its second operation and the run continues with ordinary launches; its one capture is the 3 % on the MI325X’s short rdw window. The Kerr cases do not use the region. A first version on the H200, whose trial is bound by the host, was slower (he_rdw_averaged 1.10, argon_blueshift_averaged 1.25). A verifying trial walks the whole source before its one launch, and that walk copied every operation it skipped. Comparing keys in place and copying an operation only when it is recorded removed that cost.

Device transforms of the four-mode capillary and the radial filament

rdw and radial_filament_256 gained least on devices from the batched rows and the graph replay. An NVIDIA A40 trace of one rdw trial (34 kernels, 1.30 ms of kernels, gaps of about 0.9 µs) showed why: the trial is bound by its transform shapes, not by its launches.

• Coarse rows in the forward pack. Each of the stage’s four real transforms (fine and coarse projection, absorber, held currents) took about 96 µs on 16 of the A40’s 84 multiprocessors, whatever its rows: the held pack’s 16 rows cost the same as the projection’s 4. The coarse projection’s rows now join the device forward pack, so the stage’s 28 rows are one transform (packs_forward). A device never defers the coarse quadrature, so those rows are projected at every stage.

• Few-mode transforms through the custom kernels. cuBLAS took 71 µs for the stacked projection (8 rows, 36 nodes, 10240 samples; about 84 GFLOP/s) and 27 µs for the synthesis. Below eight modes the device’s hybrid guard now keeps the custom kernels, which form the fine and coarse rows in one pass: 19 µs and 8 µs. With no vendor GEMM left, the whole source replays as a graph (158 replays to one capture).

• CUDA trial reductions in blocks of 128. Kokkos’s CUDA range reduction gave the prediction and correction 256-thread blocks of one element each: rdw’s 10496 elements were 41 blocks. Blocks of 128 occupy every multiprocessor. On the A40 the prediction went from 92.5 to 55.4 µs and the correction from 49.3 to 41.0 µs. The reductions’ order changes, so the trial norms move by rounding.

• An AMPERE86 build’s real GEMMs. cuBLAS runs a compute capability 8.6 device’s DGEMM with the FP64 tensor-op kernels of 8.0. A probe on the A40 measured 316 GFLOP/s for the radial filament’s 128 × 256 × 866 synthesis, about 54 % of the FP64 peak, at any row pitch. device_dense_gemm gives a thread four rows of one column and takes 132 µs where cuBLAS took 180. Eight rows a thread were 5 % slower. The radial source then replays as a graph too (70 replays to one capture).

All A/Bs use ab --device, three rounds, against the build before these changes.

case	NVIDIA A40	AMD Instinct MI325X
rdw (0.1)	0.799 (0.795–0.810)	0.962 (0.912–0.974)
radial_filament_256 (0.5)	0.930 (0.928–0.934)	1.000 (0.995–1.005)
he_rdw_averaged (0.3)	0.934 (0.927–0.937)	0.996 (0.990–1.004)
argon_blueshift_averaged (0.1)	0.940 (0.938–0.956)	0.994 (0.990–1.007)
h2_raman_averaged (0.4)	0.952 (0.947–0.954)	0.999 (0.989–1.025)

Every case took the same steps. Measured one change at a time on the A40, the coarse rows gave rdw 0.892 and the custom transforms 0.923, both bitwise. The CUDA reductions gave rdw 0.981 (0.968–1.008), and they are the whole of the single-mode runs’ gain. The AMPERE86 GEMM gave radial_filament_256 0.933, bitwise. The reductions’ order moves the norms by rounding, and through the step sizes it moves the final fields: by 2e-15 to 6e-14, and by 1.0e-11 on he_rdw_averaged (410 steps). The HIP reduction’s block did the same before (1.4e-11). On the MI325X only the first two changes apply. There, rdw’s field moves by 2.5e-15 (the custom transforms against rocBLAS) and the others are bitwise. The MI325X’s rdw source replays as a graph after four captures, while its radial source still holds rocBLAS.

Three leads were dropped. Eight rows a thread in the AMPERE86 GEMM ran at 1.054 against four. Padding the radial GEMM’s row pitch changed nothing: cuBLAS ran at 316 GFLOP/s with 866, 872, 896 and 1024 doubles. Recording vendor GEMMs in the graph was not tried. The trace puts its saving at about 1 % of a radial trial, and KokkosBlas sets cuBLAS’s stream at every call, which resets cuBLAS’s workspace to its default pool.

What is left on the A40’s rdw trial (0.77 ms of kernels): the inverse transform (96 µs on 4 multiprocessors, one block a row), the packed forward transform (96 µs), the band column’s four kernels (235 µs), the three projections (57 µs), the prediction (55 µs), the correction (41 µs) and the quadrature partials (39 µs on 4 teams). A timed run is also one long host gap at the first saved plane, in solprop’s output layer.

Split real-row transforms on compute capability 8.6

After the device transforms of the four-mode capillary, an NVIDIA A40 trace of one rdw trial still spent a quarter of its 0.77 ms of kernels in two cuFFT kernels: the inverse transform of the fine spectrum (4 rows of 10240 samples, 96 µs) and the packed forward transform (28 rows, 96 µs). cuFFT gives a real row of a few thousand samples one thread block, and a GA10x device (compute capability 8.6) runs FP64 at a sixty-fourth of its FP32 rate, so each block is FP64-bound on its multiprocessor while the rest of the device waits. A single-mode run’s one-row transforms used one multiprocessor of 84.

An AMPERE86 build now transforms up to 64 contiguous real rows of n = 2m samples as m-point complex rows (sample pairs as one point, the standard half-length real transform) by the four-step method (RealRowPlan, fft_plan.hpp): m = n1 n2, with n2 the largest power of two up to 512 that divides m. Our kernels pack the rows into n2 × n1 columns, cuFFT transforms n1 points per column, a kernel applies the twiddle exp(∓2πi j2 k1 / m) and transposes, cuFFT transforms n2 points per row, and a last kernel forms the real spectrum (or unpacks the samples). The roots are rounded once from long double. The inverse drops the imaginary parts of the DC and Nyquist bins, as cuFFT’s and FFTW’s do. A probe on the A40 timed 4 × 10240 at 20 µs against cuFFT’s 102, and 28 rows at 60 against 114, both agreeing with cuFFT to 2e-15 relative. Of the splits of 5120 points, 10 × 512 and 40 × 128 were fastest; those with an 80-point pass took half as long again. In the engine the two transforms of an rdw trial take about 78 µs where they took 193.

An NVIDIA H200 runs the same cuFFT kernels in 12 to 13 µs, faster than the split (12.5 µs for 4 rows and 17.8 for 28), and on it the single-mode mode-averaged runs lost 4 to 5 % with the split, so only AMPERE86 builds split their rows, as they form their own real GEMMs.

case	NVIDIA A40
rdw (0.1)	0.932 (0.924–0.939)
he_rdw_averaged (0.3)	0.751 (0.750–0.752)
argon_blueshift_averaged (0.1)	0.757 (0.756–0.758)
h2_raman_averaged (0.4)	0.770 (0.768–0.771)
radial_filament_256 (0.5)	0.996 (0.985–0.997)

ab --device, three rounds, against the build before the change. Every case took the same steps; the radial filament’s rows are too many to split and its field is bitwise. The fields agree to 1.2e-13, except he_rdw_averaged (1.1e-11): the split moves each transform by rounding, the trial norms follow, and through the step sizes the final field, as the CUDA reductions’ order did (1.0e-11). On the H200 the build with the split off agrees with the build before bit for bit and takes 0.99 to 1.01 of its time.

Two other levers were measured and dropped. A team per 128-bin chunk of a column’s quadrature partials took them from 39 to 7.8 µs, but the score’s sum over the chunks took most of that back, by one thread (2.5 to 7.9 µs) or by a team: 0.989 on rdw, 0.994 and 1.005 on the single-mode runs. The band column’s four tiled kernels (235 µs of the rdw trial) moved neither with teams of half a tile (0.985 on rdw, 1.002 on he_rdw_averaged) nor with 512-sample tiles (1.010, 1.019). The A40 we measured on gave no access to its performance counters, so what bounds them is not known. A single-mode trial’s prediction and correction (84 of its 205 µs of kernels) ran on 21 blocks of 128 threads; blocks of 64 gave 0.985, 0.980 and 0.987 on the mode-averaged runs and 1.000 on rdw, and blocks of 32 regressed the radial filament (1.029).

Retained radial synthesis and Raman modal products

Two work reductions from the CPU experiment’s handoff, each measured with python -m benchmarks ab (fresh interleaved processes, passive-bound threads, SOLPROP_FFTW_PLAN=measure, three rounds, one repeat; medians of the per-round pairs, range over rounds). Each case exercises one of them: the radial filament has no Raman term, and the capillary is not a square radial map. Work counts are equal in every pair and the final fields agree to 4–10 × 10⁻¹⁵ relative.

• Retained-spectrum radial synthesis (square radial analytic maps, backends.md): the synthesis product covers 2 n_omega real columns instead of 2 n_to. radial_filament_256, complete guide: on an Intel Xeon Gold 6240 (one socket, 18 cores) 0.928 [0.925, 0.945] at one thread and 0.954 [0.936, 0.968] at 16; on 16 cores of an AMD EPYC 7763, 0.869 [0.868, 0.891] at one thread, and at --length 0.05 0.866 [0.850, 0.868] at one thread and 0.825 [0.653, 0.894] at 16. The EPYC’s 16-thread ratios moved by about ±15 % between allocations on a shared node (a complete-guide repeat gave 0.999 [0.735, 1.020]); the Intel Xeon Gold 6240’s are the resolved ones. A 16-thread label profile on the EPYC put the synthesis GEMMs at 0.20 s against 0.38 s, of which the gather and pad, as first written per chunk and over every fine bin, took back 0.14 s; they now gather once per evaluation and write only the retained bins.

• Raman modal products (three-mode ctc, 6 product rows against 33 node drives): --length 0.05 at 16 threads 0.884 [0.881, 0.899] on a Xeon Gold 6240; on the EPYC 0.833 [0.823, 0.837] at one thread and 0.921 [0.861, 0.928] at 16 (--length 0.02: 0.850 and 0.884). The first reconstruction looped over the products per sample and cost ten times the node path’s finish (0.46 s against 0.046 s in a one-thread window); it now accumulates each product over contiguous blocks of 512 samples.

Both add FFT plans (the node transforms; the product rows), measured at construction: on a Xeon Gold 6240 a CTC worker’s wall time, construction included, rose from 16.3 s to 17.7 s.

Not adopted, on the same EPYC allocation:

• Routing every real product with a left matrix of at most 256 entries through native_order_gemm instead of the vendor BLAS: RDW 1.012 / 0.881, argon 1.057 / 1.099, CTC 1.070 / 1.377 at one / 16 threads (--length 0.05). The vendor library wins on the serial transform spaces the operation-size policy gives these shapes.

• Plain row-block loops for one-mode synthesis and projection: argon 0.996 at one thread; at 16 the loops took 0.65–0.69 s of a short window where the vendor GEMMs took 0.55 s. The vendor pool still costs 0.55 s against 0.23 s at one thread there.

• Planning the product rows at the full Raman chunk’s width rather than their own: CTC 1.001 / 0.980 [0.952, 1.012].

Host stage entry, free-space crossings and statistics

Host-only changes found by profiling the tracking set at one and sixteen threads (Kokkos label timer and perf, Intel Xeon Gold 6240, GCC 16.2 -march=cascadelake, SLEEF’s AVX2 functions, OpenMP OpenBLAS). Each was measured with python -m benchmarks ab (fresh interleaved processes on 18 cores of one socket, passive bound threads, SOLPROP_FFTW_PLAN=measure, three rounds of one repeat at 0.1 of each case’s length; medians of the per-round pairs, range over rounds). Every pair did the same work.

• The stage entry’s finish and combination as lanes (the modal routes’ host entry kernel). The deferred finish and the combination ran each bin through scalar Kokkos::complex code: on gnlse_raman_shock (16384 bins, DOPRI5, noise, Raman, a temporal collar) the entry took 2.8 s of a 7.1 s one-thread window, 21 % of it in the lambda’s own scalar code. A 64-bin block now runs as real and imaginary lanes, each bin the same operations in the same order (only the crop’s gathers stay per bin, and an identity crop reads and writes contiguously), and a lossy problem’s block whose phases all have zero real parts takes both magnitudes as exactly one without the two exp passes (only 42 of the 256 blocks of gnlse_raman_shock carry loss). One trial step from identical input differs by at most 3.4e-18 of the field’s peak.

  case	one thread	sixteen threads
  gnlse_raman_shock	0.836 (0.826–0.840)	0.962 (0.960–0.987)
  gnlse_soliton	0.955 (0.937–1.042)	0.939 (0.934–0.959)
  argon_blueshift_envelope	1.011 (0.968–1.019)	0.984 (0.981–1.013)
  ctc	1.012 (0.982–1.012)	1.022 (1.001–1.027)

  Both arms of these sixteen-thread runs predate the hold of NumPy’s BLAS.

• The separable low-storage fold in the scratch’s layout. The fold read the column-major state, accumulator, stage register and flow table 32 columns apart at each bin beside the transverse scratch, divided each column by the polarisation count and called the element’s full spatial value out of line: 63 ns an element at one thread, 27.7 s of kerr_beam_384’s 69.8 s window. The host’s registers now sit at the scratch cell’s own position; a work item takes 32 adjacent columns over every bin, loads each 16-bin tile of their state a contiguous run at a time and folds it in the scratch’s order; a tabulated fold reads only the prefactor. The flow table is formed in the same order. Fields are bitwise equal.

  case	one thread	sixteen threads
  kerr_beam_384	0.789 (0.788–0.790): 68.6 → 54.2 s	0.899 (0.895–0.903)
  kerr_beam_192	0.742 (0.732–0.746)	0.902 (0.900–0.918)

• Crossings without divisions, the monitor’s runs, the radial statistics’ GEMMs. Every host crossing between state columns and the transverse scratch located its cell with a run-time % n_pol and / n_pol; in a microbenchmark of kerr_beam_384’s chunk (4468 rows of 750 bins) at one thread they doubled the crop’s time (17 to 35 ms), and visiting the pack’s tile rows outermost took it from 36 to 22 ms. The boundary monitor’s reading divided each sample’s flat index by the row length (1.76 s of kerr_beam_384’s one-thread window); a host work item now sums a run of 512 samples of one row. A radial run’s local flux statistics synthesised each node’s field and closure with a scalar 256-mode sum inside the pointwise kernel (0.57 s of radial_filament_256’s 3.65 s window, five planes); they are now two GEMMs over the time fields’ interleaved real columns into buffers the statistics already hold, bitwise equal on the host. Against the fold above:

  case	one thread	sixteen threads
  radial_filament_256	0.834 (0.827–0.840)	0.945 (0.933–0.954)
  kerr_beam_384	0.830 (0.826–0.856): 53.7 → 44.9 s	0.943 (0.943–1.190)
  kerr_beam_192	0.855 (0.849–0.861)	0.951 (0.946–0.973)

  Together with the fold, kerr_beam_192 takes 0.635 of its time at one thread (8.04 → 5.11 s) and kerr_beam_384 0.847 at sixteen (9.36 → 7.93 s); at sixteen threads the 384² fold is now bound by memory bandwidth (about 52 GB/s), and about 5 s of that window is solprop’s own host work at the saved plane.

On 16 cores of an AMD EPYC 7763 (GCC 16.2 -march=znver3, the three changes together, same protocol), one / sixteen threads: gnlse_raman_shock 0.841 (0.836–0.842) / 0.972, gnlse_soliton 0.941 / 0.942, radial_filament_256 0.838 (0.837–0.850) / 0.938, kerr_beam_384 0.953 (0.947–0.958) / 1.003, kerr_beam_192 0.956 / 0.995, ctc 0.994 / 0.977, argon_blueshift_envelope 0.981 / 1.001. The free-space gains are small there: Zen 3 ran the 384² window in 30 s at one thread where a Xeon Gold 6240 took 69 s, its 64-bit divisions and strided streams far cheaper.

Not adopted:

• The Kerr-only free-space stage through one rank-3 FFTW transform, as a device takes it: for 384 × 192² at one thread an estimated plan took 1.03 s a transform and a measured one 0.31 s, the same as the batched planes and a strided temporal transform, after 17 s of planning; native’s transform pair with the chunks’ pack and crop is no slower.

• The pack and crop through a 32 × 32 tile buffer (as the fold’s state): 0.55 → 0.68 s and 0.52 → 0.56 s on kerr_beam_192 at one thread.

• A one-element GEMM (one mode, one node) as an elementwise product: no faster than OpenBLAS’s small kernel; both are memory-bound.

• Transparent huge pages through glibc’s glibc.malloc.hugetlb=1 (the workstation’s setting is madvise): within ±3 % on kerr_beam_192 at one thread.

Persistent FFTW wisdom

At one thread, the default plans every row transform with FFTW’s measured search, which is most of a short run’s backend construction: about 1.3 s for averaged helium’s 6400-sample real pair alone. SOLPROP_FFTW_WISDOM now keeps those plans for the next process (tuning). Fresh processes at one thread, SOLPROP_FFTW_PLAN=measure, passive-bound, through the benchmark worker’s own timed_run (--length 0.1; ctc, kerr_beam_192 and grin_nonlinear on the EPYC and ctc on a Xeon Gold 6240 --length 0.05). For each case one process stores the wisdom, then three rounds alternate a process without the variable (cold) and one with it (warm): medians of the backend’s construction, cold → warm, and the warm/cold window ratio’s median and range over the rounds. Step counts are equal in every pair. One core of an Intel Xeon Gold 6240 (an otherwise idle socket) and of a 16-core allocation of an AMD EPYC 7763:

case	Xeon construction	Xeon window	EPYC construction	EPYC window
he_rdw_averaged	1.398 → 0.033 s	0.985 [0.945, 0.988]	1.777 → 0.049 s	1.002 [0.973, 1.004]
argon_blueshift_averaged	1.330 → 0.032 s	1.006 [0.987, 1.009]	1.708 → 0.040 s	0.987 [0.961, 1.001]
h2_raman_averaged	4.003 → 0.029 s	1.003 [0.996, 1.025]	5.224 → 0.040 s	0.999 [0.987, 1.009]
gnlse_raman_shock	2.097 → 0.033 s	1.002 [0.926, 1.010]	2.801 → 0.048 s	1.004 [0.987, 1.005]
gnlse_soliton	0.453 → 0.022 s	0.999 [0.995, 1.005]	0.636 → 0.035 s	0.995 [0.990, 1.019]
rdw	1.405 → 0.051 s	0.992 [0.985, 1.003]	1.827 → 0.068 s	0.993 [0.990, 1.005]
argon_blueshift	1.321 → 0.035 s	0.995 [0.985, 0.999]	1.735 → 0.048 s	0.998 [0.996, 1.002]
argon_blueshift_envelope	0.473 → 0.024 s	0.989 [0.976, 1.021]	0.678 → 0.038 s	1.007 [0.995, 1.016]
radial_filament_256	0.613 → 0.093 s	0.994 [0.986, 1.022]	0.820 → 0.107 s	1.008 [1.004, 1.011]
ctc	5.033 → 0.056 s	0.985 [0.983, 1.016]	6.565 → 0.068 s	1.002 [1.001, 1.004]
kerr_beam_192	1.248 → 0.657 s	0.993 [0.972, 1.006]	1.079 → 0.255 s	0.995 [0.973, 1.126]
grin_nonlinear			1.984 → 0.309 s	1.015 [1.000, 1.040]

Warm windows run as fast as cold ones: the plans are those a measurement chose, in this process or an earlier one. A warm process given the seeding process’s wisdom reproduced its final field bit for bit (averaged helium), where two cold processes differ by measured plans (5.8 × 10⁻¹⁵ relative). What remains of a warm construction is allocation and setup: tens of milliseconds, a few tenths of a second for the free-space Kerr beam’s 1.4 GB of state. A job of a one-thread scan still pays solprop’s own import, case construction and translation (0.7 s and 3.0 s for averaged helium on a Xeon Gold 6240), which the native backend does not own.

Band column passes and the general route’s real signal

Host-only changes found by profiling ctc and grin_nonlinear (Kokkos label timer and perf, Intel Xeon Gold 6240, GCC 16.2 -march=cascadelake, SLEEF’s AVX2 functions). In a one-thread ctc window (0.01 of the guide) the band plasma column took 1.00 of 1.99 s inside Kokkos, 25 ns a sample, and the band statistics another 0.29 s; in a grin_nonlinear window (0.05) the complex inverse transform took 2.53 of 8.02 s, the projection 1.81 s and the synthesis 1.57 s.

• The band column’s passes as block loops. In the column’s own code a third of the samples went to the integrand pass (an out-of-line strict call per sample for the ionised fraction, a scalar division in the loss guard, a live-window note per sample), two fifths to the scans (Kokkos’ host team scan runs each closure twice, the second time from the same zero), and a fifth to the rates’ classification and scatter. A rate block of 256 samples that is mostly live now takes every sample through the vector-maths passes, a quiet one’s magnitude replaced by the table’s lower limit, and selects the limits after; a mostly quiet block gathers as before. The fractions come from one strict block loop, the guard’s branch is a selection that divides by one where it is not taken, the live window is searched from a block’s ends, and a one-thread team scans in one pass. Every value is unchanged to the bit. Label timer, column: 1.00 → 0.69 s at one thread, 0.147 → 0.100 s at sixteen.

• The filter one tap pair at a time, on x86-64. GCC had vectorised each sample’s nine-tap sum across the taps, a serial chain of additions per sample. A thread now takes blocks of 512 samples and runs the interior as one loop per tap pair; each sample keeps its order of terms, and only the contraction of its multiply-adds can move it (one ctc step: 2.2 × 10⁻¹⁶ of the field’s peak). Column: 0.69 → 0.54 s at one thread, 0.100 → 0.082 s at sixteen. Other hosts keep the per-sample traversal: the earlier tap-first filter slowed an Apple M-series argon window (see the rejected leads below).

• The rate table’s cubic vectorised out of line. Called from the blocks, the table’s cubic ran scalar, a division per sample: the compiled-out-of-line loop stored through the pointer it read while it gathered from the table, which GCC could not prove apart. It now writes a local block. Bitwise; column 0.54 → 0.51 s.

• The quiet-block skip. A rate block with no live sample takes only the table’s limits, without the gather, and the single gas’s density is read once per node, not after each store of the current (which the compiler had to assume might alias it). Bitwise; argon_blueshift column 0.72 → 0.64 s of a 0.05 window at one thread, ctc 0.51 → 0.49 s.

• The band statistics in the same blocks. The statistics scan took each sample’s rate and diagnostic fraction through the system libm and its whole-row prefix through Kokkos’ two-pass scan; it now uses the column’s rate blocks, block fractions and the one-pass scan. ctc: 22.8 → 14.0 ms a call at one thread. The statistics’ rates now carry the vector maths’ rounding: the peak ionisation rate moves by 3.6 × 10⁻¹⁴ relative (a one-ulp change of log |E| through the cubic’s slope), every other statistic not at all.

• The real form’s field through a complex-to-real transform. A real-form cubic reads only the real part of the analytic signal, which is the c2r transform of the one-sided spectrum with its doubled bins halved (exactly) and the imaginary parts of the zero and Nyquist bins dropped. The host now synthesises that half spectrum and transforms it c2r; the buffers are the complex transform’s, aliased, so the ledger and solprop’s memory estimate are unchanged. Not with SEP 6’s dispersive response, and not on devices. One grin_nonlinear step: 4.7 × 10⁻¹⁷ of the field’s peak. Inverse transform 2.53 → 1.21 s, respond 0.68 → 0.50 s at one thread.

• The projection’s block of bins up to the whole row. A host projection team walked the chunk’s rows over at most 64 bins, 1 KiB of each map row in turn; the block now widens to 2048 bins and to the first power of two covering the bins, so at one thread a team streams a mode’s map rows contiguously. No bin’s sum changes order. Project 1.81 → 1.17 s at one thread (2048 bins), 0.20 → 0.17 s at sixteen (512 bins).

• The synthesis summed in local blocks. The host synthesis zeroed each row and added every mode’s products into it in memory; it now zeroes only the bins outside the crop and sums blocks of 256 bins in two local arrays before one store, each bin’s sum in the same order. Synthesise 1.47 → 1.31 s at one thread. At one thread the synthesis and the projection each stream a 122 MB map per evaluation at 8–9 GB/s, about what one core of this socket draws from memory.

Each case measured with python -m benchmarks ab against the parent of these changes (fresh interleaved processes on 18 cores of one socket of a Xeon Gold 6240, passive bound threads, SOLPROP_FFTW_PLAN=measure, --length 0.05, three rounds of one repeat; medians of the per-round pairs, range over rounds). Every pair did the same work, and the final fields agree within the measured-plan floor (2 × 10⁻¹⁵ to 3 × 10⁻¹⁴):

case	one thread	sixteen threads
ctc	0.703 (0.694–0.705): 7.74 → 5.46 s	0.846 (0.829–0.868)
argon_blueshift	0.641 (0.627–0.649)	0.853 (0.844–0.868)
rdw	0.865 (0.859–0.886)	0.982 (0.949–0.984)
grin_nonlinear	0.701 (0.687–0.705): 8.04 → 5.67 s	0.802 (0.800–0.805)
argon_blueshift_averaged	0.913 (0.911–0.947)	0.919 (0.901–0.922)
he_rdw_averaged	0.945 (0.773–0.951)	0.958 (0.947–1.019)

The two averaged cases were measured without the quiet-block skip, which affects the resolved cases only (a separate run).

On 16 cores of an AMD EPYC 7763 (GCC 16.2 -march=znver3, the same protocol on a fragmented allocation; all but the quiet-block skip), one / sixteen threads:

case	one thread	sixteen threads
ctc	0.760 (0.754–0.770): 6.37 → 4.84 s	0.927 (0.924–0.929)
argon_blueshift	0.755 (0.744–0.765)	0.937 (0.904–0.940)
rdw	0.884 (0.874–0.885)	0.975 (0.963–1.020)
grin_nonlinear	0.750 (0.750–0.766): 5.53 → 4.16 s	0.872 (0.845–0.881)
argon_blueshift_averaged	0.918 (0.912–0.927)	0.922 (0.906–0.929)
he_rdw_averaged	0.985 (0.943–1.056)	0.967 (0.941–0.975)

With the quiet-block skip as well (another allocation, the same CPUs): argon_blueshift 0.720 (0.715–0.724) / 0.937, ctc 0.754 (0.751–0.755) / 0.907, rdw 0.883 / 0.989, grin_nonlinear 0.749 / 0.918 (0.835–0.918).

Not adopted: summing each block of eight scan increments locally before the running prefix, so that only block totals form a dependent chain (ctc column 0.538 → 0.534 s, within noise, for a change of association).

The native performance push: main at 740b24d against its result

The native performance push started from main at 740b24d, the merge of SEP 39’s native Adams march, and changed only the native backend and its build. It left the benchmark cases and stepper settings unchanged. The sections above record each adopted change with its own A/B. This section compares the two ends.

Complete runs (ab --length 1.0, three interleaved rounds, fresh processes) on a dual Intel Xeon Gold 6240, timed on the 18 physical cores of one socket. Every case took the same steps.

Case	One thread (field difference)	Sixteen threads (field difference)
gnlse_raman_shock	0.841 (\(3.6\times10^{-13}\))	0.765 (\(1.0\times10^{-12}\))
gnlse_soliton	0.802 (\(9.4\times10^{-14}\))	0.796 (\(7.9\times10^{-14}\))
rdw	0.659 (\(7.9\times10^{-14}\))	0.704 (\(1.7\times10^{-13}\))
he_rdw_averaged	0.608 (\(2.4\times10^{-10}\))	0.505 (\(2.1\times10^{-10}\))
argon_blueshift	0.563 (\(5.4\times10^{-14}\))	0.366 (\(1.1\times10^{-13}\))
argon_blueshift_averaged	0.605 (\(5.3\times10^{-14}\))	0.697 (\(1.3\times10^{-13}\))
argon_blueshift_envelope	0.724 (\(2.0\times10^{-12}\))	0.818 (\(1.7\times10^{-12}\))
h2_raman_averaged	0.645 (\(1.4\times10^{-13}\))	0.387 (\(9.5\times10^{-14}\))
radial_filament_256	0.623 (\(4.7\times10^{-15}\))	0.627 (\(4.5\times10^{-15}\))
grin_nonlinear	0.684 (\(1.7\times10^{-12}\))	0.792 (\(6.9\times10^{-13}\))
kerr_beam_384	0.559 (\(1.1\times10^{-15}\))	0.640 (\(1.1\times10^{-15}\))
kerr_beam_192	0.545 (\(9.4\times10^{-16}\))	0.659 (\(9.5\times10^{-16}\))
ctc	not run (hours)	0.679 (\(2.0\times10^{-10}\))
geometric mean	0.649	0.627 (without ctc)

The geometric means are 0.649 at one thread and 0.627 at sixteen. Short windows (--length 0.1) gave 0.663 and 0.613, but individual cases differ. rdw at sixteen threads is 0.337 in a window and 0.704 over a complete run, because the window’s baseline is dominated by per-run costs that the push removed (an idle NumPy BLAS pool spinning beside the native team).

Short windows on an exclusive AMD EPYC 7763 node (--length 0.1, three rounds, 128 cores, timing at one and sixteen threads). Every case took the same steps, with fields within \(2\times10^{-13}\).

Case	One thread	Sixteen threads
gnlse_raman_shock	0.839	0.784
gnlse_soliton	0.886	0.882
rdw	0.685	0.934
he_rdw_averaged	0.663	0.785
argon_blueshift	0.662	0.448
argon_blueshift_averaged	0.595	0.718
argon_blueshift_envelope	0.863	0.928
h2_raman_averaged	0.618	0.351
ctc	0.574	0.748
radial_filament_256	0.528	0.751
grin_nonlinear	0.763	0.811
kerr_beam_384	0.817	0.858
kerr_beam_192	0.824	0.875
geometric mean	0.707	0.735

Complete runs on an NVIDIA A40 (ab --device --length 1.0, three rounds). Every case took the same steps.

Case	main s	branch s	Ratio	Field difference
he_rdw_averaged	1.330	0.476	0.358	\(1.2\times10^{-10}\)
argon_blueshift_averaged	6.691	2.731	0.408	\(1.7\times10^{-13}\)
h2_raman_averaged	10.119	6.171	0.610	\(1.1\times10^{-13}\)
rdw	3.302	1.377	0.417	\(7.9\times10^{-14}\)
radial_filament_256	0.701	0.528	0.753	\(3.6\times10^{-15}\)
kerr_beam_384	6.455	5.036	0.791	\(1.7\times10^{-15}\)

On this card main runs some cases at one of two speeds from process to process. For rdw and radial_filament_256, the main column is therefore the faster of two complete-run jobs, and the ratio is the branch’s median over it, which makes it conservative. The final tree’s own job measured 0.403 and 0.392 directly. The other ratios are that job’s direct medians of paired rounds.

On the A40 the complete-run gains are larger than the short windows’. A window spends a larger share of its time in solprop’s per-run host work, which both arms share.

The other wheel platforms (ab --length 0.2, three interleaved rounds, each case at one thread and at the machine’s physical cores; the same steps everywhere, fields within \(2.1\times10^{-12}\)). The push was tuned on Linux x86-64; a first pass on these platforms found two regressions at full width, both on one-column problems: argon_blueshift 1.31 on macOS and gnlse_raman_shock 1.13 on Linux aarch64. Their causes and fixes are in a stable kernel team, compact CPU teams and one-mode products. The table compares main with the branch with those fixes: the CI wheels of both on macOS and Windows; on aarch64 main’s CI wheel against a local DISTRIBUTION build of the same commit.

Case	macOS arm64, Apple M5 Max, 1 / 18 threads	Linux aarch64, Google Cloud C4A (Neoverse V2), 1 / 16	Windows x86-64, Google Cloud C4D (AMD EPYC 9B45), 1 / 8
gnlse_raman_shock	0.969 / 0.769	0.982 / 0.704	0.773 / 0.941
gnlse_soliton	0.983 / 0.979	0.978 / 0.982	0.889 / 0.893
rdw	0.952 / 0.995	0.973 / 0.518	0.795 / 0.858
he_rdw_averaged	0.989 / 1.010	1.005 / 0.569	0.774 / 0.877
argon_blueshift	0.874 / 0.975	0.880 / 0.619	0.861 / 0.648
argon_blueshift_averaged	0.954 / 0.974	0.968 / 0.953	0.714 / 0.826
argon_blueshift_envelope	0.981 / 1.004	1.004 / 1.014	0.925 / 0.970
h2_raman_averaged	1.010 / 0.858	0.991 / 0.593	0.686 / 0.770
ctc	– / 0.754	– / 0.878	– / 0.717
radial_filament_256	0.715 / 0.846	0.722 / 0.586	0.564 / 0.764
grin_nonlinear	– / 0.872	– / 0.896	– / 0.764
kerr_beam_384	– / 0.954	– / 0.726	– / 0.827
kerr_beam_192	0.902 / –	0.820 / –	0.725 / –
geometric mean	0.929 / 0.911	0.928 / 0.733	0.764 / 0.816

The first-pass tables before the fixes, the bisection and the sizing records are kept with the project’s development records. The M5 Max ran with about 1.5 cores of background load (each round’s ranges stay within a few per cent); on an Intel Xeon Gold 6240 the fixes leave every case within the A/A spread (geometric mean 0.999 against the unfixed branch; main against the fixed branch 0.646 over both widths).

Accuracy. Every A/B took the same steps.

• On a Xeon Gold 6240, complete runs of he_rdw_averaged and ctc moved the final field by about \(2\times10^{-10}\), against \(10^{-13}\) in short windows. main against itself, over the same complete runs, gives \(6\times10^{-11}\) to \(2\times10^{-10}\): measured FFTW plans differ between processes at about \(10^{-16}\), and the adaptive controller amplifies the difference over a long run. The branch’s differences are at that floor.

• On the A40, fields are reproducible bitwise between processes. The device’s contraction of the Adams passes, its reduction order and the split transforms each change a step by rounding. Together they move he_rdw_averaged by \(10^{-11}\) in a window and by \(1.2\times10^{-10}\) over a complete run, with the same steps: the adaptive controller amplifies the rounding, as on the CPU.

• Every other case agrees to within about \(2\times10^{-12}\).

Measured and not adopted.

• Compiler options: -fno-math-errno -fno-trapping-math -fcx-fortran-rules, -O2, IPO off, 512-bit vectors, SLEEF’s AVX-512 functions on a Xeon Gold 6240, branch padding.

• Threading: spinning OpenMP waits, binding threads by default (on an exclusive EPYC allocation, OMP_PROC_BIND=close OMP_PLACES=cores is worth about 11 % on average; it was neutral on a Xeon Gold 6240), a BLAS width per product, serial or full-width medium operations.

• Sixteen-thread Raman and one-column team kernels: they gained on one CPU and lost on the other.

• Hand-written skinny products, Adams lane blocks of 16 and 64, band statistics through vector maths.

• On devices: a whole-transform node chunk, eight rows a thread in the compute capability 8.6 GEMM, other device reduction block sizes, band column tilings, and quadrature partials in chunks.

Outside the native backend. On complete runs, solprop’s host work outside the steps is a fixed cost per run, not per step: 0.1–0.3 s on the capillary and radial cases, about 6 s on kerr_beam_384. That is 2–9 % of a complete one-thread capillary run, and up to a third of the complete radial_filament_256 (sixteen threads) and kerr_beam_384 runs, which take one to twenty seconds. The configuration hash is computed once, at the first output write: 0.07–0.27 s, or 1.7 s for kerr_beam_384’s configuration. The cost that repeats is a saved field plane of a large Cartesian field. For kerr_beam_384’s 1.77 GB field that is about 3 s per plane:

• the saved plane’s copy, 1.1 s;

• the conversion to the optics convention, 1.0 s;

• the output store’s copy, 0.6 s;

• the native observation, 0.35 s.

A capillary or radial plane costs 7–25 ms.

The run record per plane

#427 bounds the controller’s norm history (a deque of the controller’s lag plus two, on the per-step Runge–Kutta and batch path) and makes the per-plane run record incremental: one RunOutcome and one validation of its record per plane, CallLedger.orders run-length encoded as it grows, and /meta/outcome rewritten in place. Before: main at 76eb2e57; after: the same with the pull request merged. One thread, numba_cpu, on an Apple M5 Max (18 cores, 128 GB, macOS 27.0.1, Python 3.14.5, NumPy 2.5.2, numba 0.67.0, h5py 3.16.0), a machine in ordinary use (load 4–10): a labelled A/B, to be read as ratios, not times.

python -m benchmarks ab --length 0.02 --rounds 4 --repeats 7 (the harness’s MemoryOutput, its 11 planes inside the window, so a plane every one to four steps):

case	stepper	steps (rejected)	before → after	ratio	A/A
gnlse_soliton	rk	12 (0)	14 → 14 ms	0.962 [0.906, 1.127]	1.009 [0.994, 1.057]
gnlse_soliton	multistep	18 (1)	16 → 15 ms	0.955 [0.887, 0.987]	0.991 [0.982, 1.034]
he_rdw_averaged	rk	12 (0)	46 → 46 ms	1.002 [0.992, 1.012]	0.990 [0.924, 1.036]
he_rdw_averaged	multistep	37 (1)	54 → 54 ms	0.990 [0.882, 1.026]	0.992 [0.986, 1.034]

Geometric mean 0.977 (A/A 0.996); the same steps and final fields in both builds. The harness writes no file, so the in-place /meta/outcome was timed by a hand script with the same alternation (fresh processes, four rounds of five repeats after a warm-up): the first centimetre of gnlse_soliton’s fibre (1536 bins) with 40 save planes, through run_problem into an HDF5Output:

stepper	steps	before → after	ratio	A/A	file
Multistep()	243	101.2 → 95.0 ms	0.947 [0.930, 0.968]	0.988 [0.837, 1.012]	1.50 → 1.35 MB
RungeKutta()	48	48.6 → 48.2 ms	0.987 [0.903, 1.001]	0.955 [0.714, 1.025]	1.21 → 1.13 MB

Nothing measured is slower. The Multistep file run, a plane every six steps, is 5 % faster, outside its A/A range; the others sit inside theirs. The record’s old cost grew with the steps times the planes, so these windows show only its start.

A run’s planes without a copy

The fix of #448 hands an output each plane without copying its field, by adopting the observation’s arrays, and adds one branch to StepDriver.observe for the held plane of #345. Before: f004a662; after: 643db910 on the pull request’s branch. An exclusive node, AMD EPYC 7763, 16 cores in one NUMA domain, --threading passive-bound, SOLPROP_FFTW_PLAN=measure.

python -m benchmarks ab --length 0.1 --rounds 5 --repeats 1, as the median of the per-round paired medians (after / before), with its range over the rounds:

case	threads	ratio	range
rdw	1	0.9911	0.980 to 0.998
rdw	16	0.9982	0.990 to 1.007
far_field_stats	1	1.0017	0.995 to 1.012
far_field_stats	16	0.9980	0.544 to 0.999

Geometric mean 0.9972, the fields within 3.6e-15 and the same work in both builds. The range of far_field_stats at 16 threads has one outlier round, in which the baseline took 16.2 s. This shows no regression. It does not measure the saving at large fields: these cases save few planes at a field size where the copy is small.

The saving is in memory, and is measured separately, on a loaded machine and so as field volumes rather than times. A delay scan through SimulationRunner with a handler that reduces each plane, against a hand loop over the backend protocol, on native CPU: a low-storage Kerr beam at 128 x 128 x 80 (20 MiB a volume), a window of 6 steps and 3 planes per input, the peak above the level before each input after a warm first one, on an Intel Xeon Gold 6240. The change is this fix together with #435 and #447, which release the input:

measure	hand loop	runner, before	runner, after
peak by tracemalloc, volumes	1.006	3.003	1.063
peak resident, volumes	0.94 to 0.99	2.86 to 2.94	0.93 to 1.00

numpy_ref, whose stages allocate per pass, went from 14.50 to 12.50 volumes by tracemalloc, the hand loop’s figure. Not measured at full size, and not extrapolated to it.

Single-precision collar and batches

#428 keeps the reference backend’s split collar (apply_temporal_rate) in binary32 under single precision, as the correction a + F(F^-1(a)(f - 1)); before, the product promoted to binary64 and was narrowed. Its single-precision batches now run through the plain LockstepRunner, each member its own binary32 trial; before, NumpyBatchBackend computed them in binary64 and narrowed the result, so the before column is the wrong arithmetic’s cost and the fields are not compared. ab selects neither a precision nor a batch, so both were timed by the hand script above, same machine and builds (main at 76eb2e57, then #428 merged), one thread, numpy_ref, on the same soliton.

The collar, 5 mm of RungeKutta() (43 accepted steps, one collar pass each), absorber_scheme="split_step", four rounds of five repeats:

bins	before → after	ratio	A/A
1536	34.0 → 34.0 ms	1.003 [0.999, 1.169]	0.993 [0.957, 1.003]
16384 (8 ps, 300 nm–4 µm)	183.1 → 180.2 ms	0.989 [0.971, 1.019]	0.992 [0.960, 1.020]

The pass alone, same process (best of seven times 200 calls): 18.3 → 15.9 µs at 1536 bins and 191 → 125 µs at 16384, 2–5 % of a step, which the runs do not resolve.

Batches, ten fixed steps over 2 cm with no collar, members 1 % apart in peak power, the stepping alone timed (propagate_to; the members’ construction is not):

bins	members	rounds	before → after	ratio	A/A
1536	4	4	29.8 → 30.1 ms	1.003 [0.963, 1.029]	1.025 [1.011, 1.059]
1536	16	4	116.5 → 122.3 ms	1.048 [1.036, 1.075]	1.028 [0.958, 1.072]
1536	16	8	136.5 → 141.5 ms	1.059 [0.984, 1.119]	0.981 [0.930, 1.030]
16384	8	8	496.6 → 389.7 ms	0.779 [0.676, 0.852]	0.971 [0.890, 1.017]

On a small grid with many members the member-by-member trials cost 3–6 %, at the edge of the A/A spread; on a 16384-bin grid they are 22 % faster, binary32 transforms replacing the binary64 ones. What vectorising over the members is worth was read from the binary64 batch, NumpyBatchBackend against LockstepRunner on the same members (unpaired medians of the same runs): 108.7 against 117.5 ms (16 members, 1536 bins) and 469 against 491 ms (8 members, 16384 bins), 4–8 %. A vectorised binary32 batch would recover at most that, so the reference keeps the plain runner.

A spectral walk-off absorber against the temporal collar

What the spectral walk-off absorber saves is the temporal pass, which a zero temporal rate compiles away; this entry says what that is worth, in place of a figure measured on a backend since retired. A hand script, on the one grid, one backend and one build: a field grid over 40 µm of fused silica, 143 nm to 2 µm, trange 250 fs at the default collar_fraction, so a 326 fs window of 2048 samples; envelope Kerr, an 800 nm 30 fs launch; ten fixed RungeKutta() steps of 1 µm, each case translated once beforehand and the runs alone timed. The cases are the default collar under each scheme, and custom_absorber(grid, rate, collar=absorber(grid, temporal=False)). Native CPU, eight threads bound to eight cores of a dual Intel Xeon Gold 6240 (18 cores a socket), built for -march=cascadelake. The machine was shared and in use (load 16–58 of 72), so these are to be read as ratios of medians, and each is given over its separate sessions, cycling the three cases within a session:

grid	sessions (repeats)	spectral / split_step	spectral / in-stage
line, 256 columns	2 (5, 9)	0.92, 0.90	0.58, 0.53
Cartesian, 64 × 64	2 (5, 5)	0.91, 0.89	0.47, 0.53

Single runs carried outliers of up to three times the median, which the median sets aside. The ratio against the split repeats to within 3 %; the one against the in-stage term less well, between 0.47 and 0.58. Over these ten steps the three runs’ fields agree to 1e-15 of the on-axis column’s norm inside the request, because no light has reached a collar: the entry times the pass and compares no absorption. It is a window, and the temporal pass’s share of a step moves with the response’s: a heavier response gives a smaller saving.

Reference tables

These tables back the indicative results above.

Kokkos CPU sizing thresholds

Counts are elements or scalar product work, never allocated bytes. Stages and reductions use the field’s element count; batch updates use the fused batch’s total. Packing uses matrix entries and temporal copies use output samples, without multiplying by the dense transform’s inner dimension. Responses use active columns × components × fine-grid samples. Dense products use output rows × inner dimension × samples; coupled matrix functions instead use active spectral bins × modes³, and the exponential Adams trial’s passes the field’s elements × history rows. FFT and Raman choices use each plan’s actual real-side batch dimensions, including a kernel transform and the last partial chunk.

Operation	Serial below	Full pool from	Intermediate rule
Stage/update/reduction	8,192	8,192	none
Packing/temporal copy	65,536	65,536	none
Response	32,768	196,608	up to four threads; at least 40 component columns or 65,536 fine samples per column also enables the full pool
Native/custom GEMM	1,048,576	1,048,576	none
General-modal synthesis/projection/maps	65,536	65,536	none
Coupled matrix functions	32,768	32,768	none
Exponential Adams trial passes	8,192	65,536	up to four threads
Modal/spatial FFT	32,768	131,072	up to four threads
Raman FFT	65,536	262,144	up to four threads

A response column is one active node and one component. Every width is capped by the global thread setting. A kernel pool that already holds a full-width instance gives it to the medium-width kernels too (a stable kernel team); the FFT pools keep their medium widths. A problem of one column on one node takes at most eight threads in every pool, its transforms on its kernels’ team (compact CPU teams). Explicit per-pool requests bypass this table, including the short-FFT rule. build_info() exposes the constants and override mask; backend metadata exposes effective widths and sizing inputs. Non-OpenMP host builds report that per-operation sizing is unsupported and keep their existing execution policy.

Resident Adams allocations

F = N² × 600 × 16 bytes is one complex128 field of the reduced kerr_beam_384 above and C = 5592 × 600 × 16 bytes one node chunk (a 1 GiB node budget). Each pair is numba_cpu / Kokkos capacity in GiB; arrows show the 29 → 18 volume change. The 384² capacities are calculated. These are owned allocations, not resident pages.

Allocation	Field count, numba_cpu / Kokkos	128²	256²	384²
Adams history ring	7 / 7	1.025 / 1.025	4.102 / 4.102	9.229 / 9.229
Adams moments	8 / 7	1.172 / 1.025	4.688 / 4.102	10.547 / 9.229
Adams flow / propagator	1 / 1	0.146 / 0.146	0.586 / 0.586	1.318 / 1.318
Native Adams source, exponent, saved state	0 / 3	0 / 0.439	0 / 1.758	0 / 3.955
Dense slopes	0 / 7 → 0	0 / 1.025 → 0	0 / 4.102 → 0	0 / 9.229 → 0
Dense start, end, exponent, held observation	0 / 4 → 0	0 / 0.586 → 0	0 / 2.344 → 0	0 / 5.273 → 0
Current state	1 / 1	0.146 / 0.146	0.586 / 0.586	1.318 / 1.318
RK stages	7 / 7	1.025 / 1.025	4.102 / 4.102	9.229 / 9.229
Solution, error, advanced, recycled stage	4 / 4	0.586 / 0.586	2.344 / 2.344	5.273 / 5.273
RK exponent / weighted residual	1 / 1	0.146 / 0.146	0.586 / 0.586	1.318 / 1.318
RK propagators	2 / 1	0.293 / 0.146	1.172 / 0.586	2.637 / 1.318
Transverse FFT scratch	1 / 1	0.146 / 0.146	0.586 / 0.586	1.318 / 1.318
Field-volume subtotal	32 / 44 → 33	4.688 / 6.445 → 4.834	18.750 / 25.781 → 19.336	42.188 / 58.008 → 43.506
Temporal FFT buffers (4C each)	fixed chunk	0.200 / 0.200	0.200 / 0.200	0.200 / 0.200
Absorber scratch (0 / 2C)	fixed chunk	0 / 0.100	0 / 0.100	0 / 0.100
Statistics spectrum, time, power (0 / 2.5C)	fixed chunk	0 / 0.125	0 / 0.125	0 / 0.125

The remaining ledger entries (reduction partials, axes, spatial coefficients, statistics outputs, host problem tables, masks, weights and plan heap) total a few MiB. Outside the ledgers, the input field adds one F; Kokkos keeps another F of host state, and field 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 an Apple M5 Max and 3.14.7 on a Xeon Gold 6240; 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 M5 Max’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.

The run record’s cost

#309 measured how a run’s record grows with its grid. Moving the record’s large arrays beside its text (RECORD_ARRAY_LIMIT) takes the hash from the record’s text to its arrays’ bytes.

Measured on an Intel Xeon Gold 6240, pinned to one socket with 16 threads, on a shared machine with a load average of about 17, so treat small differences as noise. One run each, built and not stepped, on numpy_ref, with #309’s Cartesian problem of 96 frequency bins. The figures are for the record alone, not for a run.

Grid		First hash, s	to_json, s	HDF5Output.open, s	Reading the file’s record, s	Hashing a read record, s	JSON, MB	File, MB
256²	before	0.442	0.163	0.171	0.177	0.174	5.61	5.62
	after	0.010	0.001	0.006	0.021	0.001	0.03	3.19
512²	before	1.763	0.652	0.711	0.755	0.777	22.66	22.67
	after	0.032	0.001	0.011	0.065	0.001	0.05	12.64
768²	before	4.250	1.588	1.711	2.155	1.735	55.59	55.60
	after	0.075	0.002	0.017	0.144	0.002	0.07	28.40

A default radial run of 256 nodes and 8192 frequency bins, which carries input noise, holds 83.21 MB of JSON before and 0.04 MB after, and its first hash goes from 2.603 s to 0.077 s. A 64-node grid is about the same before and after, since its arrays sit close to the limit. A run into MemoryOutput computes no hash at all.