chore: bump version to 1.3.9

- Replace 4*NC8 per-element -FLT_MAX comparisons with 2 row-level pn guards; masked entries naturally underflow via expf(-FLT_MAX - nm) ≈ 0; pn only guards all-masked-row edge where nm == -FLT_MAX (exp(0)=1 not 0); ~1-3% speedup on prefill (verified via standalone CUDA bench); correctness verified: prefill 4/4, decode 3/3, paged decode 13/13
- Fix stale comments: remove false 'pre-scale Q' claim, correct occupancy numbers, remove phantom sQ from smem description

.github/workflows/release.yml

@@ -0,0 +1,71 @@
+name: Release
+
+on:
+  push:
+    tags:
+      - "v*"
+
+jobs:
+  build-pure:
+    name: Build pure-Python wheel
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.12"
+
+      - name: Build wheel (no CUDA)
+        run: |
+          pip wheel . --no-deps -w dist/
+
+      - uses: actions/upload-artifact@v4
+        with:
+          name: pure-wheel
+          path: dist/*.whl
+
+  build-cuda-linux:
+    name: Build CUDA wheel (Linux)
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.12"
+
+      - name: Install torch (CUDA 12.8)
+        run: |
+          pip install torch --index-url https://download.pytorch.org/whl/cu128
+
+      - name: Setup CUDA
+        uses: Jimver/cuda-toolkit@v0.2.35
+        with:
+          cuda: "12.8.0"
+
+      - name: Build wheel (with CUDA kernels)
+        run: |
+          CSRC_KERNELS=true pip wheel . --no-deps --no-build-isolation -w dist/
+
+      - uses: actions/upload-artifact@v4
+        with:
+          name: cuda-wheel-linux
+          path: dist/*.whl
+
+  release:
+    name: Attach wheels to release
+    needs: [build-pure, build-cuda-linux]
+    runs-on: ubuntu-latest
+    permissions:
+      contents: write
+    steps:
+      - uses: actions/download-artifact@v4
+        with:
+          pattern: "*-wheel"
+          merge-multiple: true
+
+      - name: Create release & upload assets
+        uses: softprops/action-gh-release@v2
+        with:
+          files: ./*.whl
+          tag_name: ${{ github.ref_name }}
+          generate_release_notes: true

.gitignore

@@ -11,6 +11,7 @@
 !csrc/**/*.py
 
 !csrc/**/*.cu
+!csrc/**/*.h
 !csrc/**/*.cuh
 
 !scripts/**/*.sh

astrai/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "1.3.8"
+__version__ = "1.3.9"
 __author__ = "ViperEkura"
 
 from astrai.config import (

astrai/extension/__init__.py

@@ -1,19 +1,21 @@
 """CUDA attention kernel wrappers with torch fallback.
 
 Public API:
-    - ``gqa_decode_attn`` — single-query decode attention
-    - ``gqa_prefill_attn`` — multi-query prefill attention
+    - ``attn_decode`` — single-query decode attention
+    - ``attn_prefill`` — multi-query prefill attention
+    - ``attn_paged_decode`` — paged decode attention (direct page-table access)
 
-Each wrapper dispatches to its compiled CUDA kernel (``astrai.extension.gqa_*``)
+Each wrapper dispatches to its compiled CUDA kernel (``astrai.extension.attn_*``)
 when available, otherwise falls back to ``torch.nn.functional.scaled_dot_product_attention``.
 """
 
 from astrai.extension.loader import KERNEL_NAMES, is_available
-from astrai.extension.ops import gqa_decode_attn, gqa_prefill_attn
+from astrai.extension.ops import attn_decode, attn_paged_decode, attn_prefill
 
 __all__ = [
-    "gqa_decode_attn",
-    "gqa_prefill_attn",
+    "attn_decode",
+    "attn_paged_decode",
+    "attn_prefill",
     "is_available",
     "KERNEL_NAMES",
 ]

astrai/extension/loader.py

@@ -11,7 +11,7 @@ import logging
 
 logger = logging.getLogger(__name__)
 
-KERNEL_NAMES = ["gqa_decode_attn", "gqa_prefill_attn"]
+KERNEL_NAMES = ["attn_decode", "attn_prefill", "attn_paged_decode"]
 
 _available: dict[str, bool] = {}
 _modules: dict[str, object] = {}

astrai/extension/ops.py

@@ -42,7 +42,41 @@ def _torch_fallback(
     )
 
 
-def gqa_decode_attn(
+def _gather_kv_from_pages(
+    page_table: torch.Tensor,
+    k_cache: torch.Tensor,
+    v_cache: torch.Tensor,
+    page_size: int,
+    kv_len: int,
+) -> tuple[torch.Tensor, torch.Tensor]:
+    """Gather contiguous K/V from paged cache for torch SDPA fallback.
+
+    Shapes:
+        page_table : [batch, max_pages] (int64)
+        k_cache    : [n_pages, page_size, n_kv_heads, head_dim]
+        v_cache    : same as k_cache
+    Returns:
+        k, v : [batch, n_kv_heads, kv_len, head_dim]
+    """
+    batch, max_pages = page_table.shape
+    n_pages, ps, n_kv_heads, head_dim = k_cache.shape
+    if ps != page_size:
+        raise ValueError(f"k_cache page_size mismatch: {ps} vs {page_size}")
+
+    k = k_cache.new_empty(batch, n_kv_heads, kv_len, head_dim)
+    v = v_cache.new_empty(batch, n_kv_heads, kv_len, head_dim)
+
+    for b in range(batch):
+        for pos in range(kv_len):
+            log_pg = pos // page_size
+            pg_off = pos % page_size
+            phys = int(page_table[b, log_pg].item())
+            k[b, :, pos, :] = k_cache[phys, pg_off, :, :]
+            v[b, :, pos, :] = v_cache[phys, pg_off, :, :]
+    return k, v
+
+
+def attn_decode(
     q: torch.Tensor,
     k: torch.Tensor,
     v: torch.Tensor,
@@ -51,8 +85,8 @@ def gqa_decode_attn(
     causal_offset: int = 0,
     scale: float | None = None,
 ) -> torch.Tensor:
-    if _available["gqa_decode_attn"]:
-        return _modules["gqa_decode_attn"].gqa_decode_attn(
+    if _available["attn_decode"]:
+        return _modules["attn_decode"].attn_decode(
             q,
             k,
             v,
@@ -64,7 +98,7 @@ def gqa_decode_attn(
     return _torch_fallback(q, k, v, mask, is_causal, scale)
 
 
-def gqa_prefill_attn(
+def attn_prefill(
     q: torch.Tensor,
     k: torch.Tensor,
     v: torch.Tensor,
@@ -73,8 +107,8 @@ def gqa_prefill_attn(
     causal_offset: int = 0,
     scale: float | None = None,
 ) -> torch.Tensor:
-    if _available["gqa_prefill_attn"]:
-        return _modules["gqa_prefill_attn"].gqa_prefill_attn(
+    if _available["attn_prefill"]:
+        return _modules["attn_prefill"].attn_prefill(
             q,
             k,
             v,
@@ -84,3 +118,32 @@ def gqa_prefill_attn(
             scale=scale,
         )
     return _torch_fallback(q, k, v, mask, is_causal, scale)
+
+
+def attn_paged_decode(
+    q: torch.Tensor,
+    page_table: torch.Tensor,
+    k_cache: torch.Tensor,
+    v_cache: torch.Tensor,
+    page_size: int,
+    kv_len: int,
+    mask: torch.Tensor | None = None,
+    is_causal: bool = False,
+    causal_offset: int = 0,
+    scale: float | None = None,
+) -> torch.Tensor:
+    if _available["attn_paged_decode"]:
+        return _modules["attn_paged_decode"].attn_paged_decode(
+            q,
+            page_table,
+            k_cache,
+            v_cache,
+            page_size,
+            kv_len,
+            mask=mask,
+            is_causal=is_causal,
+            causal_offset=causal_offset,
+            scale=scale,
+        )
+    k, v = _gather_kv_from_pages(page_table, k_cache, v_cache, page_size, kv_len)
+    return _torch_fallback(q, k, v, mask, is_causal, scale)

astrai/parallel/setup.py

@@ -1,14 +1,21 @@
 import os
+import socket
 from abc import ABC, abstractmethod
 from contextlib import contextmanager
 from functools import wraps
-from typing import Callable
+from typing import Callable, Optional
 
 import torch
 import torch.distributed as dist
 import torch.multiprocessing as mp
 
 
+def find_free_port() -> str:
+    with socket.socket(socket.AF_INET, socket.SOCK_STREAM) as s:
+        s.bind(("", 0))
+        return str(s.getsockname()[1])
+
+
 def get_current_device():
     return os.environ["LOCAL_DEVICE"]
 
@@ -217,11 +224,13 @@ def spawn_parallel_fn(
     world_size: int,
     backend: str = "nccl",
     master_addr: str = "localhost",
-    master_port: str = "29500",
+    master_port: Optional[str] = None,
     device_type: str = "cuda",
     start_method: str = "spawn",
     **kwargs,
 ):
+    if master_port is None:
+        master_port = find_free_port()
     launcher = _detect_launcher()
     if launcher in ("torchelastic", "torchrun", "external"):
         strategy = TorchrunStrategy(

csrc/build.py

@@ -20,13 +20,14 @@ def _arch_flags() -> list[str]:
 _kernels_dir = Path("csrc/kernels")
 REGISTRY: dict[str, dict] = {}
 
-CXX_FLAGS = ["-O3", "-march=native", "-funroll-loops"]
+CXX_FLAGS = ["-O3", "-funroll-loops"]
 NVCC_FLAGS = [
     "-O3",
     "--expt-relaxed-constexpr",
     "--use_fast_math",
     "--ptxas-options=-O3,-v",
     "--extra-device-vectorization",
+    "--threads=8",
 ]
 
 
@@ -42,5 +43,6 @@ def register(name: str, sources: list[str] | None = None, **kwargs):
     }
 
 
-register("gqa_decode_attn")
-register("gqa_prefill_attn")
+register("attn_decode")
+register("attn_prefill")
+register("attn_paged_decode")

csrc/kernels/attn_common.h

@@ -0,0 +1,54 @@
+#pragma once
+
+
+template<typename T, typename AT = float>
+struct AttentionParams {
+    int batch;
+    int q_head;
+    int kv_head;
+    int q_len;
+    int kv_len;
+    int head_dim;
+    int use_mask;
+    int is_causal;
+    int causal_offset;
+    int num_splits;
+    float scale;
+
+    const T* __restrict__ q;
+    const T* __restrict__ k;
+    const T* __restrict__ v;
+    const bool* __restrict__ mask;
+
+    T* __restrict__ o;
+    AT* __restrict__ o_part;
+    AT* __restrict__ ml_part;
+};
+
+template<typename T, typename AT = float>
+struct PagedAttentionParams {
+    int batch;
+    int q_head;
+    int kv_head;
+    int q_len;
+    int kv_len;
+    int head_dim;
+    int use_mask;
+    int is_causal;
+    int causal_offset;
+    float scale;
+
+    int num_splits;
+    int page_size;
+    int max_pages;
+
+    const T* __restrict__ q;
+    const T* __restrict__ k_cache;
+    const T* __restrict__ v_cache;
+    const bool* __restrict__ mask;
+    const int64_t* __restrict__ page_table;
+    
+    T* __restrict__ o;
+    AT* __restrict__ o_part;
+    AT* __restrict__ ml_part;
+};

csrc/kernels/attn_decode.cu

@@ -0,0 +1,113 @@
+#include "attn_decode_split_kv.cuh"
+#include "attn_entry_utils.cuh"
+
+#ifndef ASTRAI_NO_MMA
+#include "attn_decode_split_kv_mma.cuh"
+#endif
+
+static int decode_num_splits(int base_blocks, int tiles_total) {
+    int sm_count = 0;
+    cudaDeviceGetAttribute(&sm_count, cudaDevAttrMultiProcessorCount, 0);
+    int n = (2 * sm_count + base_blocks - 1) / base_blocks;
+    return std::max(1, std::min(n, std::min(tiles_total, 32)));
+}
+
+// Scalar fallback: one warp per query head, split-KV across grid.z.
+static void launch_scalar_decode(AttentionParams<bf16>& p) {
+    int group_size = p.q_head / p.kv_head;
+    int chunks_total = (p.kv_len + DC_CHUNK - 1) / DC_CHUNK;
+    p.num_splits = decode_num_splits(p.batch * p.kv_head, chunks_total);
+
+    auto fopt = torch::TensorOptions().dtype(torch::kFloat32).device(torch::kCUDA);
+    auto o_part = torch::empty({p.batch, p.q_head, p.num_splits, p.head_dim}, fopt);
+    auto ml_part = torch::empty({p.batch, p.q_head, p.num_splits, 2}, fopt);
+    p.o_part = o_part.data_ptr<float>();
+    p.ml_part = ml_part.data_ptr<float>();
+
+    size_t smem = DC_CHUNK * p.head_dim * sizeof(bf16);
+    attn_decode_split_kv_kernel<<<dim3(p.batch * p.kv_head, 1, p.num_splits), dim3(32, group_size), smem>>>(p);
+    attn_decode_combine_kernel<<<p.batch * p.q_head, p.head_dim>>>(p);
+}
+
+#ifndef ASTRAI_NO_MMA
+// MMA head-packing requires G <= 16 (BR=16 rows). sm_80+ tensor-core
+// + cp.async wins even at G=1 (decode is memory-bound, not compute-bound).
+// STAGES=2 (double-buffer) for D<=128 (smem 16 KB); STAGES=1 for D=256
+// (double-buffer would be 32 KB, near the 48 KB static cap — keep single
+// to preserve occupancy).
+template <int HEAD_DIM, int BC, int STAGES = (HEAD_DIM <= 128) ? 2 : 1>
+static void launch_mma_decode(AttentionParams<bf16>& p) {
+    int tiles_total = (p.kv_len + BC - 1) / BC;
+    p.num_splits = decode_num_splits(p.batch * p.kv_head, tiles_total);
+
+    auto fopt = torch::TensorOptions().dtype(torch::kFloat32).device(torch::kCUDA);
+    auto o_part = torch::empty({p.batch, p.q_head, p.num_splits, p.head_dim}, fopt);
+    auto ml_part = torch::empty({p.batch, p.q_head, p.num_splits, 2}, fopt);
+    p.o_part = o_part.data_ptr<float>();
+    p.ml_part = ml_part.data_ptr<float>();
+
+    attn_decode_split_kv_mma_kernel<HEAD_DIM, BC, STAGES><<<dim3(p.kv_head, p.batch, p.num_splits), 32>>>(p);
+    attn_decode_combine_kernel<<<p.batch * p.q_head, p.head_dim>>>(p);
+}
+#endif
+
+template <int HEAD_DIM>
+static void dispatch_decode(AttentionParams<bf16>& p) {
+#ifndef ASTRAI_NO_MMA
+    int G = p.q_head / p.kv_head;
+    if (!p.use_mask && G >= 1 && G <= 16) {
+        launch_mma_decode<HEAD_DIM, 32>(p);
+        return;
+    }
+#endif
+    launch_scalar_decode(p);
+}
+
+torch::Tensor attn_decode(
+    torch::Tensor q,
+    torch::Tensor k,
+    torch::Tensor v,
+    c10::optional<torch::Tensor> mask,
+    bool is_causal = false,
+    int64_t causal_offset = 0,
+    c10::optional<double> scale = c10::nullopt
+) {
+    AttentionParams<bf16> p;
+    attn_pack_params(q, k, v, mask, is_causal, causal_offset, scale, p);
+    TORCH_CHECK(p.q_len == 1, "Q seq_len must be 1");
+    TORCH_CHECK(p.head_dim % 32 == 0, "head_dim must be multiple of 32");
+
+    auto O = torch::empty_like(q);
+    p.o = (bf16*)O.data_ptr();
+
+    switch (p.head_dim) {
+        case 32:
+            dispatch_decode<32>(p);
+            break;
+        case 64:
+            dispatch_decode<64>(p);
+            break;
+        case 128:
+            dispatch_decode<128>(p);
+            break;
+        case 256:
+            dispatch_decode<256>(p);
+            break;
+        default:
+            TORCH_CHECK(false, "decode: unsupported head_dim ", p.head_dim,
+                        " (supported: 32, 64, 128, 256)");
+    }
+    return O;
+}
+
+PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
+    m.def("attn_decode", &attn_decode,
+        py::arg("q"),
+        py::arg("k"),
+        py::arg("v"),
+        py::arg("mask") = py::none(),
+        py::arg("is_causal") = false,
+        py::arg("causal_offset") = 0,
+        py::arg("scale") = py::none(),
+        "GQA decode (tensor-core head-packing on sm_80+, scalar fallback)");
+}

csrc/kernels/attn_decode_split_kv.cuh

@@ -0,0 +1,116 @@
+#pragma once
+#include <cuda_bf16.h>
+#include <float.h>
+#include "attn_common.h"
+
+using bf16 = __nv_bfloat16;
+constexpr int DC_CHUNK = 64;
+
+__device__ inline float warp_reduce_sum(float val) {
+    for (int offset = 16; offset > 0; offset >>= 1)
+        val += __shfl_xor_sync(0xFFFFFFFF, val, offset);
+    return val;
+}
+
+__global__ void attn_decode_split_kv_kernel(AttentionParams<bf16> p) {
+    int batch = blockIdx.x / p.kv_head;
+    int kv_head = blockIdx.x % p.kv_head;
+    int split = blockIdx.z;
+    int group_size = blockDim.y;
+    int q_head = kv_head * group_size + threadIdx.y;
+    int lane = threadIdx.x;
+    int hd_per_thread = p.head_dim / 32;
+
+    float q_reg[8];
+    int q_off = ((batch * p.q_head + q_head) * 1) * p.head_dim + lane * hd_per_thread;
+    for (int i = 0; i < hd_per_thread; i++)
+        q_reg[i] = __bfloat162float(p.q[q_off + i]);
+
+    int kv_base = ((batch * p.kv_head + kv_head) * p.kv_len) * p.head_dim;
+    int mask_base = batch * p.kv_len;
+
+    float m = -FLT_MAX, d = 0.0f, acc_reg[8] = {0.0f};
+
+    extern __shared__ __align__(16) bf16 k_smem[];
+
+    // Split-KV: each split processes a contiguous subset of chunks
+    int chunks_total = (p.kv_len + DC_CHUNK - 1) / DC_CHUNK;
+    int chunks_per_split = (chunks_total + p.num_splits - 1) / p.num_splits;
+    int ch_begin = split * chunks_per_split;
+    int ch_end = min(chunks_total, ch_begin + chunks_per_split);
+
+    for (int ci = ch_begin; ci < ch_end; ci++) {
+        int chunk_start = ci * DC_CHUNK;
+        int this_chunk = min(DC_CHUNK, p.kv_len - chunk_start);
+
+        int total = this_chunk * p.head_dim;
+        for (int i = threadIdx.y * 32 + lane; i < total; i += blockDim.x * blockDim.y)
+            k_smem[i] = p.k[kv_base + chunk_start * p.head_dim + i];
+        __syncthreads();
+
+        for (int s = 0; s < this_chunk; s++) {
+            float partial = 0.0f;
+            for (int i = 0; i < hd_per_thread; i++)
+                partial += q_reg[i] * __bfloat162float(k_smem[s * p.head_dim + lane * hd_per_thread + i]);
+            partial = warp_reduce_sum(partial) * p.scale;
+
+            if (p.use_mask && p.mask && !p.mask[mask_base + chunk_start + s])
+                partial = -FLT_MAX;
+            if (p.is_causal && (chunk_start + s) > p.causal_offset)
+                partial = -FLT_MAX;
+
+            float new_m = fmaxf(m, partial);
+            float alpha = expf(m - new_m);
+            float beta  = expf(partial - new_m);
+            d = d * alpha + beta;
+
+            int v_off = kv_base + (chunk_start + s) * p.head_dim + lane * hd_per_thread;
+            for (int i = 0; i < hd_per_thread; i++)
+                acc_reg[i] = acc_reg[i] * alpha + __bfloat162float(p.v[v_off + i]) * beta;
+            m = new_m;
+        }
+        __syncthreads();
+    }
+
+    // ---- write UN-normalised partials for this split ----
+    size_t bh = (size_t)batch * p.q_head + q_head;
+    size_t slot = bh * p.num_splits + split;
+    int d0 = lane * hd_per_thread;
+    for (int i = 0; i < hd_per_thread; i++) {
+        int dd = d0 + i;
+        p.o_part[slot * p.head_dim + dd] = acc_reg[i];
+    }
+    if (lane == 0) {
+        p.ml_part[slot * 2] = m;
+        p.ml_part[slot * 2 + 1] = d;
+    }
+}
+
+// Reduce split-K partials into the final bf16 output. One block per (batch,
+// q_head); each thread folds across all splits with a single-pass
+// online-rescale reduction (expf + FMA counts halved vs 3-pass original).
+__global__ void attn_decode_combine_kernel(AttentionParams<bf16> p) {
+    int bh = blockIdx.x;
+    int d = threadIdx.x;
+    if (d >= p.head_dim) return;
+
+    size_t split_base = (size_t)bh * p.num_splits;
+    const float* mlp = p.ml_part + split_base * 2;
+    const float* op = p.o_part + split_base * p.head_dim;
+
+    float m = -FLT_MAX, l = 0.0f, acc = 0.0f;
+    for (int s = 0; s < p.num_splits; s++) {
+        float mi = mlp[s * 2];
+        if (mi <= -FLT_MAX) continue;
+        float li = mlp[s * 2 + 1];
+        float nm = fmaxf(m, mi);
+        float corr = __expf(m - nm);
+        float e = __expf(mi - nm);
+        acc = acc * corr + op[s * p.head_dim + d] * e;
+        l = l * corr + li * e;
+        m = nm;
+    }
+
+    float inv = (l > 1e-20f) ? (1.0f / l) : 0.0f;
+    p.o[(size_t)bh * p.head_dim + d] = __float2bfloat16(acc * inv);
+}

csrc/kernels/attn_decode_split_kv_mma.cuh

@@ -0,0 +1,200 @@
+#pragma once
+#include <cfloat>
+#include <cuda_bf16.h>
+#include "attn_common.h"
+#include "attn_mma_utils.cuh"
+
+using bf16 = __nv_bfloat16;
+
+// Split-K (FlashDecoding) tensor-core decode via GQA head-packing.
+//
+// Decode has q_len == 1, so S = q @ K^T is a GEMV per head — no tensor-core work
+// on its own. But GQA gives us G = q_head / kv_head query heads that all share
+// one kv_head. We pack those G heads into the M=16 rows of mma.sync.m16n8k16,
+// turning G independent GEMVs into a single GEMM that reuses each loaded K/V tile
+// across all G heads (K/V load is the decode bottleneck, so the reuse is the win,
+// not the flops). The KV sequence is partitioned across gridDim.z blocks so that
+// a decode with only batch*kv_head independent tasks can fill all SMs. Each
+// (batch, kv_head, split) block computes an UN-normalised partial (Oacc, m, l)
+// over its KV slice; the combine kernel below reduces across splits. Fixes the
+// "grid too small" bottleneck (0.04 waves/SM → many blocks) for long-context,
+// small-batch decode.
+//
+// Partial layout (float, contiguous):
+//   o_part : [batch, q_head, num_splits, HEAD_DIM]
+//   ml_part: [batch, q_head, num_splits, 2]  (m, l)
+//
+// Optimizations:
+//   - cp.async global→shared for K/V (bypasses registers, cuts instruction count)
+//   - XOR swizzle (swiz_col): LD=HEAD_DIM, zero waste, no bank conflicts
+//   - Q loaded directly from global into mma A-operand registers (no sQ staging,
+//     no prologue syncwarp) — frees shared memory for double-buffering
+//   - Double-buffered KV (STAGES=2): next tile's cp.async overlaps current
+//     tile's MMA compute — hides global load latency / boosts bandwidth
+//     utilization for small-batch (low-occupancy) decode
+//   - Predicated cp.async (cp_async_16_pred) for full AND partial tiles on one
+//     uniform path — eliminates the scalar fallback branch
+//
+// Smem footprint (BC=32): STAGES=2 → 2*(sK+sV) = 2*2*32*HEAD_DIM*2 bytes.
+//   D=128: 16 KB (fits 48 KB static cap).  D=256: 32 KB (also fits).
+// STAGES=1 fallback (4/8 KB) for smem-constrained configs.
+template <int HEAD_DIM, int BC, int STAGES = 2>
+__global__ void attn_decode_split_kv_mma_kernel(AttentionParams<bf16> p) {
+    constexpr int BR = 16;
+    constexpr int KD = HEAD_DIM / 16;
+    constexpr int NC8 = BC / 8;
+    constexpr int KT2 = BC / 16;
+    constexpr int DN8 = HEAD_DIM / 8;
+    constexpr int LD = HEAD_DIM;
+    constexpr int SWIZ_MASK = (HEAD_DIM >= 64) ? 7 : (HEAD_DIM / 8 - 1);
+    constexpr int VEC = 8;
+    constexpr int TOTAL = BC * HEAD_DIM;
+
+    const int lane = threadIdx.x;
+    const int gid = lane >> 2;
+    const int tid4 = lane & 3;
+
+    const int kv_head = blockIdx.x;
+    const int batch = blockIdx.y;
+    const int split = blockIdx.z;
+    const int G = p.q_head / p.kv_head;
+    const int q_head0 = kv_head * G;
+
+    // Double-buffered shared memory for K/V (no sQ needed — Q goes direct
+    // from global to registers).
+    __shared__ __align__(16) bf16 sK[STAGES * BC * LD];
+    __shared__ __align__(16) bf16 sV[STAGES * BC * LD];
+
+    // ---- Load Q directly from global into mma A-operand registers ----
+    // Same layout as prefill: frag[0]/[2] = row gid, frag[1]/[3] = row gid+8
+    // cols kt*16 + tid4*2 + {0,1} / +{8,9}. pau[0]=cols c,c+1; pau[4]=c+8,c+9.
+    const int q_base = (batch * p.q_head + q_head0) * HEAD_DIM;
+    const int qra = gid;
+    const int qrb = gid + 8;
+    const bool va = qra < G, vb = qrb < G;
+    unsigned Qa[KD][4];
+#pragma unroll
+    for (int kt = 0; kt < KD; kt++) {
+        int c = kt * 16 + tid4 * 2;
+        const unsigned* pau = reinterpret_cast<const unsigned*>(
+            &p.q[q_base + qra * HEAD_DIM + c]);
+        const unsigned* pbu = reinterpret_cast<const unsigned*>(
+            &p.q[q_base + qrb * HEAD_DIM + c]);
+        Qa[kt][0] = va ? pau[0] : 0u;
+        Qa[kt][1] = vb ? pbu[0] : 0u;
+        Qa[kt][2] = va ? pau[4] : 0u;
+        Qa[kt][3] = vb ? pbu[4] : 0u;
+    }
+
+    float Oacc[DN8][4];
+#pragma unroll
+    for (int j = 0; j < DN8; j++)
+        Oacc[j][0] = Oacc[j][1] = Oacc[j][2] = Oacc[j][3] = 0.0f;
+    float m0 = -FLT_MAX, m1 = -FLT_MAX, l0 = 0.0f, l1 = 0.0f;
+
+    const int kv_base = (batch * p.kv_head + kv_head) * p.kv_len * HEAD_DIM;
+    const int mask_base = batch * p.kv_len;
+    const int tiles_total = (p.kv_len + BC - 1) / BC;
+    const int tiles_per_split = (tiles_total + p.num_splits - 1) / p.num_splits;
+    const int ti_begin = split * tiles_per_split;
+    const int ti_end = min(tiles_total, ti_begin + tiles_per_split);
+    const int has_mask = p.use_mask && p.mask;
+
+    // ---- Load tile lambda: predicated cp.async, unified full/partial ----
+    auto load_tile = [&](int ti, int buf) {
+        int kv0 = ti * BC;
+        bf16* dK = sK + buf * BC * LD;
+        bf16* dV = sV + buf * BC * LD;
+#pragma unroll
+        for (int i = lane * VEC; i < TOTAL; i += 32 * VEC) {
+            int r = i / HEAD_DIM, d = i % HEAD_DIM;
+            int kc = kv0 + r;
+            bool valid = kc < p.kv_len;
+            int off = r * LD + swiz_col(d, r, SWIZ_MASK);
+            cp_async_16_pred(&dK[off], &p.k[kv_base + kc * HEAD_DIM + d], valid);
+            cp_async_16_pred(&dV[off], &p.v[kv_base + kc * HEAD_DIM + d], valid);
+        }
+        cp_async_commit();
+    };
+
+    // ---- Prologue: issue first tile load ----
+    if (ti_begin < ti_end) {
+        load_tile(ti_begin, 0);
+    }
+
+    for (int ti = ti_begin; ti < ti_end; ti++) {
+        constexpr int BUF_MASK = (STAGES > 1) ? (STAGES - 1) : 0;
+        int buf = (ti - ti_begin) & BUF_MASK;
+
+        // Wait for current tile, then issue next tile's prefetch (overlaps
+        // with this tile's compute). Single syncwarp covers both hazards.
+        // When STAGES==1, no prefetch — load happens at end of prior iter.
+        cp_async_wait_group<0>();
+        __syncwarp();
+        if constexpr (STAGES > 1) {
+            if (ti + 1 < ti_end)
+                load_tile(ti + 1, (ti + 1 - ti_begin) & BUF_MASK);
+        }
+
+        const bf16* bK = sK + buf * BC * LD;
+        const bf16* bV = sV + buf * BC * LD;
+        int kv0 = ti * BC;
+
+        float Sacc[NC8][4];
+        mma_compute_scores<KD, NC8>(Qa, bK, LD, SWIZ_MASK, lane, Sacc);
+
+        #pragma unroll
+        for (int n8 = 0; n8 < NC8; n8++)
+            Sacc[n8][0] *= p.scale, Sacc[n8][1] *= p.scale,
+            Sacc[n8][2] *= p.scale, Sacc[n8][3] *= p.scale;
+
+        int maxc = p.is_causal ? min(p.kv_len, p.causal_offset + 1) : p.kv_len;
+        mma_softmax_tile<NC8, DN8>(kv0, maxc, maxc,
+                                    mask_base, p.mask, has_mask,
+                                    Sacc, Oacc, m0, m1, l0, l1, lane);
+
+        mma_pv_accumulate<DN8, KT2>(Sacc, bV, LD, SWIZ_MASK, lane, Oacc);
+        __syncwarp();
+
+        if constexpr (STAGES == 1) {
+            if (ti + 1 < ti_end)
+                load_tile(ti + 1, 0);
+        }
+    }
+
+    // ---- write UN-normalised partials for this split ----
+    auto split_slot = [&](int h) -> size_t {
+        size_t bh = (size_t)batch * p.q_head + h;
+        return bh * p.num_splits + split;
+    };
+#pragma unroll
+    for (int dn8 = 0; dn8 < DN8; dn8++) {
+        int d = dn8 * 8 + 2 * tid4;
+        int r0 = gid, r1 = gid + 8;
+        if (r0 < G) {
+            int h = q_head0 + r0;
+            float* op = p.o_part + split_slot(h) * HEAD_DIM;
+            op[d] = Oacc[dn8][0];
+            op[d + 1] = Oacc[dn8][1];
+        }
+        if (r1 < G) {
+            int h = q_head0 + r1;
+            float* op = p.o_part + split_slot(h) * HEAD_DIM;
+            op[d] = Oacc[dn8][2];
+            op[d + 1] = Oacc[dn8][3];
+        }
+    }
+    if (tid4 == 0) {
+        int r0 = gid, r1 = gid + 8;
+        if (r0 < G) {
+            int h = q_head0 + r0;
+            float* mp = p.ml_part + split_slot(h) * 2;
+            mp[0] = m0; mp[1] = l0;
+        }
+        if (r1 < G) {
+            int h = q_head0 + r1;
+            float* mp = p.ml_part + split_slot(h) * 2;
+            mp[0] = m1; mp[1] = l1;
+        }
+    }
+}

csrc/kernels/attn_entry_utils.cuh

@@ -0,0 +1,43 @@
+#pragma once
+#include <torch/extension.h>
+#include "attn_common.h"
+
+template<typename T>
+inline void attn_pack_params(
+    torch::Tensor q,
+    torch::Tensor k,
+    torch::Tensor v,
+    c10::optional<torch::Tensor> mask,
+    bool is_causal,
+    int64_t causal_offset,
+    c10::optional<double> scale,
+    AttentionParams<T>& p
+) {
+    TORCH_CHECK(q.is_cuda() && k.is_cuda() && v.is_cuda());
+    TORCH_CHECK(q.dtype() == torch::kBFloat16);
+    TORCH_CHECK(k.dtype() == torch::kBFloat16);
+    TORCH_CHECK(v.dtype() == torch::kBFloat16);
+
+    p.batch = (int)q.size(0);
+    p.q_head = (int)q.size(1);
+    p.kv_head = (int)k.size(1);
+    p.q_len = (int)q.size(2);
+    p.kv_len = (int)k.size(2);
+    p.head_dim = (int)q.size(3);
+    p.use_mask = mask.has_value() ? 1 : 0;
+    p.is_causal = is_causal ? 1 : 0;
+    p.causal_offset = (int)causal_offset;
+    p.scale = scale.has_value() ? (float)scale.value() : 1.0f / sqrtf((float)p.head_dim);
+    p.q = (const T*)q.data_ptr();
+    p.k = (const T*)k.data_ptr();
+    p.v = (const T*)v.data_ptr();
+    if (p.use_mask) {
+        TORCH_CHECK(mask.value().dtype() == torch::kBool);
+        TORCH_CHECK(mask.value().dim() == 2);
+        TORCH_CHECK(mask.value().size(0) == p.batch);
+        TORCH_CHECK(mask.value().size(1) == p.kv_len);
+        p.mask = mask.value().data_ptr<bool>();
+    } else {
+        p.mask = nullptr;
+    }
+}

csrc/kernels/attn_mma_utils.cuh

@@ -0,0 +1,252 @@
+#pragma once
+#include <cfloat>
+#include <cuda_fp16.h>
+#include <cuda_runtime.h>
+
+// Shared MMA utilities for tensor-core GQA kernels.
+// mma.sync.m16n8k16 PTX wrappers, ldmatrix helpers, and bf16 packing.
+
+// mma.sync.aligned.m16n8k16.row.col.f32.bf16.bf16.f32
+__device__ __forceinline__ void mma16816(float* d, const unsigned* a,
+                                          const unsigned* b, const float* c) {
+    asm volatile(
+        "mma.sync.aligned.m16n8k16.row.col.f32.bf16.bf16.f32 "
+        "{%0,%1,%2,%3}, {%4,%5,%6,%7}, {%8,%9}, {%10,%11,%12,%13};"
+        : "=f"(d[0]), "=f"(d[1]), "=f"(d[2]), "=f"(d[3])
+        : "r"(a[0]), "r"(a[1]), "r"(a[2]), "r"(a[3]), "r"(b[0]), "r"(b[1]),
+          "f"(c[0]), "f"(c[1]), "f"(c[2]), "f"(c[3]));
+}
+
+// read two adjacent bf16 from smem as one packed .b32 (elem0 low, elem1 high)
+__device__ __forceinline__ unsigned ld2(const bf16* p) {
+    return *reinterpret_cast<const unsigned*>(p);
+}
+
+// pack two floats into one bf16x2 as .b32
+__device__ __forceinline__ unsigned pk2(float a, float b) {
+    __nv_bfloat162 v = __floats2bfloat162_rn(a, b);
+    return *reinterpret_cast<unsigned*>(&v);
+}
+
+// pack two (non-contiguous) bf16 into one .b32
+__device__ __forceinline__ unsigned pkb(bf16 a, bf16 b) {
+    __nv_bfloat162 v;
+    v.x = a;
+    v.y = b;
+    return *reinterpret_cast<unsigned*>(&v);
+}
+
+// ldmatrix: cooperatively load mma fragments from smem (one instruction per
+// 16x16 / 16x8 tile) with the exact register layout mma expects — replaces the
+// scalar per-thread fragment packing, cutting shared-load instructions and bank
+// conflicts. Each lane supplies the shared address of one 8-wide row.
+__device__ __forceinline__ void ldmatrix_x4(unsigned* r, const bf16* p) {
+    unsigned a = __cvta_generic_to_shared(p);
+    asm volatile("ldmatrix.sync.aligned.m8n8.x4.shared.b16 {%0,%1,%2,%3}, [%4];"
+                 : "=r"(r[0]), "=r"(r[1]), "=r"(r[2]), "=r"(r[3])
+                 : "r"(a));
+}
+__device__ __forceinline__ void ldmatrix_x2(unsigned* r, const bf16* p) {
+    unsigned a = __cvta_generic_to_shared(p);
+    asm volatile("ldmatrix.sync.aligned.m8n8.x2.shared.b16 {%0,%1}, [%2];"
+                 : "=r"(r[0]), "=r"(r[1])
+                 : "r"(a));
+}
+__device__ __forceinline__ void ldmatrix_x2_trans(unsigned* r, const bf16* p) {
+    unsigned a = __cvta_generic_to_shared(p);
+    asm volatile("ldmatrix.sync.aligned.m8n8.x2.trans.shared.b16 {%0,%1}, [%2];"
+                 : "=r"(r[0]), "=r"(r[1])
+                 : "r"(a));
+}
+
+// XOR swizzle for shared-memory column at 8-bf16 chunk granularity.
+// Eliminates ldmatrix bank conflicts without LD padding: consecutive rows
+// land in distinct bank groups. swiz_col(d, r, mask) = ((d>>3)^(r&mask))<<3 | (d&7).
+// mask must cover log2(HEAD_DIM/8) chunk bits but stay within LD: use 7 for
+// HEAD_DIM>=64 (8+ chunks), 3 for HEAD_DIM=32 (4 chunks). Default 7 keeps
+// existing HEAD_DIM>=64 call sites working unchanged.
+__device__ __forceinline__ int swiz_col(int d, int r, int mask = 7) {
+    return ((d >> 3) ^ (r & mask)) << 3 | (d & 7);
+}
+
+// cp.async: copy 16 bytes (8 bf16) from global to shared memory directly,
+// bypassing registers. Eliminates shared-store bank conflicts and cuts
+// load-loop instruction count in half (1 cp.async vs 1 LDG + 1 STS).
+// Requires sm_80+.
+__device__ __forceinline__ void cp_async_16(bf16* smem_ptr, const void* gmem_ptr) {
+    unsigned smem_addr = __cvta_generic_to_shared(smem_ptr);
+    asm volatile("cp.async.ca.shared.global [%0], [%1], 16;"
+                 :: "r"(smem_addr), "l"(gmem_ptr));
+}
+
+// Predicated cp.async: copy 16 bytes when `pred`, otherwise zero-fill the
+// destination (src-size operand = 0 → no bytes read from src, so an
+// out-of-bounds src address is never dereferenced). Lets full and partial
+// tiles share one uniform async load path — no scalar fallback branch.
+__device__ __forceinline__ void cp_async_16_pred(bf16* smem_ptr,
+                                                  const void* gmem_ptr,
+                                                  bool pred) {
+    unsigned smem_addr = __cvta_generic_to_shared(smem_ptr);
+    int src_size = pred ? 16 : 0;
+    asm volatile("cp.async.ca.shared.global [%0], [%1], 16, %2;"
+                 :: "r"(smem_addr), "l"(gmem_ptr), "r"(src_size));
+}
+
+__device__ __forceinline__ void cp_async_commit() {
+    asm volatile("cp.async.commit_group;");
+}
+
+__device__ __forceinline__ void cp_async_wait_all() {
+    asm volatile("cp.async.wait_all;");
+}
+
+// Wait until at most N commit groups are still in flight. Used for
+// double-buffered pipelining: wait_group<1> lets the next tile's cp.async
+// continue while ensuring the current tile's data is ready.
+template <int N>
+__device__ __forceinline__ void cp_async_wait_group() {
+    asm volatile("cp.async.wait_group %0;" :: "n"(N));
+}
+
+// ---------------------------------------------------------------------------
+// Shared MMA compute functions — used by both decode and prefill MMA kernels.
+// Extracted because S=Q@K^T, online softmax, and P@V are structurally identical
+// between the two kernels; only the per-row causal/mask bounds differ.
+// ---------------------------------------------------------------------------
+
+// S = Q @ K^T  (Qa pre-loaded by the caller; scale applied post-mma in the
+// caller to avoid bf16 precision loss).
+// LD and SWIZ_MASK are constexpr in the calling kernel — passing them as
+// runtime ints lets the compiler fold them while keeping the signature clean.
+template <int KD, int NC8>
+__device__ inline void mma_compute_scores(
+    const unsigned Qa[KD][4],
+    const bf16* __restrict__ sK,
+    int LD, int SWIZ_MASK, int lane,
+    float Sacc[NC8][4])
+{
+    #pragma unroll
+    for (int n8 = 0; n8 < NC8; n8++) {
+        Sacc[n8][0] = Sacc[n8][1] = Sacc[n8][2] = Sacc[n8][3] = 0.0f;
+        int krow_l = n8 * 8 + (lane & 7);
+        int kcol_h = (lane & 8) ? 8 : 0;
+        #pragma unroll
+        for (int kt = 0; kt < KD; kt++) {
+            unsigned b[2];
+            ldmatrix_x2(b, &sK[krow_l * LD + swiz_col(kt * 16 + kcol_h, krow_l, SWIZ_MASK)]);
+            mma16816(Sacc[n8], Qa[kt], b, Sacc[n8]);
+        }
+    }
+}
+
+// Online softmax + Oacc rescale for one K/V tile.
+//   maxc0/maxc1: per-row KV column bounds (prefill: per-query-row causal limits;
+//     decode: same value for both rows since q_len==1).
+// Reads Sacc (Q@K^T scores), applies causal/mask, computes P = exp(S - nm),
+// rescales Oacc by exp(m_old - nm), and updates m/l — all in place.
+template <int NC8, int DN8>
+__device__ inline void mma_softmax_tile(
+    int kv0,
+    int maxc0, int maxc1,
+    int mask_base,
+    const bool* __restrict__ mask,
+    bool has_mask,
+    float Sacc[NC8][4],
+    float Oacc[DN8][4],
+    float& m0, float& m1,
+    float& l0, float& l1,
+    int lane)
+{
+    int tid4 = lane & 3;
+
+    // Mask out-of-bounds / masked columns: set -FLT_MAX so expf → 0 downstream
+    // without per-element sentinel checks. Compute tile-local row maxima.
+    float rmax0 = -FLT_MAX, rmax1 = -FLT_MAX;
+    #pragma unroll
+    for (int n8 = 0; n8 < NC8; n8++) {
+        int cc = kv0 + n8 * 8 + 2 * tid4;
+        int c1 = cc + 1;
+        bool b0 = (cc >= maxc0) || (has_mask && !mask[mask_base + cc]);
+        bool b1 = (c1 >= maxc0) || (has_mask && !mask[mask_base + c1]);
+        bool b2 = (cc >= maxc1) || (has_mask && !mask[mask_base + cc]);
+        bool b3 = (c1 >= maxc1) || (has_mask && !mask[mask_base + c1]);
+        float s0 = b0 ? -FLT_MAX : Sacc[n8][0];
+        float s1 = b1 ? -FLT_MAX : Sacc[n8][1];
+        float s2 = b2 ? -FLT_MAX : Sacc[n8][2];
+        float s3 = b3 ? -FLT_MAX : Sacc[n8][3];
+        Sacc[n8][0] = s0; Sacc[n8][1] = s1;
+        Sacc[n8][2] = s2; Sacc[n8][3] = s3;
+        rmax0 = fmaxf(rmax0, fmaxf(s0, s1));
+        rmax1 = fmaxf(rmax1, fmaxf(s2, s3));
+    }
+    // Warp-reduce row maxima across the 4-lane thread group (xor 1, xor 2).
+    rmax0 = fmaxf(rmax0, __shfl_xor_sync(0xFFFFFFFF, rmax0, 1));
+    rmax0 = fmaxf(rmax0, __shfl_xor_sync(0xFFFFFFFF, rmax0, 2));
+    rmax1 = fmaxf(rmax1, __shfl_xor_sync(0xFFFFFFFF, rmax1, 1));
+    rmax1 = fmaxf(rmax1, __shfl_xor_sync(0xFFFFFFFF, rmax1, 2));
+
+    // nm = max(running max m, tile-local max rmax) — updated running maximum.
+    float nm0 = fmaxf(m0, rmax0), nm1 = fmaxf(m1, rmax1);
+    // corr rescales Oacc and l by exp(m_old - nm). When all-masked (m == nm ==
+    // -FLT_MAX), exp(0) = 1 — correct, no guard needed.
+    float corr0 = __expf(m0 - nm0);
+    float corr1 = __expf(m1 - nm1);
+    // pn guards only the all-masked-row edge: if nm == -FLT_MAX, exp(S - nm)
+    // gives 1 not 0 for masked entries. Two scalar masks replace 4*NC8
+    // per-element comparisons.
+    float pn0 = (nm0 == -FLT_MAX) ? 0.0f : 1.0f;
+    float pn1 = (nm1 == -FLT_MAX) ? 0.0f : 1.0f;
+
+    // P = exp(S - nm) for each element. Masked entries (Sacc = -FLT_MAX) give
+    // exp(-inf) ≈ 0 naturally; pn zero-fills the all-masked-row edge.
+    float rsum0 = 0.0f, rsum1 = 0.0f;
+    #pragma unroll
+    for (int n8 = 0; n8 < NC8; n8++) {
+        float p0 = pn0 * __expf(Sacc[n8][0] - nm0);
+        float p1 = pn0 * __expf(Sacc[n8][1] - nm0);
+        float p2 = pn1 * __expf(Sacc[n8][2] - nm1);
+        float p3 = pn1 * __expf(Sacc[n8][3] - nm1);
+        Sacc[n8][0] = p0; Sacc[n8][1] = p1;
+        Sacc[n8][2] = p2; Sacc[n8][3] = p3;
+        rsum0 += p0 + p1;
+        rsum1 += p2 + p3;
+    }
+    rsum0 += __shfl_xor_sync(0xFFFFFFFF, rsum0, 1);
+    rsum0 += __shfl_xor_sync(0xFFFFFFFF, rsum0, 2);
+    rsum1 += __shfl_xor_sync(0xFFFFFFFF, rsum1, 1);
+    rsum1 += __shfl_xor_sync(0xFFFFFFFF, rsum1, 2);
+    l0 = l0 * corr0 + rsum0;
+    l1 = l1 * corr1 + rsum1;
+    m0 = nm0; m1 = nm1;
+
+    #pragma unroll
+    for (int j = 0; j < DN8; j++) {
+        Oacc[j][0] *= corr0; Oacc[j][1] *= corr0;
+        Oacc[j][2] *= corr1; Oacc[j][3] *= corr1;
+    }
+}
+
+// O += P @ V  (Sacc must contain P = attention weights after softmax).
+template <int DN8, int KT2>
+__device__ inline void mma_pv_accumulate(
+    float Sacc[][4],
+    const bf16* __restrict__ sV,
+    int LD, int SWIZ_MASK, int lane,
+    float Oacc[DN8][4])
+{
+    #pragma unroll
+    for (int kt2 = 0; kt2 < KT2; kt2++) {
+        unsigned Pa[4];
+        Pa[0] = pk2(Sacc[kt2 * 2][0], Sacc[kt2 * 2][1]);
+        Pa[1] = pk2(Sacc[kt2 * 2][2], Sacc[kt2 * 2][3]);
+        Pa[2] = pk2(Sacc[kt2 * 2 + 1][0], Sacc[kt2 * 2 + 1][1]);
+        Pa[3] = pk2(Sacc[kt2 * 2 + 1][2], Sacc[kt2 * 2 + 1][3]);
+        int vrow_l = kt2 * 16 + (lane & 15);
+        #pragma unroll
+        for (int dn8 = 0; dn8 < DN8; dn8++) {
+            unsigned b[2];
+            ldmatrix_x2_trans(b, &sV[vrow_l * LD + swiz_col(dn8 * 8, vrow_l, SWIZ_MASK)]);
+            mma16816(Oacc[dn8], Pa, b, Oacc[dn8]);
+        }
+    }
+}

csrc/kernels/attn_paged_decode.cu

@@ -0,0 +1,149 @@
+#include "attn_paged_decode_split_kv.cuh"
+#ifndef ASTRAI_NO_MMA
+#include "attn_paged_decode_split_kv_mma.cuh"
+#endif
+
+#include <torch/extension.h>
+#include <c10/cuda/CUDAGuard.h>
+
+static int paged_decode_num_splits(int base_blocks, int tiles_total) {
+    int sm_count = 0;
+    cudaDeviceGetAttribute(&sm_count, cudaDevAttrMultiProcessorCount, 0);
+    int n = (2 * sm_count + base_blocks - 1) / base_blocks;
+    return std::max(1, std::min(n, std::min(tiles_total, 32)));
+}
+
+static void launch_paged_scalar_decode(PagedAttentionParams<bf16>& p) {
+    int group_size = p.q_head / p.kv_head;
+    int chunks_total = (p.kv_len + PDC_CHUNK - 1) / PDC_CHUNK;
+    p.num_splits = paged_decode_num_splits(p.batch * p.kv_head, chunks_total);
+
+    auto fopt = torch::TensorOptions().dtype(torch::kFloat32).device(torch::kCUDA);
+    auto o_part = torch::empty({p.batch, p.q_head, p.num_splits, p.head_dim}, fopt);
+    auto ml_part = torch::empty({p.batch, p.q_head, p.num_splits, 2}, fopt);
+    p.o_part = o_part.data_ptr<float>();
+    p.ml_part = ml_part.data_ptr<float>();
+
+    size_t smem = PDC_CHUNK * p.head_dim * sizeof(bf16);
+    dim3 grid = dim3(p.batch * p.kv_head, 1, p.num_splits);
+    dim3 block = dim3(32, group_size);
+    paged_attn_decode_split_kv_kernel<<<grid, block, smem>>>(p);
+    paged_attn_decode_combine_kernel<<<p.batch * p.q_head, p.head_dim>>>(p);
+}
+
+#ifndef ASTRAI_NO_MMA
+template <int HEAD_DIM, int BC, int STAGES = (HEAD_DIM <= 128) ? 2 : 1>
+static void launch_paged_mma_decode(PagedAttentionParams<bf16>& p) {
+    int tiles_total = (p.kv_len + BC - 1) / BC;
+    p.num_splits = paged_decode_num_splits(p.batch * p.kv_head, tiles_total);
+
+    auto fopt = torch::TensorOptions().dtype(torch::kFloat32).device(torch::kCUDA);
+    auto o_part = torch::empty({p.batch, p.q_head, p.num_splits, p.head_dim}, fopt);
+    auto ml_part = torch::empty({p.batch, p.q_head, p.num_splits, 2}, fopt);
+    p.o_part = o_part.data_ptr<float>();
+    p.ml_part = ml_part.data_ptr<float>();
+
+    paged_attn_decode_split_kv_mma_kernel<HEAD_DIM, BC, STAGES><<<dim3(p.kv_head, p.batch, p.num_splits), 32>>>(p);
+    paged_attn_decode_combine_kernel<<<p.batch * p.q_head, p.head_dim>>>(p);
+}
+#endif
+
+template <int HEAD_DIM>
+static void dispatch_paged_decode(PagedAttentionParams<bf16>& p) {
+#ifndef ASTRAI_NO_MMA
+    int G = p.q_head / p.kv_head;
+    if (!p.use_mask && G >= 1 && G <= 16 && p.page_size >= 32) {
+        launch_paged_mma_decode<HEAD_DIM, 32>(p);
+        return;
+    }
+#endif
+    launch_paged_scalar_decode(p);
+}
+
+torch::Tensor attn_paged_decode(
+    torch::Tensor q,
+    torch::Tensor page_table,
+    torch::Tensor k_cache,
+    torch::Tensor v_cache,
+    int64_t page_size,
+    int64_t kv_len,
+    c10::optional<torch::Tensor> mask,
+    bool is_causal = false,
+    int64_t causal_offset = 0,
+    c10::optional<double> scale = c10::nullopt
+) {
+    const at::cuda::OptionalCUDAGuard device_guard(device_of(q));
+
+    int batch = q.size(0);
+    int q_head = q.size(1);
+    int head_dim = q.size(3);
+    int kv_head = k_cache.size(2);
+    int max_pages = page_table.size(1);
+
+    TORCH_CHECK(q.is_cuda() && page_table.is_cuda() && k_cache.is_cuda() && v_cache.is_cuda());
+    TORCH_CHECK(q.dtype() == torch::kBFloat16, "q must be bf16");
+    TORCH_CHECK(k_cache.dtype() == torch::kBFloat16, "k_cache must be bf16");
+    TORCH_CHECK(v_cache.dtype() == torch::kBFloat16, "v_cache must be bf16");
+    TORCH_CHECK(page_table.dtype() == torch::kLong, "page_table must be int64");
+    TORCH_CHECK(q.size(2) == 1, "Q seq_len must be 1 (decode)");
+    TORCH_CHECK(head_dim % 32 == 0, "head_dim must be multiple of 32");
+    TORCH_CHECK(k_cache.size(1) == page_size,
+                "k_cache dim 1 must equal page_size, got ",
+                k_cache.size(1), " vs ", page_size);
+    TORCH_CHECK(k_cache.size(0) >= 0, "k_cache must have at least 0 pages");
+
+    float scale_val = scale.has_value()
+        ? static_cast<float>(scale.value())
+        : 1.0f / std::sqrt(static_cast<float>(head_dim));
+
+    auto O = torch::empty_like(q);
+
+    PagedAttentionParams<bf16, float> p;
+    p.batch = batch;
+    p.q_head = q_head;
+    p.kv_head = kv_head;
+    p.q_len = static_cast<int>(q.size(2));
+    p.kv_len = static_cast<int>(kv_len);
+    p.head_dim = head_dim;
+    p.use_mask = (mask.has_value() && mask.value().defined()) ? 1 : 0;
+    p.is_causal = is_causal ? 1 : 0;
+    p.causal_offset = static_cast<int>(causal_offset);
+    p.scale = scale_val;
+    p.page_size = static_cast<int>(page_size);
+    p.max_pages = max_pages;
+    p.page_table = page_table.data_ptr<int64_t>();
+    p.k_cache = reinterpret_cast<const bf16*>(k_cache.data_ptr());
+    p.v_cache = reinterpret_cast<const bf16*>(v_cache.data_ptr());
+    p.q = reinterpret_cast<const bf16*>(q.data_ptr());
+    p.mask = p.use_mask ? mask.value().data_ptr<bool>() : nullptr;
+    p.o = reinterpret_cast<bf16*>(O.data_ptr());
+    p.o_part = nullptr;
+    p.ml_part = nullptr;
+
+    switch (p.head_dim) {
+        case 32:  dispatch_paged_decode<32>(p);  break;
+        case 64:  dispatch_paged_decode<64>(p);  break;
+        case 128: dispatch_paged_decode<128>(p); break;
+        case 256: dispatch_paged_decode<256>(p); break;
+        default:
+            TORCH_CHECK(false, "paged_decode: unsupported head_dim ", p.head_dim,
+                        " (supported: 32, 64, 128, 256)");
+    }
+
+    return O;
+}
+
+PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
+    m.def("attn_paged_decode", &attn_paged_decode,
+        py::arg("q"),
+        py::arg("page_table"),
+        py::arg("k_cache"),
+        py::arg("v_cache"),
+        py::arg("page_size"),
+        py::arg("kv_len"),
+        py::arg("mask") = py::none(),
+        py::arg("is_causal") = false,
+        py::arg("causal_offset") = 0,
+        py::arg("scale") = py::none(),
+        "Paged GQA decode — split-KV with direct page-table access.");
+}

csrc/kernels/attn_paged_decode_split_kv.cuh

@@ -0,0 +1,140 @@
+#pragma once
+#include <cuda_bf16.h>
+#include <float.h>
+#include "attn_common.h"
+
+using bf16 = __nv_bfloat16;
+constexpr int PDC_CHUNK = 64;
+
+__device__ inline float paged_warp_reduce_sum(float val) {
+    for (int offset = 16; offset > 0; offset >>= 1)
+        val += __shfl_xor_sync(0xFFFFFFFF, val, offset);
+    return val;
+}
+
+// Split-KV scalar decode: one warp per query head, grid.z partitions KV.
+__global__ void paged_attn_decode_split_kv_kernel(PagedAttentionParams<bf16> p) {
+    int batch = blockIdx.x / p.kv_head;
+    int kv_head = blockIdx.x % p.kv_head;
+    int split = blockIdx.z;
+    int group_size = blockDim.y;
+    int q_head = kv_head * group_size + threadIdx.y;
+    int lane = threadIdx.x;
+    int hd_per_thread = p.head_dim / 32;
+
+    float q_reg[8];
+    int q_off = ((batch * p.q_head + q_head) * 1) * p.head_dim + lane * hd_per_thread;
+    #pragma unroll
+    for (int i = 0; i < hd_per_thread; i++)
+        q_reg[i] = __bfloat162float(p.q[q_off + i]);
+
+    float m = -FLT_MAX, d = 0.0f, acc_reg[8] = {0.0f};
+
+    extern __shared__ __align__(16) bf16 k_smem[];
+
+    int chunks_total = (p.kv_len + PDC_CHUNK - 1) / PDC_CHUNK;
+    int chunks_per_split = (chunks_total + p.num_splits - 1) / p.num_splits;
+    int ch_begin = split * chunks_per_split;
+    int ch_end = min(chunks_total, ch_begin + chunks_per_split);
+
+    const int mask_base = batch * p.kv_len;
+
+    for (int ci = ch_begin; ci < ch_end; ci++) {
+        int chunk_start = ci * PDC_CHUNK;
+        int this_chunk = min(PDC_CHUNK, p.kv_len - chunk_start);
+
+        int total = this_chunk * p.head_dim;
+        for (int i = threadIdx.y * 32 + lane; i < total; i += blockDim.x * blockDim.y) {
+            int s = i / p.head_dim;
+            int d_dim = i % p.head_dim;
+            int pos = chunk_start + s;
+            int logical_page = pos / p.page_size;
+            int page_offset = pos % p.page_size;
+            int phys_page = p.page_table[batch * p.max_pages + logical_page];
+            if (phys_page >= 0) {
+                int64_t off = (int64_t)phys_page * p.page_size * p.kv_head * p.head_dim
+                            + (int64_t)page_offset * p.kv_head * p.head_dim
+                            + (int64_t)kv_head * p.head_dim
+                            + d_dim;
+                k_smem[i] = p.k_cache[off];
+            } else {
+                k_smem[i] = __float2bfloat16(0.0f);
+            }
+        }
+        __syncthreads();
+
+        for (int s = 0; s < this_chunk; s++) {
+            float partial = 0.0f;
+            #pragma unroll
+            for (int i = 0; i < hd_per_thread; i++)
+                partial += q_reg[i] * __bfloat162float(k_smem[s * p.head_dim + lane * hd_per_thread + i]);
+            partial = paged_warp_reduce_sum(partial) * p.scale;
+
+            if (p.use_mask && p.mask && !p.mask[mask_base + chunk_start + s])
+                partial = -FLT_MAX;
+            if (p.is_causal && (chunk_start + s) > p.causal_offset)
+                partial = -FLT_MAX;
+
+            float new_m = fmaxf(m, partial);
+            float alpha = expf(m - new_m);
+            float beta  = expf(partial - new_m);
+            d = d * alpha + beta;
+
+            int pos = chunk_start + s;
+            int logical_page = pos / p.page_size;
+            int page_offset = pos % p.page_size;
+            int phys_page = p.page_table[batch * p.max_pages + logical_page];
+            if (phys_page >= 0) {
+                int64_t v_base = (int64_t)phys_page * p.page_size * p.kv_head * p.head_dim
+                               + (int64_t)page_offset * p.kv_head * p.head_dim
+                               + (int64_t)kv_head * p.head_dim;
+                #pragma unroll
+                for (int i = 0; i < hd_per_thread; i++)
+                    acc_reg[i] = acc_reg[i] * alpha + __bfloat162float(p.v_cache[v_base + lane * hd_per_thread + i]) * beta;
+            } else {
+                #pragma unroll
+                for (int i = 0; i < hd_per_thread; i++)
+                    acc_reg[i] = acc_reg[i] * alpha + 0.0f * beta;
+            }
+            m = new_m;
+        }
+        __syncthreads();
+    }
+
+    size_t bh = (size_t)batch * p.q_head + q_head;
+    size_t slot = bh * p.num_splits + split;
+    int d0 = lane * hd_per_thread;
+    #pragma unroll
+    for (int i = 0; i < hd_per_thread; i++)
+        p.o_part[slot * p.head_dim + (d0 + i)] = acc_reg[i];
+    if (lane == 0) {
+        p.ml_part[slot * 2] = m;
+        p.ml_part[slot * 2 + 1] = d;
+    }
+}
+
+__global__ void paged_attn_decode_combine_kernel(PagedAttentionParams<bf16> p) {
+    int bh = blockIdx.x;
+    int d = threadIdx.x;
+    if (d >= p.head_dim) return;
+
+    size_t split_base = (size_t)bh * p.num_splits;
+    const float* mlp = p.ml_part + split_base * 2;
+    const float* op = p.o_part + split_base * p.head_dim;
+
+    float m = -FLT_MAX, l = 0.0f, acc = 0.0f;
+    for (int s = 0; s < p.num_splits; s++) {
+        float mi = mlp[s * 2];
+        if (mi <= -FLT_MAX) continue;
+        float li = mlp[s * 2 + 1];
+        float nm = fmaxf(m, mi);
+        float corr = __expf(m - nm);
+        float e = __expf(mi - nm);
+        acc = acc * corr + op[s * p.head_dim + d] * e;
+        l = l * corr + li * e;
+        m = nm;
+    }
+
+    float inv = (l > 1e-20f) ? (1.0f / l) : 0.0f;
+    p.o[(size_t)bh * p.head_dim + d] = __float2bfloat16(acc * inv);
+}

csrc/kernels/attn_paged_decode_split_kv_mma.cuh

@@ -0,0 +1,182 @@
+#pragma once
+#include <cfloat>
+#include <cuda_bf16.h>
+#include "attn_common.h"
+#include "attn_mma_utils.cuh"
+
+using bf16 = __nv_bfloat16;
+
+// Paged split-KV tensor-core decode via GQA head-packing.
+// Identical algorithm to attn_decode_split_kv_mma_kernel but reads K/V
+// directly from the page pool through a page table, eliminating the gather
+// copy.  Each tile (BC=32) fits within a single page (page_size >= 32), so
+// the page-table lookup happens once per tile for cp.async.
+//
+// Optimizations mirror attn_decode_split_kv_mma_kernel:
+//   - Q loaded directly from global into mma A-operand registers (no sQ)
+//   - Double-buffered KV (STAGES=2) for D<=128, single-buffer for D=256
+//   - Predicated cp.async for unified full/partial tile path
+template <int HEAD_DIM, int BC, int STAGES = (HEAD_DIM <= 128) ? 2 : 1>
+__global__ void paged_attn_decode_split_kv_mma_kernel(PagedAttentionParams<bf16> p) {
+    constexpr int BR = 16;
+    constexpr int KD = HEAD_DIM / 16;
+    constexpr int NC8 = BC / 8;
+    constexpr int KT2 = BC / 16;
+    constexpr int DN8 = HEAD_DIM / 8;
+    constexpr int LD = HEAD_DIM;
+    constexpr int SWIZ_MASK = (HEAD_DIM >= 64) ? 7 : (HEAD_DIM / 8 - 1);
+    constexpr int VEC = 8;
+    constexpr int TOTAL = BC * HEAD_DIM;
+
+    const int lane = threadIdx.x;
+    const int gid = lane >> 2;
+    const int tid4 = lane & 3;
+
+    const int kv_head_idx = blockIdx.x;
+    const int batch = blockIdx.y;
+    const int split = blockIdx.z;
+    const int G = p.q_head / p.kv_head;
+    const int q_head0 = kv_head_idx * G;
+
+    __shared__ __align__(16) bf16 sK[STAGES * BC * LD];
+    __shared__ __align__(16) bf16 sV[STAGES * BC * LD];
+
+    // ---- Load Q directly from global into mma A-operand registers ----
+    const int q_base = (batch * p.q_head + q_head0) * HEAD_DIM;
+    const int qra = gid;
+    const int qrb = gid + 8;
+    const bool va = qra < G, vb = qrb < G;
+    unsigned Qa[KD][4];
+#pragma unroll
+    for (int kt = 0; kt < KD; kt++) {
+        int c = kt * 16 + tid4 * 2;
+        const unsigned* pau = reinterpret_cast<const unsigned*>(
+            &p.q[q_base + qra * HEAD_DIM + c]);
+        const unsigned* pbu = reinterpret_cast<const unsigned*>(
+            &p.q[q_base + qrb * HEAD_DIM + c]);
+        Qa[kt][0] = va ? pau[0] : 0u;
+        Qa[kt][1] = vb ? pbu[0] : 0u;
+        Qa[kt][2] = va ? pau[4] : 0u;
+        Qa[kt][3] = vb ? pbu[4] : 0u;
+    }
+
+    float Oacc[DN8][4];
+#pragma unroll
+    for (int j = 0; j < DN8; j++)
+        Oacc[j][0] = Oacc[j][1] = Oacc[j][2] = Oacc[j][3] = 0.0f;
+    float m0 = -FLT_MAX, m1 = -FLT_MAX, l0 = 0.0f, l1 = 0.0f;
+
+    const int mask_base = batch * p.kv_len;
+    const int tiles_total = (p.kv_len + BC - 1) / BC;
+    const int tiles_per_split = (tiles_total + p.num_splits - 1) / p.num_splits;
+    const int ti_begin = split * tiles_per_split;
+    const int ti_end = min(tiles_total, ti_begin + tiles_per_split);
+    const int has_mask = p.use_mask && p.mask;
+
+    // Paged strides (constant for the block)
+    const int64_t page_stride = (int64_t)p.page_size * p.kv_head * HEAD_DIM;
+    const int64_t pos_stride  = (int64_t)p.kv_head * HEAD_DIM;
+    const int64_t head_off    = (int64_t)kv_head_idx * HEAD_DIM;
+
+    // ---- Load tile lambda: predicated cp.async, paged addressing ----
+    auto load_tile = [&](int ti, int buf) {
+        int kv0 = ti * BC;
+        bf16* dK = sK + buf * BC * LD;
+        bf16* dV = sV + buf * BC * LD;
+        int logical_page = kv0 / p.page_size;
+        int phys_page = p.page_table[batch * p.max_pages + logical_page];
+        bool page_valid = (phys_page >= 0);
+#pragma unroll
+        for (int i = lane * VEC; i < TOTAL; i += 32 * VEC) {
+            int r = i / HEAD_DIM, d = i % HEAD_DIM;
+            int kc = kv0 + r;
+            bool valid = (kc < p.kv_len) && page_valid;
+            int page_off = kc % p.page_size;
+            int64_t gmem_base = (int64_t)phys_page * page_stride
+                              + (int64_t)page_off * pos_stride
+                              + head_off;
+            int off = r * LD + swiz_col(d, r, SWIZ_MASK);
+            cp_async_16_pred(&dK[off], &p.k_cache[gmem_base + d], valid);
+            cp_async_16_pred(&dV[off], &p.v_cache[gmem_base + d], valid);
+        }
+        cp_async_commit();
+    };
+
+    // ---- Prologue: issue first tile load ----
+    if (ti_begin < ti_end) {
+        load_tile(ti_begin, 0);
+    }
+
+    for (int ti = ti_begin; ti < ti_end; ti++) {
+        constexpr int BUF_MASK = (STAGES > 1) ? (STAGES - 1) : 0;
+        int buf = (ti - ti_begin) & BUF_MASK;
+
+        cp_async_wait_group<0>();
+        __syncwarp();
+        if constexpr (STAGES > 1) {
+            if (ti + 1 < ti_end)
+                load_tile(ti + 1, (ti + 1 - ti_begin) & BUF_MASK);
+        }
+
+        const bf16* bK = sK + buf * BC * LD;
+        const bf16* bV = sV + buf * BC * LD;
+        int kv0 = ti * BC;
+
+        float Sacc[NC8][4];
+        mma_compute_scores<KD, NC8>(Qa, bK, LD, SWIZ_MASK, lane, Sacc);
+
+#pragma unroll
+        for (int n8 = 0; n8 < NC8; n8++)
+            Sacc[n8][0] *= p.scale, Sacc[n8][1] *= p.scale,
+            Sacc[n8][2] *= p.scale, Sacc[n8][3] *= p.scale;
+
+        int maxc = p.is_causal ? min(p.kv_len, p.causal_offset + 1) : p.kv_len;
+        mma_softmax_tile<NC8, DN8>(kv0, maxc, maxc,
+                                    mask_base, p.mask, has_mask,
+                                    Sacc, Oacc, m0, m1, l0, l1, lane);
+
+        mma_pv_accumulate<DN8, KT2>(Sacc, bV, LD, SWIZ_MASK, lane, Oacc);
+        __syncwarp();
+
+        if constexpr (STAGES == 1) {
+            if (ti + 1 < ti_end)
+                load_tile(ti + 1, 0);
+        }
+    }
+
+    // ---- write UN-normalised partials for this split ----
+    auto split_slot = [&](int h) -> size_t {
+        size_t bh = (size_t)batch * p.q_head + h;
+        return bh * p.num_splits + split;
+    };
+#pragma unroll
+    for (int dn8 = 0; dn8 < DN8; dn8++) {
+        int d = dn8 * 8 + 2 * tid4;
+        int r0 = gid, r1 = gid + 8;
+        if (r0 < G) {
+            int h = q_head0 + r0;
+            float* op = p.o_part + split_slot(h) * HEAD_DIM;
+            op[d] = Oacc[dn8][0];
+            op[d + 1] = Oacc[dn8][1];
+        }
+        if (r1 < G) {
+            int h = q_head0 + r1;
+            float* op = p.o_part + split_slot(h) * HEAD_DIM;
+            op[d] = Oacc[dn8][2];
+            op[d + 1] = Oacc[dn8][3];
+        }
+    }
+    if (tid4 == 0) {
+        int r0 = gid, r1 = gid + 8;
+        if (r0 < G) {
+            int h = q_head0 + r0;
+            float* mp = p.ml_part + split_slot(h) * 2;
+            mp[0] = m0; mp[1] = l0;
+        }
+        if (r1 < G) {
+            int h = q_head0 + r1;
+            float* mp = p.ml_part + split_slot(h) * 2;
+            mp[0] = m1; mp[1] = l1;
+        }
+    }
+}

csrc/kernels/attn_prefill.cu

@@ -0,0 +1,79 @@
+#include "attn_prefill_split_q.cuh"
+#include "attn_entry_utils.cuh"
+
+#ifndef ASTRAI_NO_MMA
+#include "attn_prefill_split_q_mma.cuh"
+#endif
+
+template <int HEAD_DIM>
+static void dispatch_prefill(AttentionParams<bf16>& p) {
+#ifndef ASTRAI_NO_MMA
+    constexpr int WARPS = 4, BR = 16;
+    // KV tile: bigger tiles amortize the per-tile cp.async wait + barrier +
+    // loop overhead over more tensor-core work (this kernel is latency-bound,
+    // not compute/bandwidth-bound), so BC=32 wins ~6-8% over BC=16 for
+    // D<=128. D=256 stays at 16: BC=32 double-buffered would need 64KB smem,
+    // over the 48KB static cap.
+    constexpr int BC = (HEAD_DIM <= 128) ? 32 : 16;
+    dim3 grid((p.q_len + BR * WARPS - 1) / (BR * WARPS), p.q_head, p.batch);
+    dim3 block(WARPS * 32, 1, 1);
+    // Static shared memory — double-buffered K/V only (no sQ: Q goes direct
+    // to registers). 2*BC*LD bf16 each for sK and sV → 4*BC*HEAD_DIM*2 bytes.
+    // Occupancy is smem-capped: D=64→3 blocks/SM (16KB), D=128→1 (32KB),
+    // D=256→1 (32KB, BC=16).
+    attn_prefill_split_q_mma_kernel<HEAD_DIM, WARPS, BC><<<grid, block>>>(p);
+#else
+    constexpr int G = 8, ROWS = 32, P_BC = 32;
+    dim3 grid((p.q_len + ROWS - 1) / ROWS, p.q_head, p.batch);
+    dim3 block(G, ROWS, 1);
+    attn_prefill_split_q_kernel_t<HEAD_DIM, G, ROWS, P_BC><<<grid, block>>>(p);
+#endif
+}
+
+torch::Tensor attn_prefill(
+    torch::Tensor q,
+    torch::Tensor k,
+    torch::Tensor v,
+    c10::optional<torch::Tensor> mask,
+    bool is_causal = false,
+    int64_t causal_offset = 0,
+    c10::optional<double> scale = c10::nullopt
+) {
+    AttentionParams<bf16> p;
+    attn_pack_params(q, k, v, mask, is_causal, causal_offset, scale, p);
+    TORCH_CHECK(p.head_dim % 16 == 0, "head_dim must be multiple of 16");
+
+    auto O = torch::empty_like(q);
+    p.o = (bf16*)O.data_ptr();
+
+    switch (p.head_dim) {
+        case 32:
+            dispatch_prefill<32>(p);
+            break;
+        case 64:
+            dispatch_prefill<64>(p);
+            break;
+        case 128:
+            dispatch_prefill<128>(p);
+            break;
+        case 256:
+            dispatch_prefill<256>(p);
+            break;
+        default:
+            TORCH_CHECK(false, "prefill: unsupported head_dim ", p.head_dim,
+                        " (supported: 32,64,128,256)");
+    }
+    return O;
+}
+
+PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
+    m.def("attn_prefill", &attn_prefill,
+        py::arg("q"),
+        py::arg("k"),
+        py::arg("v"),
+        py::arg("mask") = py::none(),
+        py::arg("is_causal") = false,
+        py::arg("causal_offset") = 0,
+        py::arg("scale") = py::none(),
+        "GQA prefill (tensor-core mma on sm_80+, scalar fallback)");
+}

csrc/kernels/attn_prefill_split_q.cuh

@@ -1,5 +1,9 @@
 #pragma once
-#include "gqa_common.cuh"
+#include <cfloat>
+#include <cuda_bf16.h>
+#include "attn_common.h"
+
+using bf16 = __nv_bfloat16;
 
 // v9: group-split register blocking. G threads cooperate on one query row,
 // each owning HEAD_DIM/G dims of qreg[]/acc[]. Small per-thread footprint keeps
@@ -31,7 +35,7 @@ __device__ __forceinline__ void ld8(const bf16* p, float* o) {
 }
 
 template <int HEAD_DIM, int G, int ROWS, int P_BC>
-__global__ void gqa_prefill_attn_kernel_t(GQAParams p) {
+__global__ void attn_prefill_split_q_kernel_t(AttentionParams<bf16> p) {
     constexpr int DPT = HEAD_DIM / G;
 
     int q_tile = blockIdx.x;
@@ -43,9 +47,8 @@ __global__ void gqa_prefill_attn_kernel_t(GQAParams p) {
 
     int kv_head = q_head / (p.q_head / p.kv_head);
 
-    extern __shared__ __align__(16) bf16 smem[];
-    bf16* sK = smem;
-    bf16* sV = sK + P_BC * HEAD_DIM;
+    __shared__ __align__(16) bf16 sK[P_BC * HEAD_DIM];
+    __shared__ __align__(16) bf16 sV[P_BC * HEAD_DIM];
 
     float qreg[DPT];
     if (q_row < p.q_len) {

csrc/kernels/attn_prefill_split_q_mma.cuh

@@ -0,0 +1,204 @@
+#pragma once
+#include <cfloat>
+#include <cuda_bf16.h>
+#include "attn_common.h"
+#include "attn_mma_utils.cuh"
+
+using bf16 = __nv_bfloat16;
+
+// Tensor-core prefill flash attention (raw mma.sync PTX).
+// One warp owns BR=16 query rows. S = Q@K^T and O = P@V run on bf16 tensor
+// cores via mma.sync.m16n8k16 (f32 accumulate). Q fragments are loaded once
+// straight from global into the mma A-operand layout (no smem staging) and
+// kept resident in registers across the tile loop. S, O, and the online-softmax
+// stats (m, l) also live in registers.
+// Shared memory is statically sized via template parameters — no dynamic
+// allocation. The mma fragment layout is used directly: the S accumulator
+// (f32) maps element-for-element onto the P matrix_a (bf16) operand, so
+// softmax needs no shuffle repack; row reductions fold across the 4-lane
+// thread group. Templated on <HEAD_DIM, WARPS, BC> with BC a multiple of 16.
+//
+// Software pipeline: K/V are double-buffered and loaded via cp.async one tile
+// ahead, so the next tile streams from global memory while the current tile's
+// tensor-core math runs — hiding load latency (long_scoreboard). A single
+// __syncthreads per tile both publishes the freshly loaded tile cross-warp and
+// (because it runs before the next prefetch) guards the buffer being refilled,
+// so no second barrier is needed. Predicated cp.async (cp_async_16_pred)
+// zero-fills rows past kv_len, unifying full and partial tiles on one path.
+// BC=32 (D<=128) amortizes the per-tile wait+barrier+loop overhead over more
+// tensor-core work — this kernel is latency-bound (low occupancy from high
+// register pressure), so fewer, larger tiles beat many tiny ones.
+//
+// Optimizations: load Q fragments directly from global in mma A-operand layout
+// (no sQ staging, no prologue barriers); post-multiply scale in float after
+// S=Q@K^T to avoid bf16 precision loss; packed bf16x2 output stores;
+// causal tile skipping (block-level prefetch bound + warp-level compute skip);
+// XOR swizzle (swiz_col) → eliminates ldmatrix bank conflicts without LD
+// padding (LD=HEAD_DIM).
+
+template <int HEAD_DIM, int WARPS, int BC>
+__global__ void attn_prefill_split_q_mma_kernel(AttentionParams<bf16> p) {
+    constexpr int BR = 16;
+    constexpr int KD = HEAD_DIM / 16;  // Q/K k-tiles
+    constexpr int NC8 = BC / 8;        // S n-tiles (N=8 each)
+    constexpr int KT2 = BC / 16;       // P k-tiles (K=16 each)
+    constexpr int DN8 = HEAD_DIM / 8;  // O n-tiles (N=8 each)
+    constexpr int LD = HEAD_DIM;   // XOR swizzle (swiz_col) handles bank conflicts
+    constexpr int SWIZ_MASK = (HEAD_DIM >= 64) ? 7 : (HEAD_DIM / 8 - 1);  // chunk bits, stay within LD
+
+    const int warp = threadIdx.x / 32;
+    const int lane = threadIdx.x % 32;
+    const int gid = lane >> 2;   // 0..7  → rows gid, gid+8
+    const int tid4 = lane & 3;   // 0..3
+    const int nthreads = WARPS * 32;
+
+    const int q_head = blockIdx.y;
+    const int batch = blockIdx.z;
+    const int kv_head = q_head / (p.q_head / p.kv_head);
+    const int qrow0 = (blockIdx.x * WARPS + warp) * BR;
+
+    // Static shared memory — sized by template parameters at compile time.
+    // K/V are double-buffered (STAGES=2): the next tile's cp.async load runs
+    // while the current tile's tensor-core math executes, hiding global-load
+    // latency (FA2-style software pipeline). No dynamic smem / carveout opt-in.
+    constexpr int STAGES = 2;
+    __shared__ __align__(16) bf16 sK[STAGES * BC * LD];
+    __shared__ __align__(16) bf16 sV[STAGES * BC * LD];
+
+    // Load the Q fragments straight from global into the mma A-operand layout
+    // (m16n8k16, row-major): no sQ staging area and no serialized per-warp
+    // prologue barriers. Each lane reads exactly the 8 Q elements ldmatrix
+    // would have produced, pre-scaled by the attention scale. Kept resident in
+    // registers across the tile loop.
+    //   frag[0]/[2]: row = qrow0 + gid ;  frag[1]/[3]: row = qrow0 + gid + 8
+    //   frag[0]/[1]: cols kt*16 + tid4*2 + {0,1} ; frag[2]/[3]: + 8
+    const int q_base = ((batch * p.q_head + q_head) * p.q_len) * HEAD_DIM;
+    const int qra = qrow0 + gid;
+    const int qrb = qrow0 + gid + 8;
+    const bool va = qra < p.q_len, vb = qrb < p.q_len;
+    unsigned Qa[KD][4];
+#pragma unroll
+    for (int kt = 0; kt < KD; kt++) {
+        int c = kt * 16 + tid4 * 2;
+        const unsigned* pau = reinterpret_cast<const unsigned*>(
+            &p.q[q_base + qra * HEAD_DIM + c]);
+        const unsigned* pbu = reinterpret_cast<const unsigned*>(
+            &p.q[q_base + qrb * HEAD_DIM + c]);
+        Qa[kt][0] = va ? pau[0] : 0u;
+        Qa[kt][1] = vb ? pbu[0] : 0u;
+        Qa[kt][2] = va ? pau[4] : 0u;
+        Qa[kt][3] = vb ? pbu[4] : 0u;
+    }
+
+    float Oacc[DN8][4];
+#pragma unroll
+    for (int j = 0; j < DN8; j++)
+        Oacc[j][0] = Oacc[j][1] = Oacc[j][2] = Oacc[j][3] = 0.0f;
+    float m0 = -FLT_MAX, m1 = -FLT_MAX, l0 = 0.0f, l1 = 0.0f;
+
+    const int kv_base = ((batch * p.kv_head + kv_head) * p.kv_len) * HEAD_DIM;
+    const int tiles = (p.kv_len + BC - 1) / BC;
+    const int qr0 = qrow0 + gid;      // row for c0/c1
+    const int qr1 = qrow0 + gid + 8;  // row for c2/c3
+
+    // Causal tile-skip bounds (no-op when is_causal == 0)
+    const int use_skip = p.is_causal;
+    const int max_kv = qrow0 + BR - 1 + p.causal_offset;
+    const int block_max_kv =
+        blockIdx.x * WARPS * BR + WARPS * BR - 1 + p.causal_offset;
+    const int has_mask = p.use_mask && p.mask;
+    const int mb = batch * p.kv_len;
+
+    // Last active tile: block-level causal bound (all warps in the block share
+    // the K/V load, so the prefetch range is the block max, not per-warp).
+    int t_end = tiles - 1;
+    if (use_skip) {
+        int bt = block_max_kv / BC;
+        if (bt < t_end) t_end = bt;
+    }
+
+    constexpr int VEC = 8;  // bf16 per cp.async unit (16 bytes)
+    constexpr int TOTAL = BC * HEAD_DIM;
+
+    // Issue cp.async loads for tile `ti` into shared buffer `buf`. Predicated
+    // loads zero-fill rows past kv_len, so partial tiles need no scalar path.
+    auto load_tile = [&](int ti, int buf) {
+        int kv0 = ti * BC;
+        bf16* dK = sK + buf * BC * LD;
+        bf16* dV = sV + buf * BC * LD;
+#pragma unroll
+        for (int i = threadIdx.x * VEC; i < TOTAL; i += nthreads * VEC) {
+            int r = i / HEAD_DIM, d = i % HEAD_DIM;
+            int kc = kv0 + r;
+            bool valid = kc < p.kv_len;
+            int off = r * LD + swiz_col(d, r, SWIZ_MASK);
+            cp_async_16_pred(&dK[off], &p.k[kv_base + kc * HEAD_DIM + d], valid);
+            cp_async_16_pred(&dV[off], &p.v[kv_base + kc * HEAD_DIM + d], valid);
+        }
+        cp_async_commit();
+    };
+
+    // Prologue: kick off the first tile's load.
+    load_tile(0, 0);
+
+    for (int ti = 0; ti <= t_end; ti++) {
+        int buf = ti & 1;
+
+        // Wait for the current tile's async copies, then a single barrier: it
+        // both publishes this tile's data cross-warp AND guarantees the prior
+        // compute on the buffer we are about to refill has finished. Issuing
+        // the next tile's load *after* this barrier lets one barrier cover both
+        // hazards (vs two), while the load still overlaps this tile's math.
+        cp_async_wait_group<0>();
+        __syncthreads();
+        if (ti < t_end) load_tile(ti + 1, (ti + 1) & 1);
+
+        const bf16* bK = sK + buf * BC * LD;
+        const bf16* bV = sV + buf * BC * LD;
+        int kv0 = ti * BC;
+
+        // Warp-level causal skip
+        if (!use_skip || kv0 <= max_kv) {
+
+        // S = Q @ K^T + scale + online softmax + O += P @ V
+        float Sacc[NC8][4];
+        mma_compute_scores<KD, NC8>(Qa, bK, LD, SWIZ_MASK, lane, Sacc);
+
+        // post-multiply scale in float (no bf16 precision loss from pre-scaling Q)
+        #pragma unroll
+        for (int n8 = 0; n8 < NC8; n8++)
+            Sacc[n8][0] *= p.scale, Sacc[n8][1] *= p.scale,
+            Sacc[n8][2] *= p.scale, Sacc[n8][3] *= p.scale;
+
+        int maxc0 = p.is_causal ? min(p.kv_len, qr0 + p.causal_offset + 1)
+                                : p.kv_len;
+        int maxc1 = p.is_causal ? min(p.kv_len, qr1 + p.causal_offset + 1)
+                                : p.kv_len;
+        mma_softmax_tile<NC8, DN8>(kv0, maxc0, maxc1,
+                                    mb, p.mask, has_mask,
+                                    Sacc, Oacc, m0, m1, l0, l1, lane);
+
+        mma_pv_accumulate<DN8, KT2>(Sacc, bV, LD, SWIZ_MASK, lane, Oacc);
+        }  // if active (warp-level causal skip)
+    }
+
+    // ---- write output ---- (packed bf16x2 stores: one 32-bit STG per pair,
+    // halves store count and removes the uncoalesced scalar-store penalty)
+    float rl0 = (l0 > 1e-20f) ? (1.0f / l0) : 0.0f;
+    float rl1 = (l1 > 1e-20f) ? (1.0f / l1) : 0.0f;
+    const int o_base = ((batch * p.q_head + q_head) * p.q_len) * HEAD_DIM;
+#pragma unroll
+    for (int dn8 = 0; dn8 < DN8; dn8++) {
+        int d = dn8 * 8 + 2 * tid4;
+        if (qr0 < p.q_len) {
+            __nv_bfloat162 v = __floats2bfloat162_rn(Oacc[dn8][0] * rl0,
+                                                     Oacc[dn8][1] * rl0);
+            *reinterpret_cast<__nv_bfloat162*>(&p.o[o_base + qr0 * HEAD_DIM + d]) = v;
+        }
+        if (qr1 < p.q_len) {
+            __nv_bfloat162 v = __floats2bfloat162_rn(Oacc[dn8][2] * rl1,
+                                                     Oacc[dn8][3] * rl1);
+            *reinterpret_cast<__nv_bfloat162*>(&p.o[o_base + qr1 * HEAD_DIM + d]) = v;
+        }
+    }
+}

csrc/kernels/gqa_common.cuh

@@ -1,35 +0,0 @@
-#pragma once
-#include <cuda_bf16.h>
-#include <cuda_runtime.h>
-#include <cfloat>
-#include <algorithm>
-
-using bf16 = __nv_bfloat16;
-using std::min;
-
-constexpr int DC_CHUNK = 64;
-constexpr int Br = 32, Bc = 64;
-
-__device__ inline float warp_reduce_sum(float val) {
-    for (int offset = 16; offset > 0; offset >>= 1)
-        val += __shfl_xor_sync(0xFFFFFFFF, val, offset);
-    return val;
-}
-
-struct GQAParams {
-    int batch;
-    int q_head;
-    int kv_head;
-    int q_len;
-    int kv_len;
-    int head_dim;
-    int use_mask;
-    int is_causal;
-    int causal_offset;
-    float scale;
-    const bf16* __restrict__ q;
-    const bf16* __restrict__ k;
-    const bf16* __restrict__ v;
-    const bool* __restrict__ mask;
-    bf16* __restrict__ o;
-};

csrc/kernels/gqa_decode_attn.cu

@@ -1,114 +0,0 @@
-#include "gqa_decode_attn.cuh"
-#include <torch/extension.h>
-
-#ifndef ASTRAI_NO_MMA
-#include "gqa_decode_attn_mma.cuh"
-#endif
-
-template <int HEAD_DIM>
-static void dispatch_decode(GQAParams& p) {
-#ifndef ASTRAI_NO_MMA
-    constexpr int BC = 32, BR = 16, LD = HEAD_DIM;  // XOR swizzle → no padding
-    int G = p.q_head / p.kv_head;
-    // head-packing tensor-core path needs 1 < G <= 16 (MMA M dim) and no mask;
-    // everything else uses the scalar kernel
-    if (!p.use_mask && G > 1 && G <= 16) {
-        dim3 grid(p.kv_head, p.batch, 1);
-        dim3 block(32, 1, 1);
-        // sK + sV + sQ, each BC/BR * LD (single buffer for high occupancy)
-        int smem = (2 * BC * LD + BR * LD) * (int)sizeof(bf16);
-        cudaFuncSetAttribute(gqa_decode_attn_mma_kernel<HEAD_DIM, BC>,
-                             cudaFuncAttributeMaxDynamicSharedMemorySize, smem);
-        gqa_decode_attn_mma_kernel<HEAD_DIM, BC><<<grid, block, smem>>>(p);
-        return;
-    }
-    // scalar fallback (per-KV-head, one warp per query head)
-    int group_size = p.q_head / p.kv_head;
-    size_t smem = DC_CHUNK * p.head_dim * sizeof(bf16);
-    dim3 block(32, group_size);
-    dim3 grid(p.batch * p.kv_head);
-    gqa_decode_attn_kernel<<<grid, block, smem>>>(p);
-#else
-    // scalar fallback (per-KV-head, one warp per query head)
-    int group_size = p.q_head / p.kv_head;
-    size_t smem = DC_CHUNK * p.head_dim * sizeof(bf16);
-    dim3 block(32, group_size);
-    dim3 grid(p.batch * p.kv_head);
-    gqa_decode_attn_kernel<<<grid, block, smem>>>(p);
-#endif
-}
-
-torch::Tensor gqa_decode_attn(
-    torch::Tensor q,
-    torch::Tensor k, 
-    torch::Tensor v,
-    c10::optional<torch::Tensor> mask,
-    bool is_causal = false, 
-    int64_t causal_offset = 0,
-    c10::optional<double> scale = c10::nullopt
-) {
-    TORCH_CHECK(q.is_cuda() && k.is_cuda() && v.is_cuda());
-    TORCH_CHECK(q.dtype() == torch::kBFloat16);
-    TORCH_CHECK(k.dtype() == torch::kBFloat16);
-    TORCH_CHECK(v.dtype() == torch::kBFloat16);
-    TORCH_CHECK(q.size(2) == 1, "Q seq_len must be 1");
-
-    GQAParams p;
-    p.batch = q.size(0);
-    p.q_head = q.size(1);
-    p.kv_head = k.size(1);
-    p.q_len = 1;
-    p.kv_len = k.size(2);
-    p.head_dim = q.size(3);
-    TORCH_CHECK(p.head_dim % 32 == 0, "head_dim must be multiple of 32");
-    p.use_mask = mask.has_value();
-    p.is_causal = (int)is_causal;
-    p.causal_offset = (int)causal_offset;
-    p.scale = scale.has_value() ? (float)scale.value() : 1.0f / sqrtf((float)p.head_dim);
-    p.q = (const bf16*)q.data_ptr();
-    p.k = (const bf16*)k.data_ptr();
-    p.v = (const bf16*)v.data_ptr();
-    if (p.use_mask) {
-        TORCH_CHECK(mask.value().dtype() == torch::kBool);
-        TORCH_CHECK(mask.value().dim() == 2);
-        TORCH_CHECK(mask.value().size(0) == p.batch);
-        TORCH_CHECK(mask.value().size(1) == p.kv_len);
-        p.mask = mask.value().data_ptr<bool>();
-    } else {
-        p.mask = nullptr;
-    }
-
-    auto O = torch::empty_like(q);
-    p.o = (bf16*)O.data_ptr();
-
-    switch (p.head_dim) {
-        case 32:
-            dispatch_decode<32>(p);
-            break;
-        case 64:
-            dispatch_decode<64>(p);
-            break;
-        case 128:
-            dispatch_decode<128>(p);
-            break;
-        case 256:
-            dispatch_decode<256>(p);
-            break;
-        default:
-            TORCH_CHECK(false, "decode: unsupported head_dim ", p.head_dim,
-                        " (supported: 32, 64, 128, 256)");
-    }
-    return O;
-}
-
-PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
-    m.def("gqa_decode_attn", &gqa_decode_attn,
-        py::arg("q"),
-        py::arg("k"),
-        py::arg("v"),
-        py::arg("mask") = py::none(),
-        py::arg("is_causal") = false,
-        py::arg("causal_offset") = 0,
-        py::arg("scale") = py::none(),
-        "GQA decode (tensor-core head-packing on sm_80+, scalar fallback)");
-}

csrc/kernels/gqa_decode_attn.cuh

@@ -1,59 +0,0 @@
-#pragma once
-#include "gqa_common.cuh"
-
-__global__ void gqa_decode_attn_kernel(GQAParams p) {
-    int batch = blockIdx.x / p.kv_head;
-    int kv_head = blockIdx.x % p.kv_head;
-    int group_size = blockDim.y;
-    int q_head = kv_head * group_size + threadIdx.y;
-    int lane = threadIdx.x;
-    int hd_per_thread = p.head_dim / 32;
-
-    float q_reg[8];
-    int q_off = ((batch * p.q_head + q_head) * 1) * p.head_dim + lane * hd_per_thread;
-    for (int i = 0; i < hd_per_thread; i++)
-        q_reg[i] = __bfloat162float(p.q[q_off + i]);
-
-    int kv_base = ((batch * p.kv_head + kv_head) * p.kv_len) * p.head_dim;
-    int mask_base = batch * p.kv_len;
-
-    float m = -FLT_MAX, d = 0.0f, acc_reg[8] = {0.0f};
-
-    extern __shared__ __align__(16) bf16 k_smem[];
-
-    for (int chunk_start = 0; chunk_start < p.kv_len; chunk_start += DC_CHUNK) {
-        int this_chunk = min(DC_CHUNK, p.kv_len - chunk_start);
-
-        int total = this_chunk * p.head_dim;
-        for (int i = threadIdx.y * 32 + lane; i < total; i += blockDim.x * blockDim.y)
-            k_smem[i] = p.k[kv_base + chunk_start * p.head_dim + i];
-        __syncthreads();
-
-        for (int s = 0; s < this_chunk; s++) {
-            float partial = 0.0f;
-            for (int i = 0; i < hd_per_thread; i++)
-                partial += q_reg[i] * __bfloat162float(k_smem[s * p.head_dim + lane * hd_per_thread + i]);
-            partial = warp_reduce_sum(partial) * p.scale;
-
-            if (p.use_mask && p.mask && !p.mask[mask_base + chunk_start + s])
-                partial = -FLT_MAX;
-            if (p.is_causal && (chunk_start + s) > p.causal_offset)
-                partial = -FLT_MAX;
-
-            float new_m = fmaxf(m, partial);
-            float alpha = expf(m - new_m);
-            float beta  = expf(partial - new_m);
-            d = d * alpha + beta;
-
-            int v_off = kv_base + (chunk_start + s) * p.head_dim + lane * hd_per_thread;
-            for (int i = 0; i < hd_per_thread; i++)
-                acc_reg[i] = acc_reg[i] * alpha + __bfloat162float(p.v[v_off + i]) * beta;
-            m = new_m;
-        }
-        __syncthreads();
-    }
-
-    int out_off = ((batch * p.q_head + q_head) * 1) * p.head_dim + lane * hd_per_thread;
-    for (int i = 0; i < hd_per_thread; i++)
-        p.o[out_off + i] = __float2bfloat16(acc_reg[i] / d);
-}

csrc/kernels/gqa_decode_attn_mma.cuh

@@ -1,219 +0,0 @@
-#pragma once
-#include "gqa_common.cuh"
-#include "gqa_mma_utils.cuh"
-
-// Tensor-core decode via GQA head-packing with cp.async loads.
-//
-// Decode has q_len == 1, so S = q @ K^T is a GEMV per head — no tensor-core work
-// on its own. But GQA gives us G = q_head / kv_head query heads that all share
-// one kv_head. We pack those G heads into the M=16 rows of mma.sync.m16n8k16,
-// turning G independent GEMVs into a single GEMM that reuses each loaded K/V tile
-// across all G heads (K/V load is the decode bottleneck, so the reuse is the win,
-// not the flops). Fragment layout is identical to the prefill mma kernel; the
-// only differences are (1) the M rows come from different heads at position 0
-// instead of different sequence positions of one head, and (2) causal masking is
-// a single scalar bound shared by every row. One warp owns one (batch, kv_head);
-// requires G <= 16.
-//
-// Optimizations:
-//   - cp.async global→shared for K/V (bypasses registers, cuts instruction count)
-//   - XOR swizzle (swiz_col): LD=HEAD_DIM, zero waste, no bank conflicts
-//   - pre-scaled Q: Q scaled during load, softmax skips per-tile multiply
-//   - single-buffer: keeps smem small for high occupancy
-
-template <int HEAD_DIM, int BC>
-__global__ void gqa_decode_attn_mma_kernel(GQAParams p) {
-    constexpr int BR = 16;
-    constexpr int KD = HEAD_DIM / 16;  // Q/K k-tiles
-    constexpr int NC8 = BC / 8;        // S n-tiles (N=8 each)
-    constexpr int KT2 = BC / 16;       // P k-tiles (K=16 each)
-    constexpr int DN8 = HEAD_DIM / 8;  // O n-tiles (N=8 each)
-    constexpr int LD = HEAD_DIM;       // XOR swizzle handles bank conflicts, zero waste
-    constexpr int SWIZ_MASK = (HEAD_DIM >= 64) ? 7 : (HEAD_DIM / 8 - 1);
-
-    const int lane = threadIdx.x;  // single warp
-    const int gid = lane >> 2;     // 0..7 → rows gid, gid+8
-    const int tid4 = lane & 3;
-
-    const int kv_head = blockIdx.x;
-    const int batch = blockIdx.y;
-    const int G = p.q_head / p.kv_head;
-    const int q_head0 = kv_head * G;
-
-    extern __shared__ __align__(16) bf16 smem[];
-    bf16* sK = smem;                  // [BC][LD]
-    bf16* sV = sK + BC * LD;          // [BC][LD]
-    bf16* sQ = sV + BC * LD;          // [BR][LD]
-
-    // ---- stage Q into shared (pre-scaled, swizzled) ----
-    bf16 scale_bf16 = __float2bfloat16(p.scale);
-    for (int i = lane; i < BR * HEAD_DIM; i += 32) {
-        int r = i / HEAD_DIM, d = i % HEAD_DIM;
-        bf16 val = __float2bfloat16(0.0f);
-        if (r < G) {
-            int qh = q_head0 + r;
-            val = p.q[(batch * p.q_head + qh) * HEAD_DIM + d];  // q_len == 1
-        }
-        sQ[r * LD + swiz_col(d, r, SWIZ_MASK)] = __hmul(val, scale_bf16);
-    }
-    __syncwarp();
-
-    // Q resident A-fragments
-    unsigned Qa[KD][4];
-    int qrow_l = (lane & 7) + (lane & 8);
-    int qcol_l = (lane & 16) ? 8 : 0;
-#pragma unroll
-    for (int kt = 0; kt < KD; kt++)
-        ldmatrix_x4(Qa[kt], &sQ[qrow_l * LD + swiz_col(kt * 16 + qcol_l, qrow_l, SWIZ_MASK)]);
-
-    float Oacc[DN8][4];
-#pragma unroll
-    for (int j = 0; j < DN8; j++)
-        Oacc[j][0] = Oacc[j][1] = Oacc[j][2] = Oacc[j][3] = 0.0f;
-    float m0 = -FLT_MAX, m1 = -FLT_MAX, l0 = 0.0f, l1 = 0.0f;
-
-    const int kv_base = (batch * p.kv_head + kv_head) * p.kv_len * HEAD_DIM;
-    const int mask_base = batch * p.kv_len;
-    const int tiles = (p.kv_len + BC - 1) / BC;
-    const int has_mask = p.use_mask && p.mask;
-
-    for (int ti = 0; ti < tiles; ti++) {
-        int kv0 = ti * BC;
-
-        // ---- load K/V tile to shared (cp.async on full tiles) ----
-        bool full_tile = (kv0 + BC <= p.kv_len);
-        if (full_tile) {
-            constexpr int VEC = 8;  // 8 bf16 = 16 bytes per cp.async
-            int total = BC * HEAD_DIM;
-#pragma unroll
-            for (int i = lane * VEC; i < total; i += 32 * VEC) {
-                int r = i / HEAD_DIM, d = i % HEAD_DIM;
-                int kc = kv0 + r;
-                cp_async_16(&sK[r * LD + swiz_col(d, r, SWIZ_MASK)],
-                            &p.k[kv_base + kc * HEAD_DIM + d]);
-                cp_async_16(&sV[r * LD + swiz_col(d, r, SWIZ_MASK)],
-                            &p.v[kv_base + kc * HEAD_DIM + d]);
-            }
-            cp_async_commit();
-            cp_async_wait_all();
-        } else {
-            for (int i = lane; i < BC * HEAD_DIM; i += 32) {
-                int r = i / HEAD_DIM, d = i % HEAD_DIM;
-                int kc = kv0 + r;
-                bf16 z = __float2bfloat16(0.0f);
-                sK[r * LD + swiz_col(d, r, SWIZ_MASK)] =
-                    (kc < p.kv_len) ? p.k[kv_base + kc * HEAD_DIM + d] : z;
-                sV[r * LD + swiz_col(d, r, SWIZ_MASK)] =
-                    (kc < p.kv_len) ? p.v[kv_base + kc * HEAD_DIM + d] : z;
-            }
-        }
-        __syncwarp();
-
-        // S = Q @ K^T  (Q already pre-scaled, so Sacc includes scale)
-        float Sacc[NC8][4];
-#pragma unroll
-        for (int n8 = 0; n8 < NC8; n8++) {
-            Sacc[n8][0] = Sacc[n8][1] = Sacc[n8][2] = Sacc[n8][3] = 0.0f;
-            int krow_l = n8 * 8 + (lane & 7);
-            int kcol_h = (lane & 8) ? 8 : 0;
-#pragma unroll
-            for (int kt = 0; kt < KD; kt++) {
-                unsigned b[2];
-                ldmatrix_x2(b, &sK[krow_l * LD + swiz_col(kt * 16 + kcol_h, krow_l, SWIZ_MASK)]);
-                mma16816(Sacc[n8], Qa[kt], b, Sacc[n8]);
-            }
-        }
-
-        // ---- online softmax (Q pre-scaled → no per-tile scale multiply) ----
-        float rmax0 = -FLT_MAX, rmax1 = -FLT_MAX;
-#pragma unroll
-        for (int n8 = 0; n8 < NC8; n8++) {
-            int cc = kv0 + n8 * 8 + 2 * tid4;
-            bool bc0 = (cc >= p.kv_len) ||
-                       (has_mask && !p.mask[mask_base + cc]);
-            bool bc1 = (cc + 1 >= p.kv_len) ||
-                       (has_mask && !p.mask[mask_base + cc + 1]);
-            bool cz = p.is_causal;
-            int off = p.causal_offset;
-            bool bad0 = bc0 || (cz && cc > off);
-            bool bad1 = bc1 || (cz && (cc + 1) > off);
-            float s0 = bad0 ? -FLT_MAX : Sacc[n8][0];
-            float s1 = bad1 ? -FLT_MAX : Sacc[n8][1];
-            float s2 = bad0 ? -FLT_MAX : Sacc[n8][2];
-            float s3 = bad1 ? -FLT_MAX : Sacc[n8][3];
-            Sacc[n8][0] = s0; Sacc[n8][1] = s1; Sacc[n8][2] = s2; Sacc[n8][3] = s3;
-            rmax0 = fmaxf(rmax0, fmaxf(s0, s1));
-            rmax1 = fmaxf(rmax1, fmaxf(s2, s3));
-        }
-        rmax0 = fmaxf(rmax0, __shfl_xor_sync(0xFFFFFFFF, rmax0, 1));
-        rmax0 = fmaxf(rmax0, __shfl_xor_sync(0xFFFFFFFF, rmax0, 2));
-        rmax1 = fmaxf(rmax1, __shfl_xor_sync(0xFFFFFFFF, rmax1, 1));
-        rmax1 = fmaxf(rmax1, __shfl_xor_sync(0xFFFFFFFF, rmax1, 2));
-
-        float nm0 = fmaxf(m0, rmax0), nm1 = fmaxf(m1, rmax1);
-        float corr0 = (nm0 == -FLT_MAX) ? 1.0f : __expf(m0 - nm0);
-        float corr1 = (nm1 == -FLT_MAX) ? 1.0f : __expf(m1 - nm1);
-
-        float rsum0 = 0.0f, rsum1 = 0.0f;
-#pragma unroll
-        for (int n8 = 0; n8 < NC8; n8++) {
-            float p0 = (Sacc[n8][0] == -FLT_MAX) ? 0.0f : __expf(Sacc[n8][0] - nm0);
-            float p1 = (Sacc[n8][1] == -FLT_MAX) ? 0.0f : __expf(Sacc[n8][1] - nm0);
-            float p2 = (Sacc[n8][2] == -FLT_MAX) ? 0.0f : __expf(Sacc[n8][2] - nm1);
-            float p3 = (Sacc[n8][3] == -FLT_MAX) ? 0.0f : __expf(Sacc[n8][3] - nm1);
-            Sacc[n8][0] = p0; Sacc[n8][1] = p1; Sacc[n8][2] = p2; Sacc[n8][3] = p3;
-            rsum0 += p0 + p1;
-            rsum1 += p2 + p3;
-        }
-        rsum0 += __shfl_xor_sync(0xFFFFFFFF, rsum0, 1);
-        rsum0 += __shfl_xor_sync(0xFFFFFFFF, rsum0, 2);
-        rsum1 += __shfl_xor_sync(0xFFFFFFFF, rsum1, 1);
-        rsum1 += __shfl_xor_sync(0xFFFFFFFF, rsum1, 2);
-        l0 = l0 * corr0 + rsum0;
-        l1 = l1 * corr1 + rsum1;
-        m0 = nm0; m1 = nm1;
-
-#pragma unroll
-        for (int j = 0; j < DN8; j++) {
-            Oacc[j][0] *= corr0; Oacc[j][1] *= corr0;
-            Oacc[j][2] *= corr1; Oacc[j][3] *= corr1;
-        }
-
-        // O += P @ V
-#pragma unroll
-        for (int kt2 = 0; kt2 < KT2; kt2++) {
-            unsigned Pa[4];
-            Pa[0] = pk2(Sacc[kt2 * 2][0], Sacc[kt2 * 2][1]);
-            Pa[1] = pk2(Sacc[kt2 * 2][2], Sacc[kt2 * 2][3]);
-            Pa[2] = pk2(Sacc[kt2 * 2 + 1][0], Sacc[kt2 * 2 + 1][1]);
-            Pa[3] = pk2(Sacc[kt2 * 2 + 1][2], Sacc[kt2 * 2 + 1][3]);
-            int vrow_l = kt2 * 16 + (lane & 15);
-#pragma unroll
-            for (int dn8 = 0; dn8 < DN8; dn8++) {
-                unsigned b[2];
-                ldmatrix_x2_trans(b, &sV[vrow_l * LD + swiz_col(dn8 * 8, vrow_l, SWIZ_MASK)]);
-                mma16816(Oacc[dn8], Pa, b, Oacc[dn8]);
-            }
-        }
-        __syncwarp();  // sK/sV reused next tile
-    }
-
-    // ---- write output ----
-    float rl0 = (l0 > 1e-20f) ? (1.0f / l0) : 0.0f;
-    float rl1 = (l1 > 1e-20f) ? (1.0f / l1) : 0.0f;
-#pragma unroll
-    for (int dn8 = 0; dn8 < DN8; dn8++) {
-        int d = dn8 * 8 + 2 * tid4;
-        int r0 = gid, r1 = gid + 8;
-        if (r0 < G) {
-            int o_off = (batch * p.q_head + q_head0 + r0) * HEAD_DIM + d;
-            p.o[o_off] = __float2bfloat16(Oacc[dn8][0] * rl0);
-            p.o[o_off + 1] = __float2bfloat16(Oacc[dn8][1] * rl0);
-        }
-        if (r1 < G) {
-            int o_off = (batch * p.q_head + q_head0 + r1) * HEAD_DIM + d;
-            p.o[o_off] = __float2bfloat16(Oacc[dn8][2] * rl1);
-            p.o[o_off + 1] = __float2bfloat16(Oacc[dn8][3] * rl1);
-        }
-    }
-}

csrc/kernels/gqa_mma_utils.cuh

@@ -1,93 +0,0 @@
-#pragma once
-
-// Shared MMA utilities for tensor-core GQA kernels.
-// mma.sync.m16n8k16 PTX wrappers, ldmatrix helpers, and bf16 packing.
-
-// mma.sync.aligned.m16n8k16.row.col.f32.bf16.bf16.f32
-__device__ __forceinline__ void mma16816(float* d, const unsigned* a,
-                                          const unsigned* b, const float* c) {
-    asm volatile(
-        "mma.sync.aligned.m16n8k16.row.col.f32.bf16.bf16.f32 "
-        "{%0,%1,%2,%3}, {%4,%5,%6,%7}, {%8,%9}, {%10,%11,%12,%13};"
-        : "=f"(d[0]), "=f"(d[1]), "=f"(d[2]), "=f"(d[3])
-        : "r"(a[0]), "r"(a[1]), "r"(a[2]), "r"(a[3]), "r"(b[0]), "r"(b[1]),
-          "f"(c[0]), "f"(c[1]), "f"(c[2]), "f"(c[3]));
-}
-
-// read two adjacent bf16 from smem as one packed .b32 (elem0 low, elem1 high)
-__device__ __forceinline__ unsigned ld2(const bf16* p) {
-    return *reinterpret_cast<const unsigned*>(p);
-}
-
-// pack two floats into one bf16x2 as .b32
-__device__ __forceinline__ unsigned pk2(float a, float b) {
-    __nv_bfloat162 v = __floats2bfloat162_rn(a, b);
-    return *reinterpret_cast<unsigned*>(&v);
-}
-
-// pack two (non-contiguous) bf16 into one .b32
-__device__ __forceinline__ unsigned pkb(bf16 a, bf16 b) {
-    __nv_bfloat162 v;
-    v.x = a;
-    v.y = b;
-    return *reinterpret_cast<unsigned*>(&v);
-}
-
-// ldmatrix: cooperatively load mma fragments from smem (one instruction per
-// 16x16 / 16x8 tile) with the exact register layout mma expects — replaces the
-// scalar per-thread fragment packing, cutting shared-load instructions and bank
-// conflicts. Each lane supplies the shared address of one 8-wide row.
-__device__ __forceinline__ void ldmatrix_x4(unsigned* r, const bf16* p) {
-    unsigned a = __cvta_generic_to_shared(p);
-    asm volatile("ldmatrix.sync.aligned.m8n8.x4.shared.b16 {%0,%1,%2,%3}, [%4];"
-                 : "=r"(r[0]), "=r"(r[1]), "=r"(r[2]), "=r"(r[3])
-                 : "r"(a));
-}
-__device__ __forceinline__ void ldmatrix_x2(unsigned* r, const bf16* p) {
-    unsigned a = __cvta_generic_to_shared(p);
-    asm volatile("ldmatrix.sync.aligned.m8n8.x2.shared.b16 {%0,%1}, [%2];"
-                 : "=r"(r[0]), "=r"(r[1])
-                 : "r"(a));
-}
-__device__ __forceinline__ void ldmatrix_x2_trans(unsigned* r, const bf16* p) {
-    unsigned a = __cvta_generic_to_shared(p);
-    asm volatile("ldmatrix.sync.aligned.m8n8.x2.trans.shared.b16 {%0,%1}, [%2];"
-                 : "=r"(r[0]), "=r"(r[1])
-                 : "r"(a));
-}
-
-// XOR swizzle for shared-memory column at 8-bf16 chunk granularity.
-// Eliminates ldmatrix bank conflicts without LD padding: consecutive rows
-// land in distinct bank groups. swiz_col(d, r, mask) = ((d>>3)^(r&mask))<<3 | (d&7).
-// mask must cover log2(HEAD_DIM/8) chunk bits but stay within LD: use 7 for
-// HEAD_DIM>=64 (8+ chunks), 3 for HEAD_DIM=32 (4 chunks). Default 7 keeps
-// existing HEAD_DIM>=64 call sites working unchanged.
-__device__ __forceinline__ int swiz_col(int d, int r, int mask = 7) {
-    return ((d >> 3) ^ (r & mask)) << 3 | (d & 7);
-}
-
-// cp.async: copy 16 bytes (8 bf16) from global to shared memory directly,
-// bypassing registers. Eliminates shared-store bank conflicts and cuts
-// load-loop instruction count in half (1 cp.async vs 1 LDG + 1 STS).
-// Requires sm_80+.
-__device__ __forceinline__ void cp_async_16(bf16* smem_ptr, const void* gmem_ptr) {
-    unsigned smem_addr = __cvta_generic_to_shared(smem_ptr);
-    asm volatile("cp.async.ca.shared.global [%0], [%1], 16;"
-                 :: "r"(smem_addr), "l"(gmem_ptr));
-}
-
-__device__ __forceinline__ void cp_async_commit() {
-    asm volatile("cp.async.commit_group;");
-}
-
-__device__ __forceinline__ void cp_async_wait_all() {
-    asm volatile("cp.async.wait_all;");
-}
-
-// Wait until at most N commit groups are still in flight. Used for
-// double-buffered pipelining: wait_group<1> lets the next tile's cp.async
-// continue while ensuring the current tile's data is ready.
-template <int N>
-__device__ __forceinline__ void cp_async_wait_group() {
-    asm volatile("cp.async.wait_group %0;" :: "n"(N));
-}

csrc/kernels/gqa_prefill_attn.cu

@@ -1,100 +0,0 @@
-#include "gqa_prefill_attn.cuh"
-#include <torch/extension.h>
-
-#ifndef ASTRAI_NO_MMA
-#include "gqa_prefill_attn_mma.cuh"
-#endif
-
-template <int HEAD_DIM>
-static void dispatch_prefill(GQAParams& p) {
-#ifndef ASTRAI_NO_MMA
-    constexpr int WARPS = 4, BC = 32, BR = 16, LD = HEAD_DIM;
-    dim3 grid((p.q_len + BR * WARPS - 1) / (BR * WARPS), p.q_head, p.batch);
-    dim3 block(WARPS * 32, 1, 1);
-    // sK + sV (each BC*LD) + shared sQ staging (BR*LD)
-    int smem = (2 * BC * LD + BR * LD) * (int)sizeof(bf16);
-    cudaFuncSetAttribute(gqa_prefill_attn_mma_kernel<HEAD_DIM, WARPS, BC>,
-                         cudaFuncAttributeMaxDynamicSharedMemorySize, smem);
-    gqa_prefill_attn_mma_kernel<HEAD_DIM, WARPS, BC><<<grid, block, smem>>>(p);
-#else
-    constexpr int G = 8, ROWS = 32, P_BC = 32;
-    dim3 grid((p.q_len + ROWS - 1) / ROWS, p.q_head, p.batch);
-    dim3 block(G, ROWS, 1);
-    size_t smem = 2 * P_BC * HEAD_DIM * sizeof(bf16);
-    gqa_prefill_attn_kernel_t<HEAD_DIM, G, ROWS, P_BC><<<grid, block, smem>>>(p);
-#endif
-}
-
-torch::Tensor gqa_prefill_attn(
-    torch::Tensor q,
-    torch::Tensor k,
-    torch::Tensor v,
-    c10::optional<torch::Tensor> mask,
-    bool is_causal = false,
-    int64_t causal_offset = 0,
-    c10::optional<double> scale = c10::nullopt
-) {
-    TORCH_CHECK(q.is_cuda() && k.is_cuda() && v.is_cuda());
-    TORCH_CHECK(q.dtype() == torch::kBFloat16);
-    TORCH_CHECK(k.dtype() == torch::kBFloat16);
-    TORCH_CHECK(v.dtype() == torch::kBFloat16);
-
-    GQAParams p;
-    p.batch = q.size(0);
-    p.q_head = q.size(1);
-    p.kv_head = k.size(1);
-    p.q_len = q.size(2);
-    p.kv_len = k.size(2);
-    p.head_dim = q.size(3);
-    TORCH_CHECK(p.head_dim % 16 == 0, "head_dim must be multiple of 16");
-    p.use_mask = mask.has_value();
-    p.is_causal = (int)is_causal;
-    p.causal_offset = (int)causal_offset;
-    p.scale = scale.has_value() ? (float)scale.value() : 1.0f / sqrtf((float)p.head_dim);
-    p.q = (const bf16*)q.data_ptr();
-    p.k = (const bf16*)k.data_ptr();
-    p.v = (const bf16*)v.data_ptr();
-    if (p.use_mask) {
-        TORCH_CHECK(mask.value().dtype() == torch::kBool);
-        TORCH_CHECK(mask.value().dim() == 2);
-        TORCH_CHECK(mask.value().size(0) == p.batch);
-        TORCH_CHECK(mask.value().size(1) == p.kv_len);
-        p.mask = mask.value().data_ptr<bool>();
-    } else {
-        p.mask = nullptr;
-    }
-
-    auto O = torch::empty_like(q);
-    p.o = (bf16*)O.data_ptr();
-
-    switch (p.head_dim) {
-        case 32:
-            dispatch_prefill<32>(p);
-            break;
-        case 64:
-            dispatch_prefill<64>(p);
-            break;
-        case 128:
-            dispatch_prefill<128>(p);
-            break;
-        case 256:
-            dispatch_prefill<256>(p);
-            break;
-        default:
-            TORCH_CHECK(false, "prefill: unsupported head_dim ", p.head_dim,
-                        " (supported: 32,64,128,256)");
-    }
-    return O;
-}
-
-PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
-    m.def("gqa_prefill_attn", &gqa_prefill_attn,
-        py::arg("q"),
-        py::arg("k"),
-        py::arg("v"),
-        py::arg("mask") = py::none(),
-        py::arg("is_causal") = false,
-        py::arg("causal_offset") = 0,
-        py::arg("scale") = py::none(),
-        "GQA prefill (tensor-core mma on sm_80+, scalar fallback)");
-}

csrc/kernels/gqa_prefill_attn_mma.cuh

@@ -1,246 +0,0 @@
-#pragma once
-#include "gqa_common.cuh"
-#include "gqa_mma_utils.cuh"
-
-// Tensor-core prefill, register-resident flash attention (raw mma.sync PTX).
-// One warp owns BR=16 query rows. S = Q@K^T and O = P@V run on bf16 tensor
-// cores via mma.sync.m16n8k16 (f32 accumulate). Q stays resident in registers;
-// S, O, and the online-softmax stats (m, l) live in registers too — nothing is
-// staged through shared memory except the cooperatively-loaded K/V tiles. The
-// mma fragment layout is used directly: the S accumulator (f32) maps element-
-// for-element onto the P matrix_a (bf16) operand, so softmax needs no shuffle
-// repack; row reductions fold across the 4-lane thread group. Templated on
-// <HEAD_DIM, WARPS, BC> with BC a multiple of 16.
-//
-// Optimizations: shared sQ staging (single area, serialized per-warp load)
-// → cuts smem; pre-scale Q by attention scale during Q load; cp.async global→
-// shared for K/V; scalar fallback only for the last partial tile; causal tile
-// skipping (block-level early break + warp-level skip); XOR swizzle (swiz_col)
-// → eliminates ldmatrix bank conflicts without LD padding (LD=HEAD_DIM).
-
-template <int HEAD_DIM, int WARPS, int BC>
-__global__ void gqa_prefill_attn_mma_kernel(GQAParams p) {
-    constexpr int BR = 16;
-    constexpr int KD = HEAD_DIM / 16;  // Q/K k-tiles
-    constexpr int NC8 = BC / 8;        // S n-tiles (N=8 each)
-    constexpr int KT2 = BC / 16;       // P k-tiles (K=16 each)
-    constexpr int DN8 = HEAD_DIM / 8;  // O n-tiles (N=8 each)
-    constexpr int LD = HEAD_DIM;   // XOR swizzle (swiz_col) handles bank conflicts
-    constexpr int SWIZ_MASK = (HEAD_DIM >= 64) ? 7 : (HEAD_DIM / 8 - 1);  // chunk bits, stay within LD
-
-    const int warp = threadIdx.x / 32;
-    const int lane = threadIdx.x % 32;
-    const int gid = lane >> 2;   // 0..7  → rows gid, gid+8
-    const int tid4 = lane & 3;   // 0..3
-    const int nthreads = WARPS * 32;
-
-    const int q_head = blockIdx.y;
-    const int batch = blockIdx.z;
-    const int kv_head = q_head / (p.q_head / p.kv_head);
-    const int qrow0 = (blockIdx.x * WARPS + warp) * BR;
-
-    extern __shared__ __align__(16) bf16 smem[];
-    bf16* sK = smem;                       // [BC][LD]
-    bf16* sV = sK + BC * LD;               // [BC][LD]
-    bf16* sQ = sV + BC * LD;               // shared staging [BR][LD]
-
-    // Q resident A-fragments (loaded once per warp via shared staging).
-    // Pre-scale by attention scale so softmax doesn't need to multiply later.
-    const int q_base = ((batch * p.q_head + q_head) * p.q_len) * HEAD_DIM;
-    unsigned Qa[KD][4];
-    bf16 scale_bf16 = __float2bfloat16(p.scale);
-    int qrow_l = (lane & 7) + (lane & 8);       // 0..15
-    int qcol_l = (lane & 16) ? 8 : 0;
-    for (int w = 0; w < WARPS; w++) {
-        if (warp == w) {
-            for (int i = lane; i < BR * HEAD_DIM; i += 32) {
-                int r = i / HEAD_DIM, d = i % HEAD_DIM;
-                int qr = qrow0 + r;
-                bf16 qv = (qr < p.q_len) ? p.q[q_base + qr * HEAD_DIM + d]
-                                         : __float2bfloat16(0.0f);
-                sQ[r * LD + swiz_col(d, r, SWIZ_MASK)] = __hmul(qv, scale_bf16);
-            }
-            __syncwarp();
-        #pragma unroll
-            for (int kt = 0; kt < KD; kt++)
-                ldmatrix_x4(Qa[kt], &sQ[qrow_l * LD + swiz_col(kt * 16 + qcol_l, qrow_l, SWIZ_MASK)]);
-        }
-        __syncthreads();  // prevent next warp from overwriting sQ prematurely
-    }
-
-    float Oacc[DN8][4];
-#pragma unroll
-    for (int j = 0; j < DN8; j++)
-        Oacc[j][0] = Oacc[j][1] = Oacc[j][2] = Oacc[j][3] = 0.0f;
-    float m0 = -FLT_MAX, m1 = -FLT_MAX, l0 = 0.0f, l1 = 0.0f;
-
-    const int kv_base = ((batch * p.kv_head + kv_head) * p.kv_len) * HEAD_DIM;
-    const int tiles = (p.kv_len + BC - 1) / BC;
-    const int qr0 = qrow0 + gid;      // row for c0/c1
-    const int qr1 = qrow0 + gid + 8;  // row for c2/c3
-
-    // Causal tile-skip bounds (no-op when is_causal == 0)
-    const int use_skip = p.is_causal;
-    const int max_kv = qrow0 + BR - 1 + p.causal_offset;
-    const int block_max_kv =
-        blockIdx.x * WARPS * BR + WARPS * BR - 1 + p.causal_offset;
-    const int has_mask = p.use_mask && p.mask;
-    const int mb = batch * p.kv_len;
-
-    for (int ti = 0; ti < tiles; ti++) {
-        int kv0 = ti * BC;
-
-        // Block-level causal early break
-        if (use_skip && kv0 > block_max_kv) break;
-
-        // ---- load K/V tile to shared memory (cp.async on full tiles) ----
-        bool full_tile = (kv0 + BC <= p.kv_len);
-        if (full_tile) {
-            constexpr int VEC = 8;  // bf16 per cp.async unit (16 bytes)
-            int total = BC * HEAD_DIM;
-#pragma unroll
-            for (int i = threadIdx.x * VEC; i < total; i += nthreads * VEC) {
-                int r = i / HEAD_DIM;
-                int d = i % HEAD_DIM;
-                int kc = kv0 + r;
-                cp_async_16(&sK[r * LD + swiz_col(d, r, SWIZ_MASK)], &p.k[kv_base + kc * HEAD_DIM + d]);
-                cp_async_16(&sV[r * LD + swiz_col(d, r, SWIZ_MASK)], &p.v[kv_base + kc * HEAD_DIM + d]);
-            }
-            cp_async_commit();
-            cp_async_wait_all();
-        } else {
-            for (int i = threadIdx.x; i < BC * HEAD_DIM; i += nthreads) {
-                int r = i / HEAD_DIM, d = i % HEAD_DIM;
-                int kc = kv0 + r;
-                bf16 z = __float2bfloat16(0.0f);
-                sK[r * LD + swiz_col(d, r, SWIZ_MASK)] = (kc < p.kv_len)
-                                     ? p.k[kv_base + kc * HEAD_DIM + d] : z;
-                sV[r * LD + swiz_col(d, r, SWIZ_MASK)] = (kc < p.kv_len)
-                                     ? p.v[kv_base + kc * HEAD_DIM + d] : z;
-            }
-        }
-        __syncthreads();
-
-        // Warp-level causal skip
-        if (!use_skip || kv0 <= max_kv) {
-
-        // S = Q @ K^T  → Sacc[n8][0..3]   (n8: 8 kv cols each)
-        float Sacc[NC8][4];
-#pragma unroll
-        for (int n8 = 0; n8 < NC8; n8++) {
-            Sacc[n8][0] = Sacc[n8][1] = Sacc[n8][2] = Sacc[n8][3] = 0.0f;
-            int krow_l = n8 * 8 + (lane & 7);
-            int kcol_h = (lane & 8) ? 8 : 0;
-#pragma unroll
-            for (int kt = 0; kt < KD; kt++) {
-                unsigned b[2];
-                ldmatrix_x2(b, &sK[krow_l * LD + swiz_col(kt * 16 + kcol_h, krow_l, SWIZ_MASK)]);
-                mma16816(Sacc[n8], Qa[kt], b, Sacc[n8]);
-            }
-        }
-
-        // ---- online softmax (in registers) ----
-        // Q is pre-scaled, so Sacc already includes the attention scale.
-        int maxc0 = p.is_causal ? min(p.kv_len, qr0 + p.causal_offset + 1)
-                                : p.kv_len;
-        int maxc1 = p.is_causal ? min(p.kv_len, qr1 + p.causal_offset + 1)
-                                : p.kv_len;
-        float rmax0 = -FLT_MAX, rmax1 = -FLT_MAX;
-#pragma unroll
-        for (int n8 = 0; n8 < NC8; n8++) {
-            int cc = kv0 + n8 * 8 + 2 * tid4;
-            int c1 = cc + 1;
-            bool b0 = (cc >= maxc0) || (has_mask && !p.mask[mb + cc]);
-            bool b1 = (c1 >= maxc0) || (has_mask && !p.mask[mb + c1]);
-            bool b2 = (cc >= maxc1) || (has_mask && !p.mask[mb + cc]);
-            bool b3 = (c1 >= maxc1) || (has_mask && !p.mask[mb + c1]);
-            float s0 = b0 ? -FLT_MAX : Sacc[n8][0];
-            float s1 = b1 ? -FLT_MAX : Sacc[n8][1];
-            float s2 = b2 ? -FLT_MAX : Sacc[n8][2];
-            float s3 = b3 ? -FLT_MAX : Sacc[n8][3];
-            Sacc[n8][0] = s0; Sacc[n8][1] = s1;
-            Sacc[n8][2] = s2; Sacc[n8][3] = s3;
-            rmax0 = fmaxf(rmax0, fmaxf(s0, s1));
-            rmax1 = fmaxf(rmax1, fmaxf(s2, s3));
-        }
-        rmax0 = fmaxf(rmax0, __shfl_xor_sync(0xFFFFFFFF, rmax0, 1));
-        rmax0 = fmaxf(rmax0, __shfl_xor_sync(0xFFFFFFFF, rmax0, 2));
-        rmax1 = fmaxf(rmax1, __shfl_xor_sync(0xFFFFFFFF, rmax1, 1));
-        rmax1 = fmaxf(rmax1, __shfl_xor_sync(0xFFFFFFFF, rmax1, 2));
-
-        float nm0 = fmaxf(m0, rmax0), nm1 = fmaxf(m1, rmax1);
-        float corr0 = (nm0 == -FLT_MAX) ? 1.0f : __expf(m0 - nm0);
-        float corr1 = (nm1 == -FLT_MAX) ? 1.0f : __expf(m1 - nm1);
-
-        float rsum0 = 0.0f, rsum1 = 0.0f;
-#pragma unroll
-        for (int n8 = 0; n8 < NC8; n8++) {
-            float p0 = (Sacc[n8][0] == -FLT_MAX) ? 0.0f
-                                                 : __expf(Sacc[n8][0] - nm0);
-            float p1 = (Sacc[n8][1] == -FLT_MAX) ? 0.0f
-                                                 : __expf(Sacc[n8][1] - nm0);
-            float p2 = (Sacc[n8][2] == -FLT_MAX) ? 0.0f
-                                                 : __expf(Sacc[n8][2] - nm1);
-            float p3 = (Sacc[n8][3] == -FLT_MAX) ? 0.0f
-                                                 : __expf(Sacc[n8][3] - nm1);
-            Sacc[n8][0] = p0; Sacc[n8][1] = p1;
-            Sacc[n8][2] = p2; Sacc[n8][3] = p3;
-            rsum0 += p0 + p1;
-            rsum1 += p2 + p3;
-        }
-        rsum0 += __shfl_xor_sync(0xFFFFFFFF, rsum0, 1);
-        rsum0 += __shfl_xor_sync(0xFFFFFFFF, rsum0, 2);
-        rsum1 += __shfl_xor_sync(0xFFFFFFFF, rsum1, 1);
-        rsum1 += __shfl_xor_sync(0xFFFFFFFF, rsum1, 2);
-        l0 = l0 * corr0 + rsum0;
-        l1 = l1 * corr1 + rsum1;
-        m0 = nm0; m1 = nm1;
-
-        // rescale O accumulator by per-row correction
-#pragma unroll
-        for (int j = 0; j < DN8; j++) {
-            Oacc[j][0] *= corr0; Oacc[j][1] *= corr0;
-            Oacc[j][2] *= corr1; Oacc[j][3] *= corr1;
-        }
-
-        // O += P @ V
-#pragma unroll
-        for (int kt2 = 0; kt2 < KT2; kt2++) {
-            unsigned Pa[4];
-            Pa[0] = pk2(Sacc[kt2 * 2][0], Sacc[kt2 * 2][1]);
-            Pa[1] = pk2(Sacc[kt2 * 2][2], Sacc[kt2 * 2][3]);
-            Pa[2] = pk2(Sacc[kt2 * 2 + 1][0], Sacc[kt2 * 2 + 1][1]);
-            Pa[3] = pk2(Sacc[kt2 * 2 + 1][2], Sacc[kt2 * 2 + 1][3]);
-            int vrow_l = kt2 * 16 + (lane & 15);
-#pragma unroll
-            for (int dn8 = 0; dn8 < DN8; dn8++) {
-                unsigned b[2];
-                ldmatrix_x2_trans(b, &sV[vrow_l * LD + swiz_col(dn8 * 8, vrow_l, SWIZ_MASK)]);
-                mma16816(Oacc[dn8], Pa, b, Oacc[dn8]);
-            }
-        }
-        }  // if active (warp-level causal skip)
-        __syncthreads();
-    }
-
-    // ---- write output ----
-    float rl0 = (l0 > 1e-20f) ? (1.0f / l0) : 0.0f;
-    float rl1 = (l1 > 1e-20f) ? (1.0f / l1) : 0.0f;
-    const int o_base = ((batch * p.q_head + q_head) * p.q_len) * HEAD_DIM;
-#pragma unroll
-    for (int dn8 = 0; dn8 < DN8; dn8++) {
-        int d = dn8 * 8 + 2 * tid4;
-        if (qr0 < p.q_len) {
-            p.o[o_base + qr0 * HEAD_DIM + d] =
-                __float2bfloat16(Oacc[dn8][0] * rl0);
-            p.o[o_base + qr0 * HEAD_DIM + d + 1] =
-                __float2bfloat16(Oacc[dn8][1] * rl0);
-        }
-        if (qr1 < p.q_len) {
-            p.o[o_base + qr1 * HEAD_DIM + d] =
-                __float2bfloat16(Oacc[dn8][2] * rl1);
-            p.o[o_base + qr1 * HEAD_DIM + d + 1] =
-                __float2bfloat16(Oacc[dn8][3] * rl1);
-        }
-    }
-}

csrc/tests/attn_decode_test.cu

@@ -0,0 +1,199 @@
+/*
+Pure-C test:
+nvcc -I csrc -arch=sm_89 -O3 \
+    --use_fast_math --ptxas-options=-O3 --extra-device-vectorization \
+    csrc/tests/attn_decode_test.cu -o test && ./test
+*/
+
+#include "test_utils.cuh"
+#include "../kernels/attn_decode_split_kv.cuh"
+#ifndef ASTRAI_NO_MMA
+#include "../kernels/attn_decode_split_kv_mma.cuh"
+#endif
+
+// Split-K scratch (torch-free): the production launcher allocates these from
+// torch; here we pass pre-allocated device buffers so the bench loop doesn't
+// pay a cudaMalloc per iteration. Size for the maximum split count (32).
+struct DecodeScratch {
+    float* o_part = nullptr;
+    float* ml_part = nullptr;
+};
+
+// Launch the production decode path (tensor-core head-packing MMA on sm_80+,
+// scalar fallback otherwise), mirroring dispatch_decode() in attn_decode.cu.
+#ifndef ASTRAI_NO_MMA
+static bool decode_use_mma(const AttentionParams<bf16>& p) {
+    int G = p.q_head / p.kv_head;
+    return !p.use_mask && G > 1 && G <= 16;
+}
+
+template <int HEAD_DIM, int BC, int STAGES = (HEAD_DIM <= 128) ? 2 : 1>
+static void launch_mma_decode(AttentionParams<bf16>& p, DecodeScratch& sc) {
+    int tiles_total = (p.kv_len + BC - 1) / BC;
+    p.num_splits = compute_num_splits(p.batch * p.kv_head, tiles_total);
+    p.o_part = sc.o_part;
+    p.ml_part = sc.ml_part;
+
+    attn_decode_split_kv_mma_kernel<HEAD_DIM, BC, STAGES>
+        <<<dim3(p.kv_head, p.batch, p.num_splits), 32>>>(p);
+    attn_decode_combine_kernel<<<p.batch * p.q_head, p.head_dim>>>(p);
+}
+#endif
+
+static void launch_scalar_decode(AttentionParams<bf16>& p, DecodeScratch& sc) {
+    int gs = p.q_head / p.kv_head;
+    int chunks_total = (p.kv_len + DC_CHUNK - 1) / DC_CHUNK;
+    p.num_splits = compute_num_splits(p.batch * p.kv_head, chunks_total);
+    p.o_part = sc.o_part;
+    p.ml_part = sc.ml_part;
+
+    size_t smem = DC_CHUNK * p.head_dim * sizeof(bf16);
+    attn_decode_split_kv_kernel<<<dim3(p.batch * p.kv_head, 1, p.num_splits), dim3(32, gs), smem>>>(p);
+    attn_decode_combine_kernel<<<p.batch * p.q_head, p.head_dim>>>(p);
+}
+
+template <int HEAD_DIM>
+static void dispatch_decode_t(AttentionParams<bf16>& p, DecodeScratch& sc) {
+#ifndef ASTRAI_NO_MMA
+    if (decode_use_mma(p)) { launch_mma_decode<HEAD_DIM, 32>(p, sc); return; }
+#endif
+    launch_scalar_decode(p, sc);
+}
+
+static void dispatch_decode(AttentionParams<bf16>& p, DecodeScratch& sc) {
+    dispatch_by_head_dim(p.head_dim, [&]<int D>() { dispatch_decode_t<D>(p, sc); });
+}
+
+// Warmed-up, CUDA-event timed sweep over the production decode MMA path.
+static void bench() {
+    const int cfgs[][5] = {
+        {1, 32, 4, 512, 128},    // B, Hq, Hk, kv_len, D
+        {1, 32, 4, 1024, 128},
+        {1, 32, 4, 2048, 128},
+        {1, 32, 4, 4096, 128},
+        {16, 32, 4, 2048, 128},
+        {32, 32, 4, 1024, 128},
+    };
+    const int WARMUP = 10, ITERS = 100;
+    printf("\n===== DECODE BENCH (warmup=%d iters=%d) =====\n", WARMUP, ITERS);
+    print_bench_header();
+
+    for (int ci = 0; ci < 6; ci++) {
+        int B = cfgs[ci][0], Hq = cfgs[ci][1], Hk = cfgs[ci][2];
+        int sl = cfgs[ci][3], D = cfgs[ci][4];
+        size_t nQ = (size_t)B * Hq * D;
+        size_t nKV = (size_t)B * Hk * sl * D;
+
+        bf16 *dQ, *dK, *dV, *dO;
+        cudaMalloc(&dQ, nQ*2); cudaMalloc(&dK, nKV*2);
+        cudaMalloc(&dV, nKV*2); cudaMalloc(&dO, nQ*2);
+        size_t big = nQ > nKV ? nQ : nKV; bf16* tmp = new bf16[big];
+        for (size_t i = 0; i < nQ; i++)  tmp[i] = f2bf(randf());
+        cudaMemcpy(dQ, tmp, nQ*2, cudaMemcpyHostToDevice);
+        for (size_t i = 0; i < nKV; i++) tmp[i] = f2bf(randf());
+        cudaMemcpy(dK, tmp, nKV*2, cudaMemcpyHostToDevice);
+        for (size_t i = 0; i < nKV; i++) tmp[i] = f2bf(randf());
+        cudaMemcpy(dV, tmp, nKV*2, cudaMemcpyHostToDevice);
+        delete[] tmp;
+
+        AttentionParams<bf16> p;
+        p.batch = B; p.q_head = Hq; p.kv_head = Hk; p.q_len = 1; p.kv_len = sl;
+        p.head_dim = D; p.use_mask = 0; p.is_causal = 0; p.causal_offset = 0;
+        p.scale = 1.0f / sqrtf((float)D);
+        p.q = dQ; p.k = dK; p.v = dV; p.mask = nullptr; p.o = dO;
+
+        DecodeScratch sc;
+        cudaMalloc(&sc.o_part, (size_t)B*Hq*32*D*sizeof(float));
+        cudaMalloc(&sc.ml_part, (size_t)B*Hq*32*2*sizeof(float));
+
+        auto launch = [&]() { dispatch_decode(p, sc); };
+        double flops = 4.0 * B * Hq * (double)sl * D;
+        double bytes = 2.0 * (2.0 * nKV * sizeof(bf16));
+        BenchResult r = bench_kernel(launch, WARMUP, ITERS, flops, bytes);
+
+        char cfg[64];
+        snprintf(cfg, sizeof(cfg),
+                 "B=%2d Hq=%2d Hk=%d q=%4d kv=%4d D=%3d causal=%d",
+                 B, Hq, Hk, 1, sl, D, 0);
+        print_bench_row(cfg, r);
+
+        cudaFree(dQ); cudaFree(dK); cudaFree(dV); cudaFree(dO);
+        cudaFree(sc.o_part); cudaFree(sc.ml_part);
+    }
+}
+
+int main() {
+    const int configs[][5] = {
+        {1, 2, 1, 64, 32},    // B,Hq,Hk,seq_len,D
+        {1, 32, 4, 512, 128},
+        {1, 32, 4, 1024, 128},
+    };
+    int n_cfgs = sizeof(configs) / sizeof(configs[0]);
+
+    for (int ci = 0; ci < n_cfgs; ci++) {
+        int B = configs[ci][0], Hq = configs[ci][1], Hk = configs[ci][2];
+        int sl = configs[ci][3], D = configs[ci][4], gs = Hq / Hk;
+        printf("=== B=%d Hq=%d Hk=%d seq=%d D=%d gs=%d ===\n", B,Hq,Hk,sl,D,gs);
+
+        size_t nQ = B*Hq*1*D, nKV = B*Hk*sl*D;
+        float *hQ=new float[nQ], *hK=new float[nKV], *hV=new float[nKV];
+        for (size_t i=0;i<nQ;i++) hQ[i]=randf();
+        for (size_t i=0;i<nKV;i++){hK[i]=randf();hV[i]=randf();}
+
+        bool* hMask=new bool[B*sl];
+        for (int i=0;i<B*sl;i++) hMask[i]=true;
+
+        bf16 *dQ,*dK,*dV,*dO,*tmp;
+        bool* dMask;
+        cudaMalloc(&dQ,nQ*2); cudaMalloc(&dK,nKV*2);
+        cudaMalloc(&dV,nKV*2); cudaMalloc(&dO,nQ*2);
+        cudaMalloc(&dMask,B*sl);
+
+        tmp=new bf16[max(nQ,nKV)];
+        for (size_t i=0;i<nQ;i++) tmp[i]=f2bf(hQ[i]);
+        cudaMemcpy(dQ,tmp,nQ*2,cudaMemcpyHostToDevice);
+        for (size_t i=0;i<nKV;i++) tmp[i]=f2bf(hK[i]);
+        cudaMemcpy(dK,tmp,nKV*2,cudaMemcpyHostToDevice);
+        for (size_t i=0;i<nKV;i++) tmp[i]=f2bf(hV[i]);
+        cudaMemcpy(dV,tmp,nKV*2,cudaMemcpyHostToDevice);
+        cudaMemcpy(dMask,hMask,B*sl,cudaMemcpyHostToDevice);
+
+        AttentionParams<bf16> p;
+        p.batch=B; p.q_head=Hq; p.kv_head=Hk; p.q_len=1; p.kv_len=sl; p.head_dim=D;
+        p.use_mask=0; p.is_causal=0; p.causal_offset=0;
+        p.scale=1.0f/sqrtf((float)D);
+        p.q=dQ; p.k=dK; p.v=dV; p.mask=nullptr; p.o=dO;
+
+        // Split-K scratch (max 32 splits), sized for the production MMA path.
+        DecodeScratch sc;
+        cudaMalloc(&sc.o_part, (size_t)B*Hq*32*D*sizeof(float));
+        cudaMalloc(&sc.ml_part, (size_t)B*Hq*32*2*sizeof(float));
+
+        double t0=now_ms();
+        dispatch_decode(p, sc);
+        cudaDeviceSynchronize();
+        double kms=now_ms()-t0;
+        cudaError_t err=cudaGetLastError();
+        if (err!=cudaSuccess){printf("CUDA err: %s\n",cudaGetErrorString(err));return 1;}
+
+        bf16* hOut=new bf16[nQ];
+        cudaMemcpy(hOut,dO,nQ*2,cudaMemcpyDeviceToHost);
+
+        float* ref=new float[nQ];
+        cpu_attention_ref(hQ, hK, hV, hMask, ref, B, Hq, Hk, 1, sl, D, 0, 0);
+
+        float max_err=0;
+        for (size_t i=0;i<nQ;i++){
+            float d=fabsf(bf2f(hOut[i])-ref[i]);
+            if(d>max_err) max_err=d;
+        }
+        printf("kernel: %.3f ms  max_err: %.6e\n\n",kms,max_err);
+
+        cudaFree(dQ);cudaFree(dK);cudaFree(dV);cudaFree(dO);cudaFree(dMask);
+        cudaFree(sc.o_part);cudaFree(sc.ml_part);
+        delete[]hQ;delete[]hK;delete[]hV;delete[]hMask;delete[]hOut;delete[]ref;delete[]tmp;
+    }
+    printf("All tests passed!\n");
+    bench();
+    return 0;
+}

csrc/tests/attn_paged_decode_test.cu

@@ -0,0 +1,330 @@
+// Compile:
+//   nvcc -I csrc -arch=sm_89 -O3 --use_fast_math --ptxas-options=-O3 \
+//        --extra-device-vectorization csrc/tests/attn_paged_decode_test.cu \
+//        -o /tmp/test_paged && /tmp/test_paged
+
+#include <cstring>
+#include "test_utils.cuh"
+#include "../kernels/attn_paged_decode_split_kv.cuh"
+#ifndef ASTRAI_NO_MMA
+#include "../kernels/attn_paged_decode_split_kv_mma.cuh"
+#endif
+
+// Copy contiguous K/V from page pool (reference gather)
+static void gather_kv_cpu(
+    const bf16* h_k_pool, const bf16* h_v_pool,
+    const int64_t* h_pt, int B, int Hkv, int kv_len,
+    int page_size, int head_dim,
+    bf16* h_k, bf16* h_v)
+{
+    int max_pages = (kv_len + page_size - 1) / page_size;
+    size_t page_stride = (size_t)page_size * Hkv * head_dim;
+    for (int b = 0; b < B; b++) {
+        for (int pos = 0; pos < kv_len; pos++) {
+            int log_pg = pos / page_size;
+            int pg_off = pos % page_size;
+            int phys = (int)h_pt[b * max_pages + log_pg];
+            for (int h = 0; h < Hkv; h++) {
+                size_t src_base = (size_t)phys * page_stride
+                                + (size_t)pg_off * Hkv * head_dim
+                                + h * head_dim;
+                size_t dst_base = ((size_t)b * Hkv + h) * kv_len * head_dim + (size_t)pos * head_dim;
+                memcpy(h_k + dst_base, h_k_pool + src_base, head_dim * sizeof(bf16));
+                memcpy(h_v + dst_base, h_v_pool + src_base, head_dim * sizeof(bf16));
+            }
+        }
+    }
+}
+
+template <int HEAD_DIM>
+static void launch_paged_decode(PagedAttentionParams<bf16, float>& p) {
+#ifndef ASTRAI_NO_MMA
+    int G_check = p.q_head / p.kv_head;
+    bool use_mma = !p.use_mask && G_check >= 1 && G_check <= 16 && p.page_size >= 32;
+    if (use_mma) {
+        constexpr int STAGES = (HEAD_DIM <= 128) ? 2 : 1;
+        int tiles_total = (p.kv_len + 32 - 1) / 32;
+        p.num_splits = compute_num_splits(p.batch * p.kv_head, tiles_total);
+        paged_attn_decode_split_kv_mma_kernel<HEAD_DIM, 32, STAGES>
+            <<<dim3(p.kv_head, p.batch, p.num_splits), 32>>>(p);
+    } else
+#endif
+    {
+        int group_sz = p.q_head / p.kv_head;
+        int chunks_total = (p.kv_len + PDC_CHUNK - 1) / PDC_CHUNK;
+        p.num_splits = compute_num_splits(p.batch * p.kv_head, chunks_total);
+        size_t smem = PDC_CHUNK * p.head_dim * sizeof(bf16);
+        paged_attn_decode_split_kv_kernel<<<
+            dim3(p.batch * p.kv_head, 1, p.num_splits),
+            dim3(32, group_sz), smem>>>(p);
+    }
+    paged_attn_decode_combine_kernel<<<p.batch * p.q_head, p.head_dim>>>(p);
+}
+
+template <int HEAD_DIM>
+static int run_test(int B, int Hq, int Hkv, int kv_len, int page_size, int seed) {
+    printf("B=%d Hq=%d Hkv=%d kv_len=%d page_sz=%d head_dim=%d ... ", B, Hq, Hkv, kv_len, page_size, HEAD_DIM);
+    fflush(stdout);
+
+    int max_pages = (kv_len + page_size - 1) / page_size;
+    int n_phys_pages = B * max_pages;
+
+    size_t sz_q  = (size_t)B * Hq * 1 * HEAD_DIM * sizeof(bf16);
+    size_t sz_o  = sz_q;
+    size_t sz_kv = (size_t)n_phys_pages * page_size * Hkv * HEAD_DIM * sizeof(bf16);
+    size_t sz_pt = (size_t)B * max_pages * sizeof(int64_t);
+    int max_splits = 32;
+    size_t sz_op = (size_t)B * Hq * max_splits * HEAD_DIM * sizeof(float);
+    size_t sz_ml = (size_t)B * Hq * max_splits * 2 * sizeof(float);
+
+    bf16 *d_q, *d_o_paged, *d_o_ref;
+    bf16 *d_k_pool, *d_v_pool;
+    int64_t* d_pt;
+    float *d_op, *d_ml;
+
+    cudaMalloc(&d_q, sz_q);
+    cudaMalloc(&d_o_paged, sz_o);
+    cudaMalloc(&d_o_ref, sz_o);
+    cudaMalloc(&d_k_pool, sz_kv);
+    cudaMalloc(&d_v_pool, sz_kv);
+    cudaMalloc(&d_pt, sz_pt);
+    cudaMalloc(&d_op, sz_op);
+    cudaMalloc(&d_ml, sz_ml);
+
+    srand(seed);
+    auto rnd = [&]() { return (rand() / (float)RAND_MAX) * 2.0f - 1.0f; };
+
+    bf16* h_q = (bf16*)malloc(sz_q);
+    for (int i = 0; i < B * Hq * HEAD_DIM; i++)
+        h_q[i] = __float2bfloat16(rnd());
+    cudaMemcpy(d_q, h_q, sz_q, cudaMemcpyHostToDevice);
+
+    bf16* h_k_pool = (bf16*)malloc(sz_kv);
+    bf16* h_v_pool = (bf16*)malloc(sz_kv);
+    size_t ps = (size_t)page_size * Hkv * HEAD_DIM;
+    for (int pg = 0; pg < n_phys_pages; pg++) {
+        for (int off = 0; off < page_size; off++) {
+            for (int h = 0; h < Hkv; h++) {
+                for (int d = 0; d < HEAD_DIM; d++) {
+                    float v = sinf((float)(pg * 7919 + off * 1049 + h * 331 + d));
+                    size_t idx = (size_t)pg * ps + (size_t)off * Hkv * HEAD_DIM + h * HEAD_DIM + d;
+                    h_k_pool[idx] = __float2bfloat16(v);
+                    h_v_pool[idx] = __float2bfloat16(v * 0.3f);
+                }
+            }
+        }
+    }
+    cudaMemcpy(d_k_pool, h_k_pool, sz_kv, cudaMemcpyHostToDevice);
+    cudaMemcpy(d_v_pool, h_v_pool, sz_kv, cudaMemcpyHostToDevice);
+
+    int64_t* h_pt = (int64_t*)malloc(sz_pt);
+    int next_pg = 0;
+    for (int b = 0; b < B; b++)
+        for (int p = 0; p < max_pages; p++)
+            h_pt[b * max_pages + p] = next_pg++;
+    cudaMemcpy(d_pt, h_pt, sz_pt, cudaMemcpyHostToDevice);
+
+    bf16* h_k_cont = (bf16*)malloc((size_t)B * kv_len * Hkv * HEAD_DIM * sizeof(bf16));
+    bf16* h_v_cont = (bf16*)malloc((size_t)B * kv_len * Hkv * HEAD_DIM * sizeof(bf16));
+    gather_kv_cpu(h_k_pool, h_v_pool, h_pt, B, Hkv, kv_len, page_size, HEAD_DIM, h_k_cont, h_v_cont);
+
+    float* h_q_f = (float*)malloc((size_t)B * Hq * HEAD_DIM * sizeof(float));
+    float* h_k_f = (float*)malloc((size_t)B * kv_len * Hkv * HEAD_DIM * sizeof(float));
+    float* h_v_f = (float*)malloc((size_t)B * kv_len * Hkv * HEAD_DIM * sizeof(float));
+    for (int i = 0; i < B * Hq * HEAD_DIM; i++) h_q_f[i] = bf2f(h_q[i]);
+    for (int i = 0; i < B * kv_len * Hkv * HEAD_DIM; i++) {
+        h_k_f[i] = bf2f(h_k_cont[i]);
+        h_v_f[i] = bf2f(h_v_cont[i]);
+    }
+
+    float* h_o_ref = (float*)calloc(B * Hq * HEAD_DIM, sizeof(float));
+    cpu_attention_ref(h_q_f, h_k_f, h_v_f, nullptr, h_o_ref, B, Hq, Hkv, 1, kv_len, HEAD_DIM, 0, 0);
+
+    float scale_val = 1.0f / sqrtf((float)HEAD_DIM);
+    PagedAttentionParams<bf16, float> p;
+    p.batch = B; p.q_head = Hq; p.kv_head = Hkv; p.q_len = 1;
+    p.kv_len = kv_len; p.head_dim = HEAD_DIM;
+    p.use_mask = 0; p.is_causal = 0; p.causal_offset = 0;
+    p.num_splits = 1; p.scale = scale_val;
+    p.page_size = page_size; p.max_pages = max_pages;
+    p.page_table = d_pt;
+    p.k_cache = d_k_pool; p.v_cache = d_v_pool;
+    p.q = d_q; p.mask = nullptr; p.o = d_o_paged;
+    p.o_part = d_op; p.ml_part = d_ml;
+
+    launch_paged_decode<HEAD_DIM>(p);
+    cudaDeviceSynchronize();
+
+    bf16* h_o_bf16 = (bf16*)malloc(sz_o);
+    cudaMemcpy(h_o_bf16, d_o_paged, sz_o, cudaMemcpyDeviceToHost);
+    float* h_o_paged = (float*)malloc(B * Hq * HEAD_DIM * sizeof(float));
+    for (int i = 0; i < B * Hq * HEAD_DIM; i++)
+        h_o_paged[i] = __bfloat162float(h_o_bf16[i]);
+
+    float max_err = 0.0f;
+    int bad_idx = -1;
+    for (int i = 0; i < B * Hq * HEAD_DIM; i++) {
+        float e = fabsf(h_o_paged[i] - h_o_ref[i]);
+        if (e > max_err) { max_err = e; bad_idx = i; }
+    }
+
+    bool pass = max_err < 0.02f;
+
+    if (pass) {
+        printf("PASS (max_abs_err=%.4e)\n", max_err);
+    } else {
+        int b = bad_idx / (Hq * HEAD_DIM);
+        int h = (bad_idx / HEAD_DIM) % Hq;
+        int d = bad_idx % HEAD_DIM;
+        printf("FAIL (max_abs_err=%.4e at [%d,%d,%d]: ref=%.4f got=%.4f)\n",
+               max_err, b, h, d, h_o_ref[bad_idx], h_o_paged[bad_idx]);
+        printf("  ref[0..7]:");
+        for (int i = 0; i < 8 && i < HEAD_DIM; i++)
+            printf(" %.4f", h_o_ref[i]);
+        printf("\n  got[0..7]:");
+        for (int i = 0; i < 8 && i < HEAD_DIM; i++)
+            printf(" %.4f", h_o_paged[i]);
+        printf("\n");
+    }
+
+    free(h_q); free(h_k_pool); free(h_v_pool); free(h_pt);
+    free(h_k_cont); free(h_v_cont);
+    free(h_q_f); free(h_k_f); free(h_v_f);
+    free(h_o_ref); free(h_o_bf16); free(h_o_paged);
+    cudaFree(d_q); cudaFree(d_o_paged); cudaFree(d_o_ref);
+    cudaFree(d_k_pool); cudaFree(d_v_pool); cudaFree(d_pt);
+    cudaFree(d_op); cudaFree(d_ml);
+
+    return pass ? 0 : 1;
+}
+
+struct TestCase {
+    int head_dim;
+    int B, Hq, Hkv, kv_len, page_size, seed;
+};
+
+static const TestCase TESTS[] = {
+    {128, 1, 1, 1, 8, 128, 1},
+    {128, 1, 4, 4, 128, 128, 2},
+    {128, 2, 4, 4, 256, 128, 3},
+    {128, 1, 4, 1, 64, 64, 4},
+    {128, 1, 8, 2, 64, 128, 5},
+    {128, 2, 16, 4, 128, 128, 6},
+    {64, 1, 4, 2, 32, 128, 7},
+    {256, 1, 2, 1, 16, 128, 8},
+    {32, 1, 4, 2, 32, 64, 9},
+    {128, 3, 8, 2, 256, 128, 10},
+    {128, 2, 32, 8, 512, 128, 11},
+#ifndef ASTRAI_NO_MMA
+    {128, 1, 16, 2, 256, 128, 12},
+    {128, 2, 32, 4, 512, 128, 13},
+#endif
+};
+
+static int dispatch_test(const TestCase& tc) {
+    bool matched = false;
+    int r = 0;
+    dispatch_by_head_dim(tc.head_dim, [&]<int D>() {
+        matched = true;
+        r = run_test<D>(tc.B, tc.Hq, tc.Hkv, tc.kv_len, tc.page_size, tc.seed);
+    });
+    return matched ? r : 1;
+}
+
+// Warmed-up, CUDA-event timed sweep over paged decode configs.
+// Bytes = K + V read through page table (B*Hk*kv*D each), bf16.
+template <int HEAD_DIM>
+static void bench_config(int B, int Hq, int Hkv, int kv_len, int page_size) {
+    int max_pages = (kv_len + page_size - 1) / page_size;
+    int n_phys_pages = B * max_pages;
+
+    size_t sz_q  = (size_t)B * Hq * 1 * HEAD_DIM * sizeof(bf16);
+    size_t sz_kv = (size_t)n_phys_pages * page_size * Hkv * HEAD_DIM * sizeof(bf16);
+    size_t sz_pt = (size_t)B * max_pages * sizeof(int64_t);
+    int max_splits = 32;
+    size_t sz_op = (size_t)B * Hq * max_splits * HEAD_DIM * sizeof(float);
+    size_t sz_ml = (size_t)B * Hq * max_splits * 2 * sizeof(float);
+
+    bf16 *d_q, *d_o, *d_k_pool, *d_v_pool;
+    int64_t* d_pt;
+    float *d_op, *d_ml;
+    cudaMalloc(&d_q, sz_q); cudaMalloc(&d_o, sz_q);
+    cudaMalloc(&d_k_pool, sz_kv); cudaMalloc(&d_v_pool, sz_kv);
+    cudaMalloc(&d_pt, sz_pt);
+    cudaMalloc(&d_op, sz_op); cudaMalloc(&d_ml, sz_ml);
+
+    bf16* tmp = (bf16*)malloc(sz_kv > sz_q ? sz_kv : sz_q);
+    for (size_t i = 0; i < sz_q / sizeof(bf16); i++) tmp[i] = f2bf(randf());
+    cudaMemcpy(d_q, tmp, sz_q, cudaMemcpyHostToDevice);
+    for (size_t i = 0; i < sz_kv / sizeof(bf16); i++) tmp[i] = f2bf(randf());
+    cudaMemcpy(d_k_pool, tmp, sz_kv, cudaMemcpyHostToDevice);
+    cudaMemcpy(d_v_pool, tmp, sz_kv, cudaMemcpyHostToDevice);
+
+    int64_t* h_pt = (int64_t*)malloc(sz_pt);
+    int next_pg = 0;
+    for (int b = 0; b < B; b++)
+        for (int p = 0; p < max_pages; p++)
+            h_pt[b * max_pages + p] = next_pg++;
+    cudaMemcpy(d_pt, h_pt, sz_pt, cudaMemcpyHostToDevice);
+    free(h_pt);
+
+    float scale_val = 1.0f / sqrtf((float)HEAD_DIM);
+    PagedAttentionParams<bf16, float> pa;
+    pa.batch = B; pa.q_head = Hq; pa.kv_head = Hkv; pa.q_len = 1;
+    pa.kv_len = kv_len; pa.head_dim = HEAD_DIM;
+    pa.use_mask = 0; pa.is_causal = 0; pa.causal_offset = 0;
+    pa.num_splits = 1; pa.scale = scale_val;
+    pa.page_size = page_size; pa.max_pages = max_pages;
+    pa.page_table = d_pt;
+    pa.k_cache = d_k_pool; pa.v_cache = d_v_pool;
+    pa.q = d_q; pa.mask = nullptr; pa.o = d_o;
+    pa.o_part = d_op; pa.ml_part = d_ml;
+
+    const int WARMUP = 10, ITERS = 100;
+    auto launch = [&]() { launch_paged_decode<HEAD_DIM>(pa); };
+    double flops = 4.0 * B * Hq * (double)kv_len * HEAD_DIM;
+    size_t nKV = (size_t)B * Hkv * kv_len * HEAD_DIM;
+    double bytes = 2.0 * (2.0 * nKV * sizeof(bf16));
+    BenchResult r = bench_kernel(launch, WARMUP, ITERS, flops, bytes);
+
+    char cfg[64];
+    snprintf(cfg, sizeof(cfg),
+             "B=%2d Hq=%2d Hk=%d q=%4d kv=%4d D=%3d page=%3d",
+             B, Hq, Hkv, 1, kv_len, HEAD_DIM, page_size);
+    print_bench_row(cfg, r);
+
+    free(tmp);
+    cudaFree(d_q); cudaFree(d_o);
+    cudaFree(d_k_pool); cudaFree(d_v_pool); cudaFree(d_pt);
+    cudaFree(d_op); cudaFree(d_ml);
+}
+
+static void bench() {
+    printf("\n===== PAGED DECODE BENCH =====\n");
+    print_bench_header();
+    bench_config<128>(1, 32, 4, 512, 128);
+    bench_config<128>(1, 32, 4, 1024, 128);
+    bench_config<128>(1, 32, 4, 2048, 128);
+    bench_config<128>(1, 32, 4, 4096, 128);
+    bench_config<128>(16, 32, 4, 2048, 128);
+    bench_config<128>(32, 32, 4, 1024, 128);
+}
+
+int main() {
+    int n = sizeof(TESTS) / sizeof(TESTS[0]);
+    int fail = 0;
+    printf("=== Paged Decode vs CPU reference (%d cases) ===\n\n", n);
+
+    for (int i = 0; i < n; i++) {
+        fail += dispatch_test(TESTS[i]);
+        if (fail) break;
+    }
+
+    if (fail) {
+        printf("\nFAILED (%d/%d tests failed)\n", fail, n);
+        return fail;
+    }
+    printf("\nAll %d tests passed!\n", n);
+    bench();
+    return 0;
+}

csrc/tests/attn_prefill_test.cu

@@ -0,0 +1,176 @@
+/*
+Pure-C test:
+nvcc -I csrc -arch=sm_89 -O3 \
+    --use_fast_math --ptxas-options=-O3 --extra-device-vectorization \
+    csrc/tests/attn_prefill_test.cu -o test && ./test
+*/
+
+#include "test_utils.cuh"
+#include "../kernels/attn_prefill_split_q.cuh"
+#ifndef ASTRAI_NO_MMA
+#include "../kernels/attn_prefill_split_q_mma.cuh"
+#endif
+
+// Launch the production prefill path (tensor-core MMA on sm_80+, else the
+// scalar fallback), mirroring dispatch_prefill() in attn_prefill.cu.
+template <int HEAD_DIM>
+static void launch_prefill(AttentionParams<bf16>& p) {
+#ifndef ASTRAI_NO_MMA
+    constexpr int WARPS = 4, BR = 16;
+    constexpr int BC = (HEAD_DIM <= 128) ? 32 : 16;
+    dim3 grid((p.q_len + BR * WARPS - 1) / (BR * WARPS), p.q_head, p.batch);
+    dim3 block(WARPS * 32, 1, 1);
+    attn_prefill_split_q_mma_kernel<HEAD_DIM, WARPS, BC><<<grid, block>>>(p);
+#else
+    constexpr int G = 8, ROWS = 32, P_BC = 32;
+    dim3 grid((p.q_len + ROWS - 1) / ROWS, p.q_head, p.batch);
+    dim3 block(G, ROWS, 1);
+    attn_prefill_split_q_kernel_t<HEAD_DIM, G, ROWS, P_BC><<<grid, block>>>(p);
+#endif
+}
+
+static void dispatch_prefill(AttentionParams<bf16>& p) {
+    switch (p.head_dim) {
+        case 64:  launch_prefill<64>(p);  break;
+        case 128: launch_prefill<128>(p); break;
+        default:  printf("bench: unsupported D=%d\n", p.head_dim);
+    }
+}
+
+// Warmed-up, CUDA-event timed throughput sweep over the production MMA path.
+// Reports per-call latency and effective tensor-core TFLOP/s (2 matmuls:
+// QK^T and P@V, each 2*B*Hq*ql*kl*D flops; halved for causal).
+static void bench() {
+    const int cfgs[][7] = {
+        {1,32,4,512,512,128,0},
+        {1,32,4,1024,1024,128,0},
+        {1,32,4,2048,2048,128,0},
+        {1,32,4,2048,2048,128,1},
+        {4,32,4,2048,2048,128,1},
+        {1,32,4,4096,4096,128,1},
+    };
+    int n = sizeof(cfgs)/sizeof(cfgs[0]);
+    const int WARMUP = 10, ITERS = 50;
+    printf("\n===== PREFILL BENCH (warmup=%d iters=%d) =====\n", WARMUP, ITERS);
+    printf("%-46s | %10s | %10s | %10s\n",
+           "config", "latency", "bandwidth", "throughput");
+    printf("---------------------------------------------------------------"
+           "----------------------------\n");
+
+    for (int ci = 0; ci < n; ci++) {
+        int B=cfgs[ci][0], Hq=cfgs[ci][1], Hk=cfgs[ci][2];
+        int ql=cfgs[ci][3], kl=cfgs[ci][4], D=cfgs[ci][5], causal=cfgs[ci][6];
+        size_t nQ=(size_t)B*Hq*ql*D, nKV=(size_t)B*Hk*kl*D;
+
+        bf16 *dQ,*dK,*dV,*dO,*tmp;
+        cudaMalloc(&dQ,nQ*2); cudaMalloc(&dK,nKV*2);
+        cudaMalloc(&dV,nKV*2); cudaMalloc(&dO,nQ*2);
+        size_t big = nQ>nKV?nQ:nKV; tmp=new bf16[big];
+        for (size_t i=0;i<nQ;i++)  tmp[i]=f2bf(randf());
+        cudaMemcpy(dQ,tmp,nQ*2,cudaMemcpyHostToDevice);
+        for (size_t i=0;i<nKV;i++) tmp[i]=f2bf(randf());
+        cudaMemcpy(dK,tmp,nKV*2,cudaMemcpyHostToDevice);
+        for (size_t i=0;i<nKV;i++) tmp[i]=f2bf(randf());
+        cudaMemcpy(dV,tmp,nKV*2,cudaMemcpyHostToDevice);
+
+        AttentionParams<bf16> p;
+        p.batch=B; p.q_head=Hq; p.kv_head=Hk; p.q_len=ql; p.kv_len=kl; p.head_dim=D;
+        p.use_mask=0; p.is_causal=causal; p.causal_offset=0;
+        p.scale=1.0f/sqrtf((float)D);
+        p.q=dQ; p.k=dK; p.v=dV; p.mask=nullptr; p.o=dO;
+
+        for (int i=0;i<WARMUP;i++) dispatch_prefill(p);
+        cudaDeviceSynchronize();
+        cudaError_t err=cudaGetLastError();
+        if (err!=cudaSuccess){printf("CUDA err: %s\n",cudaGetErrorString(err));return;}
+
+        cudaEvent_t s,e; cudaEventCreate(&s); cudaEventCreate(&e);
+        cudaEventRecord(s);
+        for (int i=0;i<ITERS;i++) dispatch_prefill(p);
+        cudaEventRecord(e); cudaEventSynchronize(e);
+        float ms=0; cudaEventElapsedTime(&ms,s,e); ms/=ITERS;
+
+        double flops = 4.0*B*Hq*(double)ql*kl*D;
+        if (causal) flops *= 0.5;
+        double tflops = flops/(ms*1e-3)/1e12;
+        // HBM traffic: Q + O (B*Hq*ql*D each) + K + V (B*Hk*kl*D each), bf16.
+        double bytes = 2.0 * (2.0*nQ + 2.0*nKV);
+        double gbps = bytes/(ms*1e-3)/1e9;
+
+        char cfg[64];
+        snprintf(cfg, sizeof(cfg),
+                 "B=%2d Hq=%2d Hk=%d q=%4d kv=%4d D=%3d causal=%d",
+                 B,Hq,Hk,ql,kl,D,causal);
+        printf("%-46s | %7.4f ms | %7.1f GB/s | %6.2f TFLOP/s\n",
+               cfg, ms, gbps, tflops);
+
+        cudaFree(dQ);cudaFree(dK);cudaFree(dV);cudaFree(dO);
+        delete[]tmp; cudaEventDestroy(s); cudaEventDestroy(e);
+    }
+}
+
+int main() {
+    const int configs[][7] = {
+        {1,2,1,64,128,64,0},     // tiny: B,Hq,Hk,q,kv,D,causal
+        {1,32,4,512,512,128,0},  // standard
+        {1,32,4,128,256,128,0},  // medium
+        {1,4,2,256,256,128,1},   // causal
+    };
+    int n_configs = sizeof(configs) / sizeof(configs[0]);
+
+    for (int ci = 0; ci < n_configs; ci++) {
+        int B=configs[ci][0], Hq=configs[ci][1], Hk=configs[ci][2];
+        int ql=configs[ci][3], kl=configs[ci][4], D=configs[ci][5];
+        int causal=configs[ci][6];
+        printf("=== B=%d Hq=%d Hk=%d q=%d kv=%d D=%d causal=%d ===\n",
+               B,Hq,Hk,ql,kl,D,causal);
+
+        size_t nQ = B*Hq*ql*D, nKV = B*Hk*kl*D;
+        float *hQ=new float[nQ], *hK=new float[nKV], *hV=new float[nKV];
+        for (size_t i=0;i<nQ;i++) hQ[i]=randf();
+        for (size_t i=0;i<nKV;i++){hK[i]=randf();hV[i]=randf();}
+
+        bf16 *dQ,*dK,*dV,*dO,*tmp;
+        cudaMalloc(&dQ,nQ*2); cudaMalloc(&dK,nKV*2);
+        cudaMalloc(&dV,nKV*2); cudaMalloc(&dO,nQ*2);
+        tmp=new bf16[max(nQ,nKV)];
+        for (size_t i=0;i<nQ;i++) tmp[i]=f2bf(hQ[i]);
+        cudaMemcpy(dQ,tmp,nQ*2,cudaMemcpyHostToDevice);
+        for (size_t i=0;i<nKV;i++) tmp[i]=f2bf(hK[i]);
+        cudaMemcpy(dK,tmp,nKV*2,cudaMemcpyHostToDevice);
+        for (size_t i=0;i<nKV;i++) tmp[i]=f2bf(hV[i]);
+        cudaMemcpy(dV,tmp,nKV*2,cudaMemcpyHostToDevice);
+
+        AttentionParams<bf16> p;
+        p.batch=B; p.q_head=Hq; p.kv_head=Hk; p.q_len=ql; p.kv_len=kl; p.head_dim=D;
+        p.use_mask=0; p.is_causal=causal; p.causal_offset=0;
+        p.scale=1.0f/sqrtf((float)D);
+        p.q=dQ; p.k=dK; p.v=dV; p.mask=nullptr; p.o=dO;
+
+        double t0=now_ms();
+        dispatch_prefill(p);
+        cudaDeviceSynchronize();
+        double kms=now_ms()-t0;
+        cudaError_t err=cudaGetLastError();
+        if (err!=cudaSuccess){printf("CUDA err: %s\n",cudaGetErrorString(err));return 1;}
+
+        bf16* hOut=new bf16[nQ];
+        cudaMemcpy(hOut,dO,nQ*2,cudaMemcpyDeviceToHost);
+
+        float* ref=new float[nQ];
+        cpu_attention_ref(hQ, hK, hV, nullptr, ref, B, Hq, Hk, ql, kl, D, causal, 0);
+
+        float max_err=0;
+        for (size_t i=0;i<nQ;i++) {
+            float d=fabsf(bf2f(hOut[i])-ref[i]);
+            if(d>max_err) max_err=d;
+        }
+        printf("kernel: %.3f ms  max_err: %.6e\n\n",kms,max_err);
+
+        cudaFree(dQ);cudaFree(dK);cudaFree(dV);cudaFree(dO);
+        delete[]hQ;delete[]hK;delete[]hV;delete[]hOut;delete[]ref;delete[]tmp;
+    }
+    printf("All tests passed!\n");
+    bench();
+    return 0;
+}

csrc/tests/gqa_decode_test.cu

@@ -1,130 +0,0 @@
-/*
-Pure-C test:
-nvcc -I csrc -arch=sm_89 -O3 \
-    --use_fast_math --ptxas-options=-O3 --extra-device-vectorization \
-    csrc/tests/gqa_decode_test.cu -o test && ./test
-*/
-
-#include <cstdio>
-#include <cstdlib>
-#include <cmath>
-#include <sys/time.h>
-#include "../kernels/gqa_decode_attn.cuh"
-
-static double now_ms() {
-    struct timeval tv;
-    gettimeofday(&tv, NULL);
-    return tv.tv_sec * 1000.0 + tv.tv_usec / 1000.0;
-}
-
-static void cpu_decode(const float* Q, const float* K, const float* V,
-                        const bool* mask, float* O,
-                        int B, int Hq, int Hk, int seq_len, int D) {
-    float scale = 1.0f / sqrtf((float)D);
-    int n_rep = Hq / Hk;
-    for (int b = 0; b < B; b++) {
-        for (int h = 0; h < Hq; h++) {
-            int kv_h = h / n_rep;
-            float mv = -INFINITY, sv = 0.0f;
-            float accum[256] = {0};
-            for (int s = 0; s < seq_len; s++) {
-                if (!mask[b * seq_len + s]) continue;
-                float dot = 0.0f;
-                for (int d = 0; d < D; d++)
-                    dot += Q[((b * Hq + h) * 1 + 0) * D + d]
-                         * K[((b * Hk + kv_h) * seq_len + s) * D + d];
-                dot *= scale;
-                float nm = fmaxf(mv, dot);
-                float al = expf(mv - nm);
-                float be = expf(dot - nm);
-                sv = sv * al + be;
-                for (int d = 0; d < D; d++)
-                    accum[d] = accum[d] * al
-                             + V[((b * Hk + kv_h) * seq_len + s) * D + d] * be;
-                mv = nm;
-            }
-            float inv = 1.0f / sv;
-            for (int d = 0; d < D; d++)
-                O[((b * Hq + h) * 1 + 0) * D + d] = accum[d] * inv;
-        }
-    }
-}
-
-static bf16 f2bf(float x) { return __float2bfloat16(x); }
-static float bf2f(bf16 x) { return __bfloat162float(x); }
-static float randf() { return (float)rand() / (float)RAND_MAX - 0.5f; }
-
-int main() {
-    const int configs[][5] = {
-        {1, 2, 1, 64, 32},    // B,Hq,Hk,seq_len,D
-        {1, 32, 4, 512, 128},
-        {1, 32, 4, 1024, 128},
-    };
-    int n_cfgs = sizeof(configs) / sizeof(configs[0]);
-
-    for (int ci = 0; ci < n_cfgs; ci++) {
-        int B = configs[ci][0], Hq = configs[ci][1], Hk = configs[ci][2];
-        int sl = configs[ci][3], D = configs[ci][4], gs = Hq / Hk;
-        printf("=== B=%d Hq=%d Hk=%d seq=%d D=%d gs=%d ===\n", B,Hq,Hk,sl,D,gs);
-
-        size_t nQ = B*Hq*1*D, nKV = B*Hk*sl*D;
-        float *hQ=new float[nQ], *hK=new float[nKV], *hV=new float[nKV];
-        for (size_t i=0;i<nQ;i++) hQ[i]=randf();
-        for (size_t i=0;i<nKV;i++){hK[i]=randf();hV[i]=randf();}
-
-        bool* hMask=new bool[B*sl];
-        for (int i=0;i<B*sl;i++) hMask[i]=true;
-
-        bf16 *dQ,*dK,*dV,*dO,*tmp;
-        bool* dMask;
-        cudaMalloc(&dQ,nQ*2); cudaMalloc(&dK,nKV*2);
-        cudaMalloc(&dV,nKV*2); cudaMalloc(&dO,nQ*2);
-        cudaMalloc(&dMask,B*sl);
-
-        tmp=new bf16[max(nQ,nKV)];
-        for (size_t i=0;i<nQ;i++) tmp[i]=f2bf(hQ[i]);
-        cudaMemcpy(dQ,tmp,nQ*2,cudaMemcpyHostToDevice);
-        for (size_t i=0;i<nKV;i++) tmp[i]=f2bf(hK[i]);
-        cudaMemcpy(dK,tmp,nKV*2,cudaMemcpyHostToDevice);
-        for (size_t i=0;i<nKV;i++) tmp[i]=f2bf(hV[i]);
-        cudaMemcpy(dV,tmp,nKV*2,cudaMemcpyHostToDevice);
-        cudaMemcpy(dMask,hMask,B*sl,cudaMemcpyHostToDevice);
-
-        GQAParams p;
-        p.batch=B; p.q_head=Hq; p.kv_head=Hk; p.q_len=1; p.kv_len=sl; p.head_dim=D;
-        p.use_mask=1; p.is_causal=0; p.causal_offset=0;
-        p.scale=1.0f/sqrtf((float)D);
-        p.q=dQ; p.k=dK; p.v=dV; p.mask=dMask; p.o=dO;
-
-        size_t smem=DC_CHUNK*D*sizeof(bf16);
-        dim3 block(32, gs);
-        dim3 grid(B*Hk);
-        printf("grid=(%d,1,1) block=(%d,%d,1) smem=%zu\n",
-               grid.x, block.x, block.y, smem);
-
-        double t0=now_ms();
-        gqa_decode_attn_kernel<<<grid,block,smem>>>(p);
-        cudaDeviceSynchronize();
-        double kms=now_ms()-t0;
-        cudaError_t err=cudaGetLastError();
-        if (err!=cudaSuccess){printf("CUDA err: %s\n",cudaGetErrorString(err));return 1;}
-
-        bf16* hOut=new bf16[nQ];
-        cudaMemcpy(hOut,dO,nQ*2,cudaMemcpyDeviceToHost);
-
-        float* ref=new float[nQ];
-        cpu_decode(hQ,hK,hV,hMask,ref,B,Hq,Hk,sl,D);
-
-        float max_err=0;
-        for (size_t i=0;i<nQ;i++){
-            float d=fabsf(bf2f(hOut[i])-ref[i]);
-            if(d>max_err) max_err=d;
-        }
-        printf("kernel: %.3f ms  max_err: %.6e\n\n",kms,max_err);
-
-        cudaFree(dQ);cudaFree(dK);cudaFree(dV);cudaFree(dO);cudaFree(dMask);
-        delete[]hQ;delete[]hK;delete[]hV;delete[]hMask;delete[]hOut;delete[]ref;delete[]tmp;
-    }
-    printf("All tests passed!\n");
-    return 0;
-}

csrc/tests/gqa_prefill_test.cu

@@ -1,133 +0,0 @@
-/*
-Pure-C test:
-nvcc -I csrc -arch=sm_89 -O3 \
-    --use_fast_math --ptxas-options=-O3 --extra-device-vectorization \
-    csrc/tests/gqa_prefill_test.cu -o test && ./test
-*/
-
-#include <cstdio>
-#include <cstdlib>
-#include <cmath>
-#include <sys/time.h>
-#include "../kernels/gqa_prefill_attn.cuh"
-
-static double now_ms() {
-    struct timeval tv;
-    gettimeofday(&tv, NULL);
-    return tv.tv_sec * 1000.0 + tv.tv_usec / 1000.0;
-}
-
-static void cpu_attention(const float* Q, const float* K, const float* V, float* O,
-                          int B, int Hq, int Hk, int q_len, int kv_len, int D,
-                          int is_causal, int causal_off) {
-    float scale = 1.0f / sqrtf((float)D);
-    int n_rep = Hq / Hk;
-    for (int b = 0; b < B; b++) {
-        for (int h = 0; h < Hq; h++) {
-            for (int qi = 0; qi < q_len; qi++) {
-                int kv_h = h / n_rep;
-                float mv = -INFINITY, sv = 0.0f;
-                float accum[256] = {0};
-                int lim = is_causal ? min(kv_len, qi + causal_off + 1) : kv_len;
-                for (int kj = 0; kj < lim; kj++) {
-                    float dot = 0.0f;
-                    for (int d = 0; d < D; d++)
-                        dot += Q[((b*Hq + h)*q_len + qi)*D + d]
-                             * K[((b*Hk + kv_h)*kv_len + kj)*D + d];
-                    dot *= scale;
-                    float nm = fmaxf(mv, dot);
-                    float al = expf(mv - nm);
-                    float be = expf(dot - nm);
-                    sv = sv * al + be;
-                    for (int d = 0; d < D; d++)
-                        accum[d] = accum[d] * al
-                                 + V[((b*Hk + kv_h)*kv_len + kj)*D + d] * be;
-                    mv = nm;
-                }
-                float inv = 1.0f / sv;
-                for (int d = 0; d < D; d++)
-                    O[((b*Hq + h)*q_len + qi)*D + d] = accum[d] * inv;
-            }
-        }
-    }
-}
-
-static __nv_bfloat16 f2bf(float x) { return __float2bfloat16(x); }
-static float bf2f(__nv_bfloat16 x) { return __bfloat162float(x); }
-static float randf() { return (float)rand() / (float)RAND_MAX - 0.5f; }
-
-int main() {
-    const int configs[][7] = {
-        {1,2,1,64,128,64,0},     // tiny: B,Hq,Hk,q,kv,D,causal
-        {1,32,4,512,512,128,0},  // standard
-        {1,32,4,128,256,128,0},  // medium
-        {1,4,2,256,256,128,1},   // causal
-    };
-    int n_configs = sizeof(configs) / sizeof(configs[0]);
-
-    for (int ci = 0; ci < n_configs; ci++) {
-        int B=configs[ci][0], Hq=configs[ci][1], Hk=configs[ci][2];
-        int ql=configs[ci][3], kl=configs[ci][4], D=configs[ci][5];
-        int causal=configs[ci][6];
-        printf("=== B=%d Hq=%d Hk=%d q=%d kv=%d D=%d causal=%d ===\n",
-               B,Hq,Hk,ql,kl,D,causal);
-
-        size_t nQ = B*Hq*ql*D, nKV = B*Hk*kl*D;
-        float *hQ=new float[nQ], *hK=new float[nKV], *hV=new float[nKV];
-        for (size_t i=0;i<nQ;i++) hQ[i]=randf();
-        for (size_t i=0;i<nKV;i++){hK[i]=randf();hV[i]=randf();}
-
-        bf16 *dQ,*dK,*dV,*dO,*tmp;
-        cudaMalloc(&dQ,nQ*2); cudaMalloc(&dK,nKV*2);
-        cudaMalloc(&dV,nKV*2); cudaMalloc(&dO,nQ*2);
-        tmp=new bf16[max(nQ,nKV)];
-        for (size_t i=0;i<nQ;i++) tmp[i]=f2bf(hQ[i]);
-        cudaMemcpy(dQ,tmp,nQ*2,cudaMemcpyHostToDevice);
-        for (size_t i=0;i<nKV;i++) tmp[i]=f2bf(hK[i]);
-        cudaMemcpy(dK,tmp,nKV*2,cudaMemcpyHostToDevice);
-        for (size_t i=0;i<nKV;i++) tmp[i]=f2bf(hV[i]);
-        cudaMemcpy(dV,tmp,nKV*2,cudaMemcpyHostToDevice);
-
-        GQAParams p;
-        p.batch=B; p.q_head=Hq; p.kv_head=Hk; p.q_len=ql; p.kv_len=kl; p.head_dim=D;
-        p.use_mask=0; p.is_causal=causal; p.causal_offset=0;
-        p.scale=1.0f/sqrtf((float)D);
-        p.q=dQ; p.k=dK; p.v=dV; p.mask=nullptr; p.o=dO;
-
-        constexpr int G=8, ROWS=32, P_BC=32;
-        dim3 grid((ql+ROWS-1)/ROWS, Hq, B);
-        dim3 block(G, ROWS, 1);
-        size_t smem=2*P_BC*D*sizeof(bf16);
-        printf("grid=(%d,%d,%d) block=(%d,%d,%d) smem=%zu\n",
-               grid.x,grid.y,grid.z, block.x,block.y,block.z, smem);
-
-        double t0=now_ms();
-        switch (D) {
-            case 64:  gqa_prefill_attn_kernel_t<64, G,ROWS,P_BC><<<grid,block,smem>>>(p); break;
-            case 128: gqa_prefill_attn_kernel_t<128,G,ROWS,P_BC><<<grid,block,smem>>>(p); break;
-            default: printf("unsupported D=%d\n",D); return 1;
-        }
-        cudaDeviceSynchronize();
-        double kms=now_ms()-t0;
-        cudaError_t err=cudaGetLastError();
-        if (err!=cudaSuccess){printf("CUDA err: %s\n",cudaGetErrorString(err));return 1;}
-
-        bf16* hOut=new bf16[nQ];
-        cudaMemcpy(hOut,dO,nQ*2,cudaMemcpyDeviceToHost);
-
-        float* ref=new float[nQ];
-        cpu_attention(hQ,hK,hV,ref,B,Hq,Hk,ql,kl,D,causal,0);
-
-        float max_err=0;
-        for (size_t i=0;i<nQ;i++) {
-            float d=fabsf(bf2f(hOut[i])-ref[i]);
-            if(d>max_err) max_err=d;
-        }
-        printf("kernel: %.3f ms  max_err: %.6e\n\n",kms,max_err);
-
-        cudaFree(dQ);cudaFree(dK);cudaFree(dV);cudaFree(dO);
-        delete[]hQ;delete[]hK;delete[]hV;delete[]hOut;delete[]ref;delete[]tmp;
-    }
-    printf("All tests passed!\n");
-    return 0;
-}

csrc/tests/test_utils.cuh

@@ -0,0 +1,154 @@
+#pragma once
+
+#include <cstdio>
+#include <cstdlib>
+#include <cmath>
+#include <chrono>
+#include <cuda_bf16.h>
+
+using bf16 = __nv_bfloat16;
+
+inline bf16 f2bf(float x) { return __float2bfloat16(x); }
+inline float bf2f(bf16 x) { return __bfloat162float(x); }
+
+inline float randf() { return (float)rand() / (float)RAND_MAX - 0.5f; }
+
+inline double now_ms() {
+    using namespace std::chrono;
+    return duration_cast<milliseconds>(steady_clock::now().time_since_epoch()).count();
+}
+
+inline int compute_num_splits(int base_blocks, int tiles_total) {
+    int sm_count = 0;
+    cudaDeviceGetAttribute(&sm_count, cudaDevAttrMultiProcessorCount, 0);
+    int n = (2 * sm_count + base_blocks - 1) / base_blocks;
+    if (n > tiles_total) n = tiles_total;
+    if (n > 32) n = 32;
+    if (n < 1) n = 1;
+    return n;
+}
+
+#define CUDA_CHECK(call) \
+    do { \
+        cudaError_t _e = (call); \
+        if (_e != cudaSuccess) { \
+            printf("CUDA error %s at %s:%d\n", cudaGetErrorString(_e), __FILE__, __LINE__); \
+            exit(1); \
+        } \
+    } while (0)
+
+struct BenchResult {
+    float ms;
+    double gbps;
+    double tflops;
+};
+
+template <typename Fn>
+BenchResult bench_kernel(Fn launch, int warmup, int iters,
+                         double flops, double bytes) {
+    for (int i = 0; i < warmup; i++) launch();
+    cudaDeviceSynchronize();
+    cudaError_t err = cudaGetLastError();
+    if (err != cudaSuccess) {
+        printf("CUDA error before bench: %s\n", cudaGetErrorString(err));
+        return {0, 0, 0};
+    }
+
+    cudaEvent_t s, e;
+    cudaEventCreate(&s); cudaEventCreate(&e);
+    cudaEventRecord(s);
+    for (int i = 0; i < iters; i++) launch();
+    cudaEventRecord(e); cudaEventSynchronize(e);
+    float ms = 0; cudaEventElapsedTime(&ms, s, e); ms /= iters;
+    cudaEventDestroy(s); cudaEventDestroy(e);
+
+    return {ms, bytes / (ms * 1e-3) / 1e9, flops / (ms * 1e-3) / 1e12};
+}
+
+inline void print_bench_header() {
+    printf("%-46s | %10s | %10s | %10s\n",
+           "config", "latency", "bandwidth", "throughput");
+    printf("---------------------------------------------------------------"
+           "----------------------------\n");
+}
+
+inline void print_bench_row(const char* cfg, const BenchResult& r) {
+    printf("%-46s | %7.4f ms | %7.1f GB/s | %6.2f TFLOP/s\n",
+           cfg, r.ms, r.gbps, r.tflops);
+}
+
+template <int... Ds>
+struct _HeadSwitch;
+
+template <int D>
+struct _HeadSwitch<D> {
+    template <typename Fn>
+    static void call(int hd, Fn&& fn) { if (hd == D) fn.template operator()<D>(); }
+};
+
+template <int D, int... Rest>
+struct _HeadSwitch<D, Rest...> {
+    template <typename Fn>
+    static void call(int hd, Fn&& fn) {
+        if (hd == D) fn.template operator()<D>();
+        else _HeadSwitch<Rest...>::call(hd, fn);
+    }
+};
+
+// Default set: 32, 64, 128, 256
+template <typename Fn>
+void dispatch_by_head_dim(int head_dim, Fn&& fn) {
+    _HeadSwitch<32, 64, 128, 256>::call(head_dim, fn);
+}
+
+// Generic CPU reference for multi-query / grouped-query attention.
+// Tensor shapes (all float*):
+//   Q : [B, Hq, q_len, D]
+//   K : [B, Hk, kv_len, D]
+//   V : [B, Hk, kv_len, D]
+//   O : [B, Hq, q_len, D]
+// mask: if q_len == 1, shape is [B, kv_len]; otherwise mask is not supported.
+static void cpu_attention_ref(
+    const float* Q, const float* K, const float* V, const bool* mask,
+    float* O, int B, int Hq, int Hk, int q_len, int kv_len, int D,
+    int is_causal, int causal_offset
+) {
+    float scale = 1.0f / sqrtf((float)D);
+    int n_rep = Hq / Hk;
+    for (int b = 0; b < B; b++) {
+        for (int h = 0; h < Hq; h++) {
+            int kv_h = h / n_rep;
+            for (int qi = 0; qi < q_len; qi++) {
+                float mv = -INFINITY, sv = 0.0f;
+                float accum[256] = {0.0f};
+                int lim = kv_len;
+                if (is_causal) {
+                    int c = qi + causal_offset + 1;
+                    lim = (c < kv_len) ? c : kv_len;
+                }
+                for (int kj = 0; kj < lim; kj++) {
+                    if (mask != nullptr && q_len == 1) {
+                        if (!mask[b * kv_len + kj]) continue;
+                    }
+                    float dot = 0.0f;
+                    size_t q_idx = ((size_t)b * Hq + h) * q_len + qi;
+                    size_t kv_idx = ((size_t)b * Hk + kv_h) * kv_len + kj;
+                    for (int d = 0; d < D; d++)
+                        dot += Q[q_idx * D + d] * K[kv_idx * D + d];
+                    dot *= scale;
+                    float nm = fmaxf(mv, dot);
+                    float a = expf(mv - nm);
+                    float b_exp = expf(dot - nm);
+                    sv = sv * a + b_exp;
+                    for (int d = 0; d < D; d++)
+                        accum[d] = accum[d] * a + V[kv_idx * D + d] * b_exp;
+                    mv = nm;
+                }
+                float inv = 1.0f / sv;
+                size_t o_idx = ((size_t)b * Hq + h) * q_len + qi;
+                for (int d = 0; d < D; d++)
+                    O[o_idx * D + d] = accum[d] * inv;
+            }
+        }
+    }
+}