diff --git a/tests/pytorch/test_grouped_gemm_mxfp4.py b/tests/pytorch/test_grouped_gemm_mxfp4.py
new file mode 100644
index 0000000000..18984c3c19
--- /dev/null
+++ b/tests/pytorch/test_grouped_gemm_mxfp4.py
@@ -0,0 +1,188 @@
+# Copyright (c) 2026, Advanced Micro Devices, Inc. All rights reserved.
+# License for AMD contributions = MIT. See LICENSE for more information
+
+"""Numeric correctness tests for the grouped MXFP4 Triton GEMM (gfx950).
+
+These validate the two layout assumptions of the port: that TE's plain
+E8M0/E2M1 quantizer output is interpreted the same way by
+``tl.dot_scaled(..., "e2m1", ...)``, and that zero-padded rows contribute
+nothing to the variable-K (wgrad) contraction.
+
+Each op is checked against a *precise* reference: the same packed MXFP4 operands
+the kernel consumes are dequantized here by an independent OCP E2M1/E8M0 decoder
+(``_dequant_mxfp4``) and matmul'd in fp32. Kernel and reference start from
+identical fp4 values, so a correct kernel matches to ~bf16-rounding while a
+layout / nibble-order / scale-bias / transpose bug makes them disagree grossly.
+wgrad additionally checks against the true bf16 grouped matmul, which exercises
+the per-group zero-padding (padded rows must contribute nothing).
+"""
+
+import pytest
+import torch
+
+triton = pytest.importorskip("triton")
+
+try:
+    from transformer_engine.pytorch.quantization import check_mxfp4_support
+
+    _MXFP4_OK, _MXFP4_REASON = check_mxfp4_support()
+except Exception as exc:  # pragma: no cover - import/support probe
+    _MXFP4_OK, _MXFP4_REASON = False, str(exc)
+
+from transformer_engine.pytorch.triton_kernels.grouped_gemm_mxfp4_impl import (
+    MXFP4_BLOCK,
+    _col_operand,
+    _col_operand_grouped_padded,
+    _row_operand,
+    grouped_gemm_mxfp4_dgrad,
+    grouped_gemm_mxfp4_fprop,
+    grouped_gemm_mxfp4_wgrad,
+    grouped_linear_mxfp4,
+)
+
+pytestmark = [
+    pytest.mark.skipif(not torch.cuda.is_available(), reason="CUDA/ROCm device required"),
+    pytest.mark.skipif(not _MXFP4_OK, reason=f"MXFP4 unsupported: {_MXFP4_REASON}"),
+]
+
+DTYPE = torch.bfloat16
+# Uneven, non-128-multiple group sizes exercise fprop/dgrad masking and the
+# wgrad per-group zero-padding to 128.
+M_SPLITS = [96, 128, 160, 128]
+N, K = 256, 128
+# Precise-reference bar: kernel vs dequant-of-the-same-operands differ only by
+# ~bf16 output rounding.
+_TIGHT_TOL = 3.0e-2
+# Loose bar vs the true (unquantized) bf16 matmul: ~0.16 MXFP4 noise floor for
+# these shapes, well below the ~1.4 an uncorrelated (layout-bug) output gives.
+_REL_TOL = 2.0e-1
+
+# OCP E2M1 magnitude indexed by the 3 low bits (exp2 | mantissa1); bit 3 = sign.
+_E2M1_MAG = [0.0, 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 6.0]
+
+
+@pytest.fixture(autouse=True)
+def _mxfp4_env(monkeypatch):
+    # Triton MXFP4 quantizer (no aiter dependency); seed for a deterministic
+    # quantization error.
+    monkeypatch.setenv("NVTE_USE_CAST_TRANSPOSE_TRITON", "1")
+    torch.manual_seed(0)
+
+
+def _rand(*shape):
+    return torch.randn(*shape, dtype=DTYPE, device="cuda")
+
+
+def _rel_err(out, ref):
+    out = out.float()
+    ref = ref.float()
+    return (out - ref).norm() / ref.norm().clamp_min(1e-12)
+
+
+def _dequant_mxfp4(data_u8, scale_u8, feat):
+    """Independent OCP MXFP4 dequant of packed operands -> fp32 ``[R, feat]``.
+
+    ``data_u8`` ``[R, feat/2]`` packs two E2M1 codes per byte (low nibble = even
+    index along ``feat``); ``scale_u8`` ``[R, feat/32]`` is one E8M0 scale
+    (value ``2**(x-127)``) per 1x32 block.
+    """
+    lut = torch.tensor(_E2M1_MAG, dtype=torch.float32, device=data_u8.device)
+    lo = (data_u8 & 0xF).to(torch.long)
+    hi = ((data_u8 >> 4) & 0xF).to(torch.long)
+    codes = torch.stack((lo, hi), dim=-1).reshape(data_u8.shape[0], feat)
+    mag = lut[codes & 0x7]
+    vals = torch.where((codes & 0x8).bool(), -mag, mag)
+    scale = torch.exp2(scale_u8.to(torch.float32) - 127.0)
+    return vals * scale.repeat_interleave(MXFP4_BLOCK, dim=1)
+
+
+def test_fprop_precise():
+    total_m = sum(M_SPLITS)
+    a = _rand(total_m, K)
+    weights = [_rand(N, K) for _ in M_SPLITS]
+
+    out = grouped_gemm_mxfp4_fprop(a, weights, M_SPLITS, out_dtype=DTYPE)
+
+    # Dequantize the same row-wise operands the kernel used, then grouped matmul.
+    a_deq = _dequant_mxfp4(*_row_operand(a), K)
+    ref = torch.empty((total_m, N), dtype=torch.float32, device="cuda")
+    start = 0
+    for w, m in zip(weights, M_SPLITS):
+        w_deq = _dequant_mxfp4(*_row_operand(w), K)  # [N, K]
+        ref[start : start + m] = a_deq[start : start + m] @ w_deq.t()
+        start += m
+
+    assert out.shape == (total_m, N)
+    assert _rel_err(out, ref) < _TIGHT_TOL
+
+
+def test_dgrad_precise():
+    total_m = sum(M_SPLITS)
+    grad_out = _rand(total_m, N)
+    weights = [_rand(N, K) for _ in M_SPLITS]
+
+    dgrad = grouped_gemm_mxfp4_dgrad(grad_out, weights, M_SPLITS, out_dtype=DTYPE)
+
+    # gradO row-wise, weight col-wise (the transposed operand): dA = gradO @ W.
+    go_deq = _dequant_mxfp4(*_row_operand(grad_out), N)  # [total_M, N]
+    ref = torch.empty((total_m, K), dtype=torch.float32, device="cuda")
+    start = 0
+    for w, m in zip(weights, M_SPLITS):
+        w_col_deq = _dequant_mxfp4(*_col_operand(w), N)  # [K, N] ~ W^T
+        ref[start : start + m] = go_deq[start : start + m] @ w_col_deq.t()
+        start += m
+
+    assert dgrad.shape == (total_m, K)
+    assert _rel_err(dgrad, ref) < _TIGHT_TOL
+
+
+def test_wgrad_precise_and_padding():
+    total_m = sum(M_SPLITS)
+    a = _rand(total_m, K)
+    grad_out = _rand(total_m, N)
+
+    wgrad = grouped_gemm_mxfp4_wgrad(a, grad_out, M_SPLITS, out_dtype=DTYPE)
+    assert wgrad.shape == (len(M_SPLITS), N, K)
+
+    # Precise: dequant the exact per-group padded col operands the kernel reduced
+    # over, sliced by the padded offsets. C[g] = lhs[:, g] @ rhs[:, g]^T.
+    lhs_data, lhs_scale, go_pad = _col_operand_grouped_padded(grad_out, M_SPLITS)  # [N, Mpad/2]
+    rhs_data, rhs_scale, _ = _col_operand_grouped_padded(a, M_SPLITS)  # [K, Mpad/2]
+    m_pad_total = lhs_data.shape[1] * 2
+    lhs_deq = _dequant_mxfp4(lhs_data, lhs_scale, m_pad_total)  # [N, Mpad]
+    rhs_deq = _dequant_mxfp4(rhs_data, rhs_scale, m_pad_total)  # [K, Mpad]
+    go = go_pad.tolist()
+    ref = torch.empty((len(M_SPLITS), N, K), dtype=torch.float32, device="cuda")
+    for g in range(len(M_SPLITS)):
+        s, e = go[g], go[g + 1]
+        ref[g] = lhs_deq[:, s:e] @ rhs_deq[:, s:e].t()
+    assert _rel_err(wgrad, ref) < _TIGHT_TOL
+
+    # End-to-end vs the true (unquantized) matmul: exercises the per-group
+    # zero-padding -- the padded rows must contribute nothing.
+    ref_true = torch.empty((len(M_SPLITS), N, K), dtype=torch.float32, device="cuda")
+    start = 0
+    for g, m in enumerate(M_SPLITS):
+        ref_true[g] = grad_out[start : start + m].float().t() @ a[start : start + m].float()
+        start += m
+    assert _rel_err(wgrad, ref_true) < _REL_TOL
+
+
+def test_autograd_matches_ops():
+    total_m = sum(M_SPLITS)
+    a = _rand(total_m, K).requires_grad_(True)
+    weight = torch.stack([_rand(N, K) for _ in M_SPLITS], dim=0).requires_grad_(True)
+
+    out = grouped_linear_mxfp4(a, weight, M_SPLITS)
+    assert out.shape == (total_m, N)
+
+    grad_out = _rand(total_m, N)
+    out.backward(grad_out)
+
+    # The autograd grads must equal the direct op calls (same code path).
+    ref_da = grouped_gemm_mxfp4_dgrad(
+        grad_out, list(weight.detach().unbind(0)), M_SPLITS, out_dtype=a.dtype
+    )
+    ref_dw = grouped_gemm_mxfp4_wgrad(a.detach(), grad_out, M_SPLITS, out_dtype=weight.dtype)
+    torch.testing.assert_close(a.grad, ref_da)
+    torch.testing.assert_close(weight.grad, ref_dw)
diff --git a/transformer_engine/pytorch/triton_kernels/grouped_gemm_mxfp4.py b/transformer_engine/pytorch/triton_kernels/grouped_gemm_mxfp4.py
new file mode 100644
index 0000000000..f05652cb23
--- /dev/null
+++ b/transformer_engine/pytorch/triton_kernels/grouped_gemm_mxfp4.py
@@ -0,0 +1,601 @@
+# Copyright (c) 2026, Advanced Micro Devices, Inc. All rights reserved.
+# License for AMD contributions = MIT. See LICENSE for more information
+
+"""Grouped MXFP4 (MX_BLOCKWISE) GEMM Triton persistent kernels (gfx950).
+
+Adapted for TransformerEngine from Primus-Turbo
+(https://github.com/AMD-AGI/Primus-Turbo, MIT) -- the block-scaled grouped FP4
+kernels in ``primus_turbo/triton/grouped_gemm/grouped_gemm_fp4_kernel.py``. The
+kernel bodies are kept faithful to the source; the shared scaffolding
+(``NUM_XCDS``, ``_chiplet_transform_chunked``, the AMD Triton-knob scope, the
+output padding-tail pass) is inlined here so the module has no Primus-Turbo
+dependency.
+
+The operands are E2M1 FP4 packed two-values-per-byte along the contraction (K)
+axis, fed to the hardware block-scaled MMA via ``tl.dot_scaled(..., "e2m1", ...)``
+with E8M0 (VEC_SIZE=32) scales.
+
+  - _grouped_mxfp4_persistent_gemm_kernel:      Forward (NT)  C[g] = A[g] @ B[g]^T
+  - grouped_gemm_mxfp4_triton_kernel:           Forward public API
+  - _grouped_mxfp4_variable_k_gemm_kernel:      Backward wgrad C[g] = LHS[g] @ RHS[g]^T
+  - grouped_gemm_mxfp4_variable_k_triton_kernel: Backward (variable-K) public API
+
+FP4 packing notes:
+  * data tensors store K/2 (resp. M/2) bytes along the contraction axis;
+  * scale tensors store K/32 (resp. M/32) E8M0 bytes (one per 1x32 block);
+  * scales are stored free-major (free, K/32) / (G, free, K/32); each K-iter
+    loads a (K/32, free) tile (coalesced) and transposes it back in-reg for
+    tl.dot_scaled.
+
+Only EVEN contraction (a multiple of BLOCK_SIZE_K = 128 logical elements) is
+supported; the FP4 quantizer pads K to 128 and the MX wrapper pads per-group M
+to 128, so this always holds in practice.
+"""
+
+from __future__ import annotations
+
+import contextlib
+import functools
+import os
+from typing import Optional
+
+import torch
+import triton
+import triton.language as tl
+
+
+# ===============================================================================
+# Arch helper
+# ===============================================================================
+
+
+@functools.lru_cache
+def _is_gfx950() -> bool:
+    props = torch.cuda.get_device_properties(torch.cuda.current_device())
+    return (props.major, props.minor) == (9, 5)
+
+
+# ===============================================================================
+# AMD Triton compiler knobs
+#
+# The grouped FP4 kernels want the gfx950 knob set (async_copy,
+# scalarize_packed_fops, block_pingpong) active only while they compile/launch,
+# then restored so they never leak into other Triton kernels. These knobs are
+# read at compile time and are not part of Triton's compile cache key.
+# ===============================================================================
+
+_AMD_KNOB_ATTRS = ("use_async_copy", "scalarize_packed_fops", "use_block_pingpong")
+_AMD_KNOB_ENVS = (
+    "TRITON_HIP_USE_ASYNC_COPY",
+    "AMDGCN_SCALARIZE_PACKED_FOPS",
+    "TRITON_HIP_USE_BLOCK_PINGPONG",
+)
+
+
+def _amd_knobs_available() -> bool:
+    return hasattr(triton, "knobs") and hasattr(triton.knobs, "amd")
+
+
+@contextlib.contextmanager
+def _scoped_amd_triton_knobs():
+    """Snapshot AMD Triton compiler knobs on entry, restore them on exit.
+
+    On gfx950 the full gfx950 knob set is enabled for the duration of the scope;
+    all changes are restored on exit.
+    """
+    saved_attrs = {}
+    if _amd_knobs_available():
+        amd = triton.knobs.amd
+        for name in _AMD_KNOB_ATTRS:
+            if hasattr(amd, name):
+                saved_attrs[name] = getattr(amd, name)
+    saved_env = {name: os.environ.get(name) for name in _AMD_KNOB_ENVS}
+    try:
+        if _is_gfx950():
+            if _amd_knobs_available():
+                amd = triton.knobs.amd
+                for name in _AMD_KNOB_ATTRS:
+                    if hasattr(amd, name):
+                        setattr(amd, name, True)
+            else:
+                for name in _AMD_KNOB_ENVS:
+                    os.environ[name] = "1"
+        yield
+    finally:
+        if saved_attrs:
+            amd = triton.knobs.amd
+            for name, val in saved_attrs.items():
+                setattr(amd, name, val)
+        for name, val in saved_env.items():
+            if val is None:
+                os.environ.pop(name, None)
+            else:
+                os.environ[name] = val
+
+
+def _scoped_amd_knobs(func):
+    """Decorator: run ``func`` under :func:`_scoped_amd_triton_knobs`."""
+
+    @functools.wraps(func)
+    def wrapper(*args, **kwargs):
+        with _scoped_amd_triton_knobs():
+            return func(*args, **kwargs)
+
+    return wrapper
+
+
+# ===============================================================================
+# Chiplet transform: map the AMD round-robin program_id onto XCD-contiguous
+# chunks so a band of programs shares B columns in L2.
+# ===============================================================================
+
+NUM_XCDS = 8  # gfx942/gfx950 chiplet count
+
+
+@triton.jit
+def _chiplet_transform_chunked(
+    pid,
+    NUM_SMS: tl.constexpr,
+    NUM_XCDS: tl.constexpr,
+    CHUNK_SIZE: tl.constexpr,
+):
+    if pid > (NUM_SMS // (NUM_XCDS * CHUNK_SIZE)) * (NUM_XCDS * CHUNK_SIZE):
+        return pid
+    local_pid = pid // NUM_XCDS
+    chunk_idx = local_pid // CHUNK_SIZE
+    pos_in_chunk = local_pid % CHUNK_SIZE
+    xcd = pid % NUM_XCDS
+    return chunk_idx * NUM_XCDS * CHUNK_SIZE + xcd * CHUNK_SIZE + pos_in_chunk
+
+
+# E8M0 block size: one scale per 1x32 logical-element block (same as MXFP8).
+# tl.dot_scaled consumes the "e2m1" FP4 format string (inlined in the kernels;
+# @triton.jit cannot read module globals).
+VEC_SIZE = 32
+
+
+# ###########################################################################
+#  Forward (NT):   C[g] = A[g] @ B[g]^T
+#    A (total_M, K/2) fp4-packed, A_scale (total_M, K/32) e8m0
+#    B (G, N, K/2)    fp4-packed, B_scale (G, N, K/32)   e8m0
+# ###########################################################################
+
+
+@triton.jit
+def _grouped_mxfp4_persistent_gemm_kernel(
+    A,
+    B,
+    C,
+    A_scale,  # (total_M, K//32) uint8 e8m0
+    B_scale,  # (G, N, K//32) uint8 e8m0
+    rd_offs_ptr,  # (G+1) int64 — read base along M (A / A_scale)
+    out_offs_ptr,  # (G+1) int64 — write base along M (C, tile count)
+    G,
+    N,
+    K,  # logical contraction size (multiple of BLOCK_SIZE_K)
+    stride_am,
+    stride_ak,  # along packed K bytes
+    stride_bg,
+    stride_bn,
+    stride_bk,  # along packed K bytes
+    stride_cm,
+    stride_cn,
+    stride_asm,
+    stride_ask,
+    stride_bsg,
+    stride_bsn,
+    stride_bsk,
+    BLOCK_SIZE_M: tl.constexpr,
+    BLOCK_SIZE_N: tl.constexpr,
+    BLOCK_SIZE_K: tl.constexpr,  # logical K per iter
+    GROUP_SIZE_M: tl.constexpr,
+    NUM_SMS: tl.constexpr,
+    NUM_XCDS: tl.constexpr,
+    CHUNK_SIZE: tl.constexpr,
+    CACHE_MODIFIER: tl.constexpr,
+    VEC: tl.constexpr,
+):
+    pid = tl.program_id(0)
+    if NUM_XCDS != 1:
+        pid = _chiplet_transform_chunked(pid, NUM_SMS, NUM_XCDS, CHUNK_SIZE)
+    num_pid_n = tl.cdiv(N, BLOCK_SIZE_N)
+
+    total_tiles: tl.int32 = 0
+    for _g in range(G):
+        m_g = (tl.load(out_offs_ptr + _g + 1) - tl.load(out_offs_ptr + _g)).to(tl.int32)
+        total_tiles += tl.cdiv(m_g, BLOCK_SIZE_M) * num_pid_n
+
+    tl.assume(stride_am > 0)
+    tl.assume(stride_cm > 0)
+
+    BK_PACK: tl.constexpr = BLOCK_SIZE_K // 2  # packed bytes per K-iter
+    BK_SCALE: tl.constexpr = BLOCK_SIZE_K // VEC  # scale entries per K-iter
+
+    for global_tile_id in range(pid, total_tiles, NUM_SMS):
+        # ── locate group (linear scan) ──
+        group_idx: tl.int32 = 0
+        tile_start: tl.int32 = 0
+        cumsum: tl.int32 = 0
+        for _g in range(G):
+            m_g_i = (tl.load(out_offs_ptr + _g + 1) - tl.load(out_offs_ptr + _g)).to(tl.int32)
+            tiles_g = tl.cdiv(m_g_i, BLOCK_SIZE_M) * num_pid_n
+            new_cumsum = cumsum + tiles_g
+            if global_tile_id >= new_cumsum:
+                group_idx = _g + 1
+                tile_start = new_cumsum
+            cumsum = new_cumsum
+
+        local_tile = global_tile_id - tile_start
+        m_rd = tl.load(rd_offs_ptr + group_idx)  # int64 read base (A / A_scale)
+        m_out = tl.load(out_offs_ptr + group_idx)  # int64 write base (C)
+        M_g = (tl.load(out_offs_ptr + group_idx + 1) - m_out).to(tl.int32)
+        tiles_m_g = tl.cdiv(M_g, BLOCK_SIZE_M)
+
+        # ── swizzle ──
+        num_pid_in_group = GROUP_SIZE_M * num_pid_n
+        swizzle_group = local_tile // num_pid_in_group
+        first_pid_m = swizzle_group * GROUP_SIZE_M
+        group_size_m = min(tiles_m_g - first_pid_m, GROUP_SIZE_M)
+        pid_m = first_pid_m + ((local_tile % num_pid_in_group) % group_size_m)
+        pid_n = (local_tile % num_pid_in_group) // group_size_m
+        tl.assume(pid_m >= 0)
+        tl.assume(pid_n >= 0)
+
+        rm = (pid_m * BLOCK_SIZE_M + tl.arange(0, BLOCK_SIZE_M)) % M_g
+        rn = (pid_n * BLOCK_SIZE_N + tl.arange(0, BLOCK_SIZE_N)) % N
+        rk = tl.arange(0, BK_PACK)  # packed K bytes
+        rks = tl.arange(0, BK_SCALE)  # scale entries
+        rn = tl.max_contiguous(tl.multiple_of(rn, BLOCK_SIZE_N), BLOCK_SIZE_N)
+
+        A_BASE = A + (m_rd + rm[:, None]) * stride_am + rk[None, :] * stride_ak
+        B_BASE = B + group_idx.to(tl.int64) * stride_bg + rk[:, None] * stride_bk + rn[None, :] * stride_bn
+        AS_BASE = A_scale + rks[:, None] * stride_ask + (m_rd + rm[None, :]) * stride_asm
+        BS_BASE = (
+            B_scale
+            + group_idx.to(tl.int64) * stride_bsg
+            + rks[:, None] * stride_bsk
+            + rn[None, :] * stride_bsn
+        )
+
+        acc = tl.zeros((BLOCK_SIZE_M, BLOCK_SIZE_N), dtype=tl.float32)
+        loop_k = K // BLOCK_SIZE_K
+        for ki in range(0, loop_k):
+            a = tl.load(tl.multiple_of(A_BASE, (1, 16)), cache_modifier=CACHE_MODIFIER)  # (BM, BK/2)
+            b = tl.load(tl.multiple_of(B_BASE, (16, 1)), cache_modifier=CACHE_MODIFIER)  # (BK/2, BN)
+            a_s = tl.trans(tl.load(AS_BASE))  # (BK/32, BM) -> (BM, BK/32)
+            b_s = tl.trans(tl.load(BS_BASE))  # (BK/32, BN) -> (BN, BK/32)
+            acc = tl.dot_scaled(a, a_s, "e2m1", b, b_s, "e2m1", acc)
+            A_BASE += BK_PACK * stride_ak
+            B_BASE += BK_PACK * stride_bk
+            AS_BASE += BK_SCALE * stride_ask
+            BS_BASE += BK_SCALE * stride_bsk
+
+        c = acc.to(C.type.element_ty)
+        rm_s = (pid_m * BLOCK_SIZE_M + tl.arange(0, BLOCK_SIZE_M)) % M_g
+        rn_s = (pid_n * BLOCK_SIZE_N + tl.arange(0, BLOCK_SIZE_N)) % N
+        rn_s = tl.max_contiguous(tl.multiple_of(rn_s, BLOCK_SIZE_N), BLOCK_SIZE_N)
+        c_mask = (rm_s[:, None] < M_g) & (rn_s[None, :] < N)
+        C_ = C + (m_out + rm_s[:, None]) * stride_cm + rn_s[None, :] * stride_cn
+        tl.store(C_, c, c_mask)
+
+
+@_scoped_amd_knobs
+def grouped_gemm_mxfp4_triton_kernel(
+    a,
+    a_scale,
+    b,
+    b_scale,
+    group_offs,
+    N,
+    K,
+    group_offs_out=None,
+    out_dtype=torch.bfloat16,
+    num_cu=None,
+):
+    """A(total_M, K/2) @ B(G, N, K/2)^T -> C. FP4-packed data, e8m0 uint8 scales.
+
+    group_offs:     read offsets along M for A / A_scale.
+    group_offs_out: write offsets along M for C (defaults to group_offs).
+    K is the logical contraction size (must be a multiple of BLOCK_SIZE_K=128).
+    """
+    if group_offs_out is None:
+        group_offs_out = group_offs
+    G = b.shape[0]
+    c = torch.empty((a.shape[0], N), dtype=out_dtype, device=a.device)
+    a_s = a_scale.view(torch.uint8)
+    b_s = b_scale.view(torch.uint8)
+    a_u8 = a.view(torch.uint8)
+    b_u8 = b.view(torch.uint8)
+    cu = num_cu if num_cu is not None else torch.cuda.get_device_properties(a.device).multi_processor_count
+    BM, BN, BK = 256, 256, 128
+    m_alloc = a.shape[0]
+    avg_m = max(m_alloc // max(G, 1), 1)
+    tiles_n = (N + BN - 1) // BN
+    GM = 8 if min((avg_m + BM - 1) // BM, tiles_n) < 16 else 4
+    total_tiles = ((m_alloc + BM - 1) // BM + G) * tiles_n
+    num_sms = min(total_tiles, cu)
+    chunk = 64 if num_sms >= NUM_XCDS * 64 else 32
+    grid = (num_sms,)
+    _grouped_mxfp4_persistent_gemm_kernel[grid](
+        a_u8,
+        b_u8,
+        c,
+        a_s,
+        b_s,
+        group_offs,
+        group_offs_out,
+        G,
+        N,
+        K,
+        a_u8.stride(0),
+        a_u8.stride(1),
+        b_u8.stride(0),
+        b_u8.stride(1),
+        b_u8.stride(2),
+        c.stride(0),
+        c.stride(1),
+        a_s.stride(0),  # stride_asm (M);   a_s is (total_M, K/32)
+        a_s.stride(1),  # stride_ask (K/32 contiguous)
+        b_s.stride(0),  # stride_bsg
+        b_s.stride(1),  # stride_bsn (N);   b_s is (G, N, K/32)
+        b_s.stride(2),  # stride_bsk (K/32 contiguous)
+        BLOCK_SIZE_M=BM,
+        BLOCK_SIZE_N=BN,
+        BLOCK_SIZE_K=BK,
+        GROUP_SIZE_M=GM,
+        NUM_SMS=num_sms,
+        NUM_XCDS=NUM_XCDS,
+        CHUNK_SIZE=chunk,
+        CACHE_MODIFIER=".ca",
+        VEC=VEC_SIZE,
+        num_warps=8,
+        num_stages=2,
+        waves_per_eu=0,
+        matrix_instr_nonkdim=16,
+        kpack=1,
+    )
+    return c
+
+
+# ###########################################################################
+#  Variable-K backward (wgrad): C[g] = LHS[g] @ RHS[g]^T, reduction over M_g
+#    LHS (OUT_M, M_total/2) fp4-packed, LHS_scale (OUT_M, M_total/32) e8m0
+#    RHS (OUT_N, M_total/2) fp4-packed, RHS_scale (OUT_N, M_total/32) e8m0
+#    C   (G, OUT_M, OUT_N)
+#    go_pad: padded per-group offsets along M (each M_g a multiple of 128)
+# ###########################################################################
+
+
+@triton.jit
+def _grouped_mxfp4_variable_k_gemm_kernel(
+    LHS,
+    RHS,
+    C,
+    LHS_scale,
+    RHS_scale,
+    go_pad_ptr,
+    G,
+    OUT_M,
+    OUT_N,
+    stride_lm,
+    stride_lk,  # along packed M bytes
+    stride_rm,
+    stride_rk,  # along packed M bytes
+    stride_cg,
+    stride_cm,
+    stride_cn,
+    stride_lsm,
+    stride_lsk,
+    stride_rsm,
+    stride_rsk,
+    BLOCK_SIZE_M: tl.constexpr,  # over OUT_M
+    BLOCK_SIZE_N: tl.constexpr,  # over OUT_N
+    BLOCK_SIZE_K: tl.constexpr,  # logical reduction over M_g
+    GROUP_SIZE_M: tl.constexpr,
+    NUM_SMS: tl.constexpr,
+    NUM_XCDS: tl.constexpr,
+    CHUNK_SIZE: tl.constexpr,
+    CACHE_MODIFIER: tl.constexpr,
+    VEC: tl.constexpr,
+):
+    pid = tl.program_id(0)
+    if NUM_XCDS != 1:
+        pid = _chiplet_transform_chunked(pid, NUM_SMS, NUM_XCDS, CHUNK_SIZE)
+    tiles_m = tl.cdiv(OUT_M, BLOCK_SIZE_M)
+    tiles_n = tl.cdiv(OUT_N, BLOCK_SIZE_N)
+    tiles_per_group = tiles_m * tiles_n
+    total_tiles = G * tiles_per_group
+
+    tl.assume(stride_lm > 0)
+    tl.assume(stride_rm > 0)
+    tl.assume(stride_cm > 0)
+
+    BK_PACK: tl.constexpr = BLOCK_SIZE_K // 2
+    BK_SCALE: tl.constexpr = BLOCK_SIZE_K // VEC
+
+    for global_tile in range(pid, total_tiles, NUM_SMS):
+        group_idx = global_tile // tiles_per_group
+        local_tile = global_tile - group_idx * tiles_per_group
+
+        num_pid_in_group = GROUP_SIZE_M * tiles_n
+        swizzle_group = local_tile // num_pid_in_group
+        first_pid_m = swizzle_group * GROUP_SIZE_M
+        group_size_m = min(tiles_m - first_pid_m, GROUP_SIZE_M)
+        pid_m = first_pid_m + ((local_tile % num_pid_in_group) % group_size_m)
+        pid_n = (local_tile % num_pid_in_group) // group_size_m
+        tl.assume(pid_m >= 0)
+        tl.assume(pid_n >= 0)
+
+        m_start = tl.load(go_pad_ptr + group_idx)  # int64 (logical M, multiple of 128)
+        M_g = (tl.load(go_pad_ptr + group_idx + 1) - m_start).to(tl.int32)
+
+        rm = (pid_m * BLOCK_SIZE_M + tl.arange(0, BLOCK_SIZE_M)) % OUT_M
+        rn = (pid_n * BLOCK_SIZE_N + tl.arange(0, BLOCK_SIZE_N)) % OUT_N
+        rk = tl.arange(0, BK_PACK)
+        rks = tl.arange(0, BK_SCALE)
+        rm = tl.max_contiguous(tl.multiple_of(rm, BLOCK_SIZE_M), BLOCK_SIZE_M)
+        rn = tl.max_contiguous(tl.multiple_of(rn, BLOCK_SIZE_N), BLOCK_SIZE_N)
+
+        # packed M offset = m_start // 2 (M_g multiple of 128 -> even)
+        mp0 = (m_start // 2).to(tl.int64)
+        L_BASE = LHS + rm[:, None] * stride_lm + (mp0 + rk[None, :]) * stride_lk
+        R_BASE = RHS + (mp0 + rk[:, None]) * stride_rk + rn[None, :] * stride_rm
+        sk0 = (m_start // VEC).to(tl.int32)
+        LS_BASE = LHS_scale + (sk0 + rks[:, None]) * stride_lsk + rm[None, :] * stride_lsm
+        RS_BASE = RHS_scale + (sk0 + rks[:, None]) * stride_rsk + rn[None, :] * stride_rsm
+
+        acc = tl.zeros((BLOCK_SIZE_M, BLOCK_SIZE_N), dtype=tl.float32)
+        loop_k = M_g // BLOCK_SIZE_K  # padded -> no mask
+        for _ in range(loop_k):
+            l = tl.load(tl.multiple_of(L_BASE, (1, 16)), cache_modifier=CACHE_MODIFIER)  # (BM, BK/2)
+            r = tl.load(tl.multiple_of(R_BASE, (16, 1)), cache_modifier=CACHE_MODIFIER)  # (BK/2, BN)
+            ls = tl.trans(tl.load(LS_BASE))  # (BK/32, BM) -> (BM, BK/32)
+            rs = tl.trans(tl.load(RS_BASE))
+            acc = tl.dot_scaled(l, ls, "e2m1", r, rs, "e2m1", acc)
+            L_BASE += BK_PACK * stride_lk
+            R_BASE += BK_PACK * stride_rk
+            LS_BASE += BK_SCALE * stride_lsk
+            RS_BASE += BK_SCALE * stride_rsk
+
+        c = acc.to(C.type.element_ty)
+        rm_s = pid_m * BLOCK_SIZE_M + tl.arange(0, BLOCK_SIZE_M)
+        rn_s = pid_n * BLOCK_SIZE_N + tl.arange(0, BLOCK_SIZE_N)
+        cmask = (rm_s[:, None] < OUT_M) & (rn_s[None, :] < OUT_N)
+        C_ = C + group_idx.to(tl.int64) * stride_cg + rm_s[:, None] * stride_cm + rn_s[None, :] * stride_cn
+        tl.store(C_, c, cmask)
+
+
+@_scoped_amd_knobs
+def grouped_gemm_mxfp4_variable_k_triton_kernel(
+    lhs, lhs_scale, rhs, rhs_scale, go_pad, OUT_M, OUT_N, G, out_dtype=torch.bfloat16, num_cu=None
+):
+    """C[g] (OUT_M, OUT_N) = lhs[:,g] @ rhs[:,g]^T.
+
+    lhs (OUT_M, M_total/2) fp4-packed, rhs (OUT_N, M_total/2) fp4-packed.
+    go_pad: padded per-group offsets along M (each M_g a multiple of 128).
+    """
+    c = torch.empty((G, OUT_M, OUT_N), dtype=out_dtype, device=lhs.device)
+    ls = lhs_scale.view(torch.uint8)
+    rs = rhs_scale.view(torch.uint8)
+    l_u8 = lhs.view(torch.uint8)
+    r_u8 = rhs.view(torch.uint8)
+    cu = num_cu if num_cu is not None else torch.cuda.get_device_properties(lhs.device).multi_processor_count
+    BM, BN, BK = 256, 128, 128
+    tiles_m = (OUT_M + BM - 1) // BM
+    tiles_n = (OUT_N + BN - 1) // BN
+    GM = 8 if min(tiles_m, tiles_n) < 16 else 4
+    total_tiles = G * tiles_m * tiles_n
+    num_sms = min(total_tiles, cu)
+    chunk = 64 if num_sms >= NUM_XCDS * 64 else 32
+    _grouped_mxfp4_variable_k_gemm_kernel[(num_sms,)](
+        l_u8,
+        r_u8,
+        c,
+        ls,
+        rs,
+        go_pad,
+        G,
+        OUT_M,
+        OUT_N,
+        l_u8.stride(0),
+        l_u8.stride(1),
+        r_u8.stride(0),
+        r_u8.stride(1),
+        c.stride(0),
+        c.stride(1),
+        c.stride(2),
+        ls.stride(0),  # stride_lsm (OUT_M); ls is (OUT_M, M/32)
+        ls.stride(1),  # stride_lsk (M/32 contiguous)
+        rs.stride(0),  # stride_rsm (OUT_N); rs is (OUT_N, M/32)
+        rs.stride(1),  # stride_rsk (M/32 contiguous)
+        BLOCK_SIZE_M=BM,
+        BLOCK_SIZE_N=BN,
+        BLOCK_SIZE_K=BK,
+        GROUP_SIZE_M=GM,
+        NUM_SMS=num_sms,
+        NUM_XCDS=NUM_XCDS,
+        CHUNK_SIZE=chunk,
+        CACHE_MODIFIER=".ca",
+        VEC=VEC_SIZE,
+        num_warps=8,
+        num_stages=3,
+        waves_per_eu=2,
+        matrix_instr_nonkdim=16,
+        kpack=1,
+    )
+    return c
+
+
+# ###########################################################################
+#  Output padding tail
+# ###########################################################################
+
+
+@triton.jit
+def _grouped_gemm_output_tail_kernel(
+    C,
+    group_offs_ptr,
+    G,
+    M_total,
+    N,
+    stride_cm,
+    stride_cn,
+    NUM_SMS: tl.constexpr,
+    BLOCK_SIZE_M: tl.constexpr,
+    BLOCK_SIZE_N: tl.constexpr,
+):
+    covered = tl.load(group_offs_ptr + G).to(tl.int64)
+    m_total = tl.cast(M_total, tl.int64)
+    num_pad = m_total - covered
+    if num_pad <= 0:
+        return
+
+    num_m_blocks = tl.cdiv(num_pad, BLOCK_SIZE_M)
+    num_n_blocks = tl.cdiv(N, BLOCK_SIZE_N)
+    total_tiles = num_m_blocks * num_n_blocks
+
+    zeros = tl.zeros((BLOCK_SIZE_M, BLOCK_SIZE_N), dtype=C.dtype.element_ty)
+
+    pid = tl.program_id(0)
+    for tile in range(pid, total_tiles, NUM_SMS):
+        bm = tile // num_n_blocks
+        bn = tile % num_n_blocks
+        rows = covered + bm * BLOCK_SIZE_M + tl.arange(0, BLOCK_SIZE_M).to(tl.int64)
+        cols = bn * BLOCK_SIZE_N + tl.arange(0, BLOCK_SIZE_N)
+        mask = (rows[:, None] < m_total) & (cols[None, :] < N)
+        C_ = C + rows[:, None] * stride_cm + cols[None, :] * stride_cn
+        tl.store(C_, zeros, mask)
+
+
+def grouped_gemm_output_tail_kernel(
+    out: torch.Tensor,
+    group_offs: torch.Tensor,
+) -> torch.Tensor:
+    """Zero the uncovered padding tail of a grouped GEMM output, in place.
+
+    Clears rows ``[group_offs[-1], out.shape[0])`` -- the rows the persistent
+    forward/dgrad kernel never writes when the output is over-allocated to
+    padded token counts. CPU-sync-free (the covered boundary is read on device)
+    and a no-op when there is no padding.
+    """
+    assert out.ndim == 2, f"expected 2D grouped GEMM output, got {tuple(out.shape)}"
+    M_total, N = out.shape
+    G = group_offs.shape[0] - 1
+    if G <= 0 or M_total == 0 or N == 0:
+        return out
+
+    num_sms = torch.cuda.get_device_properties(out.device).multi_processor_count
+    _grouped_gemm_output_tail_kernel[(num_sms,)](
+        out,
+        group_offs,
+        G,
+        M_total,
+        N,
+        out.stride(0),
+        out.stride(1),
+        NUM_SMS=num_sms,
+        BLOCK_SIZE_M=64,
+        BLOCK_SIZE_N=256,
+        num_warps=4,
+    )
+    return out
diff --git a/transformer_engine/pytorch/triton_kernels/grouped_gemm_mxfp4_impl.py b/transformer_engine/pytorch/triton_kernels/grouped_gemm_mxfp4_impl.py
new file mode 100644
index 0000000000..31c228bb85
--- /dev/null
+++ b/transformer_engine/pytorch/triton_kernels/grouped_gemm_mxfp4_impl.py
@@ -0,0 +1,333 @@
+# Copyright (c) 2026, Advanced Micro Devices, Inc. All rights reserved.
+# License for AMD contributions = MIT. See LICENSE for more information
+
+"""TE-side quantization bridge + dispatch for grouped MXFP4 GEMM (gfx950).
+
+This is the glue between TransformerEngine's MXFP4 quantizer and the ported
+Triton block-scaled grouped kernels in ``grouped_gemm_mxfp4.py``. It reuses
+TE's :class:`MXFP4Quantizer` (rather than porting Primus-Turbo's fused grouped
+dual-quant) to produce the packed FP4 operands the kernels consume, then builds
+the per-group offsets and calls the kernels.
+
+Three grouped ops, paired exactly like the dense MXFP4 recipe (all NT):
+
+    fprop : C   = A_row      @ W_row^T        (contract K,  fwd kernel)
+    dgrad : dA  = gradO_row  @ W_col^T        (contract N,  fwd kernel)
+    wgrad : dW  = gradO_col  @ A_col^T        (contract M,  variable-K kernel)
+
+Operand layout the kernels expect (plain, un-swizzled OCP MXFP4):
+  * data  : E2M1 packed two-per-byte along the contraction axis (``uint8``);
+  * scale : one E8M0 (``uint8``) per 1x32 logical-element block, contiguous
+            along the block axis.
+TE's ``MXFP4Quantizer`` emits exactly this in its row-wise ``_rowwise_data`` /
+``_rowwise_scale_inv`` and col-wise ``_columnwise_data`` / ``_columnwise_scale_inv``
+members when the shuffle / GEMM-swizzle flags are all off. The scale buffers are
+over-allocated (rows to 256, block axis to 8); we slice the live ``[:, :F/32]``
+region and make it contiguous so the kernel sees tight strides.
+
+RHT (Hadamard) on the wgrad operands and the ``main_grad`` beta=1 accumulate
+fusion are deferred to the autograd layer; this module runs plain
+MX (``use_hadamard=False``) and returns a fresh weight-gradient tensor.
+"""
+
+from __future__ import annotations
+
+from typing import List, Optional, Sequence, Tuple
+
+import torch
+import torch.nn.functional as F
+
+from ..tensor.mxfp4_tensor import MXFP4Quantizer
+from ..utils import round_up_to_nearest_multiple
+from .grouped_gemm_mxfp4 import (
+    _is_gfx950,
+    grouped_gemm_mxfp4_triton_kernel,
+    grouped_gemm_mxfp4_variable_k_triton_kernel,
+)
+
+# Logical contraction tile the kernels step by (see grouped_gemm_mxfp4.py).
+BLOCK_SIZE_K = 128
+# One E8M0 scale per 32-element block.
+MXFP4_BLOCK = 32
+
+
+def _quantizer(*, rowwise: bool, columnwise: bool) -> MXFP4Quantizer:
+    """MXFP4 quantizer configured for the plain (un-swizzled) Triton layout."""
+    return MXFP4Quantizer(
+        rowwise=rowwise,
+        columnwise=columnwise,
+        shuffle_rowwise_data=False,
+        shuffle_columnwise_data=False,
+        with_gemm_swizzled_scales=False,
+        use_hadamard=False,
+    )
+
+
+def _prefix_offsets(m_splits: Sequence[int], device: torch.device) -> torch.Tensor:
+    """Tight prefix-sum offsets ``[0, m0, m0+m1, ...]`` as int64 on device."""
+    offs = [0]
+    for m in m_splits:
+        offs.append(offs[-1] + int(m))
+    return torch.tensor(offs, dtype=torch.int64, device=device)
+
+
+def _check_contract(dim: int, name: str) -> None:
+    if dim % BLOCK_SIZE_K != 0:
+        raise ValueError(
+            f"grouped MXFP4 GEMM requires the contraction dim ({name}={dim}) to be a"
+            f" multiple of {BLOCK_SIZE_K}."
+        )
+
+
+def _require_gfx950() -> None:
+    # The kernels use the CDNA4 scaled-FP4 MFMA (tl.dot_scaled e2m1); no other arch has it.
+    if not _is_gfx950():
+        raise RuntimeError(
+            "grouped MXFP4 GEMM requires gfx950 (CDNA4); the current device lacks the"
+            " scaled-FP4 MFMA used by tl.dot_scaled(..., \"e2m1\", ...)."
+        )
+
+
+def _row_operand(x: torch.Tensor) -> Tuple[torch.Tensor, torch.Tensor]:
+    """Row-wise MXFP4: ``x`` [M, F] -> (data [M, F/2] u8, scale [M, F/32] u8)."""
+    M, Feat = x.shape
+    q = _quantizer(rowwise=True, columnwise=False).quantize(x.contiguous())
+    data = q._rowwise_data.view(torch.uint8)
+    scale = q._rowwise_scale_inv[:M, : Feat // MXFP4_BLOCK].contiguous()
+    return data, scale
+
+
+def _col_operand(x: torch.Tensor) -> Tuple[torch.Tensor, torch.Tensor]:
+    """Col-wise MXFP4: ``x`` [M, F] -> (data [F, M/2] u8, scale [F, M/32] u8).
+
+    Requests both directions (matching the production transpose-quant path); the
+    row-wise result is discarded. Columnwise-only quant is a less-exercised
+    config, so we take the safe route here.
+    """
+    M, Feat = x.shape
+    q = _quantizer(rowwise=True, columnwise=True).quantize(x.contiguous())
+    data = q._columnwise_data.view(torch.uint8)
+    scale = q._columnwise_scale_inv[:Feat, : M // MXFP4_BLOCK].contiguous()
+    return data, scale
+
+
+def _col_operand_grouped_padded(
+    x: torch.Tensor,
+    m_splits: Sequence[int],
+) -> Tuple[torch.Tensor, torch.Tensor, torch.Tensor]:
+    """Per-group col-wise MXFP4 with each group's M padded to BLOCK_SIZE_K.
+
+    ``x`` [total_M, F] grouped along M -> (data [F, total_pad/2] u8,
+    scale [F, total_pad/32] u8, go_pad [G+1] int64). Padding rows are zeros, so
+    they contribute nothing to the wgrad contraction over M.
+    """
+    Feat = x.shape[1]
+    datas: List[torch.Tensor] = []
+    scales: List[torch.Tensor] = []
+    offs = [0]
+    start = 0
+    for m in m_splits:
+        m = int(m)
+        xg = x[start : start + m]
+        start += m
+        m_pad = round_up_to_nearest_multiple(m, BLOCK_SIZE_K)
+        if m_pad != m:
+            xg = F.pad(xg, (0, 0, 0, m_pad - m))
+        data, scale = _col_operand(xg)  # (F, m_pad/2), (F, m_pad/32)
+        datas.append(data)
+        scales.append(scale)
+        offs.append(offs[-1] + m_pad)
+    data = torch.cat(datas, dim=1).contiguous()
+    scale = torch.cat(scales, dim=1).contiguous()
+    go_pad = torch.tensor(offs, dtype=torch.int64, device=x.device)
+    return data, scale, go_pad
+
+
+def grouped_gemm_mxfp4_fprop(
+    a: torch.Tensor,
+    weights: List[torch.Tensor],
+    m_splits: Sequence[int],
+    *,
+    out_dtype: torch.dtype = torch.bfloat16,
+    num_cu: Optional[int] = None,
+) -> torch.Tensor:
+    """Grouped MXFP4 forward: ``C[g] = A[g] @ W[g]^T`` (contract K).
+
+    Args:
+        a: [total_M, K] activations, grouped along M by ``m_splits``.
+        weights: list of G per-expert weight tensors, each [N, K].
+        m_splits: per-group token counts (len G).
+
+    Returns:
+        [total_M, N] output in ``out_dtype``.
+    """
+    K = a.shape[1]
+    N = weights[0].shape[0]
+    _require_gfx950()
+    _check_contract(K, "K")
+
+    a_data, a_scale = _row_operand(a)  # (total_M, K/2), (total_M, K/32)
+    b_datas, b_scales = [], []
+    for w in weights:
+        d, s = _row_operand(w)  # (N, K/2), (N, K/32)
+        b_datas.append(d)
+        b_scales.append(s)
+    b_data = torch.stack(b_datas, dim=0)  # (G, N, K/2)
+    b_scale = torch.stack(b_scales, dim=0)  # (G, N, K/32)
+
+    group_offs = _prefix_offsets(m_splits, a.device)
+    return grouped_gemm_mxfp4_triton_kernel(
+        a_data,
+        a_scale,
+        b_data,
+        b_scale,
+        group_offs,
+        N,
+        K,
+        group_offs_out=group_offs,
+        out_dtype=out_dtype,
+        num_cu=num_cu,
+    )
+
+
+def grouped_gemm_mxfp4_dgrad(
+    grad_out: torch.Tensor,
+    weights: List[torch.Tensor],
+    m_splits: Sequence[int],
+    *,
+    out_dtype: torch.dtype = torch.bfloat16,
+    num_cu: Optional[int] = None,
+) -> torch.Tensor:
+    """Grouped MXFP4 dgrad: ``dA[g] = gradO[g] @ W[g]`` (contract N).
+
+    Uses the forward NT kernel with the col-wise weight as the transposed
+    operand: ``B[g] = W_col`` has shape [K, N/2], so the kernel's free dim is K
+    (in_features) and its contraction is N (out_features).
+
+    Args:
+        grad_out: [total_M, N] output gradient, grouped along M by ``m_splits``.
+        weights: list of G per-expert weight tensors, each [N, K].
+
+    Returns:
+        [total_M, K] input gradient in ``out_dtype``.
+    """
+    N = grad_out.shape[1]
+    K = weights[0].shape[1]
+    _require_gfx950()
+    _check_contract(N, "N")
+
+    go_data, go_scale = _row_operand(grad_out)  # (total_M, N/2), (total_M, N/32)
+    b_datas, b_scales = [], []
+    for w in weights:
+        d, s = _col_operand(w)  # (K, N/2), (K, N/32)
+        b_datas.append(d)
+        b_scales.append(s)
+    b_data = torch.stack(b_datas, dim=0)  # (G, K, N/2)
+    b_scale = torch.stack(b_scales, dim=0)  # (G, K, N/32)
+
+    group_offs = _prefix_offsets(m_splits, grad_out.device)
+    # Kernel free dim = b.shape[-2] = K; kernel contraction = b.shape[-1]*2 = N.
+    return grouped_gemm_mxfp4_triton_kernel(
+        go_data,
+        go_scale,
+        b_data,
+        b_scale,
+        group_offs,
+        K,
+        N,
+        group_offs_out=group_offs,
+        out_dtype=out_dtype,
+        num_cu=num_cu,
+    )
+
+
+def grouped_gemm_mxfp4_wgrad(
+    a: torch.Tensor,
+    grad_out: torch.Tensor,
+    m_splits: Sequence[int],
+    *,
+    out_dtype: torch.dtype = torch.bfloat16,
+    num_cu: Optional[int] = None,
+) -> torch.Tensor:
+    """Grouped MXFP4 wgrad: ``dW[g] = gradO[g]^T @ A[g]`` (contract M, variable-K).
+
+    Both operands are col-wise, per-group padded to a BLOCK_SIZE_K-multiple M so
+    the reduction over M needs no masking.
+
+    Args:
+        a: [total_M, K] activations, grouped along M by ``m_splits``.
+        grad_out: [total_M, N] output gradient, grouped along M by ``m_splits``.
+
+    Returns:
+        [G, N, K] weight gradient in ``out_dtype``.
+    """
+    N = grad_out.shape[1]
+    K = a.shape[1]
+    G = len(m_splits)
+    _require_gfx950()
+
+    lhs_data, lhs_scale, go_pad = _col_operand_grouped_padded(grad_out, m_splits)  # (N, Mpad/2)
+    rhs_data, rhs_scale, _ = _col_operand_grouped_padded(a, m_splits)  # (K, Mpad/2)
+
+    return grouped_gemm_mxfp4_variable_k_triton_kernel(
+        lhs_data,
+        lhs_scale,
+        rhs_data,
+        rhs_scale,
+        go_pad,
+        N,
+        K,
+        G,
+        out_dtype=out_dtype,
+        num_cu=num_cu,
+    )
+
+
+class _GroupedGemmMXFP4Func(torch.autograd.Function):
+    """Autograd for grouped MXFP4 linear: ``out[g] = a[g] @ weight[g]^T``.
+
+    Forward/backward pair the MXFP4 recipe exactly (fprop / dgrad / wgrad). Each
+    op quantizes its own operands independently; sharing the quantization across
+    the three passes is a future optimization.
+    """
+
+    @staticmethod
+    def forward(ctx, a, weight, m_splits):
+        out = grouped_gemm_mxfp4_fprop(
+            a, list(torch.unbind(weight, dim=0)), m_splits, out_dtype=a.dtype
+        )
+        ctx.save_for_backward(a, weight)
+        ctx.m_splits = m_splits
+        return out
+
+    @staticmethod
+    def backward(ctx, grad_out):
+        a, weight = ctx.saved_tensors
+        m_splits = ctx.m_splits
+        grad_out = grad_out.contiguous()
+        grad_a = grad_weight = None
+        if ctx.needs_input_grad[0]:
+            grad_a = grouped_gemm_mxfp4_dgrad(
+                grad_out, list(torch.unbind(weight, dim=0)), m_splits, out_dtype=a.dtype
+            )
+        if ctx.needs_input_grad[1]:
+            grad_weight = grouped_gemm_mxfp4_wgrad(a, grad_out, m_splits, out_dtype=weight.dtype)
+        return grad_a, grad_weight, None
+
+
+def grouped_linear_mxfp4(
+    a: torch.Tensor,
+    weight: torch.Tensor,
+    m_splits: Sequence[int],
+) -> torch.Tensor:
+    """Autograd-enabled grouped MXFP4 linear.
+
+    Args:
+        a: [total_M, K] activations, grouped along M by ``m_splits``.
+        weight: [G, N, K] stacked per-expert weights.
+        m_splits: per-group token counts (len G).
+
+    Returns:
+        [total_M, N] output.
+    """
+    return _GroupedGemmMXFP4Func.apply(a, weight, tuple(int(m) for m in m_splits))