Post-mortem: a GPU stack without nvcc, cuBLAS, or LLVM — measured boundaries, and three headline numbers we retracted

[dancinlife]

Post-mortem: a GPU stack without nvcc, cuBLAS, or LLVM — and three headline numbers we had to retract

Status: no longer actively developed (last commit 2026-07-19). MIT licensed. Every number below traces to a verdict file in this repo.

We built a compiler that emits PTX without LLVM, wrote our own GEMM kernels, and trained models with zero cuBLAS calls.

The headline result: a from-scratch GEMM can reach bit-exact parity with cuBLAS — at ~93% of its roofline, with no vendor call — in one dtype and one shape band. Outside that band it costs 1.5×, and the reason is settled rather than fixable. §1 is the whole boundary in one table.

We also published three wrong headline numbers along the way — one claiming we were far faster than PyTorch, one claiming we were catastrophically slower, one fabricated outright — and our own process caught and retracted all three (§6). If you are building something similar, that section may be worth more than the measurements.

---

Scope — read this first

This does not remove the need for NVIDIA hardware. PTX targets NVIDIA GPUs. If you are short on silicon, nothing here helps.

What it removes is dependence on NVIDIA's software stack — nvcc, cuBLAS, LLVM — and measures what that costs. If you are building a compiler or runtime and want to know whether owning the GEMM is feasible, the numbers below are the answer, including the parts where it is not.

Trust note. We published three wrong headline numbers over this project's life and retracted all three. They are catalogued in §6, including one where an agent-driven campaign fabricated pass-claims outright. Read §6 before quoting anything here — it is also the reason the numbers that survived are worth something.

---

1. The boundary, up front

Everything below is measured. This table is the whole answer if you read nothing else.

question	answer	where
Can you write a GEMM at vendor parity with no cuBLAS call?	Yes — 1.08× of cuBLAS-TF32 @d=2048, ~93% of roofline, bit-exact (rel_rms = 0)	§2
Does that hold at other shapes?	No. ~1.50× slower @d=4096, and the lever family is exhausted, not backlogged	§2.1
Does it hold at other dtypes?	No. FP16/BF16 are 11.5× off — the vendor roofline doubles and yours does not	§2.2
On a consumer card?	Own kernel edges cuBLAS @d=768 (0.95×) on an RTX 5070	§2
Full training step vs PyTorch?	FP64: we tie/win (0.90–1.10×). TF32: torch wins 3–8×	§3
Will kernel fusion save you?	No — occupancy is the wall, not launches. Closed-negative, repeatedly	§4
What actually cost us the most days?	Silent wrong-arch ships, stale artifacts, a contaminated toolchain, one op with no kernel	§8

The short version: owning the GEMM is feasible in a band and expensive outside it. What a vendor BLAS really sells is not one great kernel — it is shape adaptivity. Budget for a dispatcher.

---

2. The settled result — you can own the GEMM, in one dtype

This is the finding worth taking. A from-scratch, no-vendor-call TF32 GEMM:

shape	own TFLOP/s	cuBLAS-TF32	ratio	rel_rms	parity
D=2048 (H100)	~315	~342	1.08×	0.000e+00	YES (~93% of roofline)
D=4096 (H100)	~284	~427	~1.50×	0.000e+00	NO
D=768 (RTX 5070)	—	—	0.95×	~1.3e-5	own EDGES cuBLAS

F-GPU-ROUTEA-KEEPBAND-MEASURE, F-OP54-SUMMER-OWNGEMM-TF32

rel_rms = 0.000e+00 is an equality, not a tolerance — bit-exact against the reference at every configuration. This is a device-resident GEMM callable inline from a persistent megakernel, which cuBLAS (a host API) fundamentally cannot be.

2.1 Where it stops, and why — shape-rigidity

At D=4096 the own kernel falls to ~1.50× slower, and this is settled as unrecoverable, not a tuning backlog:

Our kernel is one fixed 128×128 tile at every D. cuBLAS is shape-adaptive — at D=4096 it scales up +24.6% (better single-pass tile + CTA swizzle) while our fixed tile scales down −9.9% (2× K-loop drain). Register spill, occupancy drop, and a ptxas ceiling were all statically excluded as causes.

Two GPU builds then exhausted the remaining lever:

• CTA-swizzle in isolation (MODE 9): regresses −1.6% (280.5 vs 285.1 TFLOP/s)
• A 2-CTA/SM-preserving 128×256 tile (MODE 10): regresses −7.1% (263.3 vs 283.5) — the sequential-halves schedule serializes the wgmma pipeline and doubles the K-drain

No bit-exact 256-N schedule on sm_90a is both 2 CTA/SM and non-serialized. The D=4096 gap is bit-exactness-bound. F-OP45GPU-OCCUPANCY-SWEEP, F-OP52-TF32-GAP-CLOSE, F-OP55-NEWTILE-D4096

Generalizable: what a vendor BLAS actually sells you is not one great kernel — it is shape adaptivity. A single excellent tile reaches parity in a band and loses outside it. Budget for a dispatcher, not a kernel.

2.2 The parity is dtype-scoped to TF32 — and this is the trap

We ported the same kernel to 16-bit operands. It did not transfer:

dtype	own	cuBLAS roofline	ratio	parity
TF32	315 @d=2048	342	1.08×	YES
FP16	71.6 @d=4096	827.2	11.55× off	NO
BF16	71.1 @d=4096	816.1	11.48× off	NO

F-FUSION-SM90-WGMMA-W14-FP16

The own kernel ran at the same absolute throughput in FP16 as in TF32 (~71–76 TFLOP/s — same decode/occupancy-bound design). But the cuBLAS FP16 roofline doubled (827 vs 431). The precision change moved the target 2× away without lifting our kernel, and the same-dtype ratio widened from 6.09× to 11.5×.

Presenting "own-GEMM ≈ cuBLAS parity" without the dtype qualifier would overstate the claim — the same honest-number failure class as the retired 1656× figure.

This also qualifies the "just switch precision" advice you will hear (including from our own earlier notes, where a BF16 pivot gave 9.67× over FP64 cuBLAS). Precision changes move your roofline too. If your kernel is occupancy-bound rather than arithmetic-bound, a lower precision raises the vendor's ceiling and leaves you further behind.

---

3. The full training step, matched-dtype and compiled

The authoritative comparison (F-BENCH-1, RTX 5070, same step DAG both sides, interpreter eliminated, byte-exact run-to-run at every cell):

 dtype  B | flame step/s | torch-eager (ratio) | torch-compile (ratio) | best torch÷flame
 -------+--------------+---------------------+-----------------------+-----------------
  FP64  1 |     604.84   |  514.31  (0.85×)    |  662.81  (1.10×)      |  1.10×
  FP64  2 |     357.88   |  299.73  (0.84×)    |  350.01  (0.98×)      |  0.98×  TIES
  FP64  4 |     196.70   |  163.32  (0.83×)    |  182.65  (0.93×)      |  0.93×  WINS
  FP64  8 |     103.33   |   85.28  (0.83×)    |   93.12  (0.90×)      |  0.90×  WINS
  TF32  1 |    1600.97   | 4856.04  (3.03×)    | 4589.95  (2.87×)      |  3.03×
  TF32  8 |     470.91   | 2705.13  (5.74×)    | 3708.76  (7.88×)      |  7.88×

Matched-dtype, the apparent 1656× collapses to single digits.

• FP64: we tie or win (0.90–1.10×), and post the highest FP64 throughput at B=8. Torch has no tensor-core FP64 path; a hand-fused deterministic FP64 step is genuinely competitive there.
• TF32: torch leads 3.03× → 7.88×. That is cuBLAS/inductor's tuned GEMM against our naive tiled CUDA-core GEMM, and a torch win at TF32 is expected. Its absolute TF32 throughput dwarfs ours.

One methodological caution about our own benchmark suite: a separate sweep (F-BENCH-7) showed the flame-calls-cuBLAS lane winning 7/8 cells against torch.compile, up to 11.6×. That lane is a hand-written CUDA harness (tool/bench/flame_bench_step.cu — "the CLMConvMoE step reduced to its load-bearing structure"), not the hexa-emitted framework path. It measures what the architecture achieves in C/CUDA. Do not read it as the language's number; F-BENCH-1 is the one we consider authoritative.

---

4. Fusion was not the lever — measured repeatedly

The intuition is that a training step is launch-bound, so fuse it and utilization rises. We tested every version. It is not.

Host-removal, idle H100, D=1536: async single-stream 10.3–12.5% util (byte-eq broken by a race), piecewise CUDA graph 13.17%, whole-step graph 13.19%. And it got worse with a bigger model: 11.94% at D=1536 → 10.39% at D=2560.

CLOSED-NEGATIVE. The binding constraint is not host launch overhead — it is occupancy. The step is a chain of many small, largely sequential kernels that under-fill the H100: median util 2%, peak 71–100% only during the brief GEMMs. The H100 is too big for this model. F-FUSION-OCCUPANCY-WALL

A cooperative megakernel moved utilization −0.08 pp. Attention fusion variants: single-CTA tensor-core fusion 9.4–15.5× slower (it uses 1 of 48 SMs); cross-layer fused block ~588× slower; the FP64 mega-kernel 1.8–4.4× slower. Batch-fill gives a real but capped self-speedup — ~3× and then it flattens (F-FUSION-BATCHFILL).

The same workload hit ~90% util on a smaller GPU (RTX 5070 at D=1536). Utilization is GPU-relative. Before optimizing, check whether your GPU is simply too large for your model.

For calibration: published fusion work (FlashFuser, arXiv 2512.12949) reports 1.24× end-to-end. That is the order of magnitude fusion buys.

---

5. Compiling clean is not being correct

We built a fused single-token decode kernel. It emitted valid PTX, passed ptxas, launched, ran, and produced finite plausible floats.

Against an f64 CPU reference: max relative error 1.61 against a 1e-2 tolerance. Reading the emitted PTX back:

• the online softmax over QK^T was omitted entirely; attention-weighted PV used inv_sqrt_hd as a "deterministic-finite stand-in" weight
• the output GEMV read attn_out from a shared-memory region never written with real attention output ("bogus reuse to mute warnings")
• RMSNorm accumulated via atom.shared.add into a shared slot not zeroed at kernel entry

The 17-kernel eager path it was benchmarked against was also wrong (rel ~1.9): placeholder unit-scale RMSNorm, raw ex2.approx softmax with no max-subtract. Both arms were numerically hollow and both compiled cleanly. It was also 40–42× slower. F-FUSION-AUTOREGRESSIVE-DECODE-TIMED

Gate on numerics against a CPU reference from your first kernel, or you will benchmark stubs.

---

6. How we fooled ourselves three times

Publish these next to your results, or you will fool yourself the same way.

6.1 "2.95× faster than PyTorch eager" — RETRACTED (unit mismatch)

It divided a 2,500-step PyTorch run wall (336.85 s) by one step of ours (114 s). The perf gate built on it was voided. stdlib/flame/PERF.md:327-331

6.2 "1,656× / 2,207× SLOWER than PyTorch" — RETIRED (three compounding errors)

Our own attempt to correct 1.1 overcorrected into a different wrong number. F-FUSION-VS-PYTORCH reported flame at 0.167 step/s vs torch eager 276.7 / torch.compile 368.5 on an H100. It was wrong for three reasons at once:

1. Unfair dtype — it compared flame FP64 against torch TF32. Torch was on a tensor-core path our FP64 could not use.
2. Interpreted, not compiled — the 5.98 s/step was the interpreted trainer's per-step host glue (per-window token copy, a ~28-call eager AdamW tail, CE/softmax-grad host glue) with the GPU at 0% between GEMMs. Not the step kernels.
3. 2-point linear extrapolation — of that ~6 s/step artifact, not a measured throughput curve.

And the follow-up probe refuted our own diagnosis of it: F-FUSION-INTERP-ELIM showed that AOT-compiling the glue was ~1.0×. The bytecode interpreter was never the wall either.

6.3 A fabricated verdict — RETRACTED

F-HEXA-GPU-ROOFLINE-MS1-1B-N-PRESERVE-LANDED: "this verdict's pass-claims were FABRICATED and are void." An agent-driven campaign produced pass-claims for measurements that never ran. The repo's own audit caught it and the verdict file now carries the retraction in place of the claim.

The discipline that came out of this (FLAME+FORGE-vs-PYTORCH+CUBLAS.md): always compare matched dtype, always compare compiled-to-compiled, never extrapolate from 2 points, and keep the retraction in the same file as the original claim.

---

7. What we were actually building — the axis that did work

Speed was never the differentiator, and pretending otherwise produced all three retractions. The measured differentiator is cross-machine bit-exact training:

The same fixed-seed step produces byte-identical weights, gradients, and loss across 6 environments / 4 architecture-libc combinations — arm64-macOS (Darwin libm), x86_64-linux (glibc, 2 hosts), arm64-linux (glibc, Pi 5), x86_64-musl (Alpine).

PyTorch cannot do this, for reasons that are structural rather than fixable:

hole	why
libm transcendentals	not correctly-rounded; glibc ≠ Darwin in the last ULP
reductions	tree/warp order-dependent; atomic-scatter races
FMA fusion	clang fuses a*b+c to one rounding on arm64, two on x86
cuBLAS accumulation	DMMA order is vendor-"unspecified", drifts across GPU generations

We closed each by construction: every exp/erf/ln/sqrt on the step path is a fixed-iteration + − × ÷ routine (dt_exp, dt_erf via A&S 7.1.26 branchless, Newton sqrt) — no libm on the step path. Whole-step run-to-run max|Δ| = 0 over 17 weights + m + v + loss, with a distinct-seed negative control at 0.344 confirming the zero is a real pass. Comptime-folded float constants are serialized as bit-exact C99 hex-float literals and locked by a CI gate so a formatter change cannot silently alter trained bits.

The honest boundary: the bit-exact identity is FP64 self-determinism. TF32 fast-mode is byte-eq run-to-run at TF32 and loss-tracks FP64 to ~1e-7, but TF32 weights are not bit-equal to FP64 weights. dt_erf sits 1.38e-7 from any one platform's libm by design — we trade "matches one platform" for "matches across all platforms." Cross-platform byte measurements were on CPU; a cross-GPU-architecture byte measurement is not part of the result.

If you need maximum TFLOP/s or the ecosystem, use PyTorch — it is 2–8× faster at TF32 and that is the right tool. If you need a training run that reproduces bit-for-bit on different hardware years later, that column is empty and we filled it.

---

8. Failure modes that cost us days — steal these

The silent wrong-arch ship. The released CUDA asset was built -gencode arch=compute_80,code=sm_80 only — SASS with no PTX. On any newer GPU (RTX 5070 = sm_120) every kernel launch failed cudaErrorNoKernelImageForDevice silently: cuda_available()==1, [OWN-GEMM-FIRED] banners printed, device buffers stayed zero-filled, results were garbage (argmax dev=0 vs host=46, max|Δ|=551). Rebuilt with the PTX gencode added: argmax 32==32==32, max|Δ|=2.5e-06. The release comment had claimed sm_80 PTX would JIT forward — the flag never embedded the PTX.
→ Ship -gencode arch=compute_$SM,code=sm_$SM -gencode arch=compute_$SM,code=compute_$SM and assert cuobjdump -lptx finds a PTX image in CI.

The stale artifact left in place. _hx_cuda_farr_silu_gate_gpu was called at runtime.c:9972 with its only extern at :10553 — after the use. Under clang C99 an implicit declaration is fatal, so the host compile failed → the host object was never built → the ar fold was skipped → the previous wrong-arch archive stayed and shipped. The device-side compile had succeeded, making it worse.
→ A build step that fails and leaves the previous artifact in place is worse than one that fails loudly. Make fold/cache steps delete their target first.

The contaminated toolchain that invalidated a measurement. A benchmark reported 75,331 ms/step. The real number was 882 ms. A stock/contaminated compiler was not builtin-aware and mis-lowered forge_dispatch_* into a slow value-dispatch CPU fallback. A "4750× slower than PyTorch" figure derived from this family was labeled phantom.
→ If a benchmark is anomalously slow, suspect your toolchain before your code. Build the compiler from current source on the measuring host.

One op with no GPU kernel at all. nn_conv1d_bwd was a host-scalar 4-nested loop — T·Cout·Cin·K ≈ 201M iterations per call, 4 calls per step, no GPU kernel. Replaced with an im2col + GEMM + col2im path.
→ Before optimizing kernels, verify every op in the step actually has one.

Cache keys that omit the toolchain. Ours did, so an artifact built with a different compiler served a "GREEN" result in ~5 ms without rebuilding. → Include CUDA version and target arch in any GPU artifact cache key.

A link guard that inverted. 37 forge_dispatch_* host dispatchers were wrapped in #ifndef HEXA_CUDA, assuming a GPU-side file would supply them under -DHEXA_CUDA. That file was never wired in, so every CUDA-variant link failed. Fixed with #if 1 + __attribute__((weak)) so a strong GPU override still wins, and locked by a build-free, GPU-free static gate (tool/forge_dispatch_cuda_link_gate.sh).

---

9. What we never finished

• gpu_launch host-side lowering and cubin embedding. PTX comes out correctly; wiring host code to launch those kernels end-to-end does not exist. Kernels were fired via dispatch scripts. Main blocker for anyone picking this up.
• sm_120 reaches only one of two entry points. hexa build --target=nvptx64-nvidia-cuda-sm120 works and runs the real codegen. The compiler binary's own dispatch (compiler/main.hexa:1066-1072) has arms for sm90/sm80 only and exits 2 — even though NVPTX_TARGET_SM120 is defined and handled downstream (nvptx_target.hexa:92,5250, PTX ISA 8.7 at :5289). Small fix; we never made it.
• The language substrate never reached the kernels' quality. Ordinary hexa scalar/array code runs 2.9–23× slower than gcc -O2 (geomean 8.6×; array-map worst at 23.2×) because every scalar op is a boxed tag-dispatch runtime call — no register allocation, no unboxed i64, no strength reduction. The hand-tuned GEMM path reaches BLIS/cuBLAS parity, but ordinary user code never goes down that path. An unboxing campaign later landed a 15–18% win, byte-eq verified across 3 targets. That gap is the honest reason the framework layer never matched the kernel layer.
• Docs drifted behind code, consistently under-claiming. Comments describing shipped code as "stub" were stale by months. We were misled by our own comments while writing this post-mortem. Read the code.

---

10. Where things are

FLAME+FORGE-vs-PYTORCH+CUBLAS.md      ★ the reconciliation SSOT — start here
.verdicts/                            ★ every verdict verbatim, incl. retractions
ARCHITECTURE.json  convergence.records  107 anti-recurrence records (§7 came from here)
compiler/codegen/nvptx_target.hexa    PTX backend (~8,200 lines, no LLVM)
compiler/codegen/{metal,rocm}_target.hexa   other backends (less exercised)
gpu/SPEC.md                           the @gpu subset
self/native/hxqwen14b_cuda.cu         own GEMM kernels (WMMA2 / tiled / split-K)
docs/forge-routea-shape-adaptive.md   the shape-adaptive selector design
self/forge/PARADIGM.md                the May paradigm campaign (superseded in parts)
stdlib/flame/PERF.md                  wall-time ledger incl. the first retraction
docs/flame-machine-independent-training.md   the cross-machine bit-exactness result

---

11. If you are picking this up

1. Compare matched-dtype, compiled-to-compiled, and never extrapolate from 2 points. All three of our retractions came from violating one of these. §1.
2. Owning a bit-exact GEMM at vendor parity is achievable — in one dtype, in one shape band. ~93% of the cuBLAS-TF32 roofline at D=2048, bit-exact, no vendor call. Budget for a shape dispatcher, because that is what the vendor library actually sells. §2.
3. Precision changes move your roofline too. Going to FP16 doubled cuBLAS's ceiling and left our occupancy-bound kernel where it was. §2.2.
4. Gate on numerics against a CPU reference from your first kernel. Clean codegen proves nothing. §5.
5. Fusion is probably not your lever; occupancy is. Check whether your GPU is too big for your model before optimizing. §4.
6. Verify every op in your step has a GPU kernel before tuning any of them. §7.
7. Decide early whether you are competing on speed. We were not, and the three retractions came from acting as if we were. The capability we actually had — byte-identical training across machines and architectures — was measurable, defensible, and unavailable elsewhere.

MIT licensed. Fork it, lift pieces, or just take the measurements.

Where we were wrong, the correction sits in the same file as the original claim. That convention is worth more than any single number in this document.