diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 96e95be..2651031 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,57 @@ 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 で有効化。 + # WebGPU (wasm) backend. Unlike cuda-build this does not stop at building — + # it runs the suite. + # + # It can, because what needs verifying is the WGSL and the C++, not Chrome: + # any runtime with navigator.gpu and JSPI will do. Deno has both, and on a + # GPU-less runner it still returns an adapter, through lavapipe (Mesa's + # software Vulkan). Headless Chrome under the same conditions does not expose + # navigator.gpu at all, across three attempts — including with a current + # Chrome and that same lavapipe working. So the obstacle was Chrome's own GPU + # gating, not the runner. + # + # This leaves WebGPU better covered than CUDA rather than worse. nvcc -ptx + # checks CUDA kernels at build time; WGSL is compiled at run time by the + # browser or runtime, so until now a shader error was invisible to the build. + # Running the suite compiles them. + # + # Both gpu and auto run against the ref oracle. auto is included because the + # TENSORLIB_WEBGPU arm of types.h's auto_threshold_() is unreachable from + # native ctest. Note the suite stays green even if every op falls back to + # CPU, so main_wasm.cpp asserts that each kernel family dispatched at least + # once and fails the run otherwise. + wasm: + name: wasm run (webgpu) / ubuntu-latest + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - uses: mymindstorm/setup-emsdk@v14 + with: + # Match the machine this was validated on: --use-port=emdawnwebgpu + # pulls a Dawn tied to the emsdk version, so moving this moves + # behaviour. + version: 6.0.3 + - uses: denoland/setup-deno@v2 + with: + deno-version: v2.x + - name: Install a software Vulkan driver (lavapipe) + run: | + sudo apt-get update -qq + sudo apt-get install -y -qq mesa-vulkan-drivers + # Build the census too, though nothing here runs it (it measures for + # minutes and is a local tool) — that it still compiles is the point. + - name: Build (tests + census) + run: ./test/wasm/build.sh + - name: Run the suite (gpu / auto) + run: deno run --allow-all test/wasm/deno_run.js + + # 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 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/README.md b/README.md index cf57cd2..b076a16 100644 --- a/README.md +++ b/README.md @@ -5,8 +5,8 @@ cpp-tensorlib A C++17 **header-only, cross-platform** tensor library with **zero third-party dependencies** — CPU (own SIMD kernels) and GPU (own Metal / -CUDA kernels) on macOS, Linux, and Windows, from a single `#include -`. No BLAS, no cuBLAS/cuDNN, no CUTLASS: every accelerated +CUDA / WGSL kernels) on macOS, Linux, Windows and the browser, from a single +`#include `. No BLAS, no cuBLAS/cuDNN, no CUTLASS: every accelerated kernel — GEMM, GEMV, attention, quantized (int4) matmul — is hand-written in this repo, so the only thing a consumer links against is the OS itself (and the CUDA *driver*, which is `dlopen`'d at runtime, not linked). @@ -19,8 +19,9 @@ all on the same own-kernel, zero-dependency foundation (see * **Header-only, C++17 minimum** (builds cleanly through C++23) — `#include `, nothing to link or build separately. * **Cross-platform**: macOS (Accelerate + Metal), Linux/Windows (own - BLIS-style CPU microkernels — scalar/AVX2/NEON — and own CUDA kernels). - One dispatch seam, no platform `#ifdef`s in consumer code. + BLIS-style CPU microkernels — scalar/AVX2/NEON — and own CUDA kernels), + and WebAssembly (own WGSL kernels over WebGPU). One dispatch seam, no + platform `#ifdef`s in consumer code. * **Zero third-party dependencies**: no OpenBLAS, cuBLAS, cuDNN, or CUTLASS anywhere in the accelerated paths (see [Dependency policy](#dependency-policy)). @@ -55,11 +56,139 @@ auto d = (a + b).relu(); // also lazy tl::eval(c, d); // evaluate both in one pass -tl::use_gpu(); // route to Metal (macOS) or CUDA (Linux/Windows) +tl::use_gpu(); // Metal (macOS), CUDA (Linux/Windows), WebGPU (wasm) auto e = a.dot(b); tl::use_auto(); // CPU or GPU per op, by measured size ``` +API reference +------------- + +`#include ` pulls in the whole library; everything lives in +namespace `tl`. The CMake target is `INTERFACE`, so it deliberately does not +impose a standard on you — set C++17 or newer yourself. + +`tokenizer.h` and `gguf.h` are *not* in the umbrella header; include them +directly if you want them. + +### Creating arrays + +Shapes are `tl::shape_t` (= `std::vector`), so braced literals work. +`{}` is a valid shape — it makes a rank-0 scalar. + +| | | +|---|---| +| `array::zeros(shape)` / `ones(shape)` | filled with 0 / 1 | +| `array::full(shape, v)` | filled with `v` | +| `array::empty(shape)` | allocated, uninitialized | +| `array::from(vector v)` | 1-d from host data (copies) | +| `array::from(vector v, shape)` | same, reshaped; throws if the count mismatches | +| `array()` | undefined array (`defined()` is false) | +| `tl::concat(vector)` | join along axis 0; other dims must match | + +### Inspecting + +`shape()`, `strides()`, `rank()`, `size()`, `contiguous()`, `defined()`, +`dt()`. Note it is **`rank()`, not `ndim()`**. + +Reading elements: `at({i, j})` (no bounds checking), `item()` for a +single-element array, `raw()` for a `const float*`, `data()` for a mutable +`float*`. All of these evaluate the graph first, and flush any pending GPU +work — mixing them with lazy ops is safe. + +**There is no printing.** No `operator<<`, no `to_string` — format arrays +yourself from `raw()`. + +For tests: `allclose(a, b, rtol = 1e-5f, atol = 1e-6f)` and `array_equal(a, b)` +(the same thing with zero tolerance). Both return false on a shape mismatch. + +### Operations + +Everything below is lazy unless noted. + +| Group | API | +|---|---| +| Arithmetic | `+ - * /` in array⊗array, array⊗float and float⊗array forms; `pow(a, b)`, `pow(a, s)` | +| Comparison | `> < >= <= == !=`, same three forms — result is an f32 mask of 1.0/0.0; `where(cond, a, b)` | +| Unary | `.exp() .log() .sqrt() .sigmoid() .relu()` | +| Softmax | `.softmax()` — **last axis only**, numerically stable, throws on rank 0 | +| Matmul | `a.dot(b)` — **rank 1 or 2 only**. 2d@2d→`{M,N}`, 2d@1d→`{M}`, 1d@2d→`{N}`, 1d@1d→scalar | +| Axis reductions | `.sum(axis, keepdims=false)`, `.mean(...)`, `.max(...)`, `.argmax(...)` — any axis, negatives count from the end; `argmax` returns indices as f32 | +| Scalar reductions | `float sum()`, `float max()`, `float mean()`, `int64_t argmax()` — **eager**, return host values | +| Views | `.transpose()`, `.transpose({axes})`, `.reshape(shape)`, `.slice(start, count)` (**axis 0 only**), `.clone()` | +| Broadcast VJP | `.sum_to(shape)` — reduce back to a shape that broadcasts to this one | +| In-place | `.add_(b)` — **eager**, mutates shared storage | +| Model ops | `array::attn_decode(q, K, V, scale)`, `array::rope(x, pos, base=10000.0f)`, `array::rmsnorm(x, weight, eps=1e-5f)`, `array::silu(x)`, `array::swiglu(gate, up)` | + +Broadcasting follows numpy's trailing-dimension rules; an incompatible pair +throws. `attn_decode` takes q `[H,D]`, K/V `[H,ctx,D]` and returns `[H,D]`. + +### Lazy evaluation + +Ops build a graph. Nothing runs until something forces it: + +```cpp +auto c = (a.dot(b) * 0.5f + 1.0f).relu(); // nothing computed yet +c.eval(); // now it runs +tl::eval(c, d, e); // or: several roots, one pass +``` + +`.eval()` (member) is the idiomatic form; the free `tl::eval(...)` is +variadic and evaluates multiple roots in a single topological pass, which is +what you want when their graphs share subexpressions. + +Evaluation also happens implicitly on `data()`, `raw()`, `item()`, `at()`, +the scalar reductions, `add_()`, and when a view is taken of a lazy source. + +### Choosing a backend + +```cpp +tl::use_cpu(); // default +tl::use_gpu(); // Metal / CUDA / WebGPU, per platform +tl::use_auto(); // per op, by measured size thresholds +bool ok = tl::gpu_available(); // compiled in AND a device is present +``` + +Selection is a process-wide switch, and which GPU backend you get is decided +at compile time. When no device is present, `use_gpu()` and `use_auto()` +silently run on the CPU. The `auto` thresholds are measured per kernel class +and per backend; `TL_BATCH_MATMUL_BIAS` overrides the batched-matmul one at +runtime. + +### Storage dtypes + +Compute and results are **always f32**. `bf16` and `q4` are weight-container +formats for the bandwidth-bound decode path, not general dtypes: + +```cpp +auto w = weights.to_bf16(); // f32 -> bf16 (round-to-nearest-even) +auto q = weights.to_q4(); // f32 [K,N] -> grouped int4; needs K % 32 == 0 +auto f = q.to_f32(); // back to f32 (a no-op if already f32) +``` + +Ops other than the native decode GEMV widen transparently. **Direct element +access requires f32** — `raw()`, `data()`, `at()` and `item()` throw on a +bf16/q4 array, so call `to_f32()` first. + +### Things that will bite you + +* **`data()` requires a contiguous array; `raw()` does not.** Calling `data()` + on a transposed view throws — `clone()` it first. +* **Views alias their base.** `transpose`/`reshape`/`slice` share storage, so a + write through `data()` is visible through every view of it — including from + unevaluated graphs still holding that array. +* **`add_()` mutates shared storage** and is not safe on an array that is still + feeding a lazy graph. +* **Evaluation is single-threaded** by design. There is no documented support + for sharing an array across threads. + +### Opt-in headers + +`#include ` gives `tl::tokenizer(gguf_path)` with `encode()`, +`decode()`, `bos_id()`, `eos_id()`. `#include ` gives +`tl::gguf::model(path)` — an mmap'd reader with `tensor(name)`, `tensors()`, +`kv(key)`, `metadata()`. Both are used by the LLM pipeline below. + Backends -------- @@ -67,6 +196,7 @@ Backends |----------|-----|-----| | macOS | Accelerate (vDSP / vForce / CBLAS) ✅ | Metal / MSL (`#embed` JIT), STEEL SGEMM ✅ | | Linux / Windows | own BLIS-style microkernels (scalar / AVX2 / NEON, runtime-dispatched) ✅ | own CUDA kernels, PTX JIT'd by a `dlopen`'d driver ✅ | +| WebAssembly | own scalar microkernels ✅ | own WGSL kernels over WebGPU ✅ | | any | reference strided implementation (oracle + fallback) ✅ | — | Everything reduces to strides through one index walker, so views, broadcast @@ -74,6 +204,31 @@ and transposed operands share a single code path. Accelerated backends replace it per-op at the `graph::eval_one` dispatch seam; the seam carries no platform `#ifdef`s (non-Apple builds get inline stubs). +All three GPU backends cover the array surface — GEMM, elementwise +unary/binary, rank-2 broadcast, row reductions (softmax / row sum / row max) +— each with the same fused affine epilogue, and each validated against the +reference implementation as an oracle. The LLM decode kernels (GEMV, +attention, RoPE, int4 dequant-matmul) are CUDA-only; on Metal and WebGPU +those ops decline at the seam and run on the CPU. + +**The WebGPU backend** is built with `emcc --use-port=emdawnwebgpu +-DTENSORLIB_WEBGPU`, and differs from the native two in two ways worth +knowing before you target it: + +* **Chrome in practice.** It needs JSPI to keep `flush()` and `sync_to_host()` + synchronous. Firefox has it behind a flag and Safari has not shipped it; + there, device acquisition fails, `available()` stays false, and every op + routes to CPU — correct, just not accelerated. That is also the fallback if + the page hands in no device. +* **No STEEL-equivalent tiling.** The auto-mode thresholds are calibrated + against measurements, but the WGSL SGEMM is a straightforward tiled kernel; + macOS gets a considerably more tuned one. + +CI runs the wasm suite on every push — `test/wasm/` builds the same oracle +suite as the native backends and runs it headless via Deno (a browser is not +required, locally either), asserting a per-kernel dispatch census so a silent +all-CPU fallback fails the build rather than passing green. + Dependency policy ------------------ @@ -90,6 +245,10 @@ platform links against: the CUDA driver is `dlopen`'d (`LoadLibrary` on Windows) at *run* time, so a binary built with `-DTENSORLIB_CUDA=ON` still runs — falling back to CPU — on a machine with no NVIDIA driver at all. +* **WebAssembly** — nothing linked. The WGSL kernels are committed as a C + string and the WebGPU entry points come from Emscripten's bundled + `emdawnwebgpu` port, so the toolchain is the only requirement and it is a + build-time one. * **Tests only** — `doctest` is vendored (not a runtime dependency of the library itself). @@ -118,7 +277,15 @@ Requires **C++17 or newer** — the headers use inline variables, compiles the tests at C++23, but consuming the headers only needs C++17. Exception: the **macOS/Metal** backend `#embed`s its shader source, so *building on macOS* additionally needs a `#embed`-capable compiler (Clang 19+); -non-Apple builds never reach that `#embed`. +non-Apple builds never reach that `#embed`. The WebGPU backend has the same +problem and solves it the other way — its WGSL is committed as a generated C +string (`kernels/tensorlib_webgpu_wgsl.inc`, regenerated by +`kernels/gen_wgsl_inc.sh` when the `.wgsl` changes), so it needs no `#embed` +and no build step. + +For a WebAssembly build, see `test/wasm/build.sh` — a flat `emcc` line with +`--use-port=emdawnwebgpu -DTENSORLIB_WEBGPU`, plus `-sJSPI=1`, which is what +lets `flush()` keep a synchronous signature. LLM inference ------------- 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. 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..874e8a5 --- /dev/null +++ b/include/gpu.h @@ -0,0 +1,38 @@ +#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__, 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 +// adding a browser GPU meant editing the CUDA header. + +#include "cuda.h" +#include "metal.h" +#include "webgpu.h" + +namespace tl { + +// 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; +#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 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/include/webgpu.h b/include/webgpu.h new file mode 100644 index 0000000..656df5e --- /dev/null +++ b/include/webgpu.h @@ -0,0 +1,702 @@ +#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 +#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. 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, 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 kUniformSlotBytes = 256; +// One flush can batch this many dispatches; past it, encode_ forces a blocking +// flush mid-graph. At 96 bytes of payload per 256-byte slot the ring is pure +// device memory (1 MB here) with no binding-size implication, so it is sized to +// put that forced stall well beyond any graph the backend is aimed at. +constexpr uint32_t kUniformSlotCount = 4096; + +// Every WGSL entry point. The context prebuilds a pipeline for each and the +// browser harness asserts each one dispatched — both need the same list, and a +// new kernel missing from either loses a guarantee silently. +inline constexpr const char* kEntryPoints[] = { + "sgemm", "ew_binary", "ew_unary", "ew_bcast", "softmax", "row_reduce"}; + +struct context { + wgpu::Instance instance; + wgpu::Device device; + wgpu::Queue queue; + wgpu::BindGroupLayout bgl; + wgpu::PipelineLayout play; + wgpu::ShaderModule mod; + wgpu::Buffer uniforms; // ring of kUniformSlots x kUniformSlot bytes + bool ready = false; + bool pending = false; + + // The open encoder and its compute pass. Created lazily on the first + // dispatch after a flush, so an idle backend submits nothing. One pass spans + // the whole batch, as metal.h keeps one MTLComputeCommandEncoder: WebGPU + // orders dispatches within a pass and makes each one's writes visible to the + // next, so chained ops stay correct, while a pass boundary per dispatch would + // serialize independent ones and give back what the batching is for. + wgpu::CommandEncoder enc; + wgpu::ComputePassEncoder pass; + 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; + + // Compute pipelines, keyed by WGSL entry point. Every one is built in the + // constructor and they all share a layout, so by the time encode_ runs this + // is a pure lookup and pipeline_'s create branch is unreachable. + 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; + } + + 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; + 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. + 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(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; + play = device.CreatePipelineLayout(&pld); + if (!play) return; + + wgpu::BufferDescriptor ud = {}; + ud.size = kUniformSlotBytes * kUniformSlotCount; + ud.usage = wgpu::BufferUsage::Uniform | wgpu::BufferUsage::CopyDst; + 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 : kEntryPoints) { + 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; + } + + 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; + } + + // 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 >= kUniformSlotCount) flush_(); // ring exhausted; new batch + + const uint32_t off = slot++ * (uint32_t)kUniformSlotBytes; + 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(); + pass = enc.BeginComputePass(); + } + pass.SetPipeline(pipe); + pass.SetBindGroup(0, bg, 1, &off); + pass.DispatchWorkgroups((uint32_t)gx, (uint32_t)gy, 1); + + 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()) { + 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; + c.pass.End(); + c.pass = nullptr; + 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) {})); +} + +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 +// 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); + m->where = context::BOTH; + } + stg.Unmap(); + } + // Only a completed memcpy makes the host copy current. Declaring BOTH on a + // failed readback would leave stale bytes permanently believed live, and no + // later sync_to_host would retry — so say so loudly instead, as the WGSL + // compile failure above does. + if (m->where == context::DEVICE) { + std::fprintf(stderr, "tensorlib webgpu: readback of %zu bytes failed\n", + m->bytes); + } + c.staging_pool[m->bytes].push_back(stg); + } + if (for_write) m->where = context::HOST; +} + +// 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 || 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; + + 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.ta = ta ? 1u : 0u; + p.tb = tb ? 1u : 0u; + p.scale = scale; + p.offset = offset; + return c.encode_("sgemm", ma, mb, mo, p, (n + 63) / 64, (m + 63) / 64); +} + +// 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 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); +} + +// 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); +} + +// 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, + 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..d72214a --- /dev/null +++ b/kernels/tensorlib_webgpu.wgsl @@ -0,0 +1,305 @@ +// 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.) + +// 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, // gemm rows | elementwise element count | reduce rows + N : u32, // gemm cols | reduce cols + K : u32, + lda : u32, + ldb : u32, + ldc : u32, + a_off : u32, + b_off : u32, + 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; +@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; + } + } +} + +// ---- 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; + +fn reduce_op(op : u32, acc : f32, v : f32) -> f32 { + if (op == OP_ROW_MAX) { return max(acc, v); } + return acc + v; +} + +// Tree-reduce one value per invocation through `scratch` and broadcast the +// result to the whole workgroup. Call only from uniform control flow — it +// barriers. The trailing barrier is what makes `scratch` safe to reuse for a +// second reduction in the same entry point (softmax's max phase, then its sum +// phase): without it, a fast invocation could overwrite scratch[0] before a +// slow one has read it. +fn tree_reduce(op : u32, lid : u32, v : f32) -> f32 { + scratch[lid] = v; + workgroupBarrier(); + for (var s : u32 = T / 2u; s > 0u; s = s >> 1u) { + if (lid < s) { scratch[lid] = reduce_op(op, scratch[lid], scratch[lid + s]); } + workgroupBarrier(); + } + let r = scratch[0]; + workgroupBarrier(); + return r; +} + +// 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]); } + let row_max = tree_reduce(OP_ROW_MAX, lid, m); + + var sum = 0.0; + for (var c : u32 = lid; c < p.N; c = c + T) { sum = sum + exp(A[src + c] - row_max); } + let inv = 1.0 / tree_reduce(OP_ROW_SUM, lid, sum); + + for (var c : u32 = lid; c < p.N; c = c + T) { + C[dst + c] = exp(A[src + c] - row_max) * inv; + } +} + +// 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]); + } + let r = tree_reduce(p.op, lid, acc); + if (lid == 0u) { C[p.c_off + row] = r * p.scale + p.offset; } +} diff --git a/kernels/tensorlib_webgpu_wgsl.inc b/kernels/tensorlib_webgpu_wgsl.inc new file mode 100644 index 0000000..47b3914 --- /dev/null +++ b/kernels/tensorlib_webgpu_wgsl.inc @@ -0,0 +1,308 @@ +// 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.) + +// 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, // gemm rows | elementwise element count | reduce rows + N : u32, // gemm cols | reduce cols + K : u32, + lda : u32, + ldb : u32, + ldc : u32, + a_off : u32, + b_off : u32, + 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; +@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; + } + } +} + +// ---- 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; + +fn reduce_op(op : u32, acc : f32, v : f32) -> f32 { + if (op == OP_ROW_MAX) { return max(acc, v); } + return acc + v; +} + +// Tree-reduce one value per invocation through `scratch` and broadcast the +// result to the whole workgroup. Call only from uniform control flow — it +// barriers. The trailing barrier is what makes `scratch` safe to reuse for a +// second reduction in the same entry point (softmax's max phase, then its sum +// phase): without it, a fast invocation could overwrite scratch[0] before a +// slow one has read it. +fn tree_reduce(op : u32, lid : u32, v : f32) -> f32 { + scratch[lid] = v; + workgroupBarrier(); + for (var s : u32 = T / 2u; s > 0u; s = s >> 1u) { + if (lid < s) { scratch[lid] = reduce_op(op, scratch[lid], scratch[lid + s]); } + workgroupBarrier(); + } + let r = scratch[0]; + workgroupBarrier(); + return r; +} + +// 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]); } + let row_max = tree_reduce(OP_ROW_MAX, lid, m); + + var sum = 0.0; + for (var c : u32 = lid; c < p.N; c = c + T) { sum = sum + exp(A[src + c] - row_max); } + let inv = 1.0 / tree_reduce(OP_ROW_SUM, lid, sum); + + for (var c : u32 = lid; c < p.N; c = c + T) { + C[dst + c] = exp(A[src + c] - row_max) * inv; + } +} + +// 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]); + } + let r = tree_reduce(p.op, lid, acc); + if (lid == 0u) { C[p.c_off + row] = r * p.scale + p.offset; } +} +)WGSL" 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) }));
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..519a179
--- /dev/null
+++ b/test/wasm/build.sh
@@ -0,0 +1,71 @@
+#!/usr/bin/env bash
+# Build the browser harnesses for the WebGPU backend.
+#
+#   ./build.sh && python3 -m http.server 8731
+#   # 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}"
+
+# Local dev sources emsdk from a checkout; CI (setup-emsdk) already has emcc on
+# PATH and has no emsdk_env.sh to source, so this must not be fatal there.
+EMSDK_ENV="${EMSDK_DIR:-$HOME/Projects/emsdk}/emsdk_env.sh"
+if [[ -f "$EMSDK_ENV" ]]; then
+  source "$EMSDK_ENV" >/dev/null 2>&1
+fi
+command -v emcc >/dev/null || { echo "emcc not found (set EMSDK_DIR or activate emsdk)" >&2; exit 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.
+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:  $(wc -c < "$OUT/tests.wasm") bytes"
+fi
+
+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: $(wc -c < "$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/deno_run.js b/test/wasm/deno_run.js
new file mode 100644
index 0000000..2fad2f4
--- /dev/null
+++ b/test/wasm/deno_run.js
@@ -0,0 +1,57 @@
+// Run the WebGPU test suite under Deno instead of a browser.
+//
+// The backend ships to browsers, but what CI needs to verify is the WGSL and
+// the C++ — not Chrome. Deno exposes navigator.gpu natively (wgpu) and its V8
+// has JSPI, so it can host this wasm with no browser, no display and no
+// puppeteer. That matters because headless Chrome on a hosted runner refuses
+// to expose WebGPU at all, even given a working software Vulkan driver; Deno
+// has none of Chrome's GPU sandboxing and blocklist machinery in the way.
+//
+// This mirrors test/wasm/worker.js: acquire the device here and hand it to
+// wasm fully formed, so the C++ side stays synchronous.
+//
+//   deno run --allow-all test/wasm/deno_run.js [path/to/tests.js]
+
+const modPath = Deno.args[0] ?? new URL("./site/tests.js", import.meta.url).href;
+
+if (!navigator.gpu) {
+  console.log("[deno] NO WEBGPU — navigator.gpu missing");
+  Deno.exit(2);
+}
+const adapter = await navigator.gpu.requestAdapter();
+if (!adapter) {
+  console.log("[deno] NO ADAPTER — requestAdapter returned null");
+  Deno.exit(2);
+}
+const device = await adapter.requestDevice();
+console.log("[deno] adapter:", JSON.stringify(adapter.info ?? {}));
+
+// One invalid dispatch discards the whole command buffer it batched into, which
+// presents as every op returning zeros. Surface that rather than debug blind.
+const gpuErrors = [];
+device.addEventListener?.("uncapturederror", (e) => gpuErrors.push(String(e.error?.message ?? e)));
+
+const createTests = (await import(modPath)).default;
+let out = "";
+const mod = await createTests({
+  preinitializedWebGPUDevice: device,
+  print: (s) => { out += s + "\n"; console.log(s); },
+  printErr: (s) => { out += s + "\n"; console.log(s); },
+});
+
+// gpu then auto, as the page does: gpu proves the kernels match the oracle,
+// auto proves the size-based routing in types.h does too.
+let ok = true;
+for (const [name, mode] of [["gpu", 1], ["auto", 2]]) {
+  out = "";
+  const t0 = performance.now();
+  const r = await mod.ccall("run_tests", "number", ["number"], [mode], { async: true });
+  console.log(`[deno] ${name} ${r ? "OK" : "FAILED"} in ${Math.round(performance.now() - t0)}ms`);
+  ok = ok && !!r;
+}
+if (gpuErrors.length) {
+  console.log("[deno] WEBGPU ERRORS:\n" + gpuErrors.join("\n"));
+  ok = false;
+}
+console.log("[deno] " + (ok ? "OK" : "FAILED"));
+Deno.exit(ok ? 0 : 1);
diff --git a/test/wasm/index.html b/test/wasm/index.html
new file mode 100644
index 0000000..62ddad8
--- /dev/null
+++ b/test/wasm/index.html
@@ -0,0 +1,60 @@
+
+
+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..a2adbe3
--- /dev/null
+++ b/test/wasm/main_wasm.cpp
@@ -0,0 +1,89 @@
+// 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.
+// mode: 0 = cpu, 1 = gpu, 2 = auto. auto is here because the TENSORLIB_WEBGPU
+// arm of auto_threshold_() (types.h) is unreachable from native ctest, so this
+// is the only place its size-based routing gets exercised against the oracle.
+EMSCRIPTEN_KEEPALIVE int run_tests(int mode) {
+  const bool gpu_mode = mode != 0;
+  const char* name = mode == 0 ? "cpu" : (mode == 1 ? "gpu" : "auto");
+  if (mode == 0) {
+    tl::use_cpu();
+  } else if (mode == 1) {
+    tl::use_gpu();
+  } else {
+    tl::use_auto();
+  }
+  std::printf("[tensorlib_wasm] mode=%s gpu_available=%d\n", name,
+              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;
+  }
+
+  // dispatch_counts is cumulative for the life of the module, so snapshot it:
+  // with more than one mode run per page load, the raw totals would credit
+  // this run with the previous one's dispatches.
+  auto before = tl::webgpu::context::get().dispatch_counts;
+
+  doctest::Context ctx;
+  int failed = ctx.run();
+
+  // A green suite does not prove the backend ran: an op the backend declines
+  // falls to CPU and its assertions pass either way. So the dispatch census is
+  // not just reported, it is asserted — every family that has a kernel must
+  // have dispatched at least once. Without this the whole suite stays green
+  // when the WGSL fails to compile, which is precisely the Phase 3 failure
+  // mode (a bad shader module makes every kernel decline, silently).
+  //
+  // gpu mode only. auto routes by size, so a family legitimately reaching zero
+  // there is a threshold decision, not a regression.
+  auto& counts = tl::webgpu::context::get().dispatch_counts;
+  for (const char* f : tl::webgpu::kEntryPoints) {
+    auto now = counts.find(f);
+    auto was = before.find(f);
+    long n = (now == counts.end() ? 0 : now->second) -
+             (was == before.end() ? 0 : was->second);
+    if (n) std::printf("[tensorlib_wasm] dispatch %s=%ld\n", f, n);
+    if (mode == 1 && n <= 0) {
+      std::printf("[tensorlib_wasm] FAIL: no %s dispatch — the backend "
+                  "declined every op in this family\n", f);
+      failed = 1;
+    }
+  }
+
+  std::printf("[tensorlib_wasm] %s %s\n", name, 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..577cfd1
--- /dev/null
+++ b/test/wasm/worker.js
@@ -0,0 +1,62 @@
+// 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();
+
+  // 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
+  // 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 ?? ""}`,
+  });
+
+  const MODES = { cpu: 0, gpu: 1, auto: 2 };
+  onmessage = async (e) => {
+    const gpuMode = MODES[e.data] ?? 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";
+    }
+    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 });
+  };
+}
+
+main().catch((e) => postMessage({ type: "fail", msg: String(e && e.stack || e) }));