cuda: add a fused candidate-scoring kernel for CLM and Kev

[xiaoyu-xyz]

One CUDA kernel that turns a query vector and a matrix of candidate vectors into a probability distribution, without materialising the similarities in between.

Two models need this readout and neither has it:

Model	Similarity
CLM (#9)	"exp(logit_scale) * cos(state_head(s), action_head(c))", answered by "the softmax over the question's candidates"
Kev (#27)	"scale * (k(h_opts) @ q(h_decide)) with scale = 1/sqrt(256), softmaxed over that question's candidates"

Neither is in #19's gemm.cu, which is a bf16 cuBLASLt GEMM. The two differ in the similarity and in whether a projection happens first; the kernel takes the shape they share and leaves the projections to the caller, which is where #6 puts kernel selection.

Core code: 123 lines (candidate_scoring.cu). The oracle, the simulation, the parity harness, the benchmark and the tests are separate.

What is fused, precisely

Not merely "fewer launches". When the similarity is a cosine — CLM's case, not Kev's, which is a plain dot — each candidate's norm and its dot with the query need the same elements, so both accumulate in a single read of C. The unfused version composes a normalisation pass with a scoring pass and reads C twice. Beyond that, one thread block per question with the softmax over shared memory means no second kernel, no atomics and no global round trip for the similarities; K ≤ 255 is what makes the in-block softmax possible at all, and a batch is a grid over questions rather than a loop of launches.

What it is worth

bench.cu is the unfused baseline: the same block with the row norm loaded from a staged array instead of accumulated beside the dot, plus a separate norm pass. The protocol and both result sets are in docs/benchmarks/scoring-fusion/; the short version is that the benchmark found a defect the correctness checks could not.

At the 256 threads this kernel was first written with, the fusion lost: 0.82x at 255 candidates of 512 and 0.57x on Kev's shape. Achieved bandwidth was 1 to 190 GB/s against a device peak near 1000, so nothing was bandwidth-bound and there was no read to save — and a single question and a batch of eight took the same 53 µs, which is only possible if each question sits on its own SM. One thread block per question means a one-question request, the common serving case, occupies one of 128 SMs, with the norm accumulation on that block's critical path.

Raising the block to 1024 threads, the maximum on sm_89, reverses every result:

shape	similarity	K	D	questions	fused	unfused	ratio
clm-cosine-5x512	cosine	5	512	1	9.2 µs	11.3 µs	1.22x
clm-cosine-64x512	cosine	64	512	1	12.1 µs	18.3 µs	1.51x
clm-cosine-255x512	cosine	255	512	1	23.6 µs	42.2 µs	1.79x
clm-cosine-batch8-255x512	cosine	255	512	8	23.8 µs	43.8 µs	1.84x
kev-dot-255x2560	dot	255	2560	1	55.3 µs	100.3 µs	1.81x

Kev is the control: it has no norm to fuse, so its gain is entirely the second kernel and the global round trip it does not need. That it lands at the same ~1.8x is the useful part.

Through a CUDA Graph. The ABI is written for capture, so the same launch was also replayed from a graph. It saves a flat ~3 µs — 1.49x of the 9 µs smallest shape and 1.06x of the 52 µs largest — which explains the launch-bound floor at K = 5 but does not change which configuration wins at any shape measured. A graph and the fused kernel are complementary; the two effects multiply.

Verified on hardware

RTX 4090 (sm_89), CUDA 13.0 (V13.0.88), on two driver stacks, build.sh ./out 89:

driver	worst absolute difference (tolerance 1e-4)
595.58.03	1.835e-07
580.76.05	1.835e-07

gpu_parity.py runs four groups and all pass:

• 13 fixed cases, including the degenerate ones — a single candidate, identical candidates, a zero candidate — and the one whose raw logits overflow a naive softmax.
• The batch path: 1, 8, 4 and 3 questions with K up to 255 and D up to 256, each row compared with the reference and checked to sum to 1. The single-question call goes through cs_score_candidates and the rest through cs_score_candidates_batch, so the two entry points have to agree.
• Rejected arguments: K above 255, K of zero, D of zero, a zero scale, a negative temperature and a null query each return an error rather than a scored answer. K above 255 is the one that matters — the shared array is sized MAX_K, so scoring the first 255 of a larger set would be a confident wrong answer instead of a failure.
• CUDA Graph capture and replay, which the ABI was written for (it queues on the caller's stream and never synchronizes) but had never been checked: exactly one node is captured, and the replay matches a direct call exactly, not merely within tolerance.

compute-sanitizer --tool memcheck reports ERROR SUMMARY: 0 errors. test_scoring.py is 17/17 on CPU. gpu_parity.py and bench.py skip rather than fail when there is no library or no device, so a machine without a GPU reports that it could not check rather than that it passed.

Defects the GPU found, in order

1. The parity harness passed host pointers. A numpy buffer is host memory, so the kernel dereferenced it as device memory — cudaErrorIllegalAddress, not a wrong answer. compute-sanitizer located it as an invalid global read with consecutive lanes on consecutive addresses, which is the coalesced load pattern.
2. The query norm was counted warps times. The inner loop d = lane; d < D; d += WARP already covers every element of D across one warp's 32 lanes, but every warp computed the same partial and the partials were then summed across warps — so the norm came out sqrt(8) too large on the 8-warp block of the time. Every unnormalized case sat at float32 rounding while the normalize cases were off by ~1e-2, which is what pointed at it.
3. The occupancy defect above. Found by the benchmark, invisible to every correctness check, and fixed by a one-line change.

The second is worth recording: the simulation did not catch it. kernel_simulation.py computed the query norm correctly because it was written from the intent rather than transcribed from the .cu. A simulation is only as good as its fidelity, and the hardware is what settles it. The third is the same lesson one level up: being correct says nothing about being fast, and only a measurement says which.

The manifest, and CI

scoring.backend.json declares status: validated, abi_version: 1 and max_abs: 1e-4 — 545x of headroom over the measured worst case, so a failure is a real failure rather than tolerance noise.

It is written to the contract #25 defines, which is not on main yet. Checked against that branch's checker and drivers: check_contract.py reports 0 errors and 0 warnings, discover_buildable.py picks the backend up for the compile matrix, and run_gpu_parity.py selects it by status: validated and runs the reference.entrypoint the manifest declares. So the CUDA workflow in that stack covers this backend with no changes here; until it merges, these checks are run by hand and the repository's ci.yml, which is Rust-only, does not run them.