spark_model/layers/glm5next_mlp/
forward.rs

1// SPDX-License-Identifier: AGPL-3.0-only
2
3//! The GLM MLP decode forward β€” dense FFN and routed MoE, one token.
4//!
5//! Launch geometry is lifted verbatim from the two gated microtests
6//! (`examples/glm5next_{ffn,moe}_microtest.rs`, Slice-10 gates 3/4/6/7), which measured this
7//! exact sequence against HF 5.16.1 on real layer-0 and layer-3 weights. Nothing here
8//! re-derives the equations.
9
10use anyhow::{Result, bail};
11use spark_runtime::gpu::{DevicePtr, GpuBackend, KernelHandle};
12use spark_runtime::kernel_args::KernelLaunch;
13
14use super::weights::{Glm5NextDenseMlpWeights, Glm5NextMoeWeights, Nvfp4Proj};
15use super::{Glm5NextMlpConfig, Glm5NextMlpKernels};
16
17const W4_TILE: u32 = 64;
18const ACT_BLOCK: u32 = 256;
19
20/// `C[M, N] = A[M, K] @ B[N, K]^T`, BF16 in and out.
21#[allow(clippy::too_many_arguments)]
22fn gemm(
23    gpu: &dyn GpuBackend,
24    k: KernelHandle,
25    gemv: KernelHandle,
26    batchm: KernelHandle,
27    a: DevicePtr,
28    b: DevicePtr,
29    c: DevicePtr,
30    m: usize,
31    n: usize,
32    kk: usize,
33    stream: u64,
34) -> Result<()> {
35    // M=1 decode -> GEMV; M=2..8 (a K-token verify sweep) -> ONE weight read for all rows;
36    // wider -> the tile GEMM. `ops::dense_mm_bf16` owns the policy and the grid coupling.
37    //
38    // πŸ”΄ The router is the worst tile-GEMM case in the whole stack: N = 288 tiles to **18
39    // blocks**, measured 6.8 GB/s. It has no FP32-out batchm twin, so it stays on gemv/tile.
40    crate::layers::ops::dense_mm_bf16(
41        gpu,
42        &crate::layers::ops::DenseMmKernels {
43            gemm: k,
44            gemv,
45            batchm,
46        },
47        a,
48        b,
49        c,
50        m,
51        n,
52        kk,
53        stream,
54    )
55}
56
57/// `C[1, N] = A[1, K] @ dequant(B)[N, K]^T` β€” the M=1 decode kernel.
58///
59/// Same NVFP4 operand triple as [`w4a16`], one output row. Used for every routed-expert
60/// projection because they are all M=1 and the tile GEMM measured 9.7 GB/s there.
61fn w4a16_gemv(
62    gpu: &dyn GpuBackend,
63    k: KernelHandle,
64    k_sw: KernelHandle,
65    a: DevicePtr,
66    w: &Nvfp4Proj,
67    c: DevicePtr,
68    n: usize,
69    kk: usize,
70    stream: u64,
71) -> Result<()> {
72    // Prefer the single-warp sibling when the target carries it. BIT-IDENTICAL β€” the
73    // per-orig-lane partial is the same function, the shuffle tree is the same tree, and
74    // the final combine is the same two-term FP32 add; only the block packing differs.
75    // This site launched the base kernel directly and so had never picked up the SW win
76    // that `ops::w4a16_decode_gemv` has been handing every other decode GEMV.
77    if k_sw.0 != 0 {
78        return crate::layers::ops::w4a16_gemv_sw_raw(
79            gpu, k_sw, a, w.packed, w.scale, w.scale_2, c, n as u32, kk as u32, stream,
80        );
81    }
82    KernelLaunch::new(gpu, k)
83        .grid([crate::layers::ops::w4a16_gemv_grid_x(n as u32), 1, 1])
84        .block([256, 1, 1])
85        .arg_ptr(a)
86        .arg_ptr(w.packed)
87        .arg_ptr(w.scale)
88        // πŸͺ€ by VALUE, as in `w4a16`.
89        .arg_f32(w.scale_2)
90        .arg_ptr(c)
91        .arg_u32(n as u32)
92        .arg_u32(kk as u32)
93        .launch(stream)?;
94    Ok(())
95}
96
97/// `C[M, N] = A[M, K] @ dequant(B)[N, K]^T` β€” NVFP4 weight, BF16 activation and output.
98#[allow(clippy::too_many_arguments)]
99#[allow(dead_code)]
100fn w4a16(
101    gpu: &dyn GpuBackend,
102    k: KernelHandle,
103    a: DevicePtr,
104    w: &Nvfp4Proj,
105    c: DevicePtr,
106    m: usize,
107    n: usize,
108    kk: usize,
109    stream: u64,
110) -> Result<()> {
111    KernelLaunch::new(gpu, k)
112        .grid([
113            (n as u32).div_ceil(W4_TILE),
114            (m as u32).div_ceil(W4_TILE),
115            1,
116        ])
117        .block([128, 1, 1])
118        .arg_ptr(a)
119        .arg_ptr(w.packed)
120        .arg_ptr(w.scale)
121        // πŸͺ€ by VALUE. `weight_scale_2` is a scalar argument, not a pointer.
122        .arg_f32(w.scale_2)
123        .arg_ptr(c)
124        .arg_u32(m as u32)
125        .arg_u32(n as u32)
126        .arg_u32(kk as u32)
127        .launch(stream)?;
128    Ok(())
129}
130
131/// `out = silu(min(gate, limit)) * clamp(up, -limit, limit)` over `n` elements.
132fn swiglu(
133    gpu: &dyn GpuBackend,
134    k: KernelHandle,
135    gate: DevicePtr,
136    up: DevicePtr,
137    out: DevicePtr,
138    n: usize,
139    limit: f32,
140    stream: u64,
141) -> Result<()> {
142    KernelLaunch::new(gpu, k)
143        .grid([(n as u32).div_ceil(ACT_BLOCK), 1, 1])
144        .block([ACT_BLOCK, 1, 1])
145        .arg_ptr(gate)
146        .arg_ptr(up)
147        .arg_ptr(out)
148        .arg_u32(n as u32)
149        .arg_f32(limit)
150        .launch(stream)?;
151    Ok(())
152}
153
154/// Scratch for one MLP site, allocated once and reused every decode step.
155///
156/// Sized for the widest thing this site can run: a dense layer's `intermediate_size`, a routed
157/// layer's `moe_intermediate_size`, and the shared expert's width.
158pub struct Glm5NextMlpWorkspace {
159    /// `[rows, max_inter]` BF16 Γ—3 β€” gate, up, activated. Shared by every projection pair.
160    a_gate: DevicePtr,
161    a_up: DevicePtr,
162    a_act: DevicePtr,
163    /// `[rows, num_experts]` F32 router logits.
164    logits: DevicePtr,
165    /// `[rows, top_k]` I32 / F32 selection.
166    ids: DevicePtr,
167    wts: DevicePtr,
168    /// `[rows, top_k, hidden]` BF16. πŸͺ€ Fully written every step β€” remote and invalid slots are
169    /// memset to zero before the loop, never left stale.
170    expert_out: DevicePtr,
171    /// `[rows, hidden]` BF16 shared-expert output.
172    shared_out: DevicePtr,
173    /// `[rows * top_k]` I32 union expert ids, `-1` = entry unused. Row-batched MoE only.
174    u_eid: DevicePtr,
175    /// `[rows * top_k, rows]` I32 slot per union entry per row, `-1` = row absent.
176    u_slot: DevicePtr,
177    max_inter: usize,
178    /// Widest verify this scratch can serve. `1` on the serial decode path.
179    max_rows: usize,
180}
181
182impl Glm5NextMlpWorkspace {
183    pub fn new(gpu: &dyn GpuBackend, cfg: &Glm5NextMlpConfig, max_rows: usize) -> Result<Self> {
184        let rows = max_rows.max(1);
185        let max_inter = cfg
186            .local_dense_intermediate
187            .max(cfg.moe_intermediate)
188            .max(cfg.local_shared_intermediate)
189            .max(1);
190        // The grouped MoE path activates all `top_k` slots in one launch, so the three
191        // activation buffers are slot-major and `top_k` times a routed expert's width. The
192        // dense path still writes only the first `inter` elements β€” `max_inter` stays the
193        // guard for it.
194        // The dense/shared arm is now `[rows, inter]`; the grouped MoE arm is still one row's
195        // `top_k` slots at a time. Both share these buffers, so take the wider.
196        // The row-batched MoE arm computes every (row, slot) pair in ONE launch, so its
197        // activations are `[rows, top_k, moe_intermediate]` β€” wider than either of the above.
198        let act_elems = (rows * max_inter)
199            .max(rows * cfg.top_k * cfg.moe_intermediate)
200            .max(1);
201        Ok(Self {
202            a_gate: gpu.alloc(act_elems * 2)?,
203            a_up: gpu.alloc(act_elems * 2)?,
204            a_act: gpu.alloc(act_elems * 2)?,
205            logits: gpu.alloc(rows * cfg.num_experts * 4)?,
206            ids: gpu.alloc(rows * cfg.top_k * 4)?,
207            wts: gpu.alloc(rows * cfg.top_k * 4)?,
208            expert_out: gpu.alloc(rows * cfg.top_k * cfg.hidden * 2)?,
209            shared_out: gpu.alloc(rows * cfg.hidden * 2)?,
210            u_eid: gpu.alloc(rows * cfg.top_k * 4)?,
211            u_slot: gpu.alloc(rows * cfg.top_k * rows * 4)?,
212            max_inter,
213            max_rows: rows,
214        })
215    }
216}
217
218/// A BF16 SwiGLU MLP of width `inter`: `down(clamped_swiglu(gate(x), up(x)))`.
219///
220/// Used for both the dense layers and the shared expert β€” identical math, different widths.
221/// With `tp_world_size > 1` the result is a **partial sum**; the caller reduces.
222#[allow(clippy::too_many_arguments)]
223pub fn forward_dense(
224    gpu: &dyn GpuBackend,
225    k: &Glm5NextMlpKernels,
226    cfg: &Glm5NextMlpConfig,
227    w: &Glm5NextDenseMlpWeights,
228    inter: usize,
229    x: DevicePtr,
230    out: DevicePtr,
231    m: usize,
232    ws: &Glm5NextMlpWorkspace,
233    stream: u64,
234) -> Result<()> {
235    if inter == 0 || inter > ws.max_inter {
236        bail!(
237            "GLM dense MLP: width {inter} does not fit a workspace built for {}",
238            ws.max_inter
239        );
240    }
241    if m == 0 || m > ws.max_rows {
242        bail!(
243            "GLM dense MLP: {m} rows do not fit a workspace built for {}",
244            ws.max_rows
245        );
246    }
247    gemm(
248        gpu,
249        k.gemm,
250        k.gemv,
251        k.gemv_batchm,
252        x,
253        w.gate_proj,
254        ws.a_gate,
255        m,
256        inter,
257        cfg.hidden,
258        stream,
259    )?;
260    gemm(
261        gpu,
262        k.gemm,
263        k.gemv,
264        k.gemv_batchm,
265        x,
266        w.up_proj,
267        ws.a_up,
268        m,
269        inter,
270        cfg.hidden,
271        stream,
272    )?;
273    swiglu(
274        gpu,
275        k.swiglu,
276        ws.a_gate,
277        ws.a_up,
278        ws.a_act,
279        // Elementwise over the whole `[m, inter]` block β€” bit-identical to m launches of
280        // `inter`, because every output element depends only on its own gate/up pair.
281        m * inter,
282        cfg.swiglu_limit,
283        stream,
284    )?;
285    gemm(
286        gpu,
287        k.gemm,
288        k.gemv,
289        k.gemv_batchm,
290        ws.a_act,
291        w.down_proj,
292        out,
293        m,
294        cfg.hidden,
295        inter,
296        stream,
297    )
298}
299
300/// `C[top_k, N] = A @ dequant(expert[ids[slot]])^T` β€” every routed slot in ONE launch.
301///
302/// Bit-identical per slot to the [`w4a16_gemv`] loop it replaces; see the kernel's header.
303/// Slots whose expert this rank does not own are skipped, so `c` must already hold whatever
304/// those rows should contribute (zero, for the routed sum).
305#[allow(clippy::too_many_arguments)]
306fn w4a16_gemv_moe(
307    gpu: &dyn GpuBackend,
308    k: KernelHandle,
309    a: DevicePtr,
310    t: &super::weights::Glm5NextExpertPtrTable,
311    c: DevicePtr,
312    ids: DevicePtr,
313    n: usize,
314    kk: usize,
315    top_k: usize,
316    num_experts: usize,
317    input_stride: usize,
318    stream: u64,
319) -> Result<()> {
320    KernelLaunch::new(gpu, k)
321        // πŸͺ€ COUPLED to the kernel's `N_PER_BLOCK_SW` = 8, and grid.y IS the slot.
322        .grid([
323            crate::layers::ops::w4a16_gemv_sw_grid_x(n as u32),
324            top_k as u32,
325            1,
326        ])
327        .block([256, 1, 1])
328        .arg_ptr(a)
329        .arg_ptr(t.packed_ptrs)
330        .arg_ptr(t.scale_ptrs)
331        .arg_ptr(t.scale2_vals)
332        .arg_ptr(c)
333        .arg_ptr(ids)
334        .arg_u32(n as u32)
335        .arg_u32(kk as u32)
336        .arg_u32(num_experts as u32)
337        .arg_u32(input_stride as u32)
338        .launch(stream)
339}
340
341/// `C[rows, top_k, N]` β€” the UNION of the rows' selected experts, each expert swept ONCE.
342///
343/// Bit-identical per (row, slot) to [`w4a16_gemv_moe`]; see the kernel header. `u_eid` /
344/// `u_slot` come from `glm5next_moe_row_union` and stay on device, so this is capturable.
345///
346/// πŸͺ€ grid.y is `rows * top_k` β€” the UNION extent, not `top_k`. Entries the routing did not
347/// fill retire immediately on `u_eid < 0`.
348#[allow(clippy::too_many_arguments)]
349fn w4a16_gemv_moe_batchm(
350    gpu: &dyn GpuBackend,
351    k: KernelHandle,
352    a: DevicePtr,
353    t: &super::weights::Glm5NextExpertPtrTable,
354    c: DevicePtr,
355    u_eid: DevicePtr,
356    u_slot: DevicePtr,
357    n: usize,
358    kk: usize,
359    rows: usize,
360    top_k: usize,
361    num_experts: usize,
362    a_row_stride: usize,
363    a_slot_stride: usize,
364    c_row_stride: usize,
365    stream: u64,
366) -> Result<()> {
367    KernelLaunch::new(gpu, k)
368        // πŸͺ€ COUPLED to the kernel's `N_PER_BLOCK_SW` = 8.
369        .grid([
370            crate::layers::ops::w4a16_gemv_sw_grid_x(n as u32),
371            (rows * top_k) as u32,
372            1,
373        ])
374        .block([256, 1, 1])
375        .arg_ptr(a)
376        .arg_ptr(t.packed_ptrs)
377        .arg_ptr(t.scale_ptrs)
378        .arg_ptr(t.scale2_vals)
379        .arg_ptr(c)
380        .arg_ptr(u_eid)
381        .arg_ptr(u_slot)
382        .arg_u32(n as u32)
383        .arg_u32(kk as u32)
384        .arg_u32(num_experts as u32)
385        .arg_u32(a_row_stride as u32)
386        .arg_u32(a_slot_stride as u32)
387        .arg_u32(c_row_stride as u32)
388        .launch(stream)
389}
390
391/// Kill switch for the row-batched routed path: `ATLAS_NO_GLM_MOE_ROW_BATCH=1` restores the
392/// one-launch-per-row grouped dispatch. Read once β€” this sits on the per-layer decode path.
393fn row_batch_disabled() -> bool {
394    static F: std::sync::OnceLock<bool> = std::sync::OnceLock::new();
395    *F.get_or_init(|| std::env::var("ATLAS_NO_GLM_MOE_ROW_BATCH").as_deref() == Ok("1"))
396}
397
398/// Widest tier `w4a16_gemv_sw_moe_batchm` may be dispatched at, `ATLAS_GLM_MOE_ROW_BATCH_MAX`.
399///
400/// πŸ”¬ An A/B lever, not a tuning knob: `=4` restores the pre-2026-08-31 cap exactly, so the
401/// width extension can be measured against itself in ONE image instead of one image per arm β€”
402/// the same lever that found A65's real defect. Clamped to the compiled tier family.
403pub(crate) fn row_batch_max() -> usize {
404    static M: std::sync::OnceLock<usize> = std::sync::OnceLock::new();
405    *M.get_or_init(|| {
406        let m = std::env::var("ATLAS_GLM_MOE_ROW_BATCH_MAX")
407            .ok()
408            .and_then(|v| v.parse::<usize>().ok())
409            .unwrap_or(MOE_ROW_BATCH_MAX_ROWS)
410            .clamp(1, MOE_ROW_BATCH_MAX_ROWS);
411        if m != MOE_ROW_BATCH_MAX_ROWS {
412            tracing::warn!(
413                "GLM MoE row-batch width capped at {m} (default {MOE_ROW_BATCH_MAX_ROWS})"
414            );
415        }
416        m
417    })
418}
419
420/// Widest compiled `w4a16_gemv_sw_moe_batchm_mR` tier. Mirror of the
421/// `ATLAS_MOE_BATCHM_ENTRY` list in `kernels/gb10/common/w4a16_gemv.cu` and of the
422/// `[KernelHandle; 7]` in `Glm5NextMlpKernels`.
423///
424/// πŸ”΄ Since 2026-09-02 this is a **sub-group width, not a caller contract**. The prefill sub-chunk
425/// is 16 rows wide (`glm5next_layer::PREFILL_ROWS`) because the DENSE tier widened to 16; the
426/// routed experts did not follow, so [`forward_moe`] splits any wider row group into even
427/// sub-groups of at most this and sweeps each one. Callers may pass any `rows` their workspace
428/// holds.
429pub const MOE_ROW_BATCH_MAX_ROWS: usize = 8;
430
431/// Split `rows` into consecutive `(start, width)` sub-groups of at most `cap`, as evenly as the
432/// count allows.
433///
434/// πŸͺ€ Even, not greedy. A greedy split of 9 rows at cap 8 leaves a trailing group of ONE, and
435/// there is no `w4a16_gemv_sw_moe_batchm_m1` tier β€” the array starts at m2. Balancing gives 5 + 4,
436/// and at the shipping cap of `MOE_ROW_BATCH_MAX_ROWS` every group is >= 2 for every `rows >= 2`.
437///
438/// πŸͺ€ A width-1 group is still REACHABLE at a small cap, where it is arithmetically forced (3 rows
439/// at cap 2 has no all->=2 split). That is not a correctness hole β€” the caller's gate requires
440/// every group to have a tier, so such a call simply runs the per-row arm β€” but it does mean
441/// `ATLAS_GLM_MOE_ROW_BATCH_MAX=2` silently disables row batching at odd widths. Only the A/B
442/// lever can reach it.
443///
444/// πŸ”΄ Splitting is EXACT. A row's routed output is the sum over ITS OWN top-k slots, each slot a
445/// single expert's GEMV whose accumulation order (`w4a16_gemv_partial_rows`) does not depend on
446/// `R` or on which rows share the sweep; the union table only decides which experts get swept and
447/// in what order the sweeps are issued, never what any row adds. So `forward_moe(rows)` returns
448/// the same bits however it is grouped. What splitting costs is amortization, not accuracy: two
449/// 8-row sweeps visit the union of 8 rows twice instead of the (smaller) union of 16 once.
450fn moe_row_groups(rows: usize, cap: usize) -> Vec<(usize, usize)> {
451    let n = rows.div_ceil(cap.max(1)).max(1);
452    let mut out = Vec::with_capacity(n);
453    let mut start = 0usize;
454    for i in 0..n {
455        let w = (rows - start).div_ceil(n - i);
456        out.push((start, w));
457        start += w;
458    }
459    out
460}
461
462/// πŸͺ€ `glm5next_moe_row_union` is ONE block of `rows * top_k` threads. A CUDA block is capped
463/// at 1024 threads, but this kernel's own scans are `O(T^2)`/`O(T^3)` over that extent and the
464/// tier family was sized around 64, so 64 is the contract. Threads past a block never run:
465/// exceeding it would SILENTLY drop union entries, so the dispatch refuses instead.
466pub const MOE_ROW_UNION_MAX_IDS: usize = 64;
467
468/// Say once PER ROW COUNT whether a verify shares its expert sweeps across its rows.
469///
470/// πŸͺ€ A plain `Once` here is a trap: the first MoE forward of a run is the prefill/K=1 step at
471/// `rows = 1`, which can never batch, so a single announcement reports "per-row" for the whole
472/// process and the batched path looks like it never engaged. Latch one bit per row count.
473fn announce_row_batch(batched: bool, rows: usize) {
474    use std::sync::atomic::{AtomicU16, Ordering};
475    static SEEN: AtomicU16 = AtomicU16::new(0);
476    let bit = 1u16 << rows.min(15);
477    if SEEN.fetch_or(bit, Ordering::Relaxed) & bit != 0 {
478        return;
479    }
480    if batched {
481        tracing::info!("GLM MoE: row-batched expert union ({rows} rows, one sweep each)");
482    } else {
483        tracing::info!("GLM MoE: per-row expert sweeps ({rows} rows)");
484    }
485}
486
487/// Kill switch for the grouped path: `ATLAS_GLM_MOE_HOST_DISPATCH=1` restores the
488/// read-ids-to-host expert loop. Read once β€” this sits on the per-layer decode path.
489fn host_dispatch_forced() -> bool {
490    static F: std::sync::OnceLock<bool> = std::sync::OnceLock::new();
491    *F.get_or_init(|| std::env::var("ATLAS_GLM_MOE_HOST_DISPATCH").as_deref() == Ok("1"))
492}
493
494/// Say once which expert-dispatch path this process took. A missing `w4a16_gemv_sw_moe`
495/// entry point falls back SILENTLY otherwise, and the fallback is the slow one.
496fn announce_dispatch(grouped: bool) {
497    static ONCE: std::sync::Once = std::sync::Once::new();
498    ONCE.call_once(|| {
499        if grouped {
500            tracing::info!("GLM MoE: grouped device dispatch (no per-layer D2H)");
501        } else {
502            tracing::info!("GLM MoE: host dispatch (per-layer stream sync + D2H)");
503        }
504    });
505}
506
507/// One routed MoE site, one token. Leaves a **partial sum** in `out` whenever this rank shares
508/// the experts (EP) or the shared expert (TP) with anyone else.
509///
510/// # πŸͺ€ The deviceβ†’host round trip
511///
512/// The expert loop reads the selected ids back to the host to decide which experts are local.
513/// That is a synchronising `copy_d2h` on the decode critical path, once per routed layer.
514/// It is deliberate for this slice: correctness first, and it is exactly what the gated
515/// microtest does. The upgrade path is the pointer-table grouped GEMM
516/// (`layers::moe::ptr_table_build`), which keeps the routing on device β€” not a change to
517/// this math.
518/// Routed MoE over `rows` rows.
519///
520/// πŸ”΄ The routed experts amortize PARTIALLY over a verify's rows. Measured on the live routing
521/// trace (42 series x 406 steps), the union of selected experts over K consecutive tokens is
522/// 8.00 / 13.74 / 18.76 / 23.35 at K = 1..4, so their weight traffic grows with K however the
523/// loop is written β€” but it grows along that curve, not along 8K.
524///
525/// πŸ”΄ RETRACTED (2026-08-29): this comment used to say the routed experts "stay one row at a
526/// time" and that deduplicating the union "needs a device-side sort the dispatch kernel does
527/// not have yet". Both are wrong. `top_k * rows <= 64` ids resolve in ONE block by pairwise
528/// scan β€” no sort β€” and `w4a16_gemv_sw_moe_batchm_mR` then sweeps each union expert once.
529/// Measured on t69, K=3, six probes byte-identical either way: open512 20.25 -> 22.12 tok/s
530/// (+9.2%), a 125.4 -> 114.8 ms step. At K=2 (the serving default) 18.91 -> 19.75.
531///
532/// πŸ”΄ WIDENED (2026-08-31) from `rows <= 4` to `rows <= 8`. The 4 was the compiled tier family,
533/// not the union's limit β€” `8 * 8 == 64` fits its single block exactly. This is what makes the
534/// 8-row batched PREFILL sub-chunk (ANOMALIES A65) amortize its routed experts too; before it,
535/// prefill batched every stage EXCEPT the experts, which by then were most of the traffic left.
536///
537/// The SHARED expert is a different animal again: it is the same weights for every row, so it
538/// runs once over all of them regardless.
539#[allow(clippy::too_many_arguments)]
540pub fn forward_moe(
541    gpu: &dyn GpuBackend,
542    k: &Glm5NextMlpKernels,
543    cfg: &Glm5NextMlpConfig,
544    w: &Glm5NextMoeWeights,
545    x: DevicePtr,
546    out: DevicePtr,
547    rows: usize,
548    ws: &Glm5NextMlpWorkspace,
549    stream: u64,
550) -> Result<()> {
551    if rows == 0 || rows > ws.max_rows {
552        bail!(
553            "GLM MoE: {rows} rows do not fit a workspace built for {}",
554            ws.max_rows
555        );
556    }
557    if w.experts.len() != cfg.local_experts {
558        bail!(
559            "GLM MoE: {} bound experts but this rank owns {} of {}",
560            w.experts.len(),
561            cfg.local_experts,
562            cfg.num_experts
563        );
564    }
565
566    use crate::layers::glm5next_layer::profile;
567
568    // πŸ”΄ The routed experts DO amortize over a verify's rows β€” just not fully. The union of
569    // the selected experts over K consecutive tokens is 8.00 / 13.74 / 18.76 / 23.35 at
570    // K = 1..4 (live routing trace, 42 series x 406 steps), so K rows sweep that many experts
571    // instead of 8K. The per-row path pays 8K; this one pays the union.
572    // πŸͺ€ The route trace reads `ids` back per row, which the batched path never does β€” leave
573    // it on the per-row arm rather than reconstructing the trace from the union table.
574    // Sub-groups the routed sweep runs at. `rows` may exceed the widest tier β€” the prefill
575    // sub-chunk is 16 wide since the dense tier widened β€” so every gate below is PER GROUP.
576    let groups = moe_row_groups(rows, row_batch_max());
577    let batched = rows >= 2
578        && !host_dispatch_forced()
579        && !row_batch_disabled()
580        && !profile::trace_on()
581        && k.moe_row_union.0 != 0
582        && groups.iter().all(|&(_, w)| {
583            // πŸͺ€ The union table is one block of `w * top_k` threads; past 64 ids it would
584            // silently drop entries. GLM-5.3 is 8 x 8 = 64 exactly, so this is a live edge.
585            w >= 2
586                && w * cfg.top_k <= MOE_ROW_UNION_MAX_IDS
587                && k.w4a16_gemv_sw_moe_batchm[w - 2].0 != 0
588        });
589    announce_row_batch(batched, rows);
590
591    // ── router: FULL expert set, FP32 logits, replicated on every rank ──
592    let t = profile::start();
593    for r in 0..rows {
594        gemm(
595            gpu,
596            k.gemm_f32,
597            k.gemv_f32,
598            // No FP32-out batchm twin exists; the router stays on gemv/tile.
599            KernelHandle(0),
600            x.offset(r * cfg.hidden * 2),
601            w.router,
602            ws.logits.offset(r * cfg.num_experts * 4),
603            1,
604            cfg.num_experts,
605            cfg.hidden,
606            stream,
607        )?;
608    }
609    // ONE top-k for every row. `glm5next_router_topk` already takes the row on `blockIdx.x`
610    // and strides `logits`/`ids`/`wts` by it, so this is the identical per-row work in one
611    // launch instead of K β€” and it was a `grid [1,1,1]` launch, 120 of them per K=3 step for
612    // 1.75 ms (nsys 2026-08-29). Bit-identical: no row's arithmetic changes.
613    KernelLaunch::new(gpu, k.router)
614        .grid([rows as u32, 1, 1])
615        .block([ACT_BLOCK, 1, 1])
616        .arg_ptr(ws.logits)
617        .arg_ptr(w.router_bias)
618        .arg_ptr(ws.ids)
619        .arg_ptr(ws.wts)
620        .arg_u32(cfg.num_experts as u32)
621        .arg_u32(cfg.top_k as u32)
622        // n_group: the parser already refuses anything but 1; the kernel refuses too.
623        .arg_u32(1)
624        .arg_f32(cfg.routed_scale)
625        .arg_u32(u32::from(cfg.renormalize))
626        .arg_u32(u32::from(cfg.router_bf16_ladder))
627        .launch(stream)?;
628    profile::end(profile::MOE_ROUTER, t, gpu, stream);
629
630    // πŸͺ€ Zero FIRST. A slot this rank does not own must contribute exactly zero to the
631    // all-reduced sum; leaving the previous token's expert output there is a wrong answer
632    // that only appears at EP > 1 and only for tokens whose routing moved. `expert_out` is
633    // `[rows, top_k, hidden]` and contiguous, so one memset covers every row.
634    gpu.memset_async(ws.expert_out, 0, rows * cfg.top_k * cfg.hidden * 2, stream)?;
635
636    for r in 0..rows {
637        if batched {
638            break; // the experts run once for ALL rows, after this loop
639        }
640        let xr = x.offset(r * cfg.hidden * 2);
641        let ids_r = ws.ids.offset(r * cfg.top_k * 4);
642        let expert_out_r = ws.expert_out.offset(r * cfg.top_k * cfg.hidden * 2);
643
644        let grouped = !host_dispatch_forced() && k.w4a16_gemv_sw_moe.0 != 0;
645        announce_dispatch(grouped);
646        if grouped {
647            // ── grouped, device-dispatched: routing never leaves the GPU ──
648            //
649            // πŸ”΄ The host loop this replaces did `synchronize` + `copy_d2h(ids)` once per routed
650            // layer β€” 42 full stream drains per decode token on GLM-5.3, and the reason the
651            // decode step could not be graph-captured. It also issued one launch per LOCAL
652            // expert per projection (~16/layer); this is four, whatever the routing picks.
653            let t = profile::start();
654            let mi = cfg.moe_intermediate;
655            // gate and up: every slot reads the SAME x, so input_stride = 0.
656            w4a16_gemv_moe(
657                gpu,
658                k.w4a16_gemv_sw_moe,
659                xr,
660                &w.ptrs.gate,
661                ws.a_gate,
662                ids_r,
663                mi,
664                cfg.hidden,
665                cfg.top_k,
666                cfg.num_experts,
667                0,
668                stream,
669            )?;
670            w4a16_gemv_moe(
671                gpu,
672                k.w4a16_gemv_sw_moe,
673                xr,
674                &w.ptrs.up,
675                ws.a_up,
676                ids_r,
677                mi,
678                cfg.hidden,
679                cfg.top_k,
680                cfg.num_experts,
681                0,
682                stream,
683            )?;
684            // Elementwise over all slots at once. Remote slots activate uninitialised rows; the
685            // down projection skips them, so those rows are never read.
686            swiglu(
687                gpu,
688                k.swiglu,
689                ws.a_gate,
690                ws.a_up,
691                ws.a_act,
692                cfg.top_k * mi,
693                cfg.swiglu_limit,
694                stream,
695            )?;
696            // down: slot-major activations, so input_stride = one expert's width.
697            w4a16_gemv_moe(
698                gpu,
699                k.w4a16_gemv_sw_moe,
700                ws.a_act,
701                &w.ptrs.down,
702                expert_out_r,
703                ids_r,
704                cfg.hidden,
705                mi,
706                cfg.top_k,
707                cfg.num_experts,
708                mi,
709                stream,
710            )?;
711            profile::end(profile::MOE_EXPERTS, t, gpu, stream);
712
713            if profile::trace_on() {
714                let mut ids = vec![0u8; cfg.top_k * 4];
715                gpu.synchronize(stream)?;
716                gpu.copy_d2h(ids_r, &mut ids)?;
717                let decoded: Vec<i32> = (0..cfg.top_k)
718                    .map(|k| {
719                        i32::from_le_bytes([
720                            ids[k * 4],
721                            ids[k * 4 + 1],
722                            ids[k * 4 + 2],
723                            ids[k * 4 + 3],
724                        ])
725                    })
726                    .collect();
727                profile::stash_route(&decoded);
728            }
729        } else {
730            // 🚩 A FULL STREAM SYNC + D2H IN THE MIDDLE OF EVERY MoE LAYER. The routing decision
731            // is read back to the host so the expert GEMMs can be launched by id. Timed on its own
732            // because it is the one span here that is pure latency and scales with layer count,
733            // not with weight bytes.
734            let t = profile::start();
735            let mut ids = vec![0u8; cfg.top_k * 4];
736            gpu.synchronize(stream)?;
737            gpu.copy_d2h(ids_r, &mut ids)?;
738            profile::end(profile::MOE_HOSTSYNC, t, gpu, stream);
739
740            if profile::trace_on() {
741                let decoded: Vec<i32> = (0..cfg.top_k)
742                    .map(|k| {
743                        i32::from_le_bytes([
744                            ids[k * 4],
745                            ids[k * 4 + 1],
746                            ids[k * 4 + 2],
747                            ids[k * 4 + 3],
748                        ])
749                    })
750                    .collect();
751                profile::stash_route(&decoded);
752            }
753
754            let t = profile::start();
755            for slot in 0..cfg.top_k {
756                let id = i32::from_le_bytes([
757                    ids[slot * 4],
758                    ids[slot * 4 + 1],
759                    ids[slot * 4 + 2],
760                    ids[slot * 4 + 3],
761                ]);
762                // -1 is the kernel's "slot unfilled" sentinel; it is reachable only if top_k exceeded
763                // the expert count, which `validate` refuses. Treat it as zero rather than as an index.
764                if id < 0 {
765                    continue;
766                }
767                let id = id as usize;
768                if id >= cfg.num_experts {
769                    bail!(
770                        "GLM MoE: router selected expert {id} of {}",
771                        cfg.num_experts
772                    );
773                }
774                let Some(local) = cfg.local_slot(id) else {
775                    continue; // another rank owns it; its zero row is already in place.
776                };
777                let e = &w.experts[local];
778                let dst = expert_out_r.offset(slot * cfg.hidden * 2);
779                let mi = cfg.moe_intermediate;
780                w4a16_gemv(
781                    gpu,
782                    k.w4a16_gemv,
783                    k.w4a16_gemv_sw,
784                    xr,
785                    &e.gate_proj,
786                    ws.a_gate,
787                    mi,
788                    cfg.hidden,
789                    stream,
790                )?;
791                w4a16_gemv(
792                    gpu,
793                    k.w4a16_gemv,
794                    k.w4a16_gemv_sw,
795                    xr,
796                    &e.up_proj,
797                    ws.a_up,
798                    mi,
799                    cfg.hidden,
800                    stream,
801                )?;
802                swiglu(
803                    gpu,
804                    k.swiglu,
805                    ws.a_gate,
806                    ws.a_up,
807                    ws.a_act,
808                    mi,
809                    cfg.swiglu_limit,
810                    stream,
811                )?;
812                w4a16_gemv(
813                    gpu,
814                    k.w4a16_gemv,
815                    k.w4a16_gemv_sw,
816                    ws.a_act,
817                    &e.down_proj,
818                    dst,
819                    cfg.hidden,
820                    mi,
821                    stream,
822                )?;
823            }
824
825            profile::end(profile::MOE_EXPERTS, t, gpu, stream);
826        }
827    }
828
829    if batched {
830        let t = profile::start();
831        let mi = cfg.moe_intermediate;
832        // ONE sweep per sub-group. At `rows <= MOE_ROW_BATCH_MAX_ROWS` this is the single pass it
833        // always was; a wider prefill sub-chunk runs it twice over disjoint row ranges, which is
834        // byte-identical (see `moe_row_groups`) and keeps the routed experts at the tier width
835        // that was actually measured.
836        for &(r0, w_rows) in &groups {
837            // The union table: one block, one thread per (row, slot) id. Stays on device.
838            // πŸͺ€ Rebuilt per sub-group over that group's slice of `ids` β€” the scratch is sized for
839            // the widest group, and a later group overwrites the earlier one's table after its
840            // sweeps have been issued on the same stream.
841            KernelLaunch::new(gpu, k.moe_row_union)
842                .grid([1, 1, 1])
843                .block([(w_rows * cfg.top_k) as u32, 1, 1])
844                .arg_ptr(ws.ids.offset(r0 * cfg.top_k * 4))
845                .arg_ptr(ws.u_eid)
846                .arg_ptr(ws.u_slot)
847                .arg_u32(w_rows as u32)
848                .arg_u32(cfg.top_k as u32)
849                .launch(stream)?;
850
851            let kb = k.w4a16_gemv_sw_moe_batchm[w_rows - 2];
852            // gate and up: a row's slots all read the SAME x, so the slot stride is 0.
853            w4a16_gemv_moe_batchm(
854                gpu,
855                kb,
856                x.offset(r0 * cfg.hidden * 2),
857                &w.ptrs.gate,
858                ws.a_gate.offset(r0 * cfg.top_k * mi * 2),
859                ws.u_eid,
860                ws.u_slot,
861                mi,
862                cfg.hidden,
863                w_rows,
864                cfg.top_k,
865                cfg.num_experts,
866                cfg.hidden,
867                0,
868                cfg.top_k * mi,
869                stream,
870            )?;
871            w4a16_gemv_moe_batchm(
872                gpu,
873                kb,
874                x.offset(r0 * cfg.hidden * 2),
875                &w.ptrs.up,
876                ws.a_up.offset(r0 * cfg.top_k * mi * 2),
877                ws.u_eid,
878                ws.u_slot,
879                mi,
880                cfg.hidden,
881                w_rows,
882                cfg.top_k,
883                cfg.num_experts,
884                cfg.hidden,
885                0,
886                cfg.top_k * mi,
887                stream,
888            )?;
889            // Elementwise over every (row, slot) at once. Slots this rank does not own activate
890            // uninitialised rows; the down projection skips them, so those rows are never read.
891            swiglu(
892                gpu,
893                k.swiglu,
894                ws.a_gate.offset(r0 * cfg.top_k * mi * 2),
895                ws.a_up.offset(r0 * cfg.top_k * mi * 2),
896                ws.a_act.offset(r0 * cfg.top_k * mi * 2),
897                w_rows * cfg.top_k * mi,
898                cfg.swiglu_limit,
899                stream,
900            )?;
901            // down: slot-major activations, so the slot stride is one expert's width.
902            w4a16_gemv_moe_batchm(
903                gpu,
904                kb,
905                ws.a_act.offset(r0 * cfg.top_k * mi * 2),
906                &w.ptrs.down,
907                ws.expert_out.offset(r0 * cfg.top_k * cfg.hidden * 2),
908                ws.u_eid,
909                ws.u_slot,
910                cfg.hidden,
911                mi,
912                w_rows,
913                cfg.top_k,
914                cfg.num_experts,
915                cfg.top_k * mi,
916                mi,
917                cfg.top_k * cfg.hidden,
918                stream,
919            )?;
920        }
921        profile::end(profile::MOE_EXPERTS, t, gpu, stream);
922    }
923
924    // ── shared expert: BF16, TP-sharded, NOT routed-scaled ──
925    let t = profile::start();
926    forward_dense(
927        gpu,
928        k,
929        cfg,
930        &w.shared,
931        cfg.local_shared_intermediate,
932        x,
933        ws.shared_out,
934        rows,
935        ws,
936        stream,
937    )?;
938
939    // πŸ”΄ The combine runs BEFORE the all-reduce, so the TP-partial shared expert and the
940    // EP-partial routed sum reduce together in one collective. Adding the shared output after
941    // a reduce β€” the `layers::moe` pattern, written for a replicated shared expert β€” would
942    // keep only this rank's half of it.
943    profile::end(profile::MOE_SHARED, t, gpu, stream);
944    let t = profile::start();
945    // ONE combine for every row: `glm5next_moe_combine` takes the row on `blockIdx.x` and
946    // strides all four buffers by it. Was K `grid [1,1,1]` launches β€” 1.50 ms of a K=3 step.
947    KernelLaunch::new(gpu, k.combine)
948        .grid([rows as u32, 1, 1])
949        .block([ACT_BLOCK, 1, 1])
950        .arg_ptr(ws.expert_out)
951        .arg_ptr(ws.wts)
952        .arg_ptr(ws.shared_out)
953        .arg_ptr(out)
954        .arg_u32(cfg.hidden as u32)
955        .arg_u32(cfg.top_k as u32)
956        .launch(stream)?;
957    profile::end(profile::MOE_COMBINE, t, gpu, stream);
958    Ok(())
959}
960
961#[cfg(test)]
962mod tests {
963    use super::{MOE_ROW_BATCH_MAX_ROWS, moe_row_groups};
964
965    /// The split must cover every row exactly once, in order, never exceed the tier width,
966    /// and never emit a group of ONE β€” there is no `w4a16_gemv_sw_moe_batchm_m1`, so a
967    /// trailing single row would silently drop the whole batched arm for that sub-chunk.
968    #[test]
969    fn row_groups_cover_and_never_orphan_a_row() {
970        for cap in 1..=MOE_ROW_BATCH_MAX_ROWS {
971            for rows in 1..=64 {
972                let g = moe_row_groups(rows, cap);
973                assert_eq!(g[0].0, 0, "rows={rows} cap={cap}: does not start at 0");
974                let mut next = 0;
975                for &(start, w) in &g {
976                    assert_eq!(start, next, "rows={rows} cap={cap}: gap or overlap");
977                    assert!(w >= 1 && w <= cap, "rows={rows} cap={cap}: width {w}");
978                    // No orphan at the width that ships. At a small cap an all->=2 split can be
979                    // arithmetically impossible (3 rows at cap 2), and the caller's per-group
980                    // tier gate handles that by falling back β€” see the fn doc.
981                    if rows >= 2 && cap == MOE_ROW_BATCH_MAX_ROWS {
982                        assert!(w >= 2, "rows={rows} cap={cap}: orphaned a single row");
983                    }
984                    next += w;
985                }
986                assert_eq!(next, rows, "rows={rows} cap={cap}: {next} rows covered");
987            }
988        }
989    }
990
991    /// The two widths this actually ships at.
992    #[test]
993    fn row_groups_at_the_shipping_widths() {
994        assert_eq!(moe_row_groups(16, 8), vec![(0, 8), (8, 8)]);
995        assert_eq!(moe_row_groups(8, 8), vec![(0, 8)]);
996        assert_eq!(moe_row_groups(9, 8), vec![(0, 5), (5, 4)]);
997    }
998}