space.irbem benchmarks

The generated benchmark tables for the space.irbem module: every routine that exists, on every lane it has, against the vendored IRBEM library run as a black box. The module page describes what the routines compute; this page holds the numbers.

Three things are being measured, and they answer three different questions:

1. Per routine — how fast is each routine on the host, on the device, and in IRBEM? The last column of that table is the one to read first: it says whether a routine can reach the GPU at all, and it is read out of gpu/dispatch.hpp — the registry and its launchers — rather than asserted, so the porting gap is visible instead of implied.

2. The crossover — at what batch size does the device start to pay? This is the module's central argument and the reason the batch entry points exist at all.

3. Arithmetic intensity — why the answer to (2) differs so wildly between kernels. Same seam, same transfers, same device; only the flops per byte changed, and the verdict flips.

The ground rules

The IRBEM lane quotes the -O2 rebuild, never the shipped binary. IRBEM's own makefile passes no -O at all, and docs/ERROR_BUDGET.md measures the shipped library at 2.7× the -O2 one. Benchmarking against the shipped artefact would inflate every ratio in these tables by that factor, and the resulting claims would be about GNU Fortran's default optimization level rather than about this code. The oracle is dlopened and called through the documented C bindings — run as a black box, never read, never linked.

Every rate is real time, never CPU time. A device lane spends nearly all of its wall clock blocked rather than running, so Google Benchmark's default CPU-time accounting would credit the GPU with several times its actual throughput. Every batch benchmark is marked UseRealTime, so the host and device columns of a ratio are the same quantity.

Transfers are inside the timed region. A kernel time that excludes the bus is a number about the ALU, not about an operation a caller can perform. The device columns include staging the inputs, the dispatch, and reading the results back.

Nothing here is hand-typed. scripts/bench_run.sh runs the binary and regenerates the region between each suite's BENCH markers from that run's JSON. The "on GPU?" column and the flops/byte column are printed by the benchmark binary itself, from the kernel registry and from a measured step count — so they cannot drift away from the code they describe.

A ratio is only quoted once its device lane has been checked against its host lane. For a long time this page did not do that, and could not: a benchmark that calls launch_igrf and gets true back has learned that a dispatch was submitted, not that the kernel unpacked the coefficient buffer the way the host packed it — and a kernel that reads the normalisation table at the wrong stride is exactly as fast as one that reads it correctly. irbem_bench --verify now runs every device lane against its own host lane and reports the largest relative deviation against docs/ERROR_BUDGET.md §5, and bench_run.sh check runs it before it grades a single row. Measured here, over the same input sets the tables are measured on:

Every lane is inside budget on all three backends, and the NVIDIA figures repeat digit for digit between runs. Two of them are not comfortably inside it: the IGRF magnitude sits at 92% of the Blocal budget and llvmpipe's I at 71% of the XJ budget. That is what an fp32 kernel measured against an fp64 reference looks like, and it is why those budgets are where they are rather than an order tighter — but it also means a kernel change that costs half an ulp will be caught by this gate rather than by a user.

Never -ffast-math. The suite builds -O3 -march=native -ffp-contract=off, which is the arithmetic the correctness tests pinned. -ffast-math would let the compiler reassociate the spherical-harmonic sums and the RK4 accumulation, and the binary would then be measuring different mathematics from the one that was verified. bench_run.sh check fails the build if the flag ever appears.

Every routine, on every lane it has

Read the on GPU? column, not just the timings. Three of its four answers are failure states, and the interesting one is registered, NO LAUNCHER: a kernel can be written, compiled, correct and even measurable in isolation while remaining unreachable through the seam every caller goes through. That is a porting gap, not a success.

The column is read out of gpu/dispatch.hpp — the registry for the binding count, and the header's own launch_* functions for whether anything actually reaches the kernel. It used to be read out of the registry plus a rule written into the benchmark ("four bindings means `dispatch_batch`, `irbem_trace_i_f32` means `launch_trace`, everything else is unreachable"), which is not a derivation at all: it is a second list, in a second file, and it went stale the way second lists do. launch_igrf landed in the seam and this column went on publishing registered, NO LAUNCHER about a kernel the seam had learned to launch. It now goes and looks.

Routines with an epoch-dependent geometry appear twice. (hot) reuses one Rotations across the ephemeris — the shape this module's API is built for — while (cold) rebuilds it per point, which is the like-for-like against IRBEM, whose entry points take a date with every call and cannot amortize. Reporting only the hot number would flatter this module; reporting only the cold one would hide its actual design advantage.

Measured by bench/irbem_bench.cpp; reproduce with scripts/bench_run.sh publish irbem-per-routine.

Also measured — variants, and routines this module has not ported

Four rows that need reading carefully

The trace comparison is against MAKE_LSTAR_SHELL_SPLITTING, and that is not a detail. IRBEM's MAKE_LSTAR is the obvious counterpart and is the one this page quoted until it was checked. It is the wrong one. IRBEM's own documentation says MAKE_LSTAR "computes the L\* parameter for locally mirroring particles (local pitch angle of 90 degrees)" and points the reader at MAKE_LSTAR_SHELL_SPLITTING "to compute L\* for arbitrary pitch angles" (docs/source/api/magnetic_coordinates.rst). A locally mirroring particle mirrors at the spacecraft, so its I spans the short arc between the point and its conjugate. This suite's lines mirror at B_local / sin²α for α = 30…80° — up to four times the local field, and a far longer arc to integrate. Called at the same eight start points, the two entry points return XJ = 1.85 Re and 5.5 × 10⁻³ Re at L = 2: two different integrals, not two implementations of one. They differ by 35× in time, and quoting the 90° one turned this module's trace win into a 0.27× loss. The 90° entry point is still measured — in Also measured, labelled for what it is — and at α = 90° the shell-splitting entry point reproduces its Lm and XJ to every printed digit while costing about twice as much, so roughly a 2× constant of its own sits inside the row this page now quotes.

What survives that correction is the smaller caveat: the two sides still trace at their own default step sizes. This module's TraceOptions::steps_per_l = 50, i.e. ds = L/50, against whatever make_lstar_shell_splitting1 chooses for itself. So the trace rows are a statement about each library's defaults as much as about its inner loop, and matching the two step counts is still owed before either is quoted as an inner-loop comparison.

L* drift shell has no CPU cell at all. driftshell.hpp::make_lstar landed while this suite was being written, so the row no longer says NOT PORTED — but nothing here times it yet, and an empty CPU cell beside IRBEM's cost is the honest rendering of "ported, not yet measured". A CPU and a device lane for it are owed.

Two device cells carry a wide spread. The dipole and IGRF device lanes finish in microseconds, so they sit close to the dispatch floor and contend with whatever else holds the queue. Their medians are directionally right — the dipole loses, the IGRF wins — but neither is a number to build a budget on. Across four independent runs of this suite on this machine the dipole device lane measured 3.50, 4.18, 6.42 and 14.01 ns/point and the IGRF device lane 7.09, 7.65, 12.06 and 48.67 (the last at ±116%), and the published speedup for the IGRF batch has read anywhere between 6.35× and 43×. Neither verdict has ever changed sign. Read the verdict and the spread; do not read the median as a budget. This is also why bench_run.sh check refuses to grade any row whose own repetitions disagree by more than 15% — the gate threshold would be inside that spread, and a gate that fires on the driver's mood is a coin flip.

Storm conditions, not only quiet ones

A radiation-belt library exists for the disturbed magnetosphere, so a benchmark suite that only ever measures maginput = 0 is measuring the one condition it is not for. Until ext_t89.hpp landed there was nothing here parameterized by activity to sweep; there is now, and the T89 rows are swept across the whole published Kp envelope — Kp × 10 from 0 to 90, IRBEM's maginput slot-1 scaling (docs/source/api/general_information.rst) — on both sides, at a real ~20° dipole tilt rather than the ψ = 0 case where Tsyganenko (1989)'s sin ψ terms drop out and the model evaluates a strictly cheaper expression than it ever does in flight.

The quiet and extreme-storm ends are in the table above; the eight bins between them are in Also measured. The sweep is flat in time on both sides, which is the expected answer and is now a measured one: T89 selects one of seven coefficient sets by Kp and then evaluates the same closed form, so activity changes the numbers and not the work. A model whose cost does move with activity — a T96 or a T01 root-find — would show it in exactly this sweep, and that is what these rows are for once one lands. Dst, Pdyn and southward Bz are still not swept anywhere, because no model in this module reads them yet; that sweep is owed with the first model that does.

The crossover — where the device starts winning

This is the whole argument for the module's shape. A single field-line trace is a serial RK4 chain, so no hardware makes it faster; the parallelism is entirely across lines, and it only becomes available if the caller hands over the whole batch at once. A loop calling trace_invariant per point cannot be accelerated no matter what hardware is present, which is why trace_invariant_batch exists and why it is the entry point to call.

Below the crossover the synchronous dispatch floor dominates and the host wins — which is why prefer_gpu consults the per-kernel crossover instead of assuming the device is always right. A blanket "always use the GPU" would make small batches slower and the dipole kernel slower at every size.

Every point of this curve traces the same distribution of field lines. The input set is indexed by the base-2 radical inverse rather than swept linearly, so the first 64 lines already cover L = 2…8 evenly and so do the first 65 536; the only thing that changes along the curve is the batch size.

Measured by bench/irbem_bench.cpp; reproduce with scripts/bench_run.sh publish irbem-crossover.

Arithmetic intensity — the number that predicts the verdict

Same seam. Same transfers. Same device. The three kernels differ in one thing — how much arithmetic they do per byte moved — and that one thing decides the verdict, including the sign of it.

The flop counts are a term model, stated rather than hidden: one spherical-harmonic coefficient slot is charged five flops (G = g·cos mφ + h·sin mφ, the multiply by the Legendre value, and the accumulation), which is the conventional way a harmonic model's cost is quoted. A full instruction count of the Slang inner loop — which also carries the Legendre recursion and the two horizontal components — comes out roughly six times higher. The argument these rows make is about the ratio between them, and a ratio is unchanged by a constant factor applied to every row.

The trace row's flop count is the only one that is measured rather than counted, because it depends on how many RK4 steps a line actually takes. It is also a lower bound: FieldLine::steps reports the two mirror-half integrations, and the initial walk to the magnetic equator is on top of that.

That is why the trace row here reads lower than the ~9 400 flops/byte quoted in gpu/dispatch.hpp's registry comment. Both are the same term model; they differ in the step count they assume, and this one's is measured over this input set — field lines started on the dipole equator of L = 2…8 at 30…80° pitch, where the walk to the true minimum-B point is short. A set of lines started off-equator would trace further and this number would rise. The figure is a property of the workload as much as of the kernel, which is exactly why it is measured here rather than written down once.

Byte counts are exact, and are the compulsory traffic per point. The IGRF coefficient and Legendre tables are ~1.7 KB for the whole batch — 0.03 bytes/point at 65 536 — so they do not appear.

All three rows are measured at 65 536 points (65 536 field lines for the tracer), which is past every one of their crossovers, so each row is the verdict at a batch size where the device has whatever advantage it is going to get. The two small kernels' device lanes are the noisy ones: they are over in microseconds and contend with the driver, so read their spreads and not only their medians. The IGRF row in particular has published a speedup between 6× and 43× across runs while never once losing; the row's claim is the sign and the order of magnitude, not the digits.

Measured by bench/irbem_bench.cpp; reproduce with scripts/bench_run.sh publish irbem-intensity.

mean RK4 steps per traced line, measured over 256 lines of the input set: 63.1

Reproducing, and the regression gate

scripts/bench_run.sh build                    # Google Benchmark at the pin, then the suite
scripts/bench_run.sh verify                   # every device lane against its own host lane
scripts/bench_run.sh publish all              # run everything and regenerate all three tables
scripts/bench_run.sh check                    # verify, then gate against bench/baseline.csv

bench/baseline.csv is a committed row per measured lane, in nanoseconds per point. bench_run.sh check re-measures and fails on a regression that clears both a 1.40× ratio and a 2 ns absolute floor — a routine measured in single nanoseconds moves 30% on machine load alone, and gating on ratio alone would make the gate a coin flip. A row that trips the gate is re-measured once before it is called a regression. The gate also fails if -ffast-math ever appears in the benchmark build, if bench/CMakeLists.txt and scripts/bench_run.sh stop agreeing about the flags, or if any device lane leaves its error budget — that check runs first, because a throughput gate cannot tell a fast kernel from a fast wrong one.

All four failure modes have been exercised deliberately and then reverted: a doctored baseline.csv row (fails, 21.01×), a deleted -march=native (fails, "the two build paths have drifted"), an added -ffast-math (fails, "it is banned"), and a device IGRF launch truncated from degree 13 to 12 (fails, max_rel_dev 1.02 against a 1e-06 budget). A gate nobody has watched fail is a gate nobody knows works.

Every table is stamped with the CPU governor it was measured under. These numbers were taken with scaling enabled and the governor in powersave, which is the honest default state of the machine rather than a tuned one; a pinned-frequency run moves the host columns and leaves the ratios roughly where they are, but the stamp is there so a future run that disagrees can be told apart from a regression.

The device is selected with CHEATAH_GPU_LINALG_VK_DEVICE (NVIDIA, Intel, llvmpipe). On a machine with no device, or a build with no cheatah-gpu-linalg on the include path, the GPU rows report themselves skipped and the table shows dashes: a truthful table, not a broken build. The same is true of the IRBEM columns when SPACE_IRBEM_ORACLE points at nothing.