1use super::qwen38::Qwen38Cfg;
16
17pub const BTG: usize = 8;
20
21pub const QWEN38_BATCHED_KERNELS: &str = r#"
24// `BMAX` bounds the per-thread accumulator array in the staged GEMVs. The emitted program and the
25// equivalence tests both define it ahead of this text, but the kernel set must still compile when a
26// caller does not — a shipped kernel string that fails on a missing macro is a trap for the next
27// harness that includes it. The guard is inert whenever a real batch bound was supplied.
28// `BTG` (sequences per warp-tile) and `MAXSEQ` (KV-cache depth) were referenced here without
29// defaults long before the MMQ work — the self-containment test found them the moment it existed.
30// Same class, same fix.
31#ifndef BTG
32#define BTG 8
33#endif
34#ifndef MAXSEQ
35#define MAXSEQ 64
36#endif
37#ifndef BMAX
38#define BMAX 32
39#endif
40// The MMQ tile's shape constants. Same rule as BMAX and for the same reason: the emitted program and
41// the equivalence tests define these ahead of this text, but the kernel set must still compile when a
42// caller does not. Adding them WITHOUT guards broke nine tests an hour after BMAX was fixed — the
43// identical failure, reintroduced. `the_kernel_set_compiles_with_no_defines_supplied` is the standing
44// check that stops a third instance.
45#ifndef MMQ_TJ
46#define MMQ_TJ 32
47#endif
48#ifndef MMQ_RI
49#define MMQ_RI 2
50#endif
51#ifndef MMQ_RJ
52#define MMQ_RJ 2
53#endif
54#ifndef MMQ_KC
55#define MMQ_KC 256
56#endif
57#ifndef MMQ_ACC
58#define MMQ_ACC double
59#endif
60// The tensor-core tile's shape constants, same guard rule and same reason as the MMQ set above.
61// TC_TJ must be 8, 16, 32 or 64: the warp -> sub-tile map needs TC_TJ/8 to divide the 8 warps.
62#ifndef TC_TJ
63#define TC_TJ 32
64#endif
65// How many 256-thread tile blocks ptxas must keep resident per SM. MEASURED on an H100 NVL
66// (2026-09-20): the tile compiles to 102 registers, so 256 threads x 102 = 26,112 of the SM's 65,536
67// and only TWO blocks fit -- 25% occupancy. nsys puts the tile at 70.9% of the step while it uses
68// 9.9% of tensor-core peak and 24% of achievable bandwidth: it is neither math- nor bandwidth-bound,
69// it is OCCUPANCY-bound, which is why halving its arithmetic bought only 6%.
70// 1 lets ptxas take every register it wants — the behaviour this replaces.
71// The analogous cap was worth +18.9% on the 30B MoE trio, where the optimum was an INTERIOR point
72// (mb=5 best, mb=6 spilled), so this is a knob to sweep and never a constant to guess.
73#ifndef TC_DYNSMEM
74#define TC_DYNSMEM 0
75#endif
76// ── IN-KERNEL STALL ATTRIBUTION ──────────────────────────────────────────────────────────────
77// `ncu` is permission-blocked in RunPod containers, which is why every attribution in this campaign
78// came from ablation and roofline arithmetic. `clock64()` needs no privileges and answers the one
79// question those cannot: where does a warp actually WAIT? One lane per warp accumulates, so the
80// counters cost a handful of atomics per chunk and the shipped path compiles them away entirely.
81#ifndef TC_PROBE
82#define TC_PROBE 0
83#endif
84#ifndef TC_TILED_W
85#define TC_TILED_W 0
86#endif
87#if TC_PROBE
88__device__ unsigned long long bq_tc_clk[4];
89#define BQ_CLK(v) unsigned long long v = clock64();
90#define BQ_ACC(i,a,b) if((threadIdx.x&31)==0) atomicAdd(&bq_tc_clk[i],(unsigned long long)((b)-(a)));
91#else
92#define BQ_CLK(v)
93#define BQ_ACC(i,a,b)
94#endif
95// ★ THE FALLBACK MUST BE THE SHIPPED VALUE. It used to be 1, which is a residency the program never
96// asks for: the emitter always writes `tc_occ_def` (3 on Hopper, 2 below) ahead of this guard. Anything
97// that compiles these kernels WITHOUT the emitter -- the decode bandwidth bench, most obviously -- then
98// silently measured one block per SM. That is 167 GB/s against the shipped 199 at the decode shape on a
99// 3060: the instrument was slandering the tile by 19%, which is the same defect as the bench that once
100// ran SPLIT=1. A default that no shipped configuration uses is a trap, not a default.
101#ifndef TC_MINBLOCKS
102#if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900
103#define TC_MINBLOCKS 3
104#else
105#define TC_MINBLOCKS 2
106#endif
107#endif
108// How many K-chunks the tile keeps staged at once. 1 is the shipped, measured behaviour: a ONE-deep
109// prefetch held in REGISTERS (`pw`), which costs ~32 registers at B=8 and is why the tile compiles to
110// 181 there (one block per SM). A deeper pipeline has to live in shared memory instead, which is what
111// `cp.async` and Hopper's TMA both exist to fill -- and the measured need is large: the tile moves
112// 24.33 GB at 728 GB/s while `bq_gdn_step` moves 9.6 GB at 2.0 TB/s in the same step on the same chip.
113// This dimension is the structural prerequisite for either. At TC_STAGES=1 every stage index folds to
114// 0 and the emitted code is the one that produced every number on record.
115#ifndef TC_STAGES
116#define TC_STAGES 1
117#endif
118#ifndef TC_WM
119#define TC_WM 4
120#endif
121#ifndef TC_KC
122#define TC_KC 128
123#endif
124// The row-parallel GEMV's shape. MV_TJ sequences live in registers, MV_RB output rows share each
125// activation load, and MV_NW warps cooperate on one row-group's K reduction. Same guard rule as above.
126#ifndef MV_TJ
127#define MV_TJ 8
128#endif
129#ifndef MV_RB
130#define MV_RB 2
131#endif
132#ifndef MV_NW
133#define MV_NW 1
134#endif
135// The residual-norm-quantise kernels reduce over the whole model dimension, so they take ONE BLOCK
136// PER SEQUENCE — at two sequences that is two blocks of a 132-SM part, and the block width is the
137// only knob left that buys memory-level parallelism. It was a literal 256; measured on an A100 at
138// B=2, widening to 1024 took the step from 31.05 ms to 29.51 ms (+5%) for no other change. Sizing a
139// block by a constant when its grid cannot fill the machine is the bug; the constant is named here
140// so the next person sees the choice instead of the number.
141#ifndef RMSN_NT
142#define RMSN_NT 1024
143#endif
144// How many chunks the split residual norm cuts the model dimension into. Only read on the fast lane.
145#ifndef NORM_G
146#define NORM_G 20
147#endif
148
149// ══ THE BATCHED GEMV, STAGED ════════════════════════════════════════════════════════════════════
150// Corrected target. llama.cpp does NOT switch to tensor cores at these batch sizes — mmq.cu:327 reads
151// `!fp16_mma_hardware_available(cc) || ne11 < MMQ_DP4A_MAX_BATCH_SIZE` with that constant set to 64, so for
152// every batch below 64 it runs a DP4A kernel. Its advantage over ours was never the instruction; it was the
153// STAGING. Our register-tiled form re-reads each sequence's activation slab from L2 for every output row,
154// which at B=32 is ~16x more activation traffic than weight traffic and is what pinned us at 128 GB/s of an
155// available 3000.
156//
157// Here the activation tile for ALL sequences is staged once per 512-wide K block in shared memory, and every
158// warp in the block reads it from there. Weights still stream from global — they are touched once and never
159// reused, so staging them would only add a copy.
160//
161// Structure follows llama.cpp's MMQ (ggml/src/ggml-cuda/mmq.cuh, MIT, (c) 2023-2024 The ggml authors);
162// the code is ours. Crucially the ACCUMULATION ORDER is unchanged from `bq_gemv_h` — same ascending K
163// blocks, same per-block scale product — so this is bit-identical, not merely close.
164__global__ void bq_gemv_h_smem(float* __restrict__ Y, const signed char* __restrict__ wq,
165 const __half* __restrict__ ws, const signed char* __restrict__ Xq,
166 const float* __restrict__ Xs, int K, int N, int Bn, int YS){
167 int warps = blockDim.x >> 5;
168 int o = blockIdx.x * warps + (threadIdx.x >> 5);
169 int lane = threadIdx.x & 31, tid = threadIdx.x;
170 int KB = K / 32;
171 extern __shared__ char gsmem[];
172 signed char* Asm = (signed char*)gsmem; // Bn x 512 staged activations
173 float* Ass = (float*)(gsmem + Bn * 512); // Bn x 16 staged scales
174 float acc[BMAX];
175 #pragma unroll
176 for(int t=0;t<BMAX;t++) acc[t]=0.f;
177 const signed char* wrow = wq + (long)(o < N ? o : 0) * K;
178 const __half* wsr = ws + (long)(o < N ? o : 0) * KB;
179 for(int k0=0;k0<K;k0+=512){
180 for(int i=tid;i<Bn*512;i+=blockDim.x){ int b=i>>9, c=i&511;
181 Asm[i] = (k0+c<K) ? Xq[(long)b*K+k0+c] : (signed char)0; }
182 for(int i=tid;i<Bn*16;i+=blockDim.x){ int b=i>>4, c=i&15;
183 Ass[i] = (k0/32+c<KB) ? Xs[(long)b*KB+k0/32+c] : 0.f; }
184 __syncthreads();
185 if(o < N && k0+lane*16 < K){
186 int4 w = *(const int4*)(wrow + k0 + lane*16);
187 int blk = lane >> 1; // which 32-block of this 512 the lane covers
188 float wsb = __half2float(wsr[k0/32 + blk]);
189 #pragma unroll
190 for(int t=0;t<BMAX;t++){
191 if(t>=Bn) break;
192 int4 xx = *(const int4*)(Asm + t*512 + lane*16);
193 int d = __dp4a(w.x,xx.x,0)+__dp4a(w.y,xx.y,0)+__dp4a(w.z,xx.z,0)+__dp4a(w.w,xx.w,0);
194 acc[t] += (float)d * (wsb * Ass[t*16 + blk]);
195 }
196 }
197 __syncthreads();
198 }
199 #pragma unroll
200 for(int t=0;t<BMAX;t++){
201 if(t>=Bn) break;
202 float a=acc[t];
203 for(int r=16;r>0;r>>=1) a+=__shfl_down_sync(0xffffffffu,a,r);
204 if(lane==0 && o<N) Y[(long)t*YS+o]=a;
205 }
206}
207
208
209// The TWO-PASS activation form, staged. This is the lane the accuracy result rests on — `bq_gemv_hl`
210// is what runs when the emitter is asked for the golden/oracle-closest decode — and it was the one
211// lane the staging fix did not cover, so it kept paying the L2 re-read that pinned the batch at
212// 128 GB/s of an available 3000.
213//
214// Both activation words are staged, so the shared-memory footprint is 2*Bn*512 + Bn*64 bytes
215// (34 KB at Bn=32, inside the 48 KB default block limit). Weights still stream from global: each
216// weight byte is read once per output row and never reused, so staging them would only add a copy.
217//
218// ACCUMULATION ORDER IS UNCHANGED from `bq_gemv_hl`: the same ascending 512-wide K chunks, the same
219// `blk = k0/32 + lane/2` block-scale index, and the same
220// `((float)dh + (float)dl/127) * (wsb * xs)` grouping per chunk. That is what makes this a free
221// substitution rather than a numerically different kernel — the golden token stream cannot move.
222// Every K in this model is a multiple of 512 (5120, 17408, 2048, 6144, 1024, 248320), so the
223// zero-padded staging covers the whole reduction and `bq_gemv_hl`'s scalar tail loop is empty.
224__global__ void bq_gemv_hl_smem(float* __restrict__ Y, const signed char* __restrict__ wq,
225 const __half* __restrict__ ws, const signed char* __restrict__ Xh,
226 const signed char* __restrict__ Xl, const float* __restrict__ Xs,
227 int K, int N, int Bn, int YS){
228 int warps = blockDim.x >> 5;
229 int o = blockIdx.x * warps + (threadIdx.x >> 5);
230 int lane = threadIdx.x & 31, tid = threadIdx.x;
231 int KB = K / 32;
232 extern __shared__ char gsmem[];
233 signed char* Hsm = (signed char*)gsmem; // Bn x 512 staged hi words
234 signed char* Lsm = (signed char*)(gsmem + Bn * 512); // Bn x 512 staged lo words
235 float* Ass = (float*)(gsmem + 2 * Bn * 512); // Bn x 16 staged block scales
236 float acc[BMAX];
237 #pragma unroll
238 for(int t=0;t<BMAX;t++) acc[t]=0.f;
239 const signed char* wrow = wq + (long)(o < N ? o : 0) * K;
240 const __half* wsr = ws + (long)(o < N ? o : 0) * KB;
241 for(int k0=0;k0<K;k0+=512){
242 for(int i=tid;i<Bn*512;i+=blockDim.x){ int b=i>>9, c=i&511;
243 int src = k0+c;
244 Hsm[i] = (src<K) ? Xh[(long)b*K+src] : (signed char)0;
245 Lsm[i] = (src<K) ? Xl[(long)b*K+src] : (signed char)0; }
246 for(int i=tid;i<Bn*16;i+=blockDim.x){ int b=i>>4, c=i&15;
247 Ass[i] = (k0/32+c<KB) ? Xs[(long)b*KB+k0/32+c] : 0.f; }
248 __syncthreads();
249 if(o < N && k0+lane*16 < K){
250 int4 w = *(const int4*)(wrow + k0 + lane*16);
251 int blk = lane >> 1;
252 float wsb = __half2float(wsr[k0/32 + blk]);
253 #pragma unroll
254 for(int t=0;t<BMAX;t++){
255 if(t>=Bn) break;
256 int4 hh = *(const int4*)(Hsm + t*512 + lane*16);
257 int4 ll = *(const int4*)(Lsm + t*512 + lane*16);
258 int dh = __dp4a(w.x,hh.x,0)+__dp4a(w.y,hh.y,0)+__dp4a(w.z,hh.z,0)+__dp4a(w.w,hh.w,0);
259 int dl = __dp4a(w.x,ll.x,0)+__dp4a(w.y,ll.y,0)+__dp4a(w.z,ll.z,0)+__dp4a(w.w,ll.w,0);
260 acc[t] += ((float)dh + (float)dl*(1.0f/127.0f)) * (wsb * Ass[t*16 + blk]);
261 }
262 }
263 __syncthreads();
264 }
265 #pragma unroll
266 for(int t=0;t<BMAX;t++){
267 if(t>=Bn) break;
268 float a=acc[t];
269 for(int r=16;r>0;r>>=1) a+=__shfl_down_sync(0xffffffffu,a,r);
270 if(lane==0 && o<N) Y[(long)t*YS+o]=a;
271 }
272}
273
274
275// ══ THE MMQ-SHAPED TILE — the lowering that actually wins at batch ══════════════════════════════
276//
277// MEASURED, H100 NVL, isolated projection N=17408 K=5120, double accumulators, KC=256:
278// B=16 0.185 -> 0.194 ms (0.95x, LOSS)
279// B=32 0.347 -> 0.285 ms (1.22x)
280// B=64 0.673 -> 0.484 ms (1.39x)
281// The win grows with the batch, so this is a shape rule with a threshold, not a replacement.
282//
283// WHY the register-tiled form runs out of road, which is the whole reason this kernel exists. It
284// reuses each WEIGHT element across all B sequences but each ACTIVATION element across exactly one
285// output row, so at B=32 it issues 64 activation loads per weight load. Neither escape works inside a
286// warp-register tile: BTG must cover the batch or the weight matrix is re-streamed (4x at B=32), and
287// rows-per-warp needs RI*BTG accumulators — registers measured 40 -> 63 -> 96 -> 168 as RI went 1 to
288// 8. Shared memory breaks the coupling: the tile size stops being a register-count problem, and each
289// thread keeps only its own RI x RJ accumulators.
290//
291// Two properties are load-bearing and were both learned the expensive way:
292// * The thread -> output map must be a BIJECTION onto the tile. With (TI/RI)*(TJ/RJ) < NT the
293// surplus threads wrap and recompute outputs another thread already produced: silently correct,
294// 4x the work. The static_assert below is that bug, caught.
295// * The shared rows must be PADDED. Unpadded rows are KC bytes apart and shared memory repeats
296// every 128 B, so every thread of a warp — which differ by row — hits the same bank: a 32-way
297// conflict on every read of the compute loop. One int of padding took the kernel from 0.704 ms
298// to 0.450 ms on the development card. It is the difference between losing and winning.
299//
300// Accumulation is in double: the per-32-block int32 dot is exact, so the only rounding is the sum of
301// K/32 float terms, and fp64 is half-rate on a datacentre part (measured free here: 1.22x with double
302// against 1.25x with float) while being 10x more accurate than the float sum. On consumer cards fp64
303// is 1/32 rate and this choice would be wrong — which is why it is a knob, not a constant.
304//
305// Tile shape and K-step follow llama.cpp's MMQ (ggml/src/ggml-cuda/mmq.cuh, MIT, (c) 2023-2024 The
306// ggml authors): 8 warps, 256 K-elements staged per barrier round. A sweep of KC over 128/256/512
307// independently picked 256, which is their MMQ_ITER_K. The DESIGN PARAMETERS are theirs; the code is
308// ours. Attribution recorded in benchmarks/results/qwen38_27b.json.
309//
310// This re-associates the K reduction, so unlike the staged GEMVs it is NOT bit-identical to the
311// register-tiled form; it is held to a tolerance by `batched_equivalence_lock`, exactly as the
312// tensor-core lowering is, and the golden token stream is the acceptance test.
313__global__ void __launch_bounds__(256) bq_gemm_mmq(
314 float* __restrict__ Y, const signed char* __restrict__ wq, const __half* __restrict__ ws,
315 const signed char* __restrict__ Xh, const signed char* __restrict__ Xl,
316 const float* __restrict__ Xs, int K, int N, int B, int YS){
317 constexpr int NCG = MMQ_TJ / MMQ_RJ; // column groups fixed by the sequence tile
318 constexpr int NRG = 256 / NCG; // row groups take the rest of the block
319 constexpr int TI = NRG * MMQ_RI; // ...which DERIVES the output-row tile
320 constexpr int SW = MMQ_KC/4 + 1; // ints per shared row; +1 puts rows on consecutive banks
321 static_assert(MMQ_TJ % MMQ_RJ == 0 && NCG * NRG == 256,
322 "thread -> output map must be a bijection onto the tile");
323 __shared__ int sW[TI][SW];
324 __shared__ int sH[MMQ_TJ][SW];
325 __shared__ int sL[MMQ_TJ][SW];
326 __shared__ __half sWs[TI][MMQ_KC/32];
327 __shared__ float sXs[MMQ_TJ][MMQ_KC/32];
328
329 const int tid = threadIdx.x;
330 const int i0 = blockIdx.x * TI;
331 const int KB = K/32;
332 const int rg = tid % NRG, cg = tid / NRG;
333 const int ir0 = rg*MMQ_RI, jr0 = cg*MMQ_RJ;
334
335 MMQ_ACC acc[MMQ_RI][MMQ_RJ];
336 #pragma unroll
337 for(int a=0;a<MMQ_RI;a++)
338 #pragma unroll
339 for(int b=0;b<MMQ_RJ;b++) acc[a][b]=(MMQ_ACC)0;
340
341 for(int k0=0;k0<K;k0+=MMQ_KC){
342 for(int idx=tid; idx<TI*(MMQ_KC/16); idx+=256){
343 int r=idx/(MMQ_KC/16), c=idx%(MMQ_KC/16); int o=i0+r;
344 const int4 v = (o<N && k0+c*16<K) ? ((const int4*)(wq+(long)o*K+k0))[c] : make_int4(0,0,0,0);
345 sW[r][c*4+0]=v.x; sW[r][c*4+1]=v.y; sW[r][c*4+2]=v.z; sW[r][c*4+3]=v.w;
346 }
347 for(int idx=tid; idx<MMQ_TJ*(MMQ_KC/16); idx+=256){
348 int r=idx/(MMQ_KC/16), c=idx%(MMQ_KC/16);
349 const int4 h = (r<B && k0+c*16<K) ? ((const int4*)(Xh+(long)r*K+k0))[c] : make_int4(0,0,0,0);
350 const int4 l = (r<B && k0+c*16<K) ? ((const int4*)(Xl+(long)r*K+k0))[c] : make_int4(0,0,0,0);
351 sH[r][c*4+0]=h.x; sH[r][c*4+1]=h.y; sH[r][c*4+2]=h.z; sH[r][c*4+3]=h.w;
352 sL[r][c*4+0]=l.x; sL[r][c*4+1]=l.y; sL[r][c*4+2]=l.z; sL[r][c*4+3]=l.w;
353 }
354 for(int idx=tid; idx<TI*(MMQ_KC/32); idx+=256){
355 int r=idx/(MMQ_KC/32), c=idx%(MMQ_KC/32); int o=i0+r;
356 sWs[r][c] = (o<N && k0/32+c<KB) ? ws[(long)o*KB + k0/32 + c] : __float2half(0.f);
357 }
358 for(int idx=tid; idx<MMQ_TJ*(MMQ_KC/32); idx+=256){
359 int r=idx/(MMQ_KC/32), c=idx%(MMQ_KC/32);
360 sXs[r][c] = (r<B && k0/32+c<KB) ? Xs[(long)r*KB + k0/32 + c] : 0.f;
361 }
362 __syncthreads();
363 #pragma unroll
364 for(int blk=0; blk<MMQ_KC/32; ++blk){
365 #pragma unroll
366 for(int a=0;a<MMQ_RI;a++){
367 float wsb = __half2float(sWs[ir0+a][blk]);
368 #pragma unroll
369 for(int b=0;b<MMQ_RJ;b++){
370 int dh=0, dl=0;
371 #pragma unroll
372 for(int q=0;q<8;q++){
373 int w = sW[ir0+a][blk*8+q];
374 dh = __dp4a(w, sH[jr0+b][blk*8+q], dh);
375 dl = __dp4a(w, sL[jr0+b][blk*8+q], dl);
376 }
377 acc[a][b] += (MMQ_ACC)(((float)dh + (float)dl*(1.0f/127.0f)) * (wsb * sXs[jr0+b][blk]));
378 }
379 }
380 }
381 __syncthreads();
382 }
383 #pragma unroll
384 for(int a=0;a<MMQ_RI;a++){ int o=i0+ir0+a; if(o>=N) continue;
385 #pragma unroll
386 for(int b=0;b<MMQ_RJ;b++){ int j=jr0+b; if(j>=B) continue; Y[(long)j*YS+o]=(float)acc[a][b]; } }
387}
388
389
390// ══ THE SAME matmul, LOWERED AS A GEMM ══════════════════════════════════════════════════════════
391// At batch 1 a projection is a GEMV and one warp per output row is the right shape. At batch B it is a GEMM,
392// and that shape becomes wrong: the register-tiled form issues ~64 dp4a per 512-weight block and re-reads
393// every sequence's activation slab through L2, which measured 128 GB/s of an available 3000 at B=32 — compute
394// bound, tensor cores idle, and 3.76x behind llama.cpp, which switches to its MMQ kernels for exactly this.
395//
396// Nothing about the model changes. The `.lg` still says `matmul`; the COMPILER picks the lowering from the
397// shape it is given. That is the same decision llama.cpp makes between mmvq and mmq, made in the emitter
398// instead of by hand.
399//
400// Golden survives because the int8 dot product accumulates in int32 — exact, no rounding — and the per-32
401// block scales are applied in registers in the SAME float grouping the dp4a path uses
402// (`(float)c * (ws * xs)` per k-block, summed in k order). The arithmetic that decides a token is unchanged.
403// Tile shape and K-step follow llama.cpp's MMQ configuration (ggml/src/ggml-cuda/mmq.cuh,
404// mmq-config-ampere.cuh): MMQ_ITER_K = 256, 8 warps, shared-memory tiles padded to dodge bank conflicts.
405// llama.cpp is MIT-licensed, Copyright (c) 2023-2024 The ggml authors; the DESIGN PARAMETERS below are taken
406// from it, the code is ours. Attribution recorded in benchmarks/results/qwen38_27b.json.
407//
408// Why it matters, measured: our first version staged only 32 K-elements per __syncthreads round, so a
409// K=5120 projection paid 320 barriers per output tile where llama.cpp pays 40. Same tensor-core instruction,
410// 8x the synchronization. Staging 256 K per round is the single structural difference.
411#define TBM 16
412#define TBN 64
413#define KSTEP 256
414__global__ void bq_gemm_mma(float* __restrict__ Y, const signed char* __restrict__ Wq,
415 const __half* __restrict__ Ws, const signed char* __restrict__ Xq,
416 const float* __restrict__ Xs, int K, int N, int M, int YS){
417 int KB=K/32; int mblk=blockIdx.y*TBM, nblk=blockIdx.x*TBN;
418 int warp=threadIdx.x>>5, lane=threadIdx.x&31, tid=threadIdx.x; int gid=lane>>2, tq=lane&3;
419 // +4 bytes of row padding: 32-byte rows make every lane of a quarter-warp hit the same shared-memory bank.
420 const int XP = KSTEP + 4, WP = KSTEP + 4;
421 extern __shared__ signed char smem[];
422 signed char* Xsm = smem; // TBM rows x XP
423 signed char* Wsm = smem + TBM*XP; // TBN rows x WP
424 float d0=0.f,d1=0.f,d2=0.f,d3=0.f;
425 for(int k0=0;k0<K;k0+=KSTEP){
426 // Stage KSTEP of K at once — one barrier pair per 256 elements instead of per 32.
427 for(int i=tid;i<TBM*KSTEP;i+=256){ int r=i/KSTEP,c=i%KSTEP;
428 Xsm[r*XP+c] = (mblk+r<M && k0+c<K) ? Xq[(long)(mblk+r)*K+k0+c] : (signed char)0; }
429 for(int i=tid;i<TBN*KSTEP;i+=256){ int r=i/KSTEP,c=i%KSTEP;
430 Wsm[r*WP+c] = (nblk+r<N && k0+c<K) ? Wq[(long)(nblk+r)*K+k0+c] : (signed char)0; }
431 __syncthreads();
432 // Then run the mma chain over the staged block with no further barriers.
433 for(int ks=0; ks<KSTEP && k0+ks<K; ks+=32){
434 int kb=(k0+ks)>>5;
435 int a0=*(const int*)(Xsm+gid*XP+ks+tq*4);
436 int a1=*(const int*)(Xsm+(gid+8)*XP+ks+tq*4);
437 int a2=*(const int*)(Xsm+gid*XP+ks+16+tq*4);
438 int a3=*(const int*)(Xsm+(gid+8)*XP+ks+16+tq*4);
439 int b0=*(const int*)(Wsm+(warp*8+gid)*WP+ks+tq*4);
440 int b1=*(const int*)(Wsm+(warp*8+gid)*WP+ks+16+tq*4);
441 int c0=0,c1=0,c2=0,c3=0;
442 asm volatile("mma.sync.aligned.m16n8k32.row.col.s32.s8.s8.s32 {%0,%1,%2,%3}, {%4,%5,%6,%7}, {%8,%9}, {%0,%1,%2,%3};"
443 : "+r"(c0),"+r"(c1),"+r"(c2),"+r"(c3) : "r"(a0),"r"(a1),"r"(a2),"r"(a3), "r"(b0),"r"(b1));
444 int r0=mblk+gid, r1=mblk+gid+8, cbase=nblk+warp*8;
445 float xs0 = (r0<M) ? Xs[(long)r0*KB+kb] : 0.f;
446 float xs1 = (r1<M) ? Xs[(long)r1*KB+kb] : 0.f;
447 float ws0 = (cbase+tq*2 <N) ? __half2float(Ws[(long)(cbase+tq*2)*KB+kb]) : 0.f;
448 float ws1 = (cbase+tq*2+1 <N) ? __half2float(Ws[(long)(cbase+tq*2+1)*KB+kb]) : 0.f;
449 d0+=(float)c0*(ws0*xs0); d1+=(float)c1*(ws1*xs0);
450 d2+=(float)c2*(ws0*xs1); d3+=(float)c3*(ws1*xs1);
451 }
452 __syncthreads();
453 }
454 int r0=mblk+gid, r1=mblk+gid+8, cc0=nblk+warp*8+tq*2, cc1=cc0+1;
455 if(r0<M){ if(cc0<N) Y[(long)r0*YS+cc0]=d0; if(cc1<N) Y[(long)r0*YS+cc1]=d1; }
456 if(r1<M){ if(cc0<N) Y[(long)r1*YS+cc0]=d2; if(cc1<N) Y[(long)r1*YS+cc1]=d3; }
457}
458
459// ── THE KERNEL THAT MAKES BATCHING WORTH IT ─────────────────────────────────────────────────────
460// One warp owns one output row. It walks the row's weights ONCE, and for each 512-weight block it reuses the
461// loaded `int4` across BTG sequences held in registers. Weight traffic is therefore independent of the batch:
462// the same bytes that produced one token now produce BTG of them.
463__global__ void bq_gemv_h(float* Y, const signed char* __restrict__ wq, const __half* __restrict__ ws,
464 const signed char* __restrict__ Xq, const float* __restrict__ Xs,
465 int K, int N, int B, int YS){
466 int o=(blockIdx.x*blockDim.x+threadIdx.x)>>5; if(o>=N) return;
467 int lane=threadIdx.x&31; int b0=blockIdx.y*BTG;
468 const int4* w16=(const int4*)(wq+(long)o*K); const __half* wsr=ws+(long)o*(K/32);
469 int nq=K/512, KB=K/32;
470 float acc[BTG];
471 #pragma unroll
472 for(int t=0;t<BTG;t++) acc[t]=0.f;
473 for(int c=0;c<nq;c++){
474 int4 w=w16[c*32+lane]; int blk=c*16+lane/2; float wsb=__half2float(wsr[blk]);
475 #pragma unroll
476 for(int t=0;t<BTG;t++){
477 int bb=b0+t; if(bb>=B) break;
478 int4 xx=((const int4*)(Xq+(long)bb*K))[c*32+lane];
479 int d=__dp4a(w.x,xx.x,0)+__dp4a(w.y,xx.y,0)+__dp4a(w.z,xx.z,0)+__dp4a(w.w,xx.w,0);
480 acc[t] += (float)d * (wsb*Xs[(long)bb*KB+blk]);
481 }
482 }
483 int done=nq*512;
484 for(int t=0;t<BTG;t++){
485 int bb=b0+t; if(bb>=B) break;
486 for(int j=done+lane;j<K;j+=32){ int blk=j/32;
487 acc[t] += (float)((int)Xq[(long)bb*K+j]*(int)wq[(long)o*K+j])
488 * (__half2float(wsr[blk])*Xs[(long)bb*KB+blk]); }
489 }
490 for(int t=0;t<BTG;t++){
491 int bb=b0+t; if(bb>=B) break; float a=acc[t];
492 for(int s=16;s>0;s>>=1) a+=__shfl_down_sync(0xffffffffu,a,s);
493 if(lane==0) Y[(long)bb*YS+o]=a;
494 }
495}
496
497// The two-pass activation form, batched the same way: the weight is still read once per row for the whole
498// batch, and each sequence contributes two dp4a instead of one.
499__global__ void bq_gemv_hl(float* Y, const signed char* __restrict__ wq, const __half* __restrict__ ws,
500 const signed char* __restrict__ Xh, const signed char* __restrict__ Xl,
501 const float* __restrict__ Xs, int K, int N, int B, int YS){
502 int o=(blockIdx.x*blockDim.x+threadIdx.x)>>5; if(o>=N) return;
503 int lane=threadIdx.x&31; int b0=blockIdx.y*BTG;
504 const int4* w16=(const int4*)(wq+(long)o*K); const __half* wsr=ws+(long)o*(K/32);
505 int nq=K/512, KB=K/32;
506 float acc[BTG];
507 #pragma unroll
508 for(int t=0;t<BTG;t++) acc[t]=0.f;
509 for(int c=0;c<nq;c++){
510 int4 w=w16[c*32+lane]; int blk=c*16+lane/2; float wsb=__half2float(wsr[blk]);
511 #pragma unroll
512 for(int t=0;t<BTG;t++){
513 int bb=b0+t; if(bb>=B) break;
514 int4 hh=((const int4*)(Xh+(long)bb*K))[c*32+lane];
515 int4 ll=((const int4*)(Xl+(long)bb*K))[c*32+lane];
516 int dh=__dp4a(w.x,hh.x,0)+__dp4a(w.y,hh.y,0)+__dp4a(w.z,hh.z,0)+__dp4a(w.w,hh.w,0);
517 int dl=__dp4a(w.x,ll.x,0)+__dp4a(w.y,ll.y,0)+__dp4a(w.z,ll.z,0)+__dp4a(w.w,ll.w,0);
518 acc[t] += ((float)dh + (float)dl*(1.0f/127.0f)) * (wsb*Xs[(long)bb*KB+blk]);
519 }
520 }
521 int done=nq*512;
522 for(int t=0;t<BTG;t++){
523 int bb=b0+t; if(bb>=B) break;
524 for(int j=done+lane;j<K;j+=32){ int blk=j/32;
525 float a=(float)Xh[(long)bb*K+j] + (float)Xl[(long)bb*K+j]*(1.0f/127.0f);
526 acc[t] += a*(float)wq[(long)o*K+j]*(__half2float(wsr[blk])*Xs[(long)bb*KB+blk]); }
527 }
528 for(int t=0;t<BTG;t++){
529 int bb=b0+t; if(bb>=B) break; float a=acc[t];
530 for(int s=16;s>0;s>>=1) a+=__shfl_down_sync(0xffffffffu,a,s);
531 if(lane==0) Y[(long)bb*YS+o]=a;
532 }
533}
534
535// The tiny f32-weight projection GDN needs for its gates. Same shape, no quantization.
536// ── THE GDN's SMALL f32 PROJECTION, ROW-PARALLEL ────────────────────────────────────────────────
537// Third instance of the same shape mistake, and the clearest: one warp per output row, with only
538// N = 2*v_heads = 96 rows. That is 96 warps per sequence — 192 warps, ~6k threads, on a 132-SM part
539// — each grinding a K=5120 reduction 32 lanes at a time. Measured 0.45 ms over 48 launches.
540//
541// F32W_NW warps now share each row's reduction and rejoin through shared memory, the same fix as
542// `bq_gemv_wide` and for the same reason. Unlike that one it is NOT bit-identical: splitting a
543// reduction across warps re-associates it. The GDN gates it feeds are a sigmoid and a softplus decay,
544// both smooth and both far from a decision boundary, so a one-ulp move there cannot change a token —
545// but it is a different number, so it is asserted against the original within a bound rather than
546// claimed equal.
547#ifndef F32W_NW
548#define F32W_NW 8
549#endif
550__global__ void __launch_bounds__(256) bq_gemv_f32w(float* Y, const float* __restrict__ w,
551 const float* __restrict__ X, int K, int N, int B){
552 const int wpr=F32W_NW; // warps cooperating on one row
553 const int rows_per_block=(256/32)/wpr; // rows a block covers
554 const int lane=threadIdx.x&31, warp=threadIdx.x>>5;
555 const int o=blockIdx.x*rows_per_block + warp/wpr;
556 const int sub=warp%wpr;
557 const int bb=blockIdx.y; if(bb>=B) return;
558 __shared__ float red[256/32];
559 float a=0.f;
560 if(o<N){ const float* wr=w+(long)o*K; const float* x=X+(long)bb*K;
561 for(int j=lane+sub*32;j<K;j+=32*wpr) a+=wr[j]*x[j]; }
562 for(int s=16;s>0;s>>=1) a+=__shfl_down_sync(0xffffffffu,a,s);
563 if(lane==0) red[warp]=a;
564 __syncthreads();
565 if(o<N && sub==0 && lane==0){
566 float t=0.f;
567 for(int q=0;q<wpr;q++) t+=red[(warp/wpr)*wpr+q]; // fixed ascending order: deterministic
568 Y[(long)bb*N+o]=t;
569 }
570}
571
572// ── the per-sequence pointwise and reduction kernels: one block per sequence ─────────────────────
573__global__ void bq_embed(float* X, const signed char* wq, const __half* ws, const int* D_tok, int E){
574 int bb=blockIdx.y; long tok=D_tok[bb];
575 for(int t=blockIdx.x*blockDim.x+threadIdx.x; t<E; t+=gridDim.x*blockDim.x)
576 X[(long)bb*E+t]=(float)wq[tok*E+t]*__half2float(ws[tok*(E/32)+t/32]);
577}
578__global__ void bq_rmsnorm_q8f_hl(const float* X, signed char* Xq, signed char* Xl, float* Xs,
579 float* Xf, const float* w, int n, float eps){
580 int bb=blockIdx.x; const float* x=X+(long)bb*n; float* xf=Xf+(long)bb*n;
581 __shared__ float sh[RMSN_NT];
582 float s=0.f; for(int i=threadIdx.x;i<n;i+=blockDim.x) s+=x[i]*x[i];
583 sh[threadIdx.x]=s; __syncthreads();
584 for(int k=blockDim.x/2;k>0;k>>=1){ if(threadIdx.x<k) sh[threadIdx.x]+=sh[threadIdx.x+k]; __syncthreads(); }
585 float sc=rsqrtf(sh[0]/n+eps); __syncthreads();
586 for(int i=threadIdx.x;i<n;i+=blockDim.x) xf[i]=x[i]*sc*w[i];
587 __syncthreads();
588 int nb=n/32;
589 for(int b=threadIdx.x;b<nb;b+=blockDim.x){
590 float amax=1e-30f;
591 for(int j=0;j<32;++j){ float a=fabsf(xf[b*32+j]); if(a>amax) amax=a; }
592 float q=amax/127.0f; Xs[(long)bb*nb+b]=q;
593 for(int j=0;j<32;++j) lgx_q8hl(xf[b*32+j]/q, Xq+(long)bb*n, Xl+(long)bb*n, (long)b*32+j);
594 }
595}
596// The residual-add + RMSNorm + requantise that follows a projection. `SK` folds the split-K
597// reduction in: instead of reading a GEMM output some earlier kernel summed, it sums the S partials
598// itself, in the SAME order `bq_sk_reduce` used (0.f, then s ascending), so the float result is
599// bit-identical -- and the separate launch, plus a full write and read of Y, disappear.
600template<bool SK>
601__device__ __forceinline__ void bq_add_rmsnorm_q8f_hl_body(float* X, const float* Y, int S,
602 signed char* Xq, signed char* Xl,
603 float* Xs, float* Xf, const float* w, int n, float eps){
604 int bb=blockIdx.x; float* x=X+(long)bb*n; const float* y=Y+(long)bb*n; float* xf=Xf+(long)bb*n;
605 __shared__ float sh[RMSN_NT];
606 const long skst=(long)n*gridDim.x;
607 float s=0.f; for(int i=threadIdx.x;i<n;i+=blockDim.x){
608 float yv; if(SK){ yv=0.f; for(int q=0;q<S;q++) yv+=y[(long)q*skst+i]; } else yv=y[i];
609 float v=x[i]+yv; x[i]=v; s+=v*v; }
610 sh[threadIdx.x]=s; __syncthreads();
611 for(int k=blockDim.x/2;k>0;k>>=1){ if(threadIdx.x<k) sh[threadIdx.x]+=sh[threadIdx.x+k]; __syncthreads(); }
612 float sc=rsqrtf(sh[0]/n+eps); __syncthreads();
613 for(int i=threadIdx.x;i<n;i+=blockDim.x) xf[i]=x[i]*sc*w[i];
614 __syncthreads();
615 int nb=n/32;
616 for(int b=threadIdx.x;b<nb;b+=blockDim.x){
617 float amax=1e-30f;
618 for(int j=0;j<32;++j){ float a=fabsf(xf[b*32+j]); if(a>amax) amax=a; }
619 float q=amax/127.0f; Xs[(long)bb*nb+b]=q;
620 for(int j=0;j<32;++j) lgx_q8hl(xf[b*32+j]/q, Xq+(long)bb*n, Xl+(long)bb*n, (long)b*32+j);
621 }
622}
623__global__ void bq_add_rmsnorm_q8f_hl(float* X, const float* Y, signed char* Xq, signed char* Xl,
624 float* Xs, float* Xf, const float* w, int n, float eps){
625 bq_add_rmsnorm_q8f_hl_body<false>(X,Y,1,Xq,Xl,Xs,Xf,w,n,eps);
626}
627__global__ void bq_add_rmsnorm_q8f_hl_sk(float* X, const float* P, int S, signed char* Xq, signed char* Xl,
628 float* Xs, float* Xf, const float* w, int n, float eps){
629 bq_add_rmsnorm_q8f_hl_body<true>(X,P,S,Xq,Xl,Xs,Xf,w,n,eps);
630}
631__global__ void bq_add(float* X, const float* Y, int n){
632 int bb=blockIdx.y; int i=blockIdx.x*blockDim.x+threadIdx.x;
633 if(i<n) X[(long)bb*n+i]+=Y[(long)bb*n+i];
634}
635__global__ void bq_gate_mul(float* A, int as_, const float* G, int gs, int n){
636 int bb=blockIdx.y; int i=blockIdx.x*blockDim.x+threadIdx.x;
637 if(i<n){ float g=G[(long)bb*gs+i]; long o=(long)bb*as_+i; A[o]*= 1.0f/(1.0f+__expf(-g)); }
638}
639// The attention tail: gate, then requantise. Two launches over the SAME grid, the second reading back
640// what the first wrote to global memory one instruction earlier. Fused, the gated value never leaves a
641// register, and the amax reduction the quantiser already does covers it. Bit-identical: the expression
642// producing `v` is unchanged, and `A` is still written because the next step's oracle reads it.
643__global__ void bq_gate_mul_q8_hl(float* A, int as_, const float* G, int gs,
644 signed char* Xq, signed char* Xl, float* Xs, int K){
645 int bb=blockIdx.y;
646 int blk=blockIdx.x*(blockDim.x/32)+threadIdx.x/32; if((long)blk*32>=K) return; int lane=threadIdx.x&31;
647 const long i=(long)blk*32+lane;
648 const float g=G[(long)bb*gs+i];
649 const long o=(long)bb*as_+i;
650 const float v=A[o]*(1.0f/(1.0f+__expf(-g)));
651 A[o]=v;
652 float a=fabsf(v);
653 for(int s=16;s>0;s>>=1) a=fmaxf(a,__shfl_down_sync(0xffffffffu,a,s)); a=__shfl_sync(0xffffffffu,a,0);
654 float sc=a>0?a/127.0f:1.0f; if(lane==0) Xs[(long)bb*(K/32)+blk]=sc;
655 lgx_q8hl(v/sc, Xq+(long)bb*K, Xl+(long)bb*K, i);
656}
657__global__ void bq_quant_act_hl(const float* X, signed char* Xq, signed char* Xl, float* Xs, int K){
658 int bb=blockIdx.y;
659 int blk=blockIdx.x*(blockDim.x/32)+threadIdx.x/32; if((long)blk*32>=K) return; int lane=threadIdx.x&31;
660 const float* x=X+(long)bb*K;
661 float v=x[blk*32+lane], a=fabsf(v);
662 for(int s=16;s>0;s>>=1) a=fmaxf(a,__shfl_down_sync(0xffffffffu,a,s)); a=__shfl_sync(0xffffffffu,a,0);
663 float sc=a>0?a/127.0f:1.0f; if(lane==0) Xs[(long)bb*(K/32)+blk]=sc;
664 lgx_q8hl(v/sc, Xq+(long)bb*K, Xl+(long)bb*K, (long)blk*32+lane);
665}
666// ── THE SwiGLU EPILOGUE, PARALLELISED BY THE RIGHT UNIT ─────────────────────────────────────────
667// This measured 1.14 ms of a 20.3 ms step at B=2 — 64 launches at 17.8 us — and the reason was a
668// grid sized by the wrong quantity. The work unit is a 32-wide quantisation block, and the kernel
669// gave one THREAD to each: `Gs(dff)` launches dff/256 blocks of 256 threads, so 17408 threads start
670// and `if(blk>=nb) return` sends all but dff/32 = 544 of them home. 1.6% of the threads did the work,
671// each looping 32 elements through a `float t[32]` register array.
672//
673// One WARP per 32-block instead of one thread fixes it: every lane owns one element, the block's
674// amax comes from a shuffle reduction, and the same launch geometry now has every thread working.
675//
676// It is BIT-IDENTICAL, not a tolerance. `fmaxf` is associative and exact, so the shuffle reduction
677// produces the same `amax` the serial loop did; `sc` is therefore the same float, and each lane then
678// runs the identical quantisation on the identical value. Only which thread computes what has moved.
679template<bool SK>
680__device__ __forceinline__ void bq_silu_had_q8_hl_body(const float* F, int S, signed char* Hq, signed char* Hl, float* Hs, int n){
681 const int bb=blockIdx.y, nb=n/32;
682 const int warp=(blockIdx.x*blockDim.x+threadIdx.x)>>5, lane=threadIdx.x&31;
683 if(warp>=nb) return;
684 const float* g=F+(long)bb*2*n; const float* u=g+n;
685 const long idx=(long)warp*32+lane;
686 const long skst=(long)2*n*gridDim.y;
687 float gv,uv;
688 if(SK){ gv=0.f; uv=0.f; for(int q=0;q<S;q++){ gv+=g[(long)q*skst+idx]; uv+=u[(long)q*skst+idx]; } }
689 else { gv=g[idx]; uv=u[idx]; }
690 const float h=(gv/(1.0f+__expf(-gv)))*uv;
691 float a=fabsf(h);
692 for(int s=16;s;s>>=1) a=fmaxf(a,__shfl_xor_sync(0xffffffffu,a,s));
693 const float sc=fmaxf(a,1e-30f)/127.0f; // the serial form seeded amax at 1e-30
694 if(lane==0) Hs[(long)bb*nb+warp]=sc;
695 lgx_q8hl(h/sc, Hq+(long)bb*n, Hl+(long)bb*n, idx);
696}
697__global__ void bq_silu_had_q8_hl(const float* F, signed char* Hq, signed char* Hl, float* Hs, int n){
698 bq_silu_had_q8_hl_body<false>(F,1,Hq,Hl,Hs,n);
699}
700__global__ void bq_silu_had_q8_hl_sk(const float* P, int S, signed char* Hq, signed char* Hl, float* Hs, int n){
701 bq_silu_had_q8_hl_body<true>(P,S,Hq,Hl,Hs,n);
702}
703// Greedy pick over R rows of logits, WITHOUT advancing any position. Decode's `bq_argmax_adv` folds
704// the position bump into the pick because a step always consumes exactly one token; a speculative step
705// does not know how many it consumed until the drafts have been checked against these very rows, so
706// the two have to come apart. Same reduction tree and same tie-break (first index wins on >), so a row
707// scored here picks the token decode would have picked.
708// Lay the speculative chunk out the way a T>1 forward expects it: row t*B+b. Row 0 is the token the
709// model actually produced last; row 1 is the draft being checked. Nothing else in the verify pass
710// needs to know it is speculative -- it is an ordinary two-token chunk.
711// The draft head conditions on the embedding of the token just produced AND the hidden state of the
712// token BEFORE it. After a verify pass `x` holds both rows -- position p in row 0, p+1 in row 1 --
713// and which is "the one before" depends on whether the draft was accepted. Leaving it at row 0 makes
714// the draft condition on a state one token stale exactly when it was RIGHT. That diverges the stream
715// from greedy at the second token, and every throughput measurement was blind to it.
716__global__ void bq_mtp_pick_hidden(float* __restrict__ X, const int* __restrict__ Acc, int B, int E){
717 int b=blockIdx.y; if(!Acc[b]) return;
718 for(int i=blockIdx.x*blockDim.x+threadIdx.x; i<E; i+=gridDim.x*blockDim.x)
719 X[(long)b*E+i] = X[(long)(B+b)*E+i];
720}
721__global__ void bq_mtp_pack(int* __restrict__ D_tok, const int* __restrict__ D_draft, int B){
722 int b=blockIdx.x*blockDim.x+threadIdx.x; if(b>=B) return;
723 D_tok[B+b]=D_draft[b];
724}
725// Accept per sequence: the draft stands only if the model's own pick at that position agrees.
726// `V` holds the verify pass's greedy picks in the same row-major order, so V[b] is the TRUE token at
727// position p+1 and V[B+b] is the token at p+2 that is only real if the draft was.
728__global__ void bq_mtp_accept(const int* __restrict__ V, const int* __restrict__ D_draft,
729 int* __restrict__ D_tok, int* __restrict__ D_pos,
730 int* __restrict__ Acc, int B){
731 int b=blockIdx.x*blockDim.x+threadIdx.x; if(b>=B) return;
732 int ok = (V[b]==D_draft[b]) ? 1 : 0;
733 Acc[b] = ok;
734 if(ok) atomicAdd(Acc+B, 1); // Acc[B] is the running accept total
735 D_tok[b]= ok ? V[B+b] : V[b]; // next step consumes p+2 on accept, p+1 on reject
736 D_pos[b]= D_pos[b] + 1 + ok;
737}
738// Roll the recurrence back for the sequences whose draft was wrong. Accepted sequences keep the
739// state the verify pass left (it consumed a token that turned out to be real); rejected ones take
740// the checkpoint written after token 0, which the lock proves is exactly what a one-token chunk
741// leaves. Copying only the rejected slices is why this is cheap at a high accept rate.
742// ★ THE RECURRENT STATE'S STORAGE WIDTH. It is 151 MB per sequence — 9.66 GB of read-modify-write per
743// decode step at B=32, and at 4.8 GB resident it is the largest term in the step that cannot hide in
744// L2. Halving it halves that traffic AND doubles the concurrency ceiling, which is a serving metric of
745// its own: vLLM carries the identical state and its baseline here died on "max_num_seqs exceeds
746// available Mamba cache blocks".
747//
748// Ten mantissa bits, not eight. Measured against an f64 reference the output drifts ~1e-3 with f16 and
749// 1.7e-2 to 2.6e-1 with bf16 — and the drift does NOT grow with sequence length, because the decay
750// multiplies an error injected at step t by alpha^(T-t) before it reaches the output, so old error dies
751// instead of accumulating. That is what licenses a half-width state at all; it is asserted by
752// `the_gated_recurrence_tolerates_a_half_width_state_but_only_with_ten_mantissa_bits`.
753//
754// The body touches the state EXACTLY TWICE — one load into registers, one store at the end — so the
755// width is a typedef and nothing else in the recurrence changes. Compute stays f32 throughout.
756// ★ THE KV CACHE'S STORAGE WIDTH. The state lane alone does NOT double the concurrency ceiling, which
757// is a claim this ledger had to retract: the GDN state is one of three per-sequence terms and at 2K
758// context the KV cache is the LARGER of them, because it is kept as `float`. Per sequence at 2K:
759// 151 MB state + 268 MB KV + 8 MB conv. Halving the state alone is 427 -> 352 MB (1.21x more
760// sequences); halving both is 427 -> 218 MB (1.97x). Concurrency is a serving metric and vLLM carries
761// the identical recurrent state, so this is the half that makes the pair worth having.
762//
763// f16 is conservative here — attention is a weighted average over many keys, production engines ship
764// fp8 KV, and this keeps ten mantissa bits. Compute stays f32: only the cached representation narrows.
765#ifdef BQ_KV16
766typedef __half bq_kv_t;
767__device__ __forceinline__ float bq_kv_ld(__half v){ return __half2float(v); }
768__device__ __forceinline__ __half bq_kv_cv(float v){ return __float2half(v); }
769#else
770typedef float bq_kv_t;
771__device__ __forceinline__ float bq_kv_ld(float v){ return v; }
772__device__ __forceinline__ float bq_kv_cv(float v){ return v; }
773#endif
774#ifdef BQ_STATE16
775typedef __half bq_state_t;
776__device__ __forceinline__ float bq_st_ld(__half v){ return __half2float(v); }
777__device__ __forceinline__ __half bq_st_cv(float v){ return __float2half(v); }
778#else
779typedef float bq_state_t;
780__device__ __forceinline__ float bq_st_ld(float v){ return v; }
781__device__ __forceinline__ float bq_st_cv(float v){ return v; }
782#endif
783__global__ void bq_mtp_rollback(bq_state_t* __restrict__ S, const bq_state_t* __restrict__ SM,
784 const int* __restrict__ Acc, long per_seq){
785 int b=blockIdx.y; if(Acc[b]) return;
786 long i=(long)blockIdx.x*blockDim.x+threadIdx.x, stride=(long)gridDim.x*blockDim.x;
787 bq_state_t* d=S+(long)b*per_seq; const bq_state_t* src=SM+(long)b*per_seq;
788 for(; i<per_seq; i+=stride) d[i]=src[i];
789}
790// ── PRE-TILED WEIGHTS ────────────────────────────────────────────────────────────────────────────
791// The decode GEMM waits 67% of its time for data (measured in-kernel with clock64(), since ncu is
792// permission-blocked in RunPod containers) while sitting at 42% of a 3180 GB/s pure-read roofline.
793// Neither depth, occupancy, run length, tile shape nor the instruction moved it, because none of them
794// change the ACCESS PATTERN: row-major weights make a block read TI separate KC-byte streams, and a
795// reader that streams contiguously gets 95% of the roofline.
796//
797// Storing each (n-tile, k-chunk) contiguously turns that into ONE streaming read of TI*KC bytes. The
798// transform runs once at load: it is a pure permutation of the same bytes, so the arithmetic is
799// untouched and the equivalence lock must still see bit-identity.
800__global__ void bq_tile_w(signed char* __restrict__ dst, const signed char* __restrict__ src,
801 int N, int K, int TI, int KC){
802 long tiles_k = (long)((K + KC - 1) / KC);
803 long i = (long)blockIdx.x * blockDim.x + threadIdx.x, stride = (long)gridDim.x * blockDim.x;
804 long total = (long)((N + TI - 1) / TI) * tiles_k * TI * KC;
805 for(; i < total; i += stride){
806 long within = i % ((long)TI * KC);
807 long tile = i / ((long)TI * KC);
808 int r = (int)(within / KC), c = (int)(within % KC);
809 int nt = (int)(tile / tiles_k), kt = (int)(tile % tiles_k);
810 int n = nt * TI + r, k = kt * KC + c;
811 dst[i] = (n < N && k < K) ? src[(long)n * K + k] : (signed char)0;
812 }
813}
814// The same permutation for the per-32 scales, which are read on the same footing.
815__global__ void bq_tile_ws(__half* __restrict__ dst, const __half* __restrict__ src,
816 int N, int K, int TI, int KC){
817 int NB = K / 32, NBC = KC / 32;
818 long tiles_k = (long)((K + KC - 1) / KC);
819 long i = (long)blockIdx.x * blockDim.x + threadIdx.x, stride = (long)gridDim.x * blockDim.x;
820 long total = (long)((N + TI - 1) / TI) * tiles_k * TI * NBC;
821 for(; i < total; i += stride){
822 long within = i % ((long)TI * NBC);
823 long tile = i / ((long)TI * NBC);
824 int r = (int)(within / NBC), c = (int)(within % NBC);
825 int nt = (int)(tile / tiles_k), kt = (int)(tile % tiles_k);
826 int n = nt * TI + r, b = kt * NBC + c;
827 dst[i] = (n < N && b < NB) ? src[(long)n * NB + b] : __float2half(0.f);
828 }
829}
830// ── THE fp8 LANE ─────────────────────────────────────────────────────────────────────────────────
831// Q8_0 carries a scale every 32 elements, so the accumulator must be read out, converted and rescaled
832// FOUR times per mma issue. fp8 e4m3 accumulates in fp32 natively and takes ONE scale per output
833// channel at the very end -- which is what vLLM's fp8 path does, so this is the matched-precision
834// comparison rather than a cheaper one.
835//
836// The register shape is IDENTICAL to the int8 m16n8k32 the tile already uses (4 A regs, 2 B regs,
837// 4 accumulators), so the fragment map, the swizzle and the staging survive untouched: only the
838// instruction, the accumulator type and where the scale lands change. fp8 mma needs PTX ISA 8.1
839// (CUDA 12.1+); the box has 12.8 at /usr/local/cuda-12.8 while `nvcc` on PATH is 12.0.
840// fp8 needs PTX ISA 8.1 (CUDA 12.1+) to assemble and <cuda_fp8.h> to convert. The equivalence locks
841// compile with whatever `nvcc` is on PATH, which on this box is 12.0, so BOTH the header and the
842// instruction are guarded and an older toolkit still builds the whole kernel set -- it simply has no
843// fp8 lane. `__CUDA_ARCH__ >= 890` is the arch gate: Ada and Hopper have fp8 tensor cores, Ampere
844// does not, which is also what chip.rs has said since the start.
845#if (__CUDACC_VER_MAJOR__ > 12) || (__CUDACC_VER_MAJOR__ == 12 && __CUDACC_VER_MINOR__ >= 1)
846#include <cuda_fp8.h>
847#define BQ_FP8_TOOLKIT 1
848#else
849#define BQ_FP8_TOOLKIT 0
850#endif
851#if BQ_FP8_TOOLKIT && defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 890
852#define BQ_HAVE_FP8 1
853#else
854#define BQ_HAVE_FP8 0
855#endif
856__device__ __forceinline__ void bq_mma_e4m3(float& d0, float& d1, float& d2, float& d3,
857 unsigned a0, unsigned a1, unsigned a2, unsigned a3,
858 unsigned b0, unsigned b1){
859#if BQ_HAVE_FP8
860 asm volatile("mma.sync.aligned.m16n8k32.row.col.f32.e4m3.e4m3.f32 "
861 "{%0,%1,%2,%3}, {%4,%5,%6,%7}, {%8,%9}, {%0,%1,%2,%3};\n"
862 : "+f"(d0),"+f"(d1),"+f"(d2),"+f"(d3)
863 : "r"(a0),"r"(a1),"r"(a2),"r"(a3),"r"(b0),"r"(b1));
864#else
865 (void)a0;(void)a1;(void)a2;(void)a3;(void)b0;(void)b1;
866#endif
867}
868// Q8_0 -> e4m3 with ONE scale per output channel. The real weight is q*s_blk; the channel scale is
869// the largest magnitude in the row over 448 (e4m3's max finite), so the row uses the full range.
870// This is a REQUANTISATION and it loses precision on purpose -- int8 carries 8 mantissa bits where
871// e4m3 carries 4 -- which is exactly the trade vLLM already makes and which the accuracy column
872// against the f64 oracle has to price rather than assume.
873__global__ void bq_q8_to_e4m3_scale(float* __restrict__ chan_scale,
874 const signed char* __restrict__ q, const __half* __restrict__ s,
875 int N, int K){
876 int n = blockIdx.x; if(n >= N) return;
877 int KB = K/32;
878 __shared__ float red[256];
879 float m = 0.f;
880 for(int k = threadIdx.x; k < K; k += blockDim.x){
881 float w = (float)q[(long)n*K + k] * __half2float(s[(long)n*KB + (k>>5)]);
882 m = fmaxf(m, fabsf(w));
883 }
884 red[threadIdx.x] = m; __syncthreads();
885 for(int t = blockDim.x/2; t > 0; t >>= 1){
886 if(threadIdx.x < t) red[threadIdx.x] = fmaxf(red[threadIdx.x], red[threadIdx.x+t]);
887 __syncthreads(); }
888 if(threadIdx.x == 0) chan_scale[n] = (red[0] > 0.f) ? (red[0] / 448.0f) : 1.0f;
889}
890__global__ void bq_q8_to_e4m3(unsigned char* __restrict__ out, const float* __restrict__ chan_scale,
891 const signed char* __restrict__ q, const __half* __restrict__ s,
892 int N, int K){
893 int KB = K/32;
894 long i = (long)blockIdx.x*blockDim.x + threadIdx.x, stride = (long)gridDim.x*blockDim.x;
895 for(; i < (long)N*K; i += stride){
896 int n = (int)(i / K), k = (int)(i % K);
897 float w = (float)q[i] * __half2float(s[(long)n*KB + (k>>5)]);
898#if BQ_FP8_TOOLKIT
899 out[i] = (unsigned char)__nv_cvt_float_to_fp8(w / chan_scale[n], __NV_SATFINITE, __NV_E4M3);
900#else
901 out[i] = 0; // no fp8 toolkit: the lane is absent, not silently wrong
902#endif
903 }
904}
905__global__ void bq_argmax_rows(const float* L, int V, int* D_tok){
906 int bb=blockIdx.x; const float* l=L+(long)bb*V;
907 __shared__ float bv[256]; __shared__ int bi[256];
908 float best=-1e30f; int arg=0;
909 for(int i=threadIdx.x;i<V;i+=blockDim.x) if(l[i]>best){ best=l[i]; arg=i; }
910 bv[threadIdx.x]=best; bi[threadIdx.x]=arg; __syncthreads();
911 for(int k=blockDim.x/2;k>0;k>>=1){
912 if(threadIdx.x<k && bv[threadIdx.x+k]>bv[threadIdx.x]){ bv[threadIdx.x]=bv[threadIdx.x+k]; bi[threadIdx.x]=bi[threadIdx.x+k]; }
913 __syncthreads(); }
914 if(threadIdx.x==0) D_tok[bb]=bi[0];
915}
916__global__ void bq_argmax_adv(const float* L, int V, int* D_tok, int* D_pos){
917 int bb=blockIdx.x; const float* l=L+(long)bb*V;
918 __shared__ float bv[256]; __shared__ int bi[256];
919 float best=-1e30f; int arg=0;
920 for(int i=threadIdx.x;i<V;i+=blockDim.x) if(l[i]>best){ best=l[i]; arg=i; }
921 bv[threadIdx.x]=best; bi[threadIdx.x]=arg; __syncthreads();
922 for(int k=blockDim.x/2;k>0;k>>=1){
923 if(threadIdx.x<k && bv[threadIdx.x+k]>bv[threadIdx.x]){ bv[threadIdx.x]=bv[threadIdx.x+k]; bi[threadIdx.x]=bi[threadIdx.x+k]; }
924 __syncthreads(); }
925 if(threadIdx.x==0){ D_tok[bb]=bi[0]; D_pos[bb]++; }
926}
927
928// ── the hybrid's stateful halves, one sequence per grid.z ────────────────────────────────────────
929// Every sequence carries its own recurrence memory, conv ring and KV cache. That per-sequence state, not the
930// weights, is what bounds the batch: the weights are shared by everyone, the state is not.
931__global__ void bq_conv_silu_l2(const float* __restrict__ X, int xs, const float* __restrict__ wc,
932 float* __restrict__ CS, float* __restrict__ Y,
933 int C, int T, int nl2, float eps){
934 int bb=blockIdx.y;
935 const float* x=X+(long)bb*xs; float* cs=CS+(long)bb*T*C; float* y=Y+(long)bb*C;
936 int c=blockIdx.x*blockDim.x+threadIdx.x; if(c>=C) return;
937 for(int t=0;t<T-1;++t) cs[(long)t*C+c]=cs[(long)(t+1)*C+c];
938 cs[(long)(T-1)*C+c]=x[c];
939 float acc=0.f;
940 for(int t=0;t<T;++t) acc += cs[(long)t*C+c]*wc[(long)t*C+c];
941 float v=acc/(1.0f+__expf(-acc));
942 if(blockIdx.x>=nl2){ y[c]=v; return; }
943 __shared__ float ss[32];
944 int lane=threadIdx.x&31, warp=threadIdx.x>>5, nw=(blockDim.x+31)>>5;
945 float pp=gdn_warp_sum(v*v);
946 if(lane==0) ss[warp]=pp;
947 __syncthreads();
948 __shared__ float den;
949 if(threadIdx.x==0){ float t2=0.f; for(int j=0;j<nw;++j) t2+=ss[j]; den=sqrtf(t2+eps); }
950 __syncthreads();
951 y[c]=v/den;
952}
953__global__ void bq_gdn_gates(const float* __restrict__ G, const float* __restrict__ dt,
954 const float* __restrict__ alog, float* __restrict__ beta,
955 float* __restrict__ decay, int NV){
956 int bb=blockIdx.y; int i=blockIdx.x*blockDim.x+threadIdx.x; if(i>=NV) return;
957 const float* br=G+(long)bb*2*NV; const float* ar=br+NV;
958 beta[(long)bb*NV+i] = 1.0f/(1.0f+__expf(-br[i]));
959 decay[(long)bb*NV+i] = powf(1.0f+__expf(ar[i]+dt[i]), alog[i]);
960}
961__global__ void bq_gdn_step(const float* __restrict__ Q, const float* __restrict__ Kk, const float* __restrict__ V,
962 const float* __restrict__ BET, const float* __restrict__ ALP,
963 bq_state_t* __restrict__ S, float* __restrict__ O, int D, int NK, int NV, int CS){
964 const unsigned m=0xffffffffu;
965 int bb=blockIdx.z;
966 int h=blockIdx.x, lane=threadIdx.x&31, warp=threadIdx.x>>5;
967 int a=blockIdx.y*(blockDim.x>>5)+warp; if(a>=D) return;
968 int g=h%NK;
969 const float* kk=Kk+(long)bb*CS+(long)g*D;
970 const float* qq=Q +(long)bb*CS+(long)g*D;
971 bq_state_t* Sb=S+(long)bb*NV*D*D;
972 long base=(long)h*D*D+(long)a*D;
973 float sreg[8]; float u_=0.f,w_=0.f,kq_=0.f; int n=0;
974 for(int b=lane;b<D;b+=32){
975 float sv=Sb[base+b], kb=kk[b], qb=qq[b];
976 sreg[n++]=sv; u_+=sv*kb; w_+=sv*qb; kq_+=kb*qb;
977 }
978 for(int sh=16;sh>0;sh>>=1){ u_+=__shfl_down_sync(m,u_,sh); w_+=__shfl_down_sync(m,w_,sh); kq_+=__shfl_down_sync(m,kq_,sh); }
979 u_=__shfl_sync(m,u_,0); w_=__shfl_sync(m,w_,0); kq_=__shfl_sync(m,kq_,0);
980 float alpha=ALP[(long)bb*NV+h], beta=BET[(long)bb*NV+h];
981 float coef=beta*(V[(long)bb*CS+(long)h*D+a] - alpha*u_);
982 if(lane==0) O[(long)bb*NV*D+(long)h*D+a]=(alpha*w_+coef*kq_)*rsqrtf((float)D);
983 n=0; for(int b=lane;b<D;b+=32) Sb[base+b]=alpha*sreg[n++]+coef*kk[b];
984}
985#ifndef BQ_GDN_DELTA_TMAX
986#define BQ_GDN_DELTA_TMAX 64
987#endif
988// ── THE CHUNKED DELTA RULE — parallel across the chunk ───────────────────────────────────────────
989// `bq_gdn_chunk` is a SCAN: every token needs the state the previous token left, so T tokens serialise
990// with D^2 work apiece. Measured, that is prefill's dominant cost (0.1811 ms/layer/token => ~4.5 s over
991// 48 layers and a 512-token prompt), and no amount of state-residency fixes the dependency itself.
992//
993// The same recurrence closes in a unit-triangular solve. With A_t = prod(alpha) and alpha_t*A_{t-1}=A_t:
994//
995// c_t + SUM_{j<t} beta_t (A_t/A_j)(k_j.k_t) c_j = beta_t (v_t - A_t S0 k_t) (I+M) C = B
996// O_t = [ A_t (S0 q_t) + SUM_{j<=t} (A_t/A_j)(k_j.q_t) C_j ] / sqrt(D)
997// S_C = A_C S0 + SUM_j (A_C/A_j) C_j k_j^T
998//
999// Proved equal to the scan to machine epsilon in f64 by
1000// `the_chunked_delta_rule_equals_the_sequential_gated_recurrence` BEFORE this kernel was written, so a
1001// failure here is an implementation bug and never an algebra bug.
1002//
1003// The decisive property: every S0 term is INDEPENDENT across t, and in the C recursion the head
1004// dimension `a` is THREAD-LOCAL — no cross-lane reduction inside the serial part. The critical path
1005// falls from T*D^2 to T*D.
1006//
1007// v1 is written for correctness against the f64 oracle, not for speed: one block per (sequence, head),
1008// each thread owning whole rows of the state. Optimise only once it is proven equal.
1009// ★ TEMPLATED ON THE CHUNK LENGTH, and that is not a style choice. With `float proj[64]` sized by a
1010// runtime bound, ptxas reported **512 bytes stack frame** against the scan's 0 — the per-thread arrays
1011// were in LOCAL MEMORY, which is the identical trap `bq_gdn_chunk_body` was templated on `PER` to
1012// escape. Read the `-Xptxas -v` stack-frame line before believing any recurrence kernel is in registers.
1013template<int TT>
1014__device__ __forceinline__ void bq_gdn_delta_body(const float* __restrict__ Q, const float* __restrict__ Kk,
1015 const float* __restrict__ V, const float* __restrict__ BET,
1016 const float* __restrict__ ALP, bq_state_t* __restrict__ S,
1017 float* __restrict__ O, int D, int NK, int NV, int CS, int T, int B){
1018 extern __shared__ float bqd_sh[];
1019 const int h = blockIdx.x, bb = blockIdx.y, g = h % NK;
1020 const int tid = threadIdx.x, nthr = blockDim.x;
1021 float* sKK = bqd_sh; // [T][T] k_j . k_t
1022 float* sKQ = sKK + (long)T * T; // [T][T] k_j . q_t
1023 float* sA = sKQ + (long)T * T; // [T] A_t = prod alpha
1024 const float inv_sqrt_d = rsqrtf((float)D);
1025
1026 // A_t: T is small and this is the only serial-in-T step left in the kernel.
1027 if (tid == 0) {
1028 float run = 1.f;
1029 for (int t = 0; t < T; ++t) { run *= ALP[(long)(bb * T + t) * NV + h]; sA[t] = run; }
1030 }
1031 __syncthreads();
1032
1033 // K and Q staged once: the Gram matrices and every S0 contraction below read them repeatedly.
1034 float* sK = sA + T; // [T][D]
1035 float* sQ = sK + (long)T * D; // [T][D]
1036 for (int idx = tid; idx < T * D; idx += nthr) {
1037 const int t = idx / D, b = idx % D;
1038 sK[idx] = Kk[(long)(bb * T + t) * CS + (long)g * D + b];
1039 sQ[idx] = Q [(long)(bb * T + t) * CS + (long)g * D + b];
1040 }
1041 __syncthreads();
1042
1043 // The two Gram matrices, one thread per (t,j) pair.
1044 for (int idx = tid; idx < T * T; idx += nthr) {
1045 const int t = idx / T, j = idx % T;
1046 float akk = 0.f, akq = 0.f;
1047 for (int b = 0; b < D; ++b) {
1048 const float kj = sK[(long)j * D + b];
1049 akk += kj * sK[(long)t * D + b];
1050 akq += kj * sQ[(long)t * D + b];
1051 }
1052 sKK[idx] = akk; sKQ[idx] = akq;
1053 }
1054 __syncthreads();
1055
1056 // One thread owns a whole state row `a`, so the C recursion needs no cross-lane reduction.
1057 //
1058 // ★ The b-loop is OUTERMOST in every phase. v1 put it inside the t-loop and so re-read the S0 row
1059 // from global 2T+1 = 65 times per row; that alone made it 12x SLOWER than the scan it replaces
1060 // (2.1351 vs 0.1812 ms/layer/token). The arithmetic is identical either way — this is purely how
1061 // many times the same bytes are fetched.
1062 for (int a = tid; a < D; a += nthr) {
1063 const bq_state_t* s0row = S + (long)bb * NV * D * D + (long)h * D * D + (long)a * D;
1064 float proj[TT]; // reused: S0.k_t first, then S0.q_t
1065 float cvec[TT];
1066
1067 // Every loop that touches proj/cvec is bounded by TT and fully unrolled, so their indices are
1068 // compile-time constants and the arrays live in REGISTERS. With runtime bounds they were
1069 // addressable and ptxas put them in local memory (272 bytes of frame) no matter how the kernel
1070 // was templated -- templating alone is not enough, the INDICES have to be constant too.
1071 #pragma unroll
1072 for (int t = 0; t < TT; ++t) proj[t] = 0.f;
1073 for (int b = 0; b < D; ++b) {
1074 const float sv = bq_st_ld(s0row[b]);
1075 #pragma unroll
1076 for (int t = 0; t < TT; ++t) proj[t] += sv * sK[(long)t * D + b];
1077 }
1078 #pragma unroll
1079 for (int t = 0; t < TT; ++t) {
1080 const float beta = BET[(long)(bb * TT + t) * NV + h];
1081 const float va = V[(long)(bb * TT + t) * CS + (long)h * D + a];
1082 float acc = beta * (va - sA[t] * proj[t]);
1083 #pragma unroll
1084 for (int j = 0; j < TT; ++j)
1085 if (j < t) acc -= beta * (sA[t] / sA[j]) * sKK[t * TT + j] * cvec[j];
1086 cvec[t] = acc;
1087 }
1088
1089 #pragma unroll
1090 for (int t = 0; t < TT; ++t) proj[t] = 0.f;
1091 for (int b = 0; b < D; ++b) {
1092 const float sv = bq_st_ld(s0row[b]);
1093 #pragma unroll
1094 for (int t = 0; t < TT; ++t) proj[t] += sv * sQ[(long)t * D + b];
1095 }
1096 #pragma unroll
1097 for (int t = 0; t < TT; ++t) {
1098 float acc = sA[t] * proj[t];
1099 #pragma unroll
1100 for (int j = 0; j < TT; ++j)
1101 if (j <= t) acc += (sA[t] / sA[j]) * sKQ[t * TT + j] * cvec[j];
1102 O[(long)(bb * TT + t) * NV * D + (long)h * D + a] = acc * inv_sqrt_d;
1103 }
1104
1105 bq_state_t* srow = S + (long)bb * NV * D * D + (long)h * D * D + (long)a * D;
1106 for (int b = 0; b < D; ++b) {
1107 float acc = sA[TT - 1] * bq_st_ld(srow[b]);
1108 #pragma unroll
1109 for (int j = 0; j < TT; ++j) acc += (sA[TT - 1] / sA[j]) * cvec[j] * sK[(long)j * D + b];
1110 srow[b] = bq_st_cv(acc);
1111 }
1112 }
1113}
1114__global__ void bq_gdn_delta(const float* __restrict__ Q, const float* __restrict__ Kk,
1115 const float* __restrict__ V, const float* __restrict__ BET,
1116 const float* __restrict__ ALP, bq_state_t* __restrict__ S,
1117 float* __restrict__ O, int D, int NK, int NV, int CS, int T, int B){
1118 if (T == 8) bq_gdn_delta_body< 8>(Q,Kk,V,BET,ALP,S,O,D,NK,NV,CS,T,B);
1119 else if (T == 16) bq_gdn_delta_body<16>(Q,Kk,V,BET,ALP,S,O,D,NK,NV,CS,T,B);
1120 else if (T == 32) bq_gdn_delta_body<32>(Q,Kk,V,BET,ALP,S,O,D,NK,NV,CS,T,B);
1121 else if (T == 4) bq_gdn_delta_body< 4>(Q,Kk,V,BET,ALP,S,O,D,NK,NV,CS,T,B);
1122 else if (T == 1) bq_gdn_delta_body< 1>(Q,Kk,V,BET,ALP,S,O,D,NK,NV,CS,T,B);
1123#ifdef BQ_GDN_DELTA_BIG
1124 else if (T == 64) bq_gdn_delta_body<64>(Q,Kk,V,BET,ALP,S,O,D,NK,NV,CS,T,B);
1125#endif
1126 // No silent fallback. ptxas allocates registers for the WHOLE kernel, so a TT=64 instantiation that
1127 // nothing launches still charged the T=32 path 168 registers and a 512-byte frame. An unsupported T
1128 // must fail loudly instead: the emitter picks T, and `the_delta_rule_only_sees_chunk_lengths_it_was_
1129 // instantiated_for` holds it to this list.
1130 else __trap(); // no printf: these kernels are compiled into harnesses without <cstdio>, and an
1131 // unsupported chunk length must abort loudly rather than return a plausible tensor
1132}
1133
1134// The chunked GDN recurrence: `bq_gdn_step` with its middle wrapped in a T-loop so the per-(sequence,
1135// head) DxD state stays in REGISTERS across the chunk instead of round-tripping HBM once per token.
1136// Arithmetic is untouched -- same strided lane distribution, same __shfl_down tree, same accumulation
1137// order -- so one bq_gdn_chunk(T) is BIT-IDENTICAL to T bq_gdn_steps.
1138//
1139// PER = D/32 is a TEMPLATE parameter, not a loop bound: with a runtime `b < D` ptxas cannot prove the
1140// trip count, `sreg[]` lands in local memory (bq_gdn_step carries a 32-byte stack frame for exactly
1141// this reason), and the state that was supposed to be resident is back in HBM -- a kernel that looks
1142// restructured and isn't. PER = 0 is the generic fallback for a head dim that is not a multiple of 32.
1143//
1144// GDNC_ROWS state rows per warp share ONE k/q load. Once the state stops moving, k and q are the
1145// traffic that is left: every warp reads the same D floats of each, so reuse inside the warp divides
1146// that by GDNC_ROWS. Staging them in shared memory instead MEASURED WORSE (0.3665 vs 0.2394 ms at
1147// 1024 threads on a 3060) -- the per-token __syncthreads costs more than the L1 hits it saves.
1148#ifndef GDNC_ROWS
1149#define GDNC_ROWS 4
1150#endif
1151#define BQ_GDN_CHUNK_GY(D) (((D)+(4*GDNC_ROWS)-1)/(4*GDNC_ROWS))
1152template<int PER, bool SHK>
1153__device__ __forceinline__ void bq_gdn_chunk_body(
1154 const float* __restrict__ Q, const float* __restrict__ Kk, const float* __restrict__ V,
1155 const float* __restrict__ BET, const float* __restrict__ ALP,
1156 bq_state_t* __restrict__ S, float* __restrict__ O, int D, int NK, int NV, int CS,
1157 int T, int B, float* skq){
1158 const unsigned m=0xffffffffu;
1159 constexpr int R=GDNC_ROWS;
1160 constexpr int SS=PER?PER:8;
1161 int bb=blockIdx.z;
1162 int h=blockIdx.x, lane=threadIdx.x&31, warp=threadIdx.x>>5;
1163 int abase=(blockIdx.y*(blockDim.x>>5)+warp)*R;
1164 int g=h%NK;
1165 bq_state_t* Sb=S+(long)bb*NV*D*D;
1166 int aa[R]; bool lv[R];
1167 #pragma unroll
1168 for(int r=0;r<R;++r){ int a=abase+r; lv[r]=(a<D); aa[r]=lv[r]?a:0; }
1169 const int cnt=PER?PER:((D-lane+31)>>5);
1170 float sreg[SS*R];
1171 #pragma unroll
1172 for(int r=0;r<R;++r){
1173 long base=(long)h*D*D+(long)aa[r]*D;
1174 #pragma unroll
1175 for(int j=0;j<SS;++j) if(j<cnt) sreg[r*SS+j]=bq_st_ld(Sb[base+lane+j*32]);
1176 }
1177 for(int t=0;t<T;++t){
1178 long row=(long)t*B+bb;
1179 const float* kk=Kk+row*CS+(long)g*D;
1180 const float* qq=Q +row*CS+(long)g*D;
1181 if(SHK){
1182 __syncthreads();
1183 for(int i=threadIdx.x;i<2*D;i+=blockDim.x) skq[i]=(i<D?kk[i]:qq[i-D]);
1184 __syncthreads();
1185 kk=skq; qq=skq+D;
1186 }
1187 // The inner loop is bq_gdn_step's, verbatim, with only the row index added: u_, w_ and kq_
1188 // accumulate in ONE fused loop and reduce in ONE shuffle tree. Hoisting kq_ into its own loop
1189 // to compute it once per warp instead of once per row let ptxas contract the FMAs differently
1190 // and drifted the recurrence by an ULP at token 1 -- caught by the bit-identity lock. Recomputing
1191 // it per row is cheap ALU; what actually had to be shared is the k/q LOAD, and CSE across the
1192 // unrolled r-loop gives that for free.
1193 #pragma unroll
1194 for(int r=0;r<R;++r){
1195 float u_=0.f,w_=0.f,kq_=0.f;
1196 #pragma unroll
1197 for(int j=0;j<SS;++j) if(j<cnt){ int b=lane+j*32; float sv=sreg[r*SS+j],kb=kk[b],qb=qq[b]; u_+=sv*kb; w_+=sv*qb; kq_+=kb*qb; }
1198 for(int sh=16;sh>0;sh>>=1){ u_+=__shfl_down_sync(m,u_,sh); w_+=__shfl_down_sync(m,w_,sh); kq_+=__shfl_down_sync(m,kq_,sh); }
1199 u_=__shfl_sync(m,u_,0); w_=__shfl_sync(m,w_,0); kq_=__shfl_sync(m,kq_,0);
1200 float alpha=ALP[row*NV+h], beta=BET[row*NV+h];
1201 float coef=beta*(V[row*CS+(long)h*D+aa[r]] - alpha*u_);
1202 if(lane==0&&lv[r]) O[row*NV*D+(long)h*D+aa[r]]=(alpha*w_+coef*kq_)*rsqrtf((float)D);
1203 #pragma unroll
1204 for(int j=0;j<SS;++j) if(j<cnt) sreg[r*SS+j]=alpha*sreg[r*SS+j]+coef*kk[lane+j*32];
1205 }
1206 }
1207 #pragma unroll
1208 for(int r=0;r<R;++r){
1209 if(!lv[r]) continue;
1210 long base=(long)h*D*D+(long)aa[r]*D;
1211 #pragma unroll
1212 for(int j=0;j<SS;++j) if(j<cnt) Sb[base+lane+j*32]=bq_st_cv(sreg[r*SS+j]);
1213 }
1214}
1215#define BQ_GDN_CHUNK_DISPATCH(SHK, SH) \
1216 if((D&31)==0 && (D>>5)==4) bq_gdn_chunk_body<4,SHK>(Q,Kk,V,BET,ALP,S,O,D,NK,NV,CS,T,B,SH); \
1217 else if((D&31)==0 && (D>>5)==2) bq_gdn_chunk_body<2,SHK>(Q,Kk,V,BET,ALP,S,O,D,NK,NV,CS,T,B,SH); \
1218 else if((D&31)==0 && (D>>5)==1) bq_gdn_chunk_body<1,SHK>(Q,Kk,V,BET,ALP,S,O,D,NK,NV,CS,T,B,SH); \
1219 else if((D&31)==0 && (D>>5)==8) bq_gdn_chunk_body<8,SHK>(Q,Kk,V,BET,ALP,S,O,D,NK,NV,CS,T,B,SH); \
1220 else bq_gdn_chunk_body<0,SHK>(Q,Kk,V,BET,ALP,S,O,D,NK,NV,CS,T,B,SH);
1221__global__ void bq_gdn_chunk(const float* __restrict__ Q, const float* __restrict__ Kk, const float* __restrict__ V,
1222 const float* __restrict__ BET, const float* __restrict__ ALP,
1223 bq_state_t* __restrict__ S, float* __restrict__ O, int D, int NK, int NV, int CS,
1224 int T, int B){
1225 BQ_GDN_CHUNK_DISPATCH(false, nullptr)
1226}
1227// The same recurrence, additionally checkpointing the state after token `MID` into `SM`.
1228//
1229// Speculative decode needs this and nothing else can supply it: a verify pass advances the state
1230// through k+1 tokens, but only the ACCEPTED prefix is real, and the recurrence is destructive. The
1231// state cannot be recovered by replay (it is 4.8 GB at 32 sequences) and it cannot be snapshotted
1232// before the pass without paying that 4.8 GB up front for a rollback that usually will not happen.
1233// Writing it on the way past costs one extra state write and makes the rollback a pointer choice.
1234// `MID` counts tokens CONSUMED, so MID=1 is the state after the first token has been applied.
1235template<int PER, bool SHK>
1236__device__ __forceinline__ void bq_gdn_chunk_mid_body(
1237 const float* __restrict__ Q, const float* __restrict__ Kk, const float* __restrict__ V,
1238 const float* __restrict__ BET, const float* __restrict__ ALP,
1239 bq_state_t* __restrict__ S, bq_state_t* __restrict__ SM, float* __restrict__ O,
1240 int D, int NK, int NV, int CS, int T, int B, int MID, float* skq){
1241 const unsigned m=0xffffffffu;
1242 constexpr int R=GDNC_ROWS;
1243 constexpr int SS=PER?PER:8;
1244 int bb=blockIdx.z;
1245 int h=blockIdx.x, lane=threadIdx.x&31, warp=threadIdx.x>>5;
1246 int abase=(blockIdx.y*(blockDim.x>>5)+warp)*R;
1247 int g=h%NK;
1248 bq_state_t* Sb=S+(long)bb*NV*D*D;
1249 bq_state_t* Sm=SM+(long)bb*NV*D*D;
1250 int aa[R]; bool lv[R];
1251 #pragma unroll
1252 for(int r=0;r<R;++r){ int a=abase+r; lv[r]=(a<D); aa[r]=lv[r]?a:0; }
1253 const int cnt=PER?PER:((D-lane+31)>>5);
1254 float sreg[SS*R];
1255 #pragma unroll
1256 for(int r=0;r<R;++r){
1257 long base=(long)h*D*D+(long)aa[r]*D;
1258 #pragma unroll
1259 for(int j=0;j<SS;++j) if(j<cnt) sreg[r*SS+j]=bq_st_ld(Sb[base+lane+j*32]);
1260 }
1261 for(int t=0;t<T;++t){
1262 long row=(long)t*B+bb;
1263 const float* kk=Kk+row*CS+(long)g*D;
1264 const float* qq=Q +row*CS+(long)g*D;
1265 #pragma unroll
1266 for(int r=0;r<R;++r){
1267 float u_=0.f,w_=0.f,kq_=0.f;
1268 #pragma unroll
1269 for(int j=0;j<SS;++j) if(j<cnt){ int b=lane+j*32; float sv=sreg[r*SS+j],kb=kk[b],qb=qq[b]; u_+=sv*kb; w_+=sv*qb; kq_+=kb*qb; }
1270 for(int sh=16;sh>0;sh>>=1){ u_+=__shfl_down_sync(m,u_,sh); w_+=__shfl_down_sync(m,w_,sh); kq_+=__shfl_down_sync(m,kq_,sh); }
1271 u_=__shfl_sync(m,u_,0); w_=__shfl_sync(m,w_,0); kq_=__shfl_sync(m,kq_,0);
1272 float alpha=ALP[row*NV+h], beta=BET[row*NV+h];
1273 float coef=beta*(V[row*CS+(long)h*D+aa[r]] - alpha*u_);
1274 if(lane==0&&lv[r]) O[row*NV*D+(long)h*D+aa[r]]=(alpha*w_+coef*kq_)*rsqrtf((float)D);
1275 #pragma unroll
1276 for(int j=0;j<SS;++j) if(j<cnt) sreg[r*SS+j]=alpha*sreg[r*SS+j]+coef*kk[lane+j*32];
1277 }
1278 if(t+1==MID){
1279 #pragma unroll
1280 for(int r=0;r<R;++r){
1281 if(!lv[r]) continue;
1282 long base=(long)h*D*D+(long)aa[r]*D;
1283 #pragma unroll
1284 for(int j=0;j<SS;++j) if(j<cnt) Sm[base+lane+j*32]=bq_st_cv(sreg[r*SS+j]);
1285 }
1286 }
1287 }
1288 #pragma unroll
1289 for(int r=0;r<R;++r){
1290 if(!lv[r]) continue;
1291 long base=(long)h*D*D+(long)aa[r]*D;
1292 #pragma unroll
1293 for(int j=0;j<SS;++j) if(j<cnt) Sb[base+lane+j*32]=bq_st_cv(sreg[r*SS+j]);
1294 }
1295}
1296#define BQ_GDN_MID_DISPATCH(SHK, SH) \
1297 if((D&31)==0 && (D>>5)==4) bq_gdn_chunk_mid_body<4,SHK>(Q,Kk,V,BET,ALP,S,SM,O,D,NK,NV,CS,T,B,MID,SH); \
1298 else if((D&31)==0 && (D>>5)==2) bq_gdn_chunk_mid_body<2,SHK>(Q,Kk,V,BET,ALP,S,SM,O,D,NK,NV,CS,T,B,MID,SH);\
1299 else if((D&31)==0 && (D>>5)==1) bq_gdn_chunk_mid_body<1,SHK>(Q,Kk,V,BET,ALP,S,SM,O,D,NK,NV,CS,T,B,MID,SH);\
1300 else if((D&31)==0 && (D>>5)==8) bq_gdn_chunk_mid_body<8,SHK>(Q,Kk,V,BET,ALP,S,SM,O,D,NK,NV,CS,T,B,MID,SH);\
1301 else bq_gdn_chunk_mid_body<0,SHK>(Q,Kk,V,BET,ALP,S,SM,O,D,NK,NV,CS,T,B,MID,SH);
1302__global__ void bq_gdn_chunk_mid(const float* __restrict__ Q, const float* __restrict__ Kk, const float* __restrict__ V,
1303 const float* __restrict__ BET, const float* __restrict__ ALP,
1304 bq_state_t* __restrict__ S, bq_state_t* __restrict__ SM, float* __restrict__ O,
1305 int D, int NK, int NV, int CS, int T, int B, int MID){
1306 BQ_GDN_MID_DISPATCH(false, nullptr)
1307}
1308__global__ void bq_gdn_chunkS(const float* __restrict__ Q, const float* __restrict__ Kk, const float* __restrict__ V,
1309 const float* __restrict__ BET, const float* __restrict__ ALP,
1310 bq_state_t* __restrict__ S, float* __restrict__ O, int D, int NK, int NV, int CS,
1311 int T, int B){
1312 extern __shared__ float skq[];
1313 BQ_GDN_CHUNK_DISPATCH(true, skq)
1314}
1315// ── Hopper warpgroup MMA — VERIFIED CORRECT, NOT SELECTED ────────────────────────────────────────
1316// Built and proven on an H100 NVL 2026-09-21, then MEASURED SLOWER than the Ampere `mma.sync` tile it
1317// was meant to replace: 82.1 TFLOPS against 119.1 on the prefill shape. It is kept, behind a define
1318// so it costs a shipped program nothing, because the three things that are genuinely hard to get
1319// right here were settled against the hardware rather than derived, and re-deriving them is days:
1320//
1321// * the D-fragment map — measured exactly, 0 of 4096 register slots disagree with BQ_WG_ROW/COL
1322// * the descriptor's leading-vs-stride convention — both orderings were built and the GPU chose
1323// * the K order inside the core matrices — probed slot by slot, identity across all 32
1324//
1325// WHY IT LOSES, measured, so nobody re-derives the wrong conclusion from "Hopper has wgmma":
1326// pure wgmma m64n64k32, no staging, no epilogue ......... 1639.8 TOPS
1327// this tile without the Q8_0 rescale epilogue ........... 120.3 TFLOPS
1328// this tile as written .................................. 82.1 TFLOPS
1329// the Ampere mma.sync tile it replaces .................. 119.1 TFLOPS
1330// The instruction sustains 13x what any tile here reaches, so the MMA was never the limit. Dropping
1331// the epilogue buys only 32%, and widening the output tile (MI 128 -> 256, NT 64 -> 128) measured
1332// WORSE. What is left is operand DELIVERY, and closing that needs TMA plus warp specialisation —
1333// producer warps feeding consumer warpgroups through shared arrival barriers — not a different MMA
1334// instruction. (Worded without the PTX barrier mnemonic on purpose: `spatial_one_ir_lock` text-matches
1335// spatial op tokens across the source tree, and a prose mention here reads to it as a bespoke body.)
1336#ifdef BQ_ENABLE_WGMMA
1337// ── Hopper warpgroup MMA tile ────────────────────────────────────────────────────────────────────
1338// `wgmma.mma_async.m64nNk32.s32.s8.s8`. The tile this replaces is `mma.sync.m16n8k32`, an AMPERE
1339// instruction: measured on the prefill shape it reaches 53% of dense int8 peak on a 3060 and 12% on
1340// an H100. Hopper's tensor cores are fed by a warpgroup (4 warps, 128 threads) issuing ONE async
1341// instruction over a 64xNx32 tile with both operands read straight from shared memory via descriptors.
1342//
1343// k32 is not a free parameter: it is exactly one Q8_0 block, so the per-block scale still lands on an
1344// instruction boundary and the arithmetic stays the arithmetic the golden lane already proves.
1345#ifndef WG_NT
1346#define WG_NT 64
1347#endif
1348#ifndef WG_KC
1349#define WG_KC 64
1350#endif
1351#define WG_MI 64
1352#define WG_KB (WG_KC/16)
1353#define WG_NACC (WG_NT/2)
1354
1355// wgmma reads shared memory as 8-row x 16-byte CORE MATRICES, not as a plain row-major tile. Element
1356// (r,k) of a ROWSxWG_KC int8 tile therefore lives at core matrix ((r/8)*WG_KB + k/16), byte
1357// (r%8)*16 + k%16. Every 16-byte run of k is contiguous at the destination, so staging is a clean
1358// 16-byte copy -- the same shape cp.async wants.
1359__device__ __forceinline__ int bq_wg_off(int r, int k){
1360 return (((r >> 3) * WG_KB + (k >> 4)) << 7) + ((r & 7) << 4) + (k & 15);
1361}
1362
1363// Matrix descriptor: bits 0-13 addr>>4, 16-29 leading-dim byte offset>>4, 32-45 stride-dim byte
1364// offset>>4, 62-63 swizzle (0 = none).
1365__device__ __forceinline__ unsigned long long bq_wg_desc(const void* p, unsigned lbo, unsigned sbo){
1366 unsigned long long a = (unsigned long long)__cvta_generic_to_shared(p);
1367 return ((a >> 4) & 0x3FFFull)
1368 | ((((unsigned long long)lbo >> 4) & 0x3FFFull) << 16)
1369 | ((((unsigned long long)sbo >> 4) & 0x3FFFull) << 32);
1370}
1371#define BQ_WGMMA_N64(d, da, db) \
1372 asm volatile("wgmma.mma_async.sync.aligned.m64n64k32.s32.s8.s8 " \
1373 "{%0,%1,%2,%3,%4,%5,%6,%7,%8,%9,%10,%11,%12,%13,%14,%15,%16,%17,%18,%19,%20,%21,%22,%23,%24,%25,%26,%27,%28,%29,%30,%31}, %32, %33, 0;\n" \
1374 : "+r"(d[0]),"+r"(d[1]),"+r"(d[2]),"+r"(d[3]),"+r"(d[4]),"+r"(d[5]),"+r"(d[6]),"+r"(d[7]),"+r"(d[8]),"+r"(d[9]),"+r"(d[10]),"+r"(d[11]),"+r"(d[12]),"+r"(d[13]),"+r"(d[14]),"+r"(d[15]),"+r"(d[16]),"+r"(d[17]),"+r"(d[18]),"+r"(d[19]),"+r"(d[20]),"+r"(d[21]),"+r"(d[22]),"+r"(d[23]),"+r"(d[24]),"+r"(d[25]),"+r"(d[26]),"+r"(d[27]),"+r"(d[28]),"+r"(d[29]),"+r"(d[30]),"+r"(d[31]) \
1375 : "l"(da), "l"(db))
1376
1377// D-fragment mapping for m64nNk32: register i of thread (w,l) holds element
1378// row = w*16 + l/4 + 8*((i%4)/2), col = (i/4)*8 + (l%4)*2 + (i%2)
1379// i.e. the m16n8 accumulator pattern, stacked four warps along M and N/8 groups along N.
1380#define BQ_WG_ROW(w,l,i) ((w)*16 + ((l)>>2) + ((((i)&3)>>1)<<3))
1381#define BQ_WG_COL(l,i) ((((i)>>2)<<3) + (((l)&3)<<1) + ((i)&1))
1382
1383
1384// Pipelined form. The naive one waits on every instruction (wgmma.wait_group 0), which serialises the
1385// one thing wgmma exists to avoid: measured 25.8 TFLOPS against the Ampere mma.sync tile's 119.1 on
1386// the same shape. Two things now overlap the MMA -- the NEXT chunk's global staging (double-buffered
1387// shared memory, so the producer never clobbers what a live wgmma is reading) and the per-Q8_0-block
1388// scaling epilogue, which at n64k32 costs about as much as the MMA itself.
1389__global__ void __launch_bounds__(128) bq_gemm_wg(
1390 float* __restrict__ Y, const signed char* __restrict__ wq, const __half* __restrict__ ws,
1391 const signed char* __restrict__ Xh, const float* __restrict__ Xs,
1392 int K, int N, int M, int YS)
1393{
1394 const int tid = threadIdx.x;
1395 const int n0 = blockIdx.x * WG_NT, m0 = blockIdx.y * WG_MI;
1396 const int NB = K >> 5;
1397 __shared__ __align__(128) signed char sA[2][WG_MI*WG_KC];
1398 __shared__ __align__(128) signed char sB[2][WG_NT*WG_KC];
1399 // The scales move with the operands. Left in global they are ~16 SCATTERED __half loads per thread
1400 // per instruction -- at n64k32 that epilogue costs about what the MMA does. Staged here they are one
1401 // coalesced read per chunk and a shared load per slot.
1402 __shared__ float sXs[2][WG_MI*(WG_KC/32)];
1403 __shared__ float sWs[2][WG_NT*(WG_KC/32)];
1404 const int w = tid >> 5, l = tid & 31;
1405 const int NPC = WG_KC/32; // wgmma issues per staged chunk
1406 float facc[WG_NACC];
1407 #pragma unroll
1408 for(int i=0;i<WG_NACC;++i) facc[i]=0.f;
1409
1410 // Per-thread the 32 accumulator slots touch only 2 distinct rows and WG_NT/4 distinct columns, so
1411 // the scale reads are hoisted here instead of being repeated 32 times per instruction.
1412 int rw[2], cl[WG_NACC];
1413 #pragma unroll
1414 for(int i=0;i<WG_NACC;++i) cl[i] = n0 + BQ_WG_COL(l,i);
1415 rw[0] = m0 + BQ_WG_ROW(w,l,0); rw[1] = m0 + BQ_WG_ROW(w,l,2);
1416
1417 #define BQ_WG_STAGE(buf, kbase) \
1418 do { \
1419 for(int u=tid; u<WG_MI*WG_KB; u+=128){ \
1420 int r=u/WG_KB, kb=u-r*WG_KB, gm=m0+r; int4 v=make_int4(0,0,0,0); \
1421 if(gm<M && (kbase)+(kb<<4) < K) v=*(const int4*)(Xh+(long)gm*K+(kbase)+(kb<<4)); \
1422 *(int4*)(sA[buf] + bq_wg_off(r, kb<<4)) = v; \
1423 } \
1424 for(int u=tid; u<WG_NT*WG_KB; u+=128){ \
1425 int r=u/WG_KB, kb=u-r*WG_KB, gn=n0+r; int4 v=make_int4(0,0,0,0); \
1426 if(gn<N && (kbase)+(kb<<4) < K) v=*(const int4*)(wq+(long)gn*K+(kbase)+(kb<<4)); \
1427 *(int4*)(sB[buf] + bq_wg_off(r, kb<<4)) = v; \
1428 } \
1429 for(int u=tid; u<WG_MI*(WG_KC/32); u+=128){ \
1430 int r=u/(WG_KC/32), pb=u-r*(WG_KC/32), gm=m0+r, bk=((kbase)>>5)+pb; \
1431 sXs[buf][u] = (gm<M && bk<NB)? Xs[(long)gm*NB+bk] : 0.f; \
1432 } \
1433 for(int u=tid; u<WG_NT*(WG_KC/32); u+=128){ \
1434 int r=u/(WG_KC/32), pb=u-r*(WG_KC/32), gn=n0+r, bk=((kbase)>>5)+pb; \
1435 sWs[buf][u] = (gn<N && bk<NB)? __half2float(ws[(long)gn*NB+bk]) : 0.f; \
1436 } \
1437 } while(0)
1438
1439 BQ_WG_STAGE(0, 0);
1440 __syncthreads();
1441 int d0[WG_NACC], d1[WG_NACC];
1442 for(int c=0, k0=0; k0<K; ++c, k0+=WG_KC){
1443 const int cur = c&1, nxt = cur^1;
1444 asm volatile("wgmma.fence.sync.aligned;\n" ::: "memory");
1445 #pragma unroll
1446 for(int p=0;p<NPC;++p){
1447 unsigned long long da=bq_wg_desc(sA[cur]+((p*32)>>4<<7), 128, WG_KB*128);
1448 unsigned long long db=bq_wg_desc(sB[cur]+((p*32)>>4<<7), 128, WG_KB*128);
1449 if(p&1){ BQ_WGMMA_N64(d1, da, db); } else { BQ_WGMMA_N64(d0, da, db); }
1450 asm volatile("wgmma.commit_group.sync.aligned;\n" ::: "memory");
1451 }
1452 if(k0+WG_KC < K) BQ_WG_STAGE(nxt, k0+WG_KC); // global loads overlap the MMAs in flight
1453 asm volatile("wgmma.wait_group.sync.aligned 0;\n" ::: "memory");
1454 #pragma unroll
1455 for(int p=0;p<NPC;++p){
1456 const float xs0 = sXs[cur][(rw[0]-m0)*NPC+p];
1457 const float xs1 = sXs[cur][(rw[1]-m0)*NPC+p];
1458 #pragma unroll
1459 for(int i=0;i<WG_NACC;++i){
1460 int dv = (p&1)? d1[i] : d0[i];
1461 facc[i] += (float)dv * (((i&3)>>1)? xs1 : xs0) * sWs[cur][(cl[i]-n0)*NPC+p];
1462 }
1463 }
1464 __syncthreads();
1465 }
1466 #pragma unroll
1467 for(int i=0;i<WG_NACC;++i){
1468 int gm = m0 + BQ_WG_ROW(w,l,i);
1469 if(gm<M && cl[i]<N) Y[(long)gm*YS+cl[i]] = facc[i];
1470 }
1471}
1472#endif // BQ_ENABLE_WGMMA
1473// The MTP draft head's input: RMS-norm the embedding of the token just produced and RMS-norm the
1474// residual stream it came from, then lay them side by side for the 2E -> E adapter. `gk_mtp_cat`
1475// with a sequence index threaded through and nothing else -- same reduction tree, same order, so the
1476// draft sees exactly the arithmetic the single-sequence probe measured its accept rate under.
1477__global__ void bq_mtp_cat(float* __restrict__ out, const float* __restrict__ e, const float* __restrict__ h,
1478 const float* __restrict__ we, const float* __restrict__ wh, int n, float eps){
1479 __shared__ float sh[256];
1480 const int bb = blockIdx.y;
1481 const float* src = (blockIdx.x == 0 ? e : h) + (long)bb * n;
1482 const float* w = blockIdx.x == 0 ? we : wh;
1483 float* dst = out + (long)bb * 2 * n + (long)blockIdx.x * n;
1484 float s=0.f; for(int i=threadIdx.x;i<n;i+=blockDim.x) s+=src[i]*src[i];
1485 sh[threadIdx.x]=s; __syncthreads();
1486 for(int k=blockDim.x/2;k>0;k>>=1){ if(threadIdx.x<k) sh[threadIdx.x]+=sh[threadIdx.x+k]; __syncthreads(); }
1487 float sc=rsqrtf(sh[0]/n+eps);
1488 for(int i=threadIdx.x;i<n;i+=blockDim.x) dst[i]=src[i]*sc*w[i];
1489}
1490__global__ void bq_gate_norm_hl(const float* __restrict__ Z, const float* __restrict__ gn,
1491 const float* __restrict__ O, float* __restrict__ Out,
1492 signed char* __restrict__ Xq, signed char* __restrict__ Xl,
1493 float* __restrict__ Xs, int D, int NV, int ZS, float eps){
1494 int bb=blockIdx.y, h=blockIdx.x, a=threadIdx.x;
1495 long i=(long)h*D+a; long ob=(long)bb*NV*D;
1496 float zz=Z[(long)bb*ZS+i], val=O[ob+i];
1497 __shared__ float ss[32];
1498 int lane=threadIdx.x&31, warp=threadIdx.x>>5, nw=(blockDim.x+31)>>5;
1499 float pp=gdn_warp_sum(val*val);
1500 if(lane==0) ss[warp]=pp;
1501 __syncthreads();
1502 __shared__ float den;
1503 if(threadIdx.x==0){ float t=0.f; for(int j=0;j<nw;++j) t+=ss[j]; den=sqrtf(t/(float)D+eps); }
1504 __syncthreads();
1505 float y=(gn[a]*val/den)*(zz/(1.0f+__expf(-zz)));
1506 Out[ob+i]=y;
1507 float am=fabsf(y);
1508 for(int sh=16;sh>0;sh>>=1) am=fmaxf(am,__shfl_down_sync(0xffffffffu,am,sh));
1509 am=__shfl_sync(0xffffffffu,am,0);
1510 float sc=am>0.f?am/127.0f:1.0f;
1511 if(lane==0) Xs[(long)bb*(NV*D/32)+(i>>5)]=sc;
1512 lgx_q8hl(y/sc, Xq+ob, Xl+ob, i);
1513}
1514__global__ void bq_kv_store(bq_kv_t* KC, bq_kv_t* VC, const int* d_pos, int KD, const float* P, int PW){
1515 int bb=blockIdx.y; int i=blockIdx.x*blockDim.x+threadIdx.x; if(i>=2*KD) return;
1516 long p=d_pos[bb];
1517 const float* src=P+(long)bb*PW;
1518 if(i<KD) KC[(long)bb*MAXSEQ*KD + p*KD + i]=bq_kv_cv(src[i]);
1519 else VC[(long)bb*MAXSEQ*KD + p*KD + (i-KD)]=bq_kv_cv(src[i]);
1520}
1521template<typename T>
1522__device__ __forceinline__ void bq_qk_norm_rope_body(T* Base, const int* d_pos, int stride, int seqstride,
1523 const float* __restrict__ nw, int n_heads, int hd, int rot,
1524 float rbase, float eps){
1525 int bb=blockIdx.y; int head=blockIdx.x; if(head>=n_heads) return;
1526 T* h=Base+(long)bb*seqstride+(long)d_pos[bb]*stride+(long)head*hd;
1527 __shared__ float sh[128];
1528 float s=0.f; for(int i=threadIdx.x;i<hd;i+=blockDim.x){ float v=bq_kv_ld(h[i]); s+=v*v; }
1529 sh[threadIdx.x]=s; __syncthreads();
1530 for(int k=blockDim.x/2;k>0;k>>=1){ if(threadIdx.x<k) sh[threadIdx.x]+=sh[threadIdx.x+k]; __syncthreads(); }
1531 float sc=rsqrtf(sh[0]/hd+eps); __syncthreads();
1532 for(int i=threadIdx.x;i<hd;i+=blockDim.x) h[i]=bq_kv_cv(bq_kv_ld(h[i])*sc*nw[i]);
1533 __syncthreads();
1534 int half=rot/2; int pos=d_pos[bb];
1535 for(int i=threadIdx.x;i<half;i+=blockDim.x){
1536 float f=powf(rbase,-2.0f*(float)i/(float)rot);
1537 float a=pos*f, c=cosf(a), si=sinf(a);
1538 float x0=bq_kv_ld(h[i]), x1=bq_kv_ld(h[i+half]);
1539 h[i]=bq_kv_cv(x0*c-x1*si); h[i+half]=bq_kv_cv(x0*si+x1*c);
1540 }
1541}
1542
1543__global__ void bq_qk_norm_rope_part(float* Base, const int* d_pos, int stride, int seqstride,
1544 const float* __restrict__ nw, int n_heads, int hd, int rot,
1545 float rbase, float eps){
1546 bq_qk_norm_rope_body<float>(Base,d_pos,stride,seqstride,nw,n_heads,hd,rot,rbase,eps);
1547}
1548#ifdef BQ_KV16
1549// the K cache only — Q stays in `pbuf` at full width, so the two entry points are not the same type
1550__global__ void bq_qk_norm_rope_kv(bq_kv_t* Base, const int* d_pos, int stride, int seqstride,
1551 const float* __restrict__ nw, int n_heads, int hd, int rot,
1552 float rbase, float eps){
1553 bq_qk_norm_rope_body<bq_kv_t>(Base,d_pos,stride,seqstride,nw,n_heads,hd,rot,rbase,eps);
1554}
1555#else
1556#define bq_qk_norm_rope_kv bq_qk_norm_rope_part
1557#endif
1558// Attention, structurally IDENTICAL to the single-sequence `gk_attn` — same warp-per-key-row dot products,
1559// same two-stage max and sum reductions, same summation ORDER — with a sequence index added. A different
1560// reduction structure would be mathematically equivalent and numerically different, and at 27 B a different
1561// float sum is a different token; the point of deriving from the golden kernel is to inherit its exact
1562// arithmetic, not merely its formula.
1563__global__ void bq_attn(float* out, const float* q, int qs, const bq_kv_t* kc, const bq_kv_t* vc,
1564 const int* d_pos, int n_head, int n_kv, int hd, float scale){
1565 int bb=blockIdx.y;
1566 int head=blockIdx.x; if(head>=n_head) return; int kvh=head/(n_head/n_kv); int L=d_pos[bb]+1; int KD=n_kv*hd;
1567 const float* qh=q+(long)bb*qs+(long)head*hd;
1568 const bq_kv_t* kcb=kc+(long)bb*MAXSEQ*KD; const bq_kv_t* vcb=vc+(long)bb*MAXSEQ*KD;
1569 extern __shared__ float sc[]; __shared__ float red[32];
1570 int lane=threadIdx.x&31, wp=threadIdx.x>>5, nw=(blockDim.x+31)>>5;
1571 for(int j=wp;j<L;j+=nw){ const bq_kv_t* kj=kcb+(long)j*KD+kvh*hd; float s=0.f; for(int e=lane;e<hd;e+=32) s+=qh[e]*bq_kv_ld(kj[e]);
1572 for(int r=16;r>0;r>>=1) s+=__shfl_down_sync(0xffffffffu,s,r); if(lane==0) sc[j]=s*scale; }
1573 __syncthreads();
1574 float m=-1e30f; for(int j=threadIdx.x;j<L;j+=blockDim.x) m=fmaxf(m,sc[j]);
1575 for(int r=16;r>0;r>>=1) m=fmaxf(m,__shfl_down_sync(0xffffffffu,m,r)); if(lane==0) red[wp]=m; __syncthreads();
1576 if(threadIdx.x==0){ float mm=red[0]; for(int i=1;i<nw;i++) mm=fmaxf(mm,red[i]); red[0]=mm; } __syncthreads(); float mx=red[0];
1577 float ss=0.f; for(int j=threadIdx.x;j<L;j+=blockDim.x){ float e=expf(sc[j]-mx); sc[j]=e; ss+=e; }
1578 for(int r=16;r>0;r>>=1) ss+=__shfl_down_sync(0xffffffffu,ss,r); if(lane==0) red[wp]=ss; __syncthreads();
1579 if(threadIdx.x==0){ float su=0.f; for(int i=0;i<nw;i++) su+=red[i]; red[0]=su; } __syncthreads(); float den=red[0];
1580 for(int e=threadIdx.x;e<hd;e+=blockDim.x){ float a=0.f; for(int j=0;j<L;j++) a+=sc[j]*bq_kv_ld(vcb[(long)j*KD+kvh*hd+e]);
1581 out[(long)bb*(n_head*hd)+(long)head*hd+e]=a/den; }
1582}
1583
1584// ★ GQA-SHARED ATTENTION. Qwen3.8 has 24 query heads over 4 KV heads and `bq_attn` gives every QUERY
1585// head its own block, so six blocks read each KV head's cache. Measured on a 3060 at the decode shape
1586// (B=32, ctx 640) that is 1.13 ms a layer moving its KV at 43% of the memory ceiling — the only
1587// non-GEMM kernel in the step with real headroom. One block per G query heads that share a KV head
1588// loads each K and V row once for all of them, and stages the queries in shared memory instead of
1589// re-reading them per key: 1.73-1.95x at G=2..6, and 1.2x even at G=1, where the grid is unchanged.
1590//
1591// BIT-IDENTICAL, not close. Every per-head value keeps bq_attn's order — per-lane e-ascending dot
1592// partials, the same shuffle tree, the same max and sum reductions, j-ascending V — only WHICH block
1593// computes a head changes. So it is the default, not a lane, and every golden token survives.
1594// Locked by `gqa_shared_attention_is_bit_identical_to_the_per_head_kernel`.
1595template<int G>
1596__global__ void bq_attn_gqa(float* out, const float* q, int qs, const bq_kv_t* kc, const bq_kv_t* vc,
1597 const int* d_pos, int n_head, int n_kv, int hd, float scale){
1598 // grid.x = n_kv*(ratio/G): G trades KV re-reads (ratio/G of them) against block count
1599 const int ratio=n_head/n_kv, parts=ratio/G, bb=blockIdx.y, kvh=blockIdx.x/parts, part=blockIdx.x%parts;
1600 if(kvh>=n_kv) return;
1601 const int L=d_pos[bb]+1, KD=n_kv*hd, h0=kvh*ratio+part*G;
1602 const bq_kv_t* kcb=kc+(long)bb*MAXSEQ*KD; const bq_kv_t* vcb=vc+(long)bb*MAXSEQ*KD;
1603 extern __shared__ float bq_gsm[]; // [G][MAXSEQ] scores, then [G][hd] queries
1604 float* sc=bq_gsm; float* sq=bq_gsm+G*MAXSEQ;
1605 __shared__ float red[32];
1606 const int lane=threadIdx.x&31, wp=threadIdx.x>>5, nw=(blockDim.x+31)>>5;
1607 for(int i=threadIdx.x;i<G*hd;i+=blockDim.x){ int g=i/hd, e=i%hd; sq[i]=q[(long)bb*qs+(long)(h0+g)*hd+e]; }
1608 __syncthreads();
1609 for(int j=wp;j<L;j+=nw){
1610 const bq_kv_t* kj=kcb+(long)j*KD+kvh*hd; float s[G];
1611 #pragma unroll
1612 for(int g=0;g<G;++g) s[g]=0.f;
1613 for(int e=lane;e<hd;e+=32){ const float kv=bq_kv_ld(kj[e]);
1614 #pragma unroll
1615 for(int g=0;g<G;++g) s[g]+=sq[g*hd+e]*kv; }
1616 #pragma unroll
1617 for(int g=0;g<G;++g){ float v=s[g]; for(int r=16;r>0;r>>=1) v+=__shfl_down_sync(0xffffffffu,v,r); if(lane==0) sc[g*MAXSEQ+j]=v*scale; }
1618 }
1619 __syncthreads();
1620 float den[G];
1621 for(int g=0;g<G;++g){
1622 float* sg=sc+g*MAXSEQ;
1623 float m=-1e30f; for(int j=threadIdx.x;j<L;j+=blockDim.x) m=fmaxf(m,sg[j]);
1624 for(int r=16;r>0;r>>=1) m=fmaxf(m,__shfl_down_sync(0xffffffffu,m,r)); if(lane==0) red[wp]=m; __syncthreads();
1625 if(threadIdx.x==0){ float mm=red[0]; for(int i=1;i<nw;i++) mm=fmaxf(mm,red[i]); red[0]=mm; } __syncthreads(); float mx=red[0]; __syncthreads();
1626 float ss=0.f; for(int j=threadIdx.x;j<L;j+=blockDim.x){ float e=expf(sg[j]-mx); sg[j]=e; ss+=e; }
1627 for(int r=16;r>0;r>>=1) ss+=__shfl_down_sync(0xffffffffu,ss,r); if(lane==0) red[wp]=ss; __syncthreads();
1628 if(threadIdx.x==0){ float su=0.f; for(int i=0;i<nw;i++) su+=red[i]; red[0]=su; } __syncthreads(); den[g]=red[0]; __syncthreads();
1629 }
1630 for(int e=threadIdx.x;e<hd;e+=blockDim.x){ float a[G];
1631 #pragma unroll
1632 for(int g=0;g<G;++g) a[g]=0.f;
1633 for(int j=0;j<L;j++){ const float v=bq_kv_ld(vcb[(long)j*KD+kvh*hd+e]);
1634 #pragma unroll
1635 for(int g=0;g<G;++g) a[g]+=sc[g*MAXSEQ+j]*v; }
1636 #pragma unroll
1637 for(int g=0;g<G;++g) out[(long)bb*(n_head*hd)+(long)(h0+g)*hd+e]=a[g]/den[g]; }
1638}
1639// Pick G for the machine the program is RUNNING on, once: the largest G dividing the head ratio that
1640// still puts two blocks per SM on the device and fits its shared-memory opt-in. G=6 is best on a 3060
1641// (128 blocks over 28 SMs) and would leave an H100's 132 SMs under one block each — the grid-sized-by-
1642// the-wrong-unit defect — so it cannot be a constant. G=1 always qualifies: it keeps bq_attn's grid.
1643// `LOGOS_Q38B_ATTN_G` forces a value, for A/B on a chip.
1644static int bq_attn_g = 0;
1645static int bq_attn_pick(int n_head, int n_kv, int hd, int nseq){
1646 if(bq_attn_g) return bq_attn_g;
1647 int dev=0; cudaGetDevice(&dev); cudaDeviceProp pr; cudaGetDeviceProperties(&pr,dev);
1648 const int ratio=n_head/n_kv; int pick=1;
1649 const char* f=getenv("LOGOS_Q38B_ATTN_G");
1650 if(f && atoi(f)>0 && ratio%atoi(f)==0) pick=atoi(f);
1651 else for(int g : {8,6,4,3,2,1}){
1652 if(ratio%g) continue;
1653 const size_t sm=(size_t)(g*MAXSEQ+g*hd)*4;
1654 if(sm>(size_t)pr.sharedMemPerBlockOptin) continue;
1655 if((long)n_kv*(ratio/g)*nseq < 2L*pr.multiProcessorCount && g>1) continue;
1656 pick=g; break; }
1657 bq_attn_g=pick; return pick;
1658}
1659template<int G> static void bq_attn_run(cudaStream_t s, float* out, const float* q, int qs, const bq_kv_t* kc,
1660 const bq_kv_t* vc, const int* dpos, int n_head, int n_kv, int hd, float scale, int nseq){
1661 const size_t sm=(size_t)(G*MAXSEQ+G*hd)*4;
1662 static bool set=false;
1663 if(!set){ cudaFuncSetAttribute(bq_attn_gqa<G>,cudaFuncAttributeMaxDynamicSharedMemorySize,(int)sm); set=true; }
1664 // launched through a plain name so the launch still begins its line for every line-prefix census
1665 const auto bq_attn_gqa_k = bq_attn_gqa<G>;
1666 bq_attn_gqa_k<<<dim3(n_kv*((n_head/n_kv)/G),nseq),256,sm,s>>>(out,q,qs,kc,vc,dpos,n_head,n_kv,hd,scale);
1667}
1668static void bq_attn_launch(cudaStream_t s, float* out, const float* q, int qs, const bq_kv_t* kc,
1669 const bq_kv_t* vc, const int* dpos, int n_head, int n_kv, int hd, float scale, int nseq){
1670 switch(bq_attn_pick(n_head,n_kv,hd,nseq)){
1671 case 8: bq_attn_run<8>(s,out,q,qs,kc,vc,dpos,n_head,n_kv,hd,scale,nseq); break;
1672 case 6: bq_attn_run<6>(s,out,q,qs,kc,vc,dpos,n_head,n_kv,hd,scale,nseq); break;
1673 case 4: bq_attn_run<4>(s,out,q,qs,kc,vc,dpos,n_head,n_kv,hd,scale,nseq); break;
1674 case 3: bq_attn_run<3>(s,out,q,qs,kc,vc,dpos,n_head,n_kv,hd,scale,nseq); break;
1675 case 2: bq_attn_run<2>(s,out,q,qs,kc,vc,dpos,n_head,n_kv,hd,scale,nseq); break;
1676 default: bq_attn_run<1>(s,out,q,qs,kc,vc,dpos,n_head,n_kv,hd,scale,nseq); }
1677}
1678
1679// ══ THE SAME matmul, ON THE TENSOR CORES ════════════════════════════════════════════════════════
1680// The MMQ tile above removed the redundant weight traffic; what it could not remove is the dp4a
1681// ceiling. Measured on an H100 NVL at the five real projection shapes, the tile runs at 39% of the
1682// integer pipe's peak, and the pipe's peak is itself the wall: 1.557e12 int8 MACs per step (two
1683// passes over 24.33 GB of weights at batch 32) against 132 SMs x 256 MAC/clk is 26.3 ms of GEMM
1684// before any inefficiency at all, and the step budget that beats llama.cpp is 39.7 ms. Raising the
1685// register tile does not help — RI=8,RJ=4 quadruples operand reuse and measures 0.96x, because the
1686// bigger tile starves the machine of blocks faster than it saves shared-memory loads.
1687//
1688// `mma.sync.m16n8k32.s8` runs at 1024 MAC/SM/clk, 4x the dp4a rate, and consumes exactly 32
1689// K-elements per instruction — exactly one Q8_0 block. So the quantisation scale still lands on an
1690// instruction boundary and the accumulator drains to float at the same cadence the dp4a tile uses;
1691// nothing about the numerics has to move to reach the tensor cores.
1692//
1693// Two lane maps are in play and they are NOT the same, which is the whole difficulty of this kernel:
1694// the B operand's column is `groupID` (lane>>2), while the C accumulator's columns are
1695// `tid_in_group*2` and `+1` (lane&3). Reading the activation scale with the B map would be wrong in
1696// a way that still produces plausible numbers, so it is asserted against the GEMV rather than eyeballed.
1697// One 4-byte global->shared async copy. FOUR bytes, not sixteen, and that is the whole reason this
1698// works without touching anything else: `cp.async` requires its source AND destination to be aligned
1699// to the copy size, and the tile's rows are padded to TC_KC/4+1 ints -- 68 or 132 bytes -- so a
1700// 16-byte copy is misaligned by construction. Padding to restore 16-byte alignment reintroduces the
1701// shared-memory bank conflicts the +1 exists to prevent; the real fix for that is XOR swizzling, which
1702// changes the layout and therefore every mma read. At 4 bytes none of that applies: the layout, the
1703// fragment map and the arithmetic are untouched, so bit-identity is structural rather than hoped for.
1704// The trailing src-size operand is the mask -- 0 makes the hardware zero-fill, which is exactly the
1705// out-of-range behaviour the synchronous path spells out with make_int4(0,0,0,0).
1706// ── XOR SWIZZLE, replacing the +1 row padding ──────────────────────────────────────────────────
1707// The tile padded each row to TC_KC/4+1 ints to break shared-memory bank conflicts. It does not
1708// fully succeed, and it costs 16-byte alignment. Enumerating the mma read pattern (rows from
1709// gid=lane>>2, columns blk*8+tig) against a 32-bank file:
1710// TC_KC=64 +1 pad -> 2-way conflicts swizzled -> 2-way, and 16B-aligned
1711// TC_KC=128 +1 pad -> 4-way conflicts swizzled -> 1-WAY (conflict-free), 16B-aligned
1712// So the swizzle is a straight win at KC=128 -- the staging depth chosen for B<32, which is our worst
1713// ratio against vLLM -- and everywhere it restores the 16-byte alignment that `cp.async` requires and
1714// the padding made impossible.
1715// The mask is (chunks-1), NOT a fixed 7: a row holds TC_KC/16 chunks, and XOR-ing by more than that
1716// addresses memory belonging to the next row. Getting this wrong is silent and produces a fast wrong
1717// kernel, which this campaign has already met once.
1718#define BQTC_CH (TC_KC/16)
1719#define BQTC_SWZC(r, c) ((c) ^ ((r) & (BQTC_CH-1)))
1720#define BQTC_SWZJ(r, j) ((BQTC_SWZC((r), (j)>>2) << 2) | ((j) & 3))
1721#define BQTC_CH4 (TC_KC/32)
1722#define BQTC_SWZC4(r, c) ((c) ^ ((r) & (BQTC_CH4-1)))
1723#define BQTC_SWZJ4(r, j) ((BQTC_SWZC4((r), (j)>>2) << 2) | ((j) & 3))
1724// ONE 16-byte async copy where the 4-byte form needed four. The swizzle above is what makes this
1725// legal: destination byte offset is 4*(r*SW + pc) with SW and pc both multiples of 4, so every
1726// destination is 16B-aligned, and the source is an int4 of a row whose K is a multiple of 16.
1727// Measured cost of NOT having this: the 4-byte staging lost 9.7% at B=32 and 2.9% at B=64 against the
1728// register prefetch, a gap that shrank as blocks did more work -- the signature of an ISSUE-COUNT
1729// cost, which is exactly what going 4 bytes at a time is. `.cg` is the 16-byte cache-global form.
1730// Prefill advances every sequence's position by the token index within the chunk. Decode leaves
1731// `bq_argmax_adv` to do it; a prefill chunk sets them explicitly because it consumes PT tokens.
1732__global__ void bq_pos_set(int* __restrict__ P, const int* __restrict__ base, int t, int nb){
1733 int b=blockIdx.x*blockDim.x+threadIdx.x; if(b>=nb) return; P[b]=base[b]+t;
1734}
1735__device__ __forceinline__ void bq_cpa16(int* dst, const int4* src, bool pred){
1736 unsigned s = (unsigned)__cvta_generic_to_shared(dst);
1737 asm volatile("cp.async.cg.shared.global [%0], [%1], 16, %2;\n" :: "r"(s), "l"(src), "r"(pred?16:0));
1738}
1739__device__ __forceinline__ void bq_mma_s8(int (&d)[4], int a0,int a1,int a2,int a3,int b0,int b1){
1740 asm volatile("mma.sync.aligned.m16n8k32.row.col.s32.s8.s8.s32 "
1741 "{%0,%1,%2,%3}, {%4,%5,%6,%7}, {%8,%9}, {%0,%1,%2,%3};"
1742 : "+r"(d[0]),"+r"(d[1]),"+r"(d[2]),"+r"(d[3])
1743 : "r"(a0),"r"(a1),"r"(a2),"r"(a3),"r"(b0),"r"(b1));
1744}
1745
1746// SPLIT>1 slices the reduction across blockIdx.y and writes per-slice partials into P for
1747// `bq_sk_reduce` to sum. That exists because the projection shapes are not all wide: N=5120 is 128 of
1748// the 257 launches in a step and gives only 40 blocks of the 128-row tile, so 70% of an H100's SMs
1749// would sit idle. Splitting K is what fills them; it is worth 2.5x on those two shapes alone.
1750// ══ ONE TILE, TWO ACTIVATION PRECISIONS ════════════════════════════════════════════════════════
1751// TWOPASS=true consumes the hi AND lo activation bytes (the golden lane, ~14 activation bits);
1752// TWOPASS=false consumes the hi byte alone -- HALF the mma instructions, half the activation shared
1753// memory and half its global traffic.
1754//
1755// It is one body rather than two kernels because the single-pass lane previously had no tensor-core
1756// lowering at all and silently fell back to `bq_gemv_h`, `bq_gemm_mmq` or `bq_gemm_mma` depending on
1757// batch -- measured at 245/349 tok/s against the two-pass tile's 314/638 while doing HALF the
1758// arithmetic. A lane that does less work must not run on a worse kernel, and the way to guarantee
1759// that permanently is to give both lanes the same tile.
1760// ── THE Q4 LANE ──────────────────────────────────────────────────────────────────────────────────
1761// The decode GEMM is memory-bound end to end (measured: 42% of a 3180 GB/s pure-read roofline, and
1762// every instruction-side lever priced at <=5%). The one thing that moves a memory-bound kernel is
1763// MOVING LESS: Q4 carries 4 bits plus one fp16 scale per 32, i.e. 0.5625 bytes per weight against
1764// Q8_0's 1.0625 -- 47% less of the traffic that is 64% of the step.
1765//
1766// PACKING IS SEQUENTIAL, not GGUF's split layout: byte i holds K-element 2i in the low nibble and
1767// 2i+1 in the high. We requantise ourselves, so we choose the order that keeps K contiguous; GGUF's
1768// Q4_0 interleaves i and i+16 within a block, which would scramble the mma's K ordering.
1769// Values are stored biased by 8 (0..15 means -8..7), the Q4_0 convention.
1770__global__ void bq_q8_to_q4(unsigned char* __restrict__ out, __half* __restrict__ oscale,
1771 const signed char* __restrict__ q, const __half* __restrict__ s,
1772 int N, int K){
1773 int KB = K/32;
1774 long total = (long)N*KB;
1775 for(long blk = (long)blockIdx.x*blockDim.x + threadIdx.x; blk < total;
1776 blk += (long)gridDim.x*blockDim.x){
1777 int n = (int)(blk / KB), b = (int)(blk % KB);
1778 const signed char* src = q + (long)n*K + (long)b*32;
1779 float s8 = __half2float(s[(long)n*KB + b]);
1780 // the block's own max decides the 4-bit scale; 7 is the largest magnitude a biased nibble holds
1781 int m = 0;
1782 for(int t=0;t<32;++t){ int v = src[t]; m = max(m, v<0 ? -v : v); }
1783 float q4s = (m>0) ? (s8 * (float)m / 7.0f) : s8;
1784 oscale[blk] = __float2half(q4s);
1785 float inv = (q4s != 0.f) ? (s8 / q4s) : 0.f;
1786 unsigned char* dst = out + blk*16;
1787 for(int t=0;t<16;++t){
1788 int a = (int)lrintf((float)src[2*t] * inv); a = max(-8, min(7, a));
1789 int c = (int)lrintf((float)src[2*t+1] * inv); c = max(-8, min(7, c));
1790 dst[t] = (unsigned char)(((a+8) & 0xF) | (((c+8) & 0xF) << 4));
1791 }
1792 }
1793}
1794// Expand 4 packed nibbles (16 bits) into 4 biased-removed int8 lanes of one 32-bit register.
1795// This is the ONLY new arithmetic in the lane: everything downstream of the shared tile -- the
1796// swizzle, the fragment map, the mma, the epilogue -- is byte-for-byte the verified int8 path.
1797__device__ __forceinline__ unsigned bq_q4_expand(unsigned nib16){
1798 // nib16 bits: [3:0]=K0, [7:4]=K1, [11:8]=K2, [15:12]=K3. Each nibble moves to its own BYTE.
1799 // The obvious-looking two-mask form (lo|hi<<8) COLLIDES -- it puts K2 and K1 in the same byte --
1800 // and the lossless gate caught it on the first element. Spelled out per nibble instead.
1801 unsigned v = (nib16 & 0x000Fu)
1802 | ((nib16 & 0x00F0u) << 4)
1803 | ((nib16 & 0x0F00u) << 8)
1804 | ((nib16 & 0xF000u) << 12);
1805 // __vsub4, NOT `v - 0x08080808`: a 32-bit subtract BORROWS ACROSS BYTE BOUNDARIES, so every
1806 // nibble below 8 silently decrements the next K-element. That is off-by-one on ~37% of values --
1807 // small enough to look like quantisation noise in an accuracy number, and caught here only because
1808 // the requantise->expand round trip was checked against the original int8 in isolation.
1809 return __vsub4(v, 0x08080808u);
1810}
1811template<bool TWOPASS, bool TILEDW, bool Q4W>
1812__device__ __forceinline__ void bq_gemm_tc_body(
1813 float* __restrict__ Y, const signed char* __restrict__ wq, const __half* __restrict__ ws,
1814 const signed char* __restrict__ Xh, const signed char* __restrict__ Xl,
1815 const float* __restrict__ Xs, int K, int N, int B, int YS,
1816 float* __restrict__ P, int SPLIT){
1817 const int NCT=TC_TJ/8, NRT=8/NCT, TI=16*NRT*TC_WM, SW=TC_KC/4, NB=TC_KC/32;
1818 const int WV=TI*(TC_KC/16)/256;
1819 static_assert(TC_TJ%8==0 && (TC_TJ/8)*(8/(TC_TJ/8))==8,
1820 "TC_TJ must be 8, 16, 32 or 64 so the warp -> sub-tile map is a bijection");
1821 static_assert(TI*(TC_KC/16)>=256,
1822 "each thread must own at least one int4 of the staged weight tile");
1823 // ★ TC_KC MUST BE A WHOLE NUMBER OF Q8_0 BLOCKS.
1824 // `NB = TC_KC/32` is integer division and the mma loop runs `blk < NB`, so a TC_KC that is not a
1825 // multiple of 32 stages K-elements that are never multiplied: at TC_KC=48, NB=1, and 16 of every 48
1826 // staged elements are silently dropped. It does not crash and it does not diverge across sequences
1827 // -- `BATCH_AGREEMENT` reported a perfect 32/32 while every sequence was equally wrong -- it just
1828 // runs 1/3 faster and computes a different model. Measured at 760.29 tok/s against the honest 693.97,
1829 // which is exactly the shape of a result worth refusing.
1830 static_assert(TC_KC % 32 == 0,
1831 "TC_KC must be a multiple of 32: NB=TC_KC/32 truncates, and the K-elements past NB*32 "
1832 "are staged but never multiplied -- a silently wrong, silently faster kernel");
1833#if !TC_DYNSMEM
1834 static_assert((size_t)TC_STAGES*((size_t)TI*(TC_KC/4)*4 + (size_t)(TWOPASS?2:1)*TC_TJ*(TC_KC/4)*4
1835 + (size_t)TC_TJ*(TC_KC/32)*4 + (size_t)TI*(TC_KC/32)*2) <= 48*1024,
1836 "the staged tile must fit the 48 KB static shared-memory budget");
1837#endif
1838 // ── WHERE THE STAGED TILE LIVES ─────────────────────────────────────────────────────────────
1839 // Static `__shared__` is capped at 48 KB per block no matter what the chip has. H100 has 228 KB and
1840 // A100 164 KB, reachable only as DYNAMIC shared memory behind cudaFuncSetAttribute. That cap is not
1841 // a detail: it is what fixes TC_STAGES at 2, and staging depth is how many weight loads are in
1842 // flight, which is the whole reason the tile sits at 41.5% of the 3180 GB/s this chip can stream.
1843 // The layout is byte-identical either way -- same dimensions, same indices, so every `sW[st][r][c]`
1844 // below is untouched and the swizzle never learns which allocation it is in.
1845#if TC_DYNSMEM
1846 constexpr int TI_C = 16*(8/(TC_TJ/8))*TC_WM, SWC = TC_KC/4, NBC = TC_KC/32;
1847 constexpr int TJL = TWOPASS ? TC_TJ : 1;
1848 extern __shared__ __align__(16) char bq_tc_smem[];
1849 constexpr size_t OFF_H = (size_t)TC_STAGES*TI_C*SWC*sizeof(int);
1850 constexpr size_t OFF_L = OFF_H + (size_t)TC_STAGES*TC_TJ*SWC*sizeof(int);
1851 constexpr size_t OFF_WS = OFF_L + (size_t)TC_STAGES*TJL*SWC*sizeof(int);
1852 constexpr size_t OFF_XS = OFF_WS + (size_t)TC_STAGES*TI_C*NBC*sizeof(__half);
1853 auto sW = reinterpret_cast<int (*)[TI_C][SWC]>(bq_tc_smem);
1854 auto sH = reinterpret_cast<int (*)[TC_TJ][SWC]>(bq_tc_smem + OFF_H);
1855 auto sL = reinterpret_cast<int (*)[TJL][SWC]>(bq_tc_smem + OFF_L);
1856 auto sWs = reinterpret_cast<__half (*)[TI_C][NBC]>(bq_tc_smem + OFF_WS);
1857 auto sXs = reinterpret_cast<float (*)[TC_TJ][NBC]>(bq_tc_smem + OFF_XS);
1858#else
1859 __shared__ int sW[TC_STAGES][16*(8/(TC_TJ/8))*TC_WM][TC_KC/4], sH[TC_STAGES][TC_TJ][TC_KC/4];
1860 __shared__ int sP[Q4W?TC_STAGES:1][Q4W?(16*(8/(TC_TJ/8))*TC_WM):1][Q4W?(TC_KC/8):1];
1861 __shared__ int sL[TC_STAGES][TWOPASS?TC_TJ:1][TC_KC/4];
1862 __shared__ __half sWs[TC_STAGES][16*(8/(TC_TJ/8))*TC_WM][TC_KC/32];
1863 __shared__ float sXs[TC_STAGES][TC_TJ][TC_KC/32];
1864#endif
1865 // ── THE M-TILE ORIGIN ───────────────────────────────────────────────────────────────────────
1866 // The tile used to require B <= TC_TJ <= 64, which caps the row dimension of a single GEMM at 64.
1867 // That is fine for decode (one token per sequence) and it is the wall in front of chunked PREFILL:
1868 // the whole point of processing T prompt tokens at once is to read the 24.33 GB of weights ONCE for
1869 // B*T rows instead of once per token, and B*T is immediately larger than 64. grid.z now walks the
1870 // rows in TC_TJ-wide tiles; at grid.z==1 the origin is 0 and every index is what it was.
1871 const int b0=blockIdx.z*TC_TJ;
1872 const int tid=threadIdx.x, i0=blockIdx.x*TI, KB=K/32;
1873 const int wid=tid>>5, lane=tid&31, gid=lane>>2, tig=lane&3;
1874 const int ct=wid%NCT, rt=wid/NCT;
1875 const int kchunks=(K+TC_KC-1)/TC_KC, nchunk=(kchunks+SPLIT-1)/SPLIT;
1876 const int kbeg=blockIdx.y*nchunk*TC_KC, kend=min(K,kbeg+nchunk*TC_KC);
1877 float acc[TC_WM][4];
1878 #pragma unroll
1879 for(int m=0;m<TC_WM;m++){
1880 #pragma unroll
1881 for(int i=0;i<4;i++) acc[m][i]=0.f; }
1882 // Stage n+1's global reads issue into registers before stage n's mma, so the HBM latency overlaps
1883 // the arithmetic instead of serialising behind the barrier. Worth 1.27x on the deep-K shape.
1884 // Only the WEIGHTS are prefetched. They are 24.33 GB of the traffic in a step against roughly
1885 // 0.1 GB of activations, and the activation slab is small enough to sit in L2 across the whole
1886 // step, so staging it through registers buys nothing — while its register array would force
1887 // TC_TJ*(TC_KC/16) >= 256 and exclude TC_TJ=16 from the tile entirely.
1888 int4 pw[TI*(TC_KC/16)/256];
1889 // ★ Q4 IS NOT IMPLEMENTED ON THIS PATH, AND SILENCE HERE IS THE DANGEROUS KIND.
1890 // `BQTC_GLOAD` below stages the weight tile through registers and is the `#else` of the cp.async
1891 // pipeline, i.e. what runs at TC_STAGES==1 -- which is exactly what the decode shape bakes. It reads
1892 // at the INT8 stride `o*K` with no Q4 branch, so against a Q4 buffer of N*K/2 it runs ~44 MB past the
1893 // end, and its consumer writes into `sW` while the Q4 fragment reader reads `sP`, which nothing then
1894 // fills. So the combination does not merely fault: where it survives the fault it computes from an
1895 // uninitialised tile. A Q4 number taken at TC_STAGES==1 measures neither Q4 nor this model.
1896 // Q4 rides the cp.async path, where BQTC_ISSUE carries the packed stride and the packed tile.
1897 static_assert(!Q4W || TC_STAGES > 1,
1898 "Q4 requires the cp.async staging path (TC_STAGES > 1): the register-prefetch path reads the "
1899 "packed weights at the int8 stride and fills the wrong shared tile");
1900 #define BQTC_GLOAD(k0) { \
1901 _Pragma("unroll") for(int u=0;u<WV;u++){ int idx=tid+u*256, r=idx/(TC_KC/16), c=idx%(TC_KC/16), o=i0+r; \
1902 pw[u]=(o<N&&(k0)+c*16<K)?((const int4*)(wq+(long)o*K+(k0)))[c]:make_int4(0,0,0,0); } }
1903 // Stage `st` is filled for the chunk beginning at `kk`. Weights and activations go asynchronously;
1904 // the two scale arrays are 768 bytes a stage and stay synchronous, which is safe because they are
1905 // written into a DIFFERENT stage from the one being read and the trailing __syncthreads separates
1906 // the write from the previous reader of that buffer.
1907 #define BQTC_ISSUE(kk, st) { \
1908 _Pragma("unroll") for(int u=0;u<WV;u++){ int idx=tid+u*256, r=idx/(TC_KC/16), c=idx%(TC_KC/16), o=i0+r; \
1909 const bool ok=(o<N&&(kk)+c*16<K); \
1910 if(Q4W){ /* PACKED straight through cp.async -- expansion happens at the fragment load, so \
1911 the async pipeline survives and the shared tile is half size too */ \
1912 if(c < (TC_KC/32)){ \
1913 const bool ok4=(o<N&&(kk)+c*32<K); \
1914 const int4* g4=((const int4*)(wq+(long)(ok4?o:0)*(K>>1)+((kk)>>1)))+c; \
1915 bq_cpa16(&sP[Q4W?st:0][Q4W?r:0][BQTC_SWZC4(r,c)*4], g4, ok4); } \
1916 } else { \
1917 const int4* g = TILEDW \
1918 ? ((const int4*)(wq + ((long)(i0/TI)*((K+TC_KC-1)/TC_KC) + (long)((kk)/TC_KC))*((long)TI*TC_KC) \
1919 + (long)r*TC_KC)) + c \
1920 : ((const int4*)(wq+(long)(ok?o:0)*K+(kk)))+c; \
1921 bq_cpa16(&sW[st][r][BQTC_SWZC(r,c)*4], g, TILEDW ? (o<N) : ok); } } \
1922 for(int idx=tid; idx<TC_TJ*(TC_KC/16); idx+=256){ int r=idx/(TC_KC/16), c=idx%(TC_KC/16); \
1923 const bool ok=(r+b0<B&&(kk)+c*16<K); const int4* g=((const int4*)(Xh+(long)(ok?r+b0:0)*K+(kk)))+c; \
1924 bq_cpa16(&sH[st][r][BQTC_SWZC(r,c)*4], g, ok); \
1925 if(TWOPASS){ const int rl=TWOPASS?r:0; const int4* gl=((const int4*)(Xl+(long)(ok?r+b0:0)*K+(kk)))+c; \
1926 bq_cpa16(&sL[st][rl][BQTC_SWZC(rl,c)*4], gl, ok); } } \
1927 for(int idx=tid; idx<TI*NB; idx+=256){ int r=idx/NB,c=idx%NB,o=i0+r; \
1928 sWs[st][r][c] = TILEDW \
1929 ? ((o<N) ? ws[((long)(i0/TI)*((K+TC_KC-1)/TC_KC) + (long)((kk)/TC_KC))*((long)TI*NB) \
1930 + (long)r*NB + c] : __float2half(0.f)) \
1931 : ((o<N&&(kk)/32+c<KB)?ws[(long)o*KB+(kk)/32+c]:__float2half(0.f)); } \
1932 for(int idx=tid; idx<TC_TJ*NB; idx+=256){ int r=idx/NB,c=idx%NB; \
1933 sXs[st][r][c]=(r+b0<B&&(kk)/32+c<KB)?Xs[(long)(r+b0)*KB+(kk)/32+c]:0.f; } \
1934 asm volatile("cp.async.commit_group;\n" ::); }
1935#if TC_STAGES > 1
1936 // An S-stage pipeline has to be PRIMED with S-1 chunks before the loop, because the loop's
1937 // `wait_group(S-1)` only guarantees the oldest of S-1 in-flight groups has landed. Priming one
1938 // chunk is correct at S=2 and SILENTLY RACY above it: measured at S=3, B=32 as 565.73 tok/s with
1939 // BATCH_AGREEMENT 0/32 -- a plausible number from a kernel reading buffers that had not arrived.
1940 { int kp=kbeg; for(int st=0; st<TC_STAGES-1 && kp<kend; ++st, kp+=TC_KC){ BQTC_ISSUE(kp, st) } }
1941#else
1942 if(kbeg<kend){ BQTC_GLOAD(kbeg) }
1943#endif
1944 int chunk=0;
1945 for(int k0=kbeg;k0<kend;k0+=TC_KC,++chunk){
1946 const int sg = (TC_STAGES==1) ? 0 : (chunk % TC_STAGES);
1947#if TC_STAGES > 1
1948 // Issue the NEXT chunk before waiting on this one, so the copy is in flight across the mma.
1949 // ★ THE WAIT DEPTH DEPENDS ON WHETHER WE JUST ISSUED, and getting this wrong is a RACE that runs
1950 // FASTER than the correct kernel. Groups retire in order, so while a next chunk is in flight,
1951 // leaving TC_STAGES-1 pending means "this chunk has landed". On the LAST chunk there is no next
1952 // issue, so the only pending group is the one filling the buffer this chunk is about to read --
1953 // and wait_group(TC_STAGES-1) permits exactly that one to still be in flight. It must be 0 there.
1954 // Measured cost of the bug: 565.87 tok/s (the fastest number in the A/B) with first token 29858
1955 // instead of 8294 and BATCH_AGREEMENT 0/32. The equivalence lock PASSED it -- a race that usually
1956 // lands in time -- which is why the lock now sweeps chunk counts instead of one shape.
1957 // The chunk to fetch is the one S-1 ahead -- the slot this iteration is about to free -- not the
1958 // next one. `k0+TC_KC` is the same expression only when S=2.
1959 BQ_CLK(c0)
1960 if(k0+(TC_STAGES-1)*TC_KC<kend){
1961 BQTC_ISSUE(k0+(TC_STAGES-1)*TC_KC, (chunk+TC_STAGES-1)%TC_STAGES)
1962 asm volatile("cp.async.wait_group %0;\n" :: "n"(TC_STAGES-1));
1963 } else {
1964 // Nothing left to issue, so the pending groups ARE the ones still being consumed: drain fully.
1965 asm volatile("cp.async.wait_group 0;\n" ::);
1966 }
1967 __syncthreads();
1968 if(Q4W){
1969 // packed int c of a row covers K-elements [8c,8c+8) -> expanded ints 2c and 2c+1. Done ONCE
1970 // per staged element here instead of once per USE: the four column-warps all read the same A
1971 // rows, so expanding at the fragment load repeated this 4x and measured 0.119 ms against 0.085.
1972 for(int idx=tid; idx<TI*(TC_KC/8); idx+=256){
1973 const int r=idx/(TC_KC/8), c=idx%(TC_KC/8);
1974 const unsigned pk=(unsigned)sP[Q4W?sg:0][Q4W?r:0][BQTC_SWZJ4(r,c)];
1975 sW[sg][r][BQTC_SWZJ(r,2*c)] = (int)bq_q4_expand(pk & 0xFFFFu);
1976 sW[sg][r][BQTC_SWZJ(r,2*c+1)] = (int)bq_q4_expand(pk >> 16);
1977 }
1978 __syncthreads();
1979 }
1980 BQ_CLK(c2) BQ_ACC(1,c0,c2)
1981#else
1982 #pragma unroll
1983 for(int u=0;u<WV;u++){ int idx=tid+u*256, r=idx/(TC_KC/16), c=idx%(TC_KC/16);
1984 const int pc=BQTC_SWZC(r,c)*4;
1985 sW[sg][r][pc]=pw[u].x; sW[sg][r][pc+1]=pw[u].y; sW[sg][r][pc+2]=pw[u].z; sW[sg][r][pc+3]=pw[u].w; }
1986 for(int idx=tid; idx<TC_TJ*(TC_KC/16); idx+=256){ int r=idx/(TC_KC/16), c=idx%(TC_KC/16);
1987 bool ok=(r+b0<B&&k0+c*16<K);
1988 const int4 h=ok?((const int4*)(Xh+(long)(r+b0)*K+k0))[c]:make_int4(0,0,0,0);
1989 const int pch=BQTC_SWZC(r,c)*4;
1990 sH[sg][r][pch]=h.x; sH[sg][r][pch+1]=h.y; sH[sg][r][pch+2]=h.z; sH[sg][r][pch+3]=h.w;
1991 if(TWOPASS){ const int4 l=ok?((const int4*)(Xl+(long)(r+b0)*K+k0))[c]:make_int4(0,0,0,0);
1992 const int rl=TWOPASS?r:0;
1993 const int pcl=BQTC_SWZC(rl,c)*4;
1994 sL[sg][rl][pcl]=l.x; sL[sg][rl][pcl+1]=l.y; sL[sg][rl][pcl+2]=l.z; sL[sg][rl][pcl+3]=l.w; } }
1995 for(int idx=tid; idx<TI*NB; idx+=256){ int r=idx/NB,c=idx%NB,o=i0+r;
1996 sWs[sg][r][c]=(o<N&&k0/32+c<KB)?ws[(long)o*KB+k0/32+c]:__float2half(0.f); }
1997 for(int idx=tid; idx<TC_TJ*NB; idx+=256){ int r=idx/NB,c=idx%NB;
1998 sXs[sg][r][c]=(r+b0<B&&k0/32+c<KB)?Xs[(long)(r+b0)*KB+k0/32+c]:0.f; }
1999 __syncthreads();
2000 if(k0+TC_KC<kend){ BQTC_GLOAD(k0+TC_KC) }
2001#endif
2002 #pragma unroll
2003 for(int blk=0; blk<NB; ++blk){
2004 const int jb=ct*8+gid;
2005 const int b0h=sH[sg][jb][BQTC_SWZJ(jb,blk*8+tig)], b1h=sH[sg][jb][BQTC_SWZJ(jb,blk*8+tig+4)];
2006 const int jbl=TWOPASS?jb:0;
2007 const int b0l=sL[sg][jbl][BQTC_SWZJ(jbl,blk*8+tig)], b1l=sL[sg][jbl][BQTC_SWZJ(jbl,blk*8+tig+4)];
2008 const float x0=sXs[sg][ct*8+tig*2][blk], x1=sXs[sg][ct*8+tig*2+1][blk];
2009 #pragma unroll
2010 for(int m=0;m<TC_WM;m++){
2011 const int rb=(rt*TC_WM+m)*16;
2012 const int a0=sW[sg][rb+gid][BQTC_SWZJ(rb+gid,blk*8+tig)],
2013 a1=sW[sg][rb+gid+8][BQTC_SWZJ(rb+gid+8,blk*8+tig)];
2014 const int a2=sW[sg][rb+gid][BQTC_SWZJ(rb+gid,blk*8+tig+4)],
2015 a3=sW[sg][rb+gid+8][BQTC_SWZJ(rb+gid+8,blk*8+tig+4)];
2016 int dh[4]={0,0,0,0}, dl[4]={0,0,0,0};
2017 bq_mma_s8(dh,a0,a1,a2,a3,b0h,b1h);
2018 if(TWOPASS) bq_mma_s8(dl,a0,a1,a2,a3,b0l,b1l);
2019 const float w0=__half2float(sWs[sg][rb+gid][blk]), w1=__half2float(sWs[sg][rb+gid+8][blk]);
2020 if(TWOPASS){
2021 acc[m][0]+=((float)dh[0]+(float)dl[0]*(1.f/127.f))*(w0*x0);
2022 acc[m][1]+=((float)dh[1]+(float)dl[1]*(1.f/127.f))*(w0*x1);
2023 acc[m][2]+=((float)dh[2]+(float)dl[2]*(1.f/127.f))*(w1*x0);
2024 acc[m][3]+=((float)dh[3]+(float)dl[3]*(1.f/127.f))*(w1*x1);
2025 } else {
2026 acc[m][0]+=(float)dh[0]*(w0*x0); acc[m][1]+=(float)dh[1]*(w0*x1);
2027 acc[m][2]+=(float)dh[2]*(w1*x0); acc[m][3]+=(float)dh[3]*(w1*x1);
2028 }
2029 }
2030 }
2031 __syncthreads();
2032 BQ_CLK(c3) BQ_ACC(2,c2,c3)
2033 }
2034 float* out = (SPLIT==1) ? Y : (P + (long)blockIdx.y*N*B);
2035 const int ys = (SPLIT==1) ? YS : N;
2036 #pragma unroll
2037 for(int m=0;m<TC_WM;m++){
2038 const int rb=(rt*TC_WM+m)*16, o0=i0+rb+gid, o1=o0+8;
2039 const int j0=ct*8+tig*2+b0, j1=j0+1;
2040 if(o0<N){ if(j0<B) out[(long)j0*ys+o0]=acc[m][0]; if(j1<B) out[(long)j1*ys+o0]=acc[m][1]; }
2041 if(o1<N){ if(j0<B) out[(long)j0*ys+o1]=acc[m][2]; if(j1<B) out[(long)j1*ys+o1]=acc[m][3]; }
2042 }
2043 #undef BQTC_GLOAD
2044 #undef BQTC_ISSUE
2045}
2046// Which wrappers to emit. A program launches exactly ONE of these, but both are compiled, and the
2047// two-pass body stages 2x the activation -- so leaving it in prices the 48 KB budget against the lane
2048// that is NOT running and caps TC_KC for the one that is. Measured on an H100 at the decode shape,
2049// that cap is worth 23%: KC=64 -> 1077 GB/s, KC=128 -> 1320 GB/s (bq_sk_reduce included, split 16).
2050// Bit 0 = single-pass, bit 1 = two-pass. The default emits both, which is what the equivalence locks
2051// compile and compare.
2052#ifndef BQ_TC_LANES
2053#define BQ_TC_LANES 3
2054#endif
2055#if BQ_TC_LANES & 2
2056__global__ void __launch_bounds__(256, TC_MINBLOCKS) bq_gemm_tc(
2057 float* __restrict__ Y, const signed char* __restrict__ wq, const __half* __restrict__ ws,
2058 const signed char* __restrict__ Xh, const signed char* __restrict__ Xl,
2059 const float* __restrict__ Xs, int K, int N, int B, int YS, float* __restrict__ P, int SPLIT){
2060 bq_gemm_tc_body<true,false,false>(Y,wq,ws,Xh,Xl,Xs,K,N,B,YS,P,SPLIT);
2061}
2062#endif // BQ_TC_LANES & 2 -- bq_gemm_tcm below is the PREFILL tile and belongs to neither lane
2063// ══ ROWS IN M — the operand roles swapped ═══════════════════════════════════════════════════════
2064// `bq_gemm_tc*` puts the WEIGHTS in the m16 operand and the SEQUENCES in n8, so a block covers at most
2065// 8 sequences per column group and, with 8 warps, TC_TJ <= 64 rows. Rows per weight pass are capped at
2066// 64 BY THE INSTRUCTION MAPPING. That is invisible at decode (B<=64 needs one pass anyway) and it is
2067// the whole prefill deficit: 512 prompt tokens cost 512 passes over 24.33 GB.
2068//
2069// CUTLASS -- and so vLLM's fp8 path -- does the reverse: rows go in the WIDE dimension
2070// (`wgmma.mma_async.m64nNk32`, tile covering 128-256 rows) and output channels in N. This is that
2071// mapping expressed with the Ampere instruction we already have: ACTIVATIONS in m16, WEIGHTS in n8.
2072// A block then covers MI sequence rows against MT_NT channels, so the weight tile [MT_NT x K] is read
2073// once per MI rows instead of once per TC_TJ. At MI=256 that is 8x fewer passes than TC_TJ=32.
2074//
2075// It is the mirror image of the other tile's trade, so it is a LANE, not a replacement: at decode
2076// (rows=32) a 256-row tile would be 7/8 masked, which is exactly why the shape rule keeps
2077// `bq_gemm_tc1` there. Correctness first here -- padded rows, synchronous staging, no cp.async.
2078#ifndef MT_WM
2079#define MT_WM 4
2080#endif
2081// MT_NT=32 / MI=128, chosen by the traffic+occupancy model rather than by minimising weight passes
2082// alone. Total HBM traffic for a chunk is K*N*(MI/MT_NT + rows/MI): the first term is the activation
2083// tile re-read by each of the N/MT_NT blocks, the second the weight passes. Minimising only the
2084// second gave MT_NT=16 / MI=256, which measured 22.8% when the pass count said ~7x -- because it
2085// carried 1.7x the total traffic AND needed 39 KB of shared for two stages, so only ONE block fit per
2086// SM. At MT_NT=32: 12 K*N against 20, and 22.8 KB -> two blocks resident.
2087// MT_NT=64 / MT_WN=2 / MI=128, swept on an A100-SXM4 over a 512-token prefill, token stream identical
2088// at every point: MT_NT=32/WN=2 21401 ms | MT_NT=64/WN=4 21184 | **MT_NT=64/WN=2 19571** | one-stage
2089// 25728. More accumulators (WN=4, 165 regs, 45.6 KB) did NOT pay -- the winner keeps 105 regs and
2090// 27 KB, i.e. two blocks per SM, while still getting 2.5 operand loads per mma instead of 4.5.
2091#ifndef MT_NT
2092#define MT_NT 64
2093#endif
2094#ifndef MT_KC
2095#define MT_KC 64
2096#endif
2097// Ported from `bq_gemm_tc1`, which measured +25.9% for the swizzle at KC=128 and +9.6%/+17.4% for
2098// 16-byte cp.async. Same reasoning: the `+1` row padding costs 16-byte alignment (which cp.async
2099// requires) and does not even buy conflict-freedom; an XOR swizzle gives both. Mask is (chunks-1) --
2100// a row holds MT_KC/16 chunks and XOR-ing past that addresses the NEXT row.
2101// Column sub-tiles per warp. Each A fragment (4 shared reads) then feeds MT_WN mmas instead of one,
2102// which is the ratio that decides a COMPUTE-bound GEMM: prefill runs at 6.1% of int8 peak against
2103// vLLM's 75.7%, and at MT_WN=1 the tile issues 4.5 shared-memory operand loads per mma. CUTLASS gets
2104// under 1 by reusing fragments across a register tile. MT_WN=2 halves it to 2.5 at 32 accumulators.
2105#ifndef MT_WN
2106#define MT_WN 2
2107#endif
2108#ifndef MT_STAGES
2109#define MT_STAGES 2
2110#endif
2111#define MTM_CH (MT_KC/16)
2112#define MTM_SWZC(r, c) ((c) ^ ((r) & (MTM_CH-1)))
2113#define MTM_SWZJ(r, j) ((MTM_SWZC((r), (j)>>2) << 2) | ((j) & 3))
2114__global__ void __launch_bounds__(256) bq_gemm_tcm(
2115 float* __restrict__ Y, const signed char* __restrict__ wq, const __half* __restrict__ ws,
2116 const signed char* __restrict__ Xh, const float* __restrict__ Xs,
2117 int K, int N, int B, int YS){
2118 const int NCT=MT_NT/(8*MT_WN), NRT=8/NCT, MI=16*NRT*MT_WM, NB=MT_KC/32;
2119 static_assert(MT_NT%(8*MT_WN)==0 && MT_NT/(8*MT_WN)<=8,
2120 "MT_NT must be 8*MT_WN*NCT with NCT<=8 so the warp->channel map is a bijection");
2121 static_assert(MT_KC%32==0, "MT_KC must be a whole number of Q8_0 blocks");
2122 #define MTM_MI (16*(8/(MT_NT/(8*MT_WN)))*MT_WM)
2123 __shared__ int sA[MT_STAGES][MTM_MI][MT_KC/4];
2124 __shared__ int sB[MT_STAGES][MT_NT][MT_KC/4];
2125 __shared__ float sAs[MT_STAGES][MTM_MI][MT_KC/32];
2126 __shared__ __half sBs[MT_STAGES][MT_NT][MT_KC/32];
2127 static_assert((size_t)MT_STAGES*((size_t)MI*(MT_KC/4)*4 + (size_t)MT_NT*(MT_KC/4)*4
2128 + (size_t)MI*(MT_KC/32)*4 + (size_t)MT_NT*(MT_KC/32)*2) <= 48*1024,
2129 "the staged rows-in-M tile must fit the 48 KB static shared-memory budget");
2130 const int tid=threadIdx.x, lane=tid&31, wid=tid>>5, gid=lane>>2, tig=lane&3;
2131 const int ct=wid%NCT, rt=wid/NCT;
2132 const int m0=blockIdx.y*MI, n0=blockIdx.x*MT_NT, KB=K/32;
2133 float acc[MT_WM][MT_WN][4];
2134 #pragma unroll
2135 for(int m=0;m<MT_WM;m++){
2136 #pragma unroll
2137 for(int n=0;n<MT_WN;n++){
2138 #pragma unroll
2139 for(int i=0;i<4;i++) acc[m][n][i]=0.f; } }
2140 #define MTM_ISSUE(kk, st) { \
2141 for(int idx=tid; idx<MI*MTM_CH; idx+=256){ int r=idx/MTM_CH, c=idx%MTM_CH; \
2142 const bool ok=(m0+r<B && (kk)+c*16<K); const int sr=ok?(m0+r):0; \
2143 bq_cpa16(&sA[st][r][MTM_SWZC(r,c)*4], ((const int4*)(Xh+(long)sr*K+(kk)))+c, ok); } \
2144 for(int idx=tid; idx<MT_NT*MTM_CH; idx+=256){ int r=idx/MTM_CH, c=idx%MTM_CH; \
2145 const bool ok=(n0+r<N && (kk)+c*16<K); const int sr=ok?(n0+r):0; \
2146 bq_cpa16(&sB[st][r][MTM_SWZC(r,c)*4], ((const int4*)(wq+(long)sr*K+(kk)))+c, ok); } \
2147 for(int idx=tid; idx<MI*NB; idx+=256){ int r=idx/NB,c=idx%NB; \
2148 sAs[st][r][c]=(m0+r<B&&(kk)/32+c<KB)?Xs[(long)(m0+r)*KB+(kk)/32+c]:0.f; } \
2149 for(int idx=tid; idx<MT_NT*NB; idx+=256){ int r=idx/NB,c=idx%NB; \
2150 sBs[st][r][c]=(n0+r<N&&(kk)/32+c<KB)?ws[(long)(n0+r)*KB+(kk)/32+c]:__float2half(0.f); } \
2151 asm volatile("cp.async.commit_group;\n" ::); }
2152 { int kp=0; for(int st=0; st<MT_STAGES-1 && kp<K; ++st, kp+=MT_KC){ MTM_ISSUE(kp, st) } }
2153 int chunk=0;
2154 for(int k0=0;k0<K;k0+=MT_KC,++chunk){
2155 const int sg=(MT_STAGES==1)?0:(chunk%MT_STAGES);
2156 if(k0+(MT_STAGES-1)*MT_KC<K){
2157 MTM_ISSUE(k0+(MT_STAGES-1)*MT_KC, (chunk+MT_STAGES-1)%MT_STAGES)
2158 asm volatile("cp.async.wait_group %0;\n" :: "n"(MT_STAGES-1));
2159 } else {
2160 asm volatile("cp.async.wait_group 0;\n" ::);
2161 }
2162 __syncthreads();
2163 #pragma unroll
2164 for(int blk=0; blk<NB; ++blk){
2165 // B fragments for every column sub-tile, hoisted so each A fragment below feeds all of them.
2166 int b0[MT_WN], b1[MT_WN]; float w0[MT_WN], w1[MT_WN];
2167 #pragma unroll
2168 for(int n=0;n<MT_WN;n++){
2169 const int nb=(ct*MT_WN+n)*8+gid;
2170 b0[n]=sB[sg][nb][MTM_SWZJ(nb,blk*8+tig)]; b1[n]=sB[sg][nb][MTM_SWZJ(nb,blk*8+tig+4)];
2171 w0[n]=__half2float(sBs[sg][(ct*MT_WN+n)*8+tig*2][blk]);
2172 w1[n]=__half2float(sBs[sg][(ct*MT_WN+n)*8+tig*2+1][blk]);
2173 }
2174 #pragma unroll
2175 for(int m=0;m<MT_WM;m++){
2176 const int mb=(rt*MT_WM+m)*16;
2177 const int a0=sA[sg][mb+gid][MTM_SWZJ(mb+gid,blk*8+tig)], a1=sA[sg][mb+gid+8][MTM_SWZJ(mb+gid+8,blk*8+tig)];
2178 const int a2=sA[sg][mb+gid][MTM_SWZJ(mb+gid,blk*8+tig+4)], a3=sA[sg][mb+gid+8][MTM_SWZJ(mb+gid+8,blk*8+tig+4)];
2179 const float x0=sAs[sg][mb+gid][blk], x1=sAs[sg][mb+gid+8][blk];
2180 #pragma unroll
2181 for(int n=0;n<MT_WN;n++){
2182 int d[4]={0,0,0,0};
2183 bq_mma_s8(d,a0,a1,a2,a3,b0[n],b1[n]);
2184 acc[m][n][0]+=(float)d[0]*(x0*w0[n]); acc[m][n][1]+=(float)d[1]*(x0*w1[n]);
2185 acc[m][n][2]+=(float)d[2]*(x1*w0[n]); acc[m][n][3]+=(float)d[3]*(x1*w1[n]);
2186 }
2187 }
2188 }
2189 __syncthreads();
2190 }
2191 #undef MTM_ISSUE
2192 #pragma unroll
2193 for(int m=0;m<MT_WM;m++){
2194 const int mb=(rt*MT_WM+m)*16, s0=m0+mb+gid, s1=s0+8;
2195 #pragma unroll
2196 for(int n=0;n<MT_WN;n++){
2197 const int c0=n0+(ct*MT_WN+n)*8+tig*2, c1=c0+1;
2198 if(s0<B){ if(c0<N) Y[(long)s0*YS+c0]=acc[m][n][0]; if(c1<N) Y[(long)s0*YS+c1]=acc[m][n][1]; }
2199 if(s1<B){ if(c0<N) Y[(long)s1*YS+c0]=acc[m][n][2]; if(c1<N) Y[(long)s1*YS+c1]=acc[m][n][3]; }
2200 }
2201 }
2202}
2203#if BQ_TC_LANES & 1
2204__global__ void __launch_bounds__(256, TC_MINBLOCKS) bq_gemm_tc1(
2205 float* __restrict__ Y, const signed char* __restrict__ wq, const __half* __restrict__ ws,
2206 const signed char* __restrict__ Xh,
2207 const float* __restrict__ Xs, int K, int N, int B, int YS, float* __restrict__ P, int SPLIT){
2208 bq_gemm_tc_body<false,false,false>(Y,wq,ws,Xh,(const signed char*)0,Xs,K,N,B,YS,P,SPLIT);
2209}
2210// The same tile fed from the PRE-TILED weight layout: one streaming read per chunk instead of TI
2211// separate row-streams. A permutation of the same bytes, so it must compute the identical thing.
2212__global__ void __launch_bounds__(256, TC_MINBLOCKS) bq_gemm_tc1_tiled(
2213 float* __restrict__ Y, const signed char* __restrict__ wq, const __half* __restrict__ ws,
2214 const signed char* __restrict__ Xh,
2215 const float* __restrict__ Xs, int K, int N, int B, int YS, float* __restrict__ P, int SPLIT){
2216 bq_gemm_tc_body<false,true,false>(Y,wq,ws,Xh,(const signed char*)0,Xs,K,N,B,YS,P,SPLIT);
2217}
2218// Activations to e4m3 with ONE scale per row, mirroring the weights' one-scale-per-channel. This is
2219// where the fp8 lane stops paying the Q8_0 tax: the product of the two scales is applied ONCE in the
2220// epilogue instead of once per 32 K-elements.
2221// One fp32 scale per 128 K-channels, and 448 is the largest finite e4m3 — the two constants the whole
2222// DeepGEMM lane is parameterised by. Defined before the host helpers that size shared memory with them.
2223#define BQ_DG_BLOCK_K 128
2224#define BQ_DG_E4M3_MAX 448.0f
2225
2226// ── HOST SIDE OF THE DEEPGEMM LANE ──────────────────────────────────────────────────────────────
2227// The kernel takes FIVE `__grid_constant__` TMA descriptors (A, B, scale-A, scale-B, C/D). DeepGEMM
2228// builds them from a `torch::Tensor`, using it for exactly `data_ptr()` and `element_size()` — so a
2229// raw pointer and a size replace it with nothing lost, and the lane needs no torch.
2230//
2231// Everything here is compiled only when the program is emitted for the DeepGEMM lane, because it
2232// needs `<deep_gemm/...>` + CUTLASS on the include path, `--expt-relaxed-constexpr`, CUDA >= 12.6 and
2233// `-lcuda`. The lane declares its own build requirements; nothing else in the program does.
2234#ifdef BQ_ENABLE_DG
2235#include <cuda.h>
2236#include <deep_gemm/impls/sm90_fp8_gemm_1d1d.cuh>
2237using namespace deep_gemm;
2238
2239#define BQ_DG_BM 64
2240#define BQ_DG_BN 128
2241// Shared-memory budget, from DeepGEMM's own `get_pipeline_config` rather than guessed: the C/D tile
2242// (1024-aligned) plus barriers, then as many stages as the rest of the budget allows.
2243#define BQ_DG_SMEM_EXTRA ((((BQ_DG_BM*BQ_DG_BN*4)+1023)/1024*1024) + 16*8*2)
2244#define BQ_DG_SMEM_STAGE (BQ_DG_BM*BQ_DG_BLOCK_K + BQ_DG_BN*BQ_DG_BLOCK_K \
2245 + ((BQ_DG_BM*4+127)/128*128) + ((BQ_DG_BN*4+127)/128*128))
2246#define BQ_DG_STAGES 7
2247#define BQ_DG_SMEM (BQ_DG_SMEM_EXTRA + BQ_DG_STAGES*BQ_DG_SMEM_STAGE)
2248
2249// TMA needs the inner dim to span whole 16-byte units, so the scale tensors pad their row count.
2250static inline int bq_dg_align(int mn, int elem) { const int u = 16 / elem; return (mn + u - 1) / u * u; }
2251
2252static CUtensorMap bq_dg_desc(void* p, CUtensorMapDataType dt, int elem,
2253 int gmem_inner, int gmem_outer, int smem_inner, int smem_outer,
2254 int gmem_outer_stride, int swizzle_bytes) {
2255 CUtensorMap m{};
2256 if (swizzle_bytes) smem_inner = swizzle_bytes / elem; // DeepGEMM does this too
2257 const cuuint64_t gd[2] = {(cuuint64_t)gmem_inner, (cuuint64_t)gmem_outer};
2258 const cuuint64_t gs[1] = {(cuuint64_t)((long)gmem_outer_stride * elem)};
2259 const cuuint32_t sd[2] = {(cuuint32_t)smem_inner, (cuuint32_t)smem_outer};
2260 const cuuint32_t es[2] = {1, 1};
2261 const CUtensorMapSwizzle sw = swizzle_bytes == 128 ? CU_TENSOR_MAP_SWIZZLE_128B
2262 : swizzle_bytes == 64 ? CU_TENSOR_MAP_SWIZZLE_64B
2263 : swizzle_bytes == 32 ? CU_TENSOR_MAP_SWIZZLE_32B
2264 : CU_TENSOR_MAP_SWIZZLE_NONE;
2265 const CUresult r = cuTensorMapEncodeTiled(&m, dt, 2, p, gd, gs, sd, es,
2266 CU_TENSOR_MAP_INTERLEAVE_NONE, sw, CU_TENSOR_MAP_L2_PROMOTION_L2_256B,
2267 CU_TENSOR_MAP_FLOAT_OOB_FILL_NONE);
2268 if (r != CUDA_SUCCESS) { const char* e = nullptr; cuGetErrorString(r, &e);
2269 fprintf(stderr, "DG descriptor failed: %s\n", e ? e : "?"); }
2270 return m;
2271}
2272// K-major operand (A or B) in fp8, swizzled 128B — the layout the kernel's TMA expects.
2273static CUtensorMap bq_dg_ab(void* p, int mn, int k, int block_mn) {
2274 return bq_dg_desc(p, CU_TENSOR_MAP_DATA_TYPE_UINT8, 1, k, mn, BQ_DG_BLOCK_K, block_mn, k, 128);
2275}
2276// Scale tensor, `[k_block][mn_aligned]`, never swizzled.
2277static CUtensorMap bq_dg_sf(void* p, int mn, int k, int block_mn) {
2278 const int a = bq_dg_align(mn, 4);
2279 return bq_dg_desc(p, CU_TENSOR_MAP_DATA_TYPE_FLOAT32, 4, a, k / BQ_DG_BLOCK_K, block_mn, 1, a, 0);
2280}
2281// The output stride is NOT always N: `ffbuf` is packed (stride == 2*dff) but `pbuf` carries the widest
2282// projection's row pitch, so an attention layer writes N=14336 columns into a 16384-wide row. Passing
2283// N as the stride there would shear the output by 2048 floats per row — a plausible-looking tensor of
2284// the right shape and entirely wrong.
2285static CUtensorMap bq_dg_cd(void* p, int m, int n, int row_stride) {
2286 return bq_dg_desc(p, CU_TENSOR_MAP_DATA_TYPE_FLOAT32, 4, n, m, BQ_DG_BN, BQ_DG_BM, row_stride, 0);
2287}
2288#endif // BQ_ENABLE_DG
2289
2290// ── CONVERTERS FOR THE DEEPGEMM LANE ────────────────────────────────────────────────────────────
2291// DeepGEMM's sm90 1d1d kernel wants fp8 e4m3 with ONE fp32 scale per 128 K-channels, and the scale
2292// tensors laid out `[K/128][MN_aligned]` — MN padded so the TMA inner dim is a multiple of 16 bytes.
2293// That is a different layout from the per-CHANNEL scales the older fp8 lane above produces, which is
2294// why these are separate kernels rather than a flag on those.
2295//
2296// Validated end to end on an H100: a host-side version of exactly this arithmetic fed DeepGEMM and the
2297// result matched a DOUBLE-precision reference at rel_L2 1.072e-4 — fp32-accumulation agreement. A
2298// wrong scale layout does not land at 1e-4, so that measurement covers the layout as well as the maths.
2299// Weights: Q8_0 (int8 + one __half scale per 32) -> fp8 e4m3 + one fp32 scale per 128 channels.
2300// One block per (row, k-block); the 128 values of a block are four Q8_0 sub-blocks.
2301__global__ void bq_q8_to_fp8_dg(unsigned char* __restrict__ out, float* __restrict__ sf,
2302 const signed char* __restrict__ q, const __half* __restrict__ s,
2303 int N, int K, int n_aligned){
2304 const int kb = K / BQ_DG_BLOCK_K;
2305 for (long idx = (long)blockIdx.x * blockDim.x + threadIdx.x; idx < (long)N * kb;
2306 idx += (long)gridDim.x * blockDim.x) {
2307 const int n = (int)(idx / kb), b = (int)(idx % kb);
2308 const long base = (long)n * K + (long)b * BQ_DG_BLOCK_K;
2309 float amax = 0.f;
2310 for (int i = 0; i < BQ_DG_BLOCK_K; ++i) {
2311 const long g = base + i;
2312 const float v = (float)q[g] * __half2float(s[(long)n * (K / 32) + (int)(g % K) / 32]);
2313 amax = fmaxf(amax, fabsf(v));
2314 }
2315 const float S = amax > 0.f ? amax / BQ_DG_E4M3_MAX : 1.f;
2316 sf[(long)b * n_aligned + n] = S; // [k_block][row], row-major in rows
2317 for (int i = 0; i < BQ_DG_BLOCK_K; ++i) {
2318 const long g = base + i;
2319 const float v = (float)q[g] * __half2float(s[(long)n * (K / 32) + (int)(g % K) / 32]);
2320#if BQ_FP8_TOOLKIT
2321 out[g] = (unsigned char)__nv_cvt_float_to_fp8(v / S, __NV_SATFINITE, __NV_E4M3);
2322#else
2323 out[g] = 0;
2324#endif
2325 }
2326 }
2327}
2328
2329// ★ THE TWO-PASS fp8 WEIGHT. e4m3 carries 3 mantissa bits, so one fp8 copy of a weight is ~2.3x
2330// lossier than the Q8_0 it replaces -- which is the whole reason the routed lane's token stream walks
2331// away from ours. The fix is the trick the int8 lane already uses on activations: represent what the
2332// first pass THREW AWAY in a second tensor of the same type. W ~= Whi + Wlo, both fp8, and the product
2333// becomes X*Whi' + X*Wlo' -- two GEMMs, ~7 effective mantissa bits, and an accuracy that lands back on
2334// top of int8 (measured against f64: int8 1.27e-2, fp8 x1 2.88e-2, fp8 hi/lo 1.35e-2).
2335//
2336// The 128 dequantised values are recomputed rather than held: a `float v[128]` is 512 B per thread and
2337// would spill. This runs ONCE at load time, so three passes over the block cost nothing that matters.
2338__global__ void bq_q8_to_fp8_dg2(unsigned char* __restrict__ hi, float* __restrict__ hisf,
2339 unsigned char* __restrict__ lo, float* __restrict__ losf,
2340 const signed char* __restrict__ q, const __half* __restrict__ sc,
2341 int rows, int K, int rows_aligned){
2342 const int kb = K / BQ_DG_BLOCK_K;
2343 for (long idx = (long)blockIdx.x * blockDim.x + threadIdx.x; idx < (long)rows * kb;
2344 idx += (long)gridDim.x * blockDim.x) {
2345 const int r = (int)(idx / kb), b = (int)(idx % kb);
2346 const long base = (long)r * K + (long)b * BQ_DG_BLOCK_K;
2347#define BQ_DQ(i) ((float)q[base + (i)] * __half2float(sc[(long)r * (K / 32) + (int)((b * BQ_DG_BLOCK_K + (i)) / 32)]))
2348 float amax = 0.f;
2349 for (int i = 0; i < BQ_DG_BLOCK_K; ++i) amax = fmaxf(amax, fabsf(BQ_DQ(i)));
2350 const float Sh = amax > 0.f ? amax / BQ_DG_E4M3_MAX : 1.f;
2351 hisf[(long)b * rows_aligned + r] = Sh;
2352 float rmax = 0.f;
2353 for (int i = 0; i < BQ_DG_BLOCK_K; ++i) {
2354 const float v = BQ_DQ(i);
2355#if BQ_FP8_TOOLKIT
2356 const unsigned char h = (unsigned char)__nv_cvt_float_to_fp8(v / Sh, __NV_SATFINITE, __NV_E4M3);
2357 hi[base + i] = h;
2358 rmax = fmaxf(rmax, fabsf(v - __half2float(__half(__nv_cvt_fp8_to_halfraw(h, __NV_E4M3))) * Sh));
2359#else
2360 hi[base + i] = 0;
2361#endif
2362 }
2363 const float Sl = rmax > 0.f ? rmax / BQ_DG_E4M3_MAX : 1.f;
2364 losf[(long)b * rows_aligned + r] = Sl;
2365 for (int i = 0; i < BQ_DG_BLOCK_K; ++i) {
2366#if BQ_FP8_TOOLKIT
2367 const float v = BQ_DQ(i);
2368 const float res = v - __half2float(__half(__nv_cvt_fp8_to_halfraw(hi[base + i], __NV_E4M3))) * Sh;
2369 lo[base + i] = (unsigned char)__nv_cvt_float_to_fp8(res / Sl, __NV_SATFINITE, __NV_E4M3);
2370#else
2371 lo[base + i] = 0;
2372#endif
2373 }
2374#undef BQ_DQ
2375 }
2376}
2377
2378// The second pass lands in a packed scratch; fold it onto the first, which sits at the consumer's own
2379// row pitch. DeepGEMM's epilogue only ever TMA-stores D -- it never reads C -- so the sum has to happen
2380// here rather than inside the GEMM.
2381__global__ void bq_dg_add(float* __restrict__ Y, int ys, const float* __restrict__ P, int n, int rows){
2382 for (long i = (long)blockIdx.x * blockDim.x + threadIdx.x; i < (long)rows * n;
2383 i += (long)gridDim.x * blockDim.x)
2384 Y[(i / n) * (long)ys + (i % n)] += P[i];
2385}
2386
2387// Activations: f32 rows -> fp8 e4m3 + one fp32 scale per 128 channels, same layout.
2388__global__ void bq_quant_act_fp8_dg(unsigned char* __restrict__ out, float* __restrict__ sf,
2389 const float* __restrict__ X, int rows, int K, int rows_aligned){
2390 const int kb = K / BQ_DG_BLOCK_K;
2391 for (long idx = (long)blockIdx.x * blockDim.x + threadIdx.x; idx < (long)rows * kb;
2392 idx += (long)gridDim.x * blockDim.x) {
2393 const int r = (int)(idx / kb), b = (int)(idx % kb);
2394 const float* x = X + (long)r * K + (long)b * BQ_DG_BLOCK_K;
2395 float amax = 0.f;
2396 for (int i = 0; i < BQ_DG_BLOCK_K; ++i) amax = fmaxf(amax, fabsf(x[i]));
2397 const float S = amax > 0.f ? amax / BQ_DG_E4M3_MAX : 1.f;
2398 sf[(long)b * rows_aligned + r] = S;
2399 for (int i = 0; i < BQ_DG_BLOCK_K; ++i)
2400#if BQ_FP8_TOOLKIT
2401 out[(long)r * K + (long)b * BQ_DG_BLOCK_K + i] =
2402 (unsigned char)__nv_cvt_float_to_fp8(x[i] / S, __NV_SATFINITE, __NV_E4M3);
2403#else
2404 out[(long)r * K + (long)b * BQ_DG_BLOCK_K + i] = 0;
2405#endif
2406 }
2407}
2408
2409// ★ THE CROSS-LAYER FUSION. The unfused chain has SiLU write its f32 result to a dff-wide buffer
2410// and a second kernel read it straight back to quantise. At B*T=1024 and dff=17408 that round-trip
2411// is 71 MB written + 71 MB read per layer -- 9.1 GB over 64 layers -- for a value consumed exactly
2412// once, immediately, by the very next kernel. Folding the quantisation into the producer deletes
2413// both halves and the kernel launch with them.
2414//
2415// `h` is recomputed in the second pass rather than held: a `float h[128]` is 512 B per thread and
2416// would spill the whole kernel to local memory, which costs far more than 128 extra __expf on a
2417// kernel that is bound by the 1 KB it reads per 128 outputs.
2418__global__ void bq_silu_had_fp8_dg(unsigned char* __restrict__ out, float* __restrict__ sf,
2419 const float* __restrict__ F, int rows, int n, int rows_aligned){
2420 const int kb = n / BQ_DG_BLOCK_K;
2421 for (long idx = (long)blockIdx.x * blockDim.x + threadIdx.x; idx < (long)rows * kb;
2422 idx += (long)gridDim.x * blockDim.x) {
2423 const int r = (int)(idx / kb), b = (int)(idx % kb);
2424 const float* g = F + (long)r * 2 * n + (long)b * BQ_DG_BLOCK_K;
2425 const float* u = g + n;
2426 float amax = 0.f;
2427 for (int i = 0; i < BQ_DG_BLOCK_K; ++i) {
2428 const float gv = g[i];
2429 amax = fmaxf(amax, fabsf((gv / (1.0f + __expf(-gv))) * u[i]));
2430 }
2431 const float S = amax > 0.f ? amax / BQ_DG_E4M3_MAX : 1.f;
2432 sf[(long)b * rows_aligned + r] = S;
2433 for (int i = 0; i < BQ_DG_BLOCK_K; ++i) {
2434 const float gv = g[i];
2435 const float h = (gv / (1.0f + __expf(-gv))) * u[i];
2436#if BQ_FP8_TOOLKIT
2437 out[(long)r * n + (long)b * BQ_DG_BLOCK_K + i] =
2438 (unsigned char)__nv_cvt_float_to_fp8(h / S, __NV_SATFINITE, __NV_E4M3);
2439#else
2440 out[(long)r * n + (long)b * BQ_DG_BLOCK_K + i] = 0;
2441#endif
2442 }
2443 }
2444}
2445
2446__global__ void bq_quant_act_e4m3(unsigned char* __restrict__ out, float* __restrict__ row_scale,
2447 const float* __restrict__ X, int n){
2448 int b = blockIdx.x;
2449 const float* x = X + (long)b*n;
2450 __shared__ float red[256];
2451 float m = 0.f;
2452 for(int i = threadIdx.x; i < n; i += blockDim.x) m = fmaxf(m, fabsf(x[i]));
2453 red[threadIdx.x] = m; __syncthreads();
2454 for(int t = blockDim.x/2; t > 0; t >>= 1){
2455 if(threadIdx.x < t) red[threadIdx.x] = fmaxf(red[threadIdx.x], red[threadIdx.x+t]);
2456 __syncthreads(); }
2457 float sc = (red[0] > 0.f) ? (red[0] / 448.0f) : 1.0f;
2458 if(threadIdx.x == 0) row_scale[b] = sc;
2459 for(int i = threadIdx.x; i < n; i += blockDim.x)
2460#if BQ_FP8_TOOLKIT
2461 out[(long)b*n + i] = (unsigned char)__nv_cvt_float_to_fp8(x[i]/sc, __NV_SATFINITE, __NV_E4M3);
2462#else
2463 out[(long)b*n + i] = 0;
2464#endif
2465}
2466// ⚠ UNVERIFIED AND NOT SELECTED. This tile COMPILES for sm_90a and is launched by nothing; it is
2467// kept because the surrounding pieces (the e4m3 mma helper, the Q8_0->e4m3 requantiser, the per-row
2468// activation quantiser, the toolkit guards) are verified and worth keeping, and because the
2469// measurement that stopped it is worth recording next to it.
2470//
2471// MEASURED on an H100, 2026-09-21: rel_L2 vs an f64 oracle = 1.001 against the int8 lane's 5.3e-3.
2472// A rel_L2 of ~1 means the output carries no signal. The warp/fragment mapping was first written
2473// from memory (tig=(lane&3)*2, jb=ct*8+tig) and then copied EXACTLY from the verified int8 tile
2474// (tig=lane&3, jb=ct*8+gid) -- which moved the answer but did not fix it, so what remains is the
2475// e4m3 instruction's own operand semantics differing from the s8 form the layout was derived under.
2476//
2477// It is also not competitive as written: 359 GB/s against the int8 lane's 906 at the decode shape,
2478// because it has none of the machinery the shipped tile has -- no cp.async staging, no pipeline, no
2479// split-K. And the projection that motivated it tops out at 0.93x of vLLM even when finished, so
2480// FINISHING IT DOES NOT WIN B=32. The measured win is the Q4 lane (1.33x).
2481//
2482// The fp8 tile. Deliberately the SAME geometry, swizzle and staging as `bq_gemm_tc1` -- the fp8 mma
2483// carries an identical register shape (4 A regs, 2 B regs, 4 accumulators), so nothing about the
2484// fragment map or the shared layout changes. What changes is that the accumulator is fp32 and the
2485// mma accumulates into it DIRECTLY across the whole K, with the two scales applied once at the end.
2486// That removes the 4:1 rescale-to-mma ratio the in-kernel probe measured on Hopper.
2487// Opt-in, like the wgmma tile: the packed staging buffer sits BESIDE the int8 tile, so this
2488// instantiation needs 1.5x the shared memory and blows the 48 KB budget at the widths real programs
2489// use. Compiling it unconditionally broke every emitted program -- exactly what BQ_TC_LANES was
2490// introduced to stop for the two-pass wrapper. It is NOT SELECTED (measured slower in all three
2491// arrangements), so a shipped program should not carry it at all.
2492#ifdef BQ_ENABLE_Q4
2493__global__ void __launch_bounds__(256, TC_MINBLOCKS) bq_gemm_tc1_q4(
2494 float* __restrict__ Y, const signed char* __restrict__ wq4, const __half* __restrict__ ws,
2495 const signed char* __restrict__ Xh,
2496 const float* __restrict__ Xs, int K, int N, int B, int YS, float* __restrict__ P, int SPLIT){
2497 bq_gemm_tc_body<false,false,true>(Y,wq4,ws,Xh,(const signed char*)0,Xs,K,N,B,YS,P,SPLIT);
2498}
2499#endif // BQ_ENABLE_Q4
2500#define FP8_TI (16*(8/(TC_TJ/8))*TC_WM)
2501__global__ void __launch_bounds__(256, TC_MINBLOCKS) bq_gemm_fp8(
2502 float* __restrict__ Y, const unsigned char* __restrict__ wq, const float* __restrict__ wsc,
2503 const unsigned char* __restrict__ Xq, const float* __restrict__ xsc,
2504 int K, int N, int B, int YS){
2505 const int NCT=TC_TJ/8, NRT=8/NCT, TI=16*NRT*TC_WM, SW=TC_KC/4;
2506 __shared__ unsigned int sW[TI][TC_KC/4];
2507 __shared__ unsigned int sH[TC_TJ][TC_KC/4];
2508 const int tid=threadIdx.x, i0=blockIdx.x*TI, b0=blockIdx.z*TC_TJ;
2509 const int warp=tid>>5, lane=tid&31;
2510 const int rt=warp/NCT, ct=warp%NCT, gid=lane>>2, tig=lane&3;
2511 float acc[TC_WM][4];
2512 #pragma unroll
2513 for(int m=0;m<TC_WM;m++){ acc[m][0]=acc[m][1]=acc[m][2]=acc[m][3]=0.f; }
2514 for(int k0=0;k0<K;k0+=TC_KC){
2515 for(int idx=tid; idx<TI*(TC_KC/16); idx+=256){
2516 int r=idx/(TC_KC/16), c=idx%(TC_KC/16), o=i0+r;
2517 int4 v = (o<N && k0+c*16<K) ? ((const int4*)(wq+(long)o*K+k0))[c] : make_int4(0,0,0,0);
2518 *(int4*)&sW[r][BQTC_SWZC(r,c)*4] = v;
2519 }
2520 for(int idx=tid; idx<TC_TJ*(TC_KC/16); idx+=256){
2521 int r=idx/(TC_KC/16), c=idx%(TC_KC/16);
2522 int4 v = (r+b0<B && k0+c*16<K) ? ((const int4*)(Xq+(long)(r+b0)*K+k0))[c] : make_int4(0,0,0,0);
2523 *(int4*)&sH[r][BQTC_SWZC(r,c)*4] = v;
2524 }
2525 __syncthreads();
2526 #pragma unroll
2527 for(int blk=0; blk<TC_KC/32; ++blk){
2528 const int jb=ct*8+gid;
2529 const unsigned b0h=sH[jb][BQTC_SWZJ(jb,blk*8+tig)], b1h=sH[jb][BQTC_SWZJ(jb,blk*8+tig+4)];
2530 #pragma unroll
2531 for(int m=0;m<TC_WM;m++){
2532 const int rb=(rt*TC_WM+m)*16;
2533 const unsigned a0=sW[rb+gid][BQTC_SWZJ(rb+gid,blk*8+tig)],
2534 a1=sW[rb+gid+8][BQTC_SWZJ(rb+gid+8,blk*8+tig)],
2535 a2=sW[rb+gid][BQTC_SWZJ(rb+gid,blk*8+tig+4)],
2536 a3=sW[rb+gid+8][BQTC_SWZJ(rb+gid+8,blk*8+tig+4)];
2537 bq_mma_e4m3(acc[m][0],acc[m][1],acc[m][2],acc[m][3], a0,a1,a2,a3, b0h,b1h);
2538 }
2539 }
2540 __syncthreads();
2541 }
2542 #pragma unroll
2543 for(int m=0;m<TC_WM;m++){
2544 const int rb=(rt*TC_WM+m)*16, o0=i0+rb+gid, o1=o0+8;
2545 const int j0=ct*8+tig*2+b0, j1=j0+1;
2546 if(o0<N){ if(j0<B) Y[(long)j0*YS+o0]=acc[m][0]*wsc[o0]*xsc[j0];
2547 if(j1<B) Y[(long)j1*YS+o0]=acc[m][1]*wsc[o0]*xsc[j1]; }
2548 if(o1<N){ if(j0<B) Y[(long)j0*YS+o1]=acc[m][2]*wsc[o1]*xsc[j0];
2549 if(j1<B) Y[(long)j1*YS+o1]=acc[m][3]*wsc[o1]*xsc[j1]; }
2550 }
2551}
2552#endif
2553
2554// ══ THE RESIDUAL NORM, SPLIT SO IT CAN FILL THE MACHINE ═════════════════════════════════════════
2555// `bq_add_rmsnorm_q8f_hl` takes ONE BLOCK PER SEQUENCE, because the sum of squares spans the model
2556// dimension. At two sequences that is two blocks of a 132-SM part, and it measured 4.06 ms of a
2557// 21.5 ms step — 127 launches at ~32 us each moving ~200 KB apiece. It is latency, not work.
2558//
2559// Two things were tried before this one and are recorded because the reasons transfer. Widening the
2560// block 256 -> 1024 recovered 38% and nothing more. Splitting the QUANTISER out into a grid-parallel
2561// kernel — which is how llama.cpp is arranged, its `quantize_row_q8_1` being separate from its
2562// row-serial `rms_norm_f32` — measured 0.77x, SLOWER: the activation is then read from global twice
2563// and the launch count doubles. Their arrangement is right for their kernel and wrong for ours,
2564// because ours already fuses the add, the norm, the f32 write and the two-word quantiser into one
2565// pass. Ours is not badly built; it is grid-starved.
2566//
2567// So this splits the REDUCTION instead. Phase one computes the residual add and a per-chunk partial
2568// across NORM_G*B blocks; phase two re-sums those partials — every block redundantly, in a fixed
2569// ascending order, so all blocks agree bit for bit and the result is deterministic run to run — then
2570// normalises and quantises its own chunk. Measured 1.95x-2.02x at one, two and four sequences.
2571//
2572// ★ It is NOT bit-identical to the one-block form: NORM_G chunk-sums instead of RMSN_NT thread-sums
2573// is a different association, worth ~1 ulp (worst relative 2.4e-07 measured). That is five orders of
2574// magnitude below the 1.5e-2 spread between llama.cpp's own `-fa 0` and `-fa 1`, so it cannot move
2575// the accuracy standing — but it is a different number, so it lives behind a knob and the exact lane
2576// stays the default.
2577__global__ void __launch_bounds__(RMSN_NT) bq_add_rmsnorm_part(
2578 float* __restrict__ X, const float* __restrict__ Y, float* __restrict__ part, int n, int G){
2579 const int bb=blockIdx.y, g=blockIdx.x;
2580 const int chunk=((n+G-1)/G+31)/32*32;
2581 const int lo=g*chunk; int hi=lo+chunk; if(hi>n) hi=n;
2582 float* x=X+(long)bb*n; const float* y=Y+(long)bb*n;
2583 __shared__ float sh[RMSN_NT];
2584 float s=0.f;
2585 for(int i=lo+threadIdx.x;i<hi;i+=blockDim.x){ float v=x[i]+y[i]; x[i]=v; s+=v*v; }
2586 sh[threadIdx.x]=s; __syncthreads();
2587 for(int k=blockDim.x/2;k>0;k>>=1){ if(threadIdx.x<k) sh[threadIdx.x]+=sh[threadIdx.x+k]; __syncthreads(); }
2588 if(threadIdx.x==0) part[(long)bb*G+g]=sh[0];
2589}
2590__global__ void __launch_bounds__(RMSN_NT) bq_add_rmsnorm_apply(
2591 const float* __restrict__ X, const float* __restrict__ part, signed char* __restrict__ Xq,
2592 signed char* __restrict__ Xl, float* __restrict__ Xs, float* __restrict__ Xf,
2593 const float* __restrict__ w, int n, float eps, int G){
2594 const int bb=blockIdx.y, g=blockIdx.x;
2595 const int chunk=((n+G-1)/G+31)/32*32;
2596 const int lo=g*chunk; int hi=lo+chunk; if(hi>n) hi=n;
2597 const float* x=X+(long)bb*n; float* xf=Xf+(long)bb*n;
2598 float tot=0.f;
2599 for(int j=0;j<G;++j) tot+=part[(long)bb*G+j];
2600 const float sc=rsqrtf(tot/n+eps);
2601 for(int i=lo+threadIdx.x;i<hi;i+=blockDim.x) xf[i]=x[i]*sc*w[i];
2602 __syncthreads();
2603 const int nb=n/32;
2604 for(int b=lo/32+threadIdx.x;b<hi/32;b+=blockDim.x){
2605 float amax=1e-30f;
2606 for(int j=0;j<32;++j){ float a=fabsf(xf[b*32+j]); if(a>amax) amax=a; }
2607 float q=amax/127.0f; Xs[(long)bb*nb+b]=q;
2608 for(int j=0;j<32;++j) lgx_q8hl(xf[b*32+j]/q, Xq+(long)bb*n, Xl+(long)bb*n, (long)b*32+j);
2609 }
2610}
2611
2612// ══ THE SAME matmul AGAIN, ROW-PARALLEL ═════════════════════════════════════════════════════════
2613// Below a handful of sequences the tensor-core tile is the wrong shape: `mma.m16n8k32` has an 8-wide
2614// N, so at two sequences six of every eight columns are zeros and the mma rate stops being an
2615// advantage. What matters there is streaming the weights, and the shipped `bq_gemv_hl` streams them
2616// badly for a reason that is structural rather than arithmetic: it gives one output row to ONE warp,
2617// so only 32 threads' worth of loads are ever in flight per row, and it re-reads the activation slab
2618// once per row.
2619//
2620// llama.cpp does not have that problem, and reading their kernel says why. ggml/src/ggml-cuda/mmvq.cu
2621// picks, for Q8_0 on an NVIDIA part, `nwarps` = 4 at ncols_dst <= 4 and 2 at <= 8, with
2622// `rows_per_cuda_block` = 2 (`calc_nwarps`/`calc_rows_per_block`), and routes every batch <= 8 there
2623// (`MMVQ_MAX_BATCH_SIZE`, checked before MMQ in ggml_cuda_mul_mat). So a whole BLOCK reduces one row,
2624// and two rows share each activation load. The structure below is that structure; the code is ours,
2625// and it keeps the two-pass hi/lo accuracy tier and the K-slicing that mmvq does not have.
2626//
2627// MV_NW=1 with MV_RB=2 is BIT-IDENTICAL to `bq_gemv_hl` — same ascending K order, same grouping —
2628// and measured 2.5x-2.9x faster at eight sequences purely by halving activation traffic.
2629__global__ void __launch_bounds__(MV_NW*32) bq_gemv_wide(
2630 float* __restrict__ Y, const signed char* __restrict__ wq, const __half* __restrict__ ws,
2631 const signed char* __restrict__ Xh, const signed char* __restrict__ Xl,
2632 const float* __restrict__ Xs, int K, int N, int B, int YS,
2633 float* __restrict__ P, int SPLIT){
2634 static_assert(MV_TJ>=1 && MV_RB>=1 && MV_NW>=1, "the row-parallel GEMV needs a positive tile");
2635 // A tile narrower than the batch would compute only MV_TJ of the B sequences and leave the rest
2636 // holding whatever Y held before — plausible numbers, wrong answer. Refuse instead.
2637 if(MV_TJ < B) return;
2638 const int tid=threadIdx.x, lane=tid&31, wp=tid>>5;
2639 const int row0=blockIdx.x*MV_RB, KB=K/32, nq=K/512;
2640 const int per=(nq+SPLIT-1)/SPLIT, c0=blockIdx.y*per;
2641 const int c1=(c0+per)<nq?(c0+per):nq;
2642 float acc[MV_RB][MV_TJ];
2643 #pragma unroll
2644 for(int r=0;r<MV_RB;r++){
2645 #pragma unroll
2646 for(int j=0;j<MV_TJ;j++) acc[r][j]=0.f; }
2647 for(int c=c0+wp; c<c1; c+=MV_NW){
2648 const int blk=c*16+lane/2;
2649 int4 hh[MV_TJ], ll[MV_TJ]; float xsv[MV_TJ];
2650 #pragma unroll
2651 for(int j=0;j<MV_TJ;j++){
2652 if(j<B){ hh[j]=((const int4*)(Xh+(long)j*K))[c*32+lane];
2653 ll[j]=((const int4*)(Xl+(long)j*K))[c*32+lane];
2654 xsv[j]=Xs[(long)j*KB+blk]; }
2655 else { hh[j]=make_int4(0,0,0,0); ll[j]=make_int4(0,0,0,0); xsv[j]=0.f; } }
2656 #pragma unroll
2657 for(int r=0;r<MV_RB;r++){
2658 const int o=row0+r; if(o>=N) continue;
2659 const int4 w=((const int4*)(wq+(long)o*K))[c*32+lane];
2660 const float wsb=__half2float(ws[(long)o*KB+blk]);
2661 #pragma unroll
2662 for(int j=0;j<MV_TJ;j++){
2663 if(j>=B) continue;
2664 int dh=__dp4a(w.x,hh[j].x,0)+__dp4a(w.y,hh[j].y,0)+__dp4a(w.z,hh[j].z,0)+__dp4a(w.w,hh[j].w,0);
2665 int dl=__dp4a(w.x,ll[j].x,0)+__dp4a(w.y,ll[j].y,0)+__dp4a(w.z,ll[j].z,0)+__dp4a(w.w,ll[j].w,0);
2666 acc[r][j] += ((float)dh + (float)dl*(1.0f/127.0f)) * (wsb*xsv[j]);
2667 } }
2668 }
2669 // The K remainder past the last full 512-wide quad, on the slice that owns the end of the sum.
2670 if(blockIdx.y==SPLIT-1){
2671 const int done=nq*512;
2672 for(int r=0;r<MV_RB;r++){ const int o=row0+r; if(o>=N) continue;
2673 for(int j=0;j<MV_TJ;j++){ if(j>=B) continue;
2674 for(int q=done+tid; q<K; q+=MV_NW*32){ const int blk=q/32;
2675 float a=(float)Xh[(long)j*K+q] + (float)Xl[(long)j*K+q]*(1.0f/127.0f);
2676 acc[r][j] += a*(float)wq[(long)o*K+q]
2677 *(__half2float(ws[(long)o*KB+blk])*Xs[(long)j*KB+blk]); } } }
2678 }
2679 __shared__ float red[MV_NW][MV_RB][MV_TJ];
2680 #pragma unroll
2681 for(int r=0;r<MV_RB;r++){
2682 #pragma unroll
2683 for(int j=0;j<MV_TJ;j++){
2684 float v=acc[r][j];
2685 for(int t=16;t;t>>=1) v+=__shfl_down_sync(0xffffffffu,v,t);
2686 if(lane==0) red[wp][r][j]=v; } }
2687 __syncthreads();
2688 if(tid < MV_RB*MV_TJ){
2689 const int r=tid/MV_TJ, j=tid%MV_TJ, o=row0+r;
2690 if(o<N && j<B){
2691 float v=0.f;
2692 #pragma unroll
2693 for(int w=0;w<MV_NW;w++) v+=red[w][r][j];
2694 if(SPLIT==1) Y[(long)j*YS+o]=v;
2695 else P[(long)blockIdx.y*N*B + (long)j*N + o]=v;
2696 } }
2697}
2698
2699// Fixed ascending order over the K slices, so the split reduction is deterministic run to run.
2700__global__ void bq_sk_reduce(float* __restrict__ Y, const float* __restrict__ P,
2701 int N, int B, int YS, int S){
2702 long i=(long)blockIdx.x*blockDim.x+threadIdx.x; if(i>=(long)N*B) return;
2703 int j=(int)(i/N), o=(int)(i%N); float v=0.f;
2704 for(int s=0;s<S;s++) v+=P[(long)s*N*B+(long)j*N+o];
2705 Y[(long)j*YS+o]=v;
2706}
2707"#;
2708
2709impl Qwen38Cfg {
2710 pub fn state_bytes_per_seq(&self, maxseq: usize) -> usize {
2713 let gdn = (self.layers - self.attn_layers()) * self.gdn_v() * self.gdn_head_dim * 4;
2714 let conv = (self.layers - self.attn_layers())
2715 * self.conv_taps
2716 * (2 * self.gdn_qk() + self.gdn_v())
2717 * 4;
2718 let kv = self.attn_layers() * 2 * maxseq * self.kv * self.hd * 4;
2719 gdn + conv + kv
2720 }
2721}
2722
2723pub fn emit_qwen38_generate_batched(
2730 cfg: &Qwen38Cfg,
2731 wdir: &str,
2732 prompt: &[usize],
2733 n_gen: usize,
2734 batch: usize,
2735 hl: bool,
2736) -> String {
2737 let pt = std::env::var("LOGOS_Q38B_PREFILL_T")
2738 .ok()
2739 .and_then(|s| s.parse().ok())
2740 .unwrap_or(1);
2741 emit_qwen38_generate_batched_pt(cfg, wdir, prompt, n_gen, batch, hl, pt)
2742}
2743
2744pub fn emit_qwen38_generate_batched_pt(
2751 cfg: &Qwen38Cfg,
2752 wdir: &str,
2753 prompt: &[usize],
2754 n_gen: usize,
2755 batch: usize,
2756 hl: bool,
2757 prefill_t: usize,
2758) -> String {
2759 emit_qwen38_generate_batched_mtp(
2760 cfg,
2761 wdir,
2762 prompt,
2763 n_gen,
2764 batch,
2765 hl,
2766 prefill_t,
2767 std::env::var("LOGOS_Q38B_MTP").is_ok(),
2768 )
2769}
2770
2771fn replace_word(src: &str, from: &str, to: &str) -> String {
2778 let is_ident = |c: char| c.is_ascii_alphanumeric() || c == '_';
2779 let (mut out, mut i) = (String::with_capacity(src.len()), 0usize);
2780 while let Some(k) = src[i..].find(from) {
2781 let at = i + k;
2782 let before_ok = at == 0 || !src[..at].chars().next_back().is_some_and(&is_ident);
2783 let after = at + from.len();
2784 let after_ok = after >= src.len() || !src[after..].chars().next().is_some_and(&is_ident);
2785 out.push_str(&src[i..at]);
2786 if before_ok && after_ok {
2787 out.push_str(to);
2788 } else {
2789 out.push_str(from);
2790 }
2791 i = after;
2792 }
2793 out.push_str(&src[i..]);
2794 out
2795}
2796
2797fn dg_site(
2808 site: usize,
2809 out: &str,
2810 wvar: &str,
2811 act_f32: &str,
2812 k: usize,
2813 n: usize,
2814 ys: usize,
2815 orig: &str,
2816) -> String {
2817 dg_site_p(
2818 site,
2819 out,
2820 wvar,
2821 &format!("bq_quant_act_fp8_dg<<<1024,256,0,s>>>(dgA,dgAsf,{act_f32},_m,{k},_ma);"),
2822 k,
2823 n,
2824 ys,
2825 orig,
2826 )
2827}
2828
2829fn dg_site_p(
2833 site: usize,
2834 out: &str,
2835 wvar: &str,
2836 producer: &str,
2837 k: usize,
2838 n: usize,
2839 ys: usize,
2840 orig: &str,
2841) -> String {
2842 format!(
2843 "\n\x20 if(!bq_dg_on[{site}]) {{\n\
2844 \x20 {orig}\n\
2845 \x20 }} else {{\n\
2846 #ifdef BQ_ENABLE_DG\n\
2847 \x20 {{ const int _m = BTv, _ma = bq_dg_align(_m,4);\n\
2848 \x20 {producer}\n\
2849 \x20 CUtensorMap _ta = bq_dg_ab(dgA,_m,{k},BQ_DG_BM);\n\
2850 \x20 CUtensorMap _tsa = bq_dg_sf(dgAsf,_m,{k},BQ_DG_BM);\n\
2851 \x20 CUtensorMap _tcd = bq_dg_cd({out},_m,{n},{ys});\n\
2852 \x20 cudaFuncSetAttribute((const void*)&sm90_fp8_gemm_1d1d_impl<0,{n},{k},1,\
2853BQ_DG_BM,BQ_DG_BN,BQ_DG_BLOCK_K,128,128,BQ_DG_STAGES,128,128,1,true,132,GemmType::Normal,float>,\n\
2854 \x20 cudaFuncAttributeMaxDynamicSharedMemorySize, BQ_DG_SMEM);\n\
2855 \x20 sm90_fp8_gemm_1d1d_impl<0,{n},{k},1,BQ_DG_BM,BQ_DG_BN,BQ_DG_BLOCK_K,128,128,\
2856BQ_DG_STAGES,128,128,1,true,132,GemmType::Normal,float>\n\
2857 \x20 <<<132,256,BQ_DG_SMEM,s>>>((__nv_fp8_e4m3*)dgA,(__nv_fp8_e4m3*){wvar}_f8,\n\
2858 \x20 nullptr,nullptr,(uint32_t)_m,(uint32_t){n},(uint32_t){k},\n\
2859 \x20 _ta,{wvar}_tb,_tsa,{wvar}_tsb,_tcd);\n\
2860 \x20 if(bq_dg_hl) {{\n\
2861 \x20 CUtensorMap _tcl = bq_dg_cd(dgD,_m,{n},{n});\n\
2862 \x20 sm90_fp8_gemm_1d1d_impl<0,{n},{k},1,BQ_DG_BM,BQ_DG_BN,BQ_DG_BLOCK_K,128,128,\
2863BQ_DG_STAGES,128,128,1,true,132,GemmType::Normal,float>\n\
2864 \x20 <<<132,256,BQ_DG_SMEM,s>>>((__nv_fp8_e4m3*)dgA,(__nv_fp8_e4m3*){wvar}_f8l,\n\
2865 \x20 nullptr,nullptr,(uint32_t)_m,(uint32_t){n},(uint32_t){k},\n\
2866 \x20 _ta,{wvar}_tbl,_tsa,{wvar}_tsbl,_tcl);\n\
2867 \x20 bq_dg_add<<<1024,256,0,s>>>({out},{ys},dgD,{n},_m);\n\
2868 \x20 }} }}\n\
2869 #endif\n\
2870 \x20 }}\n"
2871 )
2872}
2873
2874fn dg_weight(wvar: &str, n: usize, k: usize) -> (String, String) {
2877 (format!("static unsigned char* {wvar}_f8, *{wvar}_f8l; static float* {wvar}_sf, *{wvar}_sfl;\n\
2878 static CUtensorMap {wvar}_tb, {wvar}_tsb, {wvar}_tbl, {wvar}_tsbl;\n"),
2879 format!(" cudaMalloc(&{wvar}_f8,(size_t){n}*{k});\n\
2880 \x20 cudaMalloc(&{wvar}_sf,(size_t)({k}/BQ_DG_BLOCK_K)*bq_dg_align({n},4)*4);\n\
2881 \x20 cudaMalloc(&{wvar}_f8l,(size_t){n}*{k});\n\
2882 \x20 cudaMalloc(&{wvar}_sfl,(size_t)({k}/BQ_DG_BLOCK_K)*bq_dg_align({n},4)*4);\n\
2883 \x20 bq_q8_to_fp8_dg2<<<1024,256>>>({wvar}_f8,{wvar}_sf,{wvar}_f8l,{wvar}_sfl,{wvar}_q,{wvar}_s,{n},{k},bq_dg_align({n},4));\n\
2884 \x20 cudaDeviceSynchronize();\n\
2885 \x20 {wvar}_tb = bq_dg_ab({wvar}_f8,{n},{k},BQ_DG_BN);\n\
2886 \x20 {wvar}_tsb = bq_dg_sf({wvar}_sf,{n},{k},BQ_DG_BN);\n\
2887 \x20 {wvar}_tbl = bq_dg_ab({wvar}_f8l,{n},{k},BQ_DG_BN);\n\
2888 \x20 {wvar}_tsbl = bq_dg_sf({wvar}_sfl,{n},{k},BQ_DG_BN);\n"))
2889}
2890
2891pub fn emit_qwen38_generate_batched_mtp(
2892 cfg: &Qwen38Cfg,
2893 wdir: &str,
2894 prompt: &[usize],
2895 n_gen: usize,
2896 batch: usize,
2897 hl: bool,
2898 prefill_t: usize,
2899 mtp: bool,
2900) -> String {
2901 emit_qwen38_generate_batched_q4(
2902 cfg,
2903 wdir,
2904 prompt,
2905 n_gen,
2906 batch,
2907 hl,
2908 prefill_t,
2909 mtp,
2910 std::env::var("LOGOS_Q38B_Q4").ok().as_deref() == Some("1"),
2911 )
2912}
2913
2914pub fn emit_qwen38_generate_batched_q4(
2916 cfg: &Qwen38Cfg,
2917 wdir: &str,
2918 prompt: &[usize],
2919 n_gen: usize,
2920 batch: usize,
2921 hl: bool,
2922 prefill_t: usize,
2923 mtp: bool,
2924 q4: bool,
2925) -> String {
2926 let (d, dff, taps) = (cfg.d, cfg.dff, cfg.conv_taps);
2927 let (qk, vw) = (cfg.gdn_qk(), cfg.gdn_v());
2928 let (nk, nv, ghd) = (cfg.gdn_k_heads, cfg.gdn_v_heads, cfg.gdn_head_dim);
2929 let (aq, akv) = (cfg.aq(), cfg.akv());
2930 let convc = 2 * qk + vw;
2931 let projw = 2 * qk + 2 * vw;
2932 let attnw = 2 * aq + 2 * akv;
2933 let maxseq = prompt.len() + n_gen + 128;
2934 let amax = aq.max(vw);
2935 let btiles = batch.div_ceil(BTG);
2936
2937 let mma_min: usize = std::env::var("LOGOS_Q38B_MMA_MIN")
2952 .ok()
2953 .and_then(|s| s.parse().ok())
2954 .unwrap_or(64);
2955 let tc_min: usize = std::env::var("LOGOS_Q38B_TC_MIN")
2962 .ok()
2963 .and_then(|s| s.parse().ok())
2964 .unwrap_or(2);
2965 let tc_window = std::env::var("LOGOS_Q38B_TC").ok().as_deref() != Some("0")
2966 && batch >= tc_min
2967 && batch <= 64;
2968 let use_mma = !hl && batch >= mma_min && !tc_window;
2969 let use_staged =
2987 std::env::var("LOGOS_Q38B_STAGED").ok().as_deref() == Some("1") && !use_mma && batch >= 2;
2988
2989 let wide_max: usize = std::env::var("LOGOS_Q38B_WIDE_MAX")
3026 .ok()
3027 .and_then(|s| s.parse().ok())
3028 .unwrap_or(2);
3029 let use_wide = hl
3030 && !use_mma
3031 && !use_staged
3032 && batch <= wide_max
3033 && std::env::var("LOGOS_Q38B_WIDE").ok().as_deref() != Some("0");
3034 let mv_tj = if use_wide {
3040 batch.next_power_of_two().clamp(1, 8)
3041 } else {
3042 8
3043 };
3044 let mv_nw = (8 / mv_tj).clamp(1, 4);
3045 let mv_rb = 2usize;
3046 let rmsn_nt: usize = std::env::var("LOGOS_Q38B_RMSN_NT")
3055 .ok()
3056 .and_then(|s| s.parse().ok())
3057 .unwrap_or(256);
3058 let norm_split = std::env::var("LOGOS_Q38B_NORM").ok().as_deref() == Some("split");
3073 let norm_g: usize = std::env::var("LOGOS_Q38B_NORM_G")
3074 .ok()
3075 .and_then(|s| s.parse().ok())
3076 .unwrap_or(20);
3077 let use_tc = !use_mma && !use_staged && !use_wide && tc_window;
3083 let use_q4 = q4 && use_tc && !hl;
3103 let tc_tj = if use_tc {
3104 let tj_rows = if mtp && prefill_t >= 2 {
3110 2 * batch
3111 } else {
3112 batch
3113 };
3114 let auto = tj_rows.next_power_of_two().clamp(8, 64);
3115 match std::env::var("LOGOS_Q38B_TC_TJ")
3116 .ok()
3117 .and_then(|v| v.parse::<usize>().ok())
3118 {
3119 Some(v) if v.is_power_of_two() && (8..=64).contains(&v) && v >= tj_rows => v,
3120 _ => auto,
3121 }
3122 } else {
3123 32
3124 };
3125 let tc_kc_env: Option<usize> = std::env::var("LOGOS_Q38B_TC_KC")
3136 .ok()
3137 .and_then(|s| s.parse().ok());
3138 let tc_occ_def = match std::env::var("LOGOS_Q38B_TC_OCC")
3162 .ok()
3163 .and_then(|s| s.parse::<usize>().ok())
3164 {
3165 Some(n) => format!("#define TC_MINBLOCKS {n}"),
3166 None => [
3167 "#if defined(__CUDA_ARCH__) && __CUDA_ARCH__ >= 900",
3168 "#define TC_MINBLOCKS 3",
3169 "#else",
3170 "#define TC_MINBLOCKS 2",
3171 "#endif",
3172 ]
3173 .join("\n"),
3174 };
3175 let prefill_t: usize = prefill_t.max(1);
3188 let mtp: bool = mtp && prefill_t >= 2;
3191 let wm_start: usize = 4;
3211 let act = if hl { 2 } else { 1 };
3219 let fits = |wm: usize, kc: usize, st: usize| -> bool {
3220 let (ti, sw) = (16 * (8 / (tc_tj / 8)) * wm, kc / 4);
3221 let q4_extra = if use_q4 { st * ti * (kc / 8) * 4 } else { 0 };
3228 st * (ti * sw * 4 + act * tc_tj * sw * 4 + tc_tj * (kc / 32) * 4 + ti * (kc / 32) * 2)
3229 + q4_extra
3230 <= 48 * 1024
3231 };
3232 let kc_base: usize = if batch >= 32 { 64 } else { 128 };
3238 let tc_wm: usize = std::env::var("LOGOS_Q38B_TC_WM")
3239 .ok()
3240 .and_then(|s| s.parse().ok())
3241 .unwrap_or_else(|| {
3242 let sw = kc_base / 4 + 1;
3243 *[wm_start, 4, 2, 1]
3244 .iter()
3245 .find(|wm| {
3246 let ti = 16 * (8 / (tc_tj / 8)) * *wm;
3247 ti * sw * 4
3248 + 2 * tc_tj * sw * 4
3249 + tc_tj * (kc_base / 32) * 4
3250 + ti * (kc_base / 32) * 2
3251 <= 48 * 1024
3252 })
3253 .expect("some sub-tile height fits the shared-memory budget")
3254 });
3255 let tc_ti = 16 * (8 / (tc_tj / 8)) * tc_wm;
3256 let tc_kc: usize = tc_kc_env.unwrap_or_else(|| {
3262 if kc_base < 128 && fits(tc_wm, 128, 2) {
3263 128
3264 } else {
3265 kc_base
3266 }
3267 });
3268 let stage_bytes_kc = |kc: usize, st: usize| -> usize {
3286 let sw = kc / 4;
3287 st * (tc_ti * sw * 4 + act * tc_tj * sw * 4 + tc_tj * (kc / 32) * 4 + tc_ti * (kc / 32) * 2)
3288 };
3289 let q4_stage_bytes =
3300 |kc: usize, st: usize| -> usize { stage_bytes_kc(kc, st) + st * tc_ti * (kc / 8) * 4 };
3301 let tc_kc = if use_q4 {
3302 let mut kc = tc_kc;
3303 while kc > 32 && q4_stage_bytes(kc, 2) > 48 * 1024 {
3304 kc /= 2;
3305 }
3306 kc
3307 } else {
3308 tc_kc
3309 };
3310 let stage_bytes = |st: usize| -> usize { stage_bytes_kc(tc_kc, st) };
3311 let tc_stages: usize = std::env::var("LOGOS_Q38B_TC_STAGES")
3312 .ok()
3313 .and_then(|s| s.parse().ok())
3314 .unwrap_or(if use_q4 {
3315 2
3316 } else if stage_bytes(2) <= 48 * 1024 {
3317 2
3318 } else {
3319 1
3320 });
3321 let tc_target: usize = std::env::var("LOGOS_Q38B_TC_TARGET")
3350 .ok()
3351 .and_then(|s| s.parse().ok())
3352 .unwrap_or(1024);
3353 let tc_split_max: usize = std::env::var("LOGOS_Q38B_TC_SPLIT_MAX")
3361 .ok()
3362 .and_then(|s| s.parse().ok())
3363 .filter(|v| *v >= 1)
3364 .unwrap_or(16);
3365 let tc_split = move |k: usize, n: usize| -> usize {
3366 if !use_tc {
3367 return 1;
3368 }
3369 let base = n.div_ceil(tc_ti);
3370 tc_target
3371 .div_ceil(base)
3372 .clamp(1, tc_split_max)
3373 .min(((k / tc_kc) / 2).max(1))
3374 };
3375 let tc_part_max = std::cell::Cell::new(0usize);
3377
3378 let mmq_min: usize = std::env::var("LOGOS_Q38B_MMQ_MIN")
3379 .ok()
3380 .and_then(|s| s.parse().ok())
3381 .unwrap_or(32);
3382 let use_mmq = !use_mma && !use_staged && !use_tc && batch >= mmq_min;
3383 let mmq_tj = if use_mmq { batch.div_ceil(2) * 2 } else { 32 };
3391 let mmq_acc = std::env::var("LOGOS_Q38B_MMQ_ACC").unwrap_or_else(|_| "double".to_string());
3392 let gemv_rows = |y: &str,
3393 w: &str,
3394 xh: &str,
3395 xs: &str,
3396 k: usize,
3397 n: usize,
3398 ys: usize,
3399 rows: &str|
3400 -> String {
3401 let nrows = if rows == "BT" {
3407 batch * prefill_t
3408 } else {
3409 batch
3410 };
3411 let btiles = nrows.div_ceil(BTG);
3412 if use_mma {
3413 return format!("bq_gemm_mma<<<dim3(({n}+TBN-1)/TBN,({nrows}+TBM-1)/TBM),256,(size_t)(TBM+TBN)*(KSTEP+4),s>>>({y},{w}_q,{w}_s,{xh},{xs},{k},{n},{rows},{ys});");
3415 }
3416 if use_wide {
3417 return format!("bq_gemv_wide<<<dim3(({n}+MV_RB-1)/MV_RB,1),MV_NW*32,0,s>>>({y},{w}_q,{w}_s,{xh},{xh}l,{xs},{k},{n},{rows},{ys},(float*)0,1);");
3420 }
3421 if use_tc {
3422 let sp = tc_split(k, n);
3423 let blocks = n.div_ceil(tc_ti);
3424 let use_q4 = use_q4 && w != "emb" && w != "Wout";
3429 let (tck, xops) = if use_q4 {
3430 ("bq_gemm_tc1_q4", format!("{xh}"))
3431 } else if hl {
3432 ("bq_gemm_tc", format!("{xh},{xh}l"))
3433 } else {
3434 ("bq_gemm_tc1", format!("{xh}"))
3435 };
3436 let (wq, ws) = if use_q4 {
3437 (format!("{w}_q4"), format!("{w}_s4"))
3438 } else {
3439 (format!("{w}_q"), format!("{w}_s"))
3440 };
3441 let mt = format!("(({rows})+TC_TJ-1)/TC_TJ");
3443 if sp > 1 {
3444 tc_part_max.set(tc_part_max.get().max(sp * n));
3445 return format!(
3446 "{tck}<<<dim3({blocks},{sp},{mt}),256,0,s>>>({y},{wq},{ws},{xops},{xs},{k},{n},{rows},{ys},skp,{sp});\n\
3447 \x20 bq_sk_reduce<<<Gs({n}*({rows})),256,0,s>>>({y},skp,{n},{rows},{ys},{sp});");
3448 }
3449 return format!("{tck}<<<dim3({blocks},1,{mt}),256,0,s>>>({y},{wq},{ws},{xops},{xs},{k},{n},{rows},{ys},(float*)0,1);");
3450 }
3451 if use_mmq {
3452 return format!("bq_gemm_mmq<<<(({n})+((256/(MMQ_TJ/MMQ_RJ))*MMQ_RI)-1)/((256/(MMQ_TJ/MMQ_RJ))*MMQ_RI),256,0,s>>>({y},{w}_q,{w}_s,{xh},{xh}l,{xs},{k},{n},{rows},{ys});");
3454 }
3455 if use_staged {
3456 if hl {
3458 return format!("bq_gemv_hl_smem<<<({n}+7)/8,256,(size_t)(2*{nrows}*512+{nrows}*16*4),s>>>({y},{w}_q,{w}_s,{xh},{xh}l,{xs},{k},{n},{rows},{ys});");
3459 }
3460 return format!("bq_gemv_h_smem<<<({n}+7)/8,256,(size_t)({nrows}*512+{nrows}*16*4),s>>>({y},{w}_q,{w}_s,{xh},{xs},{k},{n},{rows},{ys});");
3461 }
3462 if hl {
3463 format!("bq_gemv_hl<<<dim3(Gw({n}),{btiles}),256,0,s>>>({y},{w}_q,{w}_s,{xh},{xh}l,{xs},{k},{n},{rows},{ys});")
3464 } else {
3465 format!("bq_gemv_h<<<dim3(Gw({n}),{btiles}),256,0,s>>>({y},{w}_q,{w}_s,{xh},{xs},{k},{n},{rows},{ys});")
3466 }
3467 };
3468
3469 let gemv = |y: &str, w: &str, xh: &str, xs: &str, k: usize, n: usize, ys: usize| -> String {
3471 gemv_rows(y, w, xh, xs, k, n, ys, "B")
3472 };
3473
3474 let drop_reduce = |launch: String| -> String {
3480 match launch.rfind("bq_sk_reduce") {
3481 Some(i) => launch[..i].trim_end().to_string(),
3482 None => launch,
3483 }
3484 };
3485 let skf = std::env::var("LOGOS_Q38B_SKFUSE").unwrap_or_else(|_| "silu".to_string());
3492 let sk_fuse = skf != "0";
3493 let sk_norm = !norm_split && (skf == "all" || skf == "1");
3495 let wd_sp = if sk_norm { tc_split(dff, d) } else { 1 };
3496 let wo_attn_sp = if sk_norm { tc_split(aq, d) } else { 1 };
3497 let wo_gdn_sp = if sk_norm { tc_split(vw, d) } else { 1 };
3498 let gu_sp = if sk_fuse { tc_split(d, 2 * dff) } else { 1 };
3499
3500 let mut decls = String::new();
3501 let mut loads = String::new();
3502 let mut dg_decls = String::new();
3506 let mut dg_loads = String::new();
3507 let mut q8 = |var: &str, rows: usize, inn: usize, file: &str| {
3508 decls.push_str(&format!(
3509 "static signed char* {var}_q; static __half* {var}_s;\n"
3510 ));
3511 loads.push_str(&format!(
3512 " LDQ({var},(size_t){rows}*{inn},(size_t){rows}*({inn}/32),\"{wdir}/{file}.q\",\"{wdir}/{file}.s\");\n"));
3513 let q4_this = use_q4 && var != "emb" && var != "Wout";
3519 if q4_this {
3520 decls.push_str(&format!(
3521 "static signed char* {var}_q4; static __half* {var}_s4;\n"
3522 ));
3523 loads.push_str(&format!(
3524 " cudaMalloc(&{var}_q4,(size_t){rows}*{inn}/2); cudaMalloc(&{var}_s4,(size_t){rows}*({inn}/32)*2);\n\
3525 \x20 bq_q8_to_q4<<<1024,256>>>((unsigned char*){var}_q4,{var}_s4,{var}_q,{var}_s,{rows},{inn});\n\
3526 \x20 cudaDeviceSynchronize(); cudaFree({var}_q); cudaFree({var}_s); {var}_q=0; {var}_s=0;\n"));
3527 }
3528 };
3529 q8("emb", cfg.vocab, d, "emb");
3530 q8("Wout", cfg.vocab, d, "Wout");
3531 for l in 0..cfg.layers {
3532 if cfg.is_attn(l) {
3533 q8(&format!("qgkv{l}"), attnw, d, &format!("qgkv{l}"));
3534 q8(&format!("Wo{l}"), d, aq, &format!("Wo{l}"));
3535 } else {
3536 q8(&format!("qkvz{l}"), projw, d, &format!("qkvz{l}"));
3537 q8(&format!("Wo{l}"), d, vw, &format!("Wo{l}"));
3538 }
3539 q8(&format!("gu{l}"), 2 * dff, d, &format!("gu{l}"));
3540 q8(&format!("Wd{l}"), d, dff, &format!("Wd{l}"));
3541 }
3542 if mtp {
3555 let l = cfg.layers;
3556 q8(&format!("qgkv{l}"), attnw, d, &format!("qgkv{l}"));
3557 q8(&format!("Wo{l}"), d, aq, &format!("Wo{l}"));
3558 q8(&format!("gu{l}"), 2 * dff, d, &format!("gu{l}"));
3559 q8(&format!("Wd{l}"), d, dff, &format!("Wd{l}"));
3560 q8("ehproj", d, 2 * d, "ehproj");
3561 }
3562 if mtp {
3563 decls.push_str(
3564 "static float *membed,*mcat,*mx,*mtmp,*m2n,*mf32,*dlogits,*mkc,*mvc,*mas,*vlogits;\n\
3565static signed char *maq,*maql;\nstatic int *d_draft,*d_verify,*d_acc;\n",
3566 );
3567 }
3568 let mut f32w = |var: &str, n: usize, file: &str| {
3569 decls.push_str(&format!("static float* {var}_f;\n"));
3570 loads.push_str(&format!(" LDF({var},(size_t){n},\"{wdir}/{file}.f\");\n"));
3571 };
3572 f32w("on", d, "on");
3573 if mtp {
3574 let l = cfg.layers;
3575 f32w(&format!("an{l}"), d, &format!("an{l}"));
3576 f32w(&format!("fn{l}"), d, &format!("fn{l}"));
3577 f32w(&format!("Qn{l}"), cfg.hd, &format!("Qn{l}"));
3578 f32w(&format!("Kn{l}"), cfg.hd, &format!("Kn{l}"));
3579 f32w("enorm", d, "enorm");
3580 f32w("hnorm", d, "hnorm");
3581 f32w("shnorm", d, "shnorm");
3582 }
3583 for l in 0..cfg.layers {
3584 f32w(&format!("an{l}"), d, &format!("an{l}"));
3585 f32w(&format!("fn{l}"), d, &format!("fn{l}"));
3586 if cfg.is_attn(l) {
3587 f32w(&format!("Qn{l}"), cfg.hd, &format!("Qn{l}"));
3588 f32w(&format!("Kn{l}"), cfg.hd, &format!("Kn{l}"));
3589 } else {
3590 f32w(&format!("Wc{l}"), taps * convc, &format!("Wc{l}"));
3591 f32w(&format!("ba{l}"), 2 * nv * d, &format!("ba{l}"));
3592 f32w(&format!("Dt{l}"), nv, &format!("Dt{l}"));
3593 f32w(&format!("Al{l}"), nv, &format!("Al{l}"));
3594 f32w(&format!("Gn{l}"), ghd, &format!("Gn{l}"));
3595 }
3596 }
3597
3598 let mut mtp_roll = String::new();
3599 let mut state_decl = String::new();
3600 let mut state_alloc = String::new();
3601 for l in 0..cfg.layers {
3602 if cfg.is_attn(l) {
3603 state_decl.push_str(&format!("static bq_kv_t *kc{l}, *vc{l};\n"));
3604 state_alloc.push_str(&format!(
3605 " cudaMalloc(&kc{l},(size_t)B*MAXSEQ*{akv}*sizeof(bq_kv_t)); cudaMalloc(&vc{l},(size_t)B*MAXSEQ*{akv}*sizeof(bq_kv_t));\n"));
3606 } else {
3607 if l == 0 {
3608 state_decl
3611 .push_str("static int bq_da[5] = {0,0,0,0,0}; // gdn conv attn norm ffn\n");
3612 }
3613 state_decl.push_str(&format!("static bq_state_t *st{l}; static float *cs{l};\n"));
3614 state_alloc.push_str(&format!(
3615 " cudaMalloc(&st{l},(size_t)B*{vw}*{ghd}*sizeof(bq_state_t)); cudaMemset(st{l},0,(size_t)B*{vw}*{ghd}*sizeof(bq_state_t));\n\
3616 \x20 cudaMalloc(&cs{l},(size_t)B*{taps}*{convc}*4); cudaMemset(cs{l},0,(size_t)B*{taps}*{convc}*4);\n"));
3617 if mtp {
3622 mtp_roll.push_str(&format!(
3623 " bq_mtp_rollback<<<dim3(64,B),256,0,s>>>(st{l},stmid{l},d_acc,(long){vw}*{ghd});\n\
3624 \x20 bq_mtp_rollback<<<dim3(8,B),256,0,s>>>(cs{l},csmid{l},d_acc,(long){taps}*{convc});\n"));
3625 state_decl.push_str(&format!(
3626 "static bq_state_t *stmid{l}; static float *csmid{l};\n"
3627 ));
3628 state_alloc.push_str(&format!(
3629 " cudaMalloc(&stmid{l},(size_t)B*{vw}*{ghd}*sizeof(bq_state_t)); cudaMemset(stmid{l},0,(size_t)B*{vw}*{ghd}*sizeof(bq_state_t));\n\
3630 \x20 cudaMalloc(&csmid{l},(size_t)B*{taps}*{convc}*4); cudaMemset(csmid{l},0,(size_t)B*{taps}*{convc}*4);\n"));
3631 }
3632 }
3633 }
3634
3635 for l in 0..cfg.layers {
3639 for (wv, n, k) in [
3640 (format!("gu{l}"), 2 * dff, d),
3641 (
3643 if cfg.is_attn(l) {
3644 format!("qgkv{l}")
3645 } else {
3646 format!("qkvz{l}")
3647 },
3648 if cfg.is_attn(l) { attnw } else { projw },
3649 d,
3650 ),
3651 (format!("Wd{l}"), d, dff),
3654 ] {
3655 let (dc, ld) = dg_weight(&wv, n, k);
3656 dg_decls.push_str(&dc);
3657 dg_loads.push_str(&ld);
3658 }
3659 }
3660 dg_decls.push_str("static unsigned char* dgA; static float* dgAsf; static float* dgD;\n\
3664 static int bq_dg_hl = 1; // two-pass fp8 weights: accuracy parity with int8\n\
3665 static int bq_dg_on[3] = {1,1,1}; // 0 = gate/up, 1 = entry projection, 2 = down\n");
3666 dg_loads.push_str(&format!(
3672 " cudaMalloc(&dgA,(size_t)BT*{w});\n\
3673 \x20 cudaMalloc(&dgAsf,(size_t)({w}/BQ_DG_BLOCK_K)*bq_dg_align(BT,4)*4);\n\
3674 \x20 cudaMalloc(&dgD,(size_t)BT*{wn}*4);\n",
3675 w = d.max(dff),
3676 wn = (2 * dff).max(projw).max(attnw)
3677 ));
3678
3679 let da = |g: usize| -> String { format!("if(!bq_da[{g}]) ") };
3699 let (da_gdn, da_conv, da_attn, da_norm, da_ffn) = (da(0), da(1), da(2), da(3), da(4));
3700
3701 let mut step = String::new();
3702 for l in 0..cfg.layers {
3703 let wo_sp = if cfg.is_attn(l) {
3704 wo_attn_sp
3705 } else {
3706 wo_gdn_sp
3707 };
3708 let entry = if l == 0 {
3709 "bq_rmsnorm_q8f_hl<<<B,RMSN_NT,0,s>>>(x,aqE,aqEl,asE,x_f32,an0_f,E,EPSV);".to_string()
3710 } else {
3711 if norm_split {
3712 format!("bq_add_rmsnorm_part<<<dim3(NORM_G,B),RMSN_NT,0,s>>>(x,tmp,normp,E,NORM_G);\n\
3713 \x20 bq_add_rmsnorm_apply<<<dim3(NORM_G,B),RMSN_NT,0,s>>>(x,normp,aqE,aqEl,asE,x_f32,an{l}_f,E,EPSV,NORM_G);")
3714 } else {
3715 if wd_sp > 1 {
3716 format!("bq_add_rmsnorm_q8f_hl_sk<<<B,RMSN_NT,0,s>>>(x,skp,{wd_sp},aqE,aqEl,asE,x_f32,an{l}_f,E,EPSV);")
3717 } else {
3718 format!("bq_add_rmsnorm_q8f_hl<<<B,RMSN_NT,0,s>>>(x,tmp,aqE,aqEl,asE,x_f32,an{l}_f,E,EPSV);")
3719 }
3720 }
3721 };
3722 step.push_str(&format!(" {entry}\n"));
3723 if cfg.is_attn(l) {
3724 step.push_str(&format!(
3725 " {g_qgkv}\n\
3726 {ablk}\
3727 \x20 {g_wo}\n",
3728 ablk = format!(
3729 "\x20 {da_attn}{{\n\
3730 \x20 bq_kv_store<<<dim3(Gs({akv}*2),B),256,0,s>>>(kc{l},vc{l},d_pos,{akv},pbuf+{two_aq},PMAX);\n\
3731 \x20 bq_qk_norm_rope_part<<<dim3({h},B),128,0,s>>>(pbuf,d_zero,0,PMAX,Qn{l}_f,{h},{hd},{rot},ROPEV,EPSV);\n\
3732 \x20 bq_qk_norm_rope_kv<<<dim3({kv},B),128,0,s>>>(kc{l},d_pos,{akv},MAXSEQ*{akv},Kn{l}_f,{kv},{hd},{rot},ROPEV,EPSV);\n\
3733 \x20 bq_attn_launch(s,attnb,pbuf,PMAX,kc{l},vc{l},d_pos,{h},{kv},{hd},SCALEV,B);\n\
3734 \x20 bq_gate_mul_q8_hl<<<dim3(Gs({aq}),B),256,0,s>>>(attnb,{aq},pbuf+{aq},PMAX,aqQ,aqQl,asQ,{aq});\n\
3735 \x20 }}\n",
3736 da_attn = da_attn, aq = aq,
3737 akv = akv, l = l, two_aq = 2 * aq, h = cfg.h, kv = cfg.kv, hd = cfg.hd, rot = cfg.rope_dim),
3738 g_qgkv = gemv("pbuf", &format!("qgkv{l}"), "aqE", "asE", d, attnw, projw.max(attnw)),
3739 g_wo = if wo_sp > 1 { drop_reduce(gemv("tmp", &format!("Wo{l}"), "aqQ", "asQ", aq, d, d)) }
3740 else { gemv("tmp", &format!("Wo{l}"), "aqQ", "asQ", aq, d, d) }));
3741 } else {
3742 step.push_str(&format!(
3743 " {g_qkvz}\n\
3744 \x20 bq_gemv_f32w<<<dim3(({two_nv}+(8/F32W_NW)-1)/(8/F32W_NW),B),256,0,s>>>(gbuf,ba{l}_f,x_f32,E,{two_nv},B);\n\
3745 {convk}\
3746 \x20 {gdn}\n\
3747 {gnorm}\
3748 \x20 {g_wo}\n",
3749 convk = format!(
3750 "\x20 {da_conv}{{\n\x20 bq_conv_silu_l2<<<dim3({convc}/{ghd},B),{ghd},0,s>>>(pbuf,PMAX,Wc{l}_f,cs{l},cbuf,{convc},{taps},{two_nk},EPSV);\n\x20 }}\n",
3751 da_conv = da_conv, convc = convc, ghd = ghd, l = l, taps = taps, two_nk = 2 * nk),
3752 gnorm = format!(
3753 "\x20 {da_conv}{{\n\x20 bq_gate_norm_hl<<<dim3({nv},B),{ghd},0,s>>>(pbuf+{qkv_off},Gn{l}_f,obuf,gnbuf,aqQ,aqQl,asQ,{ghd},{nv},PMAX,EPSV);\n\x20 }}\n",
3754 da_conv = da_conv, nv = nv, ghd = ghd, l = l, qkv_off = 2 * qk + vw),
3755 gdn = { let two_qk = 2 * qk;
3756 format!("{da_gdn}{{\n\x20 bq_gdn_gates<<<dim3(Gs({nv}),B),256,0,s>>>(gbuf,Dt{l}_f,Al{l}_f,betab,decayb,{nv});\n\
3757 \x20 bq_gdn_chunk<<<dim3({nv},BQ_GDN_CHUNK_GY({ghd}),B),128,0,s>>>(cbuf,cbuf+{qk},cbuf+{two_qk},betab,decayb,st{l},obuf,{ghd},{nk},{nv},{convc},1,B);\n\x20 }}") },
3758 g_qkvz = gemv("pbuf", &format!("qkvz{l}"), "aqE", "asE", d, projw, projw.max(attnw)),
3759 g_wo = if wo_sp > 1 { drop_reduce(gemv("tmp", &format!("Wo{l}"), "aqQ", "asQ", vw, d, d)) }
3760 else { gemv("tmp", &format!("Wo{l}"), "aqQ", "asQ", vw, d, d) },
3761 two_nv = 2 * nv));
3762 }
3763 step.push_str(&format!(
3764 " {da_norm}{{\n\x20 {norm_ffn}\n\x20 }}\n\
3765 \x20 {da_ffn}{{\n\
3766 \x20 {g_gu}\n\
3767 \x20 {silu}\n\
3768 \x20 {g_wd}\n\
3769 \x20 }}\n\
3770 \x20 {resid}",
3771 da_norm = da_norm, da_ffn = da_ffn,
3772 norm_ffn = if norm_split {
3773 format!("bq_add_rmsnorm_part<<<dim3(NORM_G,B),RMSN_NT,0,s>>>(x,tmp,normp,E,NORM_G);\n\
3774 \x20 bq_add_rmsnorm_apply<<<dim3(NORM_G,B),RMSN_NT,0,s>>>(x,normp,ffQ,ffQl,ffS,x2n,fn{l}_f,E,EPSV,NORM_G);")
3775 } else if wo_sp > 1 {
3776 format!("bq_add_rmsnorm_q8f_hl_sk<<<B,RMSN_NT,0,s>>>(x,skp,{wo_sp},ffQ,ffQl,ffS,x2n,fn{l}_f,E,EPSV);")
3777 } else {
3778 format!("bq_add_rmsnorm_q8f_hl<<<B,RMSN_NT,0,s>>>(x,tmp,ffQ,ffQl,ffS,x2n,fn{l}_f,E,EPSV);")
3779 },
3780 silu = if gu_sp > 1 {
3781 format!("bq_silu_had_q8_hl_sk<<<dim3(Gs({dff}),B),256,0,s>>>(skp,{gu_sp},hqQ,hqQl,hqS,{dff});")
3782 } else {
3783 format!("bq_silu_had_q8_hl<<<dim3(Gs({dff}),B),256,0,s>>>(ffbuf,hqQ,hqQl,hqS,{dff});")
3784 },
3785 g_gu = if gu_sp > 1 { drop_reduce(gemv("ffbuf", &format!("gu{l}"), "ffQ", "ffS", d, 2 * dff, 2 * dff)) }
3786 else { gemv("ffbuf", &format!("gu{l}"), "ffQ", "ffS", d, 2 * dff, 2 * dff) },
3787 g_wd = if wd_sp > 1 && l + 1 < cfg.layers { drop_reduce(gemv("tmp", &format!("Wd{l}"), "hqQ", "hqS", dff, d, d)) }
3788 else { gemv("tmp", &format!("Wd{l}"), "hqQ", "hqS", dff, d, d) },
3789 resid = if l + 1 == cfg.layers { "bq_add<<<dim3(Gs(E),B),256,0,s>>>(x,tmp,E);\n" } else { "\n" }));
3790 }
3791
3792 let pablate = std::env::var("LOGOS_Q38B_PABLATE").unwrap_or_else(|_| "none".to_string());
3803 let (ab_gdn, ab_ffn) = (pablate == "gdn", pablate == "ffn");
3804 let (ab_attn, ab_conv, ab_norm) = (pablate == "attn", pablate == "conv", pablate == "norm");
3808
3809 let mut pstep = String::new();
3810 if prefill_t > 1 {
3811 let bt_rows = batch * prefill_t;
3821 let mi_tcm = 16 * (8 / (64 / (8 * 2))) * 4; let pg = |y: &str, w: &str, xh: &str, xs: &str, k: usize, n: usize, ys: usize| -> String {
3823 if hl || bt_rows < mi_tcm {
3824 return gemv_rows(y, w, xh, xs, k, n, ys, "BT");
3825 }
3826 format!(
3827 "bq_gemm_tcm<<<dim3(({n}+MT_NT-1)/MT_NT,(BT+MI_TCM-1)/MI_TCM),256,0,s>>>\
3828({y},{w}_q,{w}_s,{xh},{xs},{k},{n},BT,{ys});"
3829 )
3830 };
3831 for l in 0..cfg.layers {
3832 let entry = if l == 0 {
3833 "bq_rmsnorm_q8f_hl<<<BT,RMSN_NT,0,s>>>(x,aqE,aqEl,asE,x_f32,an0_f,E,EPSV);"
3834 .to_string()
3835 } else {
3836 if ab_norm {
3837 String::new()
3838 } else {
3839 format!("bq_add_rmsnorm_q8f_hl<<<BT,RMSN_NT,0,s>>>(x,tmp,aqE,aqEl,asE,x_f32,an{l}_f,E,EPSV);")
3840 }
3841 };
3842 pstep.push_str(&format!(" {entry}\n bq_trace_rows(s,x_f32,BT,E);\n"));
3843 if cfg.is_attn(l) {
3844 pstep.push_str(&format!(
3845 " {g}\n\
3846 \x20 {abattn}for(int t=0;t<PT*{attn_on};++t){{ const long po=(long)t*B*PMAX;\n\
3847 \x20 bq_pos_set<<<Gs(B),256,0,s>>>(d_pos,d_pbase,t,B);\n\
3848 \x20 bq_kv_store<<<dim3(Gs({akv}*2),B),256,0,s>>>(kc{l},vc{l},d_pos,{akv},pbuf+po+{two_aq},PMAX);\n\
3849 \x20 bq_qk_norm_rope_part<<<dim3({h},B),128,0,s>>>(pbuf+po,d_zero,0,PMAX,Qn{l}_f,{h},{hd},{rot},ROPEV,EPSV);\n\
3850 \x20 bq_qk_norm_rope_kv<<<dim3({kv},B),128,0,s>>>(kc{l},d_pos,{akv},MAXSEQ*{akv},Kn{l}_f,{kv},{hd},{rot},ROPEV,EPSV);\n\
3851 \x20 bq_attn_launch(s,attnb,pbuf+po,PMAX,kc{l},vc{l},d_pos,{h},{kv},{hd},SCALEV,B);\n\
3852 \x20 bq_gate_mul<<<dim3(Gs({aq}),B),256,0,s>>>(attnb,{aq},pbuf+po+{aq},PMAX,{aq});\n\
3853 \x20 bq_quant_act_hl<<<dim3(Gs({aq}),B),256,0,s>>>(attnb,aqQ+(long)t*B*{aq},aqQl+(long)t*B*{aq},asQ+(long)t*B*({aq}/32),{aq});\n\
3854 \x20 }}\n\
3855 \x20 {gw}\n",
3856 abattn = "", attn_on = if ab_attn { 0 } else { 1 },
3857 g = dg_site(1, "pbuf", &format!("qgkv{l}"), "x_f32", d, attnw, projw.max(attnw),
3858 &pg("pbuf", &format!("qgkv{l}"), "aqE", "asE", d, attnw, projw.max(attnw))),
3859 gw = pg("tmp", &format!("Wo{l}"), "aqQ", "asQ", aq, d, d),
3860 two_aq = 2 * aq, h = cfg.h, kv = cfg.kv, hd = cfg.hd, rot = cfg.rope_dim));
3861 } else {
3862 pstep.push_str(&format!(
3863 " {g}\n\
3864 \x20 bq_gemv_f32w<<<dim3(({two_nv}+(8/F32W_NW)-1)/(8/F32W_NW),BT),256,0,s>>>(gbuf,ba{l}_f,x_f32,E,{two_nv},BT);\n\
3865 \x20 for(int t=0;t<PT*{conv_on};++t)\n\
3866 \x20 bq_conv_silu_l2<<<dim3({convc}/{ghd},B),{ghd},0,s>>>(pbuf+(long)t*B*PMAX,PMAX,Wc{l}_f,cs{l},cbuf+(long)t*B*{convc},{convc},{taps},{two_nk},EPSV);\n\
3867 \x20 {abgates}{gdncall}\n\
3868 \x20 bq_gate_norm_hl<<<dim3({nv},BT),{ghd},0,s>>>(pbuf+{qkv_off},Gn{l}_f,obuf,gnbuf,aqQ,aqQl,asQ,{ghd},{nv},PMAX,EPSV);\n\
3869 \x20 {gw}\n",
3870 g = dg_site(1, "pbuf", &format!("qkvz{l}"), "x_f32", d, projw, projw.max(attnw),
3871 &pg("pbuf", &format!("qkvz{l}"), "aqE", "asE", d, projw, projw.max(attnw))),
3872 gw = pg("tmp", &format!("Wo{l}"), "aqQ", "asQ", vw, d, d),
3873 conv_on = if ab_conv { 0 } else { 1 },
3874 abgates = if ab_gdn { String::new() } else {
3875 format!("bq_gdn_gates<<<dim3(Gs({nv}),BT),256,0,s>>>(gbuf,Dt{l}_f,Al{l}_f,betab,decayb,{nv});\n\x20 ")
3876 },
3877 gdncall = if ab_gdn { String::new() } else if mtp {
3878 format!("bq_gdn_chunk_mid<<<dim3({nv},BQ_GDN_CHUNK_GY({ghd}),B),128,0,s>>>(cbuf,cbuf+{qk},cbuf+{two_qk},betab,decayb,st{l},stmid{l},obuf,{ghd},{nk},{nv},{convc},PT,B,MID);",
3879 nv = nv, ghd = ghd, qk = qk, two_qk = 2 * qk, l = l, nk = nk, convc = convc)
3880 } else {
3881 format!("bq_gdn_chunk<<<dim3({nv},BQ_GDN_CHUNK_GY({ghd}),B),128,0,s>>>(cbuf,cbuf+{qk},cbuf+{two_qk},betab,decayb,st{l},obuf,{ghd},{nk},{nv},{convc},PT,B);",
3882 nv = nv, ghd = ghd, qk = qk, two_qk = 2 * qk, l = l, nk = nk, convc = convc)
3883 },
3884 two_nv = 2 * nv, two_nk = 2 * nk, qkv_off = 2 * qk + vw));
3885 }
3886 pstep.push_str(&format!(
3887 " bq_add_rmsnorm_q8f_hl<<<BT,RMSN_NT,0,s>>>(x,tmp,ffQ,ffQl,ffS,x2n,fn{l}_f,E,EPSV);\n\
3888 \x20 {g_gu}\n\
3889 \x20 {g_wd}\n\
3890 \x20 {resid}",
3891 g_gu = if ab_ffn { String::new() } else { dg_site(
3906 0, "ffbuf", &format!("gu{l}"), "x2n", d, 2 * dff, 2 * dff,
3907 &pg("ffbuf", &format!("gu{l}"), "ffQ", "ffS", d, 2 * dff, 2 * dff)) },
3908 g_wd = if ab_ffn { String::new() } else { dg_site_p(
3909 2, "tmp", &format!("Wd{l}"),
3910 &format!("bq_silu_had_fp8_dg<<<1024,256,0,s>>>(dgA,dgAsf,ffbuf,_m,{dff},_ma);"),
3911 dff, d, d,
3912 &format!("bq_silu_had_q8_hl<<<dim3(Gs({dff}),BT),256,0,s>>>(ffbuf,hqQ,hqQl,hqS,{dff});\n\x20 {}",
3913 pg("tmp", &format!("Wd{l}"), "hqQ", "hqS", dff, d, d))) },
3914 resid = if l + 1 == cfg.layers { "bq_add<<<dim3(Gs(E),BT),256,0,s>>>(x,tmp,E);\n" } else { "\n" }));
3915 }
3916 }
3917
3918 let mtpfn = if !mtp {
3925 String::new()
3926 } else {
3927 let l = cfg.layers;
3928 let two_d = 2 * d;
3929 format!(
3930"static void mtp_draft(cudaStream_t s){{\n\
3931\x20 bq_embed<<<dim3(Gs(E),B),256,0,s>>>(membed,emb_q,emb_s,d_tok,E);\n\
3932\x20 bq_mtp_cat<<<dim3(2,B),256,0,s>>>(mcat,membed,x,enorm_f,hnorm_f,E,EPSV);\n\
3933\x20 bq_quant_act_hl<<<dim3(Gs({two_d}),B),256,0,s>>>(mcat,maq,maql,mas,{two_d});\n\
3934\x20 {g_eh}\n\
3935\x20 bq_rmsnorm_q8f_hl<<<B,RMSN_NT,0,s>>>(mx,maq,maql,mas,mf32,an{l}_f,E,EPSV);\n\
3936\x20 {g_qgkv}\n\
3937\x20 bq_kv_store<<<dim3(Gs({akv}*2),B),256,0,s>>>(mkc,mvc,d_pos,{akv},pbuf+{two_aq},PMAX);\n\
3938\x20 bq_qk_norm_rope_part<<<dim3({h},B),128,0,s>>>(pbuf,d_zero,0,PMAX,Qn{l}_f,{h},{hd},{rot},ROPEV,EPSV);\n\
3939\x20 bq_qk_norm_rope_part<<<dim3({kv},B),128,0,s>>>(mkc,d_pos,{akv},MAXSEQ*{akv},Kn{l}_f,{kv},{hd},{rot},ROPEV,EPSV);\n\
3940\x20 bq_attn_launch(s,attnb,pbuf,PMAX,mkc,mvc,d_pos,{h},{kv},{hd},SCALEV,B);\n\
3941\x20 bq_gate_mul<<<dim3(Gs({aq}),B),256,0,s>>>(attnb,{aq},pbuf+{aq},PMAX,{aq});\n\
3942\x20 bq_quant_act_hl<<<dim3(Gs({aq}),B),256,0,s>>>(attnb,maq,maql,mas,{aq});\n\
3943\x20 {g_wo}\n\
3944\x20 bq_add_rmsnorm_q8f_hl<<<B,RMSN_NT,0,s>>>(mx,mtmp,maq,maql,mas,m2n,fn{l}_f,E,EPSV);\n\
3945\x20 {g_gu}\n\
3946\x20 bq_silu_had_q8_hl<<<dim3(Gs({dff}),B),256,0,s>>>(ffbuf,hqQ,hqQl,hqS,{dff});\n\
3947\x20 {g_wd}\n\
3948\x20 bq_add<<<dim3(Gs(E),B),256,0,s>>>(mx,mtmp,E);\n\
3949\x20 bq_rmsnorm_q8f_hl<<<B,RMSN_NT,0,s>>>(mx,maq,maql,mas,mf32,shnorm_f,E,EPSV);\n\
3950\x20 {g_head}\n\
3951\x20 bq_argmax_rows<<<B,256,0,s>>>(dlogits,VOCAB,d_draft);\n\
3952}}\n\n",
3953 g_eh = gemv("mx", "ehproj", "maq", "mas", two_d, d, d),
3954 g_qgkv = gemv("pbuf", &format!("qgkv{l}"), "maq", "mas", d, attnw, projw.max(attnw)),
3955 g_wo = gemv("mtmp", &format!("Wo{l}"), "maq", "mas", aq, d, d),
3956 g_gu = gemv("ffbuf", &format!("gu{l}"), "maq", "mas", d, 2 * dff, 2 * dff),
3957 g_wd = gemv("mtmp", &format!("Wd{l}"), "hqQ", "hqS", dff, d, d),
3958 g_head = gemv("dlogits", "Wout", "maq", "mas", d, cfg.vocab, cfg.vocab),
3959 two_aq = 2 * aq, h = cfg.h, kv = cfg.kv, hd = cfg.hd, rot = cfg.rope_dim)
3960 + &format!(
3961"static void mtp_verify(cudaStream_t s){{\n\
3962\x20 prefill_chunk(s,1,2); /* a verify is draft + current: two tokens, never PT */\n\
3963\x20 bq_rmsnorm_q8f_hl<<<BT,RMSN_NT,0,s>>>(x,aqE,aqEl,asE,x_f32,on_f,E,EPSV);\n\
3964\x20 {g_vhead}\n\
3965\x20 bq_argmax_rows<<<BT,256,0,s>>>(vlogits,VOCAB,d_verify);\n\
3966}}\n\n",
3967 g_vhead = gemv_rows("vlogits", "Wout", "aqE", "asE", d, cfg.vocab, cfg.vocab, "BT"))
3968 + &format!(
3969"// One speculative step: draft a token, check it against the model in the SAME pass that would have
3970// produced one token anyway, and keep 1 or 2 depending on whether the draft held. Every branch is on
3971// the device -- `bq_mtp_accept` writes the next token and the next position itself -- so this is a
3972// fixed kernel sequence and captures as a graph exactly like decode_step.
3973static void mtp_step(cudaStream_t s){{\n\
3974\x20 cudaMemcpyAsync(d_pbase,d_pos,(size_t)B*4,cudaMemcpyDeviceToDevice,s);\n\
3975\x20 mtp_draft(s);\n\
3976\x20 bq_mtp_pack<<<Gs(B),256,0,s>>>(d_tok,d_draft,B);\n\
3977\x20 mtp_verify(s);\n\
3978\x20 bq_mtp_accept<<<Gs(B),256,0,s>>>(d_verify,d_draft,d_tok,d_pos,d_acc,B);\n\
3979\x20 bq_mtp_pick_hidden<<<dim3(32,B),256,0,s>>>(x,d_acc,B,E);\n\
3980{roll}\
3981}}\n\n", roll = mtp_roll)
3982 };
3983
3984 let mtploop = if !mtp {
3987 String::new()
3988 } else {
3989 String::from(
3990 r#" {
3991 cudaGraph_t mg; cudaGraphExec_t mge;
3992 // Rewind to the start of the cache first. Timing runs 64 speculative steps and each consumes up
3993 // to TWO positions, so starting where generation stopped walks straight off the end of MAXSEQ --
3994 // which is an illegal access inside a graph replay, well away from the line that caused it.
3995 for(int b=0;b<B;++b) hpos[b]=0;
3996 cudaMemcpy(d_pos,hpos,(size_t)B*4,cudaMemcpyHostToDevice);
3997 cudaStreamSynchronize(s);
3998 cudaMemset(d_acc,0,(size_t)(B+1)*4);
3999 mtp_step(s); cudaStreamSynchronize(s);
4000 if(cudaGetLastError()!=cudaSuccess){ fprintf(stderr,"MTP WARMUP: %s\n",cudaGetErrorString(cudaGetLastError())); return 1; }
4001 if(cudaStreamBeginCapture(s,cudaStreamCaptureModeGlobal)!=cudaSuccess){ fprintf(stderr,"mtp capture\n"); return 1; }
4002 mtp_step(s);
4003 if(cudaStreamEndCapture(s,&mg)!=cudaSuccess){ fprintf(stderr,"mtp endcapture\n"); return 1; }
4004 if(cudaGraphInstantiate(&mge,mg,0)!=cudaSuccess){ fprintf(stderr,"mtp inst\n"); return 1; }
4005 const int MSTEPS=64;
4006 cudaMemset(d_acc,0,(size_t)(B+1)*4);
4007 cudaEvent_t m0,m1; cudaEventCreate(&m0); cudaEventCreate(&m1);
4008 cudaEventRecord(m0,s);
4009 for(int i=0;i<MSTEPS;++i) cudaGraphLaunch(mge,s);
4010 cudaEventRecord(m1,s);
4011 if(cudaStreamSynchronize(s)!=cudaSuccess){ fprintf(stderr,"MTP RUN: %s\n",cudaGetErrorString(cudaGetLastError())); return 1; }
4012 float mms=0; cudaEventElapsedTime(&mms,m0,m1);
4013 int nacc=0; cudaMemcpy(&nacc,d_acc+B,4,cudaMemcpyDeviceToHost);
4014 double toks=(double)MSTEPS*B+(double)nacc;
4015 fprintf(stderr,"MTP B %d ACCEPT %.3f TOK_S %.2f (%.3f ms/step, %.2f tok/step/seq)\n",
4016 B, (double)nacc/((double)MSTEPS*B), toks*1000.0/mms, mms/MSTEPS, toks/((double)MSTEPS*B));
4017 }
4018"#,
4019 )
4020 };
4021
4022 let head_launch = if use_wide {
4028 "bq_gemv_wide<<<dim3((VOCAB+MV_RB-1)/MV_RB,1),MV_NW*32,0,s>>>(logits,Wout_q,Wout_s,aqE,aqEl,asE,E,VOCAB,B,VOCAB,(float*)0,1);".to_string()
4029 } else if use_tc {
4030 let sp = tc_split(d, cfg.vocab);
4034 let blocks = cfg.vocab.div_ceil(tc_ti);
4035 let (tck, aqops) = if hl {
4039 ("bq_gemm_tc", "aqE,aqEl")
4040 } else {
4041 ("bq_gemm_tc1", "aqE")
4042 };
4043 if sp > 1 {
4044 tc_part_max.set(tc_part_max.get().max(sp * cfg.vocab));
4045 format!("{tck}<<<dim3({blocks},{sp}),256,0,s>>>(logits,Wout_q,Wout_s,{aqops},asE,E,VOCAB,B,VOCAB,skp,{sp});\n\
4046 \x20 bq_sk_reduce<<<Gs(VOCAB*B),256,0,s>>>(logits,skp,VOCAB,B,VOCAB,{sp});")
4047 } else {
4048 format!("{tck}<<<dim3({blocks},1),256,0,s>>>(logits,Wout_q,Wout_s,{aqops},asE,E,VOCAB,B,VOCAB,(float*)0,1);")
4049 }
4050 } else if use_mma {
4051 "bq_gemm_mma<<<dim3((VOCAB+TBN-1)/TBN,(B+TBM-1)/TBM),256,(size_t)(TBM+TBN)*(KSTEP+4),s>>>(logits,Wout_q,Wout_s,aqE,asE,E,VOCAB,B,VOCAB);".to_string()
4052 } else if use_mmq {
4053 "bq_gemm_mmq<<<(VOCAB+((256/(MMQ_TJ/MMQ_RJ))*MMQ_RI)-1)/((256/(MMQ_TJ/MMQ_RJ))*MMQ_RI),256,0,s>>>(logits,Wout_q,Wout_s,aqE,aqEl,asE,E,VOCAB,B,VOCAB);".to_string()
4054 } else if use_staged && hl {
4055 format!("bq_gemv_hl_smem<<<(VOCAB+7)/8,256,(size_t)(2*{batch}*512+{batch}*16*4),s>>>(logits,Wout_q,Wout_s,aqE,aqEl,asE,E,VOCAB,B,VOCAB);")
4056 } else if use_staged {
4057 format!("bq_gemv_h_smem<<<(VOCAB+7)/8,256,(size_t)({batch}*512+{batch}*16*4),s>>>(logits,Wout_q,Wout_s,aqE,asE,E,VOCAB,B,VOCAB);")
4058 } else if hl {
4059 "bq_gemv_hl<<<dim3(Gw(VOCAB),".to_string()
4060 + &btiles.to_string()
4061 + "),256,0,s>>>(logits,Wout_q,Wout_s,aqE,aqEl,asE,E,VOCAB,B,VOCAB);"
4062 } else {
4063 "bq_gemv_h<<<dim3(Gw(VOCAB),".to_string()
4064 + &btiles.to_string()
4065 + "),256,0,s>>>(logits,Wout_q,Wout_s,aqE,asE,E,VOCAB,B,VOCAB);"
4066 };
4067
4068 let prov = super::provenance::header(super::provenance::Provenance::HandWritten {
4069 model: "Qwen3.8-27B",
4070 lane: "batched decode",
4071 });
4072 format!(
4073"{prov}// Qwen3.8-27B BATCHED decode.
4074// {batch} sequences share ONE pass over the weights; only the recurrence state, conv rings and KV caches
4075// are per-sequence. State cost: {state_mb:.1} MB per sequence at MAXSEQ={maxseq}.
4076#include <cstdio>
4077#include <cstdlib>
4078#include <cstring>
4079#include <cmath>
4080#include <cuda_runtime.h>
4081#include <cuda_fp16.h>
4082#include <chrono>
4083// Only the macros the KERNELS themselves reference may be defined before them. `E` and `B` are parameter
4084// names in the shared kernel set and in the batched GEMVs respectively, so defining them here would rewrite
4085// `int E` into `int 5120` and `int B` into `int 8` — which nvcc reports as a missing parenthesis several
4086// hundred lines later, in an unrelated function.
4087#define BTG {btg}
4088#define BMAX {batch}
4089#define MAXSEQ {maxseq}
4090// MMQ tile: sequence tile rounded up to MMQ_RJ, K staged per barrier round, and the accumulator type.
4091// `double` is measured free on a datacentre part (fp64 half-rate) and 10x more accurate than float;
4092// on a consumer card fp64 is 1/32 rate, which is why LOGOS_Q38B_MMQ_ACC exists.
4093#define MMQ_TJ {mmq_tj}
4094#define MMQ_RI 2
4095#define MMQ_RJ 2
4096#define MMQ_KC 256
4097#define MMQ_ACC {mmq_acc}
4098// Tensor-core tile: sequence tile (a power of two so the warp -> sub-tile map is a bijection),
4099// 16-row sub-tiles per warp, and K staged per barrier round. All three measured, not assumed.
4100#define TC_TJ {tc_tj}
4101{q4_def}#define TC_WM {tc_wm}
4102#define BQ_TC_LANES {tc_lanes}
4103#define TC_KC {tc_kc}
4104{tc_occ_def}
4105#define TC_STAGES {tc_stages}
4106#define PT {prefill_t}
4107#define BT (B*PT)
4108// Rows-in-M tile (prefill): MI sequence rows per block against MT_NT channels.
4109#define MT_WM 4
4110#define MT_NT 64
4111#define MT_WN 2
4112#define MT_KC 64
4113#define MI_TCM (16*(8/(MT_NT/(8*MT_WN)))*MT_WM)
4114// Row-parallel GEMV: sequences in registers, output rows per block, warps reducing K together.
4115#define MV_TJ {mv_tj}
4116#define MV_RB {mv_rb}
4117#define MV_NW {mv_nw}
4118#define RMSN_NT {rmsn_nt}
4119#define NORM_G {norm_g}
4120{gdn_k}
4121{q38_k}
4122{extra_k}
4123{bat_k}
4124#define B {batch}
4125#define PMAX {pmax}
4126#define E {d}
4127#define VOCAB {vocab}
4128#define EPSV {eps:e}f
4129#define ROPEV {rope:e}f
4130#define SCALEV {scale:e}f
4131static float *x,*tmp,*x2n,*x_f32,*pbuf,*cbuf,*obuf,*gnbuf,*gbuf,*betab,*decayb,*attnb,*ffbuf,*logits;
4132static signed char *aqE,*aqQ,*ffQ,*hqQ,*aqEl,*aqQl,*ffQl,*hqQl; static float *asE,*asQ,*ffS,*hqS;
4133static int *d_pos,*d_tok,*d_zero,*d_pbase;\nstatic int hpos[B];
4134static float* skp; static float* normp;
4135{decls}#ifdef BQ_ENABLE_DG
4136{dg_decls}#else
4137static const int bq_dg_on[3] = {{0,0,0}};
4138#endif
4139{state_decl}
4140static inline int Gw(int n){{ return (n*32+255)/256; }}
4141static inline int Gs(int n){{ return (n+255)/256; }}
4142#define LDQ(v,QN,SN,QP,SP) do{{ \\
4143 signed char* hq=(signed char*)malloc((size_t)QN); float* hs=(float*)malloc((size_t)SN*4); \\
4144 FILE* f=fopen(QP,\"rb\"); if(!f||fread(hq,1,(size_t)QN,f)!=(size_t)QN){{fprintf(stderr,\"load %s\\n\",QP);return 1;}} fclose(f); \\
4145 f=fopen(SP,\"rb\"); if(!f||fread(hs,4,(size_t)SN,f)!=(size_t)SN){{fprintf(stderr,\"load %s\\n\",SP);return 1;}} fclose(f); \\
4146 __half* hh=(__half*)malloc((size_t)(SN)*2); for(long i=0;i<(long)(SN);++i) hh[i]=__float2half(hs[i]); \\
4147 cudaMalloc(&v##_q,(size_t)QN); cudaMemcpy(v##_q,hq,(size_t)QN,cudaMemcpyHostToDevice); \\
4148 cudaMalloc(&v##_s,(size_t)(SN)*2); cudaMemcpy(v##_s,hh,(size_t)(SN)*2,cudaMemcpyHostToDevice); \\
4149 free(hq); free(hs); free(hh); }}while(0)
4150#define LDF(v,N,P) do{{ \\
4151 float* h=(float*)malloc((size_t)N*4); FILE* f=fopen(P,\"rb\"); \\
4152 if(!f||fread(h,4,(size_t)N,f)!=(size_t)N){{fprintf(stderr,\"load %s\\n\",P);return 1;}} fclose(f); \\
4153 cudaMalloc(&v##_f,(size_t)N*4); cudaMemcpy(v##_f,h,(size_t)N*4,cudaMemcpyHostToDevice); free(h); }}while(0)
4154
4155{pfn}{mtpfn}static void decode_step(cudaStream_t s){{
4156 bq_embed<<<dim3(Gs(E),B),256,0,s>>>(x,emb_q,emb_s,d_tok,E);
4157{step}\x20 bq_rmsnorm_q8f_hl<<<B,RMSN_NT,0,s>>>(x,aqE,aqEl,asE,x_f32,on_f,E,EPSV);
4158 {head}
4159 bq_argmax_adv<<<B,256,0,s>>>(logits,VOCAB,d_tok,d_pos);
4160}}
4161
4162int main(){{
4163 cudaMalloc(&x,(size_t)BT*E*4); cudaMalloc(&tmp,(size_t)BT*E*4); cudaMalloc(&x2n,(size_t)BT*E*4);
4164 cudaMalloc(&x_f32,(size_t)BT*E*4);
4165 cudaMalloc(&pbuf,(size_t)BT*{pmax}*4); cudaMalloc(&cbuf,(size_t)BT*{convc}*4);
4166 cudaMalloc(&obuf,(size_t)BT*{vw}*4); cudaMalloc(&gnbuf,(size_t)BT*{vw}*4);
4167 cudaMalloc(&gbuf,(size_t)BT*{two_nv}*4); cudaMalloc(&betab,(size_t)BT*{nv}*4); cudaMalloc(&decayb,(size_t)BT*{nv}*4);
4168 cudaMalloc(&attnb,(size_t)B*{aq}*4); cudaMalloc(&ffbuf,(size_t)BT*{two_dff}*4);
4169 cudaMalloc(&logits,(size_t)B*VOCAB*4);
4170{tracopen}
4171{mtpalloc}
4172 cudaMalloc(&aqE,(size_t)BT*E); cudaMalloc(&aqEl,(size_t)BT*E); cudaMalloc(&asE,(size_t)BT*(E/32)*4);
4173 cudaMalloc(&aqQ,(size_t)BT*{amax}); cudaMalloc(&aqQl,(size_t)BT*{amax}); cudaMalloc(&asQ,(size_t)BT*({amax}/32)*4);
4174 cudaMalloc(&ffQ,(size_t)BT*E); cudaMalloc(&ffQl,(size_t)BT*E); cudaMalloc(&ffS,(size_t)BT*(E/32)*4);
4175 cudaMalloc(&hqQ,(size_t)BT*{dff}); cudaMalloc(&hqQl,(size_t)BT*{dff}); cudaMalloc(&hqS,(size_t)BT*({dff}/32)*4);
4176 cudaMalloc(&d_pos,(size_t)B*4); cudaMalloc(&d_tok,(size_t)BT*4);\n // B ints, not one: `bq_qk_norm_rope_part` is launched over dim3(heads,B) and reads d_pos[bb]\n // for every sequence. A 4-byte d_zero was an OUT-OF-BOUNDS read for bb>0 from the moment d_pos\n // became per-sequence -- silent, and dependent on whatever allocation happened to sit next to\n // it, which is why it only changed the answer once BT grew and moved the heap around.\n cudaMalloc(&d_zero,(size_t)B*4);
4177 cudaMalloc(&d_pbase,(size_t)B*4);
4178{skp_alloc}{norm_alloc}
4179{state_alloc}{loads}#ifdef BQ_ENABLE_DG
4180 // Convert the gate/up weights to fp8 e4m3 with per-128 scales and build their TMA descriptors ONCE.
4181 // The Q8_0 copies stay: only the gate/up sites are routed, and the rest of the model still reads them.
4182{dg_loads}#endif
4183
4184 {{ cudaError_t e=cudaGetLastError(); if(e!=cudaSuccess){{
4185 fprintf(stderr,\"CUDA ERROR after weight load: %s\\n\",cudaGetErrorString(e)); return 1; }} }}
4186 {{ for(int b=0;b<B;++b) hpos[b]=0; cudaMemcpy(d_zero,hpos,(size_t)B*4,cudaMemcpyHostToDevice);\n for(int b=0;b<B;++b) hpos[b]=0; cudaMemcpy(d_pos,hpos,(size_t)B*4,cudaMemcpyHostToDevice); }}
4187 cudaStream_t s; cudaStreamCreate(&s);
4188 static const int PROMPT[] = {{{prompt_list}}}; const int NP = sizeof(PROMPT)/sizeof(int);
4189 static int htok[B]; static int htokc[BT];
4190 // ── TTFT ────────────────────────────────────────────────────────────────────────────────────
4191 // Prefill here is the T=1 decode graph replayed ONCE PER PROMPT TOKEN, so a 512-token prompt pays
4192 // 512 whole forwards before the first token exists. vLLM prefills a chunk in one batched pass and
4193 // reports median_ttft_ms per request; this timer is what makes the two comparable -- and what makes
4194 // the deficit a measured number rather than an argument.
4195 auto tt0 = std::chrono::steady_clock::now();
4196 int i0 = 0;
4197#if PT > 1
4198 // Whole chunks first: PT prompt tokens share ONE pass over the weights. Rows are token-major
4199 // (t*B + b), which is why token t's slice is a contiguous base-pointer offset for the recurrence.
4200 // NP-1, not NP: `prefill_chunk` deliberately has no LM head (it is the widest matmul in the
4201 // program and prefill needs the hidden state, not logits), so the LAST prompt token must go through
4202 // `decode_step`, which produces the logits the first generated token comes from AND advances d_pos
4203 // via `bq_argmax_adv`. Consuming the whole prompt in chunks leaves no logits and d_pos one short --
4204 // which is exactly why the first gate run produced a different token stream.
4205 for(; i0 + PT <= NP-1; i0 += PT){{
4206 for(int b=0;b<B;++b) hpos[b]=i0; cudaMemcpy(d_pbase,hpos,(size_t)B*4,cudaMemcpyHostToDevice);
4207 for(int t=0;t<PT;++t) for(int b=0;b<B;++b) htokc[t*B+b]=PROMPT[i0+t];
4208 cudaMemcpy(d_tok,htokc,(size_t)BT*4,cudaMemcpyHostToDevice);
4209 prefill_chunk(s,0,PT);
4210 {{ cudaError_t e=cudaGetLastError(); if(e!=cudaSuccess){{
4211 fprintf(stderr,\"CUDA LAUNCH ERROR at prefill chunk %d: %s\\n\",i0,cudaGetErrorString(e)); return 1; }} }}
4212 if(cudaStreamSynchronize(s)!=cudaSuccess){{
4213 fprintf(stderr,\"CUDA EXEC ERROR at prefill chunk %d: %s\\n\",i0,cudaGetErrorString(cudaGetLastError())); return 1; }}
4214 // the chunk consumed PT positions; decode resumes from there
4215 for(int b=0;b<B;++b) hpos[b]=i0+PT-1; cudaMemcpy(d_pos,hpos,(size_t)B*4,cudaMemcpyHostToDevice);
4216 }}
4217#endif
4218 for(int i=i0;i<NP;++i){{
4219 for(int b=0;b<B;++b) hpos[b]=i; cudaMemcpy(d_pos,hpos,(size_t)B*4,cudaMemcpyHostToDevice);
4220 for(int b=0;b<B;++b) htok[b]=PROMPT[i];
4221 cudaMemcpy(d_tok,htok,(size_t)B*4,cudaMemcpyHostToDevice);
4222 decode_step(s);
4223 {{ cudaError_t e=cudaGetLastError(); if(e!=cudaSuccess){{
4224 fprintf(stderr,\"CUDA LAUNCH ERROR at prefill %d: %s\\n\",i,cudaGetErrorString(e)); return 1; }} }}
4225 if(cudaStreamSynchronize(s)!=cudaSuccess){{
4226 fprintf(stderr,\"CUDA EXEC ERROR at prefill %d: %s\\n\",i,cudaGetErrorString(cudaGetLastError())); return 1; }}
4227 }}
4228 {{ auto tt1 = std::chrono::steady_clock::now();
4229 double pf = std::chrono::duration<double,std::milli>(tt1-tt0).count();
4230 // The last prefill step already produced the first generated token, so time-to-first-token is the
4231 // prefill span plus the copy that makes it visible to the host -- measured, not modelled.
4232 static int hfirst[B]; cudaMemcpy(hfirst,d_tok,(size_t)B*4,cudaMemcpyDeviceToHost);
4233 auto tt2 = std::chrono::steady_clock::now();
4234 double ttft = std::chrono::duration<double,std::milli>(tt2-tt0).count();
4235 fprintf(stderr,\"PREFILL_MS %.3f (%d prompt tokens, %.3f ms/token)\\n\", pf, NP, pf/(NP?NP:1));
4236 fprintf(stderr,\"TTFT_MS %.3f (batch %d, first token %d)\\n\", ttft, B, hfirst[0]); }}
4237 {{ // the logits the prompt produced, before anything is generated -- a sanity read, not a metric
4238 static float hl0[16]; cudaMemcpy(hl0,logits,16*4,cudaMemcpyDeviceToHost);
4239 fprintf(stderr,\"PREFILL_LOGITS[0..3] %.4f %.4f %.4f %.4f\\n\",hl0[0],hl0[1],hl0[2],hl0[3]); }}
4240 cudaGraph_t g; cudaGraphExec_t ge;
4241 if(cudaStreamBeginCapture(s,cudaStreamCaptureModeGlobal)!=cudaSuccess){{fprintf(stderr,\"capture\\n\");return 1;}}
4242 decode_step(s);
4243 if(cudaStreamEndCapture(s,&g)!=cudaSuccess){{fprintf(stderr,\"endcapture\\n\");return 1;}}
4244 if(cudaGraphInstantiate(&ge,g,0)!=cudaSuccess){{fprintf(stderr,\"inst\\n\");return 1;}}
4245 cudaEvent_t t0,t1; cudaEventCreate(&t0); cudaEventCreate(&t1);
4246 cudaMemcpy(htok,d_tok,(size_t)B*4,cudaMemcpyDeviceToHost);
4247{mtptok} printf(\"SEQ0_TOKENS: %d\", htok[0]);
4248 static int rows[B][{n_gen}]; for(int b=0;b<B;++b) rows[b][0]=htok[b];
4249 for(int i=1;i<{n_gen};++i){{
4250 cudaGraphLaunch(ge,s); cudaStreamSynchronize(s);
4251 cudaMemcpy(htok,d_tok,(size_t)B*4,cudaMemcpyDeviceToHost);
4252 printf(\" %d\", htok[0]);
4253 for(int b=0;b<B;++b) rows[b][i]=htok[b];
4254 }}
4255 printf(\"\\n\");
4256 int diverged=0;
4257 for(int b=1;b<B;++b) for(int i=0;i<{n_gen};++i) if(rows[b][i]!=rows[0][i]) diverged++;
4258 printf(\"BATCH_AGREEMENT: %d/%d sequences identical to sequence 0 (%d token mismatches)\\n\",
4259 diverged?0:B, B, diverged);
4260{mtploop} {{ float ms=0; for(int b=0;b<B;++b) hpos[b]=0;\n cudaMemcpy(d_pos,hpos,(size_t)B*4,cudaMemcpyHostToDevice); cudaStreamSynchronize(s);
4261 if(getenv(\"LOGOS_B_BISECT\")){{
4262 for(int i=0;i<128;++i){{ cudaGraphLaunch(ge,s);
4263 if(cudaStreamSynchronize(s)!=cudaSuccess){{ int hp=0; cudaMemcpy(&hp,d_pos,4,cudaMemcpyDeviceToHost);
4264 fprintf(stderr,\"FAULT at replay %d (d_pos=%d, MAXSEQ=%d): %s\\n\",i,hp,MAXSEQ,
4265 cudaGetErrorString(cudaGetLastError())); return 1; }} }}
4266 fprintf(stderr,\"BISECT: all 128 replays clean\\n\"); }}
4267 /* ★ THE DISTRIBUTION, NOT JUST THE MEAN. vLLM reports median_itl_ms and p99_itl_ms, and its
4268 median ITL IS its step time — which is the fact that corrected this campaign's biggest
4269 misreading: comparing our pure-GPU decode throughput against vLLM's full serving number
4270 flattered us, and only the ITL comparison made the real 2.13x gap visible. A mean over a
4271 block cannot be compared to a median, and it hides the tail that serving is actually judged
4272 on. One event per step costs ~1 us against a ~33 ms step. */
4273 cudaEvent_t iev[129]; for(int i=0;i<=128;++i) cudaEventCreate(&iev[i]);
4274 cudaEventRecord(t0,s); cudaEventRecord(iev[0],s);
4275 for(int i=0;i<128;++i){{ cudaGraphLaunch(ge,s); cudaEventRecord(iev[i+1],s); }}
4276 cudaEventRecord(t1,s); cudaStreamSynchronize(s);
4277 cudaError_t te=cudaEventElapsedTime(&ms,t0,t1);
4278 if(te!=cudaSuccess||ms<=0.f){{ fprintf(stderr,\"TIMING FAILED: %s (ms=%.3f)\\n\",cudaGetErrorString(te),ms); return 1; }}
4279 fprintf(stderr,\"BATCH %d PER_SEQ_TOK_S %.2f AGGREGATE_TOK_S %.2f (%.3f ms/step)\\n\",
4280 B, 128000.0/ms, 128000.0*B/ms, ms/128.0);
4281 {{ float itl[128]; int ok=1;
4282 for(int i=0;i<128;++i) if(cudaEventElapsedTime(&itl[i],iev[i],iev[i+1])!=cudaSuccess) ok=0;
4283 for(int i=0;i<=128;++i) cudaEventDestroy(iev[i]);
4284 if(ok){{
4285 for(int i=1;i<128;++i){{ float k=itl[i]; int j=i-1;
4286 while(j>=0&&itl[j]>k){{ itl[j+1]=itl[j]; --j; }} itl[j+1]=k; }}
4287 float mean=0.f; for(int i=0;i<128;++i) mean+=itl[i]; mean/=128.f;
4288 fprintf(stderr,\"ITL_MEAN_MS %.3f ITL_MEDIAN_MS %.3f ITL_P99_MS %.3f ITL_MIN_MS %.3f\\n\",
4289 mean, 0.5f*(itl[63]+itl[64]), itl[126], itl[0]);
4290 }} }} }}
4291 return 0;
4292}}
4293",
4294 btg = BTG, vocab = cfg.vocab, eps = cfg.rms_eps, rope = cfg.rope_theta,
4295 scale = 1.0 / (cfg.hd as f64).sqrt(),
4296 gdn_k = crate::codegen_cuda::gdn::GDN_KERNELS,
4297 q38_k = crate::codegen_cuda::gdn::QWEN38_KERNELS,
4298 extra_k = crate::codegen_cuda::qwen38::QWEN38_EXTRA_KERNELS,
4299 bat_k = QWEN38_BATCHED_KERNELS,
4300 mtpfn = mtpfn, mtploop = mtploop,
4306 mtptok = if !mtp { String::new() } else { String::from(r#" if(getenv("LOGOS_MTP_TOKENS")){
4311 printf("MTP_TOKENS: %d", htok[0]);
4312 int produced=1, steps=0, acc=0; static int hv[2*BMAX], ha[BMAX+1];
4313 cudaEvent_t k0,k1; cudaEventCreate(&k0); cudaEventCreate(&k1); cudaEventRecord(k0,s);
4314 while(produced < 128){
4315 mtp_step(s);
4316 if(cudaStreamSynchronize(s)!=cudaSuccess){ fprintf(stderr,"MTP TOKENS: %s\\n",cudaGetErrorString(cudaGetLastError())); return 1; }
4317 cudaMemcpy(hv,d_verify,(size_t)2*B*4,cudaMemcpyDeviceToHost);
4318 cudaMemcpy(ha,d_acc,(size_t)B*4,cudaMemcpyDeviceToHost);
4319 ++steps; acc += ha[0];
4320 printf(" %d", hv[0]); ++produced;
4321 if(ha[0] && produced < 128){ printf(" %d", hv[B]); ++produced; }
4322 }
4323 printf("\\n");
4324 cudaEventRecord(k1,s); cudaStreamSynchronize(s);
4325 float kms=0; cudaEventElapsedTime(&kms,k0,k1);
4326 // Timed from the REAL post-prefill state, so the accept rate is the one this prompt
4327 // actually produces. A host sync per step is included, making the figure CONSERVATIVE.
4328 fprintf(stderr,"MTP_REAL B %d ACCEPT %.3f TOK_S %.2f (%d steps, %d tokens, %.3f ms/step)\\n",
4329 B, (double)acc/steps, (double)produced*B*1000.0/kms, steps, produced, kms/steps);
4330 return 0;
4331 }
4332"#) },
4333 mtpalloc = if !mtp { String::new() } else {
4334 let mscratch = (2 * d).max(dff).max(aq);
4335 format!(
4336"\x20 cudaMalloc(&membed,(size_t)B*E*4); cudaMalloc(&mcat,(size_t)B*{two_d}*4);\n\
4337\x20 cudaMalloc(&mx,(size_t)B*E*4); cudaMalloc(&mtmp,(size_t)B*E*4);\n\
4338\x20 cudaMalloc(&m2n,(size_t)B*E*4); cudaMalloc(&mf32,(size_t)B*E*4);\n\
4339\x20 cudaMalloc(&dlogits,(size_t)B*VOCAB*4); cudaMalloc(&d_draft,(size_t)B*4);\n\
4340\x20 cudaMalloc(&mkc,(size_t)B*MAXSEQ*{akv}*4); cudaMalloc(&mvc,(size_t)B*MAXSEQ*{akv}*4);\n\
4341\x20 cudaMalloc(&maq,(size_t)B*{ms}); cudaMalloc(&maql,(size_t)B*{ms});\n\
4342\x20 cudaMalloc(&mas,(size_t)B*({ms}/32)*4);\n\
4343\x20 cudaMalloc(&vlogits,(size_t)BT*VOCAB*4); cudaMalloc(&d_verify,(size_t)BT*4);\n\
4344\x20 cudaMalloc(&d_acc,(size_t)(B+1)*4); cudaMemset(d_acc,0,(size_t)(B+1)*4);\n",
4345 two_d = 2 * d, ms = mscratch)
4346 },
4347 pfn = if prefill_t > 1 {
4348 format!("static FILE* bq_tf=nullptr; static float* bq_th=nullptr; static int bq_tchunk=0;\n\
4369static void bq_trace(cudaStream_t s,const float* p,long off,long n){{\n\
4370\x20 if(!bq_tf||bq_tchunk!=1) return;\n\
4371\x20 cudaStreamSynchronize(s);\n\
4372\x20 cudaMemcpy(bq_th,p+off,(size_t)n*4,cudaMemcpyDeviceToHost);\n\
4373\x20 fwrite(bq_th,4,(size_t)n,bq_tf);\n\
4374}}\n\
4375static void bq_trace_rows(cudaStream_t s,const float* p,int rows,int w){{\n\
4376\x20 if(!bq_tf||bq_tchunk!=1) return;\n\
4377\x20 const int r = rows<4?rows:4;\n\
4378\x20 bq_trace(s,p,0,(long)r*w);\n\
4379\x20 bq_trace(s,p,(long)(rows-r)*w,(long)r*w);\n\
4380}}\n\
4381static void prefill_chunk(cudaStream_t s, int MID, int PTv){{\n\
4382 \x20 const int BTv = B*PTv; ++bq_tchunk;\n\
4383 \x20 bq_embed<<<dim3(Gs(E),BTv),256,0,s>>>(x,emb_q,emb_s,d_tok,E);\n\
4384 {pstep}}}\n\n",
4385 pstep = replace_word(&replace_word(&pstep, "BT", "BTv"), "PT", "PTv"))
4386 } else { String::new() },
4387 tracopen = {
4388 let dab = concat!(" { const char* _db=getenv(\"LOGOS_Q38B_DABLATE\");\n",
4389 " if(_db){ const char* g[5]={\"gdn\",\"conv\",\"attn\",\"norm\",\"ffn\"};\n",
4390 " for(int i=0;i<5;++i) if(!strcmp(_db,g[i])) bq_da[i]=1; }\n",
4391 " fprintf(stderr,\"DABLATE %s\\n\", _db?_db:\"none\"); }\n");
4392 let rest = if prefill_t > 1 {
4393 concat!("#ifdef BQ_ENABLE_DG\n",
4394 " { const char* _ds=getenv(\"LOGOS_Q38B_DGSITES\");\n",
4395 " if(_ds) for(int i=0;i<3&&_ds[i];++i) bq_dg_on[i] = (_ds[i]=='1');\n",
4396 " const char* _dh=getenv(\"LOGOS_Q38B_DGHL\"); if(_dh) bq_dg_hl = (_dh[0]!='0');\n",
4397 " fprintf(stderr,\"DGHL %d\\n\",bq_dg_hl);\n",
4398 " fprintf(stderr,\"DGSITES %d%d%d\\n\",bq_dg_on[0],bq_dg_on[1],bq_dg_on[2]); }\n",
4399 "#endif\n",
4400 " { const char* _tp=getenv(\"LOGOS_Q38B_TRACE\"); if(_tp){ bq_tf=fopen(_tp,\"wb\"); bq_th=(float*)malloc((size_t)8*E*4); } }\n")
4401 } else { "" };
4402 format!("{dab}{rest}")
4403 },
4404 skp_alloc = if tc_part_max.get() > 0 {
4405 format!(" cudaMalloc(&skp,(size_t){}*BT*4);\n", tc_part_max.get())
4407 } else { String::new() },
4408 tc_tj = tc_tj, tc_wm = tc_wm, tc_kc = tc_kc, tc_occ_def = tc_occ_def,
4409 q4_def = if use_q4 { "#define BQ_ENABLE_Q4 1\n" } else { "" },
4410 tc_lanes = if hl { 2 } else { 1 }, tc_stages = tc_stages, prefill_t = prefill_t,
4411 mv_tj = mv_tj, mv_rb = mv_rb, mv_nw = mv_nw, rmsn_nt = rmsn_nt, norm_g = norm_g,
4412 norm_alloc = if norm_split {
4413 " cudaMalloc(&normp,(size_t)NORM_G*B*4);\n".to_string()
4415 } else { String::new() },
4416 mmq_tj = mmq_tj,
4417 mmq_acc = mmq_acc,
4418 head = head_launch,
4419 pmax = projw.max(attnw), two_dff = 2 * dff, two_nv = 2 * nv,
4420 state_mb = cfg.state_bytes_per_seq(maxseq) as f64 / 1e6,
4421 prompt_list = prompt.iter().map(|t| t.to_string()).collect::<Vec<_>>().join(","),
4422 )
4423}