Skip to content

milnor_gpu: fix >2^32-element master via cubecl 64-bit addressing - #17

Draft
JoeyBF wants to merge 102 commits into
rref-hpcfrom
milnor-gpu-u64-addressing
Draft

milnor_gpu: fix >2^32-element master via cubecl 64-bit addressing#17
JoeyBF wants to merge 102 commits into
rref-hpcfrom
milnor-gpu-u64-addressing

Conversation

@JoeyBF

@JoeyBF JoeyBF commented Jul 27, 2026

Copy link
Copy Markdown
Owner

Problem

The batched Milnor multiply silently corrupted results once the shared masks admissible-master crossed 2³² u16 elements (~8.6 GB, around S_2 stem ~160–180). The buffer length and per-R offsets were u32, so cubecl's default U32 address_type truncated them (len as u32) → out-of-range reads → dx != 0 (d²≠0) panics, e.g. deterministically at (180,92). (160,82) stayed just under the ceiling and was clean.

Diagnosis (instrumented): masks reached 4,452,195,742 u16 (8.29 GiB) > RESIDENT_MAX_CAP = 2³²−1; the ArrayArg length metadata also truncated (u32). Confirmed not the row-reduce, not the resident basis, not a cross-stream race — pure u32 addressing overflow.

Fix — idiomatic cubecl 0.10 64-bit addressing (no hand-rolled buffer splitting)

cubecl has AddressType {U32, U64}; a plain #[cube(launch)] defaults to U32. Two CUDA-backend quirks shaped the fix:

  1. address_type = "dynamic" miscompiles for us — it picks U32 for small blocks and then narrows the Array<u64> offset arrays on read (usize::cast_from(u64) under u32 usize), corrupting results. So the multiply is static "u64".
  2. checked-mode + U64 won't compile — the bounds clamp emits min(u64, u64), ambiguous for NVRTC. So it uses launch_unchecked (safe: every access is in-bounds by construction — uploaded need_* prefix + per-column j guards).

Changes (all in the batched multiply path):

  • multiply_batch_kernel#[cube(launch_unchecked, address_type = "u64")].
  • RInfo.cs_off/mk_off u32→u64; r_cs_offset/r_mk_offset bound as Array<u64>.
  • resident grow copies (copy_into_u16/_u32) → #[cube(launch_unchecked, address_type = "dynamic")].
  • drop RESIDENT_MAX_CAP + clamps; assert out_len ≤ u32::MAX loudly.

Validation (H200)

  • (160,82): rc=0 clean, 517 s vs 486 s baseline (+6.4% for static-u64 multiply, partly offset by unchecked dropping the clamp — no regression).
  • (180,92): ran 1500 s with zero dx panics (was a deterministic panic ~130–180 s), i.e. correct well past the 2³² masks boundary. A to-completion rc=0 run is in progress.

🤖 Generated with Claude Code

JoeyBF and others added 30 commits July 26, 2026 20:03
The batched Milnor multiply silently corrupted results once the shared
`masks` admissible-master crossed 2^32 u16 elements (~8.6 GB, around
S_2 stem ~160-180): the buffer length and the per-R offsets were u32, so
cubecl's default U32 `address_type` truncated them (`len as u32`), giving
out-of-range reads and `dx != 0` (d^2 != 0) panics at e.g. (180,92).
(160,82) stayed just under the ceiling and was clean.

Fix uses cubecl 0.10's first-class 64-bit addressing rather than
hand-rolled buffer splitting:

- `multiply_batch_kernel` -> `#[cube(launch_unchecked, address_type = "u64")]`.
  Static u64 (not "dynamic"): dynamic picks u32 `usize` for small blocks and
  then narrows the u64 offset arrays on read (`usize::cast_from(u64)` under a
  u32 address type), corrupting results. `launch_unchecked` because checked
  mode emits `min(u64, u64)`, which NVRTC rejects as an ambiguous overload;
  every access is in-bounds by construction (uploaded `need_*` prefix +
  per-column `j` guards).
- `RInfo.cs_off/mk_off` u32 -> u64; `r_cs_offset`/`r_mk_offset` bound as
  `Array<u64>`.
- Resident grow copies (`copy_into_u16`/`_u32`) ->
  `#[cube(launch_unchecked, address_type = "dynamic")]` (scalar usize offsets
  adapt without narrowing; dynamic keeps small copies on u32).
- Drop `RESIDENT_MAX_CAP` and its clamps (the u32 ceiling is gone); assert
  `out_len <= u32::MAX` loudly (the row-block splitter guarantees it).

Validated on H200: (160,82) rc=0 clean, 517s vs 486s baseline (+6.4% for the
static-u64 multiply, partly offset by unchecked dropping the bounds clamp);
(180,92) ran 1500s with zero dx panics (previously a deterministic panic
~130-180s), i.e. correct well past the 2^32 masks boundary.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…he resident host master

`resolve_through_stem S_2 "" 180 92` was OOM-killed at the 500 GB cgroup limit
(~280 GB anon + ~240 GB shmem). Measurement (env `NASSAU_MEM_REPORT`, plus a
pure-CPU control at ~5-7 GB) showed the resolution's own data — differentials
and module tables — is only ~1.5 GB; the ~500 GB was entirely GPU-runtime host
memory. Two fixes cut it to ~90 GB at stem 180 (dx-clean):

1. Default `NASSAU_GPU_STREAMS` 8 -> 1. cubecl gives each CUDA stream its own
   page-locked (pinned) host pool that `memory_cleanup` never trims; ≥2 streams
   ballooned pinned host memory to 140-240 GB, single-stream holds it at ~4-6 GB.
   The payoff of multi-stream is ~nil here: cubecl's server is single-threaded,
   so extra streams buy no CPU concurrency, only GPU kernel overlap the big
   saturating multiplies barely use. Override for a dedicated large-RAM node.

2. Stop retaining the resident admissible master host-side. `RESIDENT_HOST` kept
   the full `col_sums`/`masks` (~54 GB at stem 180) forever, duplicating the
   device copy, even though after upload it is never read again (offsets come
   from `index`; growth uploads only the new tail; a capacity realloc copies the
   old *device* buffer). Now it keeps only the not-yet-uploaded tail (`*_pending`,
   ~sub-GB) plus a logical length, freeing each `R`'s data the moment it reaches
   the GPU — invariant `dev.uploaded == len - pending.len()`.

Also: bound the pinned staging on resident growth to `STAGE_CHUNK` chunks, and a
gated `NASSAU_MEM_REPORT` (differentials/modules/resident heap breakdown) for
ongoing memory work.

Validated on H200: stem 180 peak RSS ~90 GB (was 500 GB OOM), dx-clean through
the >2^32 masks region and the heavy solves. Remaining growth for higher stems
is real GPU working memory (concurrent dense output matrices) + GPU device
memory, not retained host duplicates.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…unded

Even with the resident host-master freed, host RSS still grew ~linearly toward
the cgroup limit at high stems — not a retained duplicate but the GPU path
materializing whole dense output matrices. The CPU walks signatures one at a
time (small working set); the GPU offload built the entire matrix at once and
held its dense readback (`num_rows × num_limbs` u32, ~12 GB regions at stem 180)
alongside the assembled matrix. Two caps bound it:

1. `reuse_full_matrix` now only builds the all-rows matrix when
   `rows × cols <= NASSAU_GPU_REUSE_MAX_WORK` (default 1e10). Above that it falls
   back to per-signature builds (each a bounded row subset, like the CPU),
   so the peak scales with the largest single signature, not the whole bidegree.

2. `get_partial_matrix_restricted` processes rows in batches sized so the dense
   readback stays under `NASSAU_GPU_MAX_READBACK_MB` (default 1024). Products are
   built in row order, so each batch is a contiguous slice with its `row` remapped
   batch-local; results XOR back into the global rows. The readback is freed
   between batches.

Both are correctness-neutral (rows are independent): verified bit-for-bit vs the
CPU under `NASSAU_GPU_VERIFY` with a 1 MB cap (many batches/build), 0 mismatches.

Effect at stem 180 (streams=1 + resident de-dup + this): host anon goes from
growing past 85 GB to a bounded ~18-24 GB, with no throughput loss — cross-
bidegree rayon concurrency keeps the GPU fed, and the default cap only splits the
few giant high-t builds into a handful of chunks. Peak RSS ~30 GB (was 500 GB
OOM). This bounds the host working set so it scales toward higher stems; a D-deep
async launch/collect pipeline (cubecl launches are async; only read_one blocks)
can hide chunk latency if a much smaller cap is ever needed.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
The resident admissible master (col_sums/masks) grows with degree and is the
stem-300 device-memory wall: ~34 GB at stem 180, extrapolating past the H200's
143 GB before stem 300. This adds an opt-in eviction that bounds it.

Policy (from the NASSAU_R_STATS probe): a *degree* threshold, not LRU. Low-degree
R's are the stable hot core (reference span 0.99 of the run); high-degree R's are
scattered-recurring (0.81) and also the biggest matrices, so a degree cap evicts
the most bytes for the fewest references and is stem-independent (the kept set
saturates). Env NASSAU_GPU_RESIDENT_MAX_DEGREE (default i32::MAX = keep all).

Design — no kernel change. Every output row's products share one operation R
(extract_restricted), so hot (deg<=cap) and cold rows are DISJOINT. multiply_batch_on_gpu
partitions products, compacts each group to its own dense row range (so total
readback stays num_rows, not 2x), runs the resident pass on hot and a Transient
pass on cold, and scatters results back. Transient builds a per-block master with
create_from_slice, freed with the launch, so the device copy never persists. Cold
admissible data is cached host-side (COLD_HOST) so the expensive recompute is
one-shot, like the resident path. Fast path (cap==MAX or all-hot) is byte-identical
to before — zero regression by default. GpuProduct gains Clone for the partition.

Validated bit-exact: NASSAU_GPU_VERIFY full S_2 stem-110 resolution at theta=10
(transient path heavily stressed) and theta=100: mismatches=0 dx=0 panics=0.

Memory vs throughput (stem 180, ExclusivePages, streams=1): control (no eviction)
51 GB / 1349 s; theta=125 11 GB; theta=100 12 GB (master 34 -> ~1 GB), dx=0. But
eviction costs 2-4x throughput even at a high theta (stem 150: control 282 s,
theta=140 613 s) because the evicted high-degree matrices are the biggest and are
re-uploaded every launch. So this is a fit-in-memory lever for stems where the
master won't fit at all (a slow completion beats an OOM), tuned to the highest
theta that fits; it is NOT a speedup. Next optimization to cut the re-upload: a
bounded LRU device-side cold cache (upload once, reuse across a bidegree's launches).

Also retained: the R-access probe (NASSAU_R_STATS/dump_r_stats) and the device/host
[MEM] report (NASSAU_MEM_REPORT) used to characterize this — both env-gated no-ops.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
Foundational step toward generating admissible matrices (col_sums/masks) ON the
GPU into transient scratch instead of storing/uploading the resident master --
the direction that eliminates both the stem-300 device-memory wall and the
eviction re-upload cost (a given launch would enumerate its cold R's on-device,
never uploading them).

enumerate_admissible_ref reimplements AdmissibleMatrix::next using ONLY
flag-guarded control flow -- no break, continue, or early return -- because that
is the subset the cubecl DSL compiles cleanly (cf. multiply_pair, which tracks a
`rejected` flag rather than breaking). `found` replaces the odometer's `return
true`, `handled` replaces its `continue 'mid`. This validates the tricky
restructuring on the CPU, where it is fast to debug, before the hard-to-debug
#[cube] port -- which then becomes a mechanical transcription onto per-thread
local Arrays (state is tiny: rows = |p_part|, cols <= 32).

Test admissible_enum_ref_matches asserts bit-exact equivalence with
admissible_matrices over every real R up to degree 60 (4155 R's, all match).

Not yet ported to cubecl / not wired into the multiply -- next step.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
Transcribe enumerate_admissible_ref into a cubecl #[cube] kernel that
generates each R's col_sums/masks for every admissible matrix directly
into device scratch -- the on-GPU replacement for the resident/uploaded
master. One thread per distinct R, all state in fixed-size per-thread
local Arrays; flag-based control flow (while ... && !found / handled)
since the DSL has no break/continue/return.

Backend-agnostic host driver (generic over Runtime) so the identical
kernel can run on CUDA and the cpu backend. New test
admissible_enum_gpu_matches runs it on the H200 and asserts bit-exact
output (values + per-R counts) vs the CPU reference: 1055 R's, 30385
matrices, all matching.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
Wire enumerate_admissible_kernel into the Transient (evicted) cold path of
multiply_batch_block. A cold R's col_sums/masks are now GENERATED on the GPU
into transient scratch (one enumeration launch, ordered before the multiply on
the same stream) rather than built on the host and uploaded via
create_from_slice. Only the small p-parts + per-R dimensions upload; the scratch
is freed with the launch. This kills the per-launch H2D master re-upload that
made eviction a 2-4x slowdown (bench 2026-07-27), the whole point of the
in-kernel-enumeration direction.

Also replace the COLD_HOST full-array cache with COLD_COUNT, a 12-byte-per-R
(cs_len, mk_len, num_mats) shape cache. The evicted tail of the master now lives
neither on the device nor the host -- only its sizes, needed up front to lay out
the scratch offsets and the pair-count prefix sum. num_mats is counted once per
distinct R (admissible_matrices, arrays dropped) and memoized.

De-gate enumerate_admissible_kernel + ENUM_* caps + the MAX_XI_TAU import for
production. Validated bit-exact: the isolation test (1055 R's / 30385 matrices)
still passes, and NASSAU_GPU_VERIFY full S_2 stem-110 resolutions at theta=10
(nearly all R's cold -> enumeration path stressed) and theta=100 both report
mismatches=0 dx=0 panics=0 rc=0.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…refuted)

Investigation of the eviction crash + a throughput bench that together show
in-kernel enumeration cannot beat upload-based eviction:

- enumerate_admissible_kernel -> u64 addressing (was default u32); did NOT fix
  the big-block CUDA_ERROR_LAUNCH_FAILED, so the fault is elsewhere in the
  Transient wiring, not the kernel.
- admissible_enum_gpu_matches extended to degree 145, chunked per-degree:
  proves the kernel bit-exact vs the CPU reference (144903 R's, 185M matrices)
  across the full range the eviction path exercises -- so the kernel is correct.
- bench_admissible_cpu_vs_gpu (ignored): CPU admissible_matrices 1.61s vs GPU
  kernel-only 2.02s (0.8x, SLOWER) for degrees 1-130. GPU enumeration is ~3x
  slower than just transferring the same arrays (0.68s readback) it replaces:
  the odometer is sequential per R with matrices-per-R spanning 1..millions, so
  the launch bottlenecks on its few longest threads at GPU scalar speed. The
  "trade GPU integer work for PCIe bandwidth" premise is refuted -- enumeration
  worsens the eviction re-upload cost instead of curing it.

The enum wiring in multiply_batch_block's Transient path still faults at scale
and should be reverted to upload-based eviction; kept here as validated,
dormant code documenting the dead end. Default (theta=inf) path is unaffected.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
The earlier bench compared GPU enum kernel-only to the D->H READBACK (0.68s) as
a stand-in for the upload it replaces, and wrongly concluded enum is ~3x slower.
Production uses the matrices on-device (no readback), and the readback is not the
upload. Measure the actual H->D create_from_slice upload of the same host-built
arrays instead, synced via a tiny throwaway readback (upload-only):

  CPU admissible_matrices (1 core) : 1.60 s
  GPU enumerate in-kernel          : 2.02 s
  H->D upload of host-built arrays  : 1.85 s
  -> in-kernel enum is 1.09x the upload (parity), NOT 3x

So in-kernel enumeration is throughput-neutral vs uploading AND saves the host
COLD_HOST array cache (tens of GB at stem 300). The direction is viable, not
refuted; the big-block wiring crash is worth fixing.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
Prototype step 1 toward no-copy master growth, to kill the realloc-doubling
transient that pushes cubecl into its memory-corruption regime (the ~2x spike
at a 32->64 GiB master grow = ~96 GiB live, the root cause of the stem-140+
dx=nonzero cubecl ManagedMemoryDescriptor corruption).

A segmented master keeps old segments and appends a new fixed-size segment on
growth -- never copies -- so device peak is live_size + one segment, flat. The
kernel selects a segment by static branch (cubecl has no array-of-buffers):
offset o -> segment o/seg_elems, local o%seg_elems (seg_read_u16, MASTER_MAX_SEG
separate Array args). seg_gather_kernel + seg_read_matches_contiguous validate
the mechanic bit-exact (8 segments, scrambled indices) before touching the
multiply hot path.

Not yet wired into multiply_batch_kernel / resident_dev_handle -- next step.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…-doubling spike)

Replace the doubling-realloc-with-copy growth of every resident device
buffer with a single segmented, append-only, no-copy mechanism. Growth now
allocates only the new fixed-size segment(s) it needs and stage-writes the
tail into them; existing segments are never reallocated or copied, so the
device peak is `live + one_segment` instead of the ~2x transient (old+new
buffer both live) that pushed cubecl into its silent memory-corruption
regime -- the stem-140+ `dx != 0`.

One mechanism for all four resident buffers (master cs/mk, basis pp/ln):
- SegBuf + seg_grow! replace GrowBuf, resident_dev_handle!,
  basis_dev_handles!, stage_upload!, RESIDENT_REALLOC, RESIDENT_INIT_CAP
  (all deleted -- no two coexisting growth paths).
- multiply_batch_kernel binds each store as MASTER_MAX_SEG (=16) segment
  Arrays and GATHERS a thread's matrix cs/mk + its term p-part into small
  WORKING_CAP locals via seg_read_u16/seg_read_u32 (correct at any offset,
  so no layout padding), then calls the UNCHANGED pure multiply_pair with
  base 0. Only this kernel changed; multiply_pair and the test kernel are
  untouched.
- No realloc barrier: segments never change identity or get freed, so a
  reader's cloned handles stay valid across a concurrent append.
- master_seg_elems() (env NASSAU_GPU_MASTER_SEG_ELEMS, default 1<<31, < u32
  so a segment length never truncates cubecl's 32-bit metadata) => 64 GiB
  per buffer over 16 segments. Over-cap is a clean assert, not corruption.

Validated bit-exact on H200: seg_read_matches_contiguous (16-arg
primitive), multiply_batch_matches_reference across NASSAU_GPU_MASTER_SEG_ELEMS
768..4096 (multi-segment single-launch gather), and a new
multiply_batch_incremental_growth at seg_elems=8192 (cross-launch append
into the partially-filled last segment).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
Relocate every #[cfg(test)] item (the test-only cube kernels xor_f2 /
seqno_kernel / multiply_single_r_kernel / seg_gather_kernel and their host
drivers, plus enumerate_admissible_ref / _on_runtime, seg_gather_on_gpu,
check_enum_backend, enum_launch_timed) out of the production module body and
into `mod tests`, so the production path is no longer interleaved with test
scaffolding. Drops the now-redundant inner #[cfg(test)] attributes (the
module already carries it) and applies the project's nightly rustfmt (this
also formats the segmented-master rewrite from the previous commit).

Pure move; no behavior change. GPU correctness suite still green.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…per-launch churn)

Phase 1 of reducing the per-launch allocation/copy churn that drives cubecl's
allocator into CUDA_ERROR_LAUNCH_FAILED (719) at scale (memcheck showed the
fault is on cuMemcpyHtoDAsync/cuMemAllocAsync, not our kernels). Keeps the
exclusive-pages memory mode (the ~50% footprint win) and removes the churn it
caused, instead of switching to the sub-slice pool (which reuses buffers but
~2x memory -> OOM at the 90 GB master).

- g/xi: the seqno table and (constant) xi degrees are identical every launch at
  a given degree; upload them once to shared resident handles (RESIDENT_SEQNO,
  keyed by g.len()), re-uploading only on a degree bump, instead of a
  create_from_slice each launch. Synced before publish for cross-stream reads.
- out_h: reuse a persistent per-worker (= per-stream) XOR accumulator (OUT_ACCUM
  thread-local), grown only when a larger out_len appears, instead of empty()
  + free-via-memory_cleanup every launch (the single biggest churned buffer).
  Read back only the used [0,out_len) prefix via Handle::offset_end.

memory_cleanup stays for now (trims the still-churning per-R/per-product
metadata; phase 2 makes those persistent and drops it). Bit-exact vs the prior
path: multiply_batch_matches_reference, multiply_batch_incremental_growth
(5 growing launches: out_h regrow + g/xi re-upload), multiply_single_r.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
Revert the OUT_ACCUM thread-local output buffer from the previous commit: it
was per rayon worker (~100 threads), so it held ~100 persistent block-sized
buffers (~28-50 GB) and DEFEATED the existing out_h bounding — each launch's
out_h is already chunked to NASSAU_GPU_BLOCK_MB (512 MiB) and total in-flight
output is capped by the GPU_BUDGET permit (NASSAU_GPU_MEM_BUDGET_MB). Persisting
them per thread pushed device memory toward the OOM ceiling.

out_h returns to a per-launch empty() reclaimed by memory_cleanup (properly
chunked + budget-bounded). The resident g/xi seqno tables (uploaded once vs
re-uploaded per launch) are kept — a clean churn reduction. Bit-exact:
multiply_batch_matches_reference, multiply_batch_incremental_growth,
multiply_single_r.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…ply bench

When the cubecl CUDA context is poisoned by the unresolved uninit-handle bug
(CUDA_ERROR_LAUNCH_FAILED / ServerUnhealthy at the multiply readback,
milnor_gpu.rs:2074; tracel-ai/cubecl#1401), the whole context is dead — every
later launch fails, so per-call retry cannot recover. Wrap multiply_batch_on_gpu
in catch_unwind: on the first failure, set a process-wide GPU_DISABLED flag and
finish the resolution on the CPU via a new cpu_multiply_batch. The CPU output is
bit-identical to the GPU's (validated by cpu_multiply_batch_matches_gpu, including
the multi-block out_offset path a real module row uses). RREF runs on a separate
fp-cuda runtime and is intentionally not gated by this flag.

Add benches/nassau_milnor_gpu.rs: a GPU batched-multiply throughput / regression
bench (counterpart to nassau_milnor.rs) that hammers multiply_batch_on_gpu over an
output-degree sweep with no row-reduction or resolution machinery — for
`cargo bench --baseline` comparison across cubecl commits and for isolating a
multiply/allocator crash from the rest of the pipeline. Compiles to a no-op main
without the `gpu` feature.

Also pins cubecl/cubecl-common to the JoeyBF fork branch
(claude/pool-slot-map-v0.10.0) for the #1401 generational-slot-pool validation.
TEMPORARY: revert to the crates.io release before merging upstream.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…y_cleanup

The residual high-stem CUDA_ERROR_LAUNCH_FAILED in the GPU multiply is a cross-stream
pool-reclaim race: ~100 streams each calling client.memory_cleanup() every launch, one
stream's cleanup reclaiming a shared resident-master page still in flight on another
stream's kernel -> bad device pointer -> context poison (tracel-ai/cubecl#1401; the fork's
per-command retain_until_complete does not cover cross-stream shared-page reclaim).

Gate the cleanup call on NASSAU_GPU_CLEANUP_EVERY (default 1 = every launch, N = every Nth,
0 = never). Validated: with =0, S_2 stem-200 resolves fully clean on GPU (0 LAUNCH_FAILED,
0 dx) — the first clean GPU stem-200 — with device memory plateauing ~125 GB under cubecl's
internal pressure-triggered reclaim (which does NOT hit the race). A moderate throttle
(~16-32) is the likely stem-300 config: race-avoiding while bounding memory further.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…correctness guard

An #[ignore]d GPU soak: N worker threads hammer multiply_batch_on_gpu against one
growing shared resident master, across NASSAU_GPU_STREAMS streams with per-launch
memory_cleanup on — the cross-stream access pattern the stem-200 resolution drives.
Each result is checked against the bit-identical cpu_multiply_batch oracle (up to
NASSAU_SOAK_VERIFY_MAX, to keep the CPU precompute cheap); any mid-soak context
death flips GPU_DISABLED and fails the run.

Two jobs in one:
- Correctness: catches cross-stream renumber/identity races in seconds.
- #1401 reproducer (measured): clean at max_degree<=128, but at max_degree=160 with
  NASSAU_GPU_STREAMS=48 the cubecl cross-stream pool-reclaim race fires within a ~45s
  soak (never-initialized / ServerUnhealthy cascade, gpu_disabled flips) at only
  ~28 GB host / ~22 GB GPU — the genuine timing race, NOT an OOM. That's a ~1-2 min,
  low-memory stand-in for the 40-min stem-200 crash, and the gate the coming
  single-submission-thread redesign must turn GREEN.

Repro: NASSAU_GPU_STREAMS=48 NASSAU_GPU_CLEANUP_EVERY=1 NASSAU_SOAK_MAX_DEGREE=160 \
  cargo test -p algebra --release --features gpu -- --ignored --nocapture concurrent_growth_soak

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
Drop the JoeyBF/cubecl `pool-slot-map-v0.10.0` fork pin for upstream
`tracel-ai/cubecl` tag `v0.11.0-pre.1`, whose rewritten memory manager
(+1377 lines) and CUDA stream layer are the async-engine/pool subsystem
where the cross-stream reclaim race (#1401) lives.

Migrate milnor_gpu.rs to the 0.11 launch API (227 one-line replacements):
  - launchable array kernel args `&Array<T>`/`&mut Array<T>` → slices
    `&[T]`/`&mut [T]` (incl. `&mut Array<Atomic<u32>>` → `&mut [Atomic<u32>]`);
  - host-side `ArrayArg::from_raw_parts` → `BufferArg::from_raw_parts`
    (identical 2-arg signature, from `cubecl::prelude::*`);
  - local scratch `Array<T>` passed to `#[cube]` helpers now go through
    `Array::as_slice()`; non-launch helper params take `&[T]`.
Every `address_type = "u64"`/`"dynamic"` attribute and `launch_unchecked`
call is preserved verbatim — the high-stem u64-addressing correctness path
is untouched.

Effect at the harsher-than-production d=160/48-stream soak: break rate
3/5 (fork) → 1/5 (0.11-pre), independent of NASSAU_GPU_CLEANUP_EVERY, with
0 correctness mismatches across all runs. The residual fault surfaces as a
gentler CUDA_ERROR_ILLEGAL_ADDRESS at read_one (vs the old LAUNCH_FAILED
uninit-handle cascade), always caught by the CPU multiply fallback. The
race is reduced but not eliminated; the fallback remains the safety net.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…ord])

multiply_pair computed word = row_base + (out_offset + seqno)/32 and did an
unguarded atomic XOR into out[word]. When out_offset + seqno spans past the
row's num_limbs (which nassau_gpu::get_partial_matrix_restricted already
anticipates and masks on readback — bits >= target_dim are dropped), the
device write overran: into the next row (silent corruption) or past the
buffer (CUDA_ERROR_ILLEGAL_ADDRESS). compute-sanitizer on a malloc_sync
build confirmed "Invalid __global__ atomic of size 4 bytes ... out of bounds".

Thread num_limbs through multiply_batch_kernel and skip writes with
global_bit/32 >= num_limbs — a device-side mirror of the host's existing
defensive mask. Correctness-preserving; the skipped bits are exactly the
ones the host discards. The d=160/48-stream concurrent soak that reliably
broke ~3-5/5 now runs 0/5 with 0 correctness mismatches.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
The p-part of a Milnor basis element was a `Vec<u32>`, costing a heap
allocation and a pointer chase per element. At p = 2 the internal degree of
P(R) is sum_i r_i (2^i - 1) with non-negative terms, so r_i <= deg/(2^i - 1);
sizing each field by that bound packs the whole exponent sequence into 64 bits
for every degree up to 2045. At odd primes the same bound applies divided by
q = 2(p-1), so one layout serves every prime.

`MilnorBasisElement` is now 16 bytes, `Copy`, and entirely inline. Measured
over degrees 0..=300 at p = 2, `basis_table` drops from 252 MiB in 5,036,688
allocations to 77 MiB in none.

Three things fall out of the packing:

- The packed value is a canonical key, so the hand-rolled `MilnorHashMap`
  specialization for `not(odd-primes)` is gone; a plain `HashMap` now hashes a
  single word on every path. That code also assumed a degree bound of 1536
  without enforcing it. `compute_basis` now asserts the bound up front, which
  is what lets everything downstream skip range checks.
- Trailing zeros are not represented, so the "pop trailing zeros" loops after
  building a product disappear.
- `PPartMultiplier` no longer borrows its inputs, so its lifetime parameter is
  gone, and `PPartAllocation` loses the buffer it existed to recycle.

In `ext`, `MilnorSubalgebra`'s signature test becomes one masked comparison on
the packed word instead of a loop over entries, with the mask hoisted out of
`signature_mask`'s inner loop.

Two behaviour changes worth noting:

- `basis_element_from_string("P0")` and `("Sq0")` now return the identity
  rather than `None`. P(0) is the identity, and `AdemAlgebra::try_beps_pn`
  already special-cases `x == 0` this way; the old `None` came from `vec![0]`
  and `vec![]` hashing differently, an artifact of the representation.
- `increment_p_part` now carries before incrementing. The old order
  transiently stored `max[i] + 1`, which need not fit a field whose width is
  exactly saturated by `max[i]`. The enumeration is unchanged.

The observation that every Milnor exponent sequence up to degree 512 fits in 64
bits is due to Lixiong Wu; this implementation works out the widths, finds that
the same layout holds all the way to degree 2045, and carries it through the
algebra.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01XK8YCdHkCqUV9YM7D957hV
The first packing pass regressed `milnor_ppart` by up to 8% at odd primes and
mod 4, because assembling the answer went from a memcpy plus a vectorized add
to a per-entry read-modify-write through the checked `PPart::set`, and because
`PPart::get`'s range branch landed in `update`'s inner loop.

Two changes, both confined to the kernel:

- Assemble the answer in a plain `u64` and store it once. Entries are written
  in increasing index order into a value that starts at zero, so a shift and
  an `or` suffice; the range checks become debug assertions backed by
  `compute_basis`'s degree gate.
- Pad the layout tables to 16 entries so the private `PPart::entry` can mask
  its index rather than branch on it. Padded entries have width zero and so
  read as zero, which is the answer `get` would have returned anyway. The
  public `get` keeps its explicit check, since callers outside the multiplier
  index it with a q-part-derived length that is not bounded by `MAX_LEN`.

This recovers the regression (`ppart_4/a` and `ppart_3/a` back to baseline,
`ppart_4/b` -8%) and improves the Nassau regime further.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01XK8YCdHkCqUV9YM7D957hV
`basis_table` held a `MilnorBasisElement` per basis element. At p = 2 with
unstable support off, that element is exactly
`from_p(ppart_table[t][i], t)` -- the p-part again, with a q-part that is
always zero and a degree that is the index. It was a redundant copy.

Deriving it on demand costs nothing now that `MilnorBasisElement` is `Copy`
and 16 bytes: `basis_element_from_index` returns by value and builds it in
registers rather than handing out a reference into a table. The multiply
family takes the element by value for the same reason.

The table is still built at odd primes, where the q-part varies within a
degree, and when unstable support is on, where the basis is re-sorted by
excess. Neither is a re-wrapping of `ppart_table`.

Measured over degrees 0..=250 at p = 2 (1,958,958 elements), RSS growth from
`compute_basis` drops 125.0 MB -> 95.0 MB, i.e. 66.9 -> 50.8 bytes per element.
Projected to degree 500 that is 5.24 GB -> 3.95 GB. Unlike the ranker, this
needs no basis renumbering and costs nothing at lookup time.

A test verifies the derivation matches what the table used to hold, for every
element, so the redundancy is asserted rather than assumed.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01XK8YCdHkCqUV9YM7D957hV
`basis_element_to_index` runs once per term of every product, and is a hash map
storing an entry per basis element. The index it returns is a position in an
enumeration, so a canonical key alone cannot replace the map -- but the position
can be computed.

Let counts[i][d] be the number of exponent sequences of degree d using only
xi_1..xi_i. Splitting on whether r_i is zero gives the coin-change recurrence
counts[i][d] = counts[i-1][d] + counts[i][d - xi_i]. Ranking needs the number of
sequences with r_i > v, and substituting r_i -> r_i - (v+1) is a bijection onto
all sequences of degree d - (v+1)*xi_i, so that count is a single table lookup
rather than a sum. Walking the entries downward ranks a p-part in one lookup
each, against a table covering every degree at once, where the map it would
replace grows with the basis.

Whether that is worth it depends entirely on scale, which took some measuring to
see. Against the map, at p = 2:

    degree   per-degree map   hashmap    ranker    ratio
       120          0.10 MB    11.6us    26.7us    0.43x
       300          3.12 MB     792us    1384us    0.57x
       400         12.50 MB    4490us    5310us    0.85x
       500         37.50 MB   33118us   15865us    2.09x

A lookup probes only its own degree's map. While that fits in cache the map wins
easily: one hash round and one probe, against six to ten dependent table reads.
Once it does not -- the map is 37 MB in degree 500 -- every probe misses to DRAM
at ~33 ns, whereas the ranker's table is ~43 KB, stays in L1, and costs ~16 ns
regardless of degree. Benchmarking only up to degree 120, where the map is
0.1 MB, shows a 2x loss and hides all of this; the sweep here deliberately spans
the crossover. Tuning does not move the small-degree end: nested vs flat table,
a zero-padded prefix to drop the branch, and one- vs two-pass to break the
dependency chain were all measured, and the padded variant was worst, because
doubling the table pushed it out of L1.

So the two suit opposite ends of the range, and the ranker is on the right side
of the end where the algebra's memory is the problem worth solving: replacing the
map there is 3.3 GB smaller and 2x faster.

It stays off by default and unwired even when enabled, because it numbers the
basis in colex order rather than the order compute_ppart emits, which would
invalidate saved resolutions. That order is rankable in principle, but its
natural recursion has depth equal to the sum of the entries, which is worse than
hashing. The unstable path, which re-sorts each degree by excess, is not modelled
either.

Tests verify the table reproduces the algebra's own p-part counts and that the
rank is a bijection onto 0..dim in every degree, at p = 2 and p = 3, plus one
pinning down that it really does disagree with the current basis order.

Also measured and rejected: an `unrank` recovering the p-part at a given index,
which would let basis_element_from_index drop ppart_table entirely. It ran ~15x
slower than the array read it would replace, at every degree, with none of the
crossover above -- ppart_table is 8 bytes per element against ~43 for the map, so
it stays cache-resident. The bit-packing that makes rank worth having is the same
thing that makes unrank not. The likelier route, if it is ever revisited, is
enumerating the basis in index order, which is O(1) amortised and matches how
callers actually walk it.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01XK8YCdHkCqUV9YM7D957hV
Three paths computed with unvalidated input before checking it against the
packing bounds, so the intermediate arithmetic went wrong first. All three are
reachable from public, non-panicking entry points.

- `basis_element_from_string("P^s_t")` indexed the xi-degree table with `t`,
  which has exactly `MAX_LEN` entries, so `t = MAX_LEN` was out of bounds. `p^s`
  and the degree product could also overflow. Now `t` is bounded by the table
  itself and both are computed with checked arithmetic.
- `try_beps_pn` computed `q * x + e` before bounding `x`, which overflows for a
  large `x`. The bound moves above the computation.
- `MilnorSubalgebra::packed_signature` assumed the profile was no longer than
  `PPart::MAX_LEN` and that each signature entry fit its field. Neither holds:
  `SubalgebraIterator` grows a profile without limit and `from_bytes` reads
  whatever length a file gives. Out of range, `PPart::shift` returns 64 and the
  shift overflowed; an oversized entry silently spilled into the neighbouring
  field, which could select unrelated basis elements. It now returns `None` for
  a signature no element can have, and `signature_mask` yields nothing.

`basis_element_from_string` is documented as total and `try_beps_pn` is the
non-panicking half of `beps_pn`, so these were contract violations rather than
merely untidy. Tests cover each.

The signature test checks the packed mask against the per-entry comparison it
replaced, over every element up to degree 60, for profiles that are narrower
than their fields, wider than their fields, and longer than a p-part can be.

Also adds `basis_order_at_p2_is_stable`, which pins the first nine degrees to
fixed element names. The basis order is a wire format -- saved resolutions store
coefficients by index -- so it needs a guard that does not read from
`ppart_table`, which is the thing being guarded. Verified separately that the
order is unchanged from the base commit: identical for all 4156 elements in
degrees 0..=60.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01XK8YCdHkCqUV9YM7D957hV
CI lints with the nightly toolchain, where the unstable options in
`rustfmt.toml` -- `reorder_impl_items` among them -- actually take effect.
Stable rustfmt skips them with a warning, so `cargo fmt --check` passed locally
and failed in CI.

Formatting only: the constants are sorted and the blank lines between them
dropped.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01XK8YCdHkCqUV9YM7D957hV
The prior guard bounded the intra-row limb (out_offset + seqno spanning past
the row) but not the row itself. compute-sanitizer on a malloc_sync build at
higher degree (soak d=200) caught a second out-of-bounds atomic on the same
out[word]: when row_base overruns, word = row_base + limb lands past the
buffer even with a valid limb. Add the complementary `word < out.len()`
bound so no atomic write can escape the allocation by either route.

Validated at the d=200/48-stream soak: 4/4 clean, 0 crashes, 0 fallbacks,
0 correctness mismatches.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…ation

Two CUDA consumers share this GPU: the cubecl Milnor multiply and the fp-cuda
row reduction. Moving the reduction off a single global mutex onto per-thread
streams lets concurrent rayon workers overlap, but on its own it made a
stem-200 resolution die within 5-10 seconds with CUDA_ERROR_LAUNCH_FAILED,
reproducibly (3/3).

Isolating the two runtimes shows the fault needs both of them on one device:
multiply on GPU + reduction on CPU ran clean, and reduction on GPU + multiply
on CPU ran clean, while both together died every time. compute-sanitizer finds
no invalid access in either runtime — including with cubecl's allocations
forced synchronous so every buffer is tracked — so this is contention, not an
out-of-bounds bug. Giving fp-cuda its own non-primary CUDA context did not help
either (3/3 still died); the shared *device* is what matters, not the shared
context.

Overlap is also catastrophic for throughput, not just stability. The composable
reduction is a chain of thousands of tiny sequential per-column relaunches, so
sharing the device with the multiply's saturating kernels makes every launch
queue: the same reductions take 1.8-9.7 ms on an unshared GPU and 8.6-96.8 s
co-running, with nvidia-smi showing 99% SM at 10% memory utilisation (queueing,
not compute). The comment claiming this path "needs no cross-runtime exclusion"
had it backwards — being composable means it *can* overlap without deadlocking,
not that it should.

gpu_lock arbitrates: multiplies take the shared side and still overlap each
other; a large reduction takes the device exclusively for its ~10 ms. Total cost
is ~5 s of multiply pause across a whole stem-200 run. Writer preference is
required because multiplies are continuous and would starve the reduction
indefinitely.

Two properties are load-bearing and easy to get wrong:

- WHERE the shared guard is taken. Acquiring it at multiply entry deadlocks: the
  marshalling par_iter runs chunks on other workers, which steal another
  bidegree's multiply, block on the shared side behind a waiting reduction, and
  never let the original join finish. It is taken alongside the existing
  GpuPermit, past every rayon section, for exactly the reason documented there.
- Every wait is bounded, so an unforeseen cycle degrades to lost exclusivity
  rather than a hang. The bounds must exceed how long a reduction holds the
  device; at 25 ms multiplies barged back in mid-reduction and both the slowdown
  and the crashes returned.

With this, a stem-200 resolution completes on the GPU in 2h47m with 0 crashes,
where every prior attempt died. FP_CUDA_DEVICE / NASSAU_GPU_DEVICE put the two
runtimes on separate GPUs when more than one is available, which removes the
contention by construction and makes the arbitration a no-op.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…s time

The CPU fallback turned a hard GPU fault into a silent ~100x slowdown: a run
that had lost the GPU still reported "completed", so every measurement had to be
reconstructed by grepping stderr, and a crash at 5 seconds looked like a slow
success three hours later. GPU_DISABLED is still latched first so the soak test
can tell a context death from an ordinary assertion failure, but the panic now
propagates and the run dies at the fault.

The batch counters existed but take_batch_stats() had no callers, so the
dominant cost of a resolution was never attributed. Reporting them periodically
shows where multiply time actually goes, and the split matters because the
naive reading is wrong in two ways:

- The window called "marshal" spans GpuPermit::acquire and the gpu_lock
  acquisition, so it conflates host work with time parked on our own gates.
  Separating them shows the permit costs nothing measurable.
- The shares move by 2-3x as a run matures — early samples are small batches
  where fixed overhead dominates. At 2k launches it reads 29% prep / 47% lock /
  24% device; by 226k it is 14% / 9% / 77%. Only the mature numbers mean
  anything, and they say the multiply kernel itself is the bottleneck.

Reporting keys off the value fetch_add returns rather than a separate load: with
~100 workers a load races past exact multiples and the report can fire never
(observed: zero output over 12 minutes).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
The per-bidegree signature loop is the tail's critical path: at stem 200 the
bidegrees use A(3), whose signature space is 1024, and the steps run
sequentially on one thread — b=(200,6) spends 1076 s that way, 1024 steps at
~551 ms median plus one outlier at 470 s. With only 3-5 bidegrees in flight at
that point in the wavefront, parallelising this loop looked like the largest
remaining win.

It is not available. Each step reads dx.entry(v) over its own signature mask and
then writes dx with rows of the *unmasked* matrix, whose support extends past
that mask. NASSAU_PROBE_SIG_INDEP=1 snapshots dx before the loop and compares
every read against it: 6703 of 38348 reads (17.5%) differ, i.e. were perturbed
by an earlier signature, starting from bidegrees as small as (14,2). At p=2 a
differing entry flips the zero test that drives the step, so solving the
signatures against the pre-loop dx and combining would compute different
answers. The loop is a genuine forward substitution and must stay ordered.

Kept as a probe (zero cost unless the variable is set) because the alternative
is rediscovering this by writing the parallel version and getting wrong answers.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
Formatting only, no behaviour change; `just lint` runs `cargo fmt --all --check`.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
JoeyBF and others added 30 commits August 5, 2026 19:54
…caps

Three changes to `enumerate_admissible_kernel`, of which only one is a measured win.

Sort the enumeration's `R`s by matrix count before launch (+8%). The kernel walks an
odometer, so a thread's cost tracks its `R`'s matrix count, and a warp retires only when
its slowest lane does. Matrix counts are heavily skewed — the R_STATS Lorenz curve has the
hottest 1% of `R`s carrying 31% of references — so in basis order one huge `R` idles 31
lanes. Only the enumeration's own inputs are permuted: it writes each `R` at the absolute
offset it is handed, so the scratch comes out byte-identical and the multiply's
`r_cs_offset`/`prod_r_index` keep basis order untouched.

`enum_warp_utilisation` (new, CPU-only) puts the modelled headroom at 2.16x — 30.9% of
lane slots useful in basis order against 66.8% sorted. It cashed out as 1.08x measured,
because an SM keeps many warps resident and hides most of the modelled stall. The
docstring says so, so nobody sizes a decision on the model again.

`ENUM_COL_CAP` was `WORKING_CAP` (32), conflating this cap with the multiply kernel's
assembled-p-part array. `cols` is the widest bit-length of a p-part entry, bounded by the
`r_1` field width — 11. That made `matrix`, the hottest per-thread local array, 320 u32
instead of 110. Also zero only the region each `R` reaches rather than the full cap.
Both are strictly better on local-memory footprint and NEITHER is measurable end-to-end
(1.884s -> 1.872s, inside noise) — kept for correctness of the derivation, not for speed.
`enum_col_cap_bounds_real_rs` checks the derived caps against every real `R` to degree 400
(actual maxima: 8 rows, 9 cols).

For context, in-kernel enumeration already costs 0.81x what merely uploading the same
arrays costs, so there is no large multiple waiting here.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
The eviction path bound its per-block enumeration scratch as segment 0 alone, so a block
was capped at one segment — 4 GiB of u16. A stem-200 block wants 2.17e9 u16 of `masks`,
just over the 2^31 line, so `NASSAU_GPU_RESIDENT_MAX_DEGREE=125` died on a hard assert
partway through the run. That cap is why theta=125 could not complete.

The fix is not a bigger segment. `master_seg_elems()` is deliberately under `u32::MAX` so
a segment's element count cannot overflow cubecl's 32-bit array-length metadata — the same
truncation class that caused the earlier silent corruption — and raising it would buy 2x
and reopen that door. The resident path already spans MASTER_MAX_SEG segments through the
same `seg_read`, so the transient path just had to use them: the ceiling goes 4 GiB ->
64 GiB with no new truncation surface.

Layout places each `R` wholly inside one segment, with its `col_sums` and `masks` in the
same-numbered segment, advancing both cursors to the next boundary together when either
run would straddle. That keeps the enumeration as one launch per segment against the
existing single-buffer kernel signature, so no segmented-write kernel is needed and the
validated kernel is untouched. Offsets handed to each launch are segment-local; the
multiply keeps reading global offsets through `seg_read` exactly as before.

Verified by running the eviction test with NASSAU_GPU_MASTER_SEG_ELEMS shrunk to 1048576
and 262144, which forces many segments at degree 72 — all regimes still match the CPU
golden, as does the full suite at both sizes.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…and log it

`gpu_lock` arbitration was gated on `FP_CUDA_DEVICE == multiply_device()`, where
`multiply_device()` read `NASSAU_GPU_DEVICE` — a variable nothing in `algebra` reads any
more, left behind when the multiply became multi-GPU. It answered "device 0" however many
GPUs the multiply was actually saturating, so `FP_CUDA_DEVICE=2` on a 4-GPU node would
conclude "separate devices, no arbitration needed" while the multiply hammered device 2
too. The multiply shards across every visible device, so the test is
`FP_CUDA_DEVICE < multiply_devices()`.

This is a latent bug, NOT the cause of the theta=125 LAUNCH_FAILED: with default settings
both sides resolved to 0, so arbitration was already enabled. Verified by the log line
this adds, which exists because `[batch-stats] lock=` cannot distinguish "arbitration off"
from "arbitration on but uncontended" — and the answer turned out to be a third thing.

What the instrumentation actually exposed: `lock=0.0s` in every run *with arbitration
enabled*, because the multiply takes `gpu_lock::shared()` inside the SUBMISSION closure and
drops it when submission returns. Submission only enqueues; the kernels run long after. So
the multiply releases the guard while its saturating kernels are still executing, the
reduction then takes `exclusive()` believing the device is idle, and its thousands of tiny
sequential launches interleave with them — the exact overlap the lock exists to prevent.
Fixing that means holding the guard through the fence rather than the submit (the same
scope error as the earlier GpuPermit bug); left for its own change since it alters the hot
path's concurrency and invalidates the ~5 s/stem-200 cost estimate in the lock's docs.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…e multiply

`gpu_lock::exclusive()` early-returned whenever `arbitration_needed()` was false, which
conflated two independent questions:

  - does the MULTIPLY have to yield?  (only when it shares the reduction's device)
  - do REDUCTIONS have to serialize?  (always)

The second is unconditional. The row reduction's GEMM is a persistent whole-device grid —
`num_ctas = occupancy x SMs`, cluster-aligned, with cluster sync and DSMEM multicast — so
two concurrent reductions each demand the entire GPU and neither can be placed. On Hopper
that surfaces as a bare CUDA_ERROR_LAUNCH_FAILED, which compute-sanitizer does not
attribute (0 invalid accesses across a whole run: it was never a memory bug).

This is why putting the row reduction on its own GPU did not help on its own:
`FP_CUDA_DEVICE=3` with the multiply on 0..2 turned arbitration off wholesale, so
reductions stopped serializing against each other and the run still failed 63 times in
300 s. Isolating the device removes multiply contention and leaves reduction-vs-reduction
contention untouched.

With the split, an isolated reduction GPU runs clean at the FULL GEMM grid — no CTA cap,
no throughput sacrificed. Measured on the theta=125 stem-200 repro that faulted at ~105 s
in every other configuration: 400 s, 0 launch failures, 0 panics.

The alternative was capping the persistent grid (`FP_CUDA_GEMM_CTAS`). The sweep on
bench_kernel_only (16384^3, idle H200) shows why that is the wrong trade: the largest cap
that survives the workload is 32, and 32 CTAs is 2107 TOPS against 8674 at full grid —
24% of peak. 64 CTAs reaches 49% but still fails. Throughput is linear in CTA count to
~128 (97% of peak), so the kernel only needs the device it is not being given.

Not fixed here: the multiply takes `shared()` inside the SUBMISSION closure and drops it
when submission returns, while its kernels are still resident — so on a SHARED device the
yield is ineffective (`lock=0.0s` in every run). Correct fix is a dedicated fp-cuda driver
thread that fences between reductions and overlaps transfers (copy engines do not consume
SMs, so uploads can pipeline against compute). Isolating the reduction GPU sidesteps it.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
Both fp-cuda entry points launch persistent whole-device grids (`num_ctas = occupancy x
SMs`, cluster-aligned): the row reduction's trailing GEMM and the standalone `try_mul`.
Two cannot be placed at once, and on Hopper the loser does not queue — it fails with a
bare CUDA_ERROR_LAUNCH_FAILED that compute-sanitizer cannot attribute.

`gpu_lock::exclusive()` covered `row_reduce` only; `try_mul` was deliberately lock-free
("concurrent callers do not interfere"), so even a dedicated reduction GPU had two
independent whole-device consumers and was not actually owned by anything. That is why
enabling the cooperative reduction path on an isolated device still wedged after a single
reduction: the cooperative kernel spun at its grid-wide barrier for CTAs a concurrent
`try_mul` was holding.

Routing both through one thread makes single-ownership structural instead of a discipline
each new call site must remember. Jobs run to completion there, and both end in a
device-to-host download, so serialization is on COMPLETION, not submission — the
distinction that matters, and the one `gpu_lock::shared()` still gets wrong on the multiply
side (taken inside the submit closure, dropped when submission returns, which is why
`[batch-stats] lock=` reads 0.0s everywhere).

Measured on the theta=125 stem-200 repro, reduction isolated to device 3, FULL GEMM grid:
200 s, 120 GPU reductions, 0 launch failures, both with and without FP_CUDA_RR_COOP.
Default shared-device config unchanged: 147 reductions, 0 failures, max_t=245 in 200 s.

Cooperative mode is now SAFE but not a win at this workload (closed=20406 vs 20306 in
200 s — noise). The reduction path is too small a share of the resolution for its ~2x to
show. It stays off by default.

Two things this does NOT do. The driver takes no `gpu_lock` guard: serialization among
fp-cuda jobs is structural, and taking the guard there deadlocked the run (it waits for
the multiply's readers while workers block on the driver), so yielding to the multiply on
a SHARED device still needs arranging without a guard held across a blocking job.
Transfers are serialized with compute, though copy engines do not consume SMs, so a later
change can pipeline the next job's upload against the current job's kernels.

`contended_acquisition_terminates_and_writers_are_exclusive` now measures reader overlap
in its own uncontended phase: with `exclusive` unconditional, the contended phase keeps a
writer queued almost always and writer preference correctly holds readers off, so
asserting overlap there measured the scheduler rather than the lock.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
The row reduction's GEMM sized its persistent grid to full occupancy (`occupancy x SMs`,
~264 CTAs on an H200). That grid is only placeable on a GPU this process owns outright.
When anything else holds SMs the launch is not queued — it fails, as a bare
CUDA_ERROR_LAUNCH_FAILED that compute-sanitizer cannot attribute (0 invalid accesses
across a whole run: it was never a memory bug).

Arbitrating with the other runtime would make fp-cuda and algebra depend on each other's
scheduling, and would still not help a co-tenant in another process. Asking for a share
needs no such knowledge: the persistent loop already handles any multiple of CLUSTER
(fewer CTAs do more tile-iterations each), so grid size was never a correctness parameter.

Default 1/16. The safe size is NOT a sharp threshold: on the theta=125 stem-200 workload
1/16 ran a clean 240 s while 1/8 failed 74 times, and a 32-CTA cap — essentially the same
grid as 1/8 — had passed cleanly in an earlier run under the same nominal config. Where
the launch stops being placeable depends on what the other tenant is doing at that moment.
A share sharply reduces collision probability; it does not prove it to zero.

Cost, from the idle-device sweep (bench_kernel_only, 16384^3): 16 CTAs = 1062 binary TOPS,
32 = 2107, full grid = 8674. Set FP_CUDA_GEMM_DEVICE_FRAC=1 on a dedicated linear-algebra
GPU to restore full throughput.

Left unresolved, and recorded in the source: why an oversubscribed grid fails rather than
queueing. Ordinary launches schedule in waves; something here — thread-block clusters, the
dynamic shared-memory request, or both — makes placement a hard launch-time requirement.
Until that is understood no share can be called correct, only likelier to fit.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
Replaces the "UNRESOLVED: why does an oversubscribed grid fail rather than queue" note
with the answer, found in the GEMM's own history.

Phase 9 (264b411, thread-block clusters + TMA multicast of B) introduced
`__cluster_dims__`, B multicast with an all-ranks mask, and a cluster-wide empty barrier
reached via `mapa`. A cluster's CTAs are co-resident by construction — a hard placement
constraint an ordinary launch does not carry, and the reason the launch fails instead of
scheduling in waves.

Phase 8 (60b05cd) is not implicated, which is the useful half: its persistent loop
strides (`tile += gridDim.x`) and was already work-capped (`sms.min(total_tiles)`), so the
grouped rasterization that gives the kernel its L2 behaviour composes at any grid size.
Only Phase 9's cluster layer has to go.

That makes the fix the same one the row reduction already took (see `rr_coop`): a
cluster-free variant, composable by construction, needing no knowledge of any other tenant
— in this process or another. Cost is Phase 9's 2x cut in B HBM traffic, on top of the 8x
GROUP_M rasterization already delivers, against the ~8x this share currently sacrifices
(1062 TOPS at 1/16 vs 8674 at full grid).

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
The GEMM's `__cluster_dims__(CLUSTER,1,1)` makes a cluster's CTAs co-resident by
construction -- rank 0 multicasts B straight into its mates' shared memory and
every consumer arrives on their empty barriers through `mapa`. So the launch does
not queue for SMs the way an ordinary grid does; when another runtime holds them
it fails outright with CUDA_ERROR_LAUNCH_FAILED. That is the fault compute-
sanitizer could never attribute (0 invalid accesses across a whole run: it was
never a memory bug), and the reason `FP_CUDA_GEMM_DEVICE_FRAC` had to give away
~8x of the throughput to make the resolution survive a shared GPU.

Template the kernel body on the cluster width and emit two entry points. Every
cluster-dependent construct has an exact single-CTA counterpart -- rank is 0,
cluster_sync drops (the preceding __syncthreads already orders it), arrive_cluster
becomes a local arrive, multicast becomes a plain TMA load, the empty barrier
counts 1 -- so the two kernels run identical arithmetic on an identical tile
schedule and differ only in B's HBM traffic and in whether the launch demands
co-resident CTAs.

The cost of dropping multicast turns out to be nearly nil. bench_kernel_only on
an idle H200, cluster-free vs cluster, correct=true idempotent=true throughout:

     4096^3   4091 vs 4071 binary TOPS
     8192^3   6687 vs 6778
    16384^3   8501 vs 8632
    32768^3   9608 vs 9664

Within 1.5% everywhere. GROUP_M rasterization was already keeping B resident in
L2, so the second-order saving multicast adds does not show up at these shapes --
and it was never worth a hard placement constraint.

So the cluster-free kernel takes the whole machine by default (frac 1) instead of
a 1/16 share, which is a ~8x throughput gain over the stopgap on a shared GPU.
`FP_CUDA_GEMM_COOP=1` opts back into the cluster kernel for a dedicated device,
mirroring `FP_CUDA_RR_COOP` one layer down.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
… it loses

Both odometers -- `AdmissibleMatrix::next` and `enumerate_admissible_kernel` --
recompute an anti-diagonal bitsum `d` for every visited (row, col), and the same
loop already maintains `masks[row+col] = d | new_entry`. That looks like a
precomputation waiting to be exploited: collapse an O(cols) scan in the hottest
triple-nested loop down to one load, in both implementations at once.

Measured over 46,995,344 anti-diagonal computations (degree <= 120):

  d == masks           13,322,067   28.348%
  d subset-of masks    46,995,344  100.000%
  masks==0 (sound skip) 2,736,685    5.823%
  d-scan iterations    55,846,521 -> 1.19 per check

The containment is exact and universal, and it does yield a sound skip
(`masks[row+col] == 0` implies `d == 0`). But the scan it would skip averages
1.19 iterations: `(row+col+1).saturating_sub(rows)..col` is nearly always empty
or a single step for real R's, because they are short and wide relative to the
anti-diagonal. The skip fires on 5.8% of checks and avoids 6.9% of iterations of
a loop that barely runs. There is nothing here.

Keep the probe as an ignored diagnostic with the verdict in its docstring, so the
shortcut is not re-proposed. Its assert is armed: 0 containment violations.

With the earlier null result from cutting per-thread local state 3x (1.884s ->
1.872s), the odometer's arithmetic and its local state are both cleared. What is
left proportional to the work is the emit: ~20 scattered 2-byte global stores per
matrix, ~680 GB at degree 240.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…HBM bandwidth

Two null results had already cleared the enumeration arithmetic: shrinking the
per-thread local state 3x was noise (1.884s -> 1.872s), and the anti-diagonal
bitsum averages 1.19 iterations per check (ec3ed45). That left the emit as the
only candidate proportional to the work -- `cs_len + mk_len` (~20) scattered
2-byte global stores per matrix, with adjacent lanes writing to unrelated
`r_cs_out[ri]` offsets.

Add a comptime `emit` flag so the stores compile OUT rather than branching around
them, and time the identical odometer with and without. `mat` still increments,
so the enumeration and its trip count are unchanged; only the writes differ.

  degree <= 130: 89,392 R's, 75,987,575 matrices, 1.03e9 u16 stores (2.07 GB)
    emit=true   0.0858s
    emit=false  0.0491s
    stores: 42.8% of kernel time, 1.75x if free

Large, but NOT dominant -- the odometer is still 57%. The recoverable part is
what the bandwidth says: 2.07 GB in 0.0367s is ~56 GB/s against ~4.8 TB/s of HBM,
about 1% of what the device can do. That is the signature of scattered 2-byte
transactions, not of a store volume we are stuck with. Packing u16 pairs into u32
stores, or staging a matrix through shared memory for coalesced warp writes,
should recover most of the 42.8% without touching the format the multiply reads.

Production launches pass `emit = true`; only the diagnostic passes false.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
Replacing the four `Vec`s in `AdmissibleMatrix` with fixed arrays sized to
`PPart`'s structural caps (MAX_LEN = 10 rows, width(0) = 11 cols) removes four
allocations per R, of which there are ~1.6e6 at degree 240. It looked like a
free win. Measured on benches/nassau_milnor it is a net loss:

  op24xel32   6.7% slower
  op40xel1    2.2% slower
  op8xel24, op20xel24, op32xel8, op8xel8   0.3-0.6% slower
  op24xel24, op16xel32, op46xel1           0.1-0.3% faster
  rest        no significant change

The allocations were never the cost; the zeroing is. A typical R is ~4 rows x ~6
cols = 24 entries against a 110-entry cap, so `vec![0; rows * cols]` clears about
a third of what `[0; 110]` does, and the saved mallocs do not pay for the extra
stores. An end-to-end stem-110 CPU resolution agreed: 9.59s vs 9.71s, noise.

Same effect as shrinking the GPU kernel's ENUM_COL_CAP 32 -> 11, from the other
direction: over-sized fixed state costs more than the allocation it avoids. Noted
on the struct so the change is not re-proposed.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
f3391b7 measured the emit at 42.8% of the enum kernel and read its ~56 GB/s
as ~1% of HBM -- implying a bandwidth pathology from scattered 2-byte writes, to
be fixed by relayout. Two more probe modes say otherwise.

  emit=1 production                 0.0858s
  emit=0 no stores                  0.0490s
  emit=2 lane-adjacent (coalesced)  0.0829s   92.1% of the cost REMAINS
  emit=3 every other entry          0.0685s   53.1% remains

Mode 2 writes the same bytes and the same number of stores with the warp's lanes
hitting adjacent addresses -- perfect coalescing, garbage layout, timing probe
only. If transactions or bandwidth were the constraint that would have collapsed
the cost; it moved 8%. Mode 3 halves the store count and takes almost exactly
half the cost. The emit is bound by store-instruction issue.

So relayout is dead and packing is the fix, with width tracking the win directly
since cost is linear in instruction count:

  2 x u16 -> u32   half the stores      ~21% of kernel time   ~1.27x
  4 x u16 -> u64   quarter the stores   ~32%                  ~1.47x

Same bytes moved either way; what is bought is issue slots. Requires padding each
R's cs/mk stride to the packing width so a matrix's run starts word-aligned, and
teaching the multiply kernel's reader to unpack.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…head is steep

`NASSAU_R_STATS` ranked Rs by reference count, and the eviction policy was tuned
on that plus bytes. But enumeration cost is proportional to num_mats, so what a
rebuild actually costs is `count * num_mats`, and a small hot R can cost the same
as a big cold one. Record num_mats per R and report the cost-weighted
concentration, plus the transient set on its own.

S_2 (150,75), theta=125:

  total matrices enumerated (if every ref rebuilt)  385.5e9
  transient (deg>125)   98,551 Rs (57%)  43.5M refs (3%)  22% of cost
  within transient, cost coverage  top1%=21%  top5%=51%  top10%=68%  top25%=88%

Two things follow.

Transient is 22% of total enumeration work but ~99% of GPU KERNEL time (nsys,
transient-heavy config: enumerate_admissible_kernel 98.9%, multiply_batch_kernel
1.1%). No contradiction: resident Rs are enumerated once on the CPU at first
sight and then live in the master, so they never cost GPU time again. Rebuild
work and GPU time are different denominators and were being conflated.

And the cost inside the transient set is NOT flat -- 5% of its Rs carry half of
it. A partial cache has a head to exploit. Better, since bytes are ~num_mats and
cost is count*num_mats, cost-per-byte is EXACTLY count: num_mats cancels, so
admission ranked on plain reference count is optimal for any byte budget. The
degree cutoff remains the right first filter (degree proxies bytes) but discards
hot and cold alike among what it excludes; a small count-ranked cache layered on
top of theta is the missing piece, with memory bounded by what it is given
rather than by theta.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…s already optimal

The previous commit reported that transient cost is concentrated (top 5% of Rs
carry 51%) and proposed a count-ranked pinned cache. Add a per-R CSV dump and
simulate that policy against the reference stream. It does not work, and the
reason retires the idea.

  transient at (150,75): 98,551 Rs, 43.5M refs, 83.76e9 matrices, 8.2 GB total

  budget   oracle   k=1    k=2    k=4    k=8   k=16
    1%      2.7%   2.4%   2.5%   2.5%   2.5%   2.6%
    5%     11.8%  11.7%  11.7%  11.7%  11.7%  11.6%
   10%     21.6%  21.4%  21.3%  21.3%  21.2%  21.2%
   25%     44.0%  43.8%  43.8%  43.7%  43.5%  43.0%

Savings are LINEAR in bytes cached, and every admission rule lands within 0.3pp
of an oracle that ranks by cost-per-byte with full hindsight. When the oracle
cannot beat arrival order there is no structure to exploit.

The "steep head" was an artifact of measuring concentration per R rather than per
BYTE. cost/byte is exactly `count`, and the high-cost Rs are the big ones, so 5%
of Rs is nowhere near 5% of bytes. Normalised by the memory it costs -- the only
thing a cache budget cares about -- the distribution is flat.

So a pinned cache is the same linear memory-for-time dial theta already is, and
adds nothing. What the numbers do show is that theta is not a parameter to tune
but a fallback: the whole transient set is 8.2 GB at (150,75), and giving it up
costs ~99% of GPU kernel time (2412s uncapped vs 18566s at theta=125 for stem
200). Theta should be as high as device memory allows, and the levers worth
having are the ones that make it unnecessary -- a smaller master (values are
<= 2^11 but stored as u16; bit-packing is ~31% fewer bytes) and faster
enumeration for whatever still misses.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…gh as memory allows

`resident_degree_cap`'s docstring justified the threshold on a byte metric: at
theta<=125, 43% of distinct Rs stay resident and only 4% of references miss. That
is true and it is the wrong axis. A miss is not a fixed cost -- it re-enumerates
on the GPU in EVERY block touching that R (~442 times over a run at (150,75)) --
and enumerate_admissible_kernel is ~99% of GPU kernel time whenever the transient
path is live, against multiply_batch_kernel's 1.1% (nsys).

Stem 200 to max_t=310, all complete, 0 crashes:

  theta=inf   2412s
  theta=200   2865s   50.2 GB peak
  theta=125  18566s   44.9 GB peak

theta=125 gives up 6.5x in wall time for 5.3 GB. The knee is sharp and sits above
125; the memory curve is much flatter than the time curve. exec fell 16x
(26269s -> 1619s) and fence 32x (245398s -> 7596s) on identical work (pairs
5.88e13 both) -- the enumeration simply stops happening.

So the cap is a fallback for exhausting device memory, not a knob to tune down,
and the docstring now says so with the numbers. Also records the negative result
on smarter eviction: a pinned count-ranked cache saves linearly in bytes cached,
within 0.3pp of a full-hindsight oracle, because cost/byte is exactly `count` and
the per-byte distribution is flat.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
… as an OOM

Picking theta was a blind guess: master bytes are dominated by the top degrees
(stem 150: 2.2 GB at theta<=125, 9.9 GB at theta<=150), so a lower-stem run
cannot predict a higher theta -- the degrees carrying the mass are absent from
it. The only way to learn theta was to set it too high and crash hours in.

Track master bytes per R degree at admission (one add per first-sight R). The
cumulative sum IS "how big the master would be at theta = d", so any run reports
the whole curve, and a budget picks theta off it directly. Printed with the
periodic batch-stats as well as at the end, so a run that dies on an allocation
still leaves the answer behind.

  stem 140: full=4.1GB (1.0GB/GPU over 4 devices)
            theta<=100 0.3GB, <=125 1.5GB, <=141 4.1GB

Which also corrects the premise this cap was reasoned about under. The master is
a SMALL share of device memory -- 1.0 GB/GPU at stem 140, 2.6 GB/GPU at stem 150,
against a ~50 GB peak at stem 200. theta caps the master, so it is not the knob
that governs peak memory; the concurrent dense output matrices are. A theta that
fits the master is necessary, not sufficient, and lowering theta buys far less
memory than its 6.5x wall-time cost (e7c8d2f) suggests it should.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…n is empty

step1 computes `desired_image` as the kernel of the augmentation
`target_module -> cc_module` in degree t, by building a `target_dim x target_dim`
augmented identity and row-reducing it. When the target complex is empty in
degree t that map has a zero-dimensional codomain, so the kernel is the whole
space and no computation can discover otherwise. Take it directly.

`target_dim` is dim(A_t), so at high t this allocated and reduced a matrix large
enough to reach the GPU RREF path purely to rediscover the entire space.

The guard is on the codomain being empty, NOT on resolving the sphere: every
finite target module is concentrated in finitely many degrees, so past its top
cell this holds for all t -- almost the whole resolution. The sphere is only the
extreme of it, firing from t = 1.

Verified bit-identical output on S_2, C2 and C2_eta (stem 40, max_s 20) -- i.e.
on non-sphere modules, where the guard is false at low t and true above the top
cell.

NOT verified faster. At stem 120 the difference is noise (4.85s vs 4.91s): the
waste is O(dim(A_t)^2) so it only bites at high t, and the high-t A/B was killed
before finishing. Landed on correctness and on there being no reason to do the
work, not on a measured win.

Also gate the `dump_master_by_degree` call behind `feature = "gpu"` with its
neighbour -- it broke the non-GPU build.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…rd thread

seg_grow's MASTER_MAX_SEG check was an assert, and it runs on a shard's dedicated
per-launch thread. A panic there kills only that launch: the next spawns a fresh
thread and the resolution continues, having silently dropped the failed block's
products.

That is not hypothetical. On an uncapped stem-300 run all four shard threads hit
it at b=(174,26) and the run carried on for another ~900k batches (batch-stats
calls 828k -> 1740k) with GPU_DISABLED never set, because no fallback engaged --
nothing had failed from the resolution's point of view. The output would have
been wrong, not slow, and nothing in the run's own reporting would have said so.

There is no recovery from exhausting the segment table mid-master, so abort the
process. The message now points at the actual fix -- lower theta, sized off
`dump_master_by_degree` -- rather than only at raising the caps.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…ate does not fix it

The emit ablations (f3391b7, 1a7e0f9) said stores are 42.8% of this kernel
and that store COUNT matters while coalescing does not. Both true, and both
misleading about the mechanism. ncu on an idle H200:

  Compute (SM)  5.16%    Avg Active Threads/Warp  5.88 / 32
  Memory       13.92%    Active Warps/Scheduler   2.09
  DRAM          0.65%    No Eligible             76.57%
  L2 Hit       98.83%    Warp Cycles/Issued Instr 8.90

Nothing is saturated. Stores hit cache (L2 98.8%) and never reach DRAM, which is
why making them lane-adjacent bought 8%: there was no transaction pressure to
relieve. The kernel is latency-bound, and ncu attributes 45.6% of its 8.9 stall
cycles to L1TEX waits on the odometer's own state -- `matrix`/`totals`/
`col_sums`/`masks` are indexed by runtime values, so they live in local memory.

Tried the obvious fix: state in `Shared<[u32]>`, strided `elem * BLOCK + tid` so
a warp hits consecutive banks. Bit-exact (admissible_enum_gpu_matches passes),
and the mechanism is confirmed -- warp cycles per issued instruction 8.90 -> 5.92,
a 33% cut. But 152 u32/thread is 38.9 KB per block, which forces the block from
256 to 64 threads, and active warps per scheduler fall 2.09 -> 1.61. Net
109.63 ms -> 112.72 ms, a 2.8% LOSS. Reverted.

Two things worth carrying forward. The 33% is real and only paywalled behind the
shared footprint, so a smaller state (`matrix` as u16 -- values are <= 2^11) may
still collect it. And divergence is untouched by any of this: 5.88 of 32 lanes
active, ~82% of the machine idle, because one thread owns one R and warps run at
their longest member. That is now the largest single number in the profile.

Also: ncu's "Est. Local Speedup: 76.57%" is not a forecast. It is the prize if
the stall vanishes for free; here it cost occupancy and the trade was negative.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…y are not the constraint

Follow-up to the latency profile. ncu Occupancy/LaunchStats on the production
kernel (idle H200, degree <= 130):

  Registers Per Thread   40      Block Limit Registers   24 blocks
  Theoretical Occupancy  75%     Waves Per SM            0.44
  Achieved Occupancy    8.16%    Grid Size               1397 x 64 threads

Registers allow 24 blocks/SM and a 75% theoretical ceiling, so shrinking them
buys nothing -- the achieved figure is 8.16%, an order of magnitude below a
ceiling that is already generous.

Waves Per SM = 0.44 is the fact that matters. The H200 has 132 SMs x 24 block
slots = 3168; the grid supplies 1397. One thread per R over 89392 Rs is ~89k
threads against ~270k of device capacity, so the kernel cannot fill half the
machine however well it runs. Achieved occupancy then falls to 8.16% because
those few blocks drain raggedly -- the same divergence that leaves 5.88 of 32
lanes active.

Production is likely worse still: this benchmark enumerates every R to degree 130
in one launch, while a production launch covers only one block's transient Rs.

So the granularity is the bug: thread-per-R gives too FEW threads and wildly
UNEQUAL ones, and stores / local memory / shared memory / registers are all
adjustments to a kernel that has no work to hide latency with. That explains why
every attempt today measured between -3% and +8%. The changes that would matter
are structural: batch many more Rs per launch so the grid fills the device, or
replace thread-per-R with lane cooperation over dynamically pulled work.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…% of the grid

`enumerate_admissible_kernel` is ~99% of GPU kernel time, and profiling it
standalone gave Waves Per SM = 0.44: it could not fill half the device. But that
benchmark enumerates every R to degree 130 in ONE launch, whereas production
launches cover a single block-segment's transient Rs. Count what production
actually submits.

S_2 stem 150, max_s 60, theta=125:

  enum launches 10260   Rs/launch mean 1293 (max 34157)
  blocks/launch mean 6  waves/SM 0.002

Six blocks against an H200's 3168 slots. Production is 220x smaller than the
benchmark and uses 0.2% of the machine per launch, ~1293 threads out of ~270k of
capacity, 10260 times over.

That settles batching vs lane cooperation, which was the open question. Lane
cooperation (warp per R instead of thread per R) fixes divergence -- 5.88 of 32
lanes active -- and multiplies threads by 32, but the odometer is sequential so
it is Amdahl-capped near 1.6x, and the grids would still be a fraction of a
percent of the device. Batching Rs across launches attacks a factor of hundreds.

It also explains why this kernel dominates GPU time despite being simple integer
work: it is not slow, it is run 10260 times on an almost empty GPU, so the cost
is launch latency and drain. And it closes the micro-optimisation thread --
stores, local memory, shared memory, registers and coalescing all tune a kernel
that is idle by construction. The 42.8% emit share was real, but it is 42.8% of a
kernel at 0.2% grid utilisation.

The counters ride with `[batch-stats]`, so any run reports its own geometry.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…erial

Production submits ~6 blocks per enum launch (waves/SM 0.002), so the obvious
move was to spread the same work over more SMs with smaller blocks. Measured on
S_2 stem 150 / max_s 60 / theta=125:

  ENUM_THREADS  blocks/launch  waves/SM  wall
     256              6         0.002    561 s
      64             21         0.007    569 s
      32             41         0.013    542 s

6.5x the blocks for 3.4% — noise. Reverted.

The reason is the useful part. A launch's duration is set by its LONGEST SINGLE R
-- one thread, sequential odometer -- not by how many blocks it occupies, so
extra SMs idle beside the one thread still grinding. Rs/launch averages 1293 and
peaks at 34157, so the intra-launch spread is enormous.

And 10260 launches across a ~550 s run is ~50 ms apiece, which matches
`gpu_thread` running every device section on one thread and one stream: they are
strictly serial by design (that serialisation fixed a 370 s starvation bug).

This is what makes batching worth a factor of hundreds rather than a few percent:
merging launches converts a SUM into a MAX. Ten serialised 50 ms launches cost
500 ms; the same Rs in one launch cost ~50 ms, because the short Rs run beside
the long one instead of queueing behind it. Neither block size nor any per-thread
micro-optimisation can reach that, which is why every attempt has landed between
-3% and +8%.

Two candidates, recorded on the kernel: stream concurrency (cheap -- the launches
are already independent and FIFO dispatch could round-robin across N streams and
stay fair, but streams were pinned to 1 to fix a host-memory blowup from
per-stream pinned pools), or aggregating Rs across calls (needs enumeration
decoupled from the multiply that consumes its output).

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…bound

`NASSAU_GPU_STREAMS` (default 1, i.e. unchanged) lets N worker threads drain one
device's queue, each on its own `StreamId`, so independent device sections overlap
instead of running strictly one after another. This was the cheap half of the
"batch the enum kernel" plan: `enumerate_admissible_kernel` is ~99% of GPU kernel
time and runs at `Waves Per SM = 0.002`, and merging serialised launches turns a
SUM into a MAX.

Measured (S_2 stem 150, max_s 60, theta=125), and it refutes the plan:

| streams | wall  | peak host RSS | peak GPU |
|---------|-------|---------------|----------|
| 1       | 528 s | 39.1 GB       | 31.5 GB  |
| 2       | 527 s | 53.3 GB       | 43.9 GB  |
| 4       | 598 s | 75.9 GB       | 68.8 GB  |

Wall is flat at 2 and 13% worse at 4 while both memories roughly double. The memory
growth is the load-bearing part of the result: it proves the sections really did
overlap, so this is not "the streams did not engage" — concurrency happened and
bought nothing.

The conclusion is therefore about the workload, not about streams: at this
configuration the resolution is not GPU-throughput-bound. Kernel time is 99% enum,
but kernel time is not the critical path. Measured alongside: ~918% CPU on a
128-core node (~7% of the machine) with a wavefront only ~5 bidegrees wide, i.e.
the limiter is a serial dependency chain on the host.

Kept rather than reverted: a few lines, defaults to exactly the old behaviour, and
it is the control that makes the "not GPU-bound" claim falsifiable elsewhere.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
… a time

Both sides are little-endian packed F_2 bitvectors with bit i = column i, so four
of the kernel's u32 limbs ARE one of `fp`'s u64 limbs, byte for byte. The readback
can therefore be a truncating byte copy into a scratch `FpVector` plus one limb-wise
`add`, instead of walking the set bits and calling `add_basis_element` on each.

At the logged ~26% density that is cols/64 word XORs per row in place of ~0.26*cols
bounds-checked entry writes. `update_from_bytes` refills the scratch in place, so
the whole call allocates two buffers rather than one per row.

Applies to both `get_partial_matrix` and the restricted variant; the restricted one
also drops bits past `target_dim`, which is exactly what its old `col < target_dim`
guard did — whole limbs by truncation, the partial final limb by mask.

Verified with NASSAU_GPU_VERIFY over a full S_2 stem-40 max_s-30 resolution, which
builds every matrix both ways and compares nonzero column sets: 0 mismatches.

Honest timing result: no measurable win (528 s, identical to baseline at stem 150 /
theta=125). That is the same story as the streams experiment — the run is bound by a
serial dependency chain at ~7% CPU, so making CPU work cheaper does not move wall
time. Kept because it is verified-equivalent and strictly less work, and it will
matter once the serial chain is opened up.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
…, both negative

Ran the `NASSAU_PROBE_SIG_INDEP` probe that had been sitting unused: 19.2% of the
values a signature reads (27 028 of 140 614) were written by an earlier signature,
so the loop is a genuine forward substitution and the lift must stay ordered.

That only rules out reordering the lift. Everything above it reads read-only state,
so a windowed prepare stage feeding an ordered lift is legal, and it was built and
measured. It is worthless, for a reason worth recording so nobody rebuilds it:

    sum of per-bidegree TOTAL signature time  4064.7 s
    sum of per-bidegree MAX   signature time  3895.7 s   -> ceiling 1.04x

One signature is ~96% of its bidegree, so the ideal speedup is 4%. End to end a
window of 4 measured 537 s against a 528 s baseline while raising mean CPU from
918% to 1306% — 42% more CPU to lose 1.7%, exactly as the ceiling predicts. The
refactor is reverted; only the measurement is kept.

Context for what this closes off: the run uses ~7% of a 128-core node with a
wavefront ~5 bidegrees wide, and GPU stream concurrency was independently measured
to buy nothing. Three attempts to add concurrency have now failed, which locates
the limiter as latency on a serial chain rather than any throughput shortfall.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
… spin-wait

The zero signature is ~96% of its bidegree (it is unconstrained, so its mask is the
widest and its matrix the largest — reproducible in 100% of bidegrees across four
runs), and almost all of its work sat outside every span. Adds `cpu_restricted`,
`zs_masks`, `zs_select`, `zs_source_mask` and `zs_dx_init` to close that gap.

They are all small, which is itself the result — none of the CPU regions I suspected
is the cost:

    extract_restricted  1062.7 s      cpu_restricted   97.0 s
    gpu_readback         319.6 s      pair_prepass     90.1 s
    marshal_terms         54.3 s      zs_select         1.8 s
    gpu_submit             0.3 s      zs_dx_init        0.3 s

against a zero-signature `step` total of 4844.9 s, leaving ~66% unattributed. Span
accounting could not find it because `time.busy` is inclusive AND the thing being
waited on is not inside any span: `[batch-stats]` reports `fence=1024.2 s` (caller
blocked on the GPU completion fence) against a worker-side `exec=77.1 s`.

A sampling profile settles it. `perf record` on a steady-state stem-150 run, by
shared object:

    43.84%  libcuda.so      29.18%  the binary      17.40%  libc

The whole libcuda share is one tight address range (0x461860-0x461999): the driver's
spin-wait. The process burns ~44% of its CPU busy-waiting on the device.

That retro-explains every negative result in this branch of work. Extra CUDA streams,
a parallel signature-prepare window, and a cheaper matrix readback all failed for the
same reason — they add or cheapen CPU work on a path that is waiting for the GPU, and
the extra threads spin too. It also puts `enumerate_admissible_kernel` (99% of GPU
kernel time) back on the critical path as a LATENCY problem, which is a different
target from the throughput framing that streams were testing.

Top CPU symbols for whoever picks this up: `get_partial_matrix_restricted` 11.3%,
allocator (`_int_free`/`malloc_consolidate`/`malloc`/`_int_malloc`) 8.6%,
`__memcpy_ssse3` 5.9%, `__vdso_clock_gettime` 4.7%.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
… slower than uncapped

The stream sweep (and every other A/B in this session) ran at
NASSAU_GPU_RESIDENT_MAX_DEGREE=125. The same runs' own `[MASTER-BY-DEGREE]` line
reports the FULL master at stem 150 as 1.5 GB, 0.4 GB/GPU — so the cap was forcing
the transient re-enumeration path for a master that fits resident several times over.

Same binary, same session, stem 150 / max_s 60:

| theta | wall  | enum launches | mean CPU |
|-------|-------|---------------|----------|
| 125   | 575 s | 10260         | 1110%    |
| none  | 202 s | 0             | 1618%    |

2.85x, with `enumerate_admissible_kernel` not launching once. Two thirds of the wall
time being optimised all session was enum work this stem would never do.

This does not overturn the stream result (that arm is still flat-to-worse against its
own control), but it does bound what any theta-capped benchmark can say about where
time goes: at theta=125 it mostly says where the CAP puts time.

Records the caveat next to the numbers, and says to size theta from
`[MASTER-BY-DEGREE]` — cap only when the master genuinely does not fit.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
A perf profile of an UNCAPPED stem-150 run (the regime that matters — see the theta
note) is CPU-bound, not GPU-bound: libcuda falls to 2.96% of user cycles, and the
top costs are memcpy 15.3%, get_partial_matrix_restricted 14.4%, the allocator 16.1%,
resident_info 7.7%, row_reduce 7.6%, and 2.2% in DefaultHasher/RandomState.

That last one makes swapping the resident index (and COLD_COUNT, the per-launch
intern map, and the row remap) to `rustc_hash::FxHashMap` look like free money: the
keys are packed integers with no adversary. Measured, three replicates per arm:

    pre-change  205, 208, 195   mean 202.7 s
    FxHashMap   234, 210, 213   mean 219.0 s

~8% slower, so it is reverted. Fx is a weak multiply-rotate and evidently clusters on
packed `PPart` keys where SipHash spreads them.

The replicates are the other result. Run-to-run noise at stem 150 is ~3% (and ~9% at
theta=125 across sessions), which is larger than several single-sample deltas
measured earlier in this session — those should be read as "no effect", not as small
effects. Use three runs per arm here.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
The limb-wise readback wrote the GPU's `u32` limbs into the staging buffer one at a
time (`buf[o..o+n].copy_from_slice(&w.to_le_bytes()[..n])`) — a bounds-checked
4-byte copy per limb, thousands per row, for a region that is already byte-identical
on a little-endian target. A call-graph profile of an UNCAPPED stem-150 run put that
`copy_from_slice` at 12.15% of all user cycles, under
`restricted_partial_matrix_maybe_gpu`.

`fill_limb_bytes` views the row as bytes and does a single `memcpy`, keeping an
explicit big-endian fallback so a non-little-endian target stays correct rather than
silently wrong (`fp`'s own `limb::from_bytes` takes the same shortcut under the same
`cfg`).

Verified with NASSAU_GPU_VERIFY over a full S_2 stem-40 max_s-30 resolution: 0
mismatches.

Worth 5.0% on its own (199.7 s vs 210.3 s, interleaved, 3 rounds). See the following
commit for the other half.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
`term_indices` was a `Vec<usize>` written once and only ever read, but a product is
cloned twice on its way to the device: once to compact rows into a dense range per
hot/cold group, once to fan out into per-device buckets. Each clone duplicated every
term list, and the drops showed up as 5.45% of all user cycles in `_int_free` alone
(16.1% total in the allocator) in a call-graph profile of an uncapped stem-150 run.

`Arc<[usize]>` makes the clones a refcount bump and the drops O(1). `Arc<[T]>`
implements `FromIterator`, so every construction site still just `.collect()`s; only
two `for x in &prod.term_indices` loops needed an explicit `.iter()`.

MEASURED, interleaved A/B/C, 3 rounds (arms alternate so drift in machine load hits
all of them equally):

    A  control                    211, 214, 206   mean 210.3 s
    B  bulk-memcpy readback       188, 203, 208   mean 199.7 s   -5.0%
    C  B + Arc term_indices       149, 151, 135   mean 145.0 s   -31.1%

C wins every round with no overlap against A. 1.45x.

Method note, because it changed the answer: the same two changes measured
BLOCK-sequentially looked like regressions (+7%). Arms run in blocks confound the
effect with load drift, and this node drifts by more than 5%. Interleave, or do not
believe the number — an earlier `FxHashMap` verdict in this file was measured the
block way and should be re-checked before it is trusted.

Verified with NASSAU_GPU_VERIFY over a full S_2 stem-40 max_s-30 resolution: 0
mismatches.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01SmhK9czwJj2dZeemDhXtBf
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants