From 4c2e2195915b1b43e283e7e7bd895b4ba3c8b8fc Mon Sep 17 00:00:00 2001 From: yhirose Date: Mon, 20 Jul 2026 08:39:34 -0400 Subject: [PATCH 01/12] ci: write the workflow comments in English MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The rest of the repo comments in English; ci.yml was the exception. No behaviour change — comments only. --- .github/workflows/ci.yml | 24 +++++++++++++----------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 96e95be..5265cfe 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -6,10 +6,10 @@ on: pull_request: jobs: - # 全OSで cpu / gpu / auto の3モードを実行する。 - # GitHub hosted runner のGPU事情: - # - macos-15 (arm64): paravirtual GPU 経由で Metal が動く → 実GPUテスト - # - ubuntu / windows: GPUなし → gpu/auto はCPUフォールバック経路の検証 + # All three modes — cpu / gpu / auto — on every OS. + # What each hosted runner actually offers: + # - macos-15 (arm64): Metal works through a paravirtual GPU -> a real GPU test + # - ubuntu / windows: no GPU -> gpu/auto exercise the CPU fallback path test: name: test / ${{ matrix.os }} runs-on: ${{ matrix.os }} @@ -26,9 +26,10 @@ jobs: - name: Test (cpu / gpu / auto) run: ctest --test-dir build -C Release --output-on-failure - # CUDAツールチェーンのビルド検証。hosted runner にNVIDIAドライバは無いので - # カーネルはコンパイル・リンクのみ。テスト実行は「ドライバ不在環境での - # dlopenフォールバック」の検証を兼ねる（これ自体が恒久的なテスト対象）。 + # Build verification for the CUDA toolchain. Hosted runners have no NVIDIA + # driver, so the kernels are only compiled and linked. Running the tests here + # doubles as coverage of the dlopen fallback on a driver-less machine, which + # is a permanent test target in its own right. cuda-build: name: cuda build+fallback / ${{ matrix.os }} runs-on: ${{ matrix.os }} @@ -59,10 +60,11 @@ jobs: - name: Test (CPU fallback — no driver on hosted runners) run: ctest --test-dir build -C Release --output-on-failure - # 実CUDA GPUでの実行テストは CI では行わない方針（有料GPU runnerは - # Org必須+従量課金、self-hostedは管理の手間があるため断念。2026-07）。 - # CUDAの数値回帰は NVIDIA 実機で `ctest --test-dir build` を手動実行する。 - # 将来方針を変える場合はランナーを登録し HAS_GPU_RUNNER=true で有効化。 + # Running the tests on a real CUDA GPU is deliberately out of CI (2026-07): + # paid GPU runners need an Org and bill per minute, and a self-hosted one is + # upkeep we chose not to take on. CUDA's numerical regression is instead a + # manual `ctest --test-dir build` on real NVIDIA hardware. + # To change that, register a runner and set HAS_GPU_RUNNER=true to enable this. cuda-gpu: if: vars.HAS_GPU_RUNNER == 'true' name: cuda gpu / self-hosted From 35f076aa671307358fac21f740a1ff05b42cd8ca Mon Sep 17 00:00:00 2001 From: yhirose Date: Mon, 20 Jul 2026 10:12:07 -0400 Subject: [PATCH 02/12] docs: stop tracking the internal design notes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit roadmap.md, architecture.md, performance-notes.md and the x86 runbook are development notes for my own use, not documentation for users of this library — but they sat in docs/, which is where a reader looks for the latter. Move them to docs/_internal/ on disk and ignore docs/_*, leaving docs/ free for user-facing documentation to be written into. Nothing is lost: the files stay where they were being read from, one directory down, and record.md was already ignored on the same reasoning. --- .gitignore | 4 + docs/architecture.md | 239 -------- docs/performance-notes.md | 968 --------------------------------- docs/roadmap.md | 599 -------------------- docs/x86-validation-runbook.md | 146 ----- 5 files changed, 4 insertions(+), 1952 deletions(-) delete mode 100644 docs/architecture.md delete mode 100644 docs/performance-notes.md delete mode 100644 docs/roadmap.md delete mode 100644 docs/x86-validation-runbook.md diff --git a/.gitignore b/.gitignore index 8c1c070..e6f62a7 100644 --- a/.gitignore +++ b/.gitignore @@ -7,5 +7,9 @@ compile_commands.json docs/record.md docs/record.ja.md +# Internal design notes — kept on disk under docs/_internal, not tracked. +# docs/ is being reserved for user-facing documentation. +docs/_* + # Downloaded MNIST cache for the convergence gate (bench/check_mnist.cpp) mnist-data/ diff --git a/docs/architecture.md b/docs/architecture.md deleted file mode 100644 index 87b87e8..0000000 --- a/docs/architecture.md +++ /dev/null @@ -1,239 +0,0 @@ -# cpp-tensorlib architecture - -Cross-platform F32 tensor library. Informed by silarray (the macOS-only -experiment this project succeeds) but designed fresh for three platforms; the -goal is a simple, elegant implementation that outperforms PyTorch on both CPU -and GPU. Serves as the matrix foundation for culebra's Tensor type (culebra -keeps autograd/VJP; this library owns the graph, fusion, and execution). - -## Target layer map - -``` - user code / culebra tl::array, tl::eval(), free functions - │ - array.h lazy graph + peephole fusion (M2) ← today: eager - │ │ - device dispatch (per-kernel-class thresholds, M3+) - │ - ┌───────────────────┼──────────────────────┐ - macOS: Accelerate+Metal │ own CPU microkernels │ own CUDA kernels - (M3) │ BLIS-style, AVX2/512/ │ PTX #embed, driver API - │ NEON + threadpool(M5) │ via dlopen (M6) - │ - reference backend detail:: strided loops in array.h — the correctness - (M1, permanent) oracle and universal fallback -``` - -Forward-looking plan (milestones, environment constraints, open decisions) -lives in [roadmap.md](roadmap.md). This file documents how the code works -today; performance methodology and gate results are in -[performance-notes.md](performance-notes.md). - -## Current state (M7–M9: dtype surface + decode kernels) - -The consumer-local-LLM decode path (batch≈1) is the focus (see roadmap "Target -scope"). Decode is memory-bandwidth-bound, so the levers are cutting weight -bytes (bf16/int4 storage) and fusing attention — not F32 GFLOP/s. - -- **Storage dtype (M7).** `types.h` has `dtype{f32,bf16}` + RNE converters; - `storage`/`array` carry a dtype and byte-sized allocation. bf16 is a - *weight-container* type: `array::to_bf16()`/`to_f32()`; the CUDA decode GEMV - consumes bf16 natively, and every other op transparently widens bf16 inputs to - an F32 copy in `eval_one`'s input funnel — so ref/cpu/accel/Metal kernels stay - F32-only, zero changes. Compute and results are always F32. -- **Decode GEMV (M7).** `tl_gemv_f32` / `tl_gemv_bf16v8` (8 cols/thread, 16-byte - loads, split-K for small N). `eval_one`'s dot case routes M=1 to the GEMV - (the 128×128 tile wastes 127 rows at M=1) — this alone tripled the F32 decode - baseline; bf16 weights add ~1.84× on the bandwidth-bound GEMV. -- **Fused attention (M9).** `array::attn_decode(q,K,V,scale)` (op_t::attn_dec) - → `tl_attn_decode_*`: flash-attention online-softmax in one pass (scores never - materialized), split-KV across `gridDim.y` to fill the SMs (KV-bandwidth-bound, - ~844 GB/s). CPU reference fallback; Metal returns false (widen/CPU path). -- **int4 storage dtype (M8).** `dtype::q4`: group-symmetric int4 weights, one - buffer holding packed [N,K] int4 + appended f32 scales, logical shape [K,N] so - `a.dot(Wq)` type-checks like f32 and rides the bf16 widen-fallback seam (decode - → `tl_gemv_q4`, else dequant to F32). `array::to_q4()`/`to_f32()`. Correct; - 1.79× vs bf16 on the GEMV (not yet bandwidth-bound). -- **CUDA buffer pool.** `gpu::alloc`/`release` recycle device buffers via a - size-keyed free list (like Metal's), so alloc/free-heavy workloads (training) - don't fragment the driver allocator (real-MNIST gate 7.8 → 5.8 s). -- **Result (RTX 3090, llama-7B decode shapes):** 3.5 → 61.7 tok/s over M7+M9+M8 - (F32 12.2, +fused-attn 26.1, +bf16 42.5, +int4 61.7). Full census in - performance-notes.md. Dev harness: `bench_llm` (f32/bf16 + attn; kept short so - WSL2 timing stays clean), plus direct kernel benches `bench_bf16_gemv` / - `bench_attn_decode` / `bench_q4_gemv` (each flushes — the authoritative - per-kernel numbers). - -## Current state (M5 first cut: own CPU GEMM) - -- `cpu.h` + `cpu_threadpool.h` — BLIS-style SGEMM for platforms without - Accelerate: cache-blocking (MC/KC/NC) + stride-aware packing (transposed - views feed in place, no materialization) around a register-blocked - microkernel, parallelized over the M dimension by a persistent thread - pool. Microkernel: NEON (8×8 tile, native on Apple Silicon / ARM Linux) - and a portable scalar fallback; AVX2/AVX-512 are the next drop-in behind - the same interface. Raw-pointer API (no array.h dependency, like metal.h). -- Dispatch (`eval_one` dot): `metal → accel::gemm → cpu::sgemm → ref::dot`. - `cpu::enabled_` gates it (default on; below Accelerate on Apple, primary - off-Apple; oracle tests force ref with it + `use_accelerate_` off). -- **Status: correct and complete, untuned.** Matches ref across edge/ - transposed/epilogue shapes on all modes. Scales to ~235 GFLOP/s at 1024³ - on M1 Pro NEON (vs Accelerate's AMX ~8× faster — expected; the own kernel - is the *off-Apple* path). Per-call overhead (pack-buffer allocation, - thread-pool sync even for one block) loses to the naive loop below ~256³; - removing it + the OpenBLAS-90% gate + AVX are the tuning pass (roadmap). - -## Earlier state (M4: culebra integration + tiny-tensor sprint) - -- **culebra integration** (M4): `culebra::TensorImpl` is an autograd tape - node wrapping a `tl::array`; culebra keeps VJPs, this library owns the - graph/fusion/execution. Vendored into culebra as a git submodule. - Feature-gated via `TL_RUNTIME_HOOKS`: an embedder's tensor-free binaries - reference no backend symbol (allocation/eval/barrier route through hooks - installed by `install_runtime_hooks()`). -- **Tiny-tensor sprint**: contiguous flat-loop fast paths (map/clone/add_), - eager-tiny evaluation in the graph builders (materialized ≤4096-elem - operands skip the node/eval machinery), thread-local topo-sort scratch - with visit-stamp marking, memoized const-wrap, tiny-dot/tiny-ew direct - paths in `eval_one`, pooled MTLBuffer contents caching. Per-op overhead - cut 2–3×; MNIST training 2× faster than the pre-M4 executor. -- **auto thresholds** finalized from a quiet-machine census (matmul 2e9, - elementwise 2e6, reduction 2e5 on M1 Pro) — see performance-notes.md. - -## Earlier state (M3b-3: STEEL SGEMM sprint) - -- STEEL kernels (`sgemm_steel_` 64×64, `sgemm_steel_32x64_` for M < 97) - ported from silarray/MLX: float2 fragment registers, precomputed loader - pointers, serpentine MMA, threadgroup swizzle, split edge loops - (interior / M-edge / N-edge / corner), affine-epilogue store shared - between interior and edge paths. NN operands only for now. -- Dispatch ladder in `metal::gemm`: STEEL band (NN, m ≥ 16, n ≥ 48, - k ≥ 16, BM by M band) → simple-tile family (64×32 / 32×32; handles - transposed views in place via strided loaders + float4 fast path). -- 2048³ sgemm: ~1300 → **~3200 GFLOP/s** over the sprint (loaded machine; - see performance-notes.md for the step-by-step history and census). - -## Earlier state (M3b-2: Metal SGEMM, softmax, reductions) - -- `metal_kernels.metal` adds: `sgemm_` (32×32×16 simdgroup-matrix tile, - trans_a/trans_b in-place loaders, affine epilogue fused in the store, - one edge-safe store path), `softmax_` (row-per-threadgroup, stable), - `row_sum_`/`row_max_` (last-axis reductions, epilogue applied). -- `eval_one` dispatch: dot → `metal_gemm` (maps transposed views to - trans flags, same classify as accel) → accel → ref; softmax and last-axis - sum/max → `metal_row` → ref. Other reduction axes stay on CPU. -- Full MLP forward (gemm→sigmoid→gemm→softmax) runs GPU end to end, - verified against the oracle across tile-boundary and transposed shapes. -- **Perf status**: correct and complete, not yet tuned — see - performance-notes.md "Known gap". The SGEMM optimization sprint (STEEL-class - kernel vs the PyTorch gate) is a later dedicated pass. - -## Earlier state (M3b stage 1: Metal elementwise + command-buffer batching) - -- `objc.h` — minimal objc_msgSend bridge (header-only, no .mm files). -- `metal.h` — device/queue/PSO context (lazy singleton), size-keyed - MTLBuffer pool, `#embed`'d MSL JIT-compiled on first dispatch, one - long-lived command buffer: dispatches accumulate, `flush()` = end + - commit + waitUntilCompleted at the end of every `graph::run`. - Non-Apple builds get inline stubs — callers have no `#ifdef`s. -- `storage.h` — when a Metal device exists every buffer is a pooled - shared-mode MTLBuffer (unified memory); heap otherwise. -- **CPU/GPU handoff**: every CPU-side read funnels through - `array::raw()`/`data()`, which call `metal::cpu_barrier()` (flush if - pending). One choke point makes arbitrarily mixed graphs safe — tested - gpu→cpu and cpu→gpu both ways. -- Kernels: elementwise binary/unary/affine, all applying the graph's affine - epilogue in the store (`fma(op, scale, offset)`), kernel bodies shared via - macro (silarray edge-tile lesson). Metal gemm/softmax/reductions are - stage 2. -- Dispatch order in `eval_one`: `metal (gpu mode) → accel → ref`. - `auto` joins a pending GPU pipeline but never starts one — measured - thresholds land with bench data. - -## Earlier state (M3a: Accelerate backend + dispatch seam) - -- `accel::` (in array.h; graduates to its own header with Metal) — vDSP / - vForce / CBLAS fast paths tried from `graph::eval_one`, `ref::` as - fallback. Every accel function returns nullopt/false when ineligible or on - non-Apple builds, so the evaluator carries no platform conditionals — this - is the dispatch seam Metal / CPU microkernels / CUDA plug into. -- GEMM maps transposed views onto CblasTrans (no materialization) and takes - the fused epilogue scale as alpha. `tl::use_accelerate_ = false` forces - the oracle; tests compare both paths on every dispatchable op class - including edge shapes. -- `bench/shared/speed/bench_main.cpp` — interleaved A/B medians (accel vs ref). - -## Earlier state (M2) - -- `types.h` — device mode switch (`use_cpu`/`use_gpu`/`use_auto`), - `gpu_available()`. -- `storage.h` — shared flat F32 buffer; the seam where device residency and - pooled recycling attach later. -- `array.h` — `tl::array`: shapes/strides, zero-copy views (transpose, - reshape, slice), numpy broadcast rules, lazy graph + fusion + evaluation: - - **Graph**: ops build `detail::node`s; output shapes are computed (and - shape errors thrown) at build time. `tl::eval(...)` / any data access - runs one iterative topo-sort pass over all roots. - - **Fusion**: every node carries an affine epilogue - (`result = op(...) * scale + offset`). Scalar chains compose into the - producing node — including binaries, `dot` and reductions — by copying - it with a composed epilogue. Fusion never mutates an existing node, so - shared intermediates stay valid (tested). - - **`ref::` backend**: naive strided kernels over materialized arrays. - Everything reduces to strides via `detail::for_each_index`, so one loop - shape serves every op/view/broadcast. The oracle and fallback; real - backends replace it per-op in `detail::graph::eval_one`. - -## Milestones - -Status summary; scope, environment needs, and approach are in -[roadmap.md](roadmap.md). - -| | Scope | Status | -|---|---|---| -| M1 | Core skeleton + reference CPU implementation + tests + CI | ✅ | -| M2 | Lazy graph, topo-sort eval, build-time peephole fusion | ✅ | -| M3a | Accelerate CPU backend + dispatch seam + bench harness | ✅ | -| M3b-1 | Metal foundation: context, buffer pool, elementwise, batching | ✅ | -| M3b-2 | Metal SGEMM (tiled), softmax, last-axis reductions; GPU MLP fwd | ✅ | -| M3b-3 | STEEL SGEMM sprint: 2048³ ~1300 → ~3200 GFLOP/s | ✅ | -| M4 | culebra integration; tiny-tensor per-op sprint; auto thresholds | ✅ | -| M5 | Own CPU backend: threadpool + BLIS-style microkernels (AVX2/AVX-512/NEON) | 🔨 NEON + AVX2 done (OpenBLAS gate met); AVX-512 deferred (needs hw) | -| M6 | Own CUDA backend: dlopen'd driver API, PTX `#embed`, SGEMM + split-K | 🔨 correct; ~0.82 of cuBLAS (prefill foundation — see roadmap "Target scope") | -| M7 | BF16/FP16 **storage** + decode GEMV dispatch | 🔨 storage dtype + bf16 GEMV done & integrated; bf16 **compute** GEMM (fine-tuning) deferred | -| M8 | Quantized (int4/int8) dequant-fused matmul (local-LLM inference) | 🔨 int4 GEMV done & integrated as a storage dtype (decode 61.7 tok/s); bandwidth-bound tuning + int8/GGUF-format remain | -| M9 | Fused attention (flash-attn style) + KV cache | 🔨 fused decode attention done & integrated (split-KV, ~19×); KV cache / GQA / causal mask remain | - -Milestones are now driven by the **target scope** (consumer local-LLM -inference/fine-tuning on 1–few consumer GPUs — RTX 5090 / A6000); see -`roadmap.md`. The CUDA bar is tokens/sec vs llama.cpp/exllamav2, not cuBLAS -GFLOP/s — full rationale in `performance-notes.md`. - -## Conventions - -- snake_case types/functions; trailing `_` for private members and internal - helpers; `k`-prefixed compile-time constants. -- Headers follow the cpp-httplib layout: declarations first, bodies below the - `// Implementation` divider. -- Comments state constraints and measured rationale, not narration. -- Kernel bodies must be shared templates across aligned/edge tile paths - (silarray bug-class lesson); every fused-epilogue kernel gets edge-shape - numerical tests. - -## Dependency policy - -culebra links everything statically (OS-guaranteed dynamic libraries are the -only exception), so this library carries **zero third-party dependencies** -(doctest is vendored, tests only): - -- macOS: Accelerate/Metal/MPS frameworks (OS-provided). -- CPU: own microkernels; no OpenBLAS (escape hatch: `-DUSE_OPENBLAS` may be - added later if a shape class demands it). -- CUDA: own kernels; `libcuda`/`nvcuda.dll` is dlopen'd at runtime (on WSL2: - `/usr/lib/wsl/lib/libcuda.so`), so binaries run without a driver and fall - back to CPU. No cudart, no cuBLAS, no CUTLASS (escape hatch: per-shape-band - CUTLASS instances if measurement demands). - -Size budget: < 1MB added per platform (< 2MB with CUDA kernels). AOT culebra -binaries that don't use Tensor must contain none of this — no unconditional -global initializers; lazy singletons only. diff --git a/docs/performance-notes.md b/docs/performance-notes.md deleted file mode 100644 index f1f2da2..0000000 --- a/docs/performance-notes.md +++ /dev/null @@ -1,968 +0,0 @@ -# Performance notes — how to measure, and what not to retry - -Started on day one (a silarray lesson: record refuted approaches so dead -ends are never re-walked). Every dispatch gate and design choice on a hot -path must carry the measurement that set it. - -## Performance gates - -The reference target is **PyTorch/libtorch** on the same hardware for the raw -GEMM foundation (silarray benched against MLX; this project is cross-platform, so -PyTorch is the portable BLAS-level bar). - -| Backend | Foundation gate (raw GEMM) | -|---|---| -| macOS (Accelerate + Metal) | ≥ silarray (which is ≈ MLX) and ≥ PyTorch MPS | -| Own CPU GEMM | ≥ 90% of OpenBLAS single/multi-thread, ≥ PyTorch CPU | -| Own CUDA GEMM | ≥ 90% of cuBLAS FP32 — **scoped: prefill/compute-bound only** | - -Where a gate cannot be met on a shape band after honest effort, record the -band and the data here, then consider the escape hatch (OpenBLAS / CUTLASS -for that band only) — do not ship a silent loss. - -**The CUDA GEMM gate is scoped, not absolute (2026-07-04).** The target scope is -consumer local-LLM inference/fine-tuning (roadmap.md "Target scope"), which is -decode-dominated and therefore *bandwidth-bound*, not compute-bound-GEMM-bound. -So the operative CUDA bar is **tokens/sec vs llama.cpp/exllamav2 on a real -quantized model**, and the cuBLAS-90%-on-square gate applies only to the -minority prefill/large-batch band — where it is **consciously accepted at ~0.82** -(split-K), not chased. See "CUDA GEMM regime analysis" below. - -## CUDA GEMM regime analysis — why the cuBLAS gate is de-prioritized (2026-07-04) - -The decision chain that reset the CUDA bar, recorded so it isn't relitigated: - -- **cuBLAS's moat is SASS, and it only matters in one regime.** Hand-written FP32 - SIMT SGEMM tops ~0.90–0.93 of cuBLAS on sm_86 (Boehm's public kernel is the - existence proof; ours is ~0.82 with split-K and untuned tiles). The last - ~7–18% is register-bank/instruction scheduling that ptxas won't reach from - CUDA C++ — cuBLAS's ~15-year SASS tuning. This gap is real but lives **only in - large compute-bound dense GEMM**. -- **The target workload rarely hits that regime.** Local-LLM interactive - **decode** (the dominant cost) is batch≈1 GEMV / dequant-GEMV — each weight - streamed once, FLOP-intensity ~2 → **100% VRAM-bandwidth-bound**. cuBLAS has - **no advantage** on bandwidth-bound work; a hand kernel that saturates memory - ties or beats it. **Attention** is flash-attention-style fused work cuBLAS - doesn't do. Only **prefill / larger batch** is compute-bound GEMM — a one-time, - minority cost for interactive use. -- **The SOTA precedent confirms it.** llama.cpp / exllamav2 / ggml / MLC-LLM are - the dominant local runtimes and hand-write their CUDA/Metal kernels with **no - cuBLAS on the hot paths**. Own-kernels is the *correct* approach for this scope, - not a compromise. And silarray's own data (bench/README) is the same shape: - it *tied* PyTorch-MPS on large square (0.96–1.0×) but *won 1.8×* on the - DL/small-batch shape (32×4096×768) — the win was never in the big-square regime. -- **Conclusion.** Keep the ~0.82 split-K SGEMM as the prefill foundation; do not - spend the CUTLASS escape hatch or a SASS-level sprint on the square gate. The - high-leverage CUDA work is the dtype/quant/attention kernels (roadmap M7–M9), - all bandwidth- or fusion-bound, all hand-written-friendly, all cuBLAS-irrelevant. - Benchmark them as **tokens/sec vs llama.cpp/exllamav2**, not GFLOP/s vs cuBLAS. - -## Measurement discipline (inherited from silarray, verified there) - -- GPU clock state swings absolute times ±30–60%: only same-run interleaved - A/B comparisons are meaningful. Check machine load before believing ratios. -- Per-kernel census → fix the single most expensive kernel → repeat. Don't - guess at bottlenecks. -- Warm up (~300 ms) before timing GPU work after CPU-only sections. -- Scale numeric tolerances with value magnitude (fp32 summation-order noise). -- Every dispatch gate is a measured value: comment the number and the bench - that produced it next to the gate. - -## Reference hardware - -| Machine | Notes | -|---|---| -| Apple M1 Pro | macOS backend; paravirtual GPU on CI is correctness-only — never bench on CI | -| RTX 3090 (sm_86), Ubuntu on WSL2 | CUDA backend; driver at `/usr/lib/wsl/lib/libcuda.so`. WSL2 adds submission latency vs native Linux — record which side of that line a number came from | - -## Current baselines (M1 Pro, interleaved medians) - -Sprint history for Metal SGEMM at 2048³ (all measured under load avg 5–9, -so treat as lower bounds): - -| Step | 2048³ | 1024³ | -|---|---|---| -| M3b-2 basic 32×32×16 tile, scalar loads | ~1300 GFLOP/s | ~1080 | -| + float4 fast loaders (host-gated) | ~2170 | ~1060 | -| + STEEL port (frag registers, swizzle, split edge loops) | **~3200** | **~1700–1960** | - -Reference points: Accelerate (AMX) sgemm ≈ 980 GFLOP/s at 1024³ and stays -~constant; M1 Pro GPU fp32 ceiling ≈ 4–5 TFLOP/s; MLX lands ≈ 3.5–4. The -STEEL numbers were taken on a loaded machine — re-measure quiet before -comparing against PyTorch-MPS for the gate. - -**Tile-config census (2026-07-02, loaded machine, simple-tile family):** -32×32 ≈ 64×32 ≈ 32×64 all plateau at ~2000–2200 GFLOP/s at 2048³ — the -bottleneck was the loader/pipeline structure, not tile size. 64×32 had the -best small-size behavior and is the non-STEEL fallback. - -**Remaining SGEMM work:** -- Transposed operands still use the simple-tile family (STEEL is NN-only - until SteelLoaderT is ported) — backward-pass matmuls (`xᵀ@g`, `g@Wᵀ`) - are the motivating shapes. -- Function-constant specialization (fc_mn_aligned / fc_k_aligned) not yet - ported; runtime branches instead. silarray measured conv2d FCs equivalent, - but gemm may differ — measure before porting. -- Narrow-N band (n < 48) and small shapes stay CPU-bound → auto-mode - thresholds are the fix, not kernel work. - -## MNIST training gate (M4 integration, 2026-07-02) - -Workload: culebra `benchmarks/mnist/train_bench.cul` — 784→30→10 MLP, -batch=10, hand-coded backprop, 6000 steps/epoch. The small-shape stress -test where per-op overhead dominates. Measured interleaved on a LOADED -machine (load avg ~9; medians were stable across 6 rounds regardless): - -| Stack | warm s/epoch | -|---|---| -| **culebra + tensorlib (cpu, JIT)** | **0.044** | -| PyTorch CPU (same recipe) | 0.066–0.075 | -| culebra + pre-M4 direct cblas | 0.076 | -| PyTorch MPS | 0.52 | -| culebra + tensorlib (gpu) | 3.8 | - -**Gate: PASSED for this workload** — 1.5–1.7x faster than PyTorch CPU, -1.7x faster than the pre-M4 implementation. All rows converge to -accuracy=0.9079 (numerics verified across devices). - -**The fix that flipped it** (M4 initially measured 1.6x SLOWER than -pre-M4): culebra's wrapper lifted `t * 2.0` and every VJP scalar to a -rank-0 tensor and called the tensor⊙tensor path — bypassing tl's scalar -overloads, so nothing fused and `W - d.dot(x)*lr` ran a ref:: broadcast -loop over the weight matrix (23520 elems × 6000 steps; `sample` showed -map_binary at 15x the samples of anything else). Routing -materialized rank-0 operands through the scalar overloads restored -affine/GEMM-epilogue fusion. Lesson: **fusion is only as good as the -embedder's entry points — audit the wrapper's lowering, not just the -kernels.** (Also: tl-level step microbench = 23 ms/epoch vs 44 through -culebra; the remaining ~2x is interp/JIT + tape-node overhead, known and -acceptable at batch=10.) - -GPU mode loses 87x on this shape (blocking flush per step × 6000; torch -MPS hides latency with async queues and still loses 12x) — that is what -the auto thresholds are for, not kernel work. - -## Tiny-tensor sprint (2026-07-02, microgpt gate) - -Workload: culebra microgpt (transformer, n_embd=16, 16–256 element -tensors) measured 2.6x slower than the pre-M4 direct executor — the -regime where per-op allocation dominates and the lazy-graph machinery is -pure overhead. Fixes, each measured on a 256-element op microbench -(loaded machine; ratios stable): - -| Change | Effect | -|---|---| -| Contiguous flat-loop fast paths in map_unary/map_binary/clone/add_ | add_ 1082→63 ns, clone 869→136 ns | -| graph::run: thread-local scratch + visit-stamp marking (no hash set) | part of ~2x lazy-op win | -| eval_one fast paths on raw node storage (elementwise + tiny dot ≤16K MNK, strided ok, epilogue folded) | dot 814→290 ns | -| as_node memoization + evaluated-node-as-const-cache on adoption | repeat-operand allocs gone | -| Eager-tiny in graph builders (materialized operands, ≤4096 elems, CPU-side): run flat loop, skip node/shell/eval entirely | lazy add 643→204 ns | -| culebra: broadcast check without shape construction; strides mirror → method | wrapper allocs −2/op | - -End state vs the pre-M4 executor (interleaved, per-step): - -| Workload | old/new | -|---|---| -| microgpt n_embd=16 (torture corner) | 0.55 — old wins 1.8x | -| microgpt n_embd=64 | 0.73–0.86 | -| **microgpt n_embd=128** | **1.4–1.55 — new wins** | -| **MNIST train (batch=10)** | **2.0 — new wins** (was 1.7 before sprint) | - -The remaining tiny-corner deficit is the tape+graph double bookkeeping -(~4–5 allocations/op that a direct executor doesn't pay). Grinding it -out needs node pooling / small-vector fields — invasive, and n_embd=16 -is a toy regime; the crossover sits between n_embd 64 and 128. Recorded -as accepted, not refuted: revisit only if a real workload lands in it. -Large-size and GPU numbers verified unchanged after the sprint -(checksums identical; 2048³ within load noise of its baseline). - -## auto-mode thresholds (M1 Pro, quiet-machine census 2026-07-03) - -`use_gpu_`-style gate in graph::metal_mode_: never break a pending GPU -pipeline; otherwise start GPU only above a per-kernel-class size -threshold (types.h auto_threshold_). Finalized from the crossover census -(`misc/census.cpp`, load ~4, interleaved medians — CPU=Accelerate vs -GPU-total=Metal single op + flush): - -| Class | Crossover (measured) | Threshold set | Old (provisional) | -|---|---|---|---| -| matmul | 1280³=2.1e9 tie, 1536³=3.6e9 GPU | **2e9** (~1260³) | 5e8 | -| elementwise | 1M cpu / 4M GPU | **2e6** | 4e6 | -| reduction (softmax) | 65536 tie / 262144 GPU | **2e5** | 8e6 | - -M1 Pro's AMX makes standalone GPU matmul only pay off past ~1300³; -reductions flip to GPU early (CPU's per-row loop is slow). Verified auto -tracks the faster backend at all four boundaries (mm 512³→cpu, -1536³→gpu; ew 256K→cpu, 8M→gpu) and stays CPU-fast on MNIST training -(auto 0.036 = cpu 0.039 s/epoch — tiny matmuls below threshold). - -**These are M1-Pro-specific.** They must be re-measured per Mac GPU/AMX -balance, and a full set is needed once the CUDA backend lands (M6). The -matmul value is the *standalone* crossover; a chained large-matmul -pipeline amortizes the flush and would justify a lower threshold — -unmeasured, so the conservative standalone value stands. - -**culebra default device:** left at `cpu` (not flipped to `auto`). -Rationale: the thresholds are machine-specific, current culebra -workloads (MNIST/microgpt-scale) sit below every threshold so auto would -behave like cpu anyway, and cpu is the predictable default; users opt -into GPU with `Tensor.use_gpu()`/`use_auto()`. Revisit if a large-tensor -workload appears or once CUDA lands. - -## Own CPU GEMM first cut (M5, M1 Pro NEON, 2026-07-03) - -Correctness-first scaffolding (BLIS blocking + packing + 8×8 NEON -microkernel + thread pool), untuned. Numbers on a loaded machine (load -5–11), so directional only: - -| | own (NEON, 8 thr) | vs Accelerate (AMX) | vs ref (naive) | -|---|---|---|---| -| 128³ | 24 GFLOP/s | 36× slower | 0.8× (loses) | -| 256³ | 69 GFLOP/s | 17× slower | 1.1× | -| 512³ | 156 GFLOP/s | 8× slower | 3.6× | -| 1024³ | 235 GFLOP/s | 8× slower | ~15× | - -Reading: scales correctly (compute-bound win grows with size), ~235 -GFLOP/s at 1024³ ≈ 55% of M1 Pro's NEON fp32 peak — decent for an untuned -kernel. **Loses to Accelerate by ~8× because AMX ≫ NEON — expected; the -own kernel is the off-Apple path, not a Mac competitor.** Loses to the -naive loop below ~256³ because per-call overhead (a `std::vector` pack -buffer per call/task, thread-pool sync even for one M-block) dominates -small sizes — the tuning pass's first targets. The real gate (OpenBLAS -90%) is measured off-Apple; deferred with the AVX kernels to the x86 box. - -Tuning targets, in likely priority: thread-local reusable pack buffers, -skip the pool for single-block work, MR/NR + MC/KC/NC sweep, software -prefetch, K-unroll in the microkernel. → Done; see the tuning pass below. - -**Gate caveat (2026-07-03):** `bench/cpu/speed/bench_cpu_gemm.cpp` links OpenBLAS -when present (`tensorlib_bench_cpu`). But Homebrew's OpenBLAS on M1 -appears to use generic ARMV8 microkernels, not M1-tuned ones — the -*untuned* own kernel already beats it ~1.8× at 1024³ (threads matched -via `openblas_set_num_threads`), which is implausible against a -properly-tuned OpenBLAS. So on a Mac, treat the OpenBLAS number as a -loose lower-bound only; the meaningful Mac target is **% of NEON fp32 -peak** (~410 GFLOP/s on M1 Pro → own is at ~54–64%, the real headroom). -The definitive OpenBLAS-90% gate is measured on the x86 box (M6) with a -properly-tuned OpenBLAS. The tuning sprint should track % of NEON peak -on the Mac and defer the OpenBLAS-90% verdict to x86. - -## Own CPU GEMM tuning pass (M5, M1 Pro NEON, 2026-07-03) - -The tuning sprint over the first cut. All changes oracle-verified (own -== ref, full suite green). Numbers below are the **quieter census** -(load ~4, one core taken by the `ccusage` statusline that couldn't be -stopped — so figures sit a few % under a truly idle box, but were -stable across runs); the sprint itself was done under heavy load -(4–28), where every change was accepted only on interleaved same-load -A/B. - -What landed, in order of impact: - -1. **Lane-indexed FMA microkernel** (biggest win). The first cut did 8 - `vdupq_n_f32` broadcasts per k-step; now A is loaded as two vectors - and consumed with `vfmaq_laneq_f32` (16 FMA + 4 loads per k-step, - zero dups; 20/32 NEON regs). Plus K-unroll ×4 and - `__builtin_prefetch` on both panels. Peak (2048³) reached ~354 - GFLOP/s through the dispatch path, ~365 kernel-isolated ≈ **86–89% - of NEON fp32 peak** (~410 GFLOP/s on M1 Pro), up from 54–64%. -2. **Thread-local reusable pack buffers** — the per-call/per-task - `std::vector` allocs are gone. Trap for later: a thread_local named - inside the worker lambda resolves to the *worker's* instance — the - first attempt read each worker's empty `bpack` and crashed; the fix - passes the caller's pointer in. -3. **MR-panel-granular M-parallelism** (was MC-block-granular): 256³ - has only 2 MC blocks → 2 of 8 threads used. Threads now take - contiguous MR-panel ranges, grouped into MC-row chunks to keep the - L2 blocking. -4. **KC 256→512** (MC=128, NC=2048 unchanged). **Confirmed in the quiet - census** (4 interleaved sweep rounds): KC=512 wins the compute-bound - regime — 512³ ~317 vs 256's ~304, 1024³ ~350 vs ~332, 2048³ ~363 vs - ~349 GFLOP/s — by ~1–5%. KC=256 leads only at 256³ (~231 vs ~216). - Since the gate is large-matmul throughput, KC=512 stands; the 256³ - give-up is accepted. - -Net effect through the dispatch path (`bench_cpu_gemm`, quiet census -vs the pre-tuning binary; medians of interleaved runs): - -| | before | after | vs ref | -|---|---|---|---| -| 128³ | 21 GF/s (0.8× ref — lost) | 65 GF/s | 2.4× | -| 256³ | 75 GF/s | 217 GF/s | 9.0× | -| 512³ | 131 GF/s | 310 GF/s | 12.3× | -| 1024³ | 208 GF/s | 340 GF/s | 14× | -| 2048³ | 246 GF/s | 348 GF/s | — | - -The small-size loss to the naive loop is fixed (128³ now 2.4× faster -than ref). **8×12 kernel: measured and rejected.** With the machine -quiet the 8×12 variant (24 FMA : 5 loads, 29/32 regs) was mixed and -within noise — a ~4% edge only at 2048³, a ~2% *loss* at 512³, a tie at -1024³ (KC=512, 4 interleaved rounds). Not worth its messier NR=12 edge -handling and near-full register file; 8×8 (20/32 regs) stays and keeps -headroom for future unroll/prefetch work. Block-size sweep harness: -`bench/cpu/speed/bench_cpu_sweep.cpp` — calls `cpu::sgemm` directly (compiles in -~1s) with `-DTL_CPU_MC/KC/NC` overrides (the constants in cpu.h are -`#ifndef`-guarded for this). - -Remaining for the M5 gate (needs the x86 box): AVX2/AVX-512 -microkernels, runtime CPUID dispatch, and the real OpenBLAS-90% -measurement. The on-Mac tuning is done. - -## AVX2 + runtime dispatch scaffold (M5, 2026-07-03) - -Written and compile-checked on the Mac; **not executed** (Rosetta stops -at SSE4.2 — the x86 box runs and tunes it). The compile-time `#ifdef` -microkernel pick became a runtime function pointer: `select_ukernel()` -returns NEON on ARM (fixed at compile time), and on x86 probes -`__builtin_cpu_supports("avx2"/"fma")` to pick the AVX2 kernel or the -scalar fallback. `sgemm` caches the pointer in a magic-static and calls -through it — the indirect call is amortized over the kc-loop. - -The AVX2 8×8 kernel is the x86 analogue of the NEON one: NR=8 = one -`__m256` per row → 8 accumulators (of 16 ymm), each k-step does one -B load broadcast-multiplied by the 8 A values (`_mm256_broadcast_ss` + -`_mm256_fmadd_ps`), K-unrolled ×4 with prefetch. It shares the MR=8/NR=8 -packing, so no second pack path. `__attribute__((target("avx2,fma")))` -lets it live in a baseline-x86 TU and be reached only after the CPUID -check — verified by cross-compiling (`-target x86_64-apple-macos -c`) -and disassembling: the body emits `vfmadd213/231ps` + `vbroadcastss`, -confirming the attribute took. Scalar + NEON kernels validated against a -naive triple-loop oracle across full and edge (mr<8, nr<8) tiles. - -**Deferred to the x86 box:** executing/numerically-verifying AVX2, -register-tuning the tile (8×8 vs a 6×16 that better fits 16 ymm), and an -AVX-512 kernel — a real one wants NR=16 (a second packing layout), so it -is a genuine x86-side task, not a Mac compile-check. The dispatch seam is -ready: `select_ukernel` just needs an `avx512f` branch. - -## Own CPU GEMM x86 census (M5, i7-12700KF, 2026-07-03) - -First real execution of the AVX2 path (per docs/x86-validation-runbook.md). -Box: 12th-gen i7-12700KF (Alder Lake, 8 P-core + 4 E-core = 20 logical), -Ubuntu 22.04 under **WSL2**, g++ 11.4, OpenBLAS 0.3.20. AVX2+FMA present, -no AVX-512. - -**AVX2 first-run correctness:** `check_cpu_ukernel` — scalar, avx2-8×8, and -the new avx2-6×16 all match the naive oracle across full + edge tiles; the -full array suite (cpu/gpu/auto) is green with the descriptor-driven driver. -This is the first time the AVX2 kernel has executed anywhere (Rosetta stopped -at SSE4.2 on the Mac). - -**8×8 → 6×16 register blocking (the tuning win).** The original 8×8 tile uses -only 8 of 16 ymm registers as accumulators; a 6×16 tile (6 rows × 2 ymm = 12 -accumulators, the canonical Haswell+ blocking, BLIS/OpenBLAS haswell) raises -in-register reuse. Each kernel now packs to its own tile via `ukernel_desc` -(the "second packing layout" the roadmap anticipated; AVX-512 will reuse the -NR=16 layout). Single-core (`taskset -c 0`, interleaved ×3, quiet), GFLOP/s: - -| size | 6×16 | 8×8 | 6×16 gain | -|---|---|---|---| -| 512³ | 131–142 | 114–119 | +17% | -| 1024³ | **140–148** | 116–122 | +18–20% | -| 2048³ | 128–138 | 110–118 | +15% | - -Single-P-core AVX2 peak ≈ 4.9 GHz × 2 FMA × 8 fp32 × 2 ≈ **157 GFLOP/s**, so -6×16 at ~143–148 (1024³) is **~91% of single-core peak** — matches/exceeds -the NEON kernel's 86–89%. The 8×8 sat at ~77%. 6×16 is now the x86 default; -8×8 is kept behind `-DTL_CPU_AVX2_8X8` for A/B. - -**KC re-sweep on 6×16 (single-core, interleaved ×3):** KC=512 is best/tied at -every size (1024³: 147/141/148; 2048³: 137/136/138); KC=256/384 slightly -lower, 768 ~= 512 at large N but down at 1024³. **KC=512 confirmed for x86**, -same as NEON — kept. MC=128/NC=2048 kept (not exhaustively re-swept; KC is the -dominant cache-blocking param and it held). - -**OpenBLAS-90% gate.** Full machine, threads matched (`openblas_set_num_threads` -= pool size), the stable compute-bound point is **2048³: own (6×16) ≈ 106% of -OpenBLAS** (~810–890 vs ~20–23 ms), interleaved — **gate MET**. 8×8 was -~95–103% there. Confined to 8 logical CPUs, own(6×16) 2048³ throughput -(~712 GF/s median) beat own(8×8) (~609) by ~17%, mirroring single-core. - -**Caveats (WSL2, recorded per methodology).** (1) WSL2 flattens the hybrid -topology — `lscpu` shows 10×2 with no MAXMHZ, so P/E cores can't be told apart -and `taskset` can't pin to P-cores; the Windows host scheduler places threads. -(2) `std::thread::hardware_concurrency()` ignores the taskset affinity mask -(always reports 20), so the pool always spawns 20 threads. (3) Mid-size -multi-thread (512³/1024³) is noisy and often < OpenBLAS — but single-core -1024³ is 91% of peak, so this is **not** a kernel deficit. It is the thread -pool statically splitting M across all 20 threads regardless of problem size, -over-parallelizing small GEMMs on the P/E+HT topology where OpenBLAS caps -threads. See the deferred "problem-size-aware thread count" item in -roadmap.md — best measured on native Linux, not this WSL2 box. - -## Own CUDA GEMM stage-2 (M6, RTX 3090 / WSL2, 2026-07-03) - -The SGEMM tuning sprint over the correctness-first one-output-per-thread -kernel. Bench: `bench/cuda/speed/bench_cuda_gemm.cpp` links cuBLAS + the CUDA runtime as a -**measurement reference only** (like OpenBLAS on the CPU side — never a library -dependency; cuda.h still dlopen's the driver). Timing is CUDA events on the null -stream, R launches per batch, median of interleaved own-vs-cuBLAS rounds. All -numbers are **device-side event timing on the RTX 3090 under WSL2** — WSL2 adds -submission latency, but events measure device-side, so the own/cuBLAS ratio is -WSL2-insensitive (wall-clock would not be). - -### The finding that reframed the sprint: managed memory was the first blocker - -The stage-1 backend used `cuMemAllocManaged` (chosen to reuse the Apple unified- -memory seam). Standing up the census immediately exposed an **~88× cliff** — -same cuBLAS SGEMM at 2048³, only the allocator differs: - -| allocation | cuBLAS 2048³ | -|---|---| -| `cudaMalloc` (device) | **~22,900 GFLOP/s** | -| `cudaMallocManaged` (managed) | **~260 GFLOP/s** | - -Root cause: **`cudaDevAttrConcurrentManagedAccess = 0` on WSL2** — the GPU -cannot fault-migrate managed pages, so it reads them over PCIe every access (the -~260 GF/s ≈ PCIe-bound). No managed escape exists on this box: `cuMemPrefetchAsync` -and every `cudaMemAdvise` variant return *"invalid device ordinal"*, and a -device-only-initialized managed buffer (host never touches it) still runs 256 -GF/s — migration is fundamentally off, not a first-touch issue. So the roadmap's -"first lever" (prefetch hints) is a dead end here. Because PyTorch uses -`cudaMalloc` device memory, *beating PyTorch requires device memory too* — no -kernel tuning matters on managed memory. This triggered the roadmap's pre- -authorized **device-buffer pivot**: the CUDA backend now keeps a persistent -host/device **mirror** per allocation (host malloc + `cuMemAlloc` device buffer, -dirty state keyed by the device pointer in cuda.h, so views sharing a storage -share one entry), with lazy H2D before a kernel reads a host-dirty buffer and -D2H before the CPU reads a device-dirty one (`array::raw()/data()` → -`gpu::sync_to_host`). Metal stays a strict no-op (genuinely unified). This -replaced managed everywhere; the array oracle + ctest cpu/gpu/auto validate the -coherence end-to-end. Recorded in [[cuda-wsl2-managed-memory-cliff]]. - -### SGEMM kernel ladder (interleaved own/cuBLAS, device memory) - -With operands on device memory, the naive kernel is the bottleneck. The tuned -`tl_sgemm_rb` fast path (NN, contiguous, K%8==0, N%4==0, 16B-aligned; everything -else — transpose/strides/odd shapes — falls to the correctness-first `tl_sgemm`): - -| step | 2048³ own/cuB | 4096³ own/cuB | -|---|---|---| -| naive one-output-per-thread (device mem) | ~0.09 | ~0.09 | -| 128×128×8 register block, 8×8 microtile, float4 loads | 0.68 | 0.75 | -| + warp-tiling (64×32 warp tile, WNITER=2) — conflict-free smem reads | 0.71 | 0.80 | -| **+ register-staged double buffer (current)** | **0.72–0.84** | **0.79–0.83** | - -The warp-tiled thread→output map makes every As/Bs fragment read a broadcast -within a warp, removing the 4-way bank conflict the plain `tid/16,tid%16` map -hits (its stride-8 B reads serialize). `ptxas`: 127 regs, no spills, 16KB smem -(double-buffered) → 2 blocks/SM (33% occupancy, register-limited by the 64 -accumulators — normal for GEMM). Absolute: ~18,500–19,400 GF/s at 4096³ (cuBLAS -~23,400). Correct throughout — max rel err vs cuBLAS ≤ 6e-5 at 2048³, exact at -4096³; full ctest green including the fast path (check_cuda's 64×48×40 NN case). - -**Gate status: not yet met.** Own ≈ **80% of cuBLAS**, and PyTorch-FP32 (TF32 -disabled; it bundles cuBLAS 12.8 and measured ~25,900 GF/s at 4096³, *faster* -than the linked cuBLAS 13.1 ~23,400) is the tougher bar at ~0.75. The remaining -levers to close 80%→90%+ (in likely order): larger block tiles (128×256 for more -C-reuse per load — needs register/spill management), split-K for the mid-size -wave-quantization tail (2048³ is only ~2 waves at 128² tiles — the ratio dips vs -4096³'s ~6 waves), and finer instruction scheduling. Per the roadmap escape -hatch, a per-shape-band CUTLASS instance is the fallback if a band can't reach -the gate after that. - -### Split-K census (ladder ②, 2026-07-04) - -The census signature — own/cuB rising monotonically with size (1024³ 0.63, -2048³ 0.74, 4096³ 0.82) — is the classic occupancy/wave-quantization tell: 128² -tiles make ceil(m/128)·ceil(n/128) blocks, and the RTX 3090's 82 SMs hold 2 -blocks each = 164 concurrent. So 1024³ = 64 blocks fills only ~39% of the SMs -(0.39 wave), 2048³ = 256 blocks = 1.56 waves (half-empty tail), 4096³ = 1024 -blocks = 6.24 waves (tail negligible → the 0.82 there is pure kernel efficiency, -not a fill problem). Split-K partitions the K axis into S z-slices (gridDim.z=S) -so S× more blocks run at once; crucially it *partitions* K rather than -replicating it, so A/B global traffic is unchanged — the only added cost is C -written S× via `atomicAdd` into a pre-zeroed buffer (hence gated to identity -scale/offset; fused-affine GEMM keeps the single-slice store). Forced-S census -(interleaved own/cuB, maxrel ≤1.2e-4 — fp32 atomic-order noise, acceptable): - -| size | S=1 | S=2 | S=4 | S=8 | -|---|---|---|---|---| -| 1024³ | 0.63 | **0.75** | 0.75 | 0.73 | -| 2048³ | 0.74 | **0.76** | 0.76 | 0.68 | -| 4096³ | **0.82** | 0.82 | 0.79 | 0.75 | - -**S=2 is the robust optimum** for the underfilled mid sizes (1024³ +19%, 2048³ -+3%). S≥4 regresses: the S× C-atomic traffic + the shorter per-block K (fixed -smem-pipeline fill/drain amortized over fewer slabs) overtake the occupancy gain; -by S=8 both mid sizes drop below their S=1. 4096³ is already 6 waves so any split -is pure overhead → S=1 (and S=1 keeps C exact, maxrel 0.0). Kernel change is a -single `ksplit` param + `blockIdx.z` K-range + a `gridDim.z>1 ? atomicAdd : store` -epilogue on the *same* `tl_sgemm_rb` body (no copy-paste divergence; still 127 -regs / 2 blocks/SM). Auto policy: S=2 when base<512 blocks and K≥512 (each half -≥256 K, enough to amortize the pipeline), else S=1; `TL_SPLITK` forces S for the -census. Post-split state: 1024³ ~0.72–0.82, 2048³ 0.76, 4096³ 0.82–0.83. - -**Gate still open** at 2048³/4096³ (~0.76/0.82). Split-K fixed the small-size -*fill* deficit but the per-SM *efficiency* ceiling (~0.82, visible at 4096³ where -fill is a non-issue) is untouched — that is the remaining ~10-18% and needs the -kernel-efficiency lever (ladder ③: Ampere `cp.async` staged pipeline — removes -the register-staged global loads, frees registers, enables >2-stage overlap), -then the per-band CUTLASS escape hatch if it still misses. - -**Refuted (do not retry without new evidence):** - -| Approach | Verdict | Data | -|---|---|---| -| BK 8→16 (BK=16 slab, 2×float4 loads) | **Rejected** | Helps 1024³ (0.62→0.70) but *regresses* the gate sizes: 2048³ 0.72→0.62, 4096³ 0.83→0.73. The 32KB smem + 130 regs hurt where large-matmul throughput matters. | -| Register double-buffer *before* warp-tiling | **No effect** | On the conflict-bound (pre-warptile) kernel, hiding global latency did nothing (2048³ 0.68→0.67) — the kernel was smem-bound, not latency-bound. Double buffering only paid off *after* warp-tiling made the smem reads conflict-free. | -| **128×256 large tile (ladder ①)** — 256-thread block, 8×16 microtile, WNITER=4 (`tl_sgemm_rb2`), +33% arithmetic intensity (64→85) | **Rejected** | Loses at *every* size (interleaved own/cuB): 1024³ 0.69→0.41, 2048³ 0.75→0.68, 4096³ 0.81→0.75. 128 accumulators/thread → 203 regs (no spills) but **1 block/SM = 16.7% occupancy** (rb keeps 2 blocks/SM = 33%). It loses even at 4096³ (6 waves, machine full, wave-quant irrelevant), so the arithmetic-intensity gain (−25% DRAM traffic) is dwarfed by the occupancy/latency-hiding loss — same class as the BK=16 rejection. 1024³ collapses because 128×256 halves the block count (64→32) so only ~39% of the 82 SMs get work. Big spatial tiles are ceiling-capped at rb on sm_86 for this kernel: any 128×256 variant is ≥1 block/SM. Pivoted to split-K (below) for the small-size deficit that census confirms is the real problem. | -| **cp.async staged pipeline (ladder ③)** — Ampere `cp.async` global→shared, 3/4-stage smem pipeline (`tl_sgemm_cp`), replacing the register-staged loads | **Rejected** | Loses where it was meant to win: 4096³ 0.83→0.72–0.74, 2048³ 0.75–0.80→0.62–0.69 (interleaved vs rb S=1); ties at 1024³. Correct throughout (maxrel = rb, 4096³ exact). Root cause is a **structural mismatch**: A must land *transposed* in smem (As[k][m]) for the regM fragment read to stay a contiguous float4, but cp.async copies a contiguous global chunk to a contiguous smem chunk — it can't transpose. So A is forced to 4×**4-byte** cp.async scatter (one per k of the thread's float4), and 4-byte cp.async is the inefficient granularity (16-byte `.cg` is the fast path); that penalty + the extra `__syncthreads` (2/slab vs rb's 1) sink it. 4 stages don't rescue it (4096³ still 0.73). Efficient cp.async for transposed-A SIMT needs `ldmatrix` (tensor-core-only) — out of scope for the FP32 SIMT kernel. rb's float4-ldg-then-scatter stays the load path; regs dropped 127→107 with cp.async but occupancy stayed 2 blocks/SM (64 accumulators, not loads, are the register floor), so even the freed registers bought nothing. | - -## M7–M9 decode kernels (RTX 3090 / WSL2, 2026-07-04) - -The consumer-local-LLM decode path (batch≈1, bandwidth-bound). Metric is -tokens/sec on llama-7B decoder shapes (d=4096, ffn=11008, 32 heads × head_dim -128, 32 layers, vocab 32000), measured via `bench_llm`; kernel-level numbers via -the direct `bench_bf16_gemv` / `bench_attn_decode` / `bench_q4_gemv` benches -(cuda:: API, cuda::flush + steady_clock, no cudart). RTX 3090 memory-bandwidth -peak ≈ 936 GB/s — the ceiling that matters here, not GFLOP/s. - -**Whole-decode tok/s arc (F32 baseline → M7+M9):** - -| stage | tok/s | lever | -|---|---|---| -| tile-GEMM (session start) | 3.5 | M=1 dot ran the 128² tile → 127/128 rows wasted | -| + GEMV dispatch (F32) | 12.2 | route M=1 dots (incl. attention's) to `tl_gemv_f32` | -| + bf16 weights | 14.8 | `tl_gemv_bf16v8`, halved weight bytes | -| + fused attention (F32) | 26.1 | `attn_decode` — attention was 72% of token time | -| + bf16 weights + fused attn | 42.5 | both | -| + int4 weights + fused attn | **61.7** | `dtype::q4` GEMV (M8) | - -**bf16 decode GEMV** (`bench_bf16_gemv`, layer weight shapes): f32 kernels run -at 850–935 GB/s (≈ peak — bandwidth-bound). Scalar 2-byte bf16 loads under-use -that (layer sum 1.62× vs f32); `tl_gemv_bf16v8` (8 cols/thread, one 16-byte uint4 -load, `n%8==0`) closes it to **1.84×** (ffn_gate_up/lm_head at the 2× ceiling, -~890 GB/s). Small-N shapes (N=4096) lag — split-K atomic + launch overhead is a -bigger fraction of their tiny work. - -**Fused decode attention** (`bench_attn_decode`, H=32, ctx=2048): the unfused -3-op array path (q·Kᵀ, softmax, ·V, per head × layer ≈ 1024×) is -launch/materialize-bound, ~48 ms/model. `tl_attn_decode_f32` (one block/head, -online softmax) → 11.7 ms (4×) but only ~180 GB/s: grid = H = 32 blocks fills 3% -of the 82 SMs, occupancy-bound. **Split-KV** (partition ctx over `gridDim.y`, -partial states → combine pass) → **2.54 ms, 844 GB/s (~19×)**, bandwidth-bound. -Numerically exact throughout (maxrel ~1e-7). Array-integrated at 0.110 ms/layer -(13.5× the unfused path). - -**KV cache + GQA/MQA** (`bench_attn_decode`, 2026-07-08): decode-side, kernel + -direct-bench first (STEEL). The persistent cache is a `cuda::kv_cache` object (K,V -`[n_kv_heads, max_ctx, D]` device buffers + `pos`), *outside* the lazy graph — -decode is stateful/inference-only, a poor fit for the immutable node model. The -linchpin kernel change is a **`kv_stride`** param (head stride = `max_ctx*D`, -distinct from the valid length `ctx`), so `tl_attn_decode_split` reads a persistent -max_ctx allocation as its prefix `[0,pos)`; `kv_stride==ctx*D` degenerates to the -pre-cache whole-buffer path exactly (the array `attn_decode` still uses it — MHA -numbers above unchanged, re-verified maxrel ~1e-7). A **`group`** param -(`n_q_heads/n_kv_heads`) adds GQA/MQA: query head `h` reads kv head `h/group`. -`tl_kv_append` scatters one token's k,v `[n_kv_heads,D]` into the cache at `pos`. -Verified by growing the cache token-by-token to 2048 and comparing a from-scratch -CPU reference (prefix + GQA mapping) at 11 checkpoints across the split-KV boundary -(ctx≥256) and its chunk rounding — **maxrel 1.6e-7**. GQA (32 q / 8 kv, group 4) at -ctx=2048 runs at **0.044 ms/layer vs MHA's 0.081** (~1.85× faster wall-time) and 4× -less KV cache VRAM. The reported "unique KV GB/s" (383) understates real efficiency: -with group=4, four q-head blocks each stream the *same* kv head, so physical DRAM -traffic sits between the unique (HKV) and redundant (HQ) bounds, partly served by -L2 — the next tuning lever is a per-group block that reads each kv head once (not -once per q head in the group). - -**Causal prefill** (`bench_attn_decode`, 2026-07-08): the prompt regime, baseline -kernel + direct-bench. `tl_attn_prefill_f32` processes all T query positions at -once — query `p` attends keys `0..p` (**causal mask** = loop cap `i<=p`), one block -per (head, query pos) so grid = H_q×T fills the SMs with no split-KV. It reuses the -decode online-softmax verbatim (only the loop cap + query-row index change), so no -T×T scores are materialized. `tl_kv_fill` bulk-copies a prompt's k,v `[H_kv,T,D]` -into the cache `[0,T)`; `kv_cache::prefill` chains fill → attend and sets `pos=T`. -Verified vs a from-scratch causal CPU reference (T=512): **prefill maxrel 1.6e-7**, -and a follow-on decode step at pos=512 **maxrel 8.5e-8** (the prefill→decode -handoff — prefill-filled rows + one appended token — is exact). ~335 Ktok/s (T=512, -one layer). **Deliberately un-tiled**: each query re-streams its keys from DRAM -(O(T²) traffic → bandwidth-bound, not compute-optimal). The tuning pass is a -query×key tiled flash-attention (shared-mem K/V reuse across a query block), gated -on prefill proving a bottleneck — per scope, prefill is a minority one-time cost, -so the correct baseline is the right stopping point for now. T>65535 needs -gridDim.y chunking. - -**int4 dequant GEMV** (`bench_q4_gemv`, group=32 symmetric int4, [N,K] layout): -correct vs a host dequant reference (maxrel ~1e-5), 0.625 bytes/weight (0.5 -packed + 4/32 scale) vs bf16's 2.0. **1.79× vs bf16** (layer sum 0.436 vs 0.78 -ms). *Not yet bandwidth-bound* (best ~565/936 GB/s): the global-a kernel re-reads -the activation from L2 per output warp; shared-a staging (`tl_gemv_q4s`, gated to -K·4 ≤ 24 KB — larger K collapses occupancy) helps but the strided a_sh reads -bank-conflict. Integrated as `dtype::q4` (storage-dtype path; `a.dot(Wq)` rides -the bf16 widen-fallback seam) → **decode 61.7 tok/s** (vs bf16 42.5). - -> **Update 2026-07-17 (supersedes the global-a/shared-a split above):** the whole -> global-a-vs-shared-a distinction is gone. `tl_gemv_q4` was rewritten to -> one-block-per-row + K-adaptive block size (llama.cpp `mul_mat_vec_f` strategy; -> commits 12f423a / a426bb0), and `tl_gemv_q4s` was deleted. Under one-row-per- -> block every thread reads a **disjoint** K-slice, so there is no per-warp -> activation re-read for shared staging to amortize — the a-traffic problem this -> paragraph chases was an artifact of the old 8-rows/block layout, not intrinsic. -> On the RTX 3090 at real Qwen decode this fix is net-neutral for q4 (its win was -> mostly absorbed by the same fix applied to the bf16-row path); see the M9 -> decode-speed history. The analysis above is the 2026-07-04 point-in-time record. - -### WSL2 sysmem-fallback cliff — root cause of the long-run bench slowdown (2026-07-04) - -`bench_llm` (array `eval()`) is clean at the M9 length but, once an int4 section -makes the run longer, later measurements go **5–10× slow and stay slow**. Ran it -to ground with a fixed small "probe" op timed as unrelated large shapes are -processed: - -- **Not** op-count/thermal/clock: the same op repeated 300× never degrades, and - during a degraded 400× probe the SM clock is pinned at **1710 MHz with throttle - reasons 0x0** — full clock, no throttle. The slow op is compute-idle-bound. -- **Not** the buffer pool: reproduces identically with pooling disabled (always - cuMemFree). **Not** allocation fragmentation, **not** current memory level - (degrades at 2.3 GB used / 24 GB; a single 1 GB alloc+free never triggers it), - **not** the concurrent working set (freeing each weight before the next still - degrades). pow2 size-binning (more reuse, memory capped at 2.4 GB) did **not** - fix it either. -- **It is cumulative and persistent:** degrades only after ~a dozen distinct - large allocations (lm_head *first* is clean; lm_head *last* degrades), and once - triggered it never recovers (400 iterations flat at 0.63 ms vs 0.088 ms clean). - -Diagnosis: **WSL2/WDDM sysmem fallback.** WDDM reclaims cuMemFree lazily, so the -process's committed high-water climbs monotonically; once it crosses the WDDM -working-set budget (~2 GB effective here, far below the 24 GB VRAM), WDDM evicts -device buffers to system memory, and evicted buffers stay there — every kernel -touching one then reads over PCIe (~7× slower), at full clock. This is the -documented WSL2 "shared GPU memory" / sysmem-fallback behavior (see roadmap -sources). It is **WSL2-only**: native Linux/Windows OOM instead of silently -spilling. Mitigations: the Windows-side driver setting *Prefer No Sysmem -Fallback*, and keeping the committed footprint under the budget. - -**Real workloads are unaffected (measured), so there is nothing to fix there.** -The trigger is allocation *activity* — many *distinct* large allocations in -sequence — not resident level (a single 1 GB alloc, or 3.5 GB held via the pool, -never triggers it). Real inference loads its weights *once* (low churn, high -resident is fine — same as llama.cpp running a 4 GB model on WSL2). Real training -churns activations but at *repeated* shapes, so the footprint stays bounded: a -500-step train-churn loop is **flat** (4.3 ms with the pool, 5.4 ms without — -neither degrades, memory pinned at ~1.8 GB). Only diverse-shape churn — this -bench allocating ~15 distinct multi-hundred-MB sources — climbs the high-water -past the budget. So int4's tok/s is measured in isolation and `bench_llm` stays -at the shorter f32+bf16+attn scope; the direct flush-based benches -(`bench_q4_gemv` etc.) are authoritative. The pool is a pure win for real -training regardless (20% faster here, avoiding cuMemFree's sync). - -The **CUDA buffer pool** (recycle vs cuMemFree, size-keyed, exact match; Metal -already pools) is kept regardless — it does not fix the WSL2 cliff (that's not -allocator-driven) but avoids cuMemFree's synchronization on *every* backend -(cuMemFree stalls the CPU-runs-ahead pipeline — the reason PyTorch/TF cache); -real-MNIST training gate 7.8 → 5.8 s. It helps steady-state workloads (repeated -shapes reuse exactly); for WSL2 a bounded/released pool would be safer, deferred. - -*Refuted (2026-07-04): interleaved (exllama-style) repack for conflict-free -reads.* Assigning lane l the stride-32 keys k = base+l+m·32 (packed one word/lane) -makes both the a and qweight reads coalesced/conflict-free — but measured -**slower**, not faster: shared-a 488 vs the conflicted 565 GB/s, global-a 438 vs -shared 565. So the bank conflict was *not* the binding constraint (a's L2 re-read -per output warp, which shared-a addresses, dominates; the interleaved layout also -reads 8 scales/word vs 1). Reverted. Reaching bandwidth-bound needs a different -lever (cut a-traffic further / dp4a int8 path), found by profiling, not the -obvious coalescing fix — a reminder to measure before assuming. - -**Progression pattern (all three kernels).** Each landed correct-but-un-tuned, -then a second pass to bandwidth-bound: bf16 1.62→1.84×, attention 4→19×, int4 -1.24→1.79× (mid-tune). Occupancy (blocks vs 82 SMs) was the recurring first -bottleneck — split-K/split-KV the recurring fix. - -## Real-model chat: Qwen2.5-0.5B end to end (RTX 3090 / WSL2, 2026-07-10) - -First real GGUF loaded and generating text on own kernels (roadmap M9 "real-model -chat"). The census here is the **decode-speed regime finding**, which overturned -the going-in assumption. - -**Going in:** expected bf16/q4 weight storage to be the decode lever (halve weight -bytes → faster on the bandwidth-bound GEMV) and the WSL2 2GB sysmem cliff to bite -F32. **Both were wrong for a 0.5B model.** Measured decode (`bench/cuda/chat_qwen`, -greedy, all layers on GPU, RTX 3090): - -| weight storage | resident | prefill tok/s | **decode tok/s** | -|---|---|---|---| -| F32 | ~2.0 GB | 48.8 | **67.2** | -| bf16 | ~1.0 GB | 56.9 | **65.0** | -| **llama.cpp (F16, `llama-bench`)** | 1.17 GB | 9765 | **446** | - -**Reading:** (1) bf16 gives **no decode speedup** (65 vs 67) → 0.5B decode is -**per-op-overhead-bound, not bandwidth-bound**. Per-token ≈15 ms is dominated by -the array graph's per-op kernel-launch + eval + D2H across 24 layers × ~15 ops, -not weight-byte traffic. So the M7/M8 bf16/q4 levers (real on bandwidth-bound -7B-class decode) don't move a 0.5B. (2) The **WSL2 sysmem cliff did not fire** at -2 GB F32 resident — decode stayed 67 tok/s, not the ~7× collapse the churn-driven -cliff produces. Confirms [[cuda-wsl2-sysmem-fallback-cliff]]: *resident* weights -(loaded once, never reallocated) are safe; the cliff is a churn phenomenon. (An -earlier ~27 tok/s figure was a load-inclusive mis-estimate; the isolated decode -loop is 67.) (3) **Gap to llama.cpp is 6.7×**, and its cause is overhead, not -kernels: llama.cpp fuses ops and uses CUDA-graph capture to erase per-token -launch/sync. The next decode lever is therefore **per-op overhead** — GPU-side -argmax (kill the per-token 151936-float logits→host D2H that greedy needs only an -argmax from), a fused decode step, graph reuse — *not* dtype width. - -**Correctness (the target-scope bar):** `chat_qwen`'s greedy output is -**character-for-character identical to llama.cpp** on the same F16 GGUF + same GPU -("Give me a short introduction to large language models." → "Large language -models, or LLMs, are artificial intelligence systems that can generate human-like -text…"). Validated tightly in `check_qwen` against a **numpy forward reading the -same F16 weights** (greedy exact; embedding/layer-0/final-norm/top-5 maxrel -~1e-5). Oracle subtlety worth keeping: torch/HF runs the **bf16** original weights -while the GGUF is **F16**, so torch greedy diverges after ~3 tokens (a -weight-precision difference, not a bug) — the numpy-on-F16 forward is the tight -oracle, torch is only a sanity check. head_dim=64 kernels (new {64,128} -generalization) validated at maxrel 1e-7 (`ctest attn64`). - -### Decode-overhead census + the sync-free step (2026-07-10) - -Attacking the 6.7× gap. A per-region wall-clock census of one decode step -(`bench/cuda/speed/bench_qwen_decode.cpp`, `qwenmodel::StepProf`; bf16, RTX 3090/WSL2) split -the 15.9 ms/token as: **construct (host array-graph build) 35–45%**, per-layer -`.eval()` regions (launch + sync) ~58%, and the greedy last-mile (608 KB logits -D2H + 151936-elem host argmax) only **1.7%**. Two structural facts fell out: the -**GPU compute floor is ~1 ms/token** (bf16 weights ≈1 GB ÷ ~940 GB/s) — so **~94% -is overhead, not kernels** — and a big slice of that overhead is **pure host CPU** -(building ~720 array nodes/token), which no amount of GPU work removes. - -Two levers landed, each guarded by `check_qwen` (greedy stays bit-identical to -the numpy-F16 oracle) and measured min-of-rounds (WSL2 boost-clock noise is -±20 tok/s on a single run — the min pins the boosted-clock ceiling, the only way -to resolve a ~15% change): - -- **GPU argmax** (`tl_argmax`, `cuda::argmax`): block reduction over VOCAB → - 4-byte index D2H instead of the 608 KB logits copy + host scan; tie-break - matches the host `v[i] > v[bi]` loop so greedy is bit-identical. Measured - **neutral** (75.1 vs 75.4 tok/s) — the last-mile was only 1.7%, and in - isolation it even *added* a second `CtxSynchronize` (`logits.eval()` + - `argmax`) that on WSL2 costs more than the 608 KB copy it removes. Kept - anyway: it's the **prerequisite for CUDA-graph capture** (a graph can't carry - a per-token 608 KB D2H + host branch). - -- **Sync-free `realize()`** (the real win): `array::realize()` / - `graph::run_noflush` launch an array's graph and adopt its storage **without** - the terminal `gpu::flush()`, leaving kernels in flight on the null stream. The - kv_cache append/attn and the argmax are stream-ordered after their producers, - so they see the writes with **no host sync** — collapsing a step's ~98 - `CtxSynchronize` (4/layer × 24 + tail) to **~1**. Swapping the model's - per-layer `.eval()` → `.realize()` gave **66.2 → 76.7 tok/s (+16%)** on the - identical min-of-rounds harness, greedy unchanged. So decode is **not** - purely construction-bound — the syncs were costing ~2 ms/token after all; - min-of-rounds (not single-run) was needed to see it past the noise. - -**Fused imperative decode step (C1-1, the big win).** After `realize()` the -residual was **host graph construction (~45%) + ~800 kernel launches/token**. -`qwenmodel::step_imperative` runs the whole forward as **direct `cuda::` kernel -calls on a persistent `Scratch` arena** (device buffers allocated once, residual -ping-ponged) — **zero array nodes** (no host graph-build) and **zero per-step -allocation**. Two fused kernels replace array compositions: `tl_rmsnorm` -(sum-of-squares reduction + normalize-scale, matching `(x*x).mean→rsqrt→*w`) and -`tl_swiglu` (`silu(gate)*up`); everything else reuses the existing -gemv/rope/attn/argmax kernels, all on the null stream with the single sync in -`cuda::argmax`. **Decode 77 → ~200 tok/s (+156%, 2.6×)**; `chat_qwen` end-to-end -**208 tok/s**, greedy **character-for-character unchanged** (canonical LLM -answer), and array-vs-imperative agree **0/12** on a same-cache replay -(`bench_qwen_decode`). Net session: **67 → 208 tok/s (3.1×)**, gap to llama.cpp -**6.7× → 2.14×**. This confirmed the census: killing the host construction (and -collapsing the op count) was worth far more than any dtype/bandwidth lever on a -0.5B. - -**CUDA-graph capture POC — the ceiling is only +10% (C1-2, 2026-07-10).** Built -the enabling infra (a `context.stream` all launchers target, dlopen'd -`cuStreamBeginCapture`/`cuStreamEndCapture`/`cuGraphInstantiate`/`cuGraphLaunch`, -and split-K gemv's blocking `MemsetD8` → `MemsetD8Async` so the region is -capturable) and captured the imperative forward at a **fixed pos** (numerically -meaningless, but it isolates launch overhead). Replay + argmax: **221.6 tok/s vs -imperative 201.6 — only +10%.** So after the imperative rewrite decode is -**GPU-kernel-execution-bound, not host-launch-bound**: the graph erases the host -submit sliver (~0.45 ms ÷ 19 launches ≈ 24 µs each), but the GPU still runs -**456 tiny kernels/token** (19/layer × 24) sequentially — ~4.4 ms of per-kernel -GPU-side latency at batch-1, where each M=1 gemv underfills the SMs. **This -redirects the strategy:** the remaining ~2× to llama.cpp (446) is **kernel -efficiency — fewer/fused/bigger kernels** (llama.cpp fuses norm-into-matmul, -whole-attention, and runs far fewer launches), **not** launch overhead. Graph -capture's payoff *shrinks* as you fuse (fewer kernels → less launch overhead to -hide), so it's a ~10% finisher, not the lever. Crucially the POC measured this -**before** sinking effort into the risky `pos`-as-device-scalar plumbing (which -would touch the 1e-7-validated attention kernels) that correct per-token replay -needs. The capture infra is kept (tested, dormant); the next real decode lever -is kernel fusion + tuning (a larger, separate effort). Net decode this session: -**67 → 208 tok/s (3.1×), gap 6.7× → 2.14×**, all greedy bit-identical. - -**Elementwise fusion is neutral → decode is GEMV-execution-bound (2026-07-10).** -Fused the two residual adds into the following RMSNorms (`tl_add_rmsnorm`, writes -the residual sum *and* its norm) and q/k's bias-add into RoPE (`tl_rope` optional -bias) — removing **4 elementwise launches/layer (96/token, 456 → 360 kernels)**. -Result: imperative decode **unchanged (~200 tok/s)**; only the graph ceiling -nudged 221 → 230 (fewer nodes). So the ~12 elementwise kernels/layer were -already cheap; the **7 gemvs/layer (168/token) reading the weights dominate**. -Those gemvs run at **~220 GB/s effective** vs the ~936 GB/s the isolated -`bench_bf16_gemv` hits — the loss is **M=1 small-gemv inefficiency** (short -back-to-back kernels underfill the SMs; split-K adds a `MemsetD8Async` + atomic -combine). llama.cpp's 446 tok/s ≈ **~440 GB/s effective**, so its M=1 -`mul_mat_vec` is ~2× more bandwidth-efficient than ours — **that gemv-kernel -quality gap is the entire remaining ~2×**. The fusion is kept (correct, 0/12 + -oracle, cleaner, lifts the graph path); the neutral result is the point. **Next -decode lever (its own milestone): M=1 dequant-gemv efficiency** — fused QKV / -gate-up (fewer + wider gemvs → occupancy, drop split-K) and a tuned -`mul_mat_vec`-class kernel. **Banked here: 208 tok/s, 2.14× off llama.cpp.** - -**Decode gemv census + QKV/gate-up fusion → 208 → 245 tok/s (2026-07-10).** -Diagnosed *why* in-context gemv runs ~220 GB/s vs the isolated bench's ~936, with a -new isolated bench at the **real Qwen decode shapes** (`bench_qwen_gemv`) plus a -`no_splitk` diagnostic knob (`cuda::set_no_splitk`, forces gy=1). Two hypotheses -**refuted**: (a) forcing split-K OFF makes every small-N gemv **3–16× slower** -(wq 0.30×, wd 0.06×) — split-K is essential, not the culprit; (c) the back-to-back -7-gemv layer sequence is **−2.4%** vs the sum of isolated medians — inter-kernel -cadence is ~0 (matches the graph-capture +10% finding). The real cause is a -**fixed ~14µs per-launch floor** that dominates at small N: `wk` (896×128, 229KB) -takes 14.4µs = **16 GB/s** (pure floor, bytes irrelevant); `wg` (896×4864, 8.7MB) -takes 20.2µs for 38× the bytes. The floor amortizes only at large N — `lm_head` -(896×151936) hits **897 GB/s** (≈peak), and `gateup.fused` (896×9728) **598** vs -wg+wu separately (431/357). So the lever is **fewer kernels**, not a better kernel -(a warp-per-row rewrite still can't fill the SMs at N=896 — the `no_splitk` column -proves the occupancy loss is shape-intrinsic). **Fix: concat wq|wk|wv → [896,1152] -and wg|wu → [896,9728] at load (`load_w_T_cat`), one GEMV each + device-pointer -slice** (mid-buffer pointers flow through `off_`; slice reads need no host sync as -the fused GEMV already marked the base device-live). Per-column split-K is -identical (bx=1, chunk=32 in both fused and separate) → **bit-identical per -column**. Imperative path only; the array path keeps the separate loads as the -oracle, so `bench_qwen_decode`'s array-vs-imperative is a free independent guard -on the fusion. Result: **imperative 208 → 245 tok/s (+18%)**, graph ceiling 230 → -257, `chat_qwen` end-to-end ~220, **0/12 EXACT + check_qwen oracle green + 8/8 -ctests**. **Gap to llama.cpp 446: 2.14× → 1.82×.** +466MB resident (both weight -copies) stays under the ~2GB WSL2 cliff. Remaining ~1.8×: the 5 gemvs/layer still -pay the per-launch floor (QKV.fused only 107 GB/s at N=1152) — next levers are a -lower-floor `mul_mat_vec`-class kernel and/or graph capture (now a real +5%). - -**Warp-per-row [N,K] gemv (lever A) → 245 → ~330 tok/s (2026-07-10).** The -mul_mat_vec-class kernel above. The split-K kernel is column-per-thread over -`[K,N]` weights and pays, every launch, a `MemsetD8Async` prezero + an `atomicAdd` -combine across the K-splits + an underfilled grid — all fixed cost that dwarfs the -tiny byte count at small N. New `tl_gemv_bf16_row`: weights in `[N,K]` (K -contiguous per output row = **GGML-native, so `load_w_T_row`/`load_w_T_cat_row` -drop the transpose — just a widen**), **one warp per output row**, its 32 lanes -split K via 16-byte `uint4` loads (consecutive lanes → consecutive 16B → one -coalesced 512B warp transaction), shuffle-reduce, single store. **No split-K → no -memset, no atomic, one clean launch.** Only `K % 8 == 0` needed (every dim); the -partial last 256-block is handled by the per-lane `k0 < K` guard (K=896 = 3×256+128 -→ tail lanes skip; no `K % 256` requirement like q4). Diagnosis-first held: added a -`row` column to `bench_qwen_gemv` (reuses the same `B` buffer — speed is -layout-agnostic) and **row beat split-K on *every* Qwen shape**: QKV.fused **1.89×** -(133→253 GB/s), wo **1.79×** (89→160), wd **1.68×** (322→540), gateup **1.36×** -(593→804), wk **1.91×**; `lm_head` **neutral** (1.02×, both ~950 GB/s = at HBM -peak — never floor-bound). Integrated **bf16 imperative only** (`Model.row` flag + -a `gv()` dispatch lambda in `run_layers_`): `wqkv`/`wgu` become row `[N,K]` (net 0 -memory), plus `wo_row`/`wd_row` row copies (+247MB); **`lm_head` stays split-K** -(neutral + skips the +272MB `[N,K]` copy) → 1644 MiB resident, under the cliff. The -array `[K,N]` split-K path is untouched as the oracle; because the row reduction -order differs it is **greedy-equivalent, not bit-identical** — yet -`bench_qwen_decode` still reports **0/12 EXACT**. Result: **imperative 245 → ~330 -tok/s (+35%, min of 6 rounds, stable 327–334)**, graph ceiling 343–357 (now a real -+5–8% finisher), `chat_qwen` coherent, **check_qwen greedy MATCH + 8/8 ctests**, -f32 fallback path 0/12 EXACT. **Gap to llama.cpp 446: 1.82× → ~1.35×; net session -67 → 330 tok/s (4.9×).** The per-launch floor is now gone, so the next real lever -is **bandwidth**: Q4_0 → `dtype::q4` (0.625 B/weight vs bf16's 2). - -**Q4 (group-int4) bandwidth lever → +6% at the ceiling, kept opt-in (2026-07-10).** -With the per-launch floor killed by lever A, the remaining decode cost is HBM -traffic, so the next lever is the dtype: `dtype::q4` reads **~0.625 B/weight** vs -bf16's 2 (group-32 symmetric int4, `scale = amax/7`; one buffer = `[packed int4 -[N][K/2] | scales f32 [N][K/32]]`). The q4 kernel (`tl_gemv_q4`) was already the -warp-per-row `[N,K]` shape lever A landed on, so it drops straight in — the only -change was generalizing its tail from `K % 256` to `K % 32` (a per-lane `k0 < K` -guard on the uint32 read + 8-nibble loop; bit-exact on `K % 256 == 0` shapes since -the guard is always-true there, verified via `bench_q4` maxrel 6e-6). **Diagnosis -first held again**: a `q4` column added to `bench_qwen_gemv` showed the expected -**hybrid** — q4 wins the large-N bandwidth-bound gemvs (MLP gate/up/down, lm_head) -but *loses* the small-N attention projections, where lever A already erased the -floor and the dequant overhead now exceeds the byte savings. So integrated **q4 for -the MLP gemvs only** (`wgu`/`wd`), attention stays bf16-row, lm_head optional. -Measured (graph ceiling = the noise-free axis; WSL2 boost noise ±20 tok/s masks it -at the imperative level): baseline **357** → q4-MLP **379 (+6.2%)** → q4-all (+lm_head) -**379.5** — lm_head q4 is marginal (it's already at HBM peak as split-K, so cutting -its bytes barely moves it) while costing the most quality. q4-MLP also frees **~430MB** -(replaces the bf16 row copies). **Cost is quality**: q4 is lossy, so greedy diverges -from the bf16 oracle **2/12** (q4-MLP) / **3/12** (q4-all), top-logit maxrel 3.15e-2 / -3.94e-2 — `chat_qwen` stays coherent for all configs, but the token stream can -differ. **Decision: default stays bf16-row (bit-exact oracle-matching); q4 is -opt-in** (`--q4=mlp` / `lm` / `all`) for users who want the +6% and −430MB and accept -the lossiness. Array `[K,N]` split-K path untouched as the bf16 oracle → array-vs- -imperative divergence is the built-in quality gauge. Committed 478c5d4. - -**bf16 KV cache → REFUTED for Qwen 0.5B (neutral), real 1.64× at llama-7B; kept -dormant (2026-07-16).** The flagged next decode lever: store the KV cache as bf16 -(2 B) so the attention kernels stream half the K,V bytes every step (~2× the KV -floor). Built as a storage-dtype seam that leaves the f32 path byte-identical: the -attention read cores (`attn_decode_core`, `_split_core`, `attn_prefill_core`) and -the write kernels (`tl_kv_append`, `tl_kv_fill`) gained a `typename KT ∈ {float, -__nv_bfloat16}` template with a `kv_ld`/`kv_st` widen/narrow seam; new bf16 symbols -(`tl_attn_decode_bf16[_64]`, `…_split_bf16[_64]`, `…_prefill_bf16[_64]`, -`tl_kv_{append,fill}_bf16`) sit beside the untouched f32 ones, and `kv_cache` gained -a `dtype` flag (default f32) that routes to the matching instantiation — q/out/ -scratch stay f32; all indices/`kv_stride` are element counts so only the width -changes. **Measure-first (`bench_attn_decode`, bf16 column) held the line**: at the -**llama-7B shape** (H=32, D=128) bf16 KV is a real lever — 1.05× at ctx=512 growing -to **1.64× at ctx≥2048** (901→737 GB/s at half the bytes), maxrel vs the f32 CPU -ref only **~1e-4** (bf16's 8-bit mantissa + attention's averaging, ~300× tighter -than q4). But at the **REAL Qwen2.5-0.5B shape** (14 q / **2 kv** heads, D=64) it is -**NEUTRAL — 1.00× (even 0.98×)** across every realistic ctx 128–2048: per-layer attn -is a flat ~28–32 µs regardless of dtype because the 2-KV-head footprint is tiny and -the kernel is **launch/occupancy-floor-bound, not bandwidth-bound** (the same lesson -as the whole session — on 0.5B decode the weight-reading gemvs are the cost, 2-head -attention is a rounding error). So per the "土俵で測ってから積む" discipline it is -**NOT wired into the Qwen decode path** (neutral + lossy = a strict loss there). The -infra is kept **tested + dormant** (f32 path byte-identical, 8/8 ctest green, the -arena bench exercises + documents where it pays), a real lever the moment a -llama-class model (many KV heads, long ctx) is added — mirroring the graph-capture -"kept tested + dormant" precedent. - -## vs silarray (M1 Pro, 2026-07-03) - -Head-to-head with the predecessor across cpu/gpu/auto. Two separate -binaries (header coexistence blocked by objc/Metal bridge collisions), -run alternating with per-process warmup, steady-state min of 3 rounds, -load ~3. Both use blocking eval (silarray's free `sil::eval()`, not the -non-flushing member `.eval()`). `tl/sil` < 1 means tensorlib faster. - -| Workload | cpu | gpu | auto | -|---|---|---|---| -| SGEMM 1024³ | 1.01 | **0.76** | 0.92 | -| SGEMM 256³ | 1.00 | 0.98 | **0.10** | -| elementwise 4M | 0.90 | 1.04 | 0.95 | -| softmax 1024² | **0.03** | 1.42 | **0.003** | -| MLP fwd 256×784×256×10 | **0.38** | 1.67 | 0.55 | - -Reading: -- **CPU**: tensorlib parity-to-faster. matmul parity (shared AMX). The - huge softmax/MLP wins are because **silarray's CPU softmax is a naive - correctness-only fallback** (its docs say so — silarray is GPU-focused) - — discount these as "sil didn't implement it," not an engineering edge. -- **GPU**: split. tensorlib **wins matmul** (STEEL port, 1024³ 24% - faster); silarray **wins softmax (1.4×) and the fused MLP (1.7×)** — its - online-softmax and `linear_sigmoid` GPU kernels are more mature. This is - the motivation for the deferred "GPU fused kernels" item. -- **auto**: tensorlib clearly better — the finalized thresholds route - correctly (256³→cpu 10×, softmax→gpu 300×) where silarray's auto - mis-picks (256³→slow GPU, softmax→slow CPU). auto quality is the payoff - of the measured thresholds. - -Caveats: single-op micro-bench, cross-process (looser than same-run). - -## Refuted approaches — do not retry without new evidence - -| Approach | Verdict | Data | -|---|---|---| -| Naive 64×64 tile via the simple `sgemm_body_` template (full `simdgroup_matrix` locals, 16 accumulators) | **Catastrophic** | 133 GFLOP/s at 2048³ vs 2175 for 32×32 — the "16-accumulator occupancy crisis" silarray hit; STEEL's explicit float2 fragment registers are what make 16 accumulators viable | -| Host-side alignment gating of float4 loaders (require ld%4==0, 16B base) | **Unnecessary** | silarray/MLX issue vector loads at any alignment on Apple GPUs in practice (bench shapes with ldb=50); gate removed | - -Inherited from silarray (Metal; re-verify before assuming they transfer to -CUDA): spin-wait on command-buffer status (slower), split commits to overlap -CPU/GPU (slower), trace/replay compile (no payoff — graph build was 0.2% of -step time), AOT vs JIT kernel compilation (no effect on SGEMM). - -## Bug classes worth remembering - -- silarray's only two real correctness bugs were the same class: copy-pasted - kernel bodies diverging between aligned and edge tile paths. Rule: kernel - bodies are shared templates; every fused-epilogue kernel gets edge-shape - numerical tests (M-edge, N-edge, band boundaries). diff --git a/docs/roadmap.md b/docs/roadmap.md deleted file mode 100644 index 19527f3..0000000 --- a/docs/roadmap.md +++ /dev/null @@ -1,599 +0,0 @@ -# cpp-tensorlib roadmap - -Forward-looking plan: milestone scope, environment/hardware constraints, the -approach for each remaining milestone, and open decisions. Complements -[architecture.md](architecture.md) (how the code works today) and -[performance-notes.md](performance-notes.md) (measurement methodology, gate -results, refuted approaches). This is the portable source of truth for the -plan — it travels with the repo, so any checkout on any machine has it. - -## Where things stand - -M1–M4 are done: reference backend, lazy graph + fusion, the Accelerate and -Metal (STEEL SGEMM) backends behind a single dispatch seam, culebra -integration, and a tiny-tensor per-op sprint. The macOS backend is complete -and competitive (see the silarray comparison in performance-notes.md). The -non-Apple backends have landed their F32 foundation (M5 CPU met the OpenBLAS -gate; M6 CUDA has a correct, split-K SGEMM at ~0.82 of cuBLAS). - -The decode path (the target-scope money regime) is now built out (2026-07-04): -**M7** bf16 storage + decode GEMV, and **M9** fused decode attention, are done -and integrated end to end — CUDA decode went **3.5 → 41.2 tok/s** on llama-7B -shapes (RTX 3090). **M8** int4 dequant GEMV has a correct kernel (1.79× vs bf16) -but is not yet array-integrated or bandwidth-bound. What remains: M8's tuning + -integration, M7's bf16 **compute** GEMM (fine-tuning / prefill, Tensor Core), -M9's KV cache / GQA / causal-mask, and a tail of deferred perf/CI/API work. - -Design invariants that constrain everything below: - -- **Zero third-party runtime dependencies.** macOS links only OS frameworks; - CPU/CUDA use own kernels (no OpenBLAS/cuBLAS/CUTLASS); the CUDA driver is - `dlopen`'d. doctest is vendored (tests only). -- **One dispatch seam.** Backends plug in at `graph::eval_one`; each returns - nullopt/false when ineligible, so the evaluator has no platform `#ifdef`s. - A new backend is a new set of `*_mode_()` guards + kernels, nothing else. -- **The reference `ref::` backend is permanent** — correctness oracle and the - fallback for any unsupported shape/target. -- **F32 is the compute baseline**; BF16/FP16 storage *and* compute arrive in M7, - and int4/int8 as a *storage+dequant* type in M8 (compute stays F32/BF16-accum, - no F64). Widening the dtype surface is now in scope — see below. - -## Target scope & the bar - -The project targets **consumer/prosumer local AI**: running and fine-tuning -neural nets — Local LLM inference (chat/decode) and LoRA/QLoRA fine-tuning — on -**one to a few consumer GPUs** (RTX 5090 = Blackwell sm_120; A6000 = Ampere -sm_86). Large-scale multi-node LLM *training* is explicitly **out of scope**. - -This scope resets the performance bar, and it is not "beat cuBLAS at 4096³": - -- **The relevant competitor is llama.cpp / exllamav2 / ggml / MLC-LLM, not - PyTorch+cuBLAS.** Those are the SOTA local runtimes, and they are hand-written - with **no cuBLAS on the hot paths** — proof that own-kernels is the *right* - approach here, not a compromise. The operative metric is **tokens/sec on a - real quantized model**, not GFLOP/s vs cuBLAS. -- **Interactive decode (the dominant cost) is memory-bandwidth-bound** - GEMV / dequant-matmul (batch≈1, each weight streamed once), where cuBLAS has - **no moat** — a hand kernel that saturates VRAM bandwidth matches or beats it. - Attention is flash-attention-style fused work cuBLAS doesn't even do. Only - **prefill / larger batch** is compute-bound dense GEMM (the regime where the - hand SGEMM tops ~0.90 of cuBLAS — see M6), and it is a one-time, minority cost - for interactive use. -- Therefore the M6 **cuBLAS-90%-on-square gate is consciously de-prioritized** - (measured, reached ~0.82 with split-K, kept as the prefill foundation — not - chased with a CUTLASS escape hatch the scope doesn't need). The high-leverage - work is the dtype/quant/attention kernels (M7–M9), all hand-written-friendly - and cuBLAS-irrelevant. Full analysis (the decode=bandwidth-bound regime split, - the silarray/llama.cpp precedent) is in performance-notes.md. - -## Milestones - -### M5 — Own CPU backend (Linux / Windows) 🔨 in progress - -BLIS-style GEMM: cache-blocking + packing (architecture-independent C) around -a register-blocked SIMD microkernel per ISA (AVX2, AVX-512, NEON), selected -at runtime by CPUID. Plus a small threadpool for the M/N parallel split, and -own SIMD elementwise/reduction kernels. Gate: ≥90% of OpenBLAS and ≥ PyTorch -CPU on the target ISA. Escape hatch: a `-DUSE_OPENBLAS` CMake option for any -shape band that can't meet the gate. - -**Done (first cut):** `cpu.h` + `cpu_threadpool.h` — the full BLIS structure -(blocking, stride-aware packing, persistent thread pool, M-parallel), an 8×8 -NEON microkernel + a portable scalar fallback, wired into the dot dispatch -(`metal → accel → cpu → ref`) and gated by `cpu::enabled_`. Correct on every -shape/transpose/epilogue (own == ref oracle test); ~235 GFLOP/s at 1024³ on -M1 Pro NEON. Elementwise/reductions stay on array.h's autovectorized flat -loops (memory-bound), by design. - -**Done (NEON tuning pass, 2026-07-03):** lane-indexed-FMA microkernel -(K-unroll ×4 + prefetch), thread-local reusable pack buffers, MR-panel- -granular M-parallelism, KC 256→512. Quiet-census verified: **~86–89% of -NEON fp32 peak at 2048³, ~84% at 1024³** (was 54–64%); the sub-256³ loss -to the naive loop is fixed (128³ now 2.4× ref). KC=512 confirmed and an -8×12 kernel measured-and-rejected (within noise). The on-Mac NEON tuning -is complete — see performance-notes.md "tuning pass" for the census and -the rejected 8×12. - -**Done (AVX2 + runtime dispatch scaffold, 2026-07-03):** the compile-time -`#ifdef` microkernel pick is now a runtime function pointer (`select_ukernel`) -— ARM fixes to NEON at compile time; x86 probes `__builtin_cpu_supports` -and picks AVX2 (`target("avx2,fma")`) or the scalar fallback. Added an AVX2 -8×8 microkernel (`_mm256_fmadd_ps` + `_mm256_broadcast_ss`, the x86 analogue -of the NEON lane-FMA kernel, K-unroll ×4 + prefetch), sharing the MR=8/NR=8 -packed layout. Scalar + NEON validated against a naive oracle on ARM; the -AVX2 path cross-compiled (`-target x86_64 … -c`) and disassembly-confirmed -to emit `vfmadd*`/`vbroadcastss` — but **not executed** (Rosetta stops at -SSE4.2). Full suite green on native ARM. - -**Done (AVX2 execution + 6×16 tuning + gate, 2026-07-03):** first real run on -the i7-12700KF WSL2 box (see performance-notes.md "x86 census" + the runbook). -AVX2 numerically verified vs the naive oracle and the full suite. Tuned the -tile **8×8 → 6×16** (12 ymm accumulators): single-core **~91% of AVX2 peak** -(was ~77%), matching NEON; **own ≈ 106% of OpenBLAS at 2048³ — gate MET**. Each -kernel now packs to its own tile via `ukernel_desc` (the anticipated second -packing layout; AVX-512 reuses NR=16). KC=512 re-confirmed for x86; 8×8 kept -behind `-DTL_CPU_AVX2_8X8` for A/B. - -**Remaining:** -- *AVX-512 microkernel* — deferred: a real AVX-512 kernel wants NR=16 — now a - ready-made layout (the 6×16 AVX2 kernel packs NR=16), so `select_ukernel` - needs only an `avx512f` branch + the wider ukernel. This box has no AVX-512; - needs a Sapphire/Ice Lake / Zen4 box. -- *Problem-size-aware thread count* — see the deferred table: the pool statically - splits M across all 20 threads, over-parallelizing mid-size GEMM (512³/1024³ - lose to OpenBLAS, which caps threads for small problems). Single-core is 91% - of peak, so this is a thread-policy gap, not a kernel one. Best measured on - native Linux (WSL2 hides the P/E topology the heuristic needs). - -**Environment:** Apple Silicon is ARM64 = NEON, so the NEON microkernel is -fully developable, natively runnable, and perf-tunable on a Mac (it is also -the production path on ARM Linux — Graviton/Ampere). The **AVX2/AVX-512 -microkernels can be written and compile-checked on a Mac but not executed** -(Rosetta 2 stops at SSE4.2; there is no native x86). Validating the x86 path -needs x86 Linux — the WSL2 box (see M6) or x86 CI (correctness only; CI perf -is too noisy for the gate). Note the macOS *production* CPU path stays -Accelerate (AMX beats hand-written kernels here); M5 on a Mac is validated -against the `ref::` baseline, and its real payoff is off-Apple. - -**Approach:** scaffolding (blocking loops, packing, threadpool, CPUID -dispatch) and a correct NEON microkernel first — a naive-but-correct kernel -should already beat `ref::`. Then a dedicated tuning pass to reach the -OpenBLAS-90% gate (register blocking, FMA scheduling, packing layout) — this -is the performance-sensitive part; treat it like the STEEL sprint. Defer -AVX tuning to the x86 box. - -### M6 — Own CUDA backend (Linux / Windows) - -FP32 SIMT SGEMM (own kernels, no cuBLAS/CUTLASS), driver API `dlopen`'d so a -binary runs — and falls back to CPU — without a driver present. Kernels AOT -to PTX/fatbin and `#embed`'d, loaded via the driver API. Gate: ≥90% of cuBLAS -FP32 and ≥ PyTorch CUDA. Escape hatch: per-shape-band CUTLASS instances if a -band can't meet the gate. - -**Environment:** needs an NVIDIA GPU. Available: **RTX 3090 (sm_86) on Ubuntu -under WSL2.** The WSL2 CUDA driver lives at `/usr/lib/wsl/lib/libcuda.so` — -the dlopen loader's search path must include it. WSL2 adds submission latency -vs native Linux; record which side of that line a measurement came from. -Windows-native is not yet set up. Real numerical/perf validation is manual on -this box (no self-hosted or paid GPU CI — see the CI note). CI covers the -CUDA *build* + the driver-absent CPU-fallback path only. - -**Approach:** two stages. (1) Scaffolding — dlopen loader (with the WSL2 -path), PTX `#embed` + driver-API kernel launch, build integration -(`TENSORLIB_CUDA` option, only when nvcc is found), elementwise/reduction -kernels. This mirrors the Metal M3b-1 foundation. (2) The SGEMM ladder tuned -to the cuBLAS-90% gate — the performance-sensitive sprint. The x86 WSL2 box -also validates the M5 AVX path. - -**Done (stage-1 kernels + cuda.h backend, validated on GPU, 2026-07-03):** -`kernels/tensorlib_cuda.cu` (elementwise add/sub/mul/div, unary exp/log/sqrt/ -sigmoid/relu/affine, row_sum/row_max/softmax, and a correctness-first SGEMM -with trans/ld support) compiles to PTX via nvcc 13.0 for sm_86. `include/cuda.h` -is the metal.h analogue: a hand-declared driver API loaded by dlopen, primary- -context + module-from-PTX loader, managed alloc/release, and binary/unary/gemm/ -row_op launchers — same signatures as metal.h so the seam is backend-agnostic, -reusing `metal::kop`. The PTX is embedded via **bin2c** (a build-generated -`tensorlib_cuda_ptx.inc` byte array), *not* C23 `#embed`, because the off-Apple -compilers (g++ 11 / clang 14) predate `#embed`. `bench/cuda/check/check_cuda.cpp` validates -all ops on the RTX 3090 vs CPU references (binary+affine, sigmoid, GEMM NN + NT, -row_sum+affine, softmax — all ok). Stubs compile cleanly with no `TENSORLIB_CUDA` -and on Apple, so non-CUDA builds are unaffected. - -**Done (stage-1 build + seam integration, 2026-07-03):** the `TENSORLIB_CUDA` -CMake path compiles the kernels to PTX (nvcc), bin2c's them to a byte-array -`.inc` on the include path, defines `TENSORLIB_CUDA`, and links `libdl` — no -CUDA runtime link (`cmake/bin2c.cmake` is the portable, xxd-free converter). -The eval seam is now a **`gpu::` facade** (in cuda.h: `namespace gpu = cuda` -off-Apple, `= metal` on Apple / non-CUDA builds), so array.h renamed its -`metal_*` helpers to `gpu_*` and its `metal::` calls to `gpu::` — one alias is -the only platform `#ifdef`. `storage::make_device_` uses `gpu::alloc`, so CUDA -managed memory flows in with no extra branch. **`ctest` cpu/gpu/auto + a -`cuda_ukernel` oracle all pass with the real CUDA backend on the RTX 3090** -(gpu/auto route to CUDA; `cuda_ukernel` skips-as-pass with no device, matching -the driver-absent CI fallback). Non-CUDA build unchanged (3 tests green). The -old `kernels/smoke.cu` build-check is removed — the real kernel → PTX compile -now serves that role. - -**Done (stage-2 memory-model pivot, 2026-07-03):** standing up the SGEMM census -exposed that the managed-memory model (below) collapses compute-bound GEMM ~88× -on WSL2 (`ConcurrentManagedAccess=0` → no page migration; prefetch/advise both -"invalid device ordinal"). Since no kernel tuning matters on managed memory and -PyTorch uses device memory, the backend moved to the pre-authorized **explicit -device-buffer model**: a persistent host/device **mirror** per allocation (host -malloc + `cuMemAlloc` device buffer; dirty state keyed by the device pointer in -cuda.h so views sharing a storage share one entry), lazy H2D before a kernel -reads a host-dirty buffer, D2H before the CPU reads a device-dirty one -(`array::raw()/data()` → `gpu::sync_to_host`). Metal stays a strict no-op -(genuinely unified). storage.h is unchanged; the seam addition is two lines in -array.h. ctest cpu/gpu/auto + the cuda_ukernel oracle validate the coherence. -See performance-notes.md "Own CUDA GEMM stage-2" and [[cuda-wsl2-managed-memory-cliff]]. - -**Done (stage-2 SGEMM ladder, first pass, 2026-07-03):** `tl_sgemm_rb` — a -128×128×8 blocktile / 64×32 warp-tile / 8×8-microtile kernel (float4 global -loads, register-staged double buffer, conflict-free warp-broadcast smem reads), -dispatched as the fast path for NN-contiguous K%8/N%4/16B-aligned shapes; -everything else keeps the correctness-first `tl_sgemm`. Ladder (interleaved -own/cuBLAS): naive ~0.09 → register-block 0.68/0.75 → +warp-tile 0.71/0.80 → -+double-buffer **0.72–0.84 (2048³) / 0.79–0.83 (4096³)**. Correct throughout -(max rel ≤ 6e-5). **Gate not yet met:** own ≈ 80% of cuBLAS, ≈ 75% of PyTorch- -FP32 (which bundles a faster cuBLAS 12.8). BK=16 measured-and-rejected. - -**Done (stage-2 ladder ① rejected + ② split-K landed, 2026-07-04):** the -128×256 large-tile rung (①) was built (`tl_sgemm_rb2`, 8×16 microtile) and -**rejected** — 128 accumulators → 203 regs → 1 block/SM (16.7% occ vs rb's 33%), -losing at every size *including* 4096³ (0.81→0.75, where wave-quant is a -non-issue), so the arithmetic-intensity gain is dwarfed by the occupancy loss; -big spatial tiles are ceiling-capped at rb on sm_86 (see the refuted table). The -census reframed the deficit as small-size *fill* (own/cuB rising with size), so -**split-K (②) landed**: a single `ksplit`/`blockIdx.z` addition to `tl_sgemm_rb` -partitions K into S z-slices (S× more blocks), atomicAdd'ing partials into a -zeroed C (identity scale/offset only). S=2 is the measured optimum — 1024³ -0.63→0.75, 2048³ 0.74→0.76; 4096³ stays S=1 (6 waves already, split is pure -overhead). Auto: S=2 for base<512 blocks & K≥512. Post: 1024³ ~0.72–0.82, 2048³ -0.76, 4096³ 0.82–0.83. Still 127 regs / 2 blocks/SM; ctest cpu/gpu/auto green. - -**Done (stage-2 ladder ③ cp.async rejected, 2026-07-04):** the Ampere `cp.async` -staged pipeline (`tl_sgemm_cp`, 3/4-stage) was built and **rejected** — it loses -at 4096³ (0.83→0.73) and 2048³ because the transposed-A smem layout the fragment -read needs can't be filled by an efficient 16-byte cp.async (cp.async can't -transpose), forcing 4×4-byte scatter (the slow granularity) + an extra barrier. -Efficient cp.async here needs `ldmatrix` (tensor-core-only). See the refuted -table. The register-staged `tl_sgemm_rb` load path stands. - -**Remaining — the efficiency ceiling (~0.82 at 4096³) is the open gate.** All -three documented rungs are now spent: ① 128×256 large tile (rejected, occupancy -cliff), ② split-K (**landed**, fixed the small-size fill deficit), ③ cp.async -(rejected, transpose mismatch). The remaining ~8–18% to the cuBLAS-90% / beat- -PyTorch gate is per-SM kernel efficiency at large sizes, where the hand-written -SIMT kernel structurally trails cuBLAS (which uses tensor-core-adjacent -`ldmatrix`/mma paths even for its FP32 SGEMM on Ampere). Options: (a) the -pre-authorized **per-shape CUTLASS escape hatch** for the large band; (b) accept -~0.82 and document the band as a known gap; (c) further micro-scheduling with -uncertain payoff. - -**Resolved (2026-07-04): option (b) — accept ~0.82, no CUTLASS.** Per the target -scope, large-square compute-bound GEMM is *prefill only* (a minority, one-time -cost); the dominant local-LLM cost is bandwidth-bound decode where cuBLAS has no -moat. So the ~0.82 split-K SGEMM is kept as the **prefill/compute-bound -foundation** and the cuBLAS-90% chase is dropped. The escape hatch stays -*pre-authorized but unused* — if a real workload ever proves prefill-GEMM-bound, -reopen it (with the binary-size feature-flag design: CUTLASS confined to the -build-time kernel `.cu`, `-DTENSORLIB_CUDA_CUTLASS` off by default, so -non-matmul builds carry zero CUTLASS bytes). Reaching ~0.90 by hand is also -possible (Boehm hits ~0.93 on sm_86 via a full tile autotune + fragment-register -pipeline + swizzled cp.async) but is deferred as lower-leverage than M7–M9. -Also unchanged: SteelLoaderT-equivalent for transposed operands (backward -matmuls use the naive `tl_sgemm` today). All measured on the RTX 3090, WSL2 noted. - -**Done (loader probe + memory model, 2026-07-03):** the dlopen'd-driver design -is validated on the RTX 3090 box — a standalone probe declares the driver API -itself (no CUDA headers), dlopens `/usr/lib/wsl/lib/libcuda.so.1`, and reaches -the device (`cuInit` → driver API 13.1, RTX 3090 sm_86, 82 SMs) with **no CUDA -Toolkit installed**. So the host side needs only libdl; nvcc is a build-time- -only prereq for the kernel PTX (`cuda-toolkit-13-0`, matching the 13.1 driver). -**Memory model (SUPERSEDED at the stage-2 gate — see the memory-model pivot -above): originally CUDA managed memory** (`cuMemAllocManaged`) — one allocation -visible to host and device, reusing the Apple unified-memory seam, the driver -migrating pages on demand. Chosen for minimal churn + correctness-first. The -explicit device-buffer + H2D/D2H alternative was deferred *iff* managed-memory -page migration proved a bottleneck (prefetch the first lever). It did — on WSL2 -`ConcurrentManagedAccess=0` means pages never migrate and prefetch is -unavailable, so the device-buffer mirror replaced managed at stage 2. - -### M7 — BF16 / FP16 storage **and compute** 🔨 storage done; compute deferred - -BF16/FP16 as a *storage* type first (memory is 16-bit, load widens to F32, store -narrows — a bit-shift confined to load/store, compute stays F32, no dtype × ISA × -kernel explosion). Halves bandwidth on memory-bound ops — directly the decode -lever. Pairs with a language-level dtype surface in culebra. - -**Done (storage + decode GEMV, 2026-07-04):** `types.h` `dtype{f32,bf16}` + -RNE converters; `storage`/`array` carry a dtype (byte-sized alloc). bf16 is a -weight-container type (`to_bf16()`/`to_f32()`) — consumed natively by the CUDA -decode GEMV, widened to F32 for every other op in `eval_one`'s input funnel -(ref/cpu/accel/Metal untouched). Two levers landed together: routing M=1 dots to -a GEMV kernel (`tl_gemv_f32`/`tl_gemv_bf16v8`, the 128² tile wastes 127 rows at -M=1) — F32 decode 3.5 → 12.2 tok/s — and bf16 weights (~1.84× on the -bandwidth-bound GEMV) → 14.8 tok/s. Census in performance-notes.md. - -**Remaining:** bf16 **compute** GEMM with F32-accumulate — the fine-tuning / -prefill path, and the Tensor-Core option below. Deferred as lower-leverage than -M8/M9 for decode; scheduled when the fine-tuning workload is exercised. - -**The compute path is now in scope (was "later, demand-driven" — the scope -demands it).** LoRA/QLoRA fine-tuning trains in BF16/FP16 with F32 accumulation, -and a 16-bit compute GEMM halves weight bandwidth for decode too. So M7 adds a -BF16/FP16 GEMM with F32-accumulate: on CUDA the `tl_sgemm_rb` register-blocked -kernel generalizes (16-bit smem tiles, F32 accumulators — the microtile/warp -structure is unchanged); on Ampere+ a Tensor-Core (`mma.sync` bf16→f32) path is -the fast option and, unlike the FP32 SGEMM, here Tensor Cores are the *native* -datapath so the hand kernel is on equal footing with cuBLAS. FP8 (sm_120 / -Blackwell) is a forward-looking extension of the same seam. - -### M8 — Quantized (int4/int8) dequant-fused matmul 🔨 int4 done & integrated; tuning + int8/GGUF remain - -> **Update 2026-07-17:** the `tl_gemv_q4s` shared-a variant named below was -> removed; `tl_gemv_q4` was rewritten to one-block-per-row + K-adaptive block -> size (llama.cpp's `mul_mat_vec_f` strategy), which makes the shared-a staging -> pointless (each thread reads a disjoint K-slice, no intra-block `a` reuse -> left to amortize). See the M9 decode-speed notes and commits 63a9064 / -> 12f423a / a426bb0. The 2026-07-04 analysis below is kept as the point-in-time -> record. - -**Done (int4 GEMV kernel + integration, 2026-07-04):** `tl_gemv_q4`/`tl_gemv_q4s` -— group-symmetric int4 weights in [N,K] layout (groups along K contiguous, -GGUF/GPTQ convention), one warp per output, dequant in registers -(scale·(nibble−8)), F32 accumulate. Correct vs a host dequant reference (maxrel -~1e-5); ~0.625 bytes/weight vs bf16's 2.0 → **1.79× vs bf16** on the decode GEMV. -Integrated as `dtype::q4` (chose the storage-dtype path B over a separate -q4_matrix type — see performance-notes.md): one buffer = packed [N,K] int4 + -appended scales, logical shape stays [K,N] so `a.dot(Wq)` type-checks like f32 -and rides the bf16 widen-fallback seam. `array::to_q4()`/`to_f32()`. **Decode -61.7 tok/s** (vs bf16 42.5). Also landed a **CUDA buffer pool** (recycle vs -cuMemFree; real-MNIST gate 7.8→5.8 s). **Remaining:** not yet bandwidth-bound -(best ~565/936 GB/s — the interleaved-repack fix was measured *worse*, see -refuted note; needs a profiled lever), int8, and GGUF-format interop. See -performance-notes.md. - -The heart of consumer local-LLM inference: 4-bit/8-bit weights (VRAM fit + -bandwidth) with a **fused dequantize-and-matmul** kernel — read packed int4/int8 -+ per-group scales, widen to BF16/F32 in registers, accumulate. This is a -*storage+dequant* dtype (activations stay BF16/F32; only weights are quantized), -so it rides the M7 dtype surface and the existing GEMM/GEMV structure with a -widen-on-load epilogue. Decode is GEMV-shaped and bandwidth-bound, so the win is -reading 4× fewer weight bytes — a regime cuBLAS doesn't serve and where -llama.cpp/exllamav2 hand-write. GGUF-style group quantization (Q4_K etc.) is the -reference format to stay interop-friendly. Gate: **tokens/sec within reach of -llama.cpp/exllamav2** on the same quantized model + GPU (not GFLOP/s vs cuBLAS). - -### M9 — Fused attention + KV cache 🔨 decode + KV cache + GQA + causal prefill + block-wiring + multi-layer decode loop + real-model chat (Qwen2.5-0.5B) + sync-free step + fused imperative decode step (67→208 tok/s, 2.14× off llama.cpp) + graph-capture ceiling measured (+10% only → kernel-bound) done; kernel fusion/tuning + prefill tiling remain - -**Done (fused decode attention, 2026-07-04):** `array::attn_decode(q,K,V,scale)` -(op_t::attn_dec) → `tl_attn_decode_*`. Flash-attention online-softmax in one pass -(the ctx-long scores never touch memory), one block per head with warp-shuffle -score reduction, split-KV across `gridDim.y` + a combine pass so grid = H×S fills -the SMs. KV-bandwidth-bound (~844 GB/s, ~19× the unfused 3-op array path); -numerically exact (maxrel ~1e-7); CPU reference fallback. Integrated end to end: -decode 12.2 → 26.1 tok/s (F32), 41.2 with bf16 weights. - -**Done (KV cache + GQA/MQA, 2026-07-08):** decode-side, kernel + direct-bench -proven (STEEL). A persistent device-resident **`cuda::kv_cache`** object (K,V -buffers `[n_kv_heads, max_ctx, D]` + running `pos`) lives *outside* the lazy graph -— decode is inference-only and stateful, which the immutable node model doesn't -fit (chosen A-surface; B = first-class array node was rejected as high-friction). -Two surgical kernel changes: (1) attention gains a **`kv_stride`** param (head -stride = `max_ctx*D`, distinct from the valid length `ctx`) so it reads a -persistent cache as its prefix `[0,pos)` — `kv_stride==ctx*D` degenerates to the -old whole-buffer path (the array `attn_decode` still rides that, unchanged); (2) a -**`group`** param (`n_q_heads/n_kv_heads`) maps each query head to its shared kv -head for **GQA/MQA** (`group==1` = MHA). Plus `tl_kv_append` — a scatter that -writes one token's projected k,v `[n_kv_heads,D]` into the cache at `pos`. -Verified in `bench_attn_decode` by growing the cache one token at a time to 2048 -and checking a from-scratch CPU reference (prefix + GQA mapping) at 11 checkpoints -straddling the split-KV boundary — **maxrel 1.6e-7**. GQA (32 q / 8 kv) at -ctx=2048 is ~1.85× faster wall-time than MHA (0.044 vs 0.081 ms/layer) and 4× less -KV cache VRAM. - -**Done (causal prefill + cache fill, 2026-07-08):** the prompt-processing regime, -kernel + direct-bench (STEEL), correctness-first baseline. `tl_attn_prefill_f32` -processes all T query positions at once: query `p` attends keys `0..p` (**causal -masking**), one block per (head, query pos) so grid = H_q×T fills the SMs (no -split-KV needed). It **reuses the decode online-softmax verbatim** — only the key -loop cap (`i<=p`) and the query-row index differ — so no T×T scores are ever -materialized. `tl_kv_fill` bulk-copies a prompt's k,v `[H_kv,T,D]` into the cache -`[0,T)`; `kv_cache::prefill(q,k,v,out,T,...)` chains fill → causal attend and -leaves `pos=T` so decode continues. Verified against a from-scratch causal CPU -reference (T=512): **prefill maxrel 1.6e-7**, and a follow-on decode step at -pos=512 (prefill-filled rows + one appended token) **maxrel 8.5e-8** — the -prefill→decode handoff is exact. Baseline throughput ~335 Ktok/s (T=512, one -layer's attention). **This is deliberately un-tiled** — each query re-streams its -keys from DRAM (O(T²) traffic), so it's bandwidth-bound, not compute-optimal. -**Remaining:** a query×key **tiled flash-attention** prefill (shared-mem K/V -reuse across a query block) is the tuning pass, gated on prefill proving a -bottleneck (scope: prefill is a minority one-time cost); T>65535 needs gridDim.y -chunking; the array/culebra surface for the whole cache (currently a `cuda::` -object only); a per-group shared-kv-head block is the GQA tuning lever; a bf16 KV -cache (~2× the attention floor) is the next decode lever. head_dim is still fixed -at 128 (host-gated). See performance-notes.md. - -**Done (model building-block ops + verified decoder block, 2026-07-09):** the -"turn kernels into a runnable LLM" wiring, at the **pure array-op surface** (so it -rides the tuned kernels and is autograd-ready for the fine-tuning path). The one -genuinely new op is **RoPE** (`array::rope(x,pos,base)` → `op_t::rope` → -`tl_rope`, half-split / GPT-NeoX-HF convention, GPU kernel + CPU ref; x is -`[H,D]` decode / `[H,T,D]` prefill, per-row position `pos+(row%T)`). **RMSNorm, -SiLU, SwiGLU are pure compositions** of existing ops (`rmsnorm = x·rsqrt(mean(x²) -+eps)·w`; `silu = x·σ(x)`; `swiglu = silu(gate)·up`) — no new kernels, correct by -construction. A full **llama decoder block** (RMSNorm → qkv proj → RoPE(q,k) → -append → attn_decode → out proj → residual → RMSNorm → SwiGLU MLP → residual) is -assembled from these + `dot` + `attn_decode` and verified vs a from-scratch CPU -reference in `bench/shared/check/check_llm_block.cpp` (ctest `llm_block`, backend-agnostic): -**decoder-block maxrel 1.8e-7** on GPU, all six unit+block checks green on both -CUDA and CPU builds. - -**Done (multi-layer decode loop wired to the cache, 2026-07-09):** the "runnable -LLM" end to end. Added the array↔native bridge — a public **`array::native()`** -(the evaluated device-buffer handle) — so each layer's attention step hands its -`q/k/v` arrays' buffers to the stateful `cuda::kv_cache` (`append` then `attn`, -GQA-aware) and wraps the result back into the array graph, while RMSNorm/RoPE/ -SwiGLU/proj stay pure array ops. `bench/cuda/check/check_llm_decode.cpp` (ctest `llm_decode`) -drives a **3-layer, GQA (4 q / 2 kv) llama** greedily for 12 steps (6-token prompt -+ 6 generated, integer ids over a random embedding table — no tokenizer yet) and -checks every step vs a from-scratch CPU reference: **logits maxrel 1.6e-6, greedy -token sequences match exactly.** ctest 7/7. Lifetime/coherence rest on null-stream -ordering (append/attn queue before the buffers recycle) + the device mirror (the -cache's `device_write_` marks its output live, the next `dot` reads it device-side). -The prompt is fed token-by-token through the decode path (functionally identical to -causal prefill for the cache/logits; the parallel prefill kernel is a perf -optimization). **Remaining to an actual chat:** GGUF weight loading + a BPE -tokenizer + sampling (mechanical bulk — Fable-appropriate); using `attn_prefill` -for the prompt. VJPs for RoPE (training) are deferred; head_dim still fixed at 128. - -**Done (real-model chat — Qwen2.5-0.5B-Instruct end to end, 2026-07-10):** the -project now *actually chats* on a real GGUF, own kernels only, zero third-party -runtime deps. Four pieces landed: (1) **head_dim {64,128} generalization** — the -five attention/KV kernels (`tl_attn_decode`/`split`/`prefill`, `tl_kv_append`/ -`fill`) are now `template` on head_dim (NW=AD/32 warps replacing the -hard-coded 4; the copy/combine kernels read `blockDim.x`), instantiated for -{64,128} with the launchers dispatching by D; `ctest attn64` validates the D=64 -path vs the same CPU refs (maxrel 1e-7). Qwen2's head_dim is 64 (896/14), so this -was a hard prereq. (2) **GGUF v3 loader** (`include/gguf.h`, Fable) — mmap + -typed-KV + tensor directory, validated vs gguf-py (291 tensors match). (3) **Qwen2 -GPT-2 byte-level BPE tokenizer** (`include/tokenizer.h` + generated unicode tables -in `tokenizer_data.h`, Fable) — token-exact vs HF on 7 cases incl. the chat -template. (4) **the model forward** (`bench/cuda/qwen_model.h`) — F16→F32 widen, GGML -`[out,in]`→our `[K,N]` transpose at load, QKV bias, RoPE base=1e6 half-split, -RMSNorm eps=1e-6, GQA 14q/2kv, SwiGLU, composed from the array ops + the D=64 -`kv_cache`. Verified two ways: `check_qwen` matches a **numpy-on-the-same-F16- -weights** oracle (greedy tokens exact, intermediates ~1e-5); `chat_qwen` generates -coherent text whose greedy sequence is **character-for-character identical to -llama.cpp** (F16, same GPU) — the target-scope "match llama.cpp greedy" bar, met -on real hardware. **Speed finding (the census that matters):** decode is **67 -tok/s vs llama.cpp's 446** (6.7× gap) — and bf16 storage gives **65 tok/s (no -speedup)**, proving 0.5B decode is **per-op-overhead-bound, not bandwidth-bound** -(array-graph launch/eval/D2H per op × 24 layers dominates; llama.cpp wins via -fused kernels + CUDA-graph capture, not arithmetic). Also **the WSL2 sysmem cliff -did not fire** at 2GB F32 resident (decode stayed 67 tok/s) — confirming resident -weights are churn-free-safe. So the next decode lever is **per-op overhead** -(GPU-side argmax to kill the per-token 151936-float logits D2H, fused decode step, -graph reuse), NOT bf16/q4 — those pay off only on bandwidth-bound larger models. -Sampling (temp/top-k/top-p) and Q4_0→`dtype::q4` (needs the q4 GEMV's `K%256` -relaxed to `K%32` — Qwen hidden 896 isn't a 256 multiple) remain. - -**Done (decode-overhead attack — census + sync-free step, 2026-07-10):** a -per-region census (`bench/cuda/speed/bench_qwen_decode.cpp`, `StepProf`) split the 15.9 -ms/token as **host array-graph construction 35–45%**, per-layer eval (launch + -sync) ~58%, greedy last-mile only **1.7%** — and established the **~1 ms/token -GPU compute floor** (so ~94% is overhead, and a large slice is *pure host CPU* -building ~720 nodes/token). Two levers, each `check_qwen`-guarded (greedy -bit-identical) and measured **min-of-rounds** (WSL2 boost noise is ±20 tok/s -single-run): (1) **GPU argmax** (`tl_argmax`/`cuda::argmax`, 4-byte index D2H, -host-matching tie-break) — measured **neutral** (last-mile was 1.7%; in isolation -it even adds a 2nd sync) but kept as the **CUDA-graph prerequisite**; (2) **sync- -free `realize()`** (`array::realize()` / `graph::run_noflush`) launches a graph -and adopts its storage **without** the terminal `CtxSynchronize`, so null-stream -ordering feeds the kv_cache/argmax kernels — collapsing a step's ~98 syncs to ~1 -and lifting decode **66→77 tok/s (+16%)**. **Done (fused imperative decode step, -2026-07-10):** `qwenmodel::step_imperative` runs the whole forward as **direct -`cuda::` calls on a persistent `Scratch` arena** (buffers allocated once, -residual ping-ponged) — zero array nodes (kills the ~45% host construct), zero -per-step alloc, one sync/step. Two fused kernels (`tl_rmsnorm`, `tl_swiglu`) -replace the array compositions; gemv/rope/attn/argmax reused. **Decode 77 → ~200 -tok/s (+156%)**, `chat_qwen` end-to-end **208 tok/s**, greedy bit-identical -(canonical answer unchanged; array-vs-imperative 0/12 mismatch on same-cache -replay). **Net session 67 → 208 tok/s (3.1×); gap to llama.cpp 6.7× → 2.14×.** -**Done (CUDA-graph capture POC — ceiling measured, 2026-07-10):** built the infra -(a `context.stream` all launchers target, dlopen'd graph symbols, split-K gemv's -blocking `MemsetD8`→`MemsetD8Async` so the region is capturable) and captured the -imperative forward at fixed pos. Replay + argmax = **221.6 vs imperative 201.6 -tok/s — only +10%.** So decode is now **GPU-kernel-bound, not host-launch-bound**: -the graph erases the host submit sliver but the GPU still runs **456 tiny -kernels/token** (~4.4 ms of per-kernel latency at batch-1). **The remaining ~2× -to llama.cpp is kernel efficiency (fewer/fused/bigger kernels), not launch -overhead** — graph capture's payoff shrinks as you fuse. Measured this ceiling -*before* the risky `pos`-as-device-scalar plumbing (touches the 1e-7 attention -kernels) that correct per-token replay needs; capture infra kept (tested, -dormant). **Next real decode lever: kernel fusion + tuning** (fuse -norm-into-gemv, wider gemv tiles, fused attention) — a larger separate effort; -CUDA-graph capture is a ~10% finisher on top. **Done (elementwise fusion — -neutral, 2026-07-10):** fused the 2 residual adds into the following RMSNorms -(`tl_add_rmsnorm`) and q/k bias into RoPE (`tl_rope` optional bias), −4 -launches/layer (456→360 kernels/token). Imperative decode **unchanged (~200 -tok/s)** — proving decode is **GEMV-execution-bound, not launch/elementwise- -bound**: the 7 gemvs/layer reading the weights run at ~220 GB/s effective vs -~936 GB/s isolated (M=1 underfills the SMs). llama.cpp ≈ 440 GB/s effective, so -its M=1 `mul_mat_vec` is ~2× more efficient — **that gemv-kernel gap is the -entire remaining ~2×.** Fusion kept (correct, cleaner). **Banked at 208 tok/s; -the M=1 dequant-gemv efficiency milestone (fused QKV/gate-up + tuned mul_mat_vec) -is the remaining path to llama.cpp parity.** - -Flash-attention-style fused attention (tiled QKᵀ → online softmax → ·V, never -materializing the S×S scores), causal masking, and a **KV cache** (append per -decode step, GQA/MQA head sharing). cuBLAS does not do attention at all, so this -is squarely own-kernel territory (FlashAttention is itself hand-written CUDA). -Promotes the deferred "GPU fused kernels" item. Decode attention is KV-cache- -bandwidth-bound; prefill attention is compute-bound and long-context-sensitive. -This is what turns the M6/M7/M8 GEMM kernels into an actual runnable LLM -(alongside the model layer — RoPE, norm, MLP — built in culebra or an app layer). - -## Deferred / conditional work - -Tracked so it isn't re-discovered. Each is gated on a trigger, not scheduled. - -| Item | Trigger to act | -|---|---| -| **SteelLoaderT** — transposed-operand STEEL for backward matmuls (`xᵀ@g`, `g@Wᵀ`) | a real GPU workload whose backward is transpose-bound (census-driven) | -| **GPU fused kernels** — `linear_sigmoid`, online-softmax (silarray has these; tensorlib's GPU softmax/MLP trail silarray) | **promoted → M9** (fused attention is the flagship case); general fused-epilogue kernels ride the same seam | -| **Tiny-tensor corner** (n_embd≈16) — node pooling / small-vector fields | a real workload that lives in the ~16-element regime (documented as accepted, not a bug) | -| **View-crossing fusion** — eager view eval breaks fusion across reshape/transpose | a transformer workload where this shows up in a census | -| **Problem-size-aware CPU thread count** — pool statically uses all `hardware_concurrency()` threads; over-parallelizes mid-size GEMM (512³/1024³ < OpenBLAS on hybrid P/E+HT), where OpenBLAS caps threads | a native-Linux x86 measurement session (WSL2 hides P/E topology); or a workload bound by mid-size CPU GEMM | -| **culebra default → `auto`** — currently `cpu` | large-tensor workloads appear, or after CUDA lands and thresholds are re-measured | -| **`reshape` on non-contiguous**, **other-axis slice/concat** | a workload / language feature that needs it (tensorlib can already materialize) | - -## CI / infrastructure - -Current CI (GitHub Actions, free tier): 3 OSes × cpu/gpu/auto modes; -`macos-15` runs real Metal on its paravirtual GPU; ubuntu/windows run the CPU -path + a CUDA build check + the driver-absent fallback. **No GPU numerical CI** -— self-hosted and paid GPU runners were declined (cost/ceremony); CUDA numeric -regression is manual on the RTX 3090 box. To add: a binary-size job, and (once -M5/M6 land) a `just` recipe for the manual CUDA/x86 regression run, and a -Tensor-unused-binary size check on the culebra side. - -## Open decisions - -- **CUDA performance bar — RESOLVED (2026-07-04).** The bar is **tokens/sec vs - llama.cpp/exllamav2 on a real quantized model**, not cuBLAS-90% on square GEMM. - Rationale: the target scope (consumer local-LLM inference/fine-tuning) is - decode-dominated = bandwidth-bound, where cuBLAS has no moat and hand-written - is SOTA. The M6 cuBLAS gate is kept only as the prefill/compute-bound - foundation (~0.82, split-K). Full regime analysis in performance-notes.md. -- **CUTLASS escape hatch — RESOLVED (2026-07-04): pre-authorized but unused.** - Not needed for the scope (see M6 Resolved). If ever reopened, it is confined - to the build-time kernel `.cu` behind `-DTENSORLIB_CUDA_CUTLASS` (off by - default) so non-matmul / non-opted-in binaries carry zero CUTLASS bytes and - the header-only consumer surface is preserved. -- **culebra default device** — stays `cpu` for now; revisit per the deferred - table. Rationale in performance-notes.md. -- **auto thresholds are M1-Pro-specific** — must be re-measured per Mac and a - full set produced once CUDA lands. -- **Windows-native environment** — not set up; Windows is currently - build-checked via CI only. -- **MSVC support vs header-only (AVX2 on Windows)** — decided (2026-07-03). - header-only is a hard invariant and consumers are MSVC-majority, so the - tension is: MSVC has no per-function target attribute and `/arch:AVX2` is - TU-global, so header-only runtime-dispatched AVX2 (what GCC/Clang get via - `__attribute__((target)))`) cannot be replicated on MSVC. Resolution — a - three-tier model, dropping no one: - 1. *Default, header-only:* GCC/Clang do runtime dispatch (already built). - MSVC compiles clean and runs the **scalar** kernel today — the x86 SIMD - path is gated on `__x86_64__`, which MSVC (`_M_X64`) doesn't define, so - nothing MSVC-hostile is even reached. Correct, just not accelerated. - 2. *MSVC AVX2, opt-in, header-only preserved:* build with `/arch:AVX2` → - `select_ukernel` picks AVX2 at compile time (`#if defined(__AVX2__)`). - Caveat: `/arch:AVX2` is TU-wide, so such a binary `#UD`s on pre-AVX2 - CPUs with no runtime fallback — fine for source-built-for-a-known-target - or an AVX2 baseline, not for an adaptive redistributable. - 3. *MSVC runtime dispatch (adaptive redistributable), opt-in:* link a small - AVX2-only **source** file (`cpu_avx2.cpp`, not yet written) compiled - per-file with `/arch:AVX2`; baseline TUs stay AVX2-free; CPUID picks via - the existing `ukernel_fn` pointer (the TU boundary is the isolation the - target attribute gives on GCC/Clang). Ship **source, never a prebuilt - `.obj`** — an obj would need a {toolset}×{CRT: /MT,/MTd,/MD,/MDd}×{arch} - matrix (LNK2038 on CRT mismatch), can't cover future toolsets, and a - binary blob violates the zero-dep / build-from-source / auditable ethos. - A source file auto-matches the consumer's toolchain and stays auditable. - Small mechanical seams alongside this: `__builtin_prefetch` → `_mm_prefetch` - and `__builtin_cpu_supports` → `__cpuid`/`__cpuidex` under `_MSC_VER`. None - of this blocks the gates — WSL2 unblocks all x86 M5 + M6 CUDA; native-Windows - AVX2 is a scoped follow-up. - -## Notes on effort - -The remaining milestones split the same way M3b and M6 do: **scaffolding** -(dispatch seams, loaders, packing loops, build/CMake) is mechanical and -follows documented designs; the **kernels themselves** — reaching the -OpenBLAS/cuBLAS-90% and beat-PyTorch gates — are performance-sensitive sprints -where the STEEL-sprint discipline applies (per-kernel census, interleaved A/B, -record refuted approaches). Plan each milestone as *scaffold to correct, then -a dedicated tuning pass to the gate.* diff --git a/docs/x86-validation-runbook.md b/docs/x86-validation-runbook.md deleted file mode 100644 index eb94b49..0000000 --- a/docs/x86-validation-runbook.md +++ /dev/null @@ -1,146 +0,0 @@ -# x86 validation runbook (M5 AVX2 — WSL2 / Linux) - -Step-by-step for taking the M5 CPU backend from "compile-checked on the Mac" -to "executed, numerically verified, and tuned to the OpenBLAS-90% gate" on a -native x86 box (the RTX 3090 Ubuntu-on-WSL2 machine, or any x86 Linux). The -AVX2 microkernel and the CPUID dispatch were written and disassembly-checked -on Apple but **never executed** there — Rosetta stops at SSE4.2 — so this is -where they run for real. Companion to [roadmap.md](roadmap.md) (M5 scope + -the MSVC/header-only decision) and [performance-notes.md](performance-notes.md) -(measurement discipline, the NEON census this mirrors). - -All commands assume the repo root as the working directory. The repo is -portable and travels via git, so on the box: clone, then run these. - -## 0. One-time environment - -```sh -sudo apt update -sudo apt install -y build-essential cmake pkg-config libopenblas-dev -# Compiler: any g++ >= 11 or clang++ >= 13 (whatever accepts -std=c++2b) — -# the code is C++17-level compiled under the C++23 flag, so it isn't picky; -# Ubuntu 22.04's default g++-11 already works. CMake >= 3.20 is needed for -# CMAKE_CXX_STANDARD 23. libopenblas-dev is the gate baseline; pkg-config -# resolves its include/lib paths portably. Substitute clang++ for g++ in the -# bench commands below if you prefer — both are tested (the Mac builds clang). -``` - -Sanity-check the CPU actually has the ISA the kernel dispatches to: - -```sh -lscpu | grep -iE 'model name|^flags' | grep -oiE 'avx2|fma|avx512[a-z]*|model name.*' -nproc # thread count the pool will use -``` - -If `avx2` and `fma` are absent, `select_ukernel` will fall to scalar — expected -on very old CPUs, but then there is nothing to validate here. `avx512f` present -is a bonus for the deferred AVX-512 work (not built yet). - -## 1. Build + run the test suite (native Linux) - -CMake needs no Linux-specific flags — the Accelerate/Metal frameworks are added -only `if(APPLE)`, and metal.h compiles to non-Apple stubs (its `#embed` is -inside the `__APPLE__` block, so no C23 `#embed` is needed on Linux). - -```sh -cmake -S . -B build-linux -DCMAKE_BUILD_TYPE=Release -cmake --build build-linux -j"$(nproc)" -cd build-linux && ctest --output-on-failure && cd .. -``` - -Expect cpu/gpu/auto all pass (gpu/auto degrade to the CPU/ref path with no -device). This confirms the refactored function-pointer dispatch is correct on -x86 with AVX2 actually selected (not just the ARM/NEON path the Mac tested). - -## 2. AVX2 numerical validation — the first real execution - -`bench/cpu/check/check_cpu_ukernel.cpp` compares each available microkernel against a -naive triple loop over full and edge tiles. On x86 this is the **first time the -AVX2 kernel runs and is checked numerically** (on the Mac it only ran scalar + -NEON). - -```sh -g++ -std=c++2b -O2 -I include bench/cpu/check/check_cpu_ukernel.cpp -o check_uk -pthread -./check_uk # expect: scalar / avx2-8x8 / avx2-6x16: ok / ALL OK -``` - -Done 2026-07-03 on the i7-12700KF: all three ok. See performance-notes.md -"x86 census" for the full result; the notes below are the record of how it -was reached (kept as the procedure for the next box / AVX-512). - -`-pthread` is required on Linux (the thread pool uses `std::thread`). If -`avx2: ok` does not print, confirm the CPU has AVX2+FMA (step 0) — the harness -only tests kernels the build compiled in (`#ifdef TL_CPU_X86`). - -## 3. GEMM benchmark vs ref vs OpenBLAS — the gate - -Linux build of the dispatch-path bench: drop the Mac `-framework` flags, add -`-pthread`, and link OpenBLAS via pkg-config for the gate. - -```sh -# own vs ref only: -g++ -std=c++2b -O2 -I include bench/cpu/speed/bench_cpu_gemm.cpp -o bcg -pthread -./bcg - -# with the OpenBLAS gate (threads matched inside the harness): -g++ -std=c++2b -O2 -I include -DBENCH_HAS_OPENBLAS $(pkg-config --cflags openblas) \ - bench/cpu/speed/bench_cpu_gemm.cpp -o bcg_gate -pthread $(pkg-config --libs openblas) -./bcg_gate -``` - -If `pkg-config openblas` is missing on the distro, substitute the explicit -paths (`-I/usr/include/x86_64-linux-gnu -lopenblas`, adjust per distro). Unlike -the Mac — where Homebrew's OpenBLAS was an untuned ARMv8 build and only a loose -lower bound — a properly built x86 OpenBLAS here is the **real** gate baseline. - -**Gate: own ≥ 90% of OpenBLAS** (the `own/blas %` column), single- and -multi-threaded, and ≥ PyTorch CPU. This is the verdict the Mac could not give. - -## 4. Re-tune block sizes for this microarchitecture - -The MC/KC/NC and the 8×8 tile were tuned on M1 Pro NEON; the x86 cache -hierarchy and the AVX2 register file differ, so re-sweep. `bench_cpu_sweep.cpp` -calls `cpu::sgemm` directly (≈1s compile) with `-D` overrides: - -```sh -for kc in 256 384 512 768; do - g++ -std=c++2b -O2 -I include -DTL_CPU_KC=$kc bench/cpu/speed/bench_cpu_sweep.cpp \ - -o sw_$kc -pthread -done -# interleave, take per-size medians (see below), pick the KC; repeat for MC. -for i in 1 2 3 4; do for kc in 256 384 512 768; do ./sw_$kc; done; done -``` - -A/B settled 2026-07-03: **6×16 beat 8×8** by +18% single-core (91% vs 77% of -AVX2 peak) and is now the x86 default; 8×8 is kept behind `-DTL_CPU_AVX2_8X8`. -The 6×16 kernel added the NR=16 pack variant (via `ukernel_desc`), which the -AVX-512 kernel will reuse. KC=512 re-confirmed. On a **new** µarch, re-run this -A/B and the KC sweep — the cache hierarchy and register pressure differ. - -## 5. Measurement discipline (from performance-notes.md) - -- Check `uptime` first — the load-1 average should be low (< ~nproc/2) or the - numbers are noise. On WSL2, also close Windows-side CPU hogs. -- **Interleave** A/B builds and take **medians**, never a single run. -- Record which side of the WSL2 boundary a number came from (WSL2 adds - submission latency vs native Linux). -- Track **% of AVX2 fp32 peak** as the architecture-independent yardstick: - `peak_GFLOPs = cores × freq_GHz × 2 (FMA units) × 8 (fp32/ymm) × 2 (flops/FMA)`. - Fill in the box's `cores`/`freq` and compare the 2048³ number against it, the - way the NEON census tracked ~86–89% of the M1 Pro's ~410 GFLOP/s. - -## 6. Recording results - -Add an "Own CPU GEMM x86 census" section to performance-notes.md mirroring the -NEON one: the own/ref/OpenBLAS table, % of AVX2 peak, the chosen MC/KC/NC, the -8×8-vs-6×16 verdict, and any refuted approaches. Update roadmap M5: flip the -"AVX2 execution + gate" remaining item and note the OpenBLAS-90% verdict. - -## Deferred (still not this pass) - -- **AVX-512 kernel** — wants NR=16 (a second packing layout). The dispatch seam - is ready: add an `avx512f` branch to `select_ukernel` and the wider ukernel. - Only worth it if the box's CPU has AVX-512 and the AVX2 gate is met first. -- **Native Windows / MSVC AVX2** — see the roadmap Open-decisions entry - (three-tier: header-only default, `/arch:AVX2` opt-in, `cpu_avx2.cpp` source - TU for runtime dispatch). Independent of this Linux pass. From 96f8b0e8d7d1ca5d4c81ec4ebf27259d2402e34c Mon Sep 17 00:00:00 2001 From: yhirose Date: Sun, 19 Jul 2026 23:41:53 -0400 Subject: [PATCH 03/12] =?UTF-8?q?spike:=20WebGPU=20backend=20probe=20?= =?UTF-8?q?=E2=80=94=20sync=20flush()=20holds,=20WGSL=20sgemm=20~75x=20CPU?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Standalone spike under spike/webgpu/, not wired into include/. Answers the two questions that would otherwise have been guessed at. The async boundary is not the obstacle it looked like. Device acquisition needs no async on the C++ side (JS passes a ready device via Module.preinitializedWebGPUDevice), and blocking waits are supported upstream via wgpuInstanceWaitAny(timeout>0) with the TimedWaitAny instance feature. gpu::flush() and sync_to_host() keep their signatures; the suspend surface is two call sites, not the interpreter. A 64x64-tile WGSL sgemm (the shape of sgemm_64_, minus MMA intrinsics WGSL lacks) runs ~75x the scalar wasm CPU path at n=1024, 580-630 GF/s on an M1 Pro, matching the CPU oracle to 5e-06. Two findings that change the design: dispatch+wait has a ~0.6-1.0ms fixed floor, so auto_threshold_() needs its own arm rather than Metal's; and Asyncify vs JSPI is a build-shape fork (JSPI is smaller and avoids the -fwasm-exceptions conflict but is Chrome-only, so it means two builds). See spike/webgpu/README.md for measurements and what is still uncovered. --- spike/webgpu/.gitignore | 2 + spike/webgpu/README.md | 88 ++++++++++ spike/webgpu/build.sh | 53 ++++++ spike/webgpu/index.html | 42 +++++ spike/webgpu/sgemm.wgsl | 96 ++++++++++ spike/webgpu/spike_main.cpp | 340 ++++++++++++++++++++++++++++++++++++ spike/webgpu/worker.js | 44 +++++ 7 files changed, 665 insertions(+) create mode 100644 spike/webgpu/.gitignore create mode 100644 spike/webgpu/README.md create mode 100755 spike/webgpu/build.sh create mode 100644 spike/webgpu/index.html create mode 100644 spike/webgpu/sgemm.wgsl create mode 100644 spike/webgpu/spike_main.cpp create mode 100644 spike/webgpu/worker.js diff --git a/spike/webgpu/.gitignore b/spike/webgpu/.gitignore new file mode 100644 index 0000000..dc5e40a --- /dev/null +++ b/spike/webgpu/.gitignore @@ -0,0 +1,2 @@ +site/ +sgemm_wgsl.inc diff --git a/spike/webgpu/README.md b/spike/webgpu/README.md new file mode 100644 index 0000000..308de5d --- /dev/null +++ b/spike/webgpu/README.md @@ -0,0 +1,88 @@ +# WebGPU backend spike + +A standalone probe, not a backend. Nothing here is included by `include/`; the +point was to measure the things that would otherwise be guessed at before +committing to a design. + +Build and run: + +```sh +./build.sh jspi # or: asyncify +python3 -m http.server 8731 +# open http://localhost:8731/site/index.html +``` + +Measured on an M1 Pro (`apple / metal-3`), Chrome, emsdk 6.0.3, +emdawnwebgpu `v20260423.175430`. + +## What it proves + +**1. `gpu::flush()` can stay synchronous.** This was the stated hard part and it +turns out to be a solved problem upstream rather than something we have to +engineer around: + +- Device acquisition needs no async on the C++ side at all. JS does + `requestAdapter`/`requestDevice` and passes the result as + `Module.preinitializedWebGPUDevice`; C++ picks it up with a plain + `emscripten_webgpu_get_device()`. +- Blocking waits are a supported, documented path: + `wgpuInstanceWaitAny(timeoutNS > 0)` after creating the instance with the + `TimedWaitAny` feature. Internally it suspends via Asyncify/JSPI + (`webgpu.cpp:560`, "To handle timeouts, use Asyncify and proxy back into + JS"). If the link lacks Asyncify/JSPI, `CreateInstance` returns null — it + fails loudly at startup rather than deadlocking later. + +So `flush()` and `sync_to_host()` keep their existing signatures. The suspend +surface is exactly two call sites (`OnSubmittedWorkDone`, `MapAsync`), not the +whole interpreter. Uploads (`WriteBuffer`) are queued, not awaited, and cost +nothing here. + +**2. A hand-written WGSL sgemm is clearly worth it.** 64x64 threadgroup tile, +16x16 invocations, 4x4 register accumulator — the shape of `sgemm_64_` in +`metal_kernels.metal`, minus the MMA intrinsics, which WGSL has no equivalent +for. Medians of 7 (GPU) / 3 (CPU), interleaved, after warmup: + +| n | GPU | CPU (tensorlib, scalar wasm) | speedup | maxerr | +|------|------------|------------------------------|---------|---------| +| 128 | 0.6–1.0 ms (4–7 GF/s) | 0.5 ms (8.4 GF/s) | ~1x | 4.8e-07 | +| 256 | 1.0–1.1 ms (30–34 GF/s) | 4.2 ms (8.0 GF/s) | ~4x | 7.8e-07 | +| 512 | 1.7 ms (158 GF/s) | 33–37 ms (7–8 GF/s)| ~20x | 1.4e-06 | +| 1024 | 3.4–3.7 ms (580–630 GF/s)| 263–271 ms (~8 GF/s)| ~75x | 4.8e-06 | + +## What it changes about the plan + +**There is a ~0.6–1.0 ms fixed floor per dispatch+wait.** Metal's is orders of +magnitude lower. Two consequences: + +- `auto_threshold_()` cannot reuse the Metal arm. Crossover here is around + n≈128–160, i.e. work ≈ 4–8e6 in the `num_elements * k` units `gpu_gemm` uses, + against Metal's 5e8. It is *lower* than Metal's, not higher, because the wasm + CPU baseline is scalar and weak (~8 GF/s, no threads, no SIMD) — not because + WebGPU is fast. Port `misc/census.cpp` to the browser before picking numbers. +- Batching matters more than on Metal. One `flush()` per graph eval is fine; + one per op would be dominated by the floor. + +**Asyncify vs JSPI is a real fork, and it is not a performance question.** Both +work and both measure the same (an early run suggesting otherwise was a +cold-start outlier — worth repeating, given how easily this one misleads): + +| | Asyncify | JSPI | +|---------------------|----------|------| +| spike wasm size | 397 KB | 258 KB | +| `-fwasm-exceptions` | emcc warns they are incompatible; the spike links anyway, but it does not throw across a suspend — culebra's interpreter might | clean | +| browser support | universal | Chrome shipped; Firefox behind a flag; Safari not yet | +| extra link flags | — | `-sJSPI_EXPORTS=`, and JS must `await` those calls | + +The catch: `-sJSPI` produces a wasm that *requires* JSPI, so it will not run at +all on Safari. Shipping GPU via JSPI therefore means two builds plus a JS +feature check, not one build that degrades. Asyncify is one universal build but +pays size everywhere and has an unresolved interaction with the exception model +culebra already uses. + +## Not covered + +Only `gemm`, only F32, only the untransposed case, no buffer pooling, no +`storage::native` integration. The host-visible-pointer assumption in +`gpu::alloc(bytes, float** contents)` is the next thing to design against: +WebGPU storage buffers have no CPU pointer, so the backend needs the +host-shadow + `sync_to_host()` mirror that `cuda.h` already implements. diff --git a/spike/webgpu/build.sh b/spike/webgpu/build.sh new file mode 100755 index 0000000..4c494ac --- /dev/null +++ b/spike/webgpu/build.sh @@ -0,0 +1,53 @@ +#!/usr/bin/env bash +# Build the WebGPU spike. Usage: ./build.sh [asyncify|jspi|asyncify-only] +set -euo pipefail + +MODE="${1:-asyncify}" +HERE="$(cd "$(dirname "$0")" && pwd)" +ROOT="$(cd "$HERE/../.." && pwd)" +OUT="$HERE/site" + +source "${EMSDK_DIR:-$HOME/Projects/emsdk}/emsdk_env.sh" >/dev/null 2>&1 + +mkdir -p "$OUT" + +# WGSL -> C string literal. Mirrors the bin2c step the CUDA backend uses for +# PTX; #embed is macOS/C23-only in practice and this has to build under wasm. +{ + echo '// Generated by build.sh -- do not edit.' + echo 'static const char* kSgemmWgsl = R"WGSL(' + cat "$HERE/sgemm.wgsl" + echo ')WGSL";' +} > "$HERE/sgemm_wgsl.inc" + +case "$MODE" in + asyncify) ASYNC_FLAGS=(-sASYNCIFY=1 -sASYNCIFY_STACK_SIZE=65536) ;; + # Restrict instrumentation to the WaitAny call chain. If this list is + # sufficient, the cost of a synchronous flush() is bounded and does not scale + # with the size of the interpreter linked alongside it. + asyncify-only) ASYNC_FLAGS=(-sASYNCIFY=1 -sASYNCIFY_STACK_SIZE=65536 + -sASYNCIFY_ADVISE=1) ;; + # Exports that may suspend must be declared, so the JS side gets a promising + # wrapper. The C++ call chain underneath stays ordinary synchronous code. + jspi) ASYNC_FLAGS=(-sJSPI=1 + -sJSPI_EXPORTS=wgpu_spike_init,wgpu_spike_bench) ;; + *) echo "unknown mode: $MODE" >&2; exit 1 ;; +esac + +set -x +emcc -std=c++23 -O2 -fwasm-exceptions \ + --use-port=emdawnwebgpu \ + "${ASYNC_FLAGS[@]}" \ + -sSTACK_SIZE=16MB -sALLOW_MEMORY_GROWTH=1 -sINITIAL_MEMORY=64MB \ + -sMODULARIZE=1 -sEXPORT_ES6=1 -sEXPORT_NAME=createSpike \ + -sENVIRONMENT=web,worker \ + -sEXPORTED_FUNCTIONS=_wgpu_spike_init,_wgpu_spike_bench,_wgpu_spike_log,_malloc,_free \ + -sEXPORTED_RUNTIME_METHODS=ccall,cwrap,UTF8ToString \ + -I"$ROOT/include" -I"$HERE" \ + "$HERE/spike_main.cpp" \ + -o "$OUT/spike.js" +set +x + +cp "$HERE/index.html" "$HERE/worker.js" "$OUT/" +ls -l "$OUT" +echo "wasm size: $(stat -f%z "$OUT/spike.wasm") bytes (mode=$MODE)" diff --git a/spike/webgpu/index.html b/spike/webgpu/index.html new file mode 100644 index 0000000..9cd3025 --- /dev/null +++ b/spike/webgpu/index.html @@ -0,0 +1,42 @@ + + +cpp-tensorlib WebGPU spike + +

cpp-tensorlib WebGPU spike

+

booting…

+

+

+
diff --git a/spike/webgpu/sgemm.wgsl b/spike/webgpu/sgemm.wgsl
new file mode 100644
index 0000000..f453afd
--- /dev/null
+++ b/spike/webgpu/sgemm.wgsl
@@ -0,0 +1,96 @@
+// Tiled SGEMM: C(M,N) = A(M,K) @ B(K,N), row-major, no transpose.
+//
+// 64x64 threadgroup tile, 16x16 = 256 invocations, each holding a 4x4
+// register accumulator. Mirrors the shape of the Metal sgemm_64_ kernel so the
+// measured numbers are comparable; MMA intrinsics have no WGSL equivalent, so
+// the inner product is plain FMA over registers.
+
+struct Params {
+  M : u32,
+  N : u32,
+  K : u32,
+  scale : f32,
+  offset : f32,
+  _pad0 : u32,
+  _pad1 : u32,
+  _pad2 : u32,
+};
+
+@group(0) @binding(0) var       A : array;
+@group(0) @binding(1) var       B : array;
+@group(0) @binding(2) var C : array;
+@group(0) @binding(3) var             p : Params;
+
+const BM : u32 = 64u;
+const BN : u32 = 64u;
+const BK : u32 = 16u;
+const TM : u32 = 4u;
+const TN : u32 = 4u;
+const THREADS : u32 = 256u;
+
+var As : array;  // BM * BK
+var Bs : array;  // BK * BN
+
+@compute @workgroup_size(16, 16, 1)
+fn sgemm(@builtin(workgroup_id) wg : vec3,
+         @builtin(local_invocation_id) lid : vec3) {
+  let m_base = wg.y * BM;
+  let n_base = wg.x * BN;
+  let tid = lid.y * 16u + lid.x;
+
+  var acc : array;  // TM * TN, zero-initialized
+
+  let n_tiles = (p.K + BK - 1u) / BK;
+  for (var kt : u32 = 0u; kt < n_tiles; kt = kt + 1u) {
+    let k0 = kt * BK;
+
+    // Stage A tile (BM x BK) and B tile (BK x BN), 4 elements per invocation
+    // each. Out-of-range reads are zero-filled so the tail tile needs no
+    // special case in the inner loop.
+    for (var s : u32 = 0u; s < 4u; s = s + 1u) {
+      let i = tid + s * THREADS;
+
+      let ar = i / BK;
+      let ac = i % BK;
+      let am = m_base + ar;
+      let ak = k0 + ac;
+      As[i] = select(0.0, A[am * p.K + ak], am < p.M && ak < p.K);
+
+      let br = i / BN;
+      let bc = i % BN;
+      let bk = k0 + br;
+      let bn = n_base + bc;
+      Bs[i] = select(0.0, B[bk * p.N + bn], bk < p.K && bn < p.N);
+    }
+
+    workgroupBarrier();
+
+    for (var kk : u32 = 0u; kk < BK; kk = kk + 1u) {
+      var av : array;
+      var bv : array;
+      for (var i : u32 = 0u; i < TM; i = i + 1u) {
+        av[i] = As[(lid.y * TM + i) * BK + kk];
+      }
+      for (var j : u32 = 0u; j < TN; j = j + 1u) {
+        bv[j] = Bs[kk * BN + lid.x * TN + j];
+      }
+      for (var i : u32 = 0u; i < TM; i = i + 1u) {
+        for (var j : u32 = 0u; j < TN; j = j + 1u) {
+          acc[i * TN + j] = fma(av[i], bv[j], acc[i * TN + j]);
+        }
+      }
+    }
+
+    workgroupBarrier();
+  }
+
+  for (var i : u32 = 0u; i < TM; i = i + 1u) {
+    let m = m_base + lid.y * TM + i;
+    if (m >= p.M) { continue; }
+    for (var j : u32 = 0u; j < TN; j = j + 1u) {
+      let n = n_base + lid.x * TN + j;
+      if (n >= p.N) { continue; }
+      C[m * p.N + n] = acc[i * TN + j] * p.scale + p.offset;
+    }
+  }
+}
diff --git a/spike/webgpu/spike_main.cpp b/spike/webgpu/spike_main.cpp
new file mode 100644
index 0000000..33b6072
--- /dev/null
+++ b/spike/webgpu/spike_main.cpp
@@ -0,0 +1,340 @@
+// WebGPU backend spike for cpp-tensorlib.
+//
+// What this is proving, in order of how much it would hurt to get wrong:
+//
+//  1. gpu::flush() can keep its synchronous signature under wasm. Metal's
+//     flush() is commit + waitUntilCompleted; WebGPU's queue submission is
+//     fire-and-forget and readback is promise-based. The claim under test is
+//     that wgpuInstanceWaitAny(timeoutNS > 0) restores a blocking wait, and
+//     that the Asyncify instrumentation this needs is confined to the wait
+//     call chain rather than smeared over the whole interpreter.
+//  2. Device acquisition needs no async on the C++ side at all: JS hands us a
+//     ready device via Module.preinitializedWebGPUDevice.
+//  3. A hand-written WGSL sgemm is worth the trouble versus the scalar CPU
+//     path that wasm currently falls back to.
+//
+// Deliberately standalone: nothing here touches include/, so the real backend
+// can be shaped by what this measures rather than the other way round.
+
+#include 
+#include 
+
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+
+#include "sgemm_wgsl.inc"  // const char* kSgemmWgsl, generated by build.sh
+
+// The CPU baseline is tensorlib's own path, not a hand-rolled loop -- an
+// honest speedup number has to be against the code the playground runs today.
+#include 
+
+namespace {
+
+struct Gpu {
+  wgpu::Instance instance;
+  wgpu::Device device;
+  wgpu::Queue queue;
+  wgpu::ComputePipeline pipeline;
+  wgpu::BindGroupLayout bgl;
+  bool ready = false;
+} g;
+
+std::string g_log;
+
+void logf(const char* fmt, ...) {
+  char buf[1024];
+  va_list ap;
+  va_start(ap, fmt);
+  vsnprintf(buf, sizeof(buf), fmt, ap);
+  va_end(ap);
+  g_log += buf;
+  g_log += '\n';
+  printf("%s\n", buf);
+}
+
+// Blocking wait on a single future. This is the whole ballgame: every
+// synchronous-looking call below bottoms out here, and here alone.
+bool wait(wgpu::Future f) {
+  auto st = g.instance.WaitAny(f, UINT64_MAX);
+  return st == wgpu::WaitStatus::Success;
+}
+
+double now_ms() {
+  return std::chrono::duration(
+             std::chrono::steady_clock::now().time_since_epoch())
+      .count();
+}
+
+}  // namespace
+
+extern "C" {
+
+EMSCRIPTEN_KEEPALIVE int wgpu_spike_init() {
+  g_log.clear();
+
+  // TimedWaitAny is what makes WaitAny(timeout > 0) legal. emdawnwebgpu
+  // refuses to create the instance if Asyncify/JSPI is absent, so a link
+  // without -sASYNCIFY fails loudly here instead of deadlocking later.
+  wgpu::InstanceFeatureName features[] = {
+      wgpu::InstanceFeatureName::TimedWaitAny};
+  wgpu::InstanceDescriptor idesc = {};
+  idesc.requiredFeatureCount = 1;
+  idesc.requiredFeatures = features;
+
+  g.instance = wgpu::CreateInstance(&idesc);
+  if (!g.instance) {
+    logf("FAIL: CreateInstance returned null (TimedWaitAny unavailable -- "
+         "was this linked without -sASYNCIFY / -sJSPI?)");
+    return 0;
+  }
+  logf("ok: instance created with TimedWaitAny");
+
+  // No adapter/device request round-trip: JS already did it.
+  g.device = wgpu::Device::Acquire(emscripten_webgpu_get_device());
+  if (!g.device) {
+    logf("FAIL: no preinitialized device (JS must set "
+         "Module.preinitializedWebGPUDevice before calling in)");
+    return 0;
+  }
+  g.queue = g.device.GetQueue();
+  logf("ok: device acquired synchronously from JS");
+
+  wgpu::ShaderSourceWGSL wgsl = {};
+  wgsl.code = kSgemmWgsl;
+  wgpu::ShaderModuleDescriptor smd = {};
+  smd.nextInChain = &wgsl;
+  wgpu::ShaderModule module = g.device.CreateShaderModule(&smd);
+
+  wgpu::ComputePipelineDescriptor pd = {};
+  pd.compute.module = module;
+  pd.compute.entryPoint = "sgemm";
+  g.pipeline = g.device.CreateComputePipeline(&pd);
+  if (!g.pipeline) {
+    logf("FAIL: CreateComputePipeline returned null");
+    return 0;
+  }
+  g.bgl = g.pipeline.GetBindGroupLayout(0);
+  logf("ok: WGSL sgemm compiled");
+
+  g.ready = true;
+  return 1;
+}
+
+}  // extern "C"
+
+// Resources for one problem size, created once and reused across timed
+// iterations. The real backend pools buffers the way metal.h does with its
+// size-keyed free list, so folding allocation into every measurement would
+// charge the GPU for a cost the shipping code does not pay.
+struct Gemm {
+  wgpu::Buffer a, b, c, params, staging;
+  wgpu::BindGroup bg;
+  int M = 0, N = 0, K = 0;
+  uint64_t cbytes = 0;
+};
+
+namespace {
+
+wgpu::Buffer mkbuf(uint64_t size, wgpu::BufferUsage usage) {
+  wgpu::BufferDescriptor d = {};
+  d.size = size;
+  d.usage = usage;
+  return g.device.CreateBuffer(&d);
+}
+
+struct Params {
+  uint32_t M, N, K;
+  float scale, offset;
+  uint32_t pad0, pad1, pad2;
+};
+
+Gemm prepare(const float* a, const float* b, int M, int N, int K) {
+  Gemm r;
+  r.M = M; r.N = N; r.K = K;
+  const uint64_t abytes = uint64_t(M) * K * 4;
+  const uint64_t bbytes = uint64_t(K) * N * 4;
+  r.cbytes = uint64_t(M) * N * 4;
+
+  r.a = mkbuf(abytes, wgpu::BufferUsage::Storage | wgpu::BufferUsage::CopyDst);
+  r.b = mkbuf(bbytes, wgpu::BufferUsage::Storage | wgpu::BufferUsage::CopyDst);
+  r.c = mkbuf(r.cbytes, wgpu::BufferUsage::Storage | wgpu::BufferUsage::CopySrc);
+  r.params = mkbuf(sizeof(Params),
+                   wgpu::BufferUsage::Uniform | wgpu::BufferUsage::CopyDst);
+  r.staging = mkbuf(r.cbytes,
+                    wgpu::BufferUsage::MapRead | wgpu::BufferUsage::CopyDst);
+
+  Params p = {uint32_t(M), uint32_t(N), uint32_t(K), 1.0f, 0.0f, 0, 0, 0};
+  // WriteBuffer is synchronous from the caller's point of view -- the copy is
+  // queued, not awaited. No suspend cost on the upload path.
+  g.queue.WriteBuffer(r.a, 0, a, abytes);
+  g.queue.WriteBuffer(r.b, 0, b, bbytes);
+  g.queue.WriteBuffer(r.params, 0, &p, sizeof(p));
+
+  wgpu::BindGroupEntry e[4] = {};
+  e[0].binding = 0; e[0].buffer = r.a; e[0].size = abytes;
+  e[1].binding = 1; e[1].buffer = r.b; e[1].size = bbytes;
+  e[2].binding = 2; e[2].buffer = r.c; e[2].size = r.cbytes;
+  e[3].binding = 3; e[3].buffer = r.params; e[3].size = sizeof(Params);
+  wgpu::BindGroupDescriptor bgd = {};
+  bgd.layout = g.bgl;
+  bgd.entryCount = 4;
+  bgd.entries = e;
+  r.bg = g.device.CreateBindGroup(&bgd);
+  return r;
+}
+
+// One dispatch plus a blocking wait -- this is exactly the shape gpu::flush()
+// has to keep. Returns elapsed ms.
+double dispatch(const Gemm& r) {
+  double t0 = now_ms();
+  auto enc = g.device.CreateCommandEncoder();
+  auto pass = enc.BeginComputePass();
+  pass.SetPipeline(g.pipeline);
+  pass.SetBindGroup(0, r.bg);
+  pass.DispatchWorkgroups((r.N + 63) / 64, (r.M + 63) / 64, 1);
+  pass.End();
+  auto cmds = enc.Finish();
+  g.queue.Submit(1, &cmds);
+  if (!wait(g.queue.OnSubmittedWorkDone(
+          wgpu::CallbackMode::WaitAnyOnly,
+          [](wgpu::QueueWorkDoneStatus, wgpu::StringView) {}))) {
+    return -1.0;
+  }
+  return now_ms() - t0;
+}
+
+// sync_to_host(): map, copy, unmap -- all behind a synchronous signature.
+bool readback(const Gemm& r, float* out) {
+  auto enc = g.device.CreateCommandEncoder();
+  enc.CopyBufferToBuffer(r.c, 0, r.staging, 0, r.cbytes);
+  auto cmds = enc.Finish();
+  g.queue.Submit(1, &cmds);
+
+  bool ok = false;
+  if (!wait(r.staging.MapAsync(wgpu::MapMode::Read, 0, r.cbytes,
+                               wgpu::CallbackMode::WaitAnyOnly,
+                               [&ok](wgpu::MapAsyncStatus s, wgpu::StringView) {
+                                 ok = (s == wgpu::MapAsyncStatus::Success);
+                               })) ||
+      !ok) {
+    logf("FAIL: MapAsync wait failed");
+    return false;
+  }
+  const void* src = r.staging.GetConstMappedRange(0, r.cbytes);
+  if (!src) {
+    logf("FAIL: GetConstMappedRange returned null");
+    return false;
+  }
+  memcpy(out, src, r.cbytes);
+  r.staging.Unmap();
+  return true;
+}
+
+double median(std::vector v) {
+  std::sort(v.begin(), v.end());
+  return v[v.size() / 2];
+}
+
+}  // namespace
+
+extern "C" {
+
+EMSCRIPTEN_KEEPALIVE const char* wgpu_spike_log() { return g_log.c_str(); }
+
+// Correctness against tensorlib's CPU path, then an interleaved A/B at a few
+// sizes. Interleaved because a straight "all GPU then all CPU" ordering lets
+// thermal drift masquerade as a speedup.
+EMSCRIPTEN_KEEPALIVE int wgpu_spike_bench() {
+  if (!g.ready) return 0;
+
+  auto fill = [](std::vector& v, int mul, int mod) {
+    for (size_t i = 0; i < v.size(); ++i)
+      v[i] = float((i * mul) % mod) / float(mod) - 0.5f;
+  };
+
+  // Warm both paths before recording anything. The first GPU dispatch pays
+  // lazy pipeline setup and the first tensorlib eval pays its own one-time
+  // costs; without this the smallest size absorbs both and reads as the
+  // slowest, which is exactly backwards.
+  {
+    const int w = 128;
+    std::vector a(size_t(w) * w), b(size_t(w) * w), o(size_t(w) * w);
+    fill(a, 37, 101);
+    fill(b, 53, 97);
+    Gemm r = prepare(a.data(), b.data(), w, w, w);
+    dispatch(r);
+    readback(r, o.data());
+    tl::use_cpu();
+    tl::array::from(a, {w, w}).dot(tl::array::from(b, {w, w})).eval();
+    logf("(warmup done)");
+  }
+
+  logf("%-6s %22s %22s %8s %10s", "n", "gpu", "cpu (tensorlib)", "speedup",
+       "maxerr");
+
+  const int sizes[] = {128, 256, 512, 1024};
+  for (int n : sizes) {
+    std::vector a(size_t(n) * n), b(size_t(n) * n), gpu(size_t(n) * n);
+    fill(a, 37, 101);
+    fill(b, 53, 97);
+
+    Gemm r = prepare(a.data(), b.data(), n, n, n);
+
+    // The CPU path is slow enough at n=1024 (tens of seconds, scalar wasm)
+    // that repeating it is not worth the wall clock; the GPU side is cheap so
+    // it gets a proper median.
+    const int greps = 7;
+    const int creps = n >= 1024 ? 1 : 3;
+
+    tl::use_cpu();
+    auto ta = tl::array::from(a, {n, n});
+    auto tb = tl::array::from(b, {n, n});
+
+    std::vector gt, ct;
+    tl::array tc;
+    // Interleaved: alternate GPU and CPU reps so thermal drift or a background
+    // tab hits both arms rather than whichever ran second.
+    for (int i = 0; i < std::max(greps, creps); ++i) {
+      if (i < greps) {
+        double ms = dispatch(r);
+        if (ms < 0) { logf("FAIL: dispatch failed at n=%d", n); return 0; }
+        gt.push_back(ms);
+      }
+      if (i < creps) {
+        double c0 = now_ms();
+        tc = ta.dot(tb).eval();
+        ct.push_back(now_ms() - c0);
+      }
+    }
+
+    if (!readback(r, gpu.data())) return 0;
+
+    // Tolerance is loose-ish: the WGSL kernel accumulates in a different order
+    // than the CPU microkernel, and f32 sgemm drifts accordingly as K grows.
+    const float* ref = tc.raw();
+    double maxerr = 0.0;
+    for (size_t i = 0; i < gpu.size(); ++i) {
+      double e = std::abs(double(gpu[i]) - double(ref[i]));
+      if (e > maxerr) maxerr = e;
+    }
+
+    double gms = median(gt), cms = median(ct);
+    double gflop = 2.0 * n * n * n / 1e9;
+    logf("n=%-5d %8.2f ms %8.1f GF/s %8.1f ms %8.2f GF/s %7.1fx %10.2e", n,
+         gms, gflop / (gms / 1e3), cms, gflop / (cms / 1e3), cms / gms, maxerr);
+    if (maxerr > 1e-2) {
+      logf("FAIL: n=%d exceeds tolerance", n);
+      return 0;
+    }
+  }
+  return 1;
+}
+
+}  // extern "C"
diff --git a/spike/webgpu/worker.js b/spike/webgpu/worker.js
new file mode 100644
index 0000000..32a612b
--- /dev/null
+++ b/spike/webgpu/worker.js
@@ -0,0 +1,44 @@
+// Acquire the WebGPU device here, in JS, and hand it to wasm fully formed.
+// That keeps every C++-side entry point synchronous -- the async surface is
+// this file, not the library.
+import createSpike from "./spike.js";
+
+async function main() {
+  if (!navigator.gpu) {
+    postMessage({ type: "fail", msg: "navigator.gpu missing (no WebGPU in this worker)" });
+    return;
+  }
+  const adapter = await navigator.gpu.requestAdapter();
+  if (!adapter) {
+    postMessage({ type: "fail", msg: "requestAdapter returned null" });
+    return;
+  }
+  const device = await adapter.requestDevice();
+  const info = adapter.info ?? {};
+
+  const mod = await createSpike({ preinitializedWebGPUDevice: device });
+
+  postMessage({
+    type: "ready",
+    adapter: `${info.vendor ?? "?"} / ${info.architecture ?? "?"} / ${info.description ?? ""}`,
+  });
+
+  // Under JSPI these exports return promises; under Asyncify they return
+  // values directly. `await` handles both, so the same worker drives either
+  // build without change.
+  onmessage = async () => {
+    const t0 = performance.now();
+    const okInit = await mod.ccall("wgpu_spike_init", "number", [], [], { async: true });
+    const okBench = okInit
+      ? await mod.ccall("wgpu_spike_bench", "number", [], [], { async: true })
+      : 0;
+    postMessage({
+      type: "done",
+      ok: !!okBench,
+      log: mod.UTF8ToString(mod._wgpu_spike_log()),
+      ms: performance.now() - t0,
+    });
+  };
+}
+
+main().catch((e) => postMessage({ type: "fail", msg: String(e && e.stack || e) }));

From 6598d6b5b4eefbed50a98c03d9f48926869b757c Mon Sep 17 00:00:00 2001
From: yhirose 
Date: Mon, 20 Jul 2026 00:09:13 -0400
Subject: [PATCH 04/12] refactor: move the tl::gpu facade out of cuda.h into
 gpu.h

The alias and gpu_available() sat at the bottom of cuda.h only because
that was the one place both tl::metal and tl::cuda were visible. With a
third backend coming (M10, WebGPU) that would mean editing the CUDA
header to add a browser GPU, so the choice moves to its own header.

No behaviour change: gpu.h includes both backend headers and holds the
same #if chain, and storage.h includes  instead of . The
cuda benches that include "cuda.h" use tl::cuda:: directly and are
unaffected.
---
 include/cuda.h    | 14 --------------
 include/gpu.h     | 30 ++++++++++++++++++++++++++++++
 include/metal.h   |  5 ++---
 include/storage.h |  2 +-
 4 files changed, 33 insertions(+), 18 deletions(-)
 create mode 100644 include/gpu.h

diff --git a/include/cuda.h b/include/cuda.h
index bf763b2..e2935bb 100644
--- a/include/cuda.h
+++ b/include/cuda.h
@@ -1400,18 +1400,4 @@ inline void cpu_barrier() {
 
 }  // namespace cuda
 
-// The GPU-backend facade: array.h and storage.h dispatch through tl::gpu, which
-// is Metal on Apple and CUDA elsewhere (the CUDA path only does real work when
-// TENSORLIB_CUDA is set; otherwise cuda:: is stubs, matching metal:: off-Apple).
-// One alias here is the single place the platform choice lives, so the eval
-// seam carries no #ifdefs. Both namespaces expose the identical API and share
-// tl::metal::kop, so the alias is a drop-in.
-#if defined(TENSORLIB_CUDA) && !defined(__APPLE__)
-namespace gpu = cuda;
-#else
-namespace gpu = metal;
-#endif
-
-inline bool gpu_available() { return gpu::available(); }
-
 }  // namespace tl
diff --git a/include/gpu.h b/include/gpu.h
new file mode 100644
index 0000000..84ef04d
--- /dev/null
+++ b/include/gpu.h
@@ -0,0 +1,30 @@
+#pragma once
+
+// The GPU-backend facade: array.h and storage.h dispatch through tl::gpu, so
+// the eval seam carries no platform #ifdefs. Every backend header exposes the
+// identical API (available/pending/flush/alloc/release/binary/unary/gemm/
+// row_op/sync_to_host) and shares tl::metal::kop, which makes the alias below a
+// drop-in.
+//
+// Each backend compiles to stubs unless its own gate holds, so including all of
+// them is free: metal.h is real only on __APPLE__, cuda.h only on
+// TENSORLIB_CUDA && !__APPLE__. The alias picks the one that can do real work.
+//
+// This lived at the bottom of cuda.h until M10 — the only place both namespaces
+// happened to be visible. That stopped scaling at the third backend, since
+// adding a browser GPU meant editing the CUDA header.
+
+#include "cuda.h"
+#include "metal.h"
+
+namespace tl {
+
+#if defined(TENSORLIB_CUDA) && !defined(__APPLE__)
+namespace gpu = cuda;
+#else
+namespace gpu = metal;
+#endif
+
+inline bool gpu_available() { return gpu::available(); }
+
+}  // namespace tl
diff --git a/include/metal.h b/include/metal.h
index bc8a5f8..fec5be5 100644
--- a/include/metal.h
+++ b/include/metal.h
@@ -485,8 +485,7 @@ inline void sync_to_host(void*, bool) {}
 
 }  // namespace metal
 
-// tl::gpu_available() and the tl::gpu facade (metal on Apple, cuda elsewhere)
-// are defined in cuda.h, which includes this header — so the choice of GPU
-// backend lives in one place and array.h's eval seam stays #ifdef-free.
+// tl::gpu_available() and the tl::gpu facade (which backend this build uses)
+// live in gpu.h — one place, so array.h's eval seam stays #ifdef-free.
 
 }  // namespace tl
diff --git a/include/storage.h b/include/storage.h
index c359871..e28416d 100644
--- a/include/storage.h
+++ b/include/storage.h
@@ -1,6 +1,6 @@
 #pragma once
 
-#include 
+#include 
 #include 
 
 #include 

From 83776a6189395f9b467cb738f2125217edfb1e3b Mon Sep 17 00:00:00 2001
From: yhirose 
Date: Mon, 20 Jul 2026 00:19:35 -0400
Subject: [PATCH 05/12] webgpu: host layer + WGSL sgemm, validated in the
 browser (M10 phases 1-2)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The first vertical slice: gemm runs on the GPU under wasm and agrees with
the ref oracle.

include/webgpu.h follows cuda.h, not metal.h. A WebGPU storage buffer has
no CPU-dereferenceable pointer, but gpu::alloc must return one, so each
allocation is a device buffer paired with a malloc'd host buffer and a
per-allocation dirty state drives lazy H2D/D2H — exactly CUDA's mirror.
That is what lets array.h and storage.h stay untouched.

flush() and sync_to_host() keep their synchronous signatures, as the spike
predicted: the instance is created with TimedWaitAny and the only suspend
points are OnSubmittedWorkDone and MapAsync. Dispatches accumulate in one
encoder and flush() submits once, which matters more here than on Metal
given the ~1ms dispatch floor.

View offsets resolved against the plan's open question: binding offsets
need 256-byte alignment that view offsets do not have, so bindings cover
whole buffers and offsets ride in the params block as element offsets,
folded into the indexing in WGSL. The uniform ring does use a dynamic
offset, hence the explicit bind group layout.

test/wasm/ relinks test_array.cpp unchanged — doctest builds under emcc,
and matches_gpu_oracle is already backend-agnostic. 24 cases / 335
assertions green in gpu mode on Chrome, with 149 GEMM dispatches confirmed
reaching the GPU (a temporarily instrumented run; the suite passes whether
or not the backend engages, so green alone proves nothing). Native ctest
unaffected.

Everything not yet ported returns false and falls to CPU, which is also
the Safari path when JSPI or a device is missing.
---
 include/gpu.h                     |  12 +-
 include/webgpu.h                  | 474 ++++++++++++++++++++++++++++++
 kernels/gen_wgsl_inc.sh           |  19 ++
 kernels/tensorlib_webgpu.wgsl     | 130 ++++++++
 kernels/tensorlib_webgpu_wgsl.inc | 133 +++++++++
 test/wasm/.gitignore              |   1 +
 test/wasm/build.sh                |  39 +++
 test/wasm/index.html              |  45 +++
 test/wasm/main_wasm.cpp           |  53 ++++
 test/wasm/worker.js               |  50 ++++
 10 files changed, 954 insertions(+), 2 deletions(-)
 create mode 100644 include/webgpu.h
 create mode 100755 kernels/gen_wgsl_inc.sh
 create mode 100644 kernels/tensorlib_webgpu.wgsl
 create mode 100644 kernels/tensorlib_webgpu_wgsl.inc
 create mode 100644 test/wasm/.gitignore
 create mode 100755 test/wasm/build.sh
 create mode 100644 test/wasm/index.html
 create mode 100644 test/wasm/main_wasm.cpp
 create mode 100644 test/wasm/worker.js

diff --git a/include/gpu.h b/include/gpu.h
index 84ef04d..874e8a5 100644
--- a/include/gpu.h
+++ b/include/gpu.h
@@ -8,7 +8,8 @@
 //
 // Each backend compiles to stubs unless its own gate holds, so including all of
 // them is free: metal.h is real only on __APPLE__, cuda.h only on
-// TENSORLIB_CUDA && !__APPLE__. The alias picks the one that can do real work.
+// TENSORLIB_CUDA && !__APPLE__, webgpu.h only on TENSORLIB_WEBGPU &&
+// __EMSCRIPTEN__. The alias picks the one that can do real work.
 //
 // This lived at the bottom of cuda.h until M10 — the only place both namespaces
 // happened to be visible. That stopped scaling at the third backend, since
@@ -16,10 +17,17 @@
 
 #include "cuda.h"
 #include "metal.h"
+#include "webgpu.h"
 
 namespace tl {
 
-#if defined(TENSORLIB_CUDA) && !defined(__APPLE__)
+// WebGPU is checked first: a wasm build defines neither __APPLE__ nor
+// TENSORLIB_CUDA, but a host build could define TENSORLIB_WEBGPU by accident
+// and should not silently take a backend that cannot work there — webgpu::
+// is stubs unless __EMSCRIPTEN__ too, so the order is safe either way.
+#if defined(TENSORLIB_WEBGPU) && defined(__EMSCRIPTEN__)
+namespace gpu = webgpu;
+#elif defined(TENSORLIB_CUDA) && !defined(__APPLE__)
 namespace gpu = cuda;
 #else
 namespace gpu = metal;
diff --git a/include/webgpu.h b/include/webgpu.h
new file mode 100644
index 0000000..937a737
--- /dev/null
+++ b/include/webgpu.h
@@ -0,0 +1,474 @@
+#pragma once
+
+// Own WebGPU backend (M10) — the browser GPU backend, mirroring metal.h and
+// cuda.h. Kernels are hand-written WGSL (kernels/tensorlib_webgpu.wgsl),
+// compiled at first use from a committed C-string .inc. No vendor library.
+//
+// Memory: a persistent host/device MIRROR per allocation, copied wholesale
+// from cuda.h. A WebGPU storage buffer has no CPU-dereferenceable pointer, but
+// gpu::alloc(bytes, float** contents) must hand one back (storage.h keeps it as
+// storage::ptr, and every ref::/accel::/cpu:: path then treats it as ordinary
+// memory). So each allocation is a device buffer (`native`, what kernels bind)
+// paired with a malloc'd host buffer (`contents`), and a per-allocation dirty
+// state drives lazy copies: H2D before a kernel reads a host-dirty buffer, D2H
+// before the CPU reads a device-dirty one (array::raw()/data() →
+// gpu::sync_to_host). This is why array.h and storage.h need no changes.
+//
+// Async: flush() and sync_to_host() keep their SYNCHRONOUS signatures. The
+// instance is created with the TimedWaitAny feature, which makes
+// wgpuInstanceWaitAny(timeout > 0) a legal blocking wait; emdawnwebgpu
+// implements it by suspending through Asyncify/JSPI. The suspend surface is
+// exactly two call sites — OnSubmittedWorkDone and MapAsync. Uploads
+// (WriteBuffer) are queued, never awaited. See spike/webgpu/README.md.
+//
+// Batching: one command encoder accumulates dispatches and flush() submits and
+// waits once, as metal.h does. This matters more here than on Metal — a
+// dispatch+wait has a ~0.6-1.0 ms fixed floor in the browser, so a flush per op
+// would be dominated by it.
+//
+// Real implementation is gated on TENSORLIB_WEBGPU && __EMSCRIPTEN__; a plain
+// build gets the stubs below. If the link lacks JSPI, or JS handed in no
+// device, CreateInstance/device acquisition fails and available() stays false —
+// every op then routes to CPU, which is also the Safari fallback.
+
+#include 
+
+#include "metal.h"  // reuse tl::metal::kop (platform-independent op enum)
+#include "types.h"
+
+namespace tl {
+namespace webgpu {
+
+using kop = tl::metal::kop;
+
+#if defined(TENSORLIB_WEBGPU) && defined(__EMSCRIPTEN__)
+
+}  // namespace webgpu
+}  // namespace tl
+
+// emdawnwebgpu declares emscripten_webgpu_get_device() in webgpu.h itself,
+// not in emscripten/html5_webgpu.h (that is the old built-in binding's home).
+#include 
+
+#include 
+#include 
+#include 
+#include 
+
+namespace tl {
+namespace webgpu {
+
+inline const char* wgsl_source_() {
+  static const char* src =
+#include "tensorlib_webgpu_wgsl.inc"
+      ;
+  return src;
+}
+
+// Uniform params, laid out to match the WGSL Params struct (16 x 4 bytes).
+struct gemm_params {
+  uint32_t M, N, K;
+  uint32_t lda, ldb, ldc;
+  uint32_t a_off, b_off, c_off;
+  uint32_t ta, tb;
+  float scale, offset;
+  uint32_t pad0, pad1, pad2;
+};
+
+// Dynamic uniform offsets must be a multiple of the adapter's
+// minUniformBufferOffsetAlignment; 256 is the spec's guaranteed-safe maximum.
+constexpr uint64_t kUniformSlot = 256;
+constexpr uint32_t kUniformSlots = 256;  // one flush can batch this many
+
+struct context {
+  wgpu::Instance instance;
+  wgpu::Device device;
+  wgpu::Queue queue;
+  wgpu::BindGroupLayout bgl;
+  wgpu::ComputePipeline sgemm;
+  wgpu::Buffer uniforms;  // ring of kUniformSlots x kUniformSlot bytes
+  bool ready = false;
+  bool pending = false;
+
+  // The open encoder. Created lazily on the first dispatch after a flush, so
+  // an idle backend submits nothing.
+  wgpu::CommandEncoder enc;
+  uint32_t slot = 0;  // next free uniform ring slot
+
+  // Host/device mirror per allocation, keyed by the opaque handle alloc()
+  // returns as `native`. Views sharing a storage share the key, so one dirty
+  // state serves every view. `where` tracks which copy is live.
+  enum loc { HOST, DEVICE, BOTH };
+  struct mirror {
+    float* host = nullptr;   // CPU-side buffer (storage.ptr)
+    wgpu::Buffer dev;        // device buffer
+    size_t bytes = 0;
+    loc where = HOST;
+  };
+  std::unordered_map mirrors;
+
+  // Size-keyed free lists (like Metal's MTLBuffer pool and CUDA's). Repeated
+  // alloc/free of identical shapes is the common case, and per-dispatch
+  // allocation would compound the fixed dispatch floor.
+  std::unordered_map>> pool;
+  std::unordered_map> staging_pool;
+
+  static context& get() {
+    static auto* c = new context();  // leaked: outlives all storage deleters
+    return *c;
+  }
+
+  context() {
+    // TimedWaitAny is what makes WaitAny(timeout > 0) legal. emdawnwebgpu
+    // refuses to create the instance without JSPI/Asyncify, so a link missing
+    // it fails loudly here rather than deadlocking later.
+    wgpu::InstanceFeatureName features[] = {
+        wgpu::InstanceFeatureName::TimedWaitAny};
+    wgpu::InstanceDescriptor idesc = {};
+    idesc.requiredFeatureCount = 1;
+    idesc.requiredFeatures = features;
+    instance = wgpu::CreateInstance(&idesc);
+    if (!instance) return;
+
+    // No adapter/device round-trip on this side: JS already did it and passed
+    // the result as Module.preinitializedWebGPUDevice.
+    device = wgpu::Device::Acquire(emscripten_webgpu_get_device());
+    if (!device) return;
+    queue = device.GetQueue();
+
+    wgpu::ShaderSourceWGSL wgsl = {};
+    wgsl.code = wgsl_source_();
+    wgpu::ShaderModuleDescriptor smd = {};
+    smd.nextInChain = &wgsl;
+    wgpu::ShaderModule mod = device.CreateShaderModule(&smd);
+    if (!mod) return;
+
+    // Explicit layout rather than GetBindGroupLayout(0): the auto-generated
+    // one has no dynamic offset on the uniform binding, which the ring needs.
+    wgpu::BindGroupLayoutEntry be[4] = {};
+    for (int i = 0; i < 3; ++i) {
+      be[i].binding = i;
+      be[i].visibility = wgpu::ShaderStage::Compute;
+      be[i].buffer.type = i == 2 ? wgpu::BufferBindingType::Storage
+                                 : wgpu::BufferBindingType::ReadOnlyStorage;
+    }
+    be[3].binding = 3;
+    be[3].visibility = wgpu::ShaderStage::Compute;
+    be[3].buffer.type = wgpu::BufferBindingType::Uniform;
+    be[3].buffer.hasDynamicOffset = true;
+    be[3].buffer.minBindingSize = sizeof(gemm_params);
+    wgpu::BindGroupLayoutDescriptor bgld = {};
+    bgld.entryCount = 4;
+    bgld.entries = be;
+    bgl = device.CreateBindGroupLayout(&bgld);
+
+    wgpu::PipelineLayoutDescriptor pld = {};
+    pld.bindGroupLayoutCount = 1;
+    pld.bindGroupLayouts = &bgl;
+    wgpu::PipelineLayout pl = device.CreatePipelineLayout(&pld);
+
+    wgpu::ComputePipelineDescriptor pd = {};
+    pd.layout = pl;
+    pd.compute.module = mod;
+    pd.compute.entryPoint = "sgemm";
+    sgemm = device.CreateComputePipeline(&pd);
+    if (!sgemm) return;
+
+    wgpu::BufferDescriptor ud = {};
+    ud.size = kUniformSlot * kUniformSlots;
+    ud.usage = wgpu::BufferUsage::Uniform | wgpu::BufferUsage::CopyDst;
+    uniforms = device.CreateBuffer(&ud);
+    if (!uniforms) return;
+
+    ready = true;
+  }
+
+  // Blocking wait on a single future — the one place anything suspends.
+  bool wait(wgpu::Future f) {
+    return instance.WaitAny(f, UINT64_MAX) == wgpu::WaitStatus::Success;
+  }
+
+  mirror* mirror_(void* native) {
+    auto it = mirrors.find(native);
+    return it == mirrors.end() ? nullptr : &it->second;
+  }
+
+  // A kernel is about to READ this buffer: ensure the device copy is current.
+  //
+  // WriteBuffer executes in queue order, i.e. ahead of anything still sitting
+  // in the unsubmitted encoder. That is safe precisely because a buffer in
+  // HOST state has no encoded command touching it: a pending kernel write
+  // would have set DEVICE, and a pending kernel read would have come through
+  // here and set BOTH.
+  void device_read_(void* native) {
+    mirror* m = mirror_(native);
+    if (m && m->where == HOST) {
+      queue.WriteBuffer(m->dev, 0, m->host, m->bytes);
+      m->where = BOTH;
+    }
+  }
+
+  // A kernel is about to WRITE this buffer: it becomes the live copy.
+  void device_write_(void* native) {
+    if (mirror* m = mirror_(native)) m->where = DEVICE;
+  }
+
+  wgpu::Buffer staging_(size_t bytes) {
+    auto it = staging_pool.find(bytes);
+    if (it != staging_pool.end() && !it->second.empty()) {
+      wgpu::Buffer b = it->second.back();
+      it->second.pop_back();
+      return b;
+    }
+    wgpu::BufferDescriptor d = {};
+    d.size = bytes;
+    d.usage = wgpu::BufferUsage::MapRead | wgpu::BufferUsage::CopyDst;
+    return device.CreateBuffer(&d);
+  }
+};
+
+inline bool available() { return context::get().ready; }
+inline bool pending() { return context::get().pending; }
+
+// End the batch: submit the accumulated encoder and block until the GPU
+// finishes (MLX-style eval).
+inline void flush() {
+  auto& c = context::get();
+  if (!c.pending) return;
+  c.pending = false;
+  c.slot = 0;
+  if (!c.enc) return;
+  wgpu::CommandBuffer cmds = c.enc.Finish();
+  c.enc = nullptr;
+  c.queue.Submit(1, &cmds);
+  c.wait(c.queue.OnSubmittedWorkDone(
+      wgpu::CallbackMode::WaitAnyOnly,
+      [](wgpu::QueueWorkDoneStatus, wgpu::StringView) {}));
+}
+
+// Mirror allocation: a device buffer paired with a host buffer (returned via
+// `contents`). They are DISTINCT memory — the dirty state copies between them
+// on demand. The returned handle is an opaque token, not a pointer to
+// anything dereferenceable; it is only ever a key back into `mirrors`.
+inline void* alloc(int64_t bytes, float** contents) {
+  auto& c = context::get();
+  if (!c.ready) return nullptr;
+  size_t nb = bytes > 0 ? (size_t)bytes : 4;
+  nb = (nb + 3) & ~size_t(3);  // WebGPU buffer sizes must be 4-byte multiples
+
+  wgpu::Buffer dev;
+  float* host = nullptr;
+  auto it = c.pool.find(nb);  // reuse a recycled buffer of this exact size
+  if (it != c.pool.end() && !it->second.empty()) {
+    dev = it->second.back().first;
+    host = it->second.back().second;
+    it->second.pop_back();
+  } else {
+    wgpu::BufferDescriptor d = {};
+    d.size = nb;
+    d.usage = wgpu::BufferUsage::Storage | wgpu::BufferUsage::CopyDst |
+              wgpu::BufferUsage::CopySrc;
+    dev = c.device.CreateBuffer(&d);
+    if (!dev) return nullptr;
+    host = static_cast(std::malloc(nb));
+    if (!host) return nullptr;
+  }
+
+  // The token has to be unique and stable for the allocation's lifetime; the
+  // host pointer is both, and malloc will not hand out the same address twice
+  // while it is live.
+  void* token = host;
+  c.mirrors[token] = context::mirror{host, dev, nb, context::HOST};
+  if (contents) *contents = host;
+  return token;
+}
+
+inline void release(void* buf, int64_t, float*) {
+  auto& c = context::get();
+  if (!c.ready || !buf) return;
+  auto it = c.mirrors.find(buf);
+  if (it == c.mirrors.end()) return;
+  c.pool[it->second.bytes].push_back({it->second.dev, it->second.host});
+  c.mirrors.erase(it);
+}
+
+// Reconcile a buffer for a CPU access: flush pending kernels, then D2H if the
+// device holds the live copy. for_write invalidates the device copy (the host
+// is about to mutate it). No-op for heap storages / unknown handles.
+inline void sync_to_host(void* native, bool for_write) {
+  auto& c = context::get();
+  if (!c.ready || !native) return;
+  context::mirror* m = c.mirror_(native);
+  if (!m) return;
+  if (c.pending) flush();
+  if (m->where == context::DEVICE) {
+    // No CPU-visible pointer to read from: copy device -> a MapRead staging
+    // buffer, map it (the second and last suspend point), memcpy out.
+    wgpu::Buffer stg = c.staging_(m->bytes);
+    wgpu::CommandEncoder e = c.device.CreateCommandEncoder();
+    e.CopyBufferToBuffer(m->dev, 0, stg, 0, m->bytes);
+    wgpu::CommandBuffer cmds = e.Finish();
+    c.queue.Submit(1, &cmds);
+
+    bool ok = false;
+    if (c.wait(stg.MapAsync(wgpu::MapMode::Read, 0, m->bytes,
+                            wgpu::CallbackMode::WaitAnyOnly,
+                            [&ok](wgpu::MapAsyncStatus s, wgpu::StringView) {
+                              ok = (s == wgpu::MapAsyncStatus::Success);
+                            })) &&
+        ok) {
+      if (const void* src = stg.GetConstMappedRange(0, m->bytes))
+        std::memcpy(m->host, src, m->bytes);
+      stg.Unmap();
+      m->where = context::BOTH;
+    }
+    c.staging_pool[m->bytes].push_back(stg);
+  }
+  if (for_write) m->where = context::HOST;
+}
+
+// out = op(a) @ op(b) * scale + offset. Offsets are BYTE offsets into the
+// storage; the kernel wants elements, so they must be 4-aligned (they always
+// are for f32 views).
+inline bool gemm(void* a, int64_t ao, int64_t lda, bool ta, void* b, int64_t bo,
+                 int64_t ldb, bool tb, void* out, int64_t oo, int64_t m,
+                 int64_t n, int64_t k, float scale, float offset) {
+  auto& c = context::get();
+  if (!c.ready) return false;
+  if (m <= 0 || n <= 0 || k <= 0) return false;
+  if ((ao % 4) || (bo % 4) || (oo % 4)) return false;
+  // Every operand must be a tracked mirror: an untracked pointer is a heap
+  // storage with no device buffer to bind, so this op belongs on the CPU.
+  context::mirror* ma = c.mirror_(a);
+  context::mirror* mb = c.mirror_(b);
+  context::mirror* mo = c.mirror_(out);
+  if (!ma || !mb || !mo) return false;
+  if (c.slot >= kUniformSlots) flush();  // ring exhausted; start a new batch
+
+  c.device_read_(a);
+  c.device_read_(b);
+  c.device_write_(out);
+
+  gemm_params p = {};
+  p.M = (uint32_t)m;
+  p.N = (uint32_t)n;
+  p.K = (uint32_t)k;
+  p.lda = (uint32_t)lda;
+  p.ldb = (uint32_t)ldb;
+  p.ldc = (uint32_t)n;  // the eval seam always hands us a contiguous output
+  p.a_off = (uint32_t)(ao / 4);
+  p.b_off = (uint32_t)(bo / 4);
+  p.c_off = (uint32_t)(oo / 4);
+  p.ta = ta ? 1u : 0u;
+  p.tb = tb ? 1u : 0u;
+  p.scale = scale;
+  p.offset = offset;
+
+  const uint32_t off = c.slot++ * (uint32_t)kUniformSlot;
+  c.queue.WriteBuffer(c.uniforms, off, &p, sizeof(p));
+
+  wgpu::BindGroupEntry e[4] = {};
+  e[0].binding = 0; e[0].buffer = ma->dev; e[0].size = ma->bytes;
+  e[1].binding = 1; e[1].buffer = mb->dev; e[1].size = mb->bytes;
+  e[2].binding = 2; e[2].buffer = mo->dev; e[2].size = mo->bytes;
+  e[3].binding = 3; e[3].buffer = c.uniforms; e[3].size = sizeof(gemm_params);
+  wgpu::BindGroupDescriptor bgd = {};
+  bgd.layout = c.bgl;
+  bgd.entryCount = 4;
+  bgd.entries = e;
+  wgpu::BindGroup bg = c.device.CreateBindGroup(&bgd);
+
+  if (!c.enc) c.enc = c.device.CreateCommandEncoder();
+  wgpu::ComputePassEncoder pass = c.enc.BeginComputePass();
+  pass.SetPipeline(c.sgemm);
+  pass.SetBindGroup(0, bg, 1, &off);
+  pass.DispatchWorkgroups((uint32_t)((n + 63) / 64), (uint32_t)((m + 63) / 64),
+                          1);
+  pass.End();
+
+  c.pending = true;
+  return true;
+}
+
+// ---- Not ported yet (Phase 3+). Returning false routes the op to CPU, which
+// is why each phase lands in a working state.
+inline bool binary(kop, void*, int64_t, void*, int64_t, void*, int64_t, int64_t,
+                   float, float) {
+  return false;
+}
+inline bool binary_bcast(kop, void*, int64_t, int64_t, int64_t, void*, int64_t,
+                         int64_t, int64_t, void*, int64_t, int64_t, int64_t,
+                         float, float) {
+  return false;
+}
+inline bool unary(kop, void*, int64_t, void*, int64_t, int64_t, float, float) {
+  return false;
+}
+inline bool row_op(kop, void*, int64_t, void*, int64_t, int64_t, int64_t, float,
+                   float) {
+  return false;
+}
+inline bool gemv_f32(void*, void*, void*, int64_t, int64_t) { return false; }
+inline bool gemv_bf16(void*, void*, void*, int64_t, int64_t) { return false; }
+inline bool attn_decode(void*, void*, void*, void*, int64_t, int64_t, int64_t,
+                        int64_t, int64_t, float) {
+  return false;
+}
+inline bool rope(void*, void*, int64_t, int64_t, int64_t, int64_t, float) {
+  return false;
+}
+inline bool gemv_q4(void*, void*, void*, void*, int64_t, int64_t, int64_t) {
+  return false;
+}
+
+#else  // !(TENSORLIB_WEBGPU && __EMSCRIPTEN__) — stubs, as in metal.h
+
+inline bool available() { return false; }
+inline bool pending() { return false; }
+inline void flush() {}
+inline void* alloc(int64_t, float**) { return nullptr; }
+inline void release(void*, int64_t, float*) {}
+inline void sync_to_host(void*, bool) {}
+inline bool binary(kop, void*, int64_t, void*, int64_t, void*, int64_t, int64_t,
+                   float, float) {
+  return false;
+}
+inline bool binary_bcast(kop, void*, int64_t, int64_t, int64_t, void*, int64_t,
+                         int64_t, int64_t, void*, int64_t, int64_t, int64_t,
+                         float, float) {
+  return false;
+}
+inline bool unary(kop, void*, int64_t, void*, int64_t, int64_t, float, float) {
+  return false;
+}
+inline bool gemm(void*, int64_t, int64_t, bool, void*, int64_t, int64_t, bool,
+                 void*, int64_t, int64_t, int64_t, int64_t, float, float) {
+  return false;
+}
+inline bool row_op(kop, void*, int64_t, void*, int64_t, int64_t, int64_t, float,
+                   float) {
+  return false;
+}
+inline bool gemv_f32(void*, void*, void*, int64_t, int64_t) { return false; }
+inline bool gemv_bf16(void*, void*, void*, int64_t, int64_t) { return false; }
+inline bool attn_decode(void*, void*, void*, void*, int64_t, int64_t, int64_t,
+                        int64_t, int64_t, float) {
+  return false;
+}
+inline bool rope(void*, void*, int64_t, int64_t, int64_t, int64_t, float) {
+  return false;
+}
+inline bool gemv_q4(void*, void*, void*, void*, int64_t, int64_t, int64_t) {
+  return false;
+}
+
+#endif
+
+// Every CPU-side buffer read funnels through array::raw()/data(), which call
+// this: one choke point makes mixed CPU/GPU graphs safe.
+inline void cpu_barrier() {
+  if (pending()) flush();
+}
+
+}  // namespace webgpu
+}  // namespace tl
diff --git a/kernels/gen_wgsl_inc.sh b/kernels/gen_wgsl_inc.sh
new file mode 100755
index 0000000..35affd2
--- /dev/null
+++ b/kernels/gen_wgsl_inc.sh
@@ -0,0 +1,19 @@
+#!/usr/bin/env bash
+# WGSL -> C string literal, the WebGPU counterpart of the CUDA backend's
+# bin2c step. Run this after editing tensorlib_webgpu.wgsl and commit the
+# result: the wasm build is a flat emcc line that never runs CMake, so the
+# .inc cannot be a build-time artifact the way the PTX one is.
+set -euo pipefail
+
+HERE="$(cd "$(dirname "$0")" && pwd)"
+SRC="$HERE/tensorlib_webgpu.wgsl"
+OUT="$HERE/tensorlib_webgpu_wgsl.inc"
+
+{
+  echo "// Generated by gen_wgsl_inc.sh from tensorlib_webgpu.wgsl -- do not edit."
+  echo 'R"WGSL('
+  cat "$SRC"
+  echo ')WGSL"'
+} > "$OUT"
+
+echo "wrote $OUT ($(wc -l < "$OUT") lines)"
diff --git a/kernels/tensorlib_webgpu.wgsl b/kernels/tensorlib_webgpu.wgsl
new file mode 100644
index 0000000..c18398e
--- /dev/null
+++ b/kernels/tensorlib_webgpu.wgsl
@@ -0,0 +1,130 @@
+// WebGPU compute kernels for cpp-tensorlib (M10). Hand-written WGSL, no
+// vendor library — the same stance the Metal and CUDA backends take.
+//
+// Compiled at first use from the C string in tensorlib_webgpu_wgsl.inc, which
+// kernels/gen_wgsl_inc.sh generates from this file. The generated .inc is
+// committed because the wasm build is a flat emcc line that never runs CMake
+// (the CUDA backend's PTX goes through bin2c for the same reason).
+//
+// View offsets arrive as ELEMENT offsets in the params block and are folded
+// into the indexing here, rather than as bind-group binding offsets: WebGPU
+// requires those to be 256-byte aligned, which an arbitrary view offset is
+// not. So every binding covers its whole buffer. (CUDA instead folds offsets
+// host-side into the pointer it passes, which WebGPU has no equivalent of.)
+
+struct Params {
+  M : u32,
+  N : u32,
+  K : u32,
+  lda : u32,
+  ldb : u32,
+  ldc : u32,
+  a_off : u32,
+  b_off : u32,
+  c_off : u32,
+  ta : u32,
+  tb : u32,
+  scale : f32,
+  offset : f32,
+  _pad0 : u32,
+  _pad1 : u32,
+  _pad2 : u32,
+};
+
+@group(0) @binding(0) var       A : array;
+@group(0) @binding(1) var       B : array;
+@group(0) @binding(2) var C : array;
+@group(0) @binding(3) var             p : Params;
+
+// ---- sgemm: C(M,N) = op(A)(M,K) @ op(B)(K,N) * scale + offset
+//
+// 64x64 workgroup tile, 16x16 = 256 invocations, each holding a 4x4 register
+// accumulator. Mirrors the shape of the Metal sgemm_64_ kernel; MMA intrinsics
+// have no WGSL equivalent, so the inner product is plain FMA over registers.
+// Measured at ~580-630 GF/s for n=1024 on an M1 Pro (see spike/webgpu).
+
+const BM : u32 = 64u;
+const BN : u32 = 64u;
+const BK : u32 = 16u;
+const TM : u32 = 4u;
+const TN : u32 = 4u;
+const THREADS : u32 = 256u;
+
+var As : array;  // BM * BK
+var Bs : array;  // BK * BN
+
+// Row/column strides for a possibly-transposed operand: transposing swaps
+// which axis walks by the leading dimension (cf. metal_kernels.metal:119-120).
+fn a_index(m : u32, k : u32) -> u32 {
+  let rs = select(p.lda, 1u, p.ta == 1u);
+  let cs = select(1u, p.lda, p.ta == 1u);
+  return p.a_off + m * rs + k * cs;
+}
+
+fn b_index(k : u32, n : u32) -> u32 {
+  let rs = select(p.ldb, 1u, p.tb == 1u);
+  let cs = select(1u, p.ldb, p.tb == 1u);
+  return p.b_off + k * rs + n * cs;
+}
+
+@compute @workgroup_size(16, 16, 1)
+fn sgemm(@builtin(workgroup_id) wg : vec3,
+         @builtin(local_invocation_id) lid : vec3) {
+  let m_base = wg.y * BM;
+  let n_base = wg.x * BN;
+  let tid = lid.y * 16u + lid.x;
+
+  var acc : array;  // TM * TN, zero-initialized
+
+  let n_tiles = (p.K + BK - 1u) / BK;
+  for (var kt : u32 = 0u; kt < n_tiles; kt = kt + 1u) {
+    let k0 = kt * BK;
+
+    // Stage A tile (BM x BK) and B tile (BK x BN), 4 elements per invocation
+    // each. Out-of-range reads are zero-filled so the tail tile needs no
+    // special case in the inner loop.
+    for (var s : u32 = 0u; s < 4u; s = s + 1u) {
+      let i = tid + s * THREADS;
+
+      let am = m_base + i / BK;
+      let ak = k0 + i % BK;
+      let a_ok = am < p.M && ak < p.K;
+      As[i] = select(0.0, A[a_index(am, ak)], a_ok);
+
+      let bk = k0 + i / BN;
+      let bn = n_base + i % BN;
+      let b_ok = bk < p.K && bn < p.N;
+      Bs[i] = select(0.0, B[b_index(bk, bn)], b_ok);
+    }
+
+    workgroupBarrier();
+
+    for (var kk : u32 = 0u; kk < BK; kk = kk + 1u) {
+      var av : array;
+      var bv : array;
+      for (var i : u32 = 0u; i < TM; i = i + 1u) {
+        av[i] = As[(lid.y * TM + i) * BK + kk];
+      }
+      for (var j : u32 = 0u; j < TN; j = j + 1u) {
+        bv[j] = Bs[kk * BN + lid.x * TN + j];
+      }
+      for (var i : u32 = 0u; i < TM; i = i + 1u) {
+        for (var j : u32 = 0u; j < TN; j = j + 1u) {
+          acc[i * TN + j] = fma(av[i], bv[j], acc[i * TN + j]);
+        }
+      }
+    }
+
+    workgroupBarrier();
+  }
+
+  for (var i : u32 = 0u; i < TM; i = i + 1u) {
+    let m = m_base + lid.y * TM + i;
+    if (m >= p.M) { continue; }
+    for (var j : u32 = 0u; j < TN; j = j + 1u) {
+      let n = n_base + lid.x * TN + j;
+      if (n >= p.N) { continue; }
+      C[p.c_off + m * p.ldc + n] = acc[i * TN + j] * p.scale + p.offset;
+    }
+  }
+}
diff --git a/kernels/tensorlib_webgpu_wgsl.inc b/kernels/tensorlib_webgpu_wgsl.inc
new file mode 100644
index 0000000..2cd0259
--- /dev/null
+++ b/kernels/tensorlib_webgpu_wgsl.inc
@@ -0,0 +1,133 @@
+// Generated by gen_wgsl_inc.sh from tensorlib_webgpu.wgsl -- do not edit.
+R"WGSL(
+// WebGPU compute kernels for cpp-tensorlib (M10). Hand-written WGSL, no
+// vendor library — the same stance the Metal and CUDA backends take.
+//
+// Compiled at first use from the C string in tensorlib_webgpu_wgsl.inc, which
+// kernels/gen_wgsl_inc.sh generates from this file. The generated .inc is
+// committed because the wasm build is a flat emcc line that never runs CMake
+// (the CUDA backend's PTX goes through bin2c for the same reason).
+//
+// View offsets arrive as ELEMENT offsets in the params block and are folded
+// into the indexing here, rather than as bind-group binding offsets: WebGPU
+// requires those to be 256-byte aligned, which an arbitrary view offset is
+// not. So every binding covers its whole buffer. (CUDA instead folds offsets
+// host-side into the pointer it passes, which WebGPU has no equivalent of.)
+
+struct Params {
+  M : u32,
+  N : u32,
+  K : u32,
+  lda : u32,
+  ldb : u32,
+  ldc : u32,
+  a_off : u32,
+  b_off : u32,
+  c_off : u32,
+  ta : u32,
+  tb : u32,
+  scale : f32,
+  offset : f32,
+  _pad0 : u32,
+  _pad1 : u32,
+  _pad2 : u32,
+};
+
+@group(0) @binding(0) var       A : array;
+@group(0) @binding(1) var       B : array;
+@group(0) @binding(2) var C : array;
+@group(0) @binding(3) var             p : Params;
+
+// ---- sgemm: C(M,N) = op(A)(M,K) @ op(B)(K,N) * scale + offset
+//
+// 64x64 workgroup tile, 16x16 = 256 invocations, each holding a 4x4 register
+// accumulator. Mirrors the shape of the Metal sgemm_64_ kernel; MMA intrinsics
+// have no WGSL equivalent, so the inner product is plain FMA over registers.
+// Measured at ~580-630 GF/s for n=1024 on an M1 Pro (see spike/webgpu).
+
+const BM : u32 = 64u;
+const BN : u32 = 64u;
+const BK : u32 = 16u;
+const TM : u32 = 4u;
+const TN : u32 = 4u;
+const THREADS : u32 = 256u;
+
+var As : array;  // BM * BK
+var Bs : array;  // BK * BN
+
+// Row/column strides for a possibly-transposed operand: transposing swaps
+// which axis walks by the leading dimension (cf. metal_kernels.metal:119-120).
+fn a_index(m : u32, k : u32) -> u32 {
+  let rs = select(p.lda, 1u, p.ta == 1u);
+  let cs = select(1u, p.lda, p.ta == 1u);
+  return p.a_off + m * rs + k * cs;
+}
+
+fn b_index(k : u32, n : u32) -> u32 {
+  let rs = select(p.ldb, 1u, p.tb == 1u);
+  let cs = select(1u, p.ldb, p.tb == 1u);
+  return p.b_off + k * rs + n * cs;
+}
+
+@compute @workgroup_size(16, 16, 1)
+fn sgemm(@builtin(workgroup_id) wg : vec3,
+         @builtin(local_invocation_id) lid : vec3) {
+  let m_base = wg.y * BM;
+  let n_base = wg.x * BN;
+  let tid = lid.y * 16u + lid.x;
+
+  var acc : array;  // TM * TN, zero-initialized
+
+  let n_tiles = (p.K + BK - 1u) / BK;
+  for (var kt : u32 = 0u; kt < n_tiles; kt = kt + 1u) {
+    let k0 = kt * BK;
+
+    // Stage A tile (BM x BK) and B tile (BK x BN), 4 elements per invocation
+    // each. Out-of-range reads are zero-filled so the tail tile needs no
+    // special case in the inner loop.
+    for (var s : u32 = 0u; s < 4u; s = s + 1u) {
+      let i = tid + s * THREADS;
+
+      let am = m_base + i / BK;
+      let ak = k0 + i % BK;
+      let a_ok = am < p.M && ak < p.K;
+      As[i] = select(0.0, A[a_index(am, ak)], a_ok);
+
+      let bk = k0 + i / BN;
+      let bn = n_base + i % BN;
+      let b_ok = bk < p.K && bn < p.N;
+      Bs[i] = select(0.0, B[b_index(bk, bn)], b_ok);
+    }
+
+    workgroupBarrier();
+
+    for (var kk : u32 = 0u; kk < BK; kk = kk + 1u) {
+      var av : array;
+      var bv : array;
+      for (var i : u32 = 0u; i < TM; i = i + 1u) {
+        av[i] = As[(lid.y * TM + i) * BK + kk];
+      }
+      for (var j : u32 = 0u; j < TN; j = j + 1u) {
+        bv[j] = Bs[kk * BN + lid.x * TN + j];
+      }
+      for (var i : u32 = 0u; i < TM; i = i + 1u) {
+        for (var j : u32 = 0u; j < TN; j = j + 1u) {
+          acc[i * TN + j] = fma(av[i], bv[j], acc[i * TN + j]);
+        }
+      }
+    }
+
+    workgroupBarrier();
+  }
+
+  for (var i : u32 = 0u; i < TM; i = i + 1u) {
+    let m = m_base + lid.y * TM + i;
+    if (m >= p.M) { continue; }
+    for (var j : u32 = 0u; j < TN; j = j + 1u) {
+      let n = n_base + lid.x * TN + j;
+      if (n >= p.N) { continue; }
+      C[p.c_off + m * p.ldc + n] = acc[i * TN + j] * p.scale + p.offset;
+    }
+  }
+}
+)WGSL"
diff --git a/test/wasm/.gitignore b/test/wasm/.gitignore
new file mode 100644
index 0000000..45ddf0a
--- /dev/null
+++ b/test/wasm/.gitignore
@@ -0,0 +1 @@
+site/
diff --git a/test/wasm/build.sh b/test/wasm/build.sh
new file mode 100755
index 0000000..41f5e0e
--- /dev/null
+++ b/test/wasm/build.sh
@@ -0,0 +1,39 @@
+#!/usr/bin/env bash
+# Build the browser test harness for the WebGPU backend.
+#
+#   ./build.sh && python3 -m http.server 8731
+#   # open http://localhost:8731/test/wasm/site/index.html
+#
+# A flat emcc line, not CMake — which is why the WGSL .inc is committed rather
+# than generated at build time (see kernels/gen_wgsl_inc.sh).
+set -euo pipefail
+
+HERE="$(cd "$(dirname "$0")" && pwd)"
+ROOT="$(cd "$HERE/../.." && pwd)"
+OUT="$HERE/site"
+
+source "${EMSDK_DIR:-$HOME/Projects/emsdk}/emsdk_env.sh" >/dev/null 2>&1
+
+mkdir -p "$OUT"
+
+# JSPI over Asyncify: smaller, and no conflict with -fwasm-exceptions (which
+# culebra's interpreter needs). The cost is that this wasm *requires* JSPI, so
+# it runs on Chrome only — acceptable for a test harness, and the shipping
+# playground pairs it with a CPU build behind a feature check.
+set -x
+emcc -std=c++23 -O2 -fwasm-exceptions \
+  --use-port=emdawnwebgpu \
+  -DTENSORLIB_WEBGPU \
+  -sJSPI=1 -sJSPI_EXPORTS=run_tests \
+  -sSTACK_SIZE=16MB -sALLOW_MEMORY_GROWTH=1 -sINITIAL_MEMORY=64MB \
+  -sMODULARIZE=1 -sEXPORT_ES6=1 -sEXPORT_NAME=createTests \
+  -sENVIRONMENT=web,worker \
+  -sEXPORTED_FUNCTIONS=_run_tests,_malloc,_free \
+  -sEXPORTED_RUNTIME_METHODS=ccall,cwrap,UTF8ToString \
+  -I"$ROOT/include" -I"$ROOT/kernels" -I"$ROOT/test" \
+  "$HERE/main_wasm.cpp" "$ROOT/test/test_array.cpp" \
+  -o "$OUT/tests.js"
+set +x
+
+cp "$HERE/index.html" "$HERE/worker.js" "$OUT/"
+echo "wasm size: $(stat -f%z "$OUT/tests.wasm") bytes"
diff --git a/test/wasm/index.html b/test/wasm/index.html
new file mode 100644
index 0000000..e0155a9
--- /dev/null
+++ b/test/wasm/index.html
@@ -0,0 +1,45 @@
+
+
+cpp-tensorlib WebGPU tests
+
+
cpp-tensorlib WebGPU tests

+
booting…

+

+

+
diff --git a/test/wasm/main_wasm.cpp b/test/wasm/main_wasm.cpp
new file mode 100644
index 0000000..a6b415f
--- /dev/null
+++ b/test/wasm/main_wasm.cpp
@@ -0,0 +1,53 @@
+// Browser test entry point for the WebGPU backend (M10, Phase 2).
+//
+// Native ctest cannot exercise this backend at all, so the same doctest suite
+// is relinked for wasm and driven from a page. test_array.cpp compiles in
+// unchanged: its matches_gpu_oracle helper (test/test_array.cpp:53-67) already
+// runs a graph in gpu mode and again through the ref oracle, is
+// backend-agnostic, and self-skips on !gpu_available() — so it exercises
+// WebGPU exactly as it exercises Metal, with no WebGPU-specific assertions to
+// keep in sync.
+//
+// Ops the backend has not ported yet return false and fall to CPU, so the
+// whole suite is expected green from Phase 1 onward; what changes per phase is
+// how much of it actually ran on the GPU.
+//
+// The verdict goes to the console as a single OK/FAIL line so a
+// browser-automation check can read it without scraping the page.
+
+#define DOCTEST_CONFIG_IMPLEMENT
+#include "doctest.h"
+
+#include 
+#include 
+
+#include 
+
+extern "C" {
+
+// Called from JS with {async: true}: under JSPI this returns a promise,
+// because any GPU wait beneath it suspends. Nothing in C++ is async.
+EMSCRIPTEN_KEEPALIVE int run_tests(int gpu_mode) {
+  if (gpu_mode) {
+    tl::use_gpu();
+  } else {
+    tl::use_cpu();
+  }
+  std::printf("[tensorlib_wasm] mode=%s gpu_available=%d\n",
+              gpu_mode ? "gpu" : "cpu", tl::gpu_available() ? 1 : 0);
+
+  // gpu_available() false in gpu mode means the device or JSPI is missing and
+  // every GPU assertion would trivially pass — report that rather than a green
+  // run that proved nothing.
+  if (gpu_mode && !tl::gpu_available()) {
+    std::printf("[tensorlib_wasm] FAIL: gpu mode requested but no device\n");
+    return 0;
+  }
+
+  doctest::Context ctx;
+  int failed = ctx.run();
+  std::printf("[tensorlib_wasm] %s\n", failed ? "FAIL" : "OK");
+  return failed ? 0 : 1;
+}
+
+}  // extern "C"
diff --git a/test/wasm/worker.js b/test/wasm/worker.js
new file mode 100644
index 0000000..c887db9
--- /dev/null
+++ b/test/wasm/worker.js
@@ -0,0 +1,50 @@
+// Acquire the WebGPU device here, in JS, and hand it to wasm fully formed.
+// That keeps every C++-side entry point synchronous — the async surface is
+// this file, not the library.
+import createTests from "./tests.js";
+
+let log = "";
+
+async function main() {
+  if (!navigator.gpu) {
+    postMessage({ type: "fail", msg: "navigator.gpu missing (no WebGPU in this worker)" });
+    return;
+  }
+  const adapter = await navigator.gpu.requestAdapter();
+  if (!adapter) {
+    postMessage({ type: "fail", msg: "requestAdapter returned null" });
+    return;
+  }
+  const device = await adapter.requestDevice();
+  const info = adapter.info ?? {};
+
+  // doctest writes its report to stdout; capture it rather than let it vanish
+  // into the worker console.
+  const mod = await createTests({
+    preinitializedWebGPUDevice: device,
+    print: (s) => { log += s + "\n"; },
+    printErr: (s) => { log += s + "\n"; },
+  });
+
+  postMessage({
+    type: "ready",
+    adapter: `${info.vendor ?? "?"} / ${info.architecture ?? "?"} / ${info.description ?? ""}`,
+  });
+
+  onmessage = async (e) => {
+    const gpuMode = e.data === "gpu" ? 1 : 0;
+    log = "";
+    const t0 = performance.now();
+    let ok = 0;
+    try {
+      // JSPI makes this export return a promise; the GPU waits beneath it
+      // suspend. Nothing in the C++ call chain is written as async.
+      ok = await mod.ccall("run_tests", "number", ["number"], [gpuMode], { async: true });
+    } catch (err) {
+      log += "exception: " + String(err && err.stack || err) + "\n";
+    }
+    postMessage({ type: "done", mode: e.data, ok: !!ok, log, ms: performance.now() - t0 });
+  };
+}
+
+main().catch((e) => postMessage({ type: "fail", msg: String(e && e.stack || e) }));

From f027a2ca1b053caa9e3ab865e73cbb7f858fed89 Mon Sep 17 00:00:00 2001
From: yhirose 
Date: Mon, 20 Jul 2026 00:34:23 -0400
Subject: [PATCH 06/12] webgpu: elementwise, broadcast and row reductions (M10
 phase 3)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Ports the remaining kernel families from metal_kernels.metal: the five
contiguous binaries, the six unaries, the five rank-2 broadcast binaries,
softmax and the row sum/max reductions. gemv_*, attn_decode and rope stay
false stubs, as on Metal.

WGSL has no preprocessor, so the per-op variants Metal generates from a
macro would have to be copy-pasted — the edge-tile bug class that file's
header warns about. Instead each family is one entry point switching on an
op code in the uniform block. The branch is uniform across a dispatch and
these kernels are memory-bound. That in turn lets one bind group layout and
one params struct cover the backend; one-input kernels bind their input
twice, which is legal because both bindings are read-only and the output is
always a fresh allocation.

Two failure modes made this harder than it should have been, and both
presented identically — every op returning zeros, including GEMM, which had
been working:

- CreateShaderModule returns an invalid object rather than null when the
  WGSL does not compile, as does every pipeline built from it, so the
  backend reported itself available and quietly did nothing. It now checks
  GetCompilationInfo and logs the error. The trigger was a literal: WGSL
  has no -inf, and the decimal spelling of f32::lowest rounds out of range.
- Because dispatches batch into one encoder, one invalid dispatch discards
  the entire command buffer, so a new kernel's mistake breaks unrelated
  ops. The harness now installs onuncapturederror and surfaces it.

The per-entry-point dispatch census is kept rather than removed after
debugging: this backend has no test runner that fails when it is absent
(unported ops fall back to CPU and the suite passes either way), so the
counts are the only evidence the GPU did the work.

Browser: 24 cases / 335 assertions green in gpu mode on Chrome, dispatching
sgemm 149, ew_unary 50, ew_binary 36, row_reduce 18, ew_bcast 16,
softmax 8. Native ctest unaffected.
---
 include/webgpu.h                  | 352 +++++++++++++++++++++++-------
 kernels/tensorlib_webgpu.wgsl     | 178 ++++++++++++++-
 kernels/tensorlib_webgpu_wgsl.inc | 178 ++++++++++++++-
 test/wasm/index.html              |   5 +-
 test/wasm/main_wasm.cpp           |   4 +
 test/wasm/worker.js               |  11 +
 6 files changed, 648 insertions(+), 80 deletions(-)

diff --git a/include/webgpu.h b/include/webgpu.h
index 937a737..78f4481 100644
--- a/include/webgpu.h
+++ b/include/webgpu.h
@@ -51,7 +51,9 @@ using kop = tl::metal::kop;
 #include 
 
 #include 
+#include 
 #include 
+#include 
 #include 
 #include 
 
@@ -65,16 +67,54 @@ inline const char* wgsl_source_() {
   return src;
 }
 
-// Uniform params, laid out to match the WGSL Params struct (16 x 4 bytes).
-struct gemm_params {
+// Uniform params, laid out to match the WGSL Params struct. One struct serves
+// every kernel family (see the comment on Params in the .wgsl): unused fields
+// cost a few bytes of a 256-byte slot, and it keeps one uniform ring and one
+// bind group layout for the whole backend.
+struct params {
   uint32_t M, N, K;
   uint32_t lda, ldb, ldc;
   uint32_t a_off, b_off, c_off;
   uint32_t ta, tb;
+  uint32_t ars, acs, brs, bcs;
+  uint32_t op;
   float scale, offset;
-  uint32_t pad0, pad1, pad2;
+  uint32_t pad0, pad1, pad2, pad3, pad4, pad5;
 };
 
+// WGSL gives a uniform-address-space struct align 16, so Params is 96 bytes
+// there. This must agree: the bind group's minBindingSize comes from sizeof
+// here, and a short one fails validation on every dispatch.
+static_assert(sizeof(params) == 96, "params must match the WGSL Params size");
+
+// Which operation within a family, matching the OP_* constants in the WGSL.
+// Families have separate numbering, so this is only meaningful alongside the
+// entry point it is passed to.
+inline uint32_t kernel_op_(kop op) {
+  switch (op) {
+    case kop::add: case kop::badd: return 0;
+    case kop::sub: case kop::bsub: return 1;
+    case kop::mul: case kop::bmul: return 2;
+    case kop::div: case kop::bdiv: return 3;
+    case kop::pow_: case kop::bpow: return 4;
+
+    case kop::exp_: return 0;
+    case kop::log_: return 1;
+    case kop::sqrt_: return 2;
+    case kop::sigmoid: return 3;
+    case kop::relu: return 4;
+    case kop::affine: return 5;
+
+    case kop::row_sum: return 0;
+    case kop::row_max: return 1;
+    default: return 0;
+  }
+}
+
+// WebGPU guarantees maxComputeWorkgroupsPerDimension >= 65535. Anything past
+// that returns false and falls to CPU rather than silently truncating.
+constexpr int64_t kMaxWorkgroups = 65535;
+
 // Dynamic uniform offsets must be a multiple of the adapter's
 // minUniformBufferOffsetAlignment; 256 is the spec's guaranteed-safe maximum.
 constexpr uint64_t kUniformSlot = 256;
@@ -85,7 +125,8 @@ struct context {
   wgpu::Device device;
   wgpu::Queue queue;
   wgpu::BindGroupLayout bgl;
-  wgpu::ComputePipeline sgemm;
+  wgpu::PipelineLayout play;
+  wgpu::ShaderModule mod;
   wgpu::Buffer uniforms;  // ring of kUniformSlots x kUniformSlot bytes
   bool ready = false;
   bool pending = false;
@@ -113,6 +154,17 @@ struct context {
   std::unordered_map>> pool;
   std::unordered_map> staging_pool;
 
+  // Compute pipelines, built on first use and keyed by WGSL entry point.
+  // Pipeline creation is not free, and every entry point shares one layout,
+  // so this is a pure cache.
+  std::unordered_map pipelines;
+  // Per-entry-point dispatch census. Unlike the native backends, this one has
+  // no test runner that fails when it is absent: the browser suite passes
+  // whether or not the GPU engages, because every unported op falls back to
+  // CPU. So the harness reports these counts, and a family reading zero means
+  // the backend quietly stopped doing the work. Cheap enough to always keep.
+  std::unordered_map dispatch_counts;
+
   static context& get() {
     static auto* c = new context();  // leaked: outlives all storage deleters
     return *c;
@@ -140,8 +192,32 @@ struct context {
     wgsl.code = wgsl_source_();
     wgpu::ShaderModuleDescriptor smd = {};
     smd.nextInChain = &wgsl;
-    wgpu::ShaderModule mod = device.CreateShaderModule(&smd);
-    if (!mod) return;
+    mod = device.CreateShaderModule(&smd);
+    // CreateShaderModule hands back an INVALID object, not null, when the WGSL
+    // fails to compile — and so does every pipeline built from it, and every
+    // dispatch then silently does nothing. (That cost a debugging session: one
+    // bad literal made the whole suite read zeros, including ops that had been
+    // working.) Ask for the compilation log and refuse to come up ready.
+    bool compiled = false;
+    if (mod) {
+      wait(mod.GetCompilationInfo(
+          wgpu::CallbackMode::WaitAnyOnly,
+          [&compiled](wgpu::CompilationInfoRequestStatus st,
+                      const wgpu::CompilationInfo* info) {
+            compiled = st == wgpu::CompilationInfoRequestStatus::Success;
+            if (!info) return;
+            for (size_t i = 0; i < info->messageCount; ++i) {
+              const auto& m = info->messages[i];
+              if (m.type != wgpu::CompilationMessageType::Error) continue;
+              compiled = false;
+              std::fprintf(stderr, "tensorlib webgpu: WGSL error at %llu:%llu: %.*s\n",
+                           (unsigned long long)m.lineNum,
+                           (unsigned long long)m.linePos,
+                           (int)m.message.length, m.message.data);
+            }
+          }));
+    }
+    if (!compiled) return;
 
     // Explicit layout rather than GetBindGroupLayout(0): the auto-generated
     // one has no dynamic offset on the uniform binding, which the ring needs.
@@ -156,23 +232,18 @@ struct context {
     be[3].visibility = wgpu::ShaderStage::Compute;
     be[3].buffer.type = wgpu::BufferBindingType::Uniform;
     be[3].buffer.hasDynamicOffset = true;
-    be[3].buffer.minBindingSize = sizeof(gemm_params);
+    be[3].buffer.minBindingSize = sizeof(params);
     wgpu::BindGroupLayoutDescriptor bgld = {};
     bgld.entryCount = 4;
     bgld.entries = be;
     bgl = device.CreateBindGroupLayout(&bgld);
+    if (!bgl) return;
 
     wgpu::PipelineLayoutDescriptor pld = {};
     pld.bindGroupLayoutCount = 1;
     pld.bindGroupLayouts = &bgl;
-    wgpu::PipelineLayout pl = device.CreatePipelineLayout(&pld);
-
-    wgpu::ComputePipelineDescriptor pd = {};
-    pd.layout = pl;
-    pd.compute.module = mod;
-    pd.compute.entryPoint = "sgemm";
-    sgemm = device.CreateComputePipeline(&pd);
-    if (!sgemm) return;
+    play = device.CreatePipelineLayout(&pld);
+    if (!play) return;
 
     wgpu::BufferDescriptor ud = {};
     ud.size = kUniformSlot * kUniformSlots;
@@ -180,9 +251,29 @@ struct context {
     uniforms = device.CreateBuffer(&ud);
     if (!uniforms) return;
 
+    // Build every pipeline up front rather than on first use. Lazily, a WGSL
+    // error would surface as an op quietly falling back to CPU forever; here
+    // it makes available() false, which the harness reports.
+    for (const char* ep : {"sgemm", "ew_binary", "ew_unary", "ew_bcast",
+                           "softmax", "row_reduce"}) {
+      if (!pipeline_(ep)) return;
+    }
+
     ready = true;
   }
 
+  wgpu::ComputePipeline pipeline_(const char* entry) {
+    auto it = pipelines.find(entry);
+    if (it != pipelines.end()) return it->second;
+    wgpu::ComputePipelineDescriptor pd = {};
+    pd.layout = play;
+    pd.compute.module = mod;
+    pd.compute.entryPoint = entry;
+    wgpu::ComputePipeline p = device.CreateComputePipeline(&pd);
+    if (p) pipelines[entry] = p;
+    return p;
+  }
+
   // Blocking wait on a single future — the one place anything suspends.
   bool wait(wgpu::Future f) {
     return instance.WaitAny(f, UINT64_MAX) == wgpu::WaitStatus::Success;
@@ -213,6 +304,50 @@ struct context {
     if (mirror* m = mirror_(native)) m->where = DEVICE;
   }
 
+  // The one place a dispatch is encoded. Every op differs only in which
+  // pipeline, which params and what grid — keeping the bind group, uniform
+  // ring and encoder bookkeeping in a single copy is the same discipline
+  // metal_kernels.metal applies to its kernel bodies.
+  //
+  // `b` may be the same mirror as `a` (a unary or reduce kernel binds its one
+  // input twice): two read-only bindings may alias. `out` is always a fresh
+  // allocation from the evaluator, so a writable binding never does.
+  bool encode_(const char* entry, mirror* a, mirror* b, mirror* out,
+               const params& p, int64_t gx, int64_t gy) {
+    if (gx <= 0 || gy <= 0) return false;
+    if (gx > kMaxWorkgroups || gy > kMaxWorkgroups) return false;
+    wgpu::ComputePipeline pipe = pipeline_(entry);
+    if (!pipe) return false;
+    if (slot >= kUniformSlots) flush_();  // ring exhausted; start a new batch
+
+    const uint32_t off = slot++ * (uint32_t)kUniformSlot;
+    queue.WriteBuffer(uniforms, off, &p, sizeof(p));
+
+    wgpu::BindGroupEntry e[4] = {};
+    e[0].binding = 0; e[0].buffer = a->dev;   e[0].size = a->bytes;
+    e[1].binding = 1; e[1].buffer = b->dev;   e[1].size = b->bytes;
+    e[2].binding = 2; e[2].buffer = out->dev; e[2].size = out->bytes;
+    e[3].binding = 3; e[3].buffer = uniforms; e[3].size = sizeof(params);
+    wgpu::BindGroupDescriptor bgd = {};
+    bgd.layout = bgl;
+    bgd.entryCount = 4;
+    bgd.entries = e;
+    wgpu::BindGroup bg = device.CreateBindGroup(&bgd);
+
+    if (!enc) enc = device.CreateCommandEncoder();
+    wgpu::ComputePassEncoder pass = enc.BeginComputePass();
+    pass.SetPipeline(pipe);
+    pass.SetBindGroup(0, bg, 1, &off);
+    pass.DispatchWorkgroups((uint32_t)gx, (uint32_t)gy, 1);
+    pass.End();
+
+    pending = true;
+    dispatch_counts[entry]++;
+    return true;
+  }
+
+  void flush_();  // defined below, once flush() is in scope
+
   wgpu::Buffer staging_(size_t bytes) {
     auto it = staging_pool.find(bytes);
     if (it != staging_pool.end() && !it->second.empty()) {
@@ -246,6 +381,8 @@ inline void flush() {
       [](wgpu::QueueWorkDoneStatus, wgpu::StringView) {}));
 }
 
+inline void context::flush_() { flush(); }
+
 // Mirror allocation: a device buffer paired with a host buffer (returned via
 // `contents`). They are DISTINCT memory — the dirty state copies between them
 // on demand. The returned handle is an opaque token, not a pointer to
@@ -327,87 +464,152 @@ inline void sync_to_host(void* native, bool for_write) {
   if (for_write) m->where = context::HOST;
 }
 
-// out = op(a) @ op(b) * scale + offset. Offsets are BYTE offsets into the
-// storage; the kernel wants elements, so they must be 4-aligned (they always
-// are for f32 views).
+// Shared host-side prologue for every op: resolve the operand mirrors and
+// stage the lazy copies. Returns false when any operand is untracked — that
+// means a heap storage with no device buffer, so the op belongs on the CPU.
+// `b` may be null for one-input kernels, which then bind `a` twice.
+inline bool operands_(context& c, void* a, void* b, void* out,
+                      context::mirror** ma, context::mirror** mb,
+                      context::mirror** mo) {
+  *ma = c.mirror_(a);
+  *mb = b ? c.mirror_(b) : *ma;
+  *mo = c.mirror_(out);
+  if (!*ma || !*mb || !*mo) return false;
+  c.device_read_(a);
+  if (b) c.device_read_(b);
+  c.device_write_(out);
+  return true;
+}
+
+// Byte offsets must be 4-aligned to convert to the element offsets the
+// kernels index with. They always are for f32 views; anything else falls to
+// the CPU rather than silently truncating.
+inline bool elem_off_(int64_t byte_off, uint32_t* out) {
+  if (byte_off % 4) return false;
+  *out = (uint32_t)(byte_off / 4);
+  return true;
+}
+
+// out = op(a) @ op(b) * scale + offset.
 inline bool gemm(void* a, int64_t ao, int64_t lda, bool ta, void* b, int64_t bo,
                  int64_t ldb, bool tb, void* out, int64_t oo, int64_t m,
                  int64_t n, int64_t k, float scale, float offset) {
   auto& c = context::get();
-  if (!c.ready) return false;
-  if (m <= 0 || n <= 0 || k <= 0) return false;
-  if ((ao % 4) || (bo % 4) || (oo % 4)) return false;
-  // Every operand must be a tracked mirror: an untracked pointer is a heap
-  // storage with no device buffer to bind, so this op belongs on the CPU.
-  context::mirror* ma = c.mirror_(a);
-  context::mirror* mb = c.mirror_(b);
-  context::mirror* mo = c.mirror_(out);
-  if (!ma || !mb || !mo) return false;
-  if (c.slot >= kUniformSlots) flush();  // ring exhausted; start a new batch
-
-  c.device_read_(a);
-  c.device_read_(b);
-  c.device_write_(out);
+  if (!c.ready || m <= 0 || n <= 0 || k <= 0) return false;
+  params p = {};
+  if (!elem_off_(ao, &p.a_off) || !elem_off_(bo, &p.b_off) ||
+      !elem_off_(oo, &p.c_off)) {
+    return false;
+  }
+  context::mirror *ma, *mb, *mo;
+  if (!operands_(c, a, b, out, &ma, &mb, &mo)) return false;
 
-  gemm_params p = {};
   p.M = (uint32_t)m;
   p.N = (uint32_t)n;
   p.K = (uint32_t)k;
   p.lda = (uint32_t)lda;
   p.ldb = (uint32_t)ldb;
   p.ldc = (uint32_t)n;  // the eval seam always hands us a contiguous output
-  p.a_off = (uint32_t)(ao / 4);
-  p.b_off = (uint32_t)(bo / 4);
-  p.c_off = (uint32_t)(oo / 4);
   p.ta = ta ? 1u : 0u;
   p.tb = tb ? 1u : 0u;
   p.scale = scale;
   p.offset = offset;
-
-  const uint32_t off = c.slot++ * (uint32_t)kUniformSlot;
-  c.queue.WriteBuffer(c.uniforms, off, &p, sizeof(p));
-
-  wgpu::BindGroupEntry e[4] = {};
-  e[0].binding = 0; e[0].buffer = ma->dev; e[0].size = ma->bytes;
-  e[1].binding = 1; e[1].buffer = mb->dev; e[1].size = mb->bytes;
-  e[2].binding = 2; e[2].buffer = mo->dev; e[2].size = mo->bytes;
-  e[3].binding = 3; e[3].buffer = c.uniforms; e[3].size = sizeof(gemm_params);
-  wgpu::BindGroupDescriptor bgd = {};
-  bgd.layout = c.bgl;
-  bgd.entryCount = 4;
-  bgd.entries = e;
-  wgpu::BindGroup bg = c.device.CreateBindGroup(&bgd);
-
-  if (!c.enc) c.enc = c.device.CreateCommandEncoder();
-  wgpu::ComputePassEncoder pass = c.enc.BeginComputePass();
-  pass.SetPipeline(c.sgemm);
-  pass.SetBindGroup(0, bg, 1, &off);
-  pass.DispatchWorkgroups((uint32_t)((n + 63) / 64), (uint32_t)((m + 63) / 64),
-                          1);
-  pass.End();
-
-  c.pending = true;
-  return true;
+  return c.encode_("sgemm", ma, mb, mo, p, (n + 63) / 64, (m + 63) / 64);
 }
 
-// ---- Not ported yet (Phase 3+). Returning false routes the op to CPU, which
-// is why each phase lands in a working state.
-inline bool binary(kop, void*, int64_t, void*, int64_t, void*, int64_t, int64_t,
-                   float, float) {
-  return false;
+// Contiguous elementwise binary over n elements.
+inline bool binary(kop op, void* a, int64_t ao, void* b, int64_t bo, void* out,
+                   int64_t oo, int64_t n, float scale, float offset) {
+  auto& c = context::get();
+  if (!c.ready || n <= 0) return false;
+  params p = {};
+  if (!elem_off_(ao, &p.a_off) || !elem_off_(bo, &p.b_off) ||
+      !elem_off_(oo, &p.c_off)) {
+    return false;
+  }
+  context::mirror *ma, *mb, *mo;
+  if (!operands_(c, a, b, out, &ma, &mb, &mo)) return false;
+
+  p.M = (uint32_t)n;
+  p.op = kernel_op_(op);
+  p.scale = scale;
+  p.offset = offset;
+  return c.encode_("ew_binary", ma, mb, mo, p, (n + 255) / 256, 1);
 }
-inline bool binary_bcast(kop, void*, int64_t, int64_t, int64_t, void*, int64_t,
-                         int64_t, int64_t, void*, int64_t, int64_t, int64_t,
-                         float, float) {
-  return false;
+
+inline bool unary(kop op, void* a, int64_t ao, void* out, int64_t oo, int64_t n,
+                  float scale, float offset) {
+  auto& c = context::get();
+  if (!c.ready || n <= 0) return false;
+  params p = {};
+  if (!elem_off_(ao, &p.a_off) || !elem_off_(oo, &p.c_off)) return false;
+  context::mirror *ma, *mb, *mo;
+  if (!operands_(c, a, nullptr, out, &ma, &mb, &mo)) return false;
+
+  p.b_off = p.a_off;  // the kernel binds its one input twice
+  p.M = (uint32_t)n;
+  p.op = kernel_op_(op);
+  p.scale = scale;
+  p.offset = offset;
+  return c.encode_("ew_unary", ma, mb, mo, p, (n + 255) / 256, 1);
 }
-inline bool unary(kop, void*, int64_t, void*, int64_t, int64_t, float, float) {
-  return false;
+
+// Rank-2 broadcast binary: out[r,c] = f(a[r*ars + c*acs], b[r*brs + c*bcs])
+// into a contiguous [m,n] output. One stride-parameterized kernel covers every
+// rank-2 broadcast, which keeps bias/gamma/beta chains on the GPU — falling
+// back mid-graph would cost a full submit-and-wait.
+inline bool binary_bcast(kop op, void* a, int64_t ao, int64_t ars, int64_t acs,
+                         void* b, int64_t bo, int64_t brs, int64_t bcs,
+                         void* out, int64_t oo, int64_t m, int64_t n,
+                         float scale, float offset) {
+  auto& c = context::get();
+  if (!c.ready || m <= 0 || n <= 0) return false;
+  // Broadcast strides are non-negative here (broadcast_strides only ever
+  // zeroes an axis); a negative one would wrap as u32 in the kernel.
+  if (ars < 0 || acs < 0 || brs < 0 || bcs < 0) return false;
+  params p = {};
+  if (!elem_off_(ao, &p.a_off) || !elem_off_(bo, &p.b_off) ||
+      !elem_off_(oo, &p.c_off)) {
+    return false;
+  }
+  context::mirror *ma, *mb, *mo;
+  if (!operands_(c, a, b, out, &ma, &mb, &mo)) return false;
+
+  p.M = (uint32_t)m;
+  p.N = (uint32_t)n;
+  p.ars = (uint32_t)ars;
+  p.acs = (uint32_t)acs;
+  p.brs = (uint32_t)brs;
+  p.bcs = (uint32_t)bcs;
+  p.op = kernel_op_(op);
+  p.scale = scale;
+  p.offset = offset;
+  return c.encode_("ew_bcast", ma, mb, mo, p, (n + 31) / 32, (m + 7) / 8);
 }
-inline bool row_op(kop, void*, int64_t, void*, int64_t, int64_t, int64_t, float,
-                   float) {
-  return false;
+
+// Row-wise op over the last axis: softmax writes rows x cols; row_sum/row_max
+// write one value per row, with the affine epilogue. One workgroup per row.
+inline bool row_op(kop op, void* in, int64_t io, void* out, int64_t oo,
+                   int64_t rows, int64_t cols, float scale, float offset) {
+  auto& c = context::get();
+  if (!c.ready || rows <= 0 || cols <= 0) return false;
+  params p = {};
+  if (!elem_off_(io, &p.a_off) || !elem_off_(oo, &p.c_off)) return false;
+  context::mirror *ma, *mb, *mo;
+  if (!operands_(c, in, nullptr, out, &ma, &mb, &mo)) return false;
+
+  p.b_off = p.a_off;
+  p.M = (uint32_t)rows;
+  p.N = (uint32_t)cols;
+  p.op = kernel_op_(op);
+  p.scale = scale;
+  p.offset = offset;
+  const char* entry = op == kop::softmax ? "softmax" : "row_reduce";
+  return c.encode_(entry, ma, mb, mo, p, rows, 1);
 }
+
+// ---- Not ported yet. Returning false routes the op to CPU, which is why each
+// phase lands in a working state.
 inline bool gemv_f32(void*, void*, void*, int64_t, int64_t) { return false; }
 inline bool gemv_bf16(void*, void*, void*, int64_t, int64_t) { return false; }
 inline bool attn_decode(void*, void*, void*, void*, int64_t, int64_t, int64_t,
diff --git a/kernels/tensorlib_webgpu.wgsl b/kernels/tensorlib_webgpu.wgsl
index c18398e..5386ae6 100644
--- a/kernels/tensorlib_webgpu.wgsl
+++ b/kernels/tensorlib_webgpu.wgsl
@@ -12,9 +12,16 @@
 // not. So every binding covers its whole buffer. (CUDA instead folds offsets
 // host-side into the pointer it passes, which WebGPU has no equivalent of.)
 
+// One Params struct and one bind group layout serve every kernel here, rather
+// than the per-family structs metal_kernels.metal uses. WebGPU's bind group
+// ceremony is heavy enough that a second layout would buy nothing: the fields
+// each family ignores cost 4 bytes of a 256-byte uniform slot. Kernels that
+// take one input (unary, the row reductions) get A bound to B as well — two
+// read-only bindings may alias, and the output is always a fresh allocation,
+// so no writable binding ever aliases a readable one.
 struct Params {
-  M : u32,
-  N : u32,
+  M : u32,        // gemm rows | elementwise element count | reduce rows
+  N : u32,        // gemm cols | reduce cols
   K : u32,
   lda : u32,
   ldb : u32,
@@ -24,11 +31,23 @@ struct Params {
   c_off : u32,
   ta : u32,
   tb : u32,
+  ars : u32,      // broadcast: per-operand row/col strides, in elements
+  acs : u32,
+  brs : u32,
+  bcs : u32,
+  op : u32,       // which operation, within the entry point's family
   scale : f32,
   offset : f32,
+  // A uniform-address-space struct has align 16, so its size rounds up to a
+  // multiple of 16. Pad explicitly to 96 bytes so the host struct (which the
+  // bind group's minBindingSize comes from) matches exactly — a short
+  // minBindingSize fails bind group validation for every dispatch.
   _pad0 : u32,
   _pad1 : u32,
   _pad2 : u32,
+  _pad3 : u32,
+  _pad4 : u32,
+  _pad5 : u32,
 };
 
 @group(0) @binding(0) var       A : array;
@@ -128,3 +147,158 @@ fn sgemm(@builtin(workgroup_id) wg : vec3,
     }
   }
 }
+
+// ---- Elementwise, broadcast and row reductions
+//
+// WGSL has neither templates nor a preprocessor, so the per-op variants that
+// metal_kernels.metal generates from a macro would have to be copy-pasted
+// here — exactly the edge-tile bug class that file's header warns against.
+// Instead the operation is a uniform field and each family is ONE entry point
+// that switches on it. The branch is uniform across the dispatch and these
+// kernels are memory-bound, so it costs nothing measurable; what it buys is a
+// single copy of every bounds check and epilogue.
+//
+// Op codes are assigned by kernel_op_() in webgpu.h.
+const OP_ADD : u32 = 0u;
+const OP_SUB : u32 = 1u;
+const OP_MUL : u32 = 2u;
+const OP_DIV : u32 = 3u;
+const OP_POW : u32 = 4u;
+
+const OP_EXP     : u32 = 0u;
+const OP_LOG     : u32 = 1u;
+const OP_SQRT    : u32 = 2u;
+const OP_SIGMOID : u32 = 3u;
+const OP_RELU    : u32 = 4u;
+const OP_AFFINE  : u32 = 5u;
+
+const OP_ROW_SUM : u32 = 0u;
+const OP_ROW_MAX : u32 = 1u;
+
+// Identity for a max reduction. WGSL has no -inf literal, and the decimal
+// spelling of f32::lowest rounds just past the representable range ("cannot be
+// represented as 'f32'"), so this is the largest round number safely inside
+// it. Threads whose row is shorter than the workgroup contribute this.
+const NEG_HUGE : f32 = -3.4e38;
+
+// Shared by the contiguous and the broadcast binary: same five operations,
+// only the addressing differs.
+fn binary_op(op : u32, av : f32, bv : f32) -> f32 {
+  switch (op) {
+    case 1u: { return av - bv; }
+    case 2u: { return av * bv; }
+    case 3u: { return av / bv; }
+    case 4u: { return pow(av, bv); }
+    default: { return av + bv; }
+  }
+}
+
+fn unary_op(op : u32, v : f32) -> f32 {
+  switch (op) {
+    case 1u: { return log(v); }
+    case 2u: { return sqrt(v); }
+    case 3u: { return 1.0 / (1.0 + exp(-v)); }
+    case 4u: { return max(v, 0.0); }
+    case 5u: { return v; }
+    default: { return exp(v); }
+  }
+}
+
+// Contiguous elementwise binary over p.M elements.
+@compute @workgroup_size(256, 1, 1)
+fn ew_binary(@builtin(global_invocation_id) gid : vec3) {
+  let i = gid.x;
+  if (i >= p.M) { return; }
+  let v = binary_op(p.op, A[p.a_off + i], B[p.b_off + i]);
+  C[p.c_off + i] = fma(v, p.scale, p.offset);
+}
+
+@compute @workgroup_size(256, 1, 1)
+fn ew_unary(@builtin(global_invocation_id) gid : vec3) {
+  let i = gid.x;
+  if (i >= p.M) { return; }
+  let v = unary_op(p.op, A[p.a_off + i]);
+  C[p.c_off + i] = fma(v, p.scale, p.offset);
+}
+
+// Rank-2 broadcast binary: out[r,c] = f(a[r*ars + c*acs], b[r*brs + c*bcs])
+// into a contiguous [M,N] output. Per-operand strides express every rank-2
+// broadcast (row vector, column vector, scalar) in one kernel, which keeps
+// bias/gamma/beta chains on the GPU — a CPU fallback mid-graph costs a full
+// submit-and-wait, and that is dearer here than on Metal.
+@compute @workgroup_size(32, 8, 1)
+fn ew_bcast(@builtin(global_invocation_id) gid : vec3) {
+  let c = gid.x;
+  let r = gid.y;
+  if (c >= p.N || r >= p.M) { return; }
+  let av = A[p.a_off + r * p.ars + c * p.acs];
+  let bv = B[p.b_off + r * p.brs + c * p.bcs];
+  C[p.c_off + r * p.N + c] = fma(binary_op(p.op, av, bv), p.scale, p.offset);
+}
+
+// ---- Row reductions over the last axis: one workgroup per row, 256
+// invocations, workgroup-scratch tree reduction. p.N (cols) may exceed the
+// invocation count, so each thread strides over the row first.
+
+const T : u32 = 256u;
+var scratch : array;
+
+// Numerically stable softmax (subtract the row max). Applying an affine
+// epilogue to a softmax is not meaningful, so scale/offset are ignored here,
+// as they are on Metal.
+@compute @workgroup_size(256, 1, 1)
+fn softmax(@builtin(workgroup_id) wg : vec3,
+           @builtin(local_invocation_index) lid : u32) {
+  let row = wg.x;
+  let src = p.a_off + row * p.N;
+  let dst = p.c_off + row * p.N;
+
+  var m = NEG_HUGE;
+  for (var c : u32 = lid; c < p.N; c = c + T) { m = max(m, A[src + c]); }
+  scratch[lid] = m;
+  workgroupBarrier();
+  for (var s : u32 = T / 2u; s > 0u; s = s >> 1u) {
+    if (lid < s) { scratch[lid] = max(scratch[lid], scratch[lid + s]); }
+    workgroupBarrier();
+  }
+  let row_max = scratch[0];
+  workgroupBarrier();
+
+  var sum = 0.0;
+  for (var c : u32 = lid; c < p.N; c = c + T) { sum = sum + exp(A[src + c] - row_max); }
+  scratch[lid] = sum;
+  workgroupBarrier();
+  for (var s : u32 = T / 2u; s > 0u; s = s >> 1u) {
+    if (lid < s) { scratch[lid] = scratch[lid] + scratch[lid + s]; }
+    workgroupBarrier();
+  }
+  let inv = 1.0 / scratch[0];
+  for (var c : u32 = lid; c < p.N; c = c + T) {
+    C[dst + c] = exp(A[src + c] - row_max) * inv;
+  }
+}
+
+fn reduce_op(op : u32, acc : f32, v : f32) -> f32 {
+  if (op == OP_ROW_MAX) { return max(acc, v); }
+  return acc + v;
+}
+
+// One value per row, with the affine epilogue.
+@compute @workgroup_size(256, 1, 1)
+fn row_reduce(@builtin(workgroup_id) wg : vec3,
+              @builtin(local_invocation_index) lid : u32) {
+  let row = wg.x;
+  let src = p.a_off + row * p.N;
+
+  var acc = select(0.0, NEG_HUGE, p.op == OP_ROW_MAX);
+  for (var c : u32 = lid; c < p.N; c = c + T) {
+    acc = reduce_op(p.op, acc, A[src + c]);
+  }
+  scratch[lid] = acc;
+  workgroupBarrier();
+  for (var s : u32 = T / 2u; s > 0u; s = s >> 1u) {
+    if (lid < s) { scratch[lid] = reduce_op(p.op, scratch[lid], scratch[lid + s]); }
+    workgroupBarrier();
+  }
+  if (lid == 0u) { C[p.c_off + row] = scratch[0] * p.scale + p.offset; }
+}
diff --git a/kernels/tensorlib_webgpu_wgsl.inc b/kernels/tensorlib_webgpu_wgsl.inc
index 2cd0259..1cb2702 100644
--- a/kernels/tensorlib_webgpu_wgsl.inc
+++ b/kernels/tensorlib_webgpu_wgsl.inc
@@ -14,9 +14,16 @@ R"WGSL(
 // not. So every binding covers its whole buffer. (CUDA instead folds offsets
 // host-side into the pointer it passes, which WebGPU has no equivalent of.)
 
+// One Params struct and one bind group layout serve every kernel here, rather
+// than the per-family structs metal_kernels.metal uses. WebGPU's bind group
+// ceremony is heavy enough that a second layout would buy nothing: the fields
+// each family ignores cost 4 bytes of a 256-byte uniform slot. Kernels that
+// take one input (unary, the row reductions) get A bound to B as well — two
+// read-only bindings may alias, and the output is always a fresh allocation,
+// so no writable binding ever aliases a readable one.
 struct Params {
-  M : u32,
-  N : u32,
+  M : u32,        // gemm rows | elementwise element count | reduce rows
+  N : u32,        // gemm cols | reduce cols
   K : u32,
   lda : u32,
   ldb : u32,
@@ -26,11 +33,23 @@ struct Params {
   c_off : u32,
   ta : u32,
   tb : u32,
+  ars : u32,      // broadcast: per-operand row/col strides, in elements
+  acs : u32,
+  brs : u32,
+  bcs : u32,
+  op : u32,       // which operation, within the entry point's family
   scale : f32,
   offset : f32,
+  // A uniform-address-space struct has align 16, so its size rounds up to a
+  // multiple of 16. Pad explicitly to 96 bytes so the host struct (which the
+  // bind group's minBindingSize comes from) matches exactly — a short
+  // minBindingSize fails bind group validation for every dispatch.
   _pad0 : u32,
   _pad1 : u32,
   _pad2 : u32,
+  _pad3 : u32,
+  _pad4 : u32,
+  _pad5 : u32,
 };
 
 @group(0) @binding(0) var       A : array;
@@ -130,4 +149,159 @@ fn sgemm(@builtin(workgroup_id) wg : vec3,
     }
   }
 }
+
+// ---- Elementwise, broadcast and row reductions
+//
+// WGSL has neither templates nor a preprocessor, so the per-op variants that
+// metal_kernels.metal generates from a macro would have to be copy-pasted
+// here — exactly the edge-tile bug class that file's header warns against.
+// Instead the operation is a uniform field and each family is ONE entry point
+// that switches on it. The branch is uniform across the dispatch and these
+// kernels are memory-bound, so it costs nothing measurable; what it buys is a
+// single copy of every bounds check and epilogue.
+//
+// Op codes are assigned by kernel_op_() in webgpu.h.
+const OP_ADD : u32 = 0u;
+const OP_SUB : u32 = 1u;
+const OP_MUL : u32 = 2u;
+const OP_DIV : u32 = 3u;
+const OP_POW : u32 = 4u;
+
+const OP_EXP     : u32 = 0u;
+const OP_LOG     : u32 = 1u;
+const OP_SQRT    : u32 = 2u;
+const OP_SIGMOID : u32 = 3u;
+const OP_RELU    : u32 = 4u;
+const OP_AFFINE  : u32 = 5u;
+
+const OP_ROW_SUM : u32 = 0u;
+const OP_ROW_MAX : u32 = 1u;
+
+// Identity for a max reduction. WGSL has no -inf literal, and the decimal
+// spelling of f32::lowest rounds just past the representable range ("cannot be
+// represented as 'f32'"), so this is the largest round number safely inside
+// it. Threads whose row is shorter than the workgroup contribute this.
+const NEG_HUGE : f32 = -3.4e38;
+
+// Shared by the contiguous and the broadcast binary: same five operations,
+// only the addressing differs.
+fn binary_op(op : u32, av : f32, bv : f32) -> f32 {
+  switch (op) {
+    case 1u: { return av - bv; }
+    case 2u: { return av * bv; }
+    case 3u: { return av / bv; }
+    case 4u: { return pow(av, bv); }
+    default: { return av + bv; }
+  }
+}
+
+fn unary_op(op : u32, v : f32) -> f32 {
+  switch (op) {
+    case 1u: { return log(v); }
+    case 2u: { return sqrt(v); }
+    case 3u: { return 1.0 / (1.0 + exp(-v)); }
+    case 4u: { return max(v, 0.0); }
+    case 5u: { return v; }
+    default: { return exp(v); }
+  }
+}
+
+// Contiguous elementwise binary over p.M elements.
+@compute @workgroup_size(256, 1, 1)
+fn ew_binary(@builtin(global_invocation_id) gid : vec3) {
+  let i = gid.x;
+  if (i >= p.M) { return; }
+  let v = binary_op(p.op, A[p.a_off + i], B[p.b_off + i]);
+  C[p.c_off + i] = fma(v, p.scale, p.offset);
+}
+
+@compute @workgroup_size(256, 1, 1)
+fn ew_unary(@builtin(global_invocation_id) gid : vec3) {
+  let i = gid.x;
+  if (i >= p.M) { return; }
+  let v = unary_op(p.op, A[p.a_off + i]);
+  C[p.c_off + i] = fma(v, p.scale, p.offset);
+}
+
+// Rank-2 broadcast binary: out[r,c] = f(a[r*ars + c*acs], b[r*brs + c*bcs])
+// into a contiguous [M,N] output. Per-operand strides express every rank-2
+// broadcast (row vector, column vector, scalar) in one kernel, which keeps
+// bias/gamma/beta chains on the GPU — a CPU fallback mid-graph costs a full
+// submit-and-wait, and that is dearer here than on Metal.
+@compute @workgroup_size(32, 8, 1)
+fn ew_bcast(@builtin(global_invocation_id) gid : vec3) {
+  let c = gid.x;
+  let r = gid.y;
+  if (c >= p.N || r >= p.M) { return; }
+  let av = A[p.a_off + r * p.ars + c * p.acs];
+  let bv = B[p.b_off + r * p.brs + c * p.bcs];
+  C[p.c_off + r * p.N + c] = fma(binary_op(p.op, av, bv), p.scale, p.offset);
+}
+
+// ---- Row reductions over the last axis: one workgroup per row, 256
+// invocations, workgroup-scratch tree reduction. p.N (cols) may exceed the
+// invocation count, so each thread strides over the row first.
+
+const T : u32 = 256u;
+var scratch : array;
+
+// Numerically stable softmax (subtract the row max). Applying an affine
+// epilogue to a softmax is not meaningful, so scale/offset are ignored here,
+// as they are on Metal.
+@compute @workgroup_size(256, 1, 1)
+fn softmax(@builtin(workgroup_id) wg : vec3,
+           @builtin(local_invocation_index) lid : u32) {
+  let row = wg.x;
+  let src = p.a_off + row * p.N;
+  let dst = p.c_off + row * p.N;
+
+  var m = NEG_HUGE;
+  for (var c : u32 = lid; c < p.N; c = c + T) { m = max(m, A[src + c]); }
+  scratch[lid] = m;
+  workgroupBarrier();
+  for (var s : u32 = T / 2u; s > 0u; s = s >> 1u) {
+    if (lid < s) { scratch[lid] = max(scratch[lid], scratch[lid + s]); }
+    workgroupBarrier();
+  }
+  let row_max = scratch[0];
+  workgroupBarrier();
+
+  var sum = 0.0;
+  for (var c : u32 = lid; c < p.N; c = c + T) { sum = sum + exp(A[src + c] - row_max); }
+  scratch[lid] = sum;
+  workgroupBarrier();
+  for (var s : u32 = T / 2u; s > 0u; s = s >> 1u) {
+    if (lid < s) { scratch[lid] = scratch[lid] + scratch[lid + s]; }
+    workgroupBarrier();
+  }
+  let inv = 1.0 / scratch[0];
+  for (var c : u32 = lid; c < p.N; c = c + T) {
+    C[dst + c] = exp(A[src + c] - row_max) * inv;
+  }
+}
+
+fn reduce_op(op : u32, acc : f32, v : f32) -> f32 {
+  if (op == OP_ROW_MAX) { return max(acc, v); }
+  return acc + v;
+}
+
+// One value per row, with the affine epilogue.
+@compute @workgroup_size(256, 1, 1)
+fn row_reduce(@builtin(workgroup_id) wg : vec3,
+              @builtin(local_invocation_index) lid : u32) {
+  let row = wg.x;
+  let src = p.a_off + row * p.N;
+
+  var acc = select(0.0, NEG_HUGE, p.op == OP_ROW_MAX);
+  for (var c : u32 = lid; c < p.N; c = c + T) {
+    acc = reduce_op(p.op, acc, A[src + c]);
+  }
+  scratch[lid] = acc;
+  workgroupBarrier();
+  for (var s : u32 = T / 2u; s > 0u; s = s >> 1u) {
+    if (lid < s) { scratch[lid] = reduce_op(p.op, scratch[lid], scratch[lid + s]); }
+    workgroupBarrier();
+  }
+  if (lid == 0u) { C[p.c_off + row] = scratch[0] * p.scale + p.offset; }
+}
 )WGSL"
diff --git a/test/wasm/index.html b/test/wasm/index.html
index e0155a9..1c855cd 100644
--- a/test/wasm/index.html
+++ b/test/wasm/index.html
@@ -33,7 +33,10 @@ 
cpp-tensorlib WebGPU tests

     } else if (m.type === "done") {
       status.textContent = m.ok ? "OK" : "FAILED";
       out.textContent = m.log;
-      console.log("[tltest] " + (m.ok ? "OK" : "FAILED") + " in " + Math.round(m.ms) + "ms\n" + m.log);
+      // Console gets a compact one-liner (the full report is in the page):
+      // a browser-automation check reads the verdict, not the transcript.
+      const tail = m.log.split("\n").filter((l) => /test cases:|assertions:|dispatch |FAIL:|WEBGPU|ERROR:.*validat|is invalid|DEVICE LOST/.test(l));
+      console.log("[tltest] " + (m.ok ? "OK" : "FAILED") + " in " + Math.round(m.ms) + "ms | " + tail.join(" | ").replace(/\u001b\[[0-9;]*m/g, ""));
       console.log("[tltest] DONE");
     }
   };
diff --git a/test/wasm/main_wasm.cpp b/test/wasm/main_wasm.cpp
index a6b415f..687c071 100644
--- a/test/wasm/main_wasm.cpp
+++ b/test/wasm/main_wasm.cpp
@@ -46,6 +46,10 @@ EMSCRIPTEN_KEEPALIVE int run_tests(int gpu_mode) {
 
   doctest::Context ctx;
   int failed = ctx.run();
+  // A green suite does not prove the backend ran — unported ops fall back to
+  // CPU and pass either way. Report which kernels actually dispatched.
+  for (auto& kv : tl::webgpu::context::get().dispatch_counts)
+    std::printf("[tensorlib_wasm] dispatch %s=%ld\n", kv.first.c_str(), kv.second);
   std::printf("[tensorlib_wasm] %s\n", failed ? "FAIL" : "OK");
   return failed ? 0 : 1;
 }
diff --git a/test/wasm/worker.js b/test/wasm/worker.js
index c887db9..8346373 100644
--- a/test/wasm/worker.js
+++ b/test/wasm/worker.js
@@ -16,6 +16,16 @@ async function main() {
     return;
   }
   const device = await adapter.requestDevice();
+
+  // Surface WebGPU validation errors. Without this they are only warnings in
+  // the worker console, and a single invalid dispatch silently discards the
+  // whole command buffer it was batched into — which reads as "every op
+  // returned zeros", pointing nowhere near the actual mistake.
+  const gpuErrors = [];
+  device.onuncapturederror = (e) => {
+    if (gpuErrors.length < 5) gpuErrors.push(String(e.error.message).split("\n")[0]);
+  };
+  device.lost.then((info) => gpuErrors.push("DEVICE LOST: " + info.message));
   const info = adapter.info ?? {};
 
   // doctest writes its report to stdout; capture it rather than let it vanish
@@ -43,6 +53,7 @@ async function main() {
     } catch (err) {
       log += "exception: " + String(err && err.stack || err) + "\n";
     }
+    if (gpuErrors.length) log = "WEBGPU ERRORS:\n" + gpuErrors.join("\n") + "\n\n" + log;
     postMessage({ type: "done", mode: e.data, ok: !!ok, log, ms: performance.now() - t0 });
   };
 }

From 2ce5b7a5af1db5c14bcc27bd2973a74f2c754244 Mon Sep 17 00:00:00 2001
From: yhirose 
Date: Mon, 20 Jul 2026 00:54:26 -0400
Subject: [PATCH 07/12] webgpu: calibrate the auto-mode thresholds (M10 phase
 4)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Port misc/census.cpp to the browser as its own wasm and page, measure the
CPU/GPU crossovers, and fill in the TENSORLIB_WEBGPU arms of
auto_threshold_() and batch_matmul_bias_threshold_() — until now the wasm
build was using the Metal numbers, which are calibrated against an AMX CPU
rival that does not exist here.

Measured on Chrome / M1 Pro: matmul 4e6 in num_elements*k units (two orders
below Metal's 5e8), elementwise 2e6, reduction 2e5, batch bias 8e6. The
elementwise and reduction numbers coincide with Metal's because WebGPU's
much worse dispatch floor and the much weaker scalar-wasm CPU baseline
cancel; matmul is where they do not.

The census gets a whole-block section so the batch bias is measured directly
rather than inferred, and the test harness gains an auto-mode pass — native
ctest cannot reach this arm, so that pass is the only coverage the thresholds
have. Both stay green, with auto dispatching 57 of gpu mode's 149 sgemms.
---
 include/types.h            |  47 ++++++++--
 test/wasm/build.sh         |  62 ++++++++++----
 test/wasm/census-worker.js |  53 ++++++++++++
 test/wasm/census.html      |  47 ++++++++++
 test/wasm/census_wasm.cpp  | 171 +++++++++++++++++++++++++++++++++++++
 test/wasm/index.html       |  32 ++++---
 test/wasm/main_wasm.cpp    |  31 +++++--
 test/wasm/worker.js        |   3 +-
 8 files changed, 404 insertions(+), 42 deletions(-)
 create mode 100644 test/wasm/census-worker.js
 create mode 100644 test/wasm/census.html
 create mode 100644 test/wasm/census_wasm.cpp

diff --git a/include/types.h b/include/types.h
index 1f9e979..9bd8956 100644
--- a/include/types.h
+++ b/include/types.h
@@ -76,13 +76,40 @@ enum class kernel_class {
 // auto_-mode size thresholds: the measured CPU vs GPU-total (single op +
 // flush) standalone crossover per family — the size where GPU-total first
 // beats CPU. The GPU backend is fixed at build time by the gpu:: facade
-// (CUDA off-Apple when TENSORLIB_CUDA, else Metal on Apple), so the crossover
-// is a compile-time property too and the thresholds branch on the same macro.
+// (WebGPU when TENSORLIB_WEBGPU, CUDA off-Apple when TENSORLIB_CUDA, else
+// Metal on Apple), so the crossover is a compile-time property too and the
+// thresholds branch on the same macros, in the same order as gpu.h.
 // These are the conservative standalone values; pipelined graphs amortize the
 // flush, so the effective crossover is lower. Re-measure per target with
-// misc/census.cpp.
+// misc/census.cpp — or, for WebGPU, its browser port
+// test/wasm/census_wasm.cpp, since that backend only exists in a page.
 inline int64_t auto_threshold_(kernel_class kc) {
-#if defined(TENSORLIB_CUDA)
+#if defined(TENSORLIB_WEBGPU)
+  // WebGPU / Chrome on an M1 Pro (apple/metal-3), census 2026-07-20 via
+  // test/wasm/census_wasm.cpp. Two facts shape this arm, and only one of them
+  // is about the GPU:
+  //   - There is a ~0.3-0.6 ms floor per dispatch+flush, orders of magnitude
+  //     above Metal's, so nothing small can ever win.
+  //   - The wasm CPU rival is scalar and single-threaded (~8 GF/s measured at
+  //     512^3), roughly 30x weaker than the Accelerate/AMX path the Metal arm
+  //     below is calibrated against.
+  // Those pull in opposite directions and, for the memory-bound families,
+  // cancel almost exactly — elementwise and reduction land on the same numbers
+  // as Metal. Matmul does not: the CPU side loses its AMX rival while the WGSL
+  // sgemm still reaches ~160 GF/s, so the crossover falls two orders of
+  // magnitude, from 5e8 to 4e6.
+  //   matmul      128^3=2.1e6 cpu, 160^3=4.1e6 a tie in both runs,
+  //               192^3=7.1e6 GPU clearly (1.9 vs 1.1 ms) -> ~4e6.
+  //   elementwise 1M cpu/tie (0.5 vs 0.5-0.7 ms), 4M GPU (2.0 vs 1.4)
+  //               -> ~2e6.
+  //   reduction   softmax 65536 cpu at all three row widths, 262144 GPU at all
+  //               three (1.3 vs 0.8 ms) -> ~2e5.
+  switch (kc) {
+    case kernel_class::matmul: return 4'000'000;        // ~160^3
+    case kernel_class::elementwise: return 2'000'000;   // 2M elements
+    case kernel_class::reduction: return 200'000;       // ~2e5 elements
+  }
+#elif defined(TENSORLIB_CUDA)
   // RTX 3090 (sm_86) census 2026-07-04, misc/census.cpp. No AMX rival on the
   // CPU side (own BLIS), so the GPU wins far earlier than on Apple's Metal:
   //   matmul      64^3=2.6e5 already GPU (0.038 vs cpu 0.152 ms) — crossover
@@ -138,7 +165,17 @@ inline int64_t batch_matmul_bias_threshold_() {
       long long x = std::strtoll(e, &end, 10);
       if (end != e && x > 0) return static_cast(x);
     }
-#if defined(TENSORLIB_CUDA)
+#if defined(TENSORLIB_WEBGPU)
+    // WebGPU / Chrome on an M1 Pro, census 2026-07-20. Measured directly with
+    // a transformer-shaped eval batch (qkv + attention + FFN in one graph),
+    // which is what this threshold is actually about — the Metal arm below had
+    // to infer it from a separate pipelined bench. Sum of M*N*K over the block:
+    // d=32 (2.5e6) cpu, d=48 (4.8e6) a tie across runs, d=64 (7.9e6) GPU in
+    // both (2.8 vs 1.8 ms) -> 8e6. It sits just above the per-op matmul
+    // crossover, as expected: the block's FFN gemm dominates it, so a batch
+    // earns the GPU at roughly the size its largest gemm does.
+    return 8'000'000;
+#elif defined(TENSORLIB_CUDA)
     return 4'000'000;             // GPU wins early; a couple of small gemms
 #else
     // Metal / M1 Pro, calibrated 2026-07-18 on the pipelined transformer bench
diff --git a/test/wasm/build.sh b/test/wasm/build.sh
index 41f5e0e..bea305f 100755
--- a/test/wasm/build.sh
+++ b/test/wasm/build.sh
@@ -1,16 +1,23 @@
 #!/usr/bin/env bash
-# Build the browser test harness for the WebGPU backend.
+# Build the browser harnesses for the WebGPU backend.
 #
 #   ./build.sh && python3 -m http.server 8731
-#   # open http://localhost:8731/test/wasm/site/index.html
+#   # tests:  http://localhost:8731/test/wasm/site/index.html
+#   # census: http://localhost:8731/test/wasm/site/census.html
+#
+# Two separate wasms, not one with two entry points: the test suite has to stay
+# a fast pass/fail and the census runs for minutes.
 #
 # A flat emcc line, not CMake — which is why the WGSL .inc is committed rather
 # than generated at build time (see kernels/gen_wgsl_inc.sh).
+#
+# Pass `tests` or `census` to build just one.
 set -euo pipefail
 
 HERE="$(cd "$(dirname "$0")" && pwd)"
 ROOT="$(cd "$HERE/../.." && pwd)"
 OUT="$HERE/site"
+WHAT="${1:-all}"
 
 source "${EMSDK_DIR:-$HOME/Projects/emsdk}/emsdk_env.sh" >/dev/null 2>&1
 
@@ -20,20 +27,39 @@ mkdir -p "$OUT"
 # culebra's interpreter needs). The cost is that this wasm *requires* JSPI, so
 # it runs on Chrome only — acceptable for a test harness, and the shipping
 # playground pairs it with a CPU build behind a feature check.
-set -x
-emcc -std=c++23 -O2 -fwasm-exceptions \
-  --use-port=emdawnwebgpu \
-  -DTENSORLIB_WEBGPU \
-  -sJSPI=1 -sJSPI_EXPORTS=run_tests \
-  -sSTACK_SIZE=16MB -sALLOW_MEMORY_GROWTH=1 -sINITIAL_MEMORY=64MB \
-  -sMODULARIZE=1 -sEXPORT_ES6=1 -sEXPORT_NAME=createTests \
-  -sENVIRONMENT=web,worker \
-  -sEXPORTED_FUNCTIONS=_run_tests,_malloc,_free \
-  -sEXPORTED_RUNTIME_METHODS=ccall,cwrap,UTF8ToString \
-  -I"$ROOT/include" -I"$ROOT/kernels" -I"$ROOT/test" \
-  "$HERE/main_wasm.cpp" "$ROOT/test/test_array.cpp" \
-  -o "$OUT/tests.js"
-set +x
+common=(
+  -std=c++23 -O2 -fwasm-exceptions
+  --use-port=emdawnwebgpu
+  -DTENSORLIB_WEBGPU
+  -sSTACK_SIZE=16MB -sALLOW_MEMORY_GROWTH=1 -sINITIAL_MEMORY=64MB
+  -sMODULARIZE=1 -sEXPORT_ES6=1
+  -sENVIRONMENT=web,worker
+  -sEXPORTED_RUNTIME_METHODS=ccall,cwrap,UTF8ToString
+  -I"$ROOT/include" -I"$ROOT/kernels" -I"$ROOT/test"
+)
+
+if [[ "$WHAT" == all || "$WHAT" == tests ]]; then
+  set -x
+  emcc "${common[@]}" \
+    -sJSPI=1 -sJSPI_EXPORTS=run_tests \
+    -sEXPORT_NAME=createTests \
+    -sEXPORTED_FUNCTIONS=_run_tests,_malloc,_free \
+    "$HERE/main_wasm.cpp" "$ROOT/test/test_array.cpp" \
+    -o "$OUT/tests.js"
+  set +x
+  cp "$HERE/index.html" "$HERE/worker.js" "$OUT/"
+  echo "tests.wasm:  $(stat -f%z "$OUT/tests.wasm") bytes"
+fi
 
-cp "$HERE/index.html" "$HERE/worker.js" "$OUT/"
-echo "wasm size: $(stat -f%z "$OUT/tests.wasm") bytes"
+if [[ "$WHAT" == all || "$WHAT" == census ]]; then
+  set -x
+  emcc "${common[@]}" \
+    -sJSPI=1 -sJSPI_EXPORTS=run_census \
+    -sEXPORT_NAME=createCensus \
+    -sEXPORTED_FUNCTIONS=_run_census,_malloc,_free \
+    "$HERE/census_wasm.cpp" \
+    -o "$OUT/census.js"
+  set +x
+  cp "$HERE/census.html" "$HERE/census-worker.js" "$OUT/"
+  echo "census.wasm: $(stat -f%z "$OUT/census.wasm") bytes"
+fi
diff --git a/test/wasm/census-worker.js b/test/wasm/census-worker.js
new file mode 100644
index 0000000..32acb9b
--- /dev/null
+++ b/test/wasm/census-worker.js
@@ -0,0 +1,53 @@
+// Same device-acquisition shape as worker.js: JS does the async part and hands
+// wasm a ready device, so every C++ entry point stays synchronous.
+import createCensus from "./census.js";
+
+let log = "";
+
+async function main() {
+  if (!navigator.gpu) {
+    postMessage({ type: "fail", msg: "navigator.gpu missing (no WebGPU in this worker)" });
+    return;
+  }
+  const adapter = await navigator.gpu.requestAdapter();
+  if (!adapter) {
+    postMessage({ type: "fail", msg: "requestAdapter returned null" });
+    return;
+  }
+  // A census runs far longer than the test suite and its numbers are only
+  // meaningful if nothing was silently dropped, so surface validation errors
+  // the same way — one invalid dispatch discards its whole command buffer.
+  const device = await adapter.requestDevice();
+  const gpuErrors = [];
+  device.onuncapturederror = (e) => {
+    if (gpuErrors.length < 5) gpuErrors.push(String(e.error.message).split("\n")[0]);
+  };
+  device.lost.then((info) => gpuErrors.push("DEVICE LOST: " + info.message));
+  const info = adapter.info ?? {};
+
+  const mod = await createCensus({
+    preinitializedWebGPUDevice: device,
+    print: (s) => { log += s + "\n"; postMessage({ type: "line", line: s }); },
+    printErr: (s) => { log += s + "\n"; postMessage({ type: "line", line: s }); },
+  });
+
+  postMessage({
+    type: "ready",
+    adapter: `${info.vendor ?? "?"} / ${info.architecture ?? "?"} / ${info.description ?? ""}`,
+  });
+
+  onmessage = async () => {
+    log = "";
+    const t0 = performance.now();
+    let ok = 0;
+    try {
+      ok = await mod.ccall("run_census", "number", [], [], { async: true });
+    } catch (err) {
+      log += "exception: " + String((err && err.stack) || err) + "\n";
+    }
+    if (gpuErrors.length) log = "WEBGPU ERRORS:\n" + gpuErrors.join("\n") + "\n\n" + log;
+    postMessage({ type: "done", ok: !!ok, log, ms: performance.now() - t0 });
+  };
+}
+
+main().catch((e) => postMessage({ type: "fail", msg: String((e && e.stack) || e) }));
diff --git a/test/wasm/census.html b/test/wasm/census.html
new file mode 100644
index 0000000..2d519b3
--- /dev/null
+++ b/test/wasm/census.html
@@ -0,0 +1,47 @@
+
+
+cpp-tensorlib WebGPU census
+
+
cpp-tensorlib WebGPU census

+
booting…

+

+

+
diff --git a/test/wasm/census_wasm.cpp b/test/wasm/census_wasm.cpp
new file mode 100644
index 0000000..41c9d56
--- /dev/null
+++ b/test/wasm/census_wasm.cpp
@@ -0,0 +1,171 @@
+// Browser port of misc/census.cpp for the WebGPU backend (M10, Phase 4).
+//
+// Same purpose as the native census — find, per kernel class, the size where
+// GPU-total (single op + flush) first beats CPU — but it has to run in a page
+// because that is the only place this backend exists. It is a separate wasm
+// from the test harness so the two never interfere: the census needs long runs
+// and the suite needs to stay a fast pass/fail.
+//
+// Differences from misc/census.cpp, all forced by the environment:
+//
+//   - The wasm CPU path is scalar single-threaded (~8 GF/s), roughly 30x
+//     slower than the M1 Pro's Accelerate/AMX path, so the largest native
+//     sizes would take minutes. Trial counts fall off with size.
+//   - CPU and GPU alternate *which one runs first* every trial. Blocked A/B
+//     (all CPU, then all GPU) puts every one-time cost — first-touch
+//     allocation, GPU clock ramp — on whichever side ran first; culebra's
+//     Metal calibration was misled by exactly that (see its
+//     project_gpu_bench_shape_findings).
+//   - A whole-block section, absent from the native census, anchors
+//     batch_matmul_bias_threshold_(): that threshold is about a batch of ops
+//     riding one device together, which no single-op measurement can show.
+//
+// Output goes to stdout as plain text; the page just displays it.
+
+#include 
+#include 
+
+#include 
+#include 
+#include 
+#include 
+#include 
+#include 
+
+using tl::array;
+using clk = std::chrono::steady_clock;
+
+static array rnd(tl::shape_t s, unsigned seed) {
+  std::mt19937 g(seed);
+  std::uniform_real_distribution d(-1, 1);
+  size_t n = 1;
+  for (auto x : s) n *= (size_t)x;
+  std::vector v(n);
+  for (auto& x : v) x = d(g);
+  return array::from(std::move(v), std::move(s));
+}
+
+static double time1(const std::function& f) {
+  auto t0 = clk::now();
+  f();
+  return std::chrono::duration(clk::now() - t0).count();
+}
+
+static double median(std::vector v) {
+  std::sort(v.begin(), v.end());
+  return v[v.size() / 2];
+}
+
+struct ab {
+  double cpu, gpu;
+};
+
+// Interleaved A/B with the order rotated per trial, medians, warmup on both
+// devices first (the GPU one also compiles/caches the pipeline).
+static ab bench(int trials, const std::function& f) {
+  for (int i = 0; i < 3; i++) {
+    tl::use_cpu();
+    f();
+    tl::use_gpu();
+    f();
+  }
+  std::vector c, g;
+  for (int i = 0; i < trials; i++) {
+    if (i % 2 == 0) {
+      tl::use_cpu();
+      c.push_back(time1(f));
+      tl::use_gpu();
+      g.push_back(time1(f));
+    } else {
+      tl::use_gpu();
+      g.push_back(time1(f));
+      tl::use_cpu();
+      c.push_back(time1(f));
+    }
+  }
+  tl::use_cpu();
+  return {median(std::move(c)), median(std::move(g))};
+}
+
+static void row(const char* label, double n, ab t) {
+  std::printf("  %-22s n=%-12.0f cpu %9.3f ms   gpu %9.3f ms   %s %s\n", label, n,
+              t.cpu, t.gpu, t.gpu < t.cpu ? "**" : "  ", t.gpu < t.cpu ? "GPU" : "cpu");
+}
+
+extern "C" {
+
+// {async: true} from JS: the GPU waits beneath this suspend under JSPI.
+EMSCRIPTEN_KEEPALIVE int run_census() {
+  std::printf("[census] gpu_available=%d\n", tl::gpu_available() ? 1 : 0);
+  if (!tl::gpu_available()) {
+    std::printf("[census] FAIL: no GPU device\n");
+    return 0;
+  }
+
+  std::printf("== MATMUL (square, n=M*N*K) ==\n");
+  for (int64_t m : {32, 48, 64, 96, 128, 160, 192, 256, 384, 512}) {
+    auto a = rnd({m, m}, 1), b = rnd({m, m}, 2);
+    int trials = m <= 128 ? 15 : (m <= 256 ? 9 : 5);
+    auto t = bench(trials, [&] { a.dot(b).eval(); });
+    char buf[32];
+    std::snprintf(buf, sizeof buf, "%lldx%lldx%lld", (long long)m, (long long)m,
+                  (long long)m);
+    row(buf, (double)m * m * m, t);
+  }
+
+  std::printf("== ELEMENTWISE (add, n=elements) ==\n");
+  for (int64_t n : {(int64_t)1 << 12, (int64_t)1 << 14, (int64_t)1 << 16,
+                    (int64_t)1 << 18, (int64_t)1 << 20, (int64_t)1 << 22}) {
+    auto a = rnd({n}, 3), b = rnd({n}, 4);
+    auto t = bench(n <= (1 << 20) ? 15 : 7, [&] { (a + b).eval(); });
+    char buf[32];
+    std::snprintf(buf, sizeof buf, "%lld", (long long)n);
+    row(buf, (double)n, t);
+  }
+
+  std::printf("== REDUCTION (softmax rows, n=elements) ==\n");
+  for (int64_t cols : {256, 1024, 4096}) {
+    for (int64_t rows : {16, 64, 256, 1024}) {
+      auto a = rnd({rows, cols}, 5);
+      auto t = bench(rows * cols <= (1 << 18) ? 15 : 7,
+                     [&] { a.softmax().eval(); });
+      char buf[32];
+      std::snprintf(buf, sizeof buf, "%lldx%lld", (long long)rows, (long long)cols);
+      row(buf, (double)rows * cols, t);
+    }
+  }
+
+  // Whole-block anchor for batch_matmul_bias_threshold_(). One eval batch
+  // shaped like a transformer block: three seq*d*d projections, two
+  // seq*seq*d attention gemms, two seq*d*4d FFN gemms, plus the elementwise
+  // and softmax that sit between them. The question this answers is not "does
+  // the biggest gemm earn the GPU" but "does the whole batch earn it", which
+  // is what the batch bias decides.
+  std::printf("== BLOCK (transformer-ish, one eval batch; n=sum M*N*K) ==\n");
+  const int64_t seq = 128;
+  for (int64_t d : {32, 48, 64, 96, 128, 192, 256, 384, 512}) {
+    auto x = rnd({seq, d}, 6);
+    auto wq = rnd({d, d}, 7), wk = rnd({d, d}, 8), wv = rnd({d, d}, 9);
+    auto w1 = rnd({d, 4 * d}, 10), w2 = rnd({4 * d, d}, 11);
+    auto t = bench(d <= 256 ? 9 : 5, [&] {
+      auto q = x.dot(wq), k = x.dot(wk), v = x.dot(wv);
+      auto att = q.dot(k.transpose()).softmax();
+      auto h = att.dot(v) + x;
+      (h.dot(w1).relu().dot(w2) + h).eval();
+    });
+    double work = 3.0 * seq * d * d       // q, k, v projections
+                  + (double)seq * seq * d  // q @ k^T
+                  + (double)seq * d * seq  // att @ v
+                  + 2.0 * seq * d * 4 * d; // ffn
+    char buf[32];
+    std::snprintf(buf, sizeof buf, "seq%lld d%lld", (long long)seq, (long long)d);
+    row(buf, work, t);
+  }
+
+  for (auto& kv : tl::webgpu::context::get().dispatch_counts)
+    std::printf("[census] dispatch %s=%ld\n", kv.first.c_str(), kv.second);
+  std::printf("[census] DONE\n");
+  return 1;
+}
+
+}  // extern "C"
diff --git a/test/wasm/index.html b/test/wasm/index.html
index 1c855cd..62ddad8 100644
--- a/test/wasm/index.html
+++ b/test/wasm/index.html
@@ -11,13 +11,18 @@ 
cpp-tensorlib WebGPU tests