FP8 GEMM DeepSeek V4 Flash (0731)

"""FP8 e4m3 GEMM via a CUTLASS SM90 wgmma (QGMMA) tensor-core kernel.

Computes y = (x @ w.T) * weight_scale as bf16, with real fp8 x fp8 MMA
(fp8 e4m3 inputs, fp32 accumulate) and a fused per-output-channel scale in
the epilogue. Built with torch.utils.cpp_extension.load_inline (cached).
"""
import os
import datetime

# Python 3.10 compat: the harness benchmark imports `datetime.UTC` (3.11+).
# solution.py is imported before that import in benchmark.py, so patch it here.
if not hasattr(datetime, "UTC"):
    datetime.UTC = datetime.timezone.utc

import torch
import torch.nn as nn
from torch.utils.cpp_extension import load_inline

try:
    import pybind11
    _PYBIND_INC = [pybind11.get_include()]
except ImportError:
    _PYBIND_INC = []

os.environ.setdefault("TORCH_CUDA_ARCH_LIST", "9.0a")

# Locate CUTLASS headers: vendored copy next to this file first, else /tmp/cutlass.
_SOL_DIR = os.path.dirname(os.path.abspath(__file__))
_CUTLASS_INCLUDE = os.environ.get("CUTLASS_INCLUDE", "")
if not _CUTLASS_INCLUDE:
    for cand in (os.path.join(_SOL_DIR, "cutlass_local", "include"),
                 "/tmp/cutlass/include"):
        if os.path.isdir(os.path.join(cand, "cutlass")):
            _CUTLASS_INCLUDE = cand
            break
_CUTLASS_UTIL = os.path.join(os.path.dirname(_CUTLASS_INCLUDE), "tools", "util", "include")

E4M3_MAX = 448.0

CPP_SRC = r"""
#include <torch/extension.h>
#include <ATen/cuda/CUDAContext.h>
#include <cuda_runtime.h>
#include <cutlass/cutlass.h>
#include <cutlass/numeric_types.h>
#include <cutlass/gemm/device/gemm_universal_adapter.h>
#include <cutlass/gemm/kernel/gemm_universal.hpp>
#include <cutlass/gemm/collective/collective_builder.hpp>
#include <cutlass/epilogue/collective/collective_builder.hpp>
#include <cutlass/epilogue/fusion/operations.hpp>

using namespace cute;

using ElementA = cutlass::float_e4m3_t;
using ElementB = cutlass::float_e4m3_t;
using ElementAcc = float;
using ElementOutput = cutlass::bfloat16_t;
using ElementScale = float;

using LayoutA = cutlass::layout::RowMajor;
using LayoutB = cutlass::layout::ColumnMajor;   // B is (N,K) row-major weight
using LayoutD = cutlass::layout::RowMajor;

using FusionOp = cutlass::epilogue::fusion::PerColLinCombPerColBiasEltAct<
    cutlass::epilogue::thread::Identity,
    ElementOutput, float, float, ElementOutput, ElementScale>;

template <class TileShape, class ClusterShape, class Schedule, int EpiM, class EpiSchedule>
using MakeEpilogue = typename cutlass::epilogue::collective::CollectiveBuilder<
    cutlass::arch::Sm90, cutlass::arch::OpClassTensorOp,
    TileShape, ClusterShape,
    Shape<cute::Int<EpiM>, _128>,
    ElementAcc, float,
    ElementOutput, LayoutD, 16 / sizeof(ElementOutput),
    ElementOutput, LayoutD, 16 / sizeof(ElementOutput),
    EpiSchedule,
    FusionOp>::CollectiveOp;

template <class TileShape, class ClusterShape, class Schedule, class EpilogueOp, int STAGES>
using MakeMainloop = typename cutlass::gemm::collective::CollectiveBuilder<
    cutlass::arch::Sm90, cutlass::arch::OpClassTensorOp,
    ElementA, LayoutA, 16,
    ElementB, LayoutB, 16,
    ElementAcc,
    TileShape, ClusterShape,
    typename std::conditional<STAGES == 0,
        cutlass::gemm::collective::StageCountAutoCarveout<sizeof(typename EpilogueOp::SharedStorage)>,
        cutlass::gemm::collective::StageCount<STAGES>>::type,
    Schedule>::CollectiveOp;

#define DEFINE_GEMM(NAME, TS_M, TS_N, TS_K, CL_M, CL_N, CL_K, SCH, EPISCH, STAGES) \
  using NAME##_Epi = MakeEpilogue<Shape<_##TS_M, _##TS_N, _##TS_K>, Shape<_##CL_M, _##CL_N, _##CL_K>, SCH, TS_M, EPISCH>; \
  using NAME##_Main = MakeMainloop<Shape<_##TS_M, _##TS_N, _##TS_K>, Shape<_##CL_M, _##CL_N, _##CL_K>, SCH, NAME##_Epi, STAGES>; \
  using NAME = cutlass::gemm::device::GemmUniversalAdapter< \
      cutlass::gemm::kernel::GemmUniversal<Shape<int,int,int,int>, NAME##_Main, NAME##_Epi>>;

// Compute-bound shapes: 128x128x128 cooperative, 5 pipeline stages
DEFINE_GEMM(G128N128K128, 128, 128, 128, 1, 1, 1, cutlass::gemm::KernelTmaWarpSpecializedCooperative, cutlass::epilogue::TmaWarpSpecializedCooperative, 5)
// Skinny M: 64x128x128 warp-specialized
DEFINE_GEMM(G64N128K128, 64, 128, 128, 1, 1, 1, cutlass::gemm::collective::KernelScheduleAuto, cutlass::epilogue::TmaWarpSpecialized, 0)

template <class Gemm>
void run_gemm(torch::Tensor x, torch::Tensor w, torch::Tensor scale, torch::Tensor y) {
    int M = x.size(0), K = x.size(1), N = w.size(0);
    auto stream = at::cuda::getCurrentCUDAStream();
    typename Gemm::Arguments args;
    args.mode = cutlass::gemm::GemmUniversalMode::kGemm;
    args.problem_shape = {M, N, K, 1};
    args.mainloop = {
        reinterpret_cast<ElementA const*>(x.data_ptr()),
        cute::Stride<int64_t, cute::Int<1>, int64_t>{K, cute::Int<1>{}, int64_t(0)},
        reinterpret_cast<ElementB const*>(w.data_ptr()),
        cute::Stride<int64_t, cute::Int<1>, int64_t>{K, cute::Int<1>{}, int64_t(0)}
    };
    args.epilogue.thread = {
        1.0f, 0.0f,
        reinterpret_cast<ElementScale const*>(scale.data_ptr()), nullptr,
        cute::Stride<_0, bool, int64_t>{_0{}, bool(1), 0},
        cute::Stride<_0, bool, int64_t>{_0{}, bool(1), 0},
        nullptr, cute::Stride<_0, _1, int64_t>{_0{}, _1{}, 0},
        {}
    };
    args.epilogue.ptr_C = nullptr;
    args.epilogue.dC = cute::Stride<int64_t, cute::Int<1>, int64_t>{N, cute::Int<1>{}, int64_t(0)};
    args.epilogue.ptr_D = reinterpret_cast<ElementOutput const*>(y.data_ptr());
    args.epilogue.dD = cute::Stride<int64_t, cute::Int<1>, int64_t>{N, cute::Int<1>{}, int64_t(0)};
    int dev_id = at::cuda::current_device();
    args.hw_info = {dev_id, at::cuda::getCurrentDeviceProperties()->multiProcessorCount};
    args.scheduler = {};
    Gemm gemm;
    auto status = gemm.run(args, stream);
    TORCH_CHECK(status == cutlass::Status::kSuccess, "CUTLASS GEMM failed");
}

torch::Tensor fp8_gemm(torch::Tensor x, torch::Tensor w, torch::Tensor scale, int64_t config) {
    int M = x.size(0), N = w.size(0);
    auto y = torch::empty({M, N}, x.options().dtype(torch::kBFloat16));
    if (M < 128) {
        run_gemm<G64N128K128>(x, w, scale, y);
    } else {
        run_gemm<G128N128K128>(x, w, scale, y);
    }
    return y;
}
"""

CU_SRC = r"""
#include <torch/extension.h>
torch::Tensor fp8_gemm(torch::Tensor x, torch::Tensor w, torch::Tensor scale, int64_t config);
PYBIND11_MODULE(TORCH_EXTENSION_NAME, m) {
    m.def("fp8_gemm", &fp8_gemm, "fp8 gemm");
}
"""

# Triton fallback (used only if the CUTLASS extension cannot be built/loaded).
try:
    import triton
    import triton.language as tl

    @triton.jit
    def _triton_gemm(
        x_ptr, w_ptr, s_ptr, y_ptr,
        M, N, K,
        stride_xm, stride_xk,
        stride_wk, stride_wn,
        stride_ym, stride_yn,
        BLOCK_M: tl.constexpr, BLOCK_N: tl.constexpr, BLOCK_K: tl.constexpr,
        GROUP_M: tl.constexpr,
    ):
        pid = tl.program_id(0)
        num_pid_m = tl.cdiv(M, BLOCK_M)
        num_pid_n = tl.cdiv(N, BLOCK_N)
        num_pid_in_group = GROUP_M * num_pid_n
        group_id = pid // num_pid_in_group
        first_pid_m = group_id * GROUP_M
        group_size_m = min(num_pid_m - first_pid_m, GROUP_M)
        pid_m = first_pid_m + ((pid % num_pid_in_group) % group_size_m)
        pid_n = (pid % num_pid_in_group) // group_size_m
        offs_m = pid_m * BLOCK_M + tl.arange(0, BLOCK_M)
        offs_n = pid_n * BLOCK_N + tl.arange(0, BLOCK_N)
        offs_k = tl.arange(0, BLOCK_K)
        x_ptrs = x_ptr + offs_m[:, None] * stride_xm + offs_k[None, :] * stride_xk
        w_ptrs = w_ptr + offs_k[:, None] * stride_wk + offs_n[None, :] * stride_wn
        acc = tl.zeros((BLOCK_M, BLOCK_N), dtype=tl.float32)
        for k0 in tl.range(0, tl.cdiv(K, BLOCK_K)):
            kmask = (k0 * BLOCK_K + offs_k) < K
            x = tl.load(x_ptrs, mask=(offs_m[:, None] < M) & kmask[None, :], other=0.0)
            w = tl.load(w_ptrs, mask=(offs_n[None, :] < N) & kmask[:, None], other=0.0)
            acc = tl.dot(x, w, acc)
            x_ptrs += BLOCK_K * stride_xk
            w_ptrs += BLOCK_K * stride_wk
        s = tl.load(s_ptr + offs_n, mask=offs_n < N, other=0.0)
        acc = acc * s[None, :]
        y_ptrs = y_ptr + offs_m[:, None] * stride_ym + offs_n[None, :] * stride_yn
        tl.store(y_ptrs, acc.to(tl.bfloat16), mask=(offs_m[:, None] < M) & (offs_n[None, :] < N))

    def _triton_fp8_gemm(x, w, scale, config=0):
        M, K = x.shape
        N = w.shape[0]
        y = torch.empty((M, N), dtype=torch.bfloat16, device=x.device)
        BM, BN, BK = 128, 128, 64
        grid = (triton.cdiv(M, BM) * triton.cdiv(N, BN),)
        _triton_gemm[grid](
            x, w, scale, y, M, N, K,
            x.stride(0), x.stride(1), w.stride(1), w.stride(0),
            y.stride(0), y.stride(1),
            BLOCK_M=BM, BLOCK_N=BN, BLOCK_K=BK, GROUP_M=8,
            num_warps=8, num_stages=3,
        )
        return y
except Exception:
    triton = None
    _triton_fp8_gemm = None

_ext = None
_use_triton = False


def _get_ext():
    global _ext, _use_triton
    if _ext is None and not _use_triton:
        try:
            _ext = load_inline(
                name="cutlass_fp8_gemm_solution",
                cpp_sources=[CU_SRC],
                cuda_sources=[CPP_SRC],
                extra_cuda_cflags=["-O3", "--std=c++17", "-DCUTLASS_ARCH_MMA_SM90_SUPPORTED",
                                   "-I" + _CUTLASS_INCLUDE, "-I" + _CUTLASS_UTIL],
                extra_include_paths=_PYBIND_INC,
                verbose=False,
            )
        except Exception:
            if _triton_fp8_gemm is None:
                raise
            _use_triton = True
    return _ext


class Model(nn.Module):
    """y = ((x @ w.T) * weight_scale).to(bf16) with real fp8 tensor-core MMA."""

    def __init__(self, M: int, N: int, K: int):
        super().__init__()
        self.M, self.N, self.K = M, N, K
        w = torch.empty(N, K, dtype=torch.bfloat16)
        nn.init.normal_(w, std=0.02)
        s = (w.float().abs().amax(dim=1, keepdim=True) / E4M3_MAX).clamp(min=1e-12)
        w_fp8 = (w.float() / s).to(torch.float8_e4m3fn)
        self.register_buffer("weight", w_fp8)                         # (N, K) fp8
        self.register_buffer("weight_scale", s.squeeze(1).to(torch.float32))  # (N,)
        self._wp = None
        self._wver = -1
        # padded-x cache (only used when K is not a multiple of 16).
        self._xcache = {}

    def _pad_x(self, x: torch.Tensor, Kp: int):
        """Pad x's K dim to a multiple of 16 (TMA stride requirement).

        Keyed on Python object id while holding a strong reference to x, so an
        id can never be recycled to a different tensor while cached; the
        version counter catches in-place mutation of the same tensor.
        """
        key = id(x)
        ent = self._xcache.get(key)
        if ent is not None:
            xp, xref, ver = ent
            if xref is x and ver == x._version:
                return xp
        xp = torch.zeros(x.shape[0], Kp, dtype=torch.float8_e4m3fn, device=x.device)
        xp[:, : x.shape[1]] = x
        self._xcache[key] = (xp, x, x._version)
        return xp

    def forward(self, x: torch.Tensor) -> torch.Tensor:
        ext = _get_ext()
        M, K = self.M, self.K
        N = self.N
        Kp = ((K + 15) // 16) * 16
        w = self.weight
        if self._wp is None or self._wver != w._version:
            if Kp != K:
                wp = torch.zeros(N, Kp, dtype=torch.float8_e4m3fn, device=w.device)
                wp[:, :K] = w
            else:
                wp = w
            self._wp = wp
            self._wver = w._version
        w_use = self._wp
        if Kp != K:
            x_use = self._pad_x(x, Kp)
        else:
            x_use = x
        if _use_triton:
            return _triton_fp8_gemm(x_use, w_use, self.weight_scale)
        return ext.fp8_gemm(x_use, w_use, self.weight_scale, 0)


M = 4096
N = 4096
K = 4096


def get_inputs():
    x = (torch.rand(M, K) * 8 - 4).to(torch.float8_e4m3fn)
    return [x]


def get_init_inputs():
    return [M, N, K]