meta(cpu): CPU apply! path is single-threaded — Julia threads give zero speedup

[kylebeggs]

Context

The new Chmy.jl comparison benchmark (comparisons/chmy_diffusion_2d/, Apple M5 Pro, Julia 1.12.6) shows MFO's CPU operator application is completely insensitive to -t — the thread pool sits idle:

interior	MFO apply! -t 1	MFO apply! -t auto (5 threads)
126²	9.5 µs/step	9.6 µs/step
510²	148.5 µs	144.3 µs
1022²	569.4 µs	569.2 µs
2046²	2653.9 µs	2625.3 µs

Chmy's KernelAbstractions CPU backend does scale with threads, and overtakes us at DRAM-resident sizes despite moving ~40% more memory per step (7 array passes for its materialized flux form vs our fused 5): at 2046² with 5 threads, Chmy runs 2193 µs/step vs our 2654 µs. Single-threaded we win every size 2.5–16.9×, so the loss is purely about parallelization, not kernel design.

Root cause

Every leaf apply! body is a plain serial Base broadcast over a SubArray (src/operators/laplacian.jl, gradient.jl, divergence.jl, scaling.jl). The KA kernels that do exist (src/operators/forest_packed.jl, src/transfer_kernels.jl) only launch when backend isa KernelAbstractions.GPU; on CPU they fall back to the serial host loop. There is no Threads.@threads anywhere in src/.

Why it costs so much on modern hardware

One core cannot reach the package's aggregate memory bandwidth: STREAM triad on this M5 Pro measures 129 GB/s single-thread vs 246 GB/s with 5 threads. A perfect serial stencil sweep therefore leaves ≥2× bandwidth unused — this gap only grows on higher-core-count machines, and it is exactly where Chmy beats us.

Options

1. Launch leaf kernels through the KA CPU backend. The grid already carries device::KernelAbstractions.CPU, and DESIGN.md names KA @kernel as the per-operator escape hatch — this is Chmy's model, keeps device-agnosticism, and gives threading for free. Cost: leaf bodies stop being array-level broadcasts, which is a structural decision DESIGN.md must record (AD-friendliness of the leaves is a locked invariant).
2. A threaded executor at the apply! seam. Keep leaf bodies as broadcast expressions but materialize them through a chunked Threads.@threads sweep over interior rows. Preserves the array-level-leaf invariant (Enzyme rules, Reactant tracing, Adapt), smaller blast radius, CPU-only.
3. User-level threading only (what examples/niederer_benchmark.jl does over forest leaves) — doesn't help the single-grid case, which is where we lose today.

Whichever route: the Enzyme rule constraints still apply (no allocation in rule bodies, no mixed GC-pointer/inline-float rule arguments), so the executor must not disturb the _exchange_storage! / _bc_storage! seams.

Acceptance

• If option 1 (leaf kernels through the KA CPU backend) is chosen, apply! must stay safe when a Field wraps a SubArray view — KA's CPU() backend is known to race on views (the downstream case in fix(cpu): apply! is unsafe when input/output are SubArray views #64, closed as not-an-MFO-bug because today's CPU path is plain broadcast). A test with view-backed input and output fields gates this.
• The forest reference path is tracked separately in perf(forest): thread the reference BlockForest path — leaf sweeps and exchange descriptor loops #83 (sub-issue of this one); pick the executor once for both.

comparisons/chmy_diffusion_2d/run.sh shows apply! scaling with -t, and MFO ≥ Chmy at every sweep size at -t auto. With 5/7 of Chmy's per-step traffic there is no structural reason to lose any threaded row.

---

_{🤖 Beep boop — Claude filed this one while the benchmark bars were still crawling across Kyle's terminal.}

[kylebeggs]

Extended the sweep to 4094², 8190², 16382² — two refinements to the picture above:

• Serial, the ratio stabilizes at ~2.6× in MFO's favor (MFO ~1.7 Gcells/s vs Chmy ~0.63 across all DRAM-resident sizes) — that's the clean measure of the fused-stencil advantage.
• Threaded, Chmy's edge is a bounded middle band, peaking at 1.5× ahead at 4094² and vanishing by 16382² (182.6 ms vs our serial 178.9 ms): its threaded throughput sags from 1.92 to 1.47 Gcells/s at the largest size as the 7-pass flux-form traffic (~15 GB/step) saturates the memory system. Their extra traffic caps what threading can buy them; our 5-pass loop threaded should sit above their curve everywhere, which strengthens the acceptance criterion here.

Full tables in comparisons/chmy_diffusion_2d/README.md.

---

_{🤖 Beep boop — Claude, back with bigger grids and the same conclusion wearing more data.}