diff --git a/tokenspeed-kernel-amd/python/tokenspeed_kernel_amd/ops/gemm/mm_a16w16_gfx950.py b/tokenspeed-kernel-amd/python/tokenspeed_kernel_amd/ops/gemm/mm_a16w16_gfx950.py
index 7d4fb9d34..ec12b96f0 100644
--- a/tokenspeed-kernel-amd/python/tokenspeed_kernel_amd/ops/gemm/mm_a16w16_gfx950.py
+++ b/tokenspeed-kernel-amd/python/tokenspeed_kernel_amd/ops/gemm/mm_a16w16_gfx950.py
@@ -1041,7 +1041,7 @@ def _get_partial_scratch(
 def _check_supported_dense16_shape(M: int, N: int, K: int) -> None:
     if M <= 0 or N <= 0 or K <= 0:
         raise ValueError(
-            "dense16 Gluon GEMM requires positive M/N/K, got " f"M={M}, N={N}, K={K}"
+            f"dense16 Gluon GEMM requires positive M/N/K, got M={M}, N={N}, K={K}"
         )
     if M > MAX_M:
         raise ValueError(f"dense16 Gluon GEMM requires M <= {MAX_M}; got M={M}")
@@ -1057,6 +1057,29 @@ def _check_supported_dense16_shape(M: int, N: int, K: int) -> None:
         )
 
 
+def _resolve_output(
+    A: torch.Tensor,
+    M: int,
+    N: int,
+    out_dtype: torch.dtype,
+    out: torch.Tensor | None,
+    op_name: str,
+) -> torch.Tensor:
+    if out is None:
+        return torch.empty((M, N), device=A.device, dtype=out_dtype)
+    if out.shape != (M, N):
+        raise ValueError(
+            f"{op_name} output must have shape {(M, N)}, got {tuple(out.shape)}"
+        )
+    if out.dtype != out_dtype:
+        raise TypeError(f"{op_name} output dtype must be {out_dtype}, got {out.dtype}")
+    if out.device != A.device:
+        raise ValueError(f"{op_name} output must be on {A.device}, got {out.device}")
+    if out.stride(-1) != 1:
+        raise ValueError(f"{op_name} output must be contiguous in the last dimension")
+    return out
+
+
 def _use_warp_reduce_smallm(M: int, _N: int, K: int) -> bool:
     return (M in (1, 2) and K <= 2048) or (M == 4 and K <= 1024)
 
@@ -1142,6 +1165,7 @@ def gluon_mm_a16w16_warp_reduce_smallm_gfx950(
     out_dtype: torch.dtype,
     *,
     alpha: torch.Tensor | None = None,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor:
     """Compute small-M dense ``A @ B.T`` with one warp per output.
 
@@ -1180,7 +1204,7 @@ def gluon_mm_a16w16_warp_reduce_smallm_gfx950(
         )
     _check_supported_dense16_shape(M, N, K)
 
-    C = torch.empty((M, N), device=A.device, dtype=out_dtype)
+    C = _resolve_output(A, M, N, out_dtype, out, "small-M dense16 warp-reduce GEMM")
     total_outputs = M * N
     grid = (triton.cdiv(total_outputs, WARP_REDUCE_OUTPUTS),)
     _warp_reduce_smallm_kernel[grid](
@@ -1202,7 +1226,7 @@ def gluon_mm_a16w16_warp_reduce_smallm_gfx950(
     )
 
     if alpha is not None:
-        C = C * alpha.to(dtype=C.dtype)
+        C.mul_(alpha.to(device=C.device, dtype=C.dtype))
     return C
 
 
@@ -1212,6 +1236,7 @@ def gluon_mm_a16w16_mfma_lds_smallm_gfx950(
     out_dtype: torch.dtype,
     *,
     alpha: torch.Tensor | None = None,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor:
     """Compute high-K small-M dense ``A @ B.T`` with MFMA and partial-sum reduction.
 
@@ -1224,8 +1249,7 @@ def gluon_mm_a16w16_mfma_lds_smallm_gfx950(
     """
     if A.ndim != 2 or B.ndim != 2:
         raise ValueError(
-            "small-M dense16 MFMA LDS GEMM expects 2D inputs, got "
-            f"{A.ndim=} {B.ndim=}"
+            f"small-M dense16 MFMA LDS GEMM expects 2D inputs, got {A.ndim=} {B.ndim=}"
         )
     if A.dtype not in _SUPPORTED_DTYPES or B.dtype not in _SUPPORTED_DTYPES:
         raise TypeError(
@@ -1261,8 +1285,13 @@ def gluon_mm_a16w16_mfma_lds_smallm_gfx950(
             f"got M={M}, N={N}, K={K}"
         )
 
-    C_full = torch.empty((MFMA_LDS_REDUCE_M, N), device=A.device, dtype=out_dtype)
-    C = C_full[:M, :]
+    if out is not None:
+        _resolve_output(A, M, N, out_dtype, out, "small-M dense16 MFMA LDS GEMM")
+    if out is not None and M == MFMA_LDS_REDUCE_M:
+        C = out
+    else:
+        C_full = torch.empty((MFMA_LDS_REDUCE_M, N), device=A.device, dtype=out_dtype)
+        C = C_full[:M, :]
     split_k = _choose_mfma_lds_split_k(K)
     num_n_tiles = triton.cdiv(N, MFMA_LDS_BLOCK_N)
     total_work = num_n_tiles * split_k
@@ -1310,7 +1339,10 @@ def gluon_mm_a16w16_mfma_lds_smallm_gfx950(
     )
 
     if alpha is not None:
-        C = C * alpha.to(dtype=C.dtype)
+        C.mul_(alpha.to(device=C.device, dtype=C.dtype))
+    if out is not None and C is not out:
+        out.copy_(C)
+        return out
     return C
 
 
@@ -1320,6 +1352,7 @@ def gluon_mm_a16w16_mfma_lds_mediumm_gfx950(
     out_dtype: torch.dtype,
     *,
     alpha: torch.Tensor | None = None,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor:
     """Compute tuned medium-M dense ``A @ B.T`` with MFMA/LDS tiling.
 
@@ -1364,7 +1397,7 @@ def gluon_mm_a16w16_mfma_lds_mediumm_gfx950(
         )
 
     block_m, block_n, block_k, warps_m, warps_n, num_buffers = config
-    C = torch.empty((M, N), device=A.device, dtype=out_dtype)
+    C = _resolve_output(A, M, N, out_dtype, out, "medium-M dense16 MFMA LDS GEMM")
     grid = (triton.cdiv(M, block_m) * triton.cdiv(N, block_n),)
     _mfma_lds_mediumm_kernel[grid](
         A,
@@ -1391,7 +1424,7 @@ def gluon_mm_a16w16_mfma_lds_mediumm_gfx950(
     )
 
     if alpha is not None:
-        C = C * alpha.to(dtype=C.dtype)
+        C.mul_(alpha.to(device=C.device, dtype=C.dtype))
     return C
 
 
@@ -1401,6 +1434,7 @@ def gluon_mm_a16w16_gfx950(
     out_dtype: torch.dtype,
     *,
     alpha: torch.Tensor | None = None,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor | None:
     """Compute dense16 GEMM ``A @ B.T`` on gfx950 when supported.
 
@@ -1444,6 +1478,7 @@ def gluon_mm_a16w16_gfx950(
             B,
             out_dtype,
             alpha=alpha,
+            out=out,
         )
     if M > MAX_M:
         return None
@@ -1454,6 +1489,7 @@ def gluon_mm_a16w16_gfx950(
             B,
             out_dtype,
             alpha=alpha,
+            out=out,
         )
     if _use_mfma_lds_smallm(M, N, K):
         return gluon_mm_a16w16_mfma_lds_smallm_gfx950(
@@ -1461,6 +1497,7 @@ def gluon_mm_a16w16_gfx950(
             B,
             out_dtype,
             alpha=alpha,
+            out=out,
         )
     if _use_mfma_lds_mediumm(M, N, K):
         return gluon_mm_a16w16_mfma_lds_mediumm_gfx950(
@@ -1468,5 +1505,6 @@ def gluon_mm_a16w16_gfx950(
             B,
             out_dtype,
             alpha=alpha,
+            out=out,
         )
     return None
diff --git a/tokenspeed-kernel-amd/python/tokenspeed_kernel_amd/ops/gemm/mm_a16w16_largem_gfx950.py b/tokenspeed-kernel-amd/python/tokenspeed_kernel_amd/ops/gemm/mm_a16w16_largem_gfx950.py
index e101a9646..713781147 100644
--- a/tokenspeed-kernel-amd/python/tokenspeed_kernel_amd/ops/gemm/mm_a16w16_largem_gfx950.py
+++ b/tokenspeed-kernel-amd/python/tokenspeed_kernel_amd/ops/gemm/mm_a16w16_largem_gfx950.py
@@ -394,12 +394,41 @@ def _supports_largem_shape(M: int, N: int, K: int) -> bool:
     )
 
 
+def _resolve_largem_output(
+    A: torch.Tensor,
+    M: int,
+    N: int,
+    out_dtype: torch.dtype,
+    out: torch.Tensor | None,
+) -> torch.Tensor:
+    if out is None:
+        return torch.empty((M, N), device=A.device, dtype=out_dtype)
+    if out.shape != (M, N):
+        raise ValueError(
+            f"large-M dense16 output must have shape {(M, N)}, got {tuple(out.shape)}"
+        )
+    if out.dtype != out_dtype:
+        raise TypeError(
+            f"large-M dense16 output dtype must be {out_dtype}, got {out.dtype}"
+        )
+    if out.device != A.device:
+        raise ValueError(
+            f"large-M dense16 output must be on {A.device}, got {out.device}"
+        )
+    if out.stride(-1) != 1:
+        raise ValueError(
+            "large-M dense16 output must be contiguous in the last dimension"
+        )
+    return out
+
+
 def gluon_mm_a16w16_largem_gfx950(
     A: torch.Tensor,
     B: torch.Tensor,
     out_dtype: torch.dtype,
     *,
     alpha: torch.Tensor | None = None,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor | None:
     """Compute large, aligned dense16 ``A @ B.T`` with an 8-wave MFMA/LDS tile."""
     if A.ndim != 2 or B.ndim != 2:
@@ -416,7 +445,7 @@ def gluon_mm_a16w16_largem_gfx950(
     if K_b != K or not _supports_largem_shape(M, N, K):
         return None
 
-    C = torch.empty((M, N), device=A.device, dtype=out_dtype)
+    C = _resolve_largem_output(A, M, N, out_dtype, out)
     grid_m = triton.cdiv(M, LARGEM_BLOCK_M)
     grid_n = triton.cdiv(N, LARGEM_BLOCK_N)
     grid_mn = grid_m * grid_n
@@ -446,5 +475,5 @@ def gluon_mm_a16w16_largem_gfx950(
     )
 
     if alpha is not None:
-        C = C * alpha.to(dtype=C.dtype)
+        C.mul_(alpha.to(device=C.device, dtype=C.dtype))
     return C
diff --git a/tokenspeed-kernel-amd/test/ops/test_mm_a16w16_gfx950.py b/tokenspeed-kernel-amd/test/ops/test_mm_a16w16_gfx950.py
index 31367a30d..cf4f753d1 100644
--- a/tokenspeed-kernel-amd/test/ops/test_mm_a16w16_gfx950.py
+++ b/tokenspeed-kernel-amd/test/ops/test_mm_a16w16_gfx950.py
@@ -76,6 +76,39 @@ def test_dense16_kernel_variant_correctness(
     torch.testing.assert_close(out, torch.mm(a, b.T), atol=1e-2, rtol=1e-2)
 
 
+@pytest.mark.parametrize("kernel,shape", _CORRECTNESS_CASES)
+def test_dense16_kernel_variant_writes_strided_out(
+    kernel, shape: tuple[int, int, int]
+) -> None:
+    torch.manual_seed(0)
+    dtype = torch.bfloat16
+    m, n, k = shape
+    a = torch.randn((m, k), device="cuda", dtype=dtype) * 0.25
+    b = torch.randn((n, k), device="cuda", dtype=dtype) * 0.25
+    backing = torch.empty((m, n + 17), device="cuda", dtype=dtype)
+    out = backing[:, :n]
+
+    actual = kernel(a, b, dtype, out=out)
+
+    assert actual is out
+    torch.testing.assert_close(out, torch.mm(a, b.T), atol=1e-2, rtol=1e-2)
+
+
+def test_splitk_smallm_out_handles_padded_reducer_rows() -> None:
+    torch.manual_seed(0)
+    dtype = torch.bfloat16
+    m, n, k = 2, 256, 2048
+    a = torch.randn((m, k), device="cuda", dtype=dtype) * 0.25
+    b = torch.randn((n, k), device="cuda", dtype=dtype) * 0.25
+    backing = torch.empty((m, n + 17), device="cuda", dtype=dtype)
+    out = backing[:, :n]
+
+    actual = gluon_mm_a16w16_mfma_lds_smallm_gfx950(a, b, dtype, out=out)
+
+    assert actual is out
+    torch.testing.assert_close(out, torch.mm(a, b.T), atol=1e-2, rtol=1e-2)
+
+
 def test_use_warp_reduce_covers_small_k_decode_shapes() -> None:
     assert _use_warp_reduce_smallm(1, 1280, 1024)
     assert _use_warp_reduce_smallm(2, 2560, 2048)
diff --git a/tokenspeed-kernel/python/tokenspeed_kernel/__init__.py b/tokenspeed-kernel/python/tokenspeed_kernel/__init__.py
index 9b7abe688..b5e7edb19 100644
--- a/tokenspeed-kernel/python/tokenspeed_kernel/__init__.py
+++ b/tokenspeed-kernel/python/tokenspeed_kernel/__init__.py
@@ -43,7 +43,7 @@
     rel_mha_plan,
     rel_mha_prefill,
 )
-from tokenspeed_kernel.ops.gemm import mm
+from tokenspeed_kernel.ops.gemm import bmm, mm
 from tokenspeed_kernel.ops.moe import moe_apply, moe_plan, moe_process_weights
 from tokenspeed_kernel.ops.quantization import (
     quantize_fp8,
@@ -60,6 +60,7 @@
     # exceptions
     "NoKernelFoundError",
     # gemm
+    "bmm",
     "mm",
     # attention
     "mha_plan",
diff --git a/tokenspeed-kernel/python/tokenspeed_kernel/benchmark/throughput.py b/tokenspeed-kernel/python/tokenspeed_kernel/benchmark/throughput.py
index 70c111bc1..d96067590 100644
--- a/tokenspeed-kernel/python/tokenspeed_kernel/benchmark/throughput.py
+++ b/tokenspeed-kernel/python/tokenspeed_kernel/benchmark/throughput.py
@@ -54,6 +54,21 @@ def gemm_mm_bytes(M: int, N: int, K: int, dtype: torch.dtype) -> int:
         element_size = _dtype_nbytes(dtype)
         return (M * K + K * N + M * N) * element_size
 
+    @staticmethod
+    def gemm_bmm_flops(batch: int, M: int, N: int, K: int) -> int:
+        return batch * ThroughputCalculator.gemm_mm_flops(M, N, K)
+
+    @staticmethod
+    def gemm_bmm_bytes(
+        batch: int,
+        M: int,
+        N: int,
+        K: int,
+        dtype: torch.dtype,
+    ) -> int:
+        element_size = _dtype_nbytes(dtype)
+        return (batch * M * K + batch * N * K + batch * M * N) * element_size
+
     @staticmethod
     def attn_flops(
         batch: int,
@@ -105,6 +120,28 @@ def compute(
             bandwidth = bytes_moved / seconds / 1e9
             return tflops, bandwidth
 
+        if op_family == "gemm" and op_mode == "bmm":
+            batch = _shape_int(shape_params, "B", "batch")
+            M = _shape_int(shape_params, "M")
+            N = _shape_int(shape_params, "N")
+            K = _shape_int(shape_params, "K")
+            if any(value is None for value in (batch, M, N, K)):
+                return None, None
+
+            flops = ThroughputCalculator.gemm_bmm_flops(batch, M, N, K)
+            bytes_moved = ThroughputCalculator.gemm_bmm_bytes(
+                batch,
+                M,
+                N,
+                K,
+                dtype,
+            )
+            seconds = latency_us * 1e-6
+
+            tflops = flops / seconds / 1e12
+            bandwidth = bytes_moved / seconds / 1e9
+            return tflops, bandwidth
+
         if op_family in {"attention", "attn"} and op_mode == "decode":
             batch = _shape_int(
                 shape_params,
diff --git a/tokenspeed-kernel/python/tokenspeed_kernel/numerics/gemm.py b/tokenspeed-kernel/python/tokenspeed_kernel/numerics/gemm.py
index e2970412c..5a41b1646 100644
--- a/tokenspeed-kernel/python/tokenspeed_kernel/numerics/gemm.py
+++ b/tokenspeed-kernel/python/tokenspeed_kernel/numerics/gemm.py
@@ -104,6 +104,14 @@ def tolerance(
 
 
 class GemmInputGenerator(InputGenerator):
+    def _trait_value(self, name: str, default: str) -> str:
+        value = self.traits.get(name)
+        if isinstance(value, str):
+            return value
+        if isinstance(value, (frozenset, set)) and len(value) == 1:
+            return str(next(iter(value)))
+        return default
+
     def _generate_value(self, shape: tuple[int, ...], dtype) -> torch.Tensor:
         values = torch.randn(
             *shape,
@@ -142,6 +150,7 @@ def _scale_for_format(
         tensor_format: TensorFormat | None,
         role: str,
         *,
+        batch_shape: tuple[int, ...],
         M: int,
         N: int,
         K: int,
@@ -159,47 +168,76 @@ def _scale_for_format(
             block_n, block_k = block_size
             k_tiles = math.ceil(K / block_k)
             if role == "a":
-                return self._generate_scales((M, k_tiles), scale.storage_dtype)
+                return self._generate_scales(
+                    (*batch_shape, M, k_tiles), scale.storage_dtype
+                )
             if role == "b":
                 n_tiles = math.ceil(N / block_n)
-                return self._generate_scales((n_tiles, k_tiles), scale.storage_dtype)
+                return self._generate_scales(
+                    (*batch_shape, n_tiles, k_tiles), scale.storage_dtype
+                )
 
         if scale.granularity == "channel":
             return self._generate_scales(
-                (M,) if role == "a" else (N,),
+                (*batch_shape, M) if role == "a" else (*batch_shape, N),
                 scale.storage_dtype,
             )
 
         return self._generate_scales((1,), scale.storage_dtype)
 
-    def generate(
+    def _a_shape(
+        self,
+        *,
+        batch_shape: tuple[int, ...],
+        M: int,
+        K: int,
+        layout: str,
+    ) -> tuple[int, ...]:
+        if layout in {"KM", "BKM"}:
+            return (*batch_shape, K, M)
+        return (*batch_shape, M, K)
+
+    def _b_shape(
+        self,
+        *,
+        batch_shape: tuple[int, ...],
+        N: int,
+        K: int,
+        layout: str,
+    ) -> tuple[int, ...]:
+        if layout in {"KN", "BKN"}:
+            return (*batch_shape, K, N)
+        return (*batch_shape, N, K)
+
+    def _generate_for_shape(
         self,
+        *,
+        batch_shape: tuple[int, ...],
         M: int,
         N: int,
         K: int,
+        a_layout: str,
+        b_layout: str,
     ) -> dict[str, Any]:
-        a_layout = self.traits.get("a_layout")
-        b_layout = self.traits.get("b_layout")
         a_format = self._format("a")
         b_format = self._format("b")
+
         a_dtype = a_format.storage_dtype if a_format is not None else self.dtype
         b_dtype = b_format.storage_dtype if b_format is not None else self.dtype
-
-        A = (
-            self._generate_value((K, M), a_dtype)
-            if a_layout == {"KM"}
-            else self._generate_value((M, K), a_dtype)
+        A = self._generate_value(
+            self._a_shape(batch_shape=batch_shape, M=M, K=K, layout=a_layout),
+            a_dtype,
         )
-        B = (
-            self._generate_value((K, N), b_dtype)
-            if b_layout == {"KN"}
-            else self._generate_value((N, K), b_dtype)
+        B = self._generate_value(
+            self._b_shape(batch_shape=batch_shape, N=N, K=K, layout=b_layout),
+            b_dtype,
         )
 
         block_size = self._block_size(a_format, b_format)
         A_scales = self._scale_for_format(
             a_format,
             "a",
+            batch_shape=batch_shape,
             M=M,
             N=N,
             K=K,
@@ -208,27 +246,62 @@ def generate(
         B_scales = self._scale_for_format(
             b_format,
             "b",
+            batch_shape=batch_shape,
             M=M,
             N=N,
             K=K,
             block_size=block_size,
         )
 
-        out_dtype = torch.bfloat16
-        alpha = None
-
         return {
             "A": A,
             "B": B,
             "A_scales": A_scales,
             "B_scales": B_scales,
-            "out_dtype": out_dtype,
-            "alpha": alpha,
+            "out_dtype": torch.bfloat16,
+            "alpha": None,
             "block_size": block_size,
         }
 
+    def generate(
+        self,
+        M: int,
+        N: int,
+        K: int,
+    ) -> dict[str, Any]:
+        return self._generate_for_shape(
+            batch_shape=(),
+            M=M,
+            N=N,
+            K=K,
+            a_layout=self._trait_value("a_layout", "MK"),
+            b_layout=self._trait_value("b_layout", "NK"),
+        )
+
+
+class BmmInputGenerator(GemmInputGenerator):
+    def generate(
+        self,
+        M: int,
+        N: int,
+        K: int,
+        B: int | None = None,
+        batch: int | None = None,
+        **_: Any,
+    ) -> dict[str, Any]:
+        batch_size = B if B is not None else batch if batch is not None else 1
+        return self._generate_for_shape(
+            batch_shape=(batch_size,),
+            M=M,
+            N=N,
+            K=K,
+            a_layout="BMK",
+            b_layout="BNK",
+        )
+
 
 set_input_generator("gemm", "mm", GemmInputGenerator)
+set_input_generator("gemm", "bmm", BmmInputGenerator)
 
 # ---------------------------------------------------------------------------
 # Shape Presets
@@ -263,5 +336,24 @@ def generate(
 ]
 
 
+GEMM_BMM_STANDARD_SHAPES: list[dict[str, Any]] = [
+    {"batch": 1, "M": 4, "N": 32, "K": 64},
+    {"batch": 4, "M": 8, "N": 128, "K": 64},
+    {"batch": 8, "M": 2, "N": 128, "K": 128},
+    {"batch": 4, "M": 16, "N": 256, "K": 256},
+]
+
+
+GEMM_BMM_BENCHMARK_SHAPES: list[dict[str, Any]] = [
+    {"batch": 32, "M": 1, "N": 4096, "K": 4096},
+    {"batch": 8, "M": 16, "N": 4096, "K": 4096},
+    {"batch": 4, "M": 64, "N": 4096, "K": 4096},
+    {"batch": 32, "M": 1, "N": 512, "K": 512},
+    {"batch": 8, "M": 16, "N": 512, "K": 512},
+]
+
+
 set_standard_shapes("gemm", "mm", GEMM_MM_STANDARD_SHAPES)
 set_benchmark_shapes("gemm", "mm", GEMM_MM_BENCHMARK_SHAPES)
+set_standard_shapes("gemm", "bmm", GEMM_BMM_STANDARD_SHAPES)
+set_benchmark_shapes("gemm", "bmm", GEMM_BMM_BENCHMARK_SHAPES)
diff --git a/tokenspeed-kernel/python/tokenspeed_kernel/numerics/reference/gemm.py b/tokenspeed-kernel/python/tokenspeed_kernel/numerics/reference/gemm.py
index 0b256ee9e..24661aadd 100644
--- a/tokenspeed-kernel/python/tokenspeed_kernel/numerics/reference/gemm.py
+++ b/tokenspeed-kernel/python/tokenspeed_kernel/numerics/reference/gemm.py
@@ -38,17 +38,60 @@
     storage_dtype=torch.float32,
     granularity="tensor",
 )
+_FP8_CHANNEL_SCALE = ScaleFormat(
+    storage_dtype=torch.float32,
+    granularity="channel",
+)
 _MXFP8_FORMAT_SIGNATURES = format_signatures(
     ("a", "b"), "mxfp8", {fp8_dtype}, scale=_FP8_BLOCK_SCALE
 )
 _FP8_TENSOR_FORMAT_SIGNATURES = format_signatures(
     ("a", "b"), "scaled-fp8", {fp8_dtype}, scale=_FP8_TENSOR_SCALE
 )
+_FP8_SCALED_BMM_FORMAT_SIGNATURES = _FP8_TENSOR_FORMAT_SIGNATURES | format_signatures(
+    ("a", "b"), "scaled-fp8", {fp8_dtype}, scale=_FP8_CHANNEL_SCALE
+)
 _DENSE_GEMM_FORMAT_SIGNATURES = format_signatures(
     ("a", "b"), "dense", {torch.bfloat16, torch.float16, torch.float32}
 )
 
 
+def _as_baddbmm_alpha(alpha: torch.Tensor | float | int | None) -> float | int:
+    if alpha is None:
+        return 1
+    if isinstance(alpha, torch.Tensor):
+        if alpha.numel() != 1:
+            raise ValueError(
+                f"torch_bmm alpha expects a scalar tensor, got {alpha.shape}"
+            )
+        return alpha.item()
+    return alpha
+
+
+def _reference_mxfp8_quantize(
+    A: torch.Tensor,
+    *,
+    block_k: int,
+) -> tuple[torch.Tensor, torch.Tensor]:
+    k_tiles = math.ceil(A.shape[-1] / block_k)
+    qA = torch.empty(A.shape, device=A.device, dtype=fp8_dtype)
+    A_scales = torch.empty(
+        (*A.shape[:-1], k_tiles),
+        device=A.device,
+        dtype=torch.float32,
+    )
+    fp8_max = float(torch.finfo(fp8_dtype).max)
+    min_scale = torch.finfo(torch.float32).tiny
+    for tile_idx in range(k_tiles):
+        start = tile_idx * block_k
+        end = min(start + block_k, A.shape[-1])
+        tile = A[..., start:end].float()
+        scale = (tile.abs().amax(dim=-1) / fp8_max).clamp_min(min_scale)
+        A_scales[..., tile_idx] = scale
+        qA[..., start:end] = (tile / scale.unsqueeze(-1)).to(fp8_dtype)
+    return qA, A_scales
+
+
 @register_kernel(
     "gemm",
     "mm",
@@ -68,11 +111,12 @@ def torch_mm_fp8_blockscale(
     *,
     alpha: torch.Tensor | None = None,
     block_size: list[int] | None = None,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor:
     assert block_size is not None, "block_size is required for mxfp8 reference"
-    assert (
-        A_scales is not None and B_scales is not None
-    ), "A_scales and B_scales are required for mxfp8 reference"
+    if A_scales is None:
+        A, A_scales = _reference_mxfp8_quantize(A, block_k=block_size[1])
+    assert B_scales is not None, "B_scales is required for mxfp8 reference"
     assert A.ndim == 2 and B.ndim == 2, f"Expected 2D inputs, got {A.ndim=} {B.ndim=}"
 
     M, K = A.shape
@@ -101,7 +145,13 @@ def torch_mm_fp8_blockscale(
 
     if alpha is not None:
         output = output * alpha.float()
-    return output.to(out_dtype)
+    output = output.to(out_dtype)
+    if out is not None:
+        # Reference expressions materialize through fp32 dequantization first,
+        # so support pre-allocated output by copying into caller storage.
+        out.copy_(output)
+        return out
+    return output
 
 
 @register_kernel(
@@ -125,6 +175,7 @@ def torch_mm_fp8_scaled_mnk(
     *,
     alpha: torch.Tensor | None = None,
     block_size: list[int] | None = None,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor:
     assert block_size is None, "block_size is not supported for fp8 scaled reference"
     assert (
@@ -143,7 +194,13 @@ def torch_mm_fp8_scaled_mnk(
 
     if alpha is not None:
         output = output * alpha.float()
-    return output.to(out_dtype)
+    output = output.to(out_dtype)
+    if out is not None:
+        # Reference expressions materialize through fp32 dequantization first,
+        # so support pre-allocated output by copying into caller storage.
+        out.copy_(output)
+        return out
+    return output
 
 
 @register_kernel(
@@ -167,6 +224,7 @@ def torch_mm_fp8_scaled_nkm(
     *,
     alpha: torch.Tensor | None = None,
     block_size: list[int] | None = None,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor:
     assert block_size is None, "block_size is not supported for fp8 scaled reference"
     assert (
@@ -183,7 +241,13 @@ def torch_mm_fp8_scaled_nkm(
 
     if alpha is not None:
         output = output * alpha.float()
-    return output.to(out_dtype)
+    output = output.to(out_dtype)
+    if out is not None:
+        # Reference expressions materialize through fp32 dequantization first,
+        # so support pre-allocated output by copying into caller storage.
+        out.copy_(output)
+        return out
+    return output
 
 
 @register_kernel(
@@ -206,7 +270,22 @@ def torch_mm(
     alpha: torch.Tensor | None = None,
     block_size: list[int] | None = None,
     bias: torch.Tensor | None = None,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor:
+    if out is not None:
+        if out_dtype != A.dtype:
+            raise ValueError(
+                f"torch_mm out= requires out_dtype {A.dtype}, got {out_dtype}"
+            )
+        # F.linear has no portable out= form, so write the GEMM natively and
+        # apply the epilogue in place.
+        output = torch.mm(A, B.T, out=out)
+        if alpha is not None:
+            output.mul_(alpha.to(dtype=output.dtype))
+        if bias is not None:
+            output.add_(bias.to(dtype=output.dtype))
+        return output
+
     if alpha is None:
         # F.linear fuses the bias add inside the GEMM epilogue.
         output = F.linear(A, B, bias)
@@ -216,3 +295,217 @@ def torch_mm(
         if bias is not None:
             output = output + bias.to(dtype=output.dtype)
     return output.to(out_dtype)
+
+
+@register_kernel(
+    "gemm",
+    "bmm",
+    name="torch_bmm_fp8_blockscale",
+    solution="reference",
+    signatures=_MXFP8_FORMAT_SIGNATURES,
+    traits={},
+    priority=Priority.PORTABLE + 2,
+    tags={"portability"},
+)
+def torch_bmm_fp8_blockscale(
+    A: torch.Tensor,
+    B: torch.Tensor,
+    A_scales: torch.Tensor | None,
+    B_scales: torch.Tensor | None,
+    out_dtype: torch.dtype,
+    *,
+    alpha: torch.Tensor | None = None,
+    block_size: list[int] | None = None,
+    out: torch.Tensor | None = None,
+) -> torch.Tensor:
+    assert block_size is not None, "block_size is required for mxfp8 reference"
+    if A_scales is None:
+        A, A_scales = _reference_mxfp8_quantize(A, block_k=block_size[1])
+    assert B_scales is not None, "B_scales is required for mxfp8 reference"
+    assert A.ndim == 3 and B.ndim == 3, f"Expected 3D inputs, got {A.ndim=} {B.ndim=}"
+
+    batch, M, K = A.shape
+    B_batch, N, K_b = B.shape
+    assert B_batch == batch, f"Expected matching batch dims, got {A.shape=} {B.shape=}"
+    assert K_b == K, f"Expected B in [B, N, K] layout, got shape={tuple(B.shape)}"
+
+    block_n, block_k = block_size
+    k_tiles = math.ceil(K / block_k)
+    n_tiles = math.ceil(N / block_n)
+    assert A_scales.shape == (batch, M, k_tiles), (
+        f"A_scales shape mismatch: expected {(batch, M, k_tiles)}, "
+        f"got {tuple(A_scales.shape)}"
+    )
+    assert B_scales.shape == (batch, n_tiles, k_tiles), (
+        f"B_scales shape mismatch: expected {(batch, n_tiles, k_tiles)}, "
+        f"got {tuple(B_scales.shape)}"
+    )
+
+    A_scaled = A_scales.float().repeat_interleave(block_k, dim=2)[:, :, :K]
+    B_scaled = (
+        B_scales.float()
+        .repeat_interleave(block_n, dim=1)
+        .repeat_interleave(block_k, dim=2)[:, :N, :K]
+    )
+    output = torch.bmm(A.float() * A_scaled, (B.float() * B_scaled).transpose(1, 2))
+
+    if alpha is not None:
+        output = output * alpha.float()
+    output = output.to(out_dtype)
+    if out is not None:
+        # Reference expressions materialize through fp32 dequantization first,
+        # so support pre-allocated output by copying into caller storage.
+        out.copy_(output)
+        return out
+    return output
+
+
+def _bmm_scaled_fp8_scale(
+    scale: torch.Tensor | None,
+    *,
+    batch: int,
+    rows: int,
+    name: str,
+) -> float | torch.Tensor:
+    assert scale is not None, f"{name} is required for fp8 scaled reference"
+    if scale.shape == (1,):
+        return float(scale.item())
+    if scale.shape == (batch, rows):
+        return scale.float().unsqueeze(-1)
+    raise AssertionError(
+        f"{name} shape mismatch: expected (1,) or {(batch, rows)}, "
+        f"got {tuple(scale.shape)}"
+    )
+
+
+@register_kernel(
+    "gemm",
+    "bmm",
+    name="torch_bmm_fp8_scaled",
+    solution="reference",
+    signatures=_FP8_SCALED_BMM_FORMAT_SIGNATURES,
+    traits={},
+    priority=Priority.PORTABLE,
+    tags={"portability"},
+)
+def torch_bmm_fp8_scaled(
+    A: torch.Tensor,
+    B: torch.Tensor,
+    A_scales: torch.Tensor | None,
+    B_scales: torch.Tensor | None,
+    out_dtype: torch.dtype,
+    *,
+    alpha: torch.Tensor | None = None,
+    block_size: list[int] | None = None,
+    out: torch.Tensor | None = None,
+) -> torch.Tensor:
+    assert block_size is None, "block_size is not supported for fp8 scaled reference"
+    assert A.ndim == 3 and B.ndim == 3, f"Expected 3D inputs, got {A.ndim=} {B.ndim=}"
+
+    batch, M, K = A.shape
+    B_batch, N, K_b = B.shape
+    assert B_batch == batch, f"Expected matching batch dims, got {A.shape=} {B.shape=}"
+    assert K_b == K, f"Expected B in [B, N, K] layout, got shape={tuple(B.shape)}"
+
+    A_scale = _bmm_scaled_fp8_scale(A_scales, batch=batch, rows=M, name="A_scales")
+    B_scale = _bmm_scaled_fp8_scale(B_scales, batch=batch, rows=N, name="B_scales")
+    output = torch.bmm(A.float() * A_scale, (B.float() * B_scale).transpose(1, 2))
+
+    if alpha is not None:
+        output = output * alpha.float()
+    output = output.to(out_dtype)
+    if out is not None:
+        # Reference expressions materialize through fp32 dequantization first,
+        # so support pre-allocated output by copying into caller storage.
+        out.copy_(output)
+        return out
+    return output
+
+
+@register_kernel(
+    "gemm",
+    "bmm",
+    name="torch_bmm",
+    solution="reference",
+    signatures=_DENSE_GEMM_FORMAT_SIGNATURES,
+    priority=Priority.PORTABLE + 3,
+    tags={"determinism", "portability"},
+)
+def torch_bmm(
+    A: torch.Tensor,
+    B: torch.Tensor,
+    A_scales: torch.Tensor | None,
+    B_scales: torch.Tensor | None,
+    out_dtype: torch.dtype,
+    *,
+    alpha: torch.Tensor | None = None,
+    block_size: list[int] | None = None,
+    bias: torch.Tensor | None = None,
+    out: torch.Tensor | None = None,
+) -> torch.Tensor:
+    if A_scales is not None:
+        raise ValueError("A_scales are not supported for dense reference BMM")
+    if B_scales is not None:
+        raise ValueError("B_scales are not supported for dense reference BMM")
+    if block_size is not None:
+        raise ValueError("block_size is not supported for dense reference BMM")
+    if A.ndim != 3:
+        raise ValueError(f"torch_bmm expects A=[B, M, K], got {A.shape}")
+    if B.ndim != 3:
+        raise ValueError(f"torch_bmm expects B=[B, N, K], got {B.shape}")
+    if A.shape[0] != B.shape[0]:
+        raise ValueError(f"torch_bmm batch mismatch: {A.shape=} {B.shape=}")
+    if A.shape[2] != B.shape[2]:
+        raise ValueError(f"torch_bmm K mismatch: {A.shape=} {B.shape=}")
+    if out is not None and out_dtype != A.dtype:
+        raise ValueError(
+            f"torch_bmm out= requires out_dtype {A.dtype}, got {out_dtype}"
+        )
+    if bias is not None and out_dtype == A.dtype:
+        bias = bias.to(dtype=A.dtype)
+        # torch.baddbmm is the batched equivalent of the dense F.linear path:
+        # it combines bias, alpha, matmul, and out= for native output dtype.
+        if bias.ndim == 1:
+            bias_view = bias.view(1, 1, -1)
+        elif bias.ndim == 2:
+            bias_view = bias.view(bias.shape[0], 1, bias.shape[1])
+        else:
+            raise ValueError(
+                f"torch_bmm bias expects shape [N] or [B, N], got {bias.shape}"
+            )
+        if out is not None:
+            return torch.baddbmm(
+                bias_view,
+                A,
+                B.transpose(1, 2),
+                alpha=_as_baddbmm_alpha(alpha),
+                out=out,
+            )
+        return torch.baddbmm(
+            bias_view,
+            A,
+            B.transpose(1, 2),
+            alpha=_as_baddbmm_alpha(alpha),
+        )
+
+    output = torch.bmm(A, B.transpose(1, 2), out=out)
+    if alpha is not None:
+        output.mul_(alpha.to(dtype=output.dtype))
+    if bias is not None:
+        bias = bias.to(dtype=output.dtype)
+        # Match mm bias broadcasting for either a shared [N] vector or a
+        # per-batch [B, N] bias matrix.
+        if bias.ndim == 1:
+            bias_view = bias.view(1, 1, -1)
+        elif bias.ndim == 2:
+            bias_view = bias.view(bias.shape[0], 1, bias.shape[1])
+        else:
+            raise ValueError(
+                f"torch_bmm bias expects shape [N] or [B, N], got {bias.shape}"
+            )
+        output.add_(bias_view)
+    if output.dtype != out_dtype:
+        # torch.bmm only writes the input dtype; non-native output dtypes are
+        # represented as a reference cast after the native epilogue.
+        output = output.to(out_dtype)
+    return output
diff --git a/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/__init__.py b/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/__init__.py
index 14157230c..37b0e505d 100644
--- a/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/__init__.py
+++ b/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/__init__.py
@@ -32,6 +32,7 @@
 import torch
 from tokenspeed_kernel.platform import ArchVersion, Platform
 from tokenspeed_kernel.profiling import ShapeCapture, kernel_scope
+from tokenspeed_kernel.registry import KernelRegistry
 from tokenspeed_kernel.selection import select_kernel
 from tokenspeed_kernel.signature import (
     ScaleFormat,
@@ -42,16 +43,17 @@
 
 logger = logging.getLogger(__name__)
 
-__all__ = ["mm"]
+__all__ = ["bmm", "mm"]
 
 _platform = Platform.get()
 _fp8_dtype = _platform.fp8e4m3fn.dtype
 
-# Kernels that natively fuse a bias-vector add inside their GEMM kernel.
-# For any kernel not listed here, ``mm`` applies the bias with a post-GEMM
+# Kernels that accept and own bias application inside their GEMM wrapper.
+# For any kernel not listed here, dispatch applies the bias with a post-GEMM
 # add instead of passing it to the kernel.
 _KERNELS_WITH_FUSED_BIAS: frozenset[str] = frozenset(
     {
+        "torch_bmm",
         "torch_mm",
         "triton_mm_fp8_scaled",
     }
@@ -240,11 +242,34 @@ def ensure_row_major_scales(
         raise ValueError(f"No online quantization defined for kernel {kernel_name!r}")
 
 
+def _kernel_handles_online_mxfp8(kernel_name: str) -> bool:
+    spec = KernelRegistry.get().get_by_name(kernel_name)
+    return spec is not None and spec.solution == "reference"
+
+
 # ---------------------------------------------------------------------------
 # Dispatch
 # ---------------------------------------------------------------------------
 
 
+def _validate_gemm_out(
+    out: torch.Tensor,
+    *,
+    shape: tuple[int, ...],
+    dtype: torch.dtype,
+    device: torch.device,
+    op: str,
+) -> None:
+    if tuple(out.shape) != shape:
+        raise ValueError(f"{op} out expects shape {shape}, got {tuple(out.shape)}")
+    if out.dtype != dtype:
+        raise ValueError(f"{op} out expects dtype {dtype}, got {out.dtype}")
+    if out.device != device:
+        raise ValueError(f"{op} out expects device {device}, got {out.device}")
+    if out.stride(-1) != 1:
+        raise ValueError(f"{op} out must have stride(-1) == 1")
+
+
 def mm(
     A: torch.Tensor,
     B: torch.Tensor,
@@ -252,6 +277,7 @@ def mm(
     A_scales: torch.Tensor | None = None,
     B_scales: torch.Tensor | None = None,
     bias: torch.Tensor | None = None,
+    out: torch.Tensor | None = None,
     out_dtype: torch.dtype | None = None,
     alpha: torch.Tensor | None = None,
     block_size: list[int] | None = None,
@@ -276,6 +302,8 @@ def mm(
             output.  When the selected kernel supports a fused bias
             epilogue (see ``_KERNELS_WITH_FUSED_BIAS``) it is passed
             into the kernel; otherwise it is added after the GEMM.
+        out: Optional output buffer. The output may be a strided view
+            but must have contiguous rows (``stride(-1) == 1``).
         out_dtype: Output dtype (defaults to ``A.dtype``).
         alpha: Global scaling factor (nvfp4 only).
         block_size: Block size for block-wise quantization, e.g.
@@ -283,11 +311,13 @@ def mm(
         quant: Explicit quant type override.  One of ``"mxfp8"``,
             ``"fp8"``, ``"nvfp4"``, ``"mxfp4"``, ``"none"``.
             If ``None``, inferred from input dtypes and scales.
+        enable_pdl: Whether to request Programmatic Dependent Launch support
+            from kernels that accept it.
         override: Force selection of a specific kernel by name (e.g.
             ``"cublaslt_mm_nvfp4"``). Bypasses heuristic scoring.
         expected_kernel_name: Debug hint for expected kernel selection.
     """
-    out_dtype = out_dtype or A.dtype
+    out_dtype = out_dtype or (out.dtype if out is not None else A.dtype)
 
     M = A.shape[0]
     if quant == "mxfp4":
@@ -297,6 +327,15 @@ def mm(
         K = A.shape[-1]
         N = B.shape[-1] if B.shape[0] == K else B.shape[0]
 
+    if out is not None:
+        _validate_gemm_out(
+            out,
+            shape=(M, N),
+            dtype=out_dtype,
+            device=A.device,
+            op="mm",
+        )
+
     traits: dict[str, object] = {
         "n_align_16": N % 16 == 0,
         "k_align_16": K % 16 == 0,
@@ -321,14 +360,23 @@ def mm(
     )
 
     # Online activation quantization
-    if quant == "mxfp8" and A_scales is None:
+    if (
+        quant == "mxfp8"
+        and A_scales is None
+        and not _kernel_handles_online_mxfp8(kernel.name)
+    ):
         assert (
             block_size is not None
         ), "block_size is required for online activation quantization"
         A, A_scales = _online_quantize_mxfp8(A, block_size, kernel.name)
 
     kernel_args = (A, B, A_scales, B_scales, out_dtype)
-    kernel_kwargs: dict[str, object] = {"alpha": alpha, "block_size": block_size}
+    kernel_kwargs: dict[str, object] = {
+        "alpha": alpha,
+        "block_size": block_size,
+    }
+    if out is not None:
+        kernel_kwargs["out"] = out
 
     fused_bias = bias is not None and kernel.name in _KERNELS_WITH_FUSED_BIAS
     if fused_bias:
@@ -351,9 +399,167 @@ def mm(
         select_dtype,
         kernel_name=kernel.name,
         **shape_params,
+        has_out=out is not None,
+    ):
+        output = kernel(*kernel_args, **kernel_kwargs)
+
+    if bias is not None and not fused_bias:
+        if out is not None:
+            output.add_(bias.to(dtype=output.dtype))
+        else:
+            output = output + bias.to(dtype=output.dtype)
+    return output
+
+
+def bmm(
+    A: torch.Tensor,
+    B: torch.Tensor,
+    *,
+    A_scales: torch.Tensor | None = None,
+    B_scales: torch.Tensor | None = None,
+    bias: torch.Tensor | None = None,
+    out: torch.Tensor | None = None,
+    out_dtype: torch.dtype | None = None,
+    alpha: torch.Tensor | None = None,
+    block_size: list[int] | None = None,
+    quant: str | None = None,
+    enable_pdl: bool = False,
+    override: str | None = None,
+    expected_kernel_name: str | None = None,
+) -> torch.Tensor:
+    """Batched matrix multiply with automatic kernel selection.
+
+    Mirrors :func:`mm` for batched inputs with an outer batch dimension.
+    ``A`` must use ``[B, M, K]`` layout and ``B`` must use ``[B, N, K]``
+    layout. The result shape is ``[B, M, N]``.
+    If ``out`` is provided, kernels that support direct output write into that
+    buffer. The output may be a strided view but must have contiguous rows
+    (``stride(-1) == 1``).
+
+    Args:
+        A: Activation batch ``[batch, M, K]``.
+        B: Weight batch ``[batch, N, K]``.
+        A_scales: Activation scales.
+        B_scales: Weight scales (layout depends on quant type).
+        bias: Optional bias vector of shape ``[N]`` or per-batch bias matrix
+            of shape ``[batch, N]`` added to the output. When the selected
+            kernel supports a fused bias epilogue (see
+            ``_KERNELS_WITH_FUSED_BIAS``) it is passed into the kernel;
+            otherwise it is added after the BMM.
+        out: Optional output buffer. The output may be a strided view
+            but must have contiguous rows (``stride(-1) == 1``).
+        out_dtype: Output dtype (defaults to ``A.dtype``).
+        alpha: Global scaling factor (nvfp4 only).
+        block_size: Block size for block-wise quantization, e.g.
+            ``[128, 128]``.
+        quant: Explicit quant type override. One of ``"mxfp8"``,
+            ``"fp8"``, ``"nvfp4"``, ``"mxfp4"``, ``"none"``.
+            If ``None``, inferred from input dtypes and scales.
+        enable_pdl: Whether to request Programmatic Dependent Launch support
+            from kernels that accept it.
+        override: Force selection of a specific kernel by name. Bypasses
+            heuristic scoring.
+        expected_kernel_name: Debug hint for expected kernel selection.
+    """
+    out_dtype = out_dtype or (out.dtype if out is not None else A.dtype)
+    if A.ndim != 3:
+        raise ValueError(f"bmm expects A with shape [B, M, K], got {tuple(A.shape)}")
+    if B.ndim != 3:
+        raise ValueError(f"bmm expects B with shape [B, N, K], got {tuple(B.shape)}")
+
+    batch, M, A_storage_K = A.shape
+    B_batch, N, B_storage_K = B.shape
+    if B_batch != batch:
+        raise ValueError(f"bmm batch mismatch: A batch={batch}, B batch={B_batch}")
+    if B_storage_K != A_storage_K:
+        raise ValueError(f"bmm K mismatch: A K={A_storage_K}, B K={B_storage_K}")
+    K = A_storage_K * 2 if quant == "mxfp4" else A_storage_K
+
+    if out is not None:
+        _validate_gemm_out(
+            out,
+            shape=(batch, M, N),
+            dtype=out_dtype,
+            device=A.device,
+            op="bmm",
+        )
+
+    traits: dict[str, object] = {
+        "n_align_16": N % 16 == 0,
+        "k_align_16": K % 16 == 0,
+        "n_align_64": N % 64 == 0,
+        "n_align_128": N % 128 == 0,
+        "k_align_64": K % 64 == 0,
+        "k_align_128": K % 128 == 0,
+    }
+
+    signature = _gemm_format_signature(
+        A, B, A_scales, B_scales, out_dtype, quant, block_size
+    )
+    select_dtype = signature.storage_dtype_for("a") or A.dtype
+
+    kernel = select_kernel(
+        "gemm",
+        "bmm",
+        signature,
+        traits=traits,
+        override=override,
+        expected_kernel_name=expected_kernel_name,
+    )
+
+    if (
+        quant == "mxfp8"
+        and A_scales is None
+        and not _kernel_handles_online_mxfp8(kernel.name)
+    ):
+        assert (
+            block_size is not None
+        ), "block_size is required for online activation quantization"
+        A, A_scales = _online_quantize_mxfp8(A, block_size, kernel.name)
+
+    kernel_args = (A, B, A_scales, B_scales, out_dtype)
+    kernel_kwargs: dict[str, object] = {
+        "alpha": alpha,
+        "block_size": block_size,
+    }
+    if out is not None:
+        kernel_kwargs["out"] = out
+
+    fused_bias = bias is not None and kernel.name in _KERNELS_WITH_FUSED_BIAS
+    if fused_bias:
+        kernel_kwargs["bias"] = bias
+
+    if kernel.name in _KERNELS_WITH_PDL:
+        kernel_kwargs["enable_pdl"] = enable_pdl
+
+    shape_params = {"B": batch, "M": M, "N": N, "K": K}
+    ShapeCapture.get().record(
+        "gemm",
+        "bmm",
+        kernel.name,
+        select_dtype,
+        shape_params,
+    )
+    with kernel_scope(
+        "gemm",
+        "bmm",
+        select_dtype,
+        kernel_name=kernel.name,
+        **shape_params,
+        has_out=out is not None,
     ):
         output = kernel(*kernel_args, **kernel_kwargs)
 
     if bias is not None and not fused_bias:
-        output = output + bias.to(dtype=output.dtype)
+        bias = bias.to(dtype=output.dtype)
+        if bias.ndim == 1:
+            bias_view = bias.view(1, 1, -1)
+        elif bias.ndim == 2:
+            bias_view = bias.view(bias.shape[0], 1, bias.shape[1])
+        else:
+            raise ValueError(f"bmm bias expects shape [N] or [B, N], got {bias.shape}")
+        if out is not None:
+            output.add_(bias_view)
+        else:
+            output = output + bias_view
     return output
diff --git a/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/deep_gemm.py b/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/deep_gemm.py
index f484f5743..8f6980609 100644
--- a/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/deep_gemm.py
+++ b/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/deep_gemm.py
@@ -80,6 +80,7 @@ def deep_gemm_mm_fp8_blockscale(
         *,
         alpha: torch.Tensor | None = None,
         block_size: list[int] | None = None,
+        out: torch.Tensor | None = None,
     ) -> torch.Tensor:
         assert (
             A_scales is not None
@@ -89,4 +90,8 @@ def deep_gemm_mm_fp8_blockscale(
         N = B.shape[0]
         C = A.new_empty(A.shape[0], N, dtype=torch.bfloat16)
         fp8_gemm_nt((A, A_scales), (B, B_scales), C)
-        return C.to(out_dtype)
+        output = C.to(out_dtype)
+        if out is not None:
+            out.copy_(output)
+            return out
+        return output
diff --git a/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/flashinfer.py b/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/flashinfer.py
index c078e6883..8a90338f8 100644
--- a/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/flashinfer.py
+++ b/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/flashinfer.py
@@ -114,6 +114,7 @@ def flashinfer_mm_fp8_blockscale(
         *,
         alpha: torch.Tensor | None = None,
         block_size: list[int] | None = None,
+        out: torch.Tensor | None = None,
     ) -> torch.Tensor:
         assert (
             A_scales is not None
@@ -140,7 +141,11 @@ def flashinfer_mm_fp8_blockscale(
             scale_major_mode="MN",
             out_dtype=out_dtype,
         )
-        return output[:orig_m] if output.shape[0] != orig_m else output
+        output = output[:orig_m] if output.shape[0] != orig_m else output
+        if out is not None:
+            out.copy_(output)
+            return out
+        return output
 
 
 # ---- FlashInfer FP4 -----------------------------------------------------
@@ -178,9 +183,10 @@ def flashinfer_mm_nvfp4(
         alpha: torch.Tensor | None = None,
         block_size: list[int] | None = None,
         enable_pdl: bool = False,
+        out: torch.Tensor | None = None,
     ) -> torch.Tensor:
         # backend="cutlass" (not "auto") to skip flashinfer's cuDNN-graph plan compile.
-        return mm_fp4(
+        output = mm_fp4(
             A,
             B,
             A_scales,
@@ -190,3 +196,7 @@ def flashinfer_mm_nvfp4(
             backend="cutlass",
             enable_pdl=enable_pdl,
         )
+        if out is not None:
+            out.copy_(output)
+            return out
+        return output
diff --git a/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/triton.py b/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/triton.py
index ee0db163b..2a181efe4 100644
--- a/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/triton.py
+++ b/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/triton.py
@@ -25,7 +25,7 @@
 import logging
 import math
 import os
-from typing import Any, Dict, List, Optional, Tuple
+from typing import Any, Dict, List, Optional
 
 import torch
 from tokenspeed_kernel._triton import tl, triton
@@ -72,7 +72,8 @@ def prepare_block_fp8_matmul_inputs(
     Bs: torch.Tensor,
     block_size: List[int],
     output_dtype: torch.dtype = torch.float16,
-) -> Tuple[int, int, int]:
+    out: torch.Tensor | None = None,
+) -> tuple[int, int, int, torch.Tensor]:
     assert len(block_size) == 2
     block_n, block_k = block_size[0], block_size[1]
 
@@ -108,7 +109,7 @@ def prepare_block_fp8_matmul_inputs(
         raise NotImplementedError
 
     C_shape = A.shape[:-1] + (N,)
-    C = A.new_empty(C_shape, dtype=output_dtype)
+    C = out if out is not None else A.new_empty(C_shape, dtype=output_dtype)
 
     return M, N, K, C
 
@@ -415,6 +416,7 @@ def w8a8_block_fp8_matmul_triton(
     Bs: torch.Tensor,
     block_size: List[int],
     output_dtype: torch.dtype = torch.float16,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor:
     """This function performs matrix multiplication with block-wise quantization.
 
@@ -432,7 +434,9 @@ def w8a8_block_fp8_matmul_triton(
     Returns:
         torch.Tensor: The result of matmul.
     """
-    M, N, K, C = prepare_block_fp8_matmul_inputs(A, B, As, Bs, block_size, output_dtype)
+    M, N, K, C = prepare_block_fp8_matmul_inputs(
+        A, B, As, Bs, block_size, output_dtype, out=out
+    )
 
     block_n, block_k = block_size
 
@@ -648,6 +652,7 @@ def triton_scaled_mm(
     block_size_n: int = 32,
     block_size_k: int = 32,
     use_heuristic=True,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor:
     M, K = input.shape
     N = weight.shape[1]
@@ -671,7 +676,15 @@ def triton_scaled_mm(
         triton.cdiv(M, META["BLOCK_SIZE_M"]) * triton.cdiv(N, META["BLOCK_SIZE_N"]),
     )
 
-    result = torch.empty((M, N), dtype=out_dtype, device=input.device)
+    result = (
+        out
+        if out is not None
+        else torch.empty((M, N), dtype=out_dtype, device=input.device)
+    )
+    if result.dtype != out_dtype:
+        raise ValueError(
+            f"triton_scaled_mm out expects dtype {out_dtype}, got {result.dtype}"
+        )
 
     has_scalar = lambda x: x.shape[0] == 1 and x.shape[1] == 1
 
@@ -720,7 +733,7 @@ def triton_scaled_mm(
         BLOCK_SIZE_SCALE_B=block_size_sb,
     )
 
-    return result.to(out_dtype)
+    return result
 
 
 # ---- Triton block-scaled FP8 ----------------------------------------------
@@ -748,6 +761,7 @@ def triton_mm_fp8_blockscale(
     *,
     alpha: torch.Tensor | None = None,
     block_size: list[int] | None = None,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor:
     assert block_size is not None, "block_size is required for triton_mm_fp8_blockscale"
     assert (
@@ -760,6 +774,7 @@ def triton_mm_fp8_blockscale(
         B_scales,
         block_size=block_size,
         output_dtype=out_dtype,
+        out=out,
     )
 
 
@@ -855,6 +870,7 @@ def triton_mm_mxfp4(
     *,
     alpha: torch.Tensor | None = None,
     block_size: list[int] | None = None,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor:
     del alpha, block_size
     if A.dtype != torch.uint8 or B.dtype != torch.uint8:
@@ -868,7 +884,7 @@ def triton_mm_mxfp4(
     N = B.shape[0]
     if K % 32 != 0:
         raise ValueError("MXFP4 GEMM requires K divisible by 32")
-    C = torch.empty(M, N, device=A.device, dtype=out_dtype)
+    C = out if out is not None else torch.empty(M, N, device=A.device, dtype=out_dtype)
     grid = (triton.cdiv(M, 16), triton.cdiv(N, 32))
     _mxfp4_mm_kernel[grid](
         A,
@@ -926,6 +942,7 @@ def triton_mm_fp8_scaled(
     alpha: torch.Tensor | None = None,
     block_size: list[int] | None = None,
     bias: torch.Tensor | None = None,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor:
     return triton_scaled_mm(
         A,
@@ -934,4 +951,5 @@ def triton_mm_fp8_scaled(
         B_scales,
         out_dtype=out_dtype,
         bias=bias,
+        out=out,
     )
diff --git a/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/trtllm.py b/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/trtllm.py
index 3bc0faf1a..c819578e6 100644
--- a/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/trtllm.py
+++ b/tokenspeed-kernel/python/tokenspeed_kernel/ops/gemm/trtllm.py
@@ -97,9 +97,10 @@ def cublaslt_mm_nvfp4(
     *,
     alpha: torch.Tensor,
     block_size: list[int] | None = None,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor:
     runner = _get_runner(out_dtype)
-    return runner.run_gemm(
+    output = runner.run_gemm(
         A,
         B.T,
         A_scales,
@@ -108,3 +109,7 @@ def cublaslt_mm_nvfp4(
         False,
         _CUBLASLT_HEURISTIC_ALGO,
     )
+    if out is not None:
+        out.copy_(output)
+        return out
+    return output
diff --git a/tokenspeed-kernel/test/ops/test_gemm.py b/tokenspeed-kernel/test/ops/test_gemm.py
new file mode 100644
index 000000000..66779e4f0
--- /dev/null
+++ b/tokenspeed-kernel/test/ops/test_gemm.py
@@ -0,0 +1,77 @@
+# Copyright (c) 2026 LightSeek Foundation
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+from __future__ import annotations
+
+import pytest
+import tokenspeed_kernel
+import torch
+
+
+def test_mm_rejects_bad_out_layout() -> None:
+    a = torch.empty((4, 8), dtype=torch.bfloat16)
+    b = torch.empty((16, 8), dtype=torch.bfloat16)
+    out = torch.empty((16, 4), dtype=torch.bfloat16).transpose(0, 1)
+
+    with pytest.raises(ValueError, match=r"stride\(-1\) == 1"):
+        tokenspeed_kernel.mm(a, b, out=out)
+
+
+def test_mm_reference_rejects_out_dtype_mismatch() -> None:
+    a = torch.empty((4, 8), dtype=torch.float32)
+    b = torch.empty((16, 8), dtype=torch.float32)
+    out = torch.empty((4, 16), dtype=torch.bfloat16)
+
+    with pytest.raises(ValueError, match="torch_mm out= requires out_dtype"):
+        tokenspeed_kernel.mm(a, b, out=out, override="torch_mm")
+
+
+def test_bmm_rejects_batch_mismatch() -> None:
+    a = torch.empty((2, 4, 8), dtype=torch.bfloat16)
+    b = torch.empty((3, 16, 8), dtype=torch.bfloat16)
+
+    with pytest.raises(ValueError, match="batch mismatch"):
+        tokenspeed_kernel.bmm(a, b)
+
+
+def test_bmm_rejects_rank2_weights() -> None:
+    a = torch.empty((2, 4, 8), dtype=torch.bfloat16)
+    b = torch.empty((16, 8), dtype=torch.bfloat16)
+
+    with pytest.raises(ValueError, match=r"B with shape \[B, N, K\]"):
+        tokenspeed_kernel.bmm(a, b)
+
+
+def test_bmm_rejects_bad_out_layout() -> None:
+    a = torch.empty((2, 4, 8), dtype=torch.bfloat16)
+    b = torch.empty((2, 16, 8), dtype=torch.bfloat16)
+    out = torch.empty((2, 16, 4), dtype=torch.bfloat16).transpose(1, 2)
+
+    with pytest.raises(ValueError, match=r"stride\(-1\) == 1"):
+        tokenspeed_kernel.bmm(a, b, out=out)
+
+
+def test_bmm_reference_rejects_out_dtype_mismatch() -> None:
+    a = torch.empty((2, 4, 8), dtype=torch.float32)
+    b = torch.empty((2, 16, 8), dtype=torch.float32)
+    out = torch.empty((2, 4, 16), dtype=torch.bfloat16)
+
+    with pytest.raises(ValueError, match="torch_bmm out= requires out_dtype"):
+        tokenspeed_kernel.bmm(a, b, out=out, override="torch_bmm")
diff --git a/tokenspeed-kernel/test/test_benchmark.py b/tokenspeed-kernel/test/test_benchmark.py
index e5fcd6a1a..c6a538b5c 100644
--- a/tokenspeed-kernel/test/test_benchmark.py
+++ b/tokenspeed-kernel/test/test_benchmark.py
@@ -53,6 +53,7 @@ def _torch_mm(
     *,
     alpha: torch.Tensor | None = None,
     block_size: list[int] | None = None,
+    out: torch.Tensor | None = None,
 ) -> torch.Tensor:
     _ = A_scales, B_scales, alpha, block_size
     return (A @ B).to(out_dtype)
@@ -232,6 +233,36 @@ def test_export_import_roundtrip(tmp_path, setup_gemm_case):
     assert loaded[0].max_rel_diff == results[0].max_rel_diff
 
 
+def test_bmm_throughput_with_batched_weights():
+    shape = {"batch": 2, "M": 4, "N": 8, "K": 16}
+
+    tflops, bandwidth = ThroughputCalculator.compute(
+        "gemm",
+        "bmm",
+        shape,
+        latency_us=1000.0,
+        dtype=torch.float16,
+    )
+
+    assert tflops == pytest.approx(2.048e-6)
+    assert bandwidth == pytest.approx(0.000896)
+
+
+def test_bmm_throughput_missing_required_shape_returns_none():
+    shape = {"batch": 2, "M": 4, "N": 8}
+
+    tflops, bandwidth = ThroughputCalculator.compute(
+        "gemm",
+        "bmm",
+        shape,
+        latency_us=1000.0,
+        dtype=torch.float16,
+    )
+
+    assert tflops is None
+    assert bandwidth is None
+
+
 def test_attention_decode_throughput_with_explicit_kv_heads():
     shape = {
         "batch": 2,
diff --git a/tokenspeed-kernel/test/test_kernel_api_selection.py b/tokenspeed-kernel/test/test_kernel_api_selection.py
index dd72fc66a..a930a506d 100644
--- a/tokenspeed-kernel/test/test_kernel_api_selection.py
+++ b/tokenspeed-kernel/test/test_kernel_api_selection.py
@@ -274,6 +274,12 @@ def _mm_dense_cdna4_aligned() -> torch.Tensor:
     return tokenspeed_kernel.mm(a, b)
 
 
+def _bmm_dense() -> torch.Tensor:
+    a = torch.empty((4, 2, 16), dtype=torch.bfloat16)
+    b = torch.empty((4, 32, 16), dtype=torch.bfloat16)
+    return tokenspeed_kernel.bmm(a, b)
+
+
 def _mm_mxfp8() -> torch.Tensor:
     a = torch.empty((4, 128), dtype=_fp8_dtype())
     b = torch.empty((128, 128), dtype=_fp8_dtype())
@@ -317,6 +323,33 @@ def test_gemm_mxfp8_online_activation_signature_uses_quantized_storage() -> None
     assert b_format.scale.block_shape == (128, 128)
 
 
+def test_bmm_mxfp8_online_activation_signature_uses_quantized_storage() -> None:
+    a = torch.empty((2, 4, 128), dtype=torch.bfloat16)
+    b = torch.empty((2, 128, 128), dtype=_fp8_dtype())
+    b_scales = torch.empty((2, 1, 1), dtype=torch.float32)
+
+    signature = _gemm_pkg._gemm_format_signature(
+        a,
+        b,
+        None,
+        b_scales,
+        torch.bfloat16,
+        "mxfp8",
+        [128, 128],
+    )
+
+    a_format = signature.format_for("a")
+    b_format = signature.format_for("b")
+    assert a_format is not None
+    assert b_format is not None
+    assert a_format.storage_dtype == _fp8_dtype()
+    assert b_format.storage_dtype == _fp8_dtype()
+    assert a_format.scale is not None
+    assert b_format.scale is not None
+    assert a_format.scale.block_shape == (128, 128)
+    assert b_format.scale.block_shape == (128, 128)
+
+
 def test_gemm_mxfp8_online_activation_preserves_repeated_rows() -> None:
     if not torch.cuda.is_available():
         pytest.skip("CUDA is required for online mxfp8 GEMM verification")
@@ -386,6 +419,32 @@ def test_gemm_fp8_scaled_signature_uses_fp8_format_with_scale() -> None:
         assert tensor_format.scale.storage_dtype == torch.float32
 
 
+def test_bmm_fp8_scaled_signature_uses_fp8_format_with_scale() -> None:
+    a = torch.empty((2, 4, 128), dtype=_fp8_dtype())
+    b = torch.empty((2, 128, 128), dtype=_fp8_dtype())
+    a_scales = torch.empty((1,), dtype=torch.float32)
+    b_scales = torch.empty((1,), dtype=torch.float32)
+
+    signature = _gemm_pkg._gemm_format_signature(
+        a,
+        b,
+        a_scales,
+        b_scales,
+        torch.bfloat16,
+        "fp8",
+        None,
+    )
+
+    for role in ("a", "b"):
+        tensor_format = signature.format_for(role)
+        assert tensor_format is not None
+        assert tensor_format.format == "scaled-fp8"
+        assert tensor_format.storage_dtype == _fp8_dtype()
+        assert tensor_format.scale is not None
+        assert tensor_format.scale.granularity == "tensor"
+        assert tensor_format.scale.storage_dtype == torch.float32
+
+
 def test_gemm_fp8_scaled_signature_uses_channel_granularity() -> None:
     a = torch.empty((4, 128), dtype=_fp8_dtype())
     b = torch.empty((128, 128), dtype=_fp8_dtype())
@@ -409,6 +468,179 @@ def test_gemm_fp8_scaled_signature_uses_channel_granularity() -> None:
         assert tensor_format.scale.granularity == "channel"
 
 
+def test_bmm_fp8_scaled_signature_uses_channel_granularity() -> None:
+    a = torch.empty((4, 2, 128), dtype=_fp8_dtype())
+    b = torch.empty((4, 32, 128), dtype=_fp8_dtype())
+    a_scales = torch.empty((4, 2), dtype=torch.float32)
+    b_scales = torch.empty((4, 32), dtype=torch.float32)
+
+    signature = _gemm_pkg._gemm_format_signature(
+        a,
+        b,
+        a_scales,
+        b_scales,
+        torch.bfloat16,
+        "fp8",
+        None,
+    )
+
+    for role in ("a", "b"):
+        tensor_format = signature.format_for(role)
+        assert tensor_format is not None
+        assert tensor_format.scale is not None
+        assert tensor_format.scale.granularity == "channel"
+
+
+def test_gemm_quantized_reference_dispatches_fp8_inputs() -> None:
+    fp8_dtype = _fp8_dtype()
+    a = torch.zeros((4, 128), dtype=fp8_dtype)
+    a_bf16 = torch.zeros((4, 128), dtype=torch.bfloat16)
+    b = torch.zeros((128, 128), dtype=fp8_dtype)
+    tensor_scales = torch.ones((1,), dtype=torch.float32)
+    block_a_scales = torch.ones((4, 1), dtype=torch.float32)
+    block_b_scales = torch.ones((1, 1), dtype=torch.float32)
+
+    blockscale = tokenspeed_kernel.mm(
+        a,
+        b,
+        A_scales=block_a_scales,
+        B_scales=block_b_scales,
+        out_dtype=torch.bfloat16,
+        block_size=[128, 128],
+        quant="mxfp8",
+        override="torch_mm_fp8_blockscale",
+    )
+    assert blockscale.shape == (4, 128)
+    assert blockscale.dtype == torch.bfloat16
+
+    online_blockscale = tokenspeed_kernel.mm(
+        a_bf16,
+        b,
+        B_scales=block_b_scales,
+        out_dtype=torch.bfloat16,
+        block_size=[128, 128],
+        quant="mxfp8",
+        override="torch_mm_fp8_blockscale",
+    )
+    assert online_blockscale.shape == (4, 128)
+    assert online_blockscale.dtype == torch.bfloat16
+
+    tensor_scaled = tokenspeed_kernel.mm(
+        a,
+        b,
+        A_scales=tensor_scales,
+        B_scales=tensor_scales,
+        out_dtype=torch.bfloat16,
+        quant="fp8",
+        override="torch_mm_fp8_scaled_mnk",
+    )
+    assert tensor_scaled.shape == (4, 128)
+    assert tensor_scaled.dtype == torch.bfloat16
+
+
+def test_bmm_quantized_reference_dispatches_fp8_inputs() -> None:
+    fp8_dtype = _fp8_dtype()
+    a = torch.zeros((2, 4, 128), dtype=fp8_dtype)
+    a_bf16 = torch.zeros((2, 4, 128), dtype=torch.bfloat16)
+    b = torch.zeros((2, 128, 128), dtype=fp8_dtype)
+    tensor_scales = torch.ones((1,), dtype=torch.float32)
+    channel_a_scales = torch.ones((2, 4), dtype=torch.float32)
+    channel_b_scales = torch.ones((2, 128), dtype=torch.float32)
+    block_a_scales = torch.ones((2, 4, 1), dtype=torch.float32)
+    block_b_scales = torch.ones((2, 1, 1), dtype=torch.float32)
+
+    blockscale = tokenspeed_kernel.bmm(
+        a,
+        b,
+        A_scales=block_a_scales,
+        B_scales=block_b_scales,
+        out_dtype=torch.bfloat16,
+        block_size=[128, 128],
+        quant="mxfp8",
+        override="torch_bmm_fp8_blockscale",
+    )
+    assert blockscale.shape == (2, 4, 128)
+    assert blockscale.dtype == torch.bfloat16
+
+    online_blockscale = tokenspeed_kernel.bmm(
+        a_bf16,
+        b,
+        B_scales=block_b_scales,
+        out_dtype=torch.bfloat16,
+        block_size=[128, 128],
+        quant="mxfp8",
+        override="torch_bmm_fp8_blockscale",
+    )
+    assert online_blockscale.shape == (2, 4, 128)
+    assert online_blockscale.dtype == torch.bfloat16
+
+    tensor_scaled = tokenspeed_kernel.bmm(
+        a,
+        b,
+        A_scales=tensor_scales,
+        B_scales=tensor_scales,
+        out_dtype=torch.bfloat16,
+        quant="fp8",
+        override="torch_bmm_fp8_scaled",
+    )
+    assert tensor_scaled.shape == (2, 4, 128)
+    assert tensor_scaled.dtype == torch.bfloat16
+
+    channel_scaled = tokenspeed_kernel.bmm(
+        a,
+        b,
+        A_scales=channel_a_scales,
+        B_scales=channel_b_scales,
+        out_dtype=torch.bfloat16,
+        quant="fp8",
+        override="torch_bmm_fp8_scaled",
+    )
+    assert channel_scaled.shape == (2, 4, 128)
+    assert channel_scaled.dtype == torch.bfloat16
+
+
+def _copy_out_mm_kernel(
+    A: torch.Tensor,
+    B: torch.Tensor,
+    A_scales: torch.Tensor | None,
+    B_scales: torch.Tensor | None,
+    out_dtype: torch.dtype,
+    *,
+    alpha: torch.Tensor | None = None,
+    block_size: list[int] | None = None,
+    out: torch.Tensor | None = None,
+) -> torch.Tensor:
+    assert A_scales is None
+    assert B_scales is None
+    assert block_size is None
+    output = A @ B.T
+    if alpha is not None:
+        output = output * alpha.to(dtype=output.dtype)
+    output = output.to(out_dtype)
+    if out is not None:
+        out.copy_(output)
+        return out
+    return output
+
+
+def test_mm_non_native_out_kernel_copies_to_out(monkeypatch) -> None:
+    torch.manual_seed(1)
+    a = torch.randn((4, 8), dtype=torch.float32)
+    b = torch.randn((16, 8), dtype=torch.float32)
+    out = torch.empty((4, 16), dtype=torch.float32)
+
+    def select_copy_out_kernel(*args, **kwargs) -> SelectedKernel:
+        return SelectedKernel("test_mm_copy_out_kernel", _copy_out_mm_kernel)
+
+    monkeypatch.setattr(_gemm_pkg, "select_kernel", select_copy_out_kernel)
+
+    actual = tokenspeed_kernel.mm(a, b, out=out, override="test_mm_copy_out_kernel")
+    expected = a @ b.T
+
+    assert actual is out
+    torch.testing.assert_close(out, expected)
+
+
 def _mm_nvfp4() -> torch.Tensor:
     a = torch.empty((4, 64), dtype=torch.uint8)
     b = torch.empty((128, 64), dtype=torch.uint8)
@@ -449,6 +681,29 @@ def test_gemm_nvfp4_signature_uses_fixed_block_shape() -> None:
         assert tensor_format.scale.block_shape == (16,)
 
 
+def test_bmm_nvfp4_signature_uses_fixed_block_shape() -> None:
+    a = torch.empty((2, 4, 64), dtype=torch.uint8)
+    b = torch.empty((2, 128, 64), dtype=torch.uint8)
+    a_scales = torch.empty((2, 4, 1), dtype=torch.float32)
+    b_scales = torch.empty((2, 128, 1), dtype=torch.float32)
+
+    signature = _gemm_pkg._gemm_format_signature(
+        a,
+        b,
+        a_scales,
+        b_scales,
+        torch.bfloat16,
+        "nvfp4",
+        None,
+    )
+
+    for role in ("a", "b"):
+        tensor_format = signature.format_for(role)
+        assert tensor_format is not None
+        assert tensor_format.scale is not None
+        assert tensor_format.scale.block_shape == (16,)
+
+
 def _attention_prefill() -> object:
     q = torch.empty((4, 16, 64), dtype=torch.bfloat16)
     k = torch.empty((4, 8, 64), dtype=torch.bfloat16)
@@ -1687,6 +1942,14 @@ def _case(
     ),
     # GEMM API x architecture golden cases.
     _case(_is_supported_gpu, "supported-gpu", "gemm", "mm", "torch_mm", _mm_dense),
+    _case(
+        _is_supported_gpu,
+        "supported-gpu",
+        "gemm",
+        "bmm",
+        "torch_bmm",
+        _bmm_dense,
+    ),
     _case(
         _is_cdna4,
         "cdna4",
@@ -1921,6 +2184,11 @@ def fake_call(self: SelectedKernel, *args, **kwargs):
 
         if case.family == "gemm":
             a, b, _a_scales, _b_scales, out_dtype = args[:5]
+            if case.mode == "bmm":
+                n = b.shape[-2]
+                return torch.empty(
+                    (a.shape[0], a.shape[1], n), dtype=out_dtype, device=a.device
+                )
             n = b.shape[-1] if b.shape[0] == a.shape[-1] else b.shape[0]
             return torch.empty((a.shape[0], n), dtype=out_dtype, device=a.device)
 
diff --git a/tokenspeed-kernel/test/test_numerics.py b/tokenspeed-kernel/test/test_numerics.py
index f1398d4cf..a237197ec 100644
--- a/tokenspeed-kernel/test/test_numerics.py
+++ b/tokenspeed-kernel/test/test_numerics.py
@@ -24,6 +24,13 @@
 import torch
 from tokenspeed_kernel.numerics.comparison import compare_outputs, format_comparison
 from tokenspeed_kernel.numerics.inputs import get_input_generator
+from tokenspeed_kernel.numerics.reference.gemm import (
+    torch_bmm_fp8_blockscale,
+    torch_bmm_fp8_scaled,
+    torch_mm_fp8_blockscale,
+    torch_mm_fp8_scaled_mnk,
+    torch_mm_fp8_scaled_nkm,
+)
 from tokenspeed_kernel.numerics.tolerance import Tolerance
 from tokenspeed_kernel.numerics.verify import (
     _verification_signature_and_reference,
@@ -115,6 +122,146 @@ def test_gemm_input_generator_requires_mxfp8_block_shape() -> None:
         generator.generate(M=4, N=256, K=128)
 
 
+def test_quantized_reference_gemm_supports_out() -> None:
+    A = torch.empty((2, 128), dtype=_fp8_dtype)
+    B = torch.empty((256, 128), dtype=_fp8_dtype)
+    tensor_scale = torch.ones((1,), dtype=torch.float32)
+    block_a_scale = torch.ones((2, 1), dtype=torch.float32)
+    block_b_scale = torch.ones((2, 1), dtype=torch.float32)
+    out = torch.empty((2, 256), dtype=torch.bfloat16)
+
+    blockscale = torch_mm_fp8_blockscale(
+        A,
+        B,
+        block_a_scale,
+        block_b_scale,
+        torch.bfloat16,
+        block_size=[128, 128],
+        out=out,
+    )
+    assert blockscale is out
+
+    scaled_mnk = torch_mm_fp8_scaled_mnk(
+        A,
+        B,
+        tensor_scale,
+        tensor_scale,
+        torch.bfloat16,
+        out=out,
+    )
+    assert scaled_mnk is out
+
+    scaled_nkm = torch_mm_fp8_scaled_nkm(
+        A,
+        B.T,
+        tensor_scale,
+        tensor_scale,
+        torch.bfloat16,
+        out=out,
+    )
+    assert scaled_nkm is out
+
+
+def test_quantized_reference_bmm_supports_out() -> None:
+    A = torch.empty((2, 4, 128), dtype=_fp8_dtype)
+    B = torch.empty((2, 256, 128), dtype=_fp8_dtype)
+    tensor_scale = torch.ones((1,), dtype=torch.float32)
+    channel_a_scale = torch.ones((2, 4), dtype=torch.float32)
+    channel_b_scale = torch.ones((2, 256), dtype=torch.float32)
+    block_a_scale = torch.ones((2, 4, 1), dtype=torch.float32)
+    block_b_scale = torch.ones((2, 2, 1), dtype=torch.float32)
+    out = torch.empty((2, 4, 256), dtype=torch.bfloat16)
+
+    blockscale = torch_bmm_fp8_blockscale(
+        A,
+        B,
+        block_a_scale,
+        block_b_scale,
+        torch.bfloat16,
+        block_size=[128, 128],
+        out=out,
+    )
+    assert blockscale is out
+
+    tensor_scaled = torch_bmm_fp8_scaled(
+        A,
+        B,
+        tensor_scale,
+        tensor_scale,
+        torch.bfloat16,
+        out=out,
+    )
+    assert tensor_scaled is out
+
+    channel_scaled = torch_bmm_fp8_scaled(
+        A,
+        B,
+        channel_a_scale,
+        channel_b_scale,
+        torch.bfloat16,
+        out=out,
+    )
+    assert channel_scaled is out
+
+
+def test_bmm_input_generator_uses_signature_scale_metadata() -> None:
+    scale = ScaleFormat(
+        storage_dtype=torch.float32,
+        granularity="block",
+        block_shape=(128, 128),
+    )
+    signature = next(
+        iter(format_signatures(("a", "b"), "mxfp8", {_fp8_dtype}, scale=scale))
+    )
+    generator = get_input_generator(
+        "gemm",
+        "bmm",
+        dtype=_fp8_dtype,
+        traits={},
+        format_signature=signature,
+        device="cpu",
+    )
+
+    inputs = generator.generate(batch=2, M=4, N=256, K=128)
+
+    assert inputs["A"].shape == (2, 4, 128)
+    assert inputs["B"].shape == (2, 256, 128)
+    assert inputs["A"].dtype == _fp8_dtype
+    assert inputs["B"].dtype == _fp8_dtype
+    assert inputs["A_scales"].shape == (2, 4, 1)
+    assert inputs["B_scales"].shape == (2, 2, 1)
+    assert inputs["A_scales"].dtype == torch.float32
+    assert inputs["B_scales"].dtype == torch.float32
+    assert inputs["block_size"] == [128, 128]
+
+
+def test_bmm_input_generator_supports_scaled_fp8_channel_scales() -> None:
+    scale = ScaleFormat(storage_dtype=torch.float32, granularity="channel")
+    signature = next(
+        iter(format_signatures(("a", "b"), "scaled-fp8", {_fp8_dtype}, scale=scale))
+    )
+    generator = get_input_generator(
+        "gemm",
+        "bmm",
+        dtype=_fp8_dtype,
+        traits={},
+        format_signature=signature,
+        device="cpu",
+    )
+
+    inputs = generator.generate(batch=2, M=4, N=256, K=128)
+
+    assert inputs["A"].shape == (2, 4, 128)
+    assert inputs["B"].shape == (2, 256, 128)
+    assert inputs["A"].dtype == _fp8_dtype
+    assert inputs["B"].dtype == _fp8_dtype
+    assert inputs["A_scales"].shape == (2, 4)
+    assert inputs["B_scales"].shape == (2, 256)
+    assert inputs["A_scales"].dtype == torch.float32
+    assert inputs["B_scales"].dtype == torch.float32
+    assert inputs["block_size"] is None
+
+
 def test_verification_uses_signature_with_compatible_reference(fresh_registry) -> None:
     tensor_scale = ScaleFormat(storage_dtype=torch.float32, granularity="tensor")
     channel_scale = ScaleFormat(storage_dtype=torch.float32, granularity="channel")
diff --git a/tokenspeed-kernel/test/test_selection.py b/tokenspeed-kernel/test/test_selection.py
index 30044a175..67bf0bb56 100644
--- a/tokenspeed-kernel/test/test_selection.py
+++ b/tokenspeed-kernel/test/test_selection.py
@@ -1354,6 +1354,7 @@ def _impl(
             *,
             alpha: torch.Tensor | None = None,
             block_size: list[int] | None = None,
+            out: torch.Tensor | None = None,
         ) -> torch.Tensor:
             _ = A_scales, B_scales, alpha, block_size
             call_log.append(name)
@@ -1428,6 +1429,7 @@ def test_mm_wraps_triton_kernel_execution_in_scope(self, monkeypatch):
                     "M": 4,
                     "N": 6,
                     "K": 8,
+                    "has_out": False,
                 },
             )
         ]
@@ -1465,6 +1467,7 @@ def test_mm_wraps_non_triton_kernel_execution_in_scope(self, monkeypatch):
                     "M": 4,
                     "N": 6,
                     "K": 8,
+                    "has_out": False,
                 },
             )
         ]