diff --git a/.agents/specs/vt-matmul-fp8-block-ref.md b/.agents/specs/vt-matmul-fp8-block-ref.md new file mode 100644 index 000000000..32822f14b --- /dev/null +++ b/.agents/specs/vt-matmul-fp8-block-ref.md @@ -0,0 +1,386 @@ +# VT-MATMUL-FP8-BLOCK-REF — the block-scaled FP8 GEMM, CPU reference arm + +Issue: [#1189](https://github.com/mudler/vllm.cpp/issues/1189), milestone **M2**. +Row: `VT-MATMUL-FP8-BLOCK-REF`. +Pinned oracle: vLLM `5559679229bc961848b121ccdeaa8fa5d79bec98` +(`.agents/upstream-sync.md`), asserted as the HEAD of the local checkout before +any `file:line` below was read. + +## Scope + +Add one `vt` op with a CPU kernel only: + +```c++ +void MatmulFp8BlockScaled(Queue& q, Tensor& out, + const Tensor& a_fp8, const Tensor& a_scale, + const Tensor& b_fp8, const Tensor& b_scale, + int block_n, int block_k); +``` + +`a_fp8` is `[M,K]` i8 carrying raw fp8-e4m3fn bytes, as `vt::QuantFp8Group` +emits. `a_scale` is `[M, cdiv(K, block_k)]` f32, as `vt::QuantFp8Group` emits. +`b_fp8` is `[N,K]` i8, the on-disk weight. `b_scale` is +`[cdiv(N, block_n), cdiv(K, block_k)]` f32, the on-disk `weight_scale_inv`. +`out` is `[M,N]` f32 or bf16. + +This is the numerical oracle every later block-FP8 kernel is measured against. +It is milestone M2 of #1189 and it deliberately stops there. + +**Out of scope, each owned by another milestone of #1189**: `Fp8BlockWeight`, +the loader rung and the config reader (M3); `layers::Fp8BlockLinearMethod` and +the Qwen3.5 wiring (M4); the mainloop-scaled CUTLASS kernel for `sm_121a` (M5); +merged `gate_up` and QKV (M6). No CUDA arm lands here: M5 owns it, and a +device-free host cannot gate one. + +## Which implementation actually runs, and why it does not change the answer + +M1 established that upstream's Triton source is not what executes for the +activation quant, and that the two arms disagree in polarity. The same question +has to be asked here before any arithmetic is mirrored, and the executing chain +was read end to end rather than inferred from the Python. + +The dispatch decision, in order: + +| Step | Where | +|---|---| +| CUDA block-FP8 kernel priority list: FlashInfer, DeepGEMM, **CUTLASS**, Marlin, Triton, Humming | `vllm/model_executor/kernels/linear/__init__.py:355-377` | +| DeepGEMM is auto-disabled for `qwen3_5_text` on device-capability family 120 | `vllm/utils/deep_gemm.py:27-46` | +| the selected kernel's apply: `ops.cutlass_scaled_mm(A, B.T, scale_a=As, scale_b=Bs.T)` | `vllm/model_executor/kernels/linear/scaled_mm/cutlass.py:312-326` | +| `cutlass_scaled_mm` routes every `sm >= 120` device to the sm120 entry | `csrc/libtorch_stable/quantization/w8a8/cutlass/scaled_mm_entry.cu:220-226` | +| the sm120 entry hands `cutlass_scaled_mm_blockwise_sm120_fp8` to the shared dispatcher | `csrc/libtorch_stable/quantization/w8a8/cutlass/scaled_mm_c3x_sm120.cu:13-20` | +| the dispatcher takes the blockwise branch whenever a scale is not per-tensor or per-row, requires both scales **f32** and 2-D, checks the shapes with **`ceil_div`**, and refuses a bias | `csrc/libtorch_stable/quantization/w8a8/cutlass/c3x/scaled_mm_helper.hpp:15-18,39-55` | +| the blockwise kernel's accumulator and scale element type are both `float`, and the two scale pointers are **mainloop** arguments | `csrc/libtorch_stable/quantization/w8a8/cutlass/c3x/scaled_mm_blockwise_sm120_fp8_dispatch.cuh:56-58,218-235` | +| CUTLASS 4.5.0, the mainloop line itself: `accum(i) += tmp_accum(i) * tCrScaleAViewAsC(i) * tCrScaleBViewAsC(i)` | `include/cutlass/gemm/collective/sm120_mma_tma_blockwise_scaling.hpp:714-717` | + +So the executing kernel is CUTLASS, and unlike M1 it **agrees** with the +reference in both placement and polarity: a per-K-block partial product is +formed in an f32 accumulator and **multiplied** by `a_s * b_s` before it is +added to the running sum. Mirroring `native_w8a8_block_matmul` is therefore +mirroring the executing arm, and that is a measurement rather than an +assumption: upstream compares the CUTLASS arm against exactly this reference at +`rel_diff < 0.001` (`tests/kernels/quantization/test_block_fp8.py:194-200`), a +gate no polarity or placement error survives. + +One difference is real and is not a divergence. The reference forms the scale +product first, `s = As_tiles[i] * Bs[j][i]` then `partial * s` +(`tests/kernels/quant_utils.py:150-151`); CUTLASS associates left to right, +`(partial * a_s) * b_s` (`sm120_mma_tma_blockwise_scaling.hpp:715`). The two +differ by at most one f32 ULP per K-block and upstream's own 1e-3 relative gate +is what admits it. We mirror the reference's association, because this op **is** +the reference port and M5's kernel will be measured against it. + +## The constraint that decides correctness + +**The scales apply in the mainloop, once per K-block, into an f32 accumulator — +not in the epilogue.** Written out (`fp8_utils.py:826-836`, which is the Triton +arm spelling the same structure the CUTLASS mainloop implements): + +```text +accumulator = zeros(f32) +for k_block: + a_s = a_scale[m, k_block] per token, per K-group + b_s = b_scale[n_block, k_block] per (N-block, K-block) + accumulator += dot(a_tile, b_tile) * a_s * b_s +``` + +Our existing per-tensor FP8 path folds one scalar `alpha` into the epilogue +(`vt::MatmulFp8Cutlass`, `src/vt/cpu/cpu_ops.cpp` `MatmulFp8CutlassKernel`: +`alpha * acc` after the whole K reduction). **An epilogue-only application +cannot express a per-K-block scale at all**: it has exactly one degree of +freedom per output element, and the block scheme has `cdiv(K, block_k)` of them. +This is a correctness constraint, not an optimisation choice, and it is the +reason M2 is a new op rather than a parameter of `kMatmulFp8Cutlass`. + +G4 is the instrument for it, and it is constructed so that no epilogue-folded +alpha can pass: it builds a case whose two K-blocks carry different scales and +whose correct output is not `alpha * (plain fp8 GEMM)` for any `alpha`. + +## Upstream anchors + +| What | Where | +|---|---| +| the reference this op ports, whole | `tests/kernels/quant_utils.py:91-154` | +| the fp32 compute type and the `.to(output_dtype)` store | `quant_utils.py:98,109-110,153` | +| `ceil` tiling of both dimensions | `quant_utils.py:123-124` | +| the `Bs` shape assertions, `n_tiles`/`k_tiles` | `quant_utils.py:125-126` | +| the `As` last-dimension `ceil` assertion | `quant_utils.py:115-116` | +| the mainloop accumulate, scales multiplied per K-block | `quant_utils.py:145-151` | +| the production wrapper's shape contract, `cdiv` on N and K | `vllm/model_executor/layers/quantization/utils/fp8_utils.py:928-936` | +| the Triton mainloop spelling the same structure | `fp8_utils.py:826-836` | +| the ported test, its grid, its tolerance | `tests/kernels/quantization/test_block_fp8.py:48-55,123-153` | +| the ragged-shape case upstream calls out by name (DSV3 `kv_a_proj_with_mqa`) | `test_block_fp8.py:156-200` | +| the activation quant that produces `a_fp8`/`a_scale` | `.agents/specs/vt-quant-fp8-group.md`, landed `ad5f175e7` | + +## Design + +### Numerics + +```text +for each output element (m, n): + acc = 0 f32 + for kt in [0, cdiv(K, block_k)): + part = 0 f32, a SEPARATE accumulator + for k in the k-tile (ragged final tile is short): + part += f32(a_fp8[m,k]) * f32(b_fp8[n,k]) + acc += part * (a_scale[m, kt] * b_scale[n / block_n, kt]) + out[m, n] = acc stored to out's dtype +``` + +`part` is a separate register that is scaled and then folded into `acc`. That is +the whole point: collapsing it into `acc` is what an epilogue alpha does, and it +is unrepresentable. + +`n / block_n` is integer division of the **output column index**, mirroring +`offs_bsn = offs_bn // group_n` (`fp8_utils.py:823`) and the reference's tiling +at `quant_utils.py:131-143`. `block_n` and `block_k` are validated positive +before either divides anything, because `x % 0` and `x / 0` are undefined +behaviour and a zero must refuse rather than trap. + +### Ragged edges + +`cdiv`, not floor, on both dimensions. Upstream's production wrapper asserts +`triton.cdiv(N, block_n) == Bs.shape[0]` and +`triton.cdiv(K, block_k) == Bs.shape[1]` (`fp8_utils.py:935-936`), so a final +short block is legal and must work. The final K-tile runs to `K`, not to +`(kt+1)*block_k`; the final N-block is short and its scale row still exists. +`N=576` (`4*128 + 64`) and `K=3884` (`30*128 + 44`) are the shapes from +upstream's own grid that expose an integer-division bug, and both are in the +ported grid, separately and together. + +### Memory format + +`out` is f32 or bf16 and nothing here is f32 by default. The accumulator is f32 +because upstream's `ElementAccumulator` is `float` +(`scaled_mm_blockwise_sm120_fp8_dispatch.cuh:56`) and the reference's +`compute_type` is `torch.float32` (`quant_utils.py:98`); the **store** rounds to +whatever `out` carries, exactly as `.to(output_dtype)` does +(`quant_utils.py:153`). Both scale tensors are f32 and the op refuses any other +dtype, because upstream refuses it too +(`scaled_mm_helper.hpp:15-18`). This op therefore widens nothing: it is the same +polarity `.agents/porting.md` requires. + +### Structure + +`vt::MatmulNvfp4Fp4` (`include/vt/ops.h:1477`) is the closest existing +signature — a block-scaled GEMM taking two packed operands and two scale +streams — and `MatmulFp8BlockScaled` follows its shape, with the scalar `alpha` +replaced by the block geometry it cannot express. `OpId::kMatmulFp8BlockScaled` +is appended before `kCount`, the additive convention documented at +`include/vt/ops.h:363-368`, so no existing op's id shifts. + +The kernel is a **correctness reference, not a performance path**, in the same +sense as `MatmulFp8CutlassKernel` beside it: a naive nest that makes the +block-FP8 seam resolvable on a CPU queue so M3 and M4 can be gated without a +GPU. It makes no speed claim. + +## Risks + +| Risk | Control | +|---|---| +| a later reader folds the scales into the epilogue, "simplifying" the inner accumulator away | the prose above and in the header, plus G4, which is unpassable by any single-alpha form | +| a ragged N or K uses floor tiling and silently drops or misindexes a block | `cdiv` everywhere; G2 runs `N=576`, `K=3884` separately and together, and G5 refuses a floor-sized scale tensor by name | +| the `n / block_n` index is computed from the tile index rather than the column, which agrees for round N and diverges for ragged N | G2's `N=576` case, where the two disagree in the final block | +| an all-zero output passes every value comparison | G2, G3 and G4 each carry a `nonzero == numel` vacuity guard | +| a scale dtype silently widens or narrows | G5 refuses a non-f32 scale by name | +| the test's reference is the implementation rewritten, so agreement proves nothing | the G2/G3 reference accumulates in `double` with a different loop nest and an exhaustively derived fp8 decode table, and its bound is an f32 forward-error bound rather than a fudge factor | +| `block_n`/`block_k` of zero divides by zero | validated positive before use; G5 asserts the refusal | +| nothing dispatches the op, so it is dead on arrival | acknowledged and named under `## Owed`: this is the staged-slice exception of `.agents/reachability.md`, and M4 owns the wiring | + +## Tests + +`tests/vt/test_ops_matmul_fp8_block_cpu.cpp`, registered in +`tests/CMakeLists.txt`. + +- **G1** the registration itself: `OpRegistered(kMatmulFp8BlockScaled, kCPU)`, + and `OpName` is not `unknown`. A refusal case cannot stand in for this — M1 + measured a deleted registration passing every refusal test it had. +- **G2** the ported upstream case, `test_block_fp8.py:123-153`, against an + independently written `double` reference, carrying upstream's own + `rel_diff < 0.001` criterion verbatim **and** a tighter per-element f32 + forward-error bound, because our arm is the reference rather than a kernel + measured against it. Vacuity guard. +- **G3** the ragged grid, run as part of G2's table and called out here because + it is the reason the grid is what it is: `N=576`, `K=3884`, and the pair + together, plus upstream's dedicated DSV3 case `M=32, N=576, K=7168` + (`test_block_fp8.py:156-200`). +- **G4** **the mainloop constraint, made unpassable by an epilogue**. Two + K-blocks with deliberately different scales and a hand-computed expected + output; then the same operands with the two K-block scales *swapped*, which + leaves every per-tensor summary of the scales identical and changes the + correct answer. A kernel that folds one alpha produces the same output for + both and fails. Vacuity guard. +- **G5** the refusals, each by name: a non-f32 scale, a `b_scale` sized by + floor instead of `cdiv` on N and on K, an `a_scale` with the wrong number of + groups or rows, a zero or negative `block_n`/`block_k`, a rank mismatch, a + non-contiguous operand, a device mismatch, a non-i8 packed operand, and an + `out` that is neither f32 nor bf16. +- **G6** the M1 seam: `vt::QuantFp8Group` feeding `vt::MatmulFp8BlockScaled` on + a CPU queue, end to end, which is the pair a block-FP8 linear method will run. + This is a *composition* test, not a reachability claim; nothing in production + calls it yet and `## Owed` says so. + +### The adaptation of the upstream grid, and why it is unavoidable + +Upstream's grid is `itertools.product(M, N, K, ...)` with `M=[1,7,8,83,4096]`, +`N=[128,512,576,7168,13824]`, `K=[256,3884,4096,13824,16384]` +(`test_block_fp8.py:48-55`) — 125 combinations, the largest of which is +`4096x13824x16384`, about 9.3e11 multiply-accumulates. That is a GPU grid. A +naive CPU reference nest cannot run it, and running it twice (op and reference) +is worse. + +The adaptation preserves the **axes**, not the product: every value of every +axis appears at least once, the ragged values appear separately and together, +and the total is under 1e9 MACs. Which values are paired with which is the only +thing dropped, and the parameters, dtypes, tolerance and failure criterion are +preserved exactly. The grid is written out in the test with this reasoning +beside it. + +## Gates + +| Gate | Command | +|---|---| +| focused | `ctest -R test_ops_matmul_fp8_block_cpu --output-on-failure` | +| op provider totality | `ctest -R test_op_provider` | +| the M1 sibling, unchanged | `ctest -R test_ops_quant_fp8_group_cpu` | +| the per-tensor sibling, unchanged | `ctest -R test_ops_fp8_cpu` | +| record | `scripts/agent-preflight.sh --fail-on-skip` | + +No GPU lease is taken and none is needed: the op has a CPU arm only. + +## Owed + +- **Nothing reaches this op.** `vt::MatmulFp8BlockScaled` is dispatched by no + production entry point at this merge commit: `include/vllm.h` does not expose + it, no loader builds an `Fp8BlockWeight`, and `ModelRegistry::Forward` has no + block-FP8 linear method to call it from. The wiring is owned by #1189 + milestone M4 (`layers::Fp8BlockLinearMethod` and the Qwen3.5 dense forward), + which needs M3 first. This is the staged-slice exception of + `.agents/reachability.md`, named here, in the commit body, and in the pull + request body. +- **The CUDA arm.** There is none. Milestone M5 owns the mainloop-scaled CUTLASS + kernel for `sm_121a`, needs a GPU, and will be measured against this op. +- **The `swap_ab` path.** Upstream's sm120 blockwise dispatch swaps the operands + for some shapes (`scaled_mm_blockwise_sm120_fp8_dispatch.cuh:221-235`), which + changes the scale-pointer assignment but not the mathematical result. It is a + kernel-scheduling concern with no reference-arm meaning. Owed by M5. +- **The column-major and TMA-aligned activation-scale layouts** + (`fp8_utils.py:610-628`), which the CUTLASS kernel reads and this reference + does not. Already owed by M5 in `.agents/specs/vt-quant-fp8-group.md`; repeated + here because M2 is the first op that could have consumed them. +- **`bias`.** Upstream refuses a bias on the blockwise path outright + (`scaled_mm_helper.hpp:55`). This op takes none, which mirrors the refusal + rather than deferring a feature. + +## Stop conditions + +Stop and report `NEEDS_DECISION` if any of the following holds. + +- The pinned oracle's checkout is not at + `5559679229bc961848b121ccdeaa8fa5d79bec98`. Every anchor above was read at + that revision, asserted before the first read. +- G4 cannot be made to fail an epilogue-folded implementation. If a single alpha + can reproduce the mainloop result on the constructed case, the case is wrong + and the constraint is untested, and widening the case is not the fix. +- The op cannot express upstream's contract without a parameter that no caller + in this tree passes. + +Stop and report `NEEDS_CONTEXT` if the work requires a GPU lease. The row is +scoped so that it does not. + +## Evidence + +Taken on the merged tree, `origin/main` at `2d26da5a1`. + +**RED.** With the test present and the implementation reversed out — the reverse +of the implementation diff applied to the four source files, `git diff --stat` +printing `4 files changed, 200 deletions(-)` so the revert is proven to have +landed — the focused build fails, `compile_rc=1`, with 29 errors: 23 +``'MatmulFp8BlockScaled' is not a member of 'vt'`` and 6 +``'kMatmulFp8BlockScaled' is not a member of 'vt::OpId'``. Restoring the four +files returns `git status --porcelain` to empty and the build to green. + +**GREEN.** `test_ops_matmul_fp8_block_cpu` reports 6 cases, 80 assertions, 0 +failed, in 3.5 s under `ctest`. Per block, each run through a `-tc` prefix +filter that contains no comma, because doctest splits `-tc` on commas and a +name that contains one yields `0 cases ran` under a `SUCCESS!` banner: + +| Block | Cases | Assertions | +|---|---:|---:| +| G1 registration and name | 1 | 2 | +| G2 the ported grid | 1 | 22 | +| G3 the ragged edges | 1 | 13 | +| G4 the mainloop constraint | 1 | 14 | +| G5 the refusals | 1 | 25 | +| G6 the M1 seam | 1 | 4 | +| **sum** | **6** | **80** | + +The buckets sum to the whole-run count, so no block is silently empty and no +filter selected nothing. + +The other declared gates on the same tree: `test_op_provider` passed, +`test_ops_quant_fp8_group_cpu` passed (M1 unchanged), +`test_ops_fp8_cpu` passed (the per-tensor sibling unchanged). +`scripts/agent-preflight.sh --fail-on-skip` reports **All gates green** with 80 +`ok` results, no `FAIL` and no `SKIP` before the verdict line. +`test_cpu_x86_llamacpp_floor` (#618) passed on this run rather than reporting +`NO_QUIET_WINDOW`, so this row needed no pristine-baseline reproduction. + +### The ragged edge, specifically + +`N=576` is `4*128 + 64` and `K=3884` is `30*128 + 44`, the two non-round shapes +from upstream's own grid (`test_block_fp8.py:49-50`). G3 runs each separately +and both together, plus upstream's dedicated DSV3 `kv_a_proj_with_mqa` case +`M=32, N=576, K=7168` (`test_block_fp8.py:156-200`). This is load-bearing and +measured, not asserted: mutation M2 (floor `k_tiles`) and mutation M3 (a floor +N-block index) each leave **G2 entirely green** and fail only G3, because every +`N` and every `K` in G2 is a multiple of 128. A grid of round shapes passes +while being wrong. + +### Mutation results + +Every mutation printed `compile_rc` **and** `git diff --stat` before the run, +because a mutation that fails to build and a mutation that never applied both +read as a passing test, and M1 of #1189 hit each of those once. Each was +restored with `git checkout --` and the restore verified. + +| Mutation | `compile_rc` | Result | +|---|---|---| +| the scales folded into ONE epilogue alpha, `part` collapsed into `acc` | 0 | 4 of 6 cases fail, 26 of 80 assertions. G4's `got != got_swapped` fires, which is the epilogue signature exactly | +| `k_tiles = k / block_k`, floor instead of cdiv, in the kernel | 0 | **only G3 fails**, 4 assertions. G2 stays green: every K in it is a multiple of 128 | +| the b_scale row index forced to 0 | **1** | proves nothing: `-Werror=unused-parameter` on `block_n` | +| the same defect with the parameter kept live: `nb = min(col / block_n, n / block_n - 1)`, a floor tile count that makes a short final N-block reuse the previous scale row | 0 | 3 of 6 cases fail, 13 assertions: G3, G4, G6. **G2 stays green**, for the same reason | +| `a_scale` dropped from the scale product | 0 | 4 of 6 cases fail, 24 assertions | +| the CPU registration deleted, the kernel symbol kept live | 0 | all 6 cases fail after 7 assertions. G1's `REQUIRE(OpRegistered(...))` is what catches it | +| `n_tiles = n / block_n`, floor, in the wrapper's validation | 0 | 4 of 6 cases fail. G5's floor-sized `b_scale` refusal stops firing and G3's well-formed ragged calls are refused instead | +| the `block_n > 0 && block_k > 0` check removed | 0 | **SIGFPE, core dumped, `run_rc=136`**. The check converts undefined behaviour into a named refusal | +| the `a_scale` shape check widened to a rank check | 0 | G5's narrow and tall `a_scale` refusals both stop firing, 2 assertions | +| the store zeroed, `acc * 0.0F` | 0 | 4 of 6 cases fail, 29 assertions, including G4's two vacuity guards — the guards are live | + +Three results are worth keeping. + +**The epilogue mutation is the row's whole claim, and it is now measured.** +Folding the per-block scales into a single alpha compiles, runs, and produces the +identical number for a scale tensor and for that same tensor with its two +K-block entries swapped. G4 is the only block that can see that, and it does. + +**A grid of round shapes is blind to the ragged bug in both directions.** Two +different floor-vs-ceil mutations — one on the K tiling in the kernel, one on the +N-block index — left G2 completely green and were caught only by G3. That is the +argument for keeping `N=576` and `K=3884` in the grid rather than the shapes the +target checkpoint happens to use, which are all multiples of 128. + +**A refusal removed can crash rather than report.** With the block-size +positivity check deleted, the process took SIGFPE on integer division and doctest +had already printed `assertions: 52 | 52 passed | 0 failed` before it died. The +run's exit status is the verdict; its printed summary is not. + +### A false red, recorded because it cost a cycle + +Mid-run, a direct invocation of the test binary reproduced the zeroed-store +mutation's exact signature — 4 failed cases, 29 failed assertions — against a +tree whose sources were clean. The binary was the one the mutation harness had +last linked; the whole-tree build that followed had not reached the test target +before the run. Nothing was wrong with the code. The control is the one +`.agents/verification.md` already names: rebuild explicitly, require +`ninja: no work to do`, and check that the executable's mtime is newer than +every source it links, before believing any verdict about it. diff --git a/include/vt/ops.h b/include/vt/ops.h index 7905f4da4..a4142f6af 100644 --- a/include/vt/ops.h +++ b/include/vt/ops.h @@ -433,6 +433,14 @@ enum class OpId : uint8_t { // See vt::QuantFp8Group below for the contract. // Appended before kCount so no existing op's id shifts. kQuantFp8Group, + // --- Block-wise FP8 (VT-MATMUL-FP8-BLOCK-REF, #1189 milestone M2). The + // 128x128 block-scaled fp8 GEMM that consumes kQuantFp8Group's output. It is + // NOT a parameter of kMatmulFp8Cutlass and cannot be: that op folds ONE + // scalar alpha into the epilogue, which has exactly one degree of freedom per + // output element, while this scheme has cdiv(K, block_k) of them and applies + // them in the MAINLOOP. See vt::MatmulFp8BlockScaled below for the contract. + // Appended before kCount so no existing op's id shifts. + kMatmulFp8BlockScaled, kCount }; @@ -936,6 +944,13 @@ using QuantFp8StaticFn = void (*)(Queue&, Tensor&, const Tensor&, float); // and the f32 [M, K/group_size] per-group scale. using QuantFp8GroupFn = void (*)(Queue&, Tensor& /*out_fp8*/, Tensor& /*out_scale*/, const Tensor& /*x*/, int /*group_size*/); +// Two scale streams and the block geometry, where the per-tensor fp8 GEMMs above +// carry one scalar alpha. The geometry is not a convenience: it is what selects +// which scale pair each K-block multiplies by. +using MatmulFp8BlockScaledFn = void (*)(Queue&, Tensor& /*out*/, const Tensor& /*a_fp8*/, + const Tensor& /*a_scale*/, const Tensor& /*b_fp8*/, + const Tensor& /*b_scale*/, int /*block_n*/, + int /*block_k*/); using RmsNormQuantFp8Fn = void (*)(Queue&, Tensor& /*out_fp8*/, Tensor* /*out_bf16*/, const Tensor& /*x*/, const Tensor& /*weight*/, const RmsNormArgs&, Tensor* /*residual*/, float /*input_scale*/); @@ -1582,6 +1597,69 @@ void QuantFp8Static(Queue& q, Tensor& out_fp8, const Tensor& x, float input_scal void QuantFp8Group(Queue& q, Tensor& out_fp8, Tensor& out_scale, const Tensor& x, int group_size); +// MatmulFp8BlockScaled (VT-MATMUL-FP8-BLOCK-REF, #1189 M2, +// .agents/specs/vt-matmul-fp8-block-ref.md) — the 128x128 block-scaled fp8 GEMM, +// mirroring native_w8a8_block_matmul (tests/kernels/quant_utils.py:91-154): +// +// for each (m, n): +// acc = 0 f32 +// for kt in [0, cdiv(K, block_k)): +// part = 0 f32, a SEPARATE register +// for k in the k-tile: part += f8(a[m,k]) * f8(b[n,k]) +// acc += part * ( a_scale[m, kt] * b_scale[n / block_n, kt] ) +// out[m, n] = acc stored to out's dtype +// +// THE SCALES APPLY IN THE MAINLOOP, ONCE PER K-BLOCK, INTO AN F32 ACCUMULATOR — +// NOT IN THE EPILOGUE, and that is a correctness constraint rather than an +// optimisation choice. MatmulFp8Cutlass above folds one scalar alpha after the +// whole K reduction. An epilogue has exactly ONE degree of freedom per output +// element; this scheme has cdiv(K, block_k) of them. An epilogue-only +// application therefore cannot express a per-K-block scale AT ALL, which is why +// this is a separate op. DO NOT "simplify" `part` away into `acc`: that IS the +// epilogue form, and tests/vt/test_ops_matmul_fp8_block_cpu.cpp G4 is built so +// that no single-alpha implementation can pass it. +// +// WHICH UPSTREAM ARM THIS MIRRORS. The Triton kernel at fp8_utils.py:826-836 is +// not what executes on the target architecture; CUTLASS is +// (vllm/model_executor/kernels/linear/__init__.py:355-377 ranks it third, +// DeepGEMM is auto-disabled for qwen3_5_text on family 120 at +// vllm/utils/deep_gemm.py:27-46, and Marlin is excluded at cc >= 89). Unlike the +// QuantFp8Group case above, the two AGREE: csrc/.../c3x/ +// scaled_mm_blockwise_sm120_fp8_dispatch.cuh:56-58,218-235 hands both scale +// pointers to the MAINLOOP arguments over an `ElementAccumulator = float`, and +// cutlass 4.5.0's sm120_mma_tma_blockwise_scaling.hpp:714-717 is literally +// `accum(i) += tmp_accum(i) * tCrScaleAViewAsC(i) * tCrScaleBViewAsC(i)`. +// CUTLASS associates the two scale multiplies left to right where the reference +// forms their product first (quant_utils.py:150-151); the difference is at most +// one f32 ULP per K-block and upstream's own gate admits it, comparing the two +// at rel_diff < 0.001 (test_block_fp8.py:194-200). We mirror the reference's +// association, because this op IS the reference port and the CUDA kernel that +// #1189 milestone M5 lands will be measured against it. +// +// SHAPES, with CEIL on every tiling, so a ragged final block is legal and must +// work (upstream asserts exactly this at fp8_utils.py:935-936): +// a_fp8 [M,K] i8, raw fp8-e4m3fn bytes +// a_scale [M, cdiv(K, block_k)] F32 +// b_fp8 [N,K] i8, raw fp8-e4m3fn bytes +// b_scale [cdiv(N, block_n), cdiv(K, block_k)] F32 +// out [M,N] f32 or bf16 +// The scales are f32 because upstream refuses any other dtype on this path +// (csrc/.../c3x/scaled_mm_helper.hpp:15-18) and the accumulator is f32. +// a_scale's K axis is a CEIL too, so this op accepts a K that vt::QuantFp8Group +// would refuse; that asymmetry is upstream's own (fp8_utils.py:930 uses cdiv +// where fp8_utils.py:596-599 demands divisibility). +// +// No bias: upstream refuses one outright on the blockwise path +// (scaled_mm_helper.hpp:54), so this mirrors a refusal rather than deferring a +// feature. +// +// CPU only. A CORRECTNESS REFERENCE, NOT A PERFORMANCE PATH — it is the +// numerical oracle #1189 milestone M5's CUTLASS kernel is measured against, and +// it makes no speed claim. M5 owns the CUDA arm. +void MatmulFp8BlockScaled(Queue& q, Tensor& out, const Tensor& a_fp8, const Tensor& a_scale, + const Tensor& b_fp8, const Tensor& b_scale, int block_n, + int block_k); + // RmsNormQuantFp8 (fused fp8 RMSNorm -> static per-tensor activation quant). One // HBM pass mirrors vLLM's Inductor `fused_add_rms_norm_static_fp8_quant` // (vllm/compilation/passes/fusion/rms_quant_fusion.py:124) — the RMSNorm producer diff --git a/src/vt/cpu/cpu_ops.cpp b/src/vt/cpu/cpu_ops.cpp index fcb0afd21..4768b52eb 100644 --- a/src/vt/cpu/cpu_ops.cpp +++ b/src/vt/cpu/cpu_ops.cpp @@ -628,6 +628,79 @@ void QuantFp8GroupKernel(Queue&, Tensor& out_fp8, Tensor& out_scale, const Tenso }); } +// --- Block-wise FP8 (VT-MATMUL-FP8-BLOCK-REF, #1189 M2). MatmulFp8BlockScaled +// CPU kernel: the 128x128 block-scaled fp8 GEMM, mirroring +// native_w8a8_block_matmul (tests/kernels/quant_utils.py:91-154). +// +// THE SCALES APPLY IN THE MAINLOOP, ONCE PER K-BLOCK. `part` is a SEPARATE f32 +// accumulator that is scaled and only then folded into `acc`. That separation is +// the whole point of the op: MatmulFp8CutlassKernel fifty lines below folds one +// scalar alpha AFTER the whole K reduction, which has exactly one degree of +// freedom per output element while this scheme has cdiv(K, block_k) of them. An +// epilogue-only application cannot express a per-K-block scale AT ALL. Collapsing +// `part` into `acc` here IS that epilogue form, and +// tests/vt/test_ops_matmul_fp8_block_cpu.cpp G4 is constructed so it cannot pass. +// +// The scale PRODUCT is formed first, `part * (a_s * b_s)`, mirroring +// `s = As_tiles[i] * Bs[j][i]` then `c += matmul(a, b.t()) * s` +// (quant_utils.py:150-151). The kernel that executes upstream is CUTLASS, not +// Triton (vllm/model_executor/kernels/linear/__init__.py:355-377, +// vllm/utils/deep_gemm.py:27-46), and it associates left to right -- +// `accum(i) += tmp_accum(i) * tCrScaleAViewAsC(i) * tCrScaleBViewAsC(i)`, cutlass +// 4.5.0 sm120_mma_tma_blockwise_scaling.hpp:714-717. The two differ by at most one +// f32 ULP per K-block, and upstream's own gate is what admits it: it compares the +// CUTLASS arm against THIS reference at rel_diff < 0.001 +// (test_block_fp8.py:194-200). We mirror the reference, because this op is the +// reference port. +// +// `col / block_n` indexes the b_scale ROW by OUTPUT COLUMN, mirroring +// `offs_bsn = offs_bn // group_n` (fp8_utils.py:823) and the reference's tiling +// (quant_utils.py:131-143). For a round N that agrees with a tile counter; for a +// ragged N it does not, which is why N=576 is in the ported grid. +// +// CEIL tiling and a short final K-tile: `k1 = min(k0 + block_k, k)` +// (quant_utils.py:131-141). Upstream's wrapper asserts the ceil shapes +// (fp8_utils.py:935-936), so a ragged block is legal rather than tolerated. +// +// A CORRECTNESS REFERENCE, NOT A PERFORMANCE PATH, in the same sense as +// MatmulFp8CutlassKernel below: a naive nest that makes the block-fp8 seam +// resolvable on a CPU queue so #1189 milestones M3 and M4 can be gated without a +// GPU. It makes no speed claim. M5 owns the CUDA arm and will be measured +// against this one. Parallel over ROWS; each output row is independent and the +// reduction order inside a row is fixed, so the result does not depend on the +// thread count. +void MatmulFp8BlockScaledKernel(Queue&, Tensor& out, const Tensor& a_fp8, const Tensor& a_scale, + const Tensor& b_fp8, const Tensor& b_scale, int block_n, + int block_k) { + const int64_t m = a_fp8.shape[0], k = a_fp8.shape[1], n = b_fp8.shape[0]; + const int64_t k_tiles = (k + block_k - 1) / block_k; + const auto* ap = a_fp8.Ptr(); + const auto* bp = b_fp8.Ptr(); + const auto* asp = a_scale.Ptr(); + const auto* bsp = b_scale.Ptr(); + ForRows(m, [&](int64_t r0, int64_t r1) { + std::vector arow(static_cast(k)); + for (int64_t i = r0; i < r1; ++i) { + // Decode the A row once and reuse it across N, as MatmulNvfp4Fp4Kernel does. + for (int64_t kk = 0; kk < k; ++kk) + arow[static_cast(kk)] = Fp8ToF32(ap[i * k + kk]); + for (int64_t col = 0; col < n; ++col) { + const int64_t nb = col / block_n; // fp8_utils.py:823 + float acc = 0.0F; + for (int64_t kt = 0; kt < k_tiles; ++kt) { + const int64_t k0 = kt * block_k; + const int64_t k1 = std::min(k0 + static_cast(block_k), k); + float part = 0.0F; // the MAINLOOP register, kept separate + for (int64_t kk = k0; kk < k1; ++kk) + part += arow[static_cast(kk)] * Fp8ToF32(bp[col * k + kk]); + acc += part * (asp[i * k_tiles + kt] * bsp[nb * k_tiles + kt]); + } + StoreF32(out, i * n + col, acc); + } + } + }); +} + // MatmulFp8Cutlass CPU kernel: out[m,n] = alpha * Sum_k f8val(a[m,k])*f8val(b[n,k]), // f32 accumulate, ONE folded alpha (= input_scale*weight_scale — our recorded // deviation from upstream's two epilogue scalars, see include/vt/ops.h). @@ -3256,6 +3329,9 @@ struct Registrar { reinterpret_cast(static_cast(&QuantFp8GroupKernel))); RegisterOp(OpId::kMatmulFp8Cutlass, DeviceType::kCPU, reinterpret_cast(static_cast(&MatmulFp8CutlassKernel))); + RegisterOp(OpId::kMatmulFp8BlockScaled, DeviceType::kCPU, + reinterpret_cast( + static_cast(&MatmulFp8BlockScaledKernel))); RegisterOp(OpId::kSiluAndMul, DeviceType::kCPU, reinterpret_cast(static_cast(&SiluAndMulKernel))); RegisterOp(OpId::kGeluAndMul, DeviceType::kCPU, diff --git a/src/vt/op_provider.cpp b/src/vt/op_provider.cpp index bd91bd174..24f6432d7 100644 --- a/src/vt/op_provider.cpp +++ b/src/vt/op_provider.cpp @@ -491,6 +491,8 @@ const char* OpNameImpl(OpId op) { return "AttentionRelPos"; case OpId::kQuantFp8Group: return "QuantFp8Group"; + case OpId::kMatmulFp8BlockScaled: + return "MatmulFp8BlockScaled"; case OpId::kCount: break; } diff --git a/src/vt/ops.cpp b/src/vt/ops.cpp index 8a8f6aaa2..26b8d78d1 100644 --- a/src/vt/ops.cpp +++ b/src/vt/ops.cpp @@ -673,6 +673,50 @@ void QuantFp8Group(Queue& q, Tensor& out_fp8, Tensor& out_scale, const Tensor& x reinterpret_cast(GetOp(OpId::kQuantFp8Group, q.device.type))( q, out_fp8, out_scale, x, group_size); } +void MatmulFp8BlockScaled(Queue& q, Tensor& out, const Tensor& a_fp8, const Tensor& a_scale, + const Tensor& b_fp8, const Tensor& b_scale, int block_n, + int block_k) { + VT_CHECK(out.rank == 2 && a_fp8.rank == 2 && a_scale.rank == 2 && b_fp8.rank == 2 && + b_scale.rank == 2, + "matmul_fp8_block_scaled: out/a_fp8/a_scale/b_fp8/b_scale must be rank-2"); + // Validated BEFORE either one divides anything: `x / 0` and `x % 0` are + // undefined behaviour, so a zero must refuse rather than trap. + VT_CHECK(block_n > 0 && block_k > 0, + "matmul_fp8_block_scaled: block_n and block_k must be positive"); + const int64_t m = a_fp8.shape[0], k = a_fp8.shape[1]; + const int64_t n = b_fp8.shape[0]; + // Upstream: `assert A.shape[-1] == B.shape[-1]` (quant_utils.py:111). + VT_CHECK(b_fp8.shape[1] == k, + "matmul_fp8_block_scaled: a_fp8 [M,K] and b_fp8 [N,K] must share K"); + VT_CHECK(out.shape[0] == m && out.shape[1] == n, "matmul_fp8_block_scaled: out must be [M,N]"); + VT_CHECK(a_fp8.dtype == DType::kI8 && b_fp8.dtype == DType::kI8, + "matmul_fp8_block_scaled: a_fp8/b_fp8 must be i8 (raw fp8-e4m3fn bytes)"); + // f32, not the model dtype: upstream refuses any other scale dtype on this + // path (csrc/.../w8a8/cutlass/c3x/scaled_mm_helper.hpp:15-18) and the + // accumulator these multiply into is f32. + VT_CHECK(a_scale.dtype == DType::kF32 && b_scale.dtype == DType::kF32, + "matmul_fp8_block_scaled: a_scale/b_scale must be f32"); + VT_CHECK(out.dtype == DType::kF32 || out.dtype == DType::kBF16, + "matmul_fp8_block_scaled: out must be f32 or bf16"); + // CEIL on every tiling, so a ragged final block is legal: upstream asserts + // `triton.cdiv(N, block_n) == Bs.shape[0]` and + // `triton.cdiv(K, block_k) == Bs.shape[1]` (fp8_utils.py:935-936), and + // `triton.cdiv(A.shape[-1], block_k) == As.shape[-1]` (fp8_utils.py:930). + const int64_t k_tiles = (k + block_k - 1) / block_k; + const int64_t n_tiles = (n + block_n - 1) / block_n; + VT_CHECK(a_scale.shape[0] == m && a_scale.shape[1] == k_tiles, + "matmul_fp8_block_scaled: a_scale must be [M, cdiv(K, block_k)]"); + VT_CHECK(b_scale.shape[0] == n_tiles && b_scale.shape[1] == k_tiles, + "matmul_fp8_block_scaled: b_scale must be [cdiv(N, block_n), cdiv(K, block_k)]"); + VT_CHECK(out.IsContiguous() && a_fp8.IsContiguous() && a_scale.IsContiguous() && + b_fp8.IsContiguous() && b_scale.IsContiguous(), + "matmul_fp8_block_scaled: contiguous tensors required"); + VT_CHECK(out.device == q.device && a_fp8.device == q.device && a_scale.device == q.device && + b_fp8.device == q.device && b_scale.device == q.device, + "matmul_fp8_block_scaled: device mismatch (out/a_fp8/a_scale/b_fp8/b_scale/queue)"); + reinterpret_cast(GetOp(OpId::kMatmulFp8BlockScaled, q.device.type))( + q, out, a_fp8, a_scale, b_fp8, b_scale, block_n, block_k); +} void RmsNormQuantFp8(Queue& q, Tensor& out_fp8, Tensor* out_bf16, const Tensor& x, const Tensor& weight, const RmsNormArgs& args, Tensor* residual, float input_scale) { diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index d97303bc8..ffddf5fc5 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1743,6 +1743,11 @@ vllm_cpp_add_test(test_ops_fp8_cpu vt/test_ops_fp8_cpu.cpp) # quant. CPU-gateable by construction; its CPU-vs-CUDA byte-identity arm (G6) is # CUDA-gated and reports PENDING rather than skipping where no device exists. vllm_cpp_add_test(test_ops_quant_fp8_group_cpu vt/test_ops_quant_fp8_group_cpu.cpp) +# VT-MATMUL-FP8-BLOCK-REF (#1189 M2): the 128x128 block-scaled fp8 GEMM, CPU +# reference arm — the numerical oracle every later block-fp8 kernel is measured +# against. Its G4 holds the constraint the row exists for: the scales apply in +# the MAINLOOP, once per K-block, which an epilogue-folded alpha cannot express. +vllm_cpp_add_test(test_ops_matmul_fp8_block_cpu vt/test_ops_matmul_fp8_block_cpu.cpp) # Opt-in arm: run the fp8 plan-cache byte-exact case with the cache ENABLED # (VT_FP8_PLAN_CACHE=1 -> first MatmulFp8CublasLt call builds the plan fresh, # later calls hit the cache). Proves the cached-plan GEMM is BYTE-identical to the diff --git a/tests/vt/test_ops_matmul_fp8_block_cpu.cpp b/tests/vt/test_ops_matmul_fp8_block_cpu.cpp new file mode 100644 index 000000000..480b6ecf7 --- /dev/null +++ b/tests/vt/test_ops_matmul_fp8_block_cpu.cpp @@ -0,0 +1,684 @@ +// vllm.cpp — the 128x128 block-scaled FP8 GEMM, CPU reference arm. +// +// VT-MATMUL-FP8-BLOCK-REF (.agents/specs/vt-matmul-fp8-block-ref.md), issue +// #1189 milestone M2. Pinned oracle: vLLM 5559679229bc961848b121ccdeaa8fa5d79bec98, +// asserted as the local checkout's HEAD before every anchor below was read. +// +// THE CONSTRAINT THIS FILE EXISTS TO HOLD. The scales apply in the GEMM +// MAINLOOP, once per K-block, into an f32 accumulator — not in the epilogue: +// +// accumulator = zeros(f32) +// for k_block: +// accumulator += dot(a_tile, b_tile) * a_s[m,k_block] * b_s[n_block,k_block] +// +// (vllm/model_executor/layers/quantization/utils/fp8_utils.py:826-836, the +// Triton arm spelling out the structure). Our per-tensor FP8 path folds ONE +// scalar alpha into the epilogue (`MatmulFp8CutlassKernel`, `alpha * acc` after +// the whole K reduction). An epilogue has exactly one degree of freedom per +// output element and the block scheme has cdiv(K, block_k) of them, so an +// epilogue-only application CANNOT EXPRESS a per-K-block scale at all. That is a +// correctness constraint, not an optimisation choice. G4 is the instrument for +// it and no single-alpha implementation can pass G4. +// +// WHICH IMPLEMENTATION ACTUALLY RUNS, asked because M1 found upstream's Triton +// source is not what executes and that the two arms DISAGREE in polarity. Here +// the chain lands differently — the executing kernel is CUTLASS and it AGREES +// with the reference: +// vllm/model_executor/kernels/linear/__init__.py:355-377 CUDA priority list +// vllm/utils/deep_gemm.py:27-46 DeepGEMM off for +// qwen3_5_text on family 120 +// .../scaled_mm/cutlass.py:312-326 ops.cutlass_scaled_mm(A, B.T, As, Bs.T) +// .../cutlass/scaled_mm_entry.cu:220-226 every sm>=120 -> sm120 +// .../cutlass/scaled_mm_c3x_sm120.cu:13-20 -> blockwise_sm120_fp8 +// .../c3x/scaled_mm_helper.hpp:15-18,39-55 blockwise branch: both +// scales f32 and 2-D, +// ceil_div shapes, no bias +// .../c3x/scaled_mm_blockwise_sm120_fp8_dispatch.cuh:56-58,218-235 +// ElementAccumulator=float, +// both scale pointers are +// MAINLOOP arguments +// cutlass 4.5.0 include/cutlass/gemm/collective/ +// sm120_mma_tma_blockwise_scaling.hpp:714-717 +// accum(i) += tmp_accum(i) * tCrScaleAViewAsC(i) * tCrScaleBViewAsC(i) +// which is native_w8a8_block_matmul's `c += matmul(a, b.t()) * s` +// (tests/kernels/quant_utils.py:145-151) with the two scale multiplies +// associated left to right instead of as a product. The two differ by at most +// one f32 ULP per K-block, and upstream's own gate is what admits it: it +// compares the CUTLASS arm against this very reference at rel_diff < 0.001 +// (tests/kernels/quantization/test_block_fp8.py:194-200). We mirror the +// reference's association, because this op IS the reference port. +// +// G1 the registration itself, and OpName totality. +// G2 the ported upstream case (test_block_fp8.py:123-153) over an adapted +// grid, against an INDEPENDENTLY WRITTEN double reference, carrying +// upstream's rel_diff < 0.001 verbatim AND a tighter f32 forward-error +// bound. Vacuity guard. +// G3 the ragged edges, called out separately because they are why the grid +// is what it is: N=576 (4*128+64), K=3884 (30*128+44), and both at once. +// G4 THE MAINLOOP CONSTRAINT, by exact hand-computed values plus a scale +// SWAP that no epilogue-folded alpha can distinguish. +// G5 the refusals, each by name. +// G6 the M1 seam: vt::QuantFp8Group -> vt::MatmulFp8BlockScaled end to end +// on a CPU queue. A COMPOSITION test, not a reachability claim: nothing +// in production dispatches either op yet and the spec's `## Owed` says so. +// +// THE REFERENCE IS INDEPENDENT, or the gate is a tautology. It accumulates in +// `double` with a DIFFERENT loop nest — upstream's, k-tile outer and n-block +// inner over whole tiles — while the kernel walks output elements with the +// k-tiles inside. Its fp8 decode is derived from the e4m3fn field layout +// (1 sign, 4 exponent bias 7, 3 mantissa, no infinities), not taken from any +// codec in src/. Reduction-order ULPs are exactly what the f32 forward-error +// bound admits; nothing here claims bit-exactness against a double. +// +// THE UPSTREAM GRID IS ADAPTED, and the adaptation is stated rather than +// hidden. Upstream runs itertools.product(M, N, K) with M=[1,7,8,83,4096], +// N=[128,512,576,7168,13824], K=[256,3884,4096,13824,16384] +// (test_block_fp8.py:48-55) — 125 combinations whose largest is 4096x13824x16384, +// about 9.3e11 multiply-accumulates. That is a GPU grid; a naive CPU reference +// nest, run twice, cannot execute it. Every axis VALUE below appears at least +// once and the ragged values appear separately and together; only the PAIRING +// is dropped. Parameters, dtypes, tolerance and failure criterion are preserved +// exactly. +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "vt/backend.h" +#include "vt/dtype.h" +#include "vt/ops.h" + +namespace { + +using vt::Device; +using vt::DeviceType; +using vt::DType; +using vt::Queue; +using vt::Tensor; + +Device Cpu() { return Device{DeviceType::kCPU, 0}; } + +Tensor MakeTensor(void* data, DType dt, Device dev, const std::vector& shape) { + Tensor t; + t.data = data; + t.dtype = dt; + t.device = dev; + t.rank = static_cast(shape.size()); + int64_t stride = 1; + for (int i = t.rank - 1; i >= 0; --i) { + t.shape[i] = shape[static_cast(i)]; + t.stride[i] = stride; + stride *= shape[static_cast(i)]; + } + return t; +} + +int64_t CDiv(int64_t a, int64_t b) { return (a + b - 1) / b; } + +// --- the independent e4m3fn decode ---------------------------------------- +// From the format: 1 sign bit, 4 exponent bits (bias 7), 3 mantissa bits, no +// infinities, 0x7F/0xFF the only NaN encodings ("fn"). Not taken from src/. +constexpr float kFp8Max = 448.0f; // quant_utils.py:27-35 finfo(e4m3fn).max +constexpr float kFp8Min = -448.0f; + +double E4m3Exact(unsigned exp_field, unsigned mant) { + if (exp_field == 0) return std::ldexp(static_cast(mant), -9); // mant / 512 + return std::ldexp(1.0 + static_cast(mant) / 8.0, static_cast(exp_field) - 7); +} + +const std::vector& DequantTable() { + static const std::vector table = [] { + std::vector v(256); + for (unsigned b = 0; b < 256; ++b) { + const double m = E4m3Exact((b >> 3) & 0xFu, b & 0x7u); + v[b] = (b & 0x80u) != 0 ? -m : m; + } + return v; + }(); + return table; +} + +// The 127 finite magnitudes 0x00..0x7E are MONOTONIC in the byte value, so the +// nearest-value encode is a binary search over their midpoints, with an exact +// midpoint resolving to the even mantissa (table index == byte value, so +// `lo & 1` IS the mantissa's low bit). +const std::vector& Magnitudes() { + static const std::vector table = [] { + std::vector v; + v.reserve(127); + for (unsigned e = 0; e <= 15; ++e) + for (unsigned m = 0; m <= 7; ++m) { + if (e == 15 && m == 7) continue; // the NaN encoding is not a value + v.push_back(E4m3Exact(e, m)); + } + return v; + }(); + return table; +} + +uint8_t EncodeRne(float r) { + const auto sign = static_cast(std::signbit(r) ? 0x80u : 0x00u); + const double a = std::fabs(static_cast(r)); + const std::vector& mag = Magnitudes(); + const auto it = std::lower_bound(mag.begin(), mag.end(), a); + size_t hi = static_cast(it - mag.begin()); + if (hi == 0) return sign; // a == 0 + if (hi >= mag.size()) return static_cast(sign | 0x7Eu); // a >= 448 + const size_t lo = hi - 1; + const double dlo = a - mag[lo], dhi = mag[hi] - a; + size_t pick = 0; + if (dlo < dhi) { + pick = lo; + } else if (dhi < dlo) { + pick = hi; + } else { + pick = (lo & 1u) == 0u ? lo : hi; + } + return static_cast(sign | static_cast(pick)); +} + +// --- operands --------------------------------------------------------------- +// Upstream builds each operand as +// A_fp32 = (rand(M,K) - 0.5) * 2 * fp8_max +// A_fp8 = A_fp32.clamp(min=fp8_min, max=fp8_max).to(fp8) +// (test_block_fp8.py:130-136), i.e. the fp8 bytes induced by a value that is +// UNIFORM on [-448, 448]. The transform is reproduced; the RNG is not torch's, +// which is a harness fact rather than a parameter — the values are random on +// both sides and nothing here depends on a particular draw. +// +// The encode is hoisted into a 2^20-entry table of uniformly spaced values +// because the grid below needs ~4e7 operand bytes and a binary search per byte, +// at the -O0 the CPU lane builds with, dominated the case's run time. The +// resulting byte distribution is upstream's to a value resolution of 8.5e-4, +// and a uniform draw lands below that magnitude with probability ~2e-6. +const std::vector& UniformByteTable() { + static const std::vector table = [] { + constexpr int kN = 1 << 20; + std::vector v(static_cast(kN)); + for (int i = 0; i < kN; ++i) { + const float u = (static_cast(i) + 0.5f) / static_cast(kN); // (0,1) + const float x = (u - 0.5f) * 2.0f * kFp8Max; + v[static_cast(i)] = EncodeRne(std::fmin(std::fmax(x, kFp8Min), kFp8Max)); + } + return v; + }(); + return table; +} + +std::vector RandomFp8(int64_t n, uint32_t seed) { + const std::vector& tab = UniformByteTable(); + std::mt19937 rng(seed); + std::uniform_int_distribution idx(0, static_cast(tab.size() - 1)); + std::vector out(static_cast(n)); + for (int64_t i = 0; i < n; ++i) out[static_cast(i)] = tab[idx(rng)]; + return out; +} + +// As = rand(M, k_tiles) * 1e-2 ; Bs = rand(n_tiles, k_tiles) * 1e-2 +// (test_block_fp8.py:143-144, factor_for_scale = 1e-2 at :127). +std::vector RandomScales(int64_t n, uint32_t seed) { + std::mt19937 rng(seed); + std::uniform_real_distribution u(0.0f, 1.0f); + std::vector out(static_cast(n)); + for (int64_t i = 0; i < n; ++i) out[static_cast(i)] = u(rng) * 1e-2f; + return out; +} + +// --- the independent reference --------------------------------------------- +// native_w8a8_block_matmul (tests/kernels/quant_utils.py:91-154), transcribed: +// k-tile outer, n-block inner, whole tiles, f32 semantics widened to double so +// that the kernel's own f32 rounding is what the bound measures. `abs_sum` +// carries the same nest and is the forward-error bound's magnitude term. +struct RefResult { + std::vector c; // [M,N] + std::vector abs_sum; // [M,N] sum over k-tiles of |partial| * |s| +}; + +RefResult RefBlockMatmul(const std::vector& ad, const std::vector& bd, + const std::vector& as, const std::vector& bs, int64_t m, + int64_t n, int64_t k, int64_t block_n, int64_t block_k) { + const int64_t n_tiles = CDiv(n, block_n); // quant_utils.py:123-124, CEIL + const int64_t k_tiles = CDiv(k, block_k); + RefResult r; + r.c.assign(static_cast(m * n), 0.0); + r.abs_sum.assign(static_cast(m * n), 0.0); + for (int64_t i = 0; i < k_tiles; ++i) { // quant_utils.py:145 + const int64_t k0 = i * block_k; + const int64_t k1 = std::min((i + 1) * block_k, k); + for (int64_t j = 0; j < n_tiles; ++j) { // quant_utils.py:146 + const int64_t n0 = j * block_n; + const int64_t n1 = std::min((j + 1) * block_n, n); + const double b_s = static_cast(bs[static_cast(j * k_tiles + i)]); + for (int64_t row = 0; row < m; ++row) { + const double a_s = static_cast(as[static_cast(row * k_tiles + i)]); + const double s = a_s * b_s; // quant_utils.py:150, PRODUCT first + for (int64_t col = n0; col < n1; ++col) { + double part = 0.0, part_abs = 0.0; + for (int64_t kk = k0; kk < k1; ++kk) { + const double p = ad[static_cast(row * k + kk)] * + bd[static_cast(col * k + kk)]; + part += p; + part_abs += std::fabs(p); + } + r.c[static_cast(row * n + col)] += part * s; // quant_utils.py:151 + r.abs_sum[static_cast(row * n + col)] += part_abs * std::fabs(s); + } + } + } + } + return r; +} + +std::vector Decode(const std::vector& bytes) { + const std::vector& tab = DequantTable(); + std::vector out(bytes.size()); + for (size_t i = 0; i < bytes.size(); ++i) out[i] = tab[bytes[i]]; + return out; +} + +double ReadOut(const std::vector& buf, DType dt, int64_t i) { + if (dt == DType::kF32) { + float v = 0.0f; + std::memcpy(&v, buf.data() + static_cast(i) * 4, sizeof(v)); + return static_cast(v); + } + uint16_t v = 0; + std::memcpy(&v, buf.data() + static_cast(i) * 2, sizeof(v)); + return static_cast(vt::BF16ToF32(v)); +} + +// One shape of the ported case. `out_dtype` is the store width; upstream runs +// bfloat16 (test_block_fp8.py:54) and the op admits f32 as well. +void RunPorted(int64_t m, int64_t n, int64_t k, DType out_dtype, uint32_t seed) { + constexpr int64_t kBlockN = 128, kBlockK = 128; // test_block_fp8.py:53 + const int64_t n_tiles = CDiv(n, kBlockN), k_tiles = CDiv(k, kBlockK); + CAPTURE(m); + CAPTURE(n); + CAPTURE(k); + + const std::vector a = RandomFp8(m * k, seed); + const std::vector b = RandomFp8(n * k, seed + 1u); + const std::vector as = RandomScales(m * k_tiles, seed + 2u); + const std::vector bs = RandomScales(n_tiles * k_tiles, seed + 3u); + + const int64_t out_elems = m * n; + std::vector out(static_cast(out_elems) * + (out_dtype == DType::kF32 ? 4u : 2u)); + + Queue q; + q.device = Cpu(); + Tensor ta = MakeTensor(const_cast(a.data()), DType::kI8, Cpu(), {m, k}); + Tensor tas = MakeTensor(const_cast(as.data()), DType::kF32, Cpu(), {m, k_tiles}); + Tensor tb = MakeTensor(const_cast(b.data()), DType::kI8, Cpu(), {n, k}); + Tensor tbs = + MakeTensor(const_cast(bs.data()), DType::kF32, Cpu(), {n_tiles, k_tiles}); + Tensor tout = MakeTensor(out.data(), out_dtype, Cpu(), {m, n}); + vt::MatmulFp8BlockScaled(q, tout, ta, tas, tb, tbs, static_cast(kBlockN), + static_cast(kBlockK)); + + const RefResult ref = RefBlockMatmul(Decode(a), Decode(b), as, bs, m, n, k, kBlockN, kBlockK); + + // (1) upstream's own criterion, verbatim (test_block_fp8.py:150-153): + // mean|out - ref| / mean|ref| < 0.001, with BOTH sides at the store width + // so that a bf16 store is common-mode exactly as it is upstream. + // (2) a tighter per-element f32 forward-error bound, because our arm IS the + // reference rather than a kernel measured against it: the f32 recursive-sum + // bound over K terms, x4 margin, plus the store's own half-ulp for bf16. + double sum_abs_diff = 0.0, sum_abs_ref = 0.0, worst_ratio = 0.0; + int64_t bad = 0, first_bad = -1, nonzero = 0; + const double eps_f32 = static_cast(std::numeric_limits::epsilon()); + for (int64_t i = 0; i < out_elems; ++i) { + const double got = ReadOut(out, out_dtype, i); + const double want = ref.c[static_cast(i)]; + const double want_stored = + out_dtype == DType::kF32 + ? static_cast(static_cast(want)) + : static_cast(vt::BF16ToF32(vt::F32ToBF16(static_cast(want)))); + sum_abs_diff += std::fabs(got - want_stored); + sum_abs_ref += std::fabs(want_stored); + const double tol = + 4.0 * static_cast(k) * eps_f32 * ref.abs_sum[static_cast(i)] + + (out_dtype == DType::kBF16 ? std::fabs(want) * 0.004 : 0.0); + const double diff = std::fabs(got - want); + if (tol > 0.0 && diff / tol > worst_ratio) worst_ratio = diff / tol; + if (!(diff <= tol)) { + if (bad == 0) { + first_bad = i; + CAPTURE(got); + CAPTURE(want); + CAPTURE(tol); + } + ++bad; + } + if (std::fabs(want) > 0.0) ++nonzero; + } + if (bad != 0) { + CAPTURE(bad); + CAPTURE(first_bad); + } + CAPTURE(worst_ratio); + CHECK(bad == 0); + const double rel_diff = sum_abs_diff / sum_abs_ref; + CAPTURE(rel_diff); + CHECK(rel_diff < 0.001); // test_block_fp8.py:153 + // VACUITY GUARD: an all-zero reference would make any implementation pass. + CHECK(nonzero == out_elems); +} + +} // namespace + +// =========================================================================== +// G1 — THE REGISTRATION ITSELF. +// +// A refusal test cannot stand in for this and that is measured, not argued: M1 +// deleted its CPU registration and every refusal case still passed, because the +// refusals live in the wrapper's validation and fire before dispatch. Only an +// OpRegistered assertion caught it (.agents/specs/vt-quant-fp8-group.md, +// "A refusal test cannot stand in for a registration test"). +TEST_CASE("G1: MatmulFp8BlockScaled is registered on the CPU backend and named") { + REQUIRE(vt::OpRegistered(vt::OpId::kMatmulFp8BlockScaled, DeviceType::kCPU)); + CHECK(std::string(vt::OpName(vt::OpId::kMatmulFp8BlockScaled)) == "MatmulFp8BlockScaled"); +} + +// =========================================================================== +// G2/G3 — the ported upstream case over the adapted grid. +// +// Every axis value of upstream's M/N/K lists appears at least once. The three +// ragged entries are G3 and are marked: N=576 is 4*128 + 64, K=3884 is +// 30*128 + 44, and one case carries both at once, which is where a floor +// division and a ceil division give different answers in different places. +TEST_CASE("G2: MatmulFp8BlockScaled matches native_w8a8_block_matmul over the ported grid") { + REQUIRE(vt::OpRegistered(vt::OpId::kMatmulFp8BlockScaled, DeviceType::kCPU)); + // M N K out seed + RunPorted(4096, 128, 256, DType::kBF16, 1101); // M=4096 + RunPorted(83, 128, 4096, DType::kBF16, 1201); // M=83, K=4096 + RunPorted(8, 7168, 256, DType::kBF16, 1301); // N=7168 + RunPorted(7, 13824, 256, DType::kBF16, 1401); // M=7, N=13824 + RunPorted(1, 512, 13824, DType::kBF16, 1501); // M=1 decode, K=13824 + RunPorted(1, 128, 16384, DType::kBF16, 1601); // K=16384 + // the f32 store arm, which upstream lists but leaves commented out at :54 + RunPorted(83, 128, 4096, DType::kF32, 1701); +} + +TEST_CASE("G3: MatmulFp8BlockScaled handles a ragged final N-block and K-block") { + REQUIRE(vt::OpRegistered(vt::OpId::kMatmulFp8BlockScaled, DeviceType::kCPU)); + RunPorted(8, 576, 4096, DType::kBF16, 2101); // N=576 = 4*128 + 64, round K + RunPorted(32, 576, 3884, DType::kBF16, 2201); // BOTH ragged: K=3884 = 30*128 + 44 + // Upstream's own dedicated ragged case, "weight.shape % 128 != 0, like in DSV3 + // kv_a_proj_with_mqa" (test_block_fp8.py:156-200). + RunPorted(32, 576, 7168, DType::kBF16, 2301); + RunPorted(32, 576, 3884, DType::kF32, 2401); +} + +// =========================================================================== +// G4 — THE MAINLOOP CONSTRAINT, and it is constructed so that no epilogue-folded +// alpha can pass it. +// +// Every value here is exact in bf16 and in f32, so these are equalities and not +// tolerances. Written as bare `==` rather than doctest's Approx: Approx carries +// a `scale` term and a strict `<`, so `epsilon(0.0)` compares `|a-b| < 0` and +// fails on values that ARE equal, while the default epsilon admits ~1.19e-5. Case A has two K-blocks whose partial products differ (128 and +// 256) and whose scales differ (0.25 and 0.5). SWAPPING the two K-block scales +// leaves every per-tensor summary of the scale tensor identical — same set, +// same sum, same product, same max — and changes the correct answer from 160 to +// 128. A kernel that reduces the scales to one epilogue alpha returns the same +// number for both and fails here. +TEST_CASE("G4: the scales apply per K-BLOCK in the mainloop, not once in the epilogue") { + REQUIRE(vt::OpRegistered(vt::OpId::kMatmulFp8BlockScaled, DeviceType::kCPU)); + constexpr uint8_t kOne = 0x38; // e4m3fn 1.0 = exponent field 7, mantissa 0 + constexpr uint8_t kTwo = 0x40; // e4m3fn 2.0 = exponent field 8, mantissa 0 + REQUIRE(DequantTable()[kOne] == 1.0); + REQUIRE(DequantTable()[kTwo] == 2.0); + + Queue q; + q.device = Cpu(); + + SUBCASE("A: swapping the two K-block scales changes the answer") { + constexpr int64_t kM = 1, kN = 1, kK = 256; // 2 K-blocks of 128, 1 N-block + std::vector a(static_cast(kM * kK), kOne); + std::vector b(static_cast(kN * kK), kOne); + for (int64_t i = 128; i < kK; ++i) b[static_cast(i)] = kTwo; // block 1: b = 2 + // partial(block 0) = 128 * 1 * 1 = 128 ; partial(block 1) = 128 * 1 * 2 = 256 + std::vector as = {1.0f, 1.0f}; + std::vector bs = {0.25f, 0.5f}; + + Tensor ta = MakeTensor(a.data(), DType::kI8, Cpu(), {kM, kK}); + Tensor tb = MakeTensor(b.data(), DType::kI8, Cpu(), {kN, kK}); + Tensor tas = MakeTensor(as.data(), DType::kF32, Cpu(), {kM, 2}); + Tensor tbs = MakeTensor(bs.data(), DType::kF32, Cpu(), {1, 2}); + + float got = 0.0f; + Tensor tout = MakeTensor(&got, DType::kF32, Cpu(), {kM, kN}); + vt::MatmulFp8BlockScaled(q, tout, ta, tas, tb, tbs, 128, 128); + // 128*0.25 + 256*0.5 = 32 + 128 + CHECK(static_cast(got) == 160.0); + + std::swap(bs[0], bs[1]); + float got_swapped = 0.0f; + Tensor tout2 = MakeTensor(&got_swapped, DType::kF32, Cpu(), {kM, kN}); + vt::MatmulFp8BlockScaled(q, tout2, ta, tas, tb, tbs, 128, 128); + // 128*0.5 + 256*0.25 = 64 + 64 + CHECK(static_cast(got_swapped) == 128.0); + // The load-bearing statement: an epilogue alpha cannot tell these apart. + CHECK(got != got_swapped); + // VACUITY GUARD. + CHECK(got != 0.0f); + CHECK(got_swapped != 0.0f); + } + + SUBCASE("B: the b-scale row is indexed by OUTPUT COLUMN / block_n, ragged block included") { + // N = 129 -> 2 N-blocks, the second of them one column wide. a = b = 1.0 + // everywhere, so every K-block partial is exactly 128 and the output is + // decided entirely by which scale pair the kernel picked. + constexpr int64_t kM = 2, kN = 129, kK = 256; + std::vector a(static_cast(kM * kK), kOne); + std::vector b(static_cast(kN * kK), kOne); + std::vector as = {1.0f, 1.0f, 2.0f, 2.0f}; // [2,2] + std::vector bs = {0.25f, 0.5f, 1.0f, 0.125f}; // [2,2] + + Tensor ta = MakeTensor(a.data(), DType::kI8, Cpu(), {kM, kK}); + Tensor tb = MakeTensor(b.data(), DType::kI8, Cpu(), {kN, kK}); + Tensor tas = MakeTensor(as.data(), DType::kF32, Cpu(), {kM, 2}); + Tensor tbs = MakeTensor(bs.data(), DType::kF32, Cpu(), {2, 2}); + + std::vector out(static_cast(kM * kN), 0.0f); + Tensor tout = MakeTensor(out.data(), DType::kF32, Cpu(), {kM, kN}); + vt::MatmulFp8BlockScaled(q, tout, ta, tas, tb, tbs, 128, 128); + + int64_t nonzero = 0; + for (int64_t row = 0; row < kM; ++row) { + const double a_row = row == 0 ? 1.0 : 2.0; + for (int64_t col = 0; col < kN; ++col) { + const size_t nb = col < 128 ? 0u : 1u; // OUTPUT COLUMN / block_n + const double want = 128.0 * (a_row * static_cast(bs[nb * 2 + 0])) + + 128.0 * (a_row * static_cast(bs[nb * 2 + 1])); + const double got = static_cast(out[static_cast(row * kN + col)]); + if (got != want) { + CAPTURE(row); + CAPTURE(col); + CAPTURE(got); + CAPTURE(want); + REQUIRE(got == want); + } + if (got != 0.0) ++nonzero; + } + } + // VACUITY GUARD. + CHECK(nonzero == kM * kN); + // Named explicitly so the ragged column is not merely covered by the loop: + // column 128 is the one-wide second N-block and must read bs row 1. + CHECK(static_cast(out[static_cast(128)]) == 144.0); + CHECK(static_cast(out[static_cast(kN + 128)]) == 288.0); + } +} + +// =========================================================================== +// G5 — the refusals, each by name. +TEST_CASE("G5: MatmulFp8BlockScaled refuses a malformed call by name") { + constexpr int64_t kM = 2, kN = 256, kK = 256; + constexpr int kBn = 128, kBk = 128; + const int64_t nt = 2, kt = 2; + std::vector a(static_cast(kM * kK), 0x38); + std::vector b(static_cast(kN * kK), 0x38); + std::vector as(static_cast(kM * kt), 1.0f); + std::vector bs(static_cast(nt * kt), 1.0f); + std::vector out(static_cast(kM * kN), 0.0f); + + Queue q; + q.device = Cpu(); + Tensor ta = MakeTensor(a.data(), DType::kI8, Cpu(), {kM, kK}); + Tensor tb = MakeTensor(b.data(), DType::kI8, Cpu(), {kN, kK}); + Tensor tas = MakeTensor(as.data(), DType::kF32, Cpu(), {kM, kt}); + Tensor tbs = MakeTensor(bs.data(), DType::kF32, Cpu(), {nt, kt}); + Tensor tout = MakeTensor(out.data(), DType::kF32, Cpu(), {kM, kN}); + // The well-formed call this case perturbs. + REQUIRE_NOTHROW(vt::MatmulFp8BlockScaled(q, tout, ta, tas, tb, tbs, kBn, kBk)); + + SUBCASE("a zero or negative block size, validated BEFORE it divides anything") { + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, tout, ta, tas, tb, tbs, 0, kBk), + std::runtime_error); + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, tout, ta, tas, tb, tbs, kBn, 0), + std::runtime_error); + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, tout, ta, tas, tb, tbs, -128, kBk), + std::runtime_error); + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, tout, ta, tas, tb, tbs, kBn, -128), + std::runtime_error); + } + SUBCASE("a b_scale sized by FLOOR instead of cdiv, on N and on K") { + // N = 200 -> cdiv = 2, floor = 1. A floor-tiled kernel would accept this. + constexpr int64_t kNr = 200; + std::vector br(static_cast(kNr * kK), 0x38); + std::vector bsf(static_cast(1 * kt), 1.0f); + std::vector outr(static_cast(kM * kNr), 0.0f); + Tensor tbr = MakeTensor(br.data(), DType::kI8, Cpu(), {kNr, kK}); + Tensor tbsf = MakeTensor(bsf.data(), DType::kF32, Cpu(), {1, kt}); + Tensor toutr = MakeTensor(outr.data(), DType::kF32, Cpu(), {kM, kNr}); + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, toutr, ta, tas, tbr, tbsf, kBn, kBk), + std::runtime_error); + // K = 200 -> cdiv = 2, floor = 1, on the K axis of b_scale. + constexpr int64_t kKr = 200; + std::vector ar2(static_cast(kM * kKr), 0x38); + std::vector br2(static_cast(kN * kKr), 0x38); + std::vector as2(static_cast(kM * 2), 1.0f); + std::vector bs2(static_cast(nt * 1), 1.0f); + Tensor ta2 = MakeTensor(ar2.data(), DType::kI8, Cpu(), {kM, kKr}); + Tensor tb2 = MakeTensor(br2.data(), DType::kI8, Cpu(), {kN, kKr}); + Tensor tas2 = MakeTensor(as2.data(), DType::kF32, Cpu(), {kM, 2}); + Tensor tbs2 = MakeTensor(bs2.data(), DType::kF32, Cpu(), {nt, 1}); + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, tout, ta2, tas2, tb2, tbs2, kBn, kBk), + std::runtime_error); + } + SUBCASE("an a_scale with the wrong group count or the wrong row count") { + std::vector narrow(static_cast(kM * 1), 1.0f); + Tensor tnarrow = MakeTensor(narrow.data(), DType::kF32, Cpu(), {kM, 1}); + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, tout, ta, tnarrow, tb, tbs, kBn, kBk), + std::runtime_error); + std::vector tall(static_cast((kM + 1) * kt), 1.0f); + Tensor ttall = MakeTensor(tall.data(), DType::kF32, Cpu(), {kM + 1, kt}); + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, tout, ta, ttall, tb, tbs, kBn, kBk), + std::runtime_error); + } + SUBCASE("a scale that is not f32") { + std::vector as16(static_cast(kM * kt), 0); + Tensor tas16 = MakeTensor(as16.data(), DType::kBF16, Cpu(), {kM, kt}); + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, tout, ta, tas16, tb, tbs, kBn, kBk), + std::runtime_error); + std::vector bs16(static_cast(nt * kt), 0); + Tensor tbs16 = MakeTensor(bs16.data(), DType::kBF16, Cpu(), {nt, kt}); + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, tout, ta, tas, tb, tbs16, kBn, kBk), + std::runtime_error); + } + SUBCASE("a packed operand that is not i8, and a mismatched inner dimension") { + std::vector a16(static_cast(kM * kK), 0); + Tensor ta16 = MakeTensor(a16.data(), DType::kBF16, Cpu(), {kM, kK}); + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, tout, ta16, tas, tb, tbs, kBn, kBk), + std::runtime_error); + std::vector bk(static_cast(kN * 128), 0x38); + Tensor tbk = MakeTensor(bk.data(), DType::kI8, Cpu(), {kN, 128}); + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, tout, ta, tas, tbk, tbs, kBn, kBk), + std::runtime_error); + } + SUBCASE("an out that is the wrong shape, the wrong rank, or an unsupported dtype") { + std::vector small(static_cast(kM * 8), 0.0f); + Tensor tsmall = MakeTensor(small.data(), DType::kF32, Cpu(), {kM, 8}); + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, tsmall, ta, tas, tb, tbs, kBn, kBk), + std::runtime_error); + Tensor t3d = MakeTensor(out.data(), DType::kF32, Cpu(), {1, kM, kN}); + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, t3d, ta, tas, tb, tbs, kBn, kBk), + std::runtime_error); + std::vector i8out(static_cast(kM * kN), 0); + Tensor ti8 = MakeTensor(i8out.data(), DType::kI8, Cpu(), {kM, kN}); + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, ti8, ta, tas, tb, tbs, kBn, kBk), + std::runtime_error); + } + SUBCASE("a non-contiguous operand") { + Tensor gappy = ta; + gappy.stride[0] = kK + 8; // a row gap: the K run is no longer contiguous + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, tout, gappy, tas, tb, tbs, kBn, kBk), + std::runtime_error); + } + SUBCASE("a device mismatch") { + Tensor foreign = ta; + foreign.device = Device{DeviceType::kCUDA, 0}; + CHECK_THROWS_AS(vt::MatmulFp8BlockScaled(q, tout, foreign, tas, tb, tbs, kBn, kBk), + std::runtime_error); + } +} + +// =========================================================================== +// G6 — the M1 seam, end to end on a CPU queue. +// +// vt::QuantFp8Group produces exactly the a_fp8/a_scale pair this GEMM consumes, +// and this is the composition a block-FP8 linear method will run. It is NOT a +// reachability claim: nothing in production dispatches either op at this merge +// commit, milestone M4 owns the wiring, and the spec's `## Owed` records it. +TEST_CASE("G6: QuantFp8Group feeds MatmulFp8BlockScaled on a CPU queue") { + REQUIRE(vt::OpRegistered(vt::OpId::kQuantFp8Group, DeviceType::kCPU)); + REQUIRE(vt::OpRegistered(vt::OpId::kMatmulFp8BlockScaled, DeviceType::kCPU)); + constexpr int64_t kM = 5, kN = 320, kK = 256, kG = 128; + const int64_t nt = CDiv(kN, 128), kt = kK / kG; + + std::mt19937 rng(9001); + std::uniform_real_distribution ux(-3.0f, 3.0f); + std::vector x(static_cast(kM * kK)); + for (auto& v : x) v = ux(rng); + + std::vector aq(static_cast(kM * kK)); + std::vector aqs(static_cast(kM * kt)); + Queue q; + q.device = Cpu(); + Tensor tx = MakeTensor(x.data(), DType::kF32, Cpu(), {kM, kK}); + Tensor taq = MakeTensor(aq.data(), DType::kI8, Cpu(), {kM, kK}); + Tensor taqs = MakeTensor(aqs.data(), DType::kF32, Cpu(), {kM, kt}); + vt::QuantFp8Group(q, taq, taqs, tx, static_cast(kG)); + + const std::vector b = RandomFp8(kN * kK, 7001); + const std::vector bs = RandomScales(nt * kt, 7002); + std::vector out(static_cast(kM * kN), 0); + Tensor tb = MakeTensor(const_cast(b.data()), DType::kI8, Cpu(), {kN, kK}); + Tensor tbs = MakeTensor(const_cast(bs.data()), DType::kF32, Cpu(), {nt, kt}); + Tensor tout = MakeTensor(out.data(), DType::kBF16, Cpu(), {kM, kN}); + vt::MatmulFp8BlockScaled(q, tout, taq, taqs, tb, tbs, 128, static_cast(kG)); + + const RefResult ref = RefBlockMatmul(Decode(aq), Decode(b), aqs, bs, kM, kN, kK, 128, kG); + const double eps_f32 = static_cast(std::numeric_limits::epsilon()); + int64_t bad = 0, nonzero = 0; + for (int64_t i = 0; i < kM * kN; ++i) { + const double got = static_cast(vt::BF16ToF32(out[static_cast(i)])); + const double want = ref.c[static_cast(i)]; + const double tol = + 4.0 * static_cast(kK) * eps_f32 * ref.abs_sum[static_cast(i)] + + std::fabs(want) * 0.004; + if (!(std::fabs(got - want) <= tol)) ++bad; + if (std::fabs(want) > 0.0) ++nonzero; + } + CHECK(bad == 0); + // VACUITY GUARD. + CHECK(nonzero == kM * kN); +}