submission 399585

#!POPCORN leaderboard nvfp4_group_gemm
#!POPCORN gpu NVIDIA

import torch
from task import input_t, output_t
from torch.utils.cpp_extension import load_inline

CUDA_SRC = r"""
#include <cudaTypedefs.h>
#include <cuda_fp16.h>

#include <torch/library.h>
#include <ATen/core/Tensor.h>

constexpr int WARP_SIZE = 32;
constexpr int MMA_K = 64;  // 32 bytes

__device__ __host__
constexpr int cdiv(int a, int b) { return (a + b - 1) / b; }

// https://github.com/NVIDIA/cutlass/blob/v4.3.2/include/cute/arch/copy_sm90_desc.hpp#L193-L197
constexpr uint64_t EVICT_NORMAL = 0x1000000000000000;
constexpr uint64_t EVICT_FIRST = 0x12F0000000000000;
constexpr uint64_t EVICT_LAST = 0x14F0000000000000;

__device__ inline
constexpr uint64_t desc_encode(uint64_t x) { return (x & 0x3'FFFFULL) >> 4ULL; };

// https://github.com/NVIDIA/cutlass/blob/v4.2.1/include/cute/arch/cluster_sm90.hpp#L180
__device__
uint32_t elect_sync() {
  uint32_t pred = 0;
  asm volatile(
    "{\n\t"
    ".reg .pred %%px;\n\t"
    "elect.sync _|%%px, %1;\n\t"
    "@%%px mov.s32 %0, 1;\n\t"
    "}"
    : "+r"(pred)
    : "r"(0xFFFFFFFF)
  );
  return pred;
}

__device__ inline
void mbarrier_init(int mbar_addr, int count) {
  asm volatile("mbarrier.init.shared::cta.b64 [%0], %1;" :: "r"(mbar_addr), "r"(count));
}

// NOTE: using .shared::cluster
__device__ inline
void mbarrier_arrive_expect_tx(int mbar_addr, int size) {
  asm volatile("mbarrier.arrive.expect_tx.release.cta.shared::cluster.b64 _, [%0], %1;" :: "r"(mbar_addr), "r"(size) : "memory");
}

// https://github.com/NVIDIA/cutlass/blob/v4.2.1/include/cutlass/arch/barrier.h#L408
__device__
void mbarrier_wait(int mbar_addr, int phase) {
  uint32_t ticks = 0x989680;  // this is optional
  asm volatile(
    "{\n\t"
    ".reg .pred P1;\n\t"
    "LAB_WAIT:\n\t"
    "mbarrier.try_wait.parity.acquire.cta.shared::cta.b64 P1, [%0], %1, %2;\n\t"
    "@!P1 bra.uni LAB_WAIT;\n\t"
    "}"
    :: "r"(mbar_addr), "r"(phase), "r"(ticks)
  );
}

__device__ inline
void tma_prefetch(const void *src, int size, uint64_t cache_policy) {
  asm volatile("cp.async.bulk.prefetch.L2.global.L2::cache_hint [%0], %1, %2;"
              :: "l"(src), "r"(size), "l"(cache_policy) : "memory");
}

__device__ inline
void tma_1d_prefetch(const void *tmap_ptr, int x, uint64_t cache_policy) {
  asm volatile("cp.async.bulk.prefetch.tensor.1d.L2.global.L2::cache_hint [%0, {%1}], %2;"
              :: "l"(tmap_ptr), "r"(x), "l"(cache_policy) : "memory");
}

__device__ inline
void tma_2d_prefetch(const void *tmap_ptr, int x, int y, uint64_t cache_policy) {
  asm volatile("cp.async.bulk.prefetch.tensor.2d.L2.global.L2::cache_hint [%0, {%1, %2}], %3;"
              :: "l"(tmap_ptr), "r"(x), "r"(y), "l"(cache_policy) : "memory");
}

__device__ inline
void tma_3d_prefetch(const void *tmap_ptr, int x, int y, int z, uint64_t cache_policy) {
  asm volatile("cp.async.bulk.prefetch.tensor.3d.L2.global.L2::cache_hint [%0, {%1, %2, %3}], %4;"
              :: "l"(tmap_ptr), "r"(x), "r"(y), "r"(z), "l"(cache_policy) : "memory");
}

__device__ inline
void tma_gmem2smem(int dst, const void *src, int size, int mbar_addr, uint64_t cache_policy) {
  asm volatile("cp.async.bulk.shared::cta.global.mbarrier::complete_tx::bytes.L2::cache_hint [%0], [%1], %2, [%3], %4;"
              :: "r"(dst), "l"(src), "r"(size), "r"(mbar_addr), "l"(cache_policy));
}

template <int CTA_GROUP = 1>
__device__ inline
void tma_1d_gmem2smem(int dst, const void *tmap_ptr, int x, int mbar_addr, uint64_t cache_policy) {
  asm volatile("cp.async.bulk.tensor.1d.shared::cta.global.mbarrier::complete_tx::bytes.cta_group::%5.L2::cache_hint "
              "[%0], [%1, {%2}], [%3], %4;"
              :: "r"(dst), "l"(tmap_ptr), "r"(x), "r"(mbar_addr), "l"(cache_policy), "n"(CTA_GROUP)
              : "memory");
}

template <int CTA_GROUP = 1>
__device__ inline
void tma_1d_gmem2smem_mcast(int dst, const void *tmap_ptr, int x, int mbar_addr, int16_t cta_mask, uint64_t cache_policy) {
  asm volatile("cp.async.bulk.tensor.1d.shared::cluster.global.mbarrier::complete_tx::bytes.multicast::cluster.cta_group::%6.L2::cache_hint "
              "[%0], [%1, {%2}], [%3], %4, %5;"
              :: "r"(dst), "l"(tmap_ptr), "r"(x), "r"(mbar_addr), "h"(cta_mask), "l"(cache_policy), "n"(CTA_GROUP)
              : "memory");
}

template <int CTA_GROUP = 1>
__device__ inline
void tma_2d_gmem2smem(int dst, const void *tmap_ptr, int x, int y, int mbar_addr, uint64_t cache_policy) {
  asm volatile("cp.async.bulk.tensor.2d.shared::cta.global.mbarrier::complete_tx::bytes.cta_group::%6.L2::cache_hint "
              "[%0], [%1, {%2, %3}], [%4], %5;"
              :: "r"(dst), "l"(tmap_ptr), "r"(x), "r"(y), "r"(mbar_addr), "l"(cache_policy), "n"(CTA_GROUP)
              : "memory");
}

template <int CTA_GROUP = 1>
__device__ inline
void tma_2d_gmem2smem_mcast(int dst, const void *tmap_ptr, int x, int y, int mbar_addr, int16_t cta_mask, uint64_t cache_policy) {
  asm volatile("cp.async.bulk.tensor.2d.shared::cluster.global.mbarrier::complete_tx::bytes.multicast::cluster.cta_group::%7.L2::cache_hint "
              "[%0], [%1, {%2, %3}], [%4], %5, %6;"
              :: "r"(dst), "l"(tmap_ptr), "r"(x), "r"(y), "r"(mbar_addr), "h"(cta_mask), "l"(cache_policy), "n"(CTA_GROUP)
              : "memory");
}

template <int CTA_GROUP = 1>
__device__ inline
void tma_3d_gmem2smem(int dst, const void *tmap_ptr, int x, int y, int z, int mbar_addr, uint64_t cache_policy) {
  asm volatile("cp.async.bulk.tensor.3d.shared::cta.global.mbarrier::complete_tx::bytes.cta_group::%7.L2::cache_hint "
              "[%0], [%1, {%2, %3, %4}], [%5], %6;"
              :: "r"(dst), "l"(tmap_ptr), "r"(x), "r"(y), "r"(z), "r"(mbar_addr), "l"(cache_policy), "n"(CTA_GROUP)
              : "memory");
}

template <int CTA_GROUP = 1>
__device__ inline
void tma_3d_gmem2smem_mcast(int dst, const void *tmap_ptr, int x, int y, int z, int mbar_addr, int16_t cta_mask, uint64_t cache_policy) {
  asm volatile("cp.async.bulk.tensor.3d.shared::cluster.global.mbarrier::complete_tx::bytes.multicast::cluster.cta_group::%8.L2::cache_hint "
              "[%0], [%1, {%2, %3, %4}], [%5], %6, %7;"
              :: "r"(dst), "l"(tmap_ptr), "r"(x), "r"(y), "r"(z), "r"(mbar_addr), "h"(cta_mask), "l"(cache_policy), "n"(CTA_GROUP)
              : "memory");
}

template <int CTA_GROUP = 1>
__device__ inline
void tcgen05_cp_nvfp4(int taddr, uint64_t s_desc) {
  // .32x128b corresponds to (32, 16) 8-bit scale -> 1 MMA for nvfp4.
  // .warpx4 duplicates data across 32-lane groups.
  asm volatile("tcgen05.cp.cta_group::%2.32x128b.warpx4 [%0], %1;" :: "r"(taddr), "l"(s_desc), "n"(CTA_GROUP));
}

template <int CTA_GROUP = 1>
__device__ inline
void tcgen05_alloc(int addr, int size) {
  asm volatile("tcgen05.alloc.cta_group::%2.sync.aligned.shared::cta.b32 [%0], %1;" :: "r"(addr), "r"(size), "n"(CTA_GROUP));
}

template <int CTA_GROUP = 1>
__device__ inline
void tcgen05_dealloc(int addr, int size) {
  asm volatile("tcgen05.dealloc.cta_group::%2.sync.aligned.b32 %0, %1;" :: "r"(0), "r"(size), "n"(CTA_GROUP));
}

template <int CTA_GROUP = 1>
__device__ inline
void tcgen05_commit(int mbar_addr) {
  asm volatile("tcgen05.commit.cta_group::%1.mbarrier::arrive::one.shared::cluster.b64 [%0];"
              :: "r"(mbar_addr), "n"(CTA_GROUP) : "memory");
}

template <int CTA_GROUP = 1>
__device__ inline
void tcgen05_commit_mcast(int mbar_addr, uint16_t cta_mask) {
  asm volatile("tcgen05.commit.cta_group::%2.mbarrier::arrive::one.shared::cluster.multicast::cluster.b64 [%0], %1;"
              :: "r"(mbar_addr), "h"(cta_mask), "n"(CTA_GROUP) : "memory");
}

struct COLLECTOR_USAGE {
  static constexpr char NONE[]      = "";
  static constexpr char A_FILL[]    = ".collector::a::fill";
  static constexpr char A_USE[]     = ".collector::a::use";
  static constexpr char A_LASTUSE[] = ".collector::a::lastuse";
  static constexpr char A_DISCARD[] = ".collector::a::discard";
};

template <int CTA_GROUP = 1, const char *collector_usage = COLLECTOR_USAGE::NONE>
__device__ inline
void tcgen05_mma_nvfp4(
  int d_tmem,
  uint64_t a_desc,
  uint64_t b_desc,
  uint32_t i_desc,
  int scale_A_tmem,
  int scale_B_tmem,
  int enable_input_d
) {
  asm volatile(
    "{\n\t"
    ".reg .pred p;\n\t"  // predicate register enable-input-d
    "setp.ne.b32 p, %6, 0;\n\t"
    "tcgen05.mma.cta_group::%7.kind::mxf4nvf4.block_scale.block16%8 [%0], %1, %2, %3, [%4], [%5], p;\n\t"
    "}"
    :: "r"(d_tmem), "l"(a_desc), "l"(b_desc), "r"(i_desc),
       "r"(scale_A_tmem), "r"(scale_B_tmem), "r"(enable_input_d),
       "n"(CTA_GROUP), "C"(collector_usage)
  );
}

// see https://docs.nvidia.com/cuda/inline-ptx-assembly/index.html
struct SHAPE {
  static constexpr char _32x32b[]  = ".32x32b";   // 32x1 tile for each warp
  static constexpr char _16x128b[] = ".16x128b";  // 16x4 tile
  static constexpr char _16x256b[] = ".16x256b";  // 16x8 tile
};

template <int NUM_REGS, const char *SHAPE, int NUM>
__device__ inline
void tcgen05_ld(float *tmp, int row, int col) {
  int addr = (row << 16) | col;

  if constexpr (NUM_REGS == 1)
  asm volatile("tcgen05.ld.sync.aligned%3.x%4.b32 {%0}, [%1];"
              : "=f"(tmp[0]) : "r"(addr), "C"(SHAPE), "n"(NUM));
  if constexpr (NUM_REGS == 2)
  asm volatile("tcgen05.ld.sync.aligned%3.x%4.b32 {%0, %1}, [%2];"
              : "=f"(tmp[0]), "=f"(tmp[1]) : "r"(addr), "C"(SHAPE), "n"(NUM));
  if constexpr (NUM_REGS == 4)
  asm volatile("tcgen05.ld.sync.aligned%5.x%6.b32 "
              "{%0, %1, %2, %3}, [%4];"
              : "=f"(tmp[0]), "=f"(tmp[1]), "=f"(tmp[2]), "=f"(tmp[3])
              : "r"(addr), "C"(SHAPE), "n"(NUM));
  if constexpr (NUM_REGS == 8)
  asm volatile("tcgen05.ld.sync.aligned%9.x%10.b32 "
              "{ %0,  %1,  %2,  %3,  %4,  %5,  %6,  %7}, [%8];"
              : "=f"(tmp[0]), "=f"(tmp[1]), "=f"(tmp[2]), "=f"(tmp[3]), "=f"(tmp[4]), "=f"(tmp[5]), "=f"(tmp[6]), "=f"(tmp[7])
              : "r"(addr), "C"(SHAPE), "n"(NUM));
  if constexpr (NUM_REGS == 16)
  asm volatile("tcgen05.ld.sync.aligned%17.x%18.b32 "
              "{ %0,  %1,  %2,  %3,  %4,  %5,  %6,  %7, "
              "  %8,  %9, %10, %11, %12, %13, %14, %15}, [%16];"
              : "=f"(tmp[ 0]), "=f"(tmp[ 1]), "=f"(tmp[ 2]), "=f"(tmp[ 3]), "=f"(tmp[ 4]), "=f"(tmp[ 5]), "=f"(tmp[ 6]), "=f"(tmp[ 7]),
                "=f"(tmp[ 8]), "=f"(tmp[ 9]), "=f"(tmp[10]), "=f"(tmp[11]), "=f"(tmp[12]), "=f"(tmp[13]), "=f"(tmp[14]), "=f"(tmp[15])
              : "r"(addr), "C"(SHAPE), "n"(NUM));
  if constexpr (NUM_REGS == 32)
  asm volatile("tcgen05.ld.sync.aligned%33.x%34.b32 "
              "{ %0,  %1,  %2,  %3,  %4,  %5,  %6,  %7, "
              "  %8,  %9, %10, %11, %12, %13, %14, %15, "
              " %16, %17, %18, %19, %20, %21, %22, %23, "
              " %24, %25, %26, %27, %28, %29, %30, %31}, [%32];"
              : "=f"(tmp[ 0]), "=f"(tmp[ 1]), "=f"(tmp[ 2]), "=f"(tmp[ 3]), "=f"(tmp[ 4]), "=f"(tmp[ 5]), "=f"(tmp[ 6]), "=f"(tmp[ 7]),
                "=f"(tmp[ 8]), "=f"(tmp[ 9]), "=f"(tmp[10]), "=f"(tmp[11]), "=f"(tmp[12]), "=f"(tmp[13]), "=f"(tmp[14]), "=f"(tmp[15]),
                "=f"(tmp[16]), "=f"(tmp[17]), "=f"(tmp[18]), "=f"(tmp[19]), "=f"(tmp[20]), "=f"(tmp[21]), "=f"(tmp[22]), "=f"(tmp[23]),
                "=f"(tmp[24]), "=f"(tmp[25]), "=f"(tmp[26]), "=f"(tmp[27]), "=f"(tmp[28]), "=f"(tmp[29]), "=f"(tmp[30]), "=f"(tmp[31])
              : "r"(addr), "C"(SHAPE), "n"(NUM));
  if constexpr (NUM_REGS == 64)
  asm volatile("tcgen05.ld.sync.aligned%65.x%66.b32 "
              "{ %0,  %1,  %2,  %3,  %4,  %5,  %6,  %7, "
              "  %8,  %9, %10, %11, %12, %13, %14, %15, "
              " %16, %17, %18, %19, %20, %21, %22, %23, "
              " %24, %25, %26, %27, %28, %29, %30, %31, "
              " %32, %33, %34, %35, %36, %37, %38, %39, "
              " %40, %41, %42, %43, %44, %45, %46, %47, "
              " %48, %49, %50, %51, %52, %53, %54, %55, "
              " %56, %57, %58, %59, %60, %61, %62, %63}, [%64];"
              : "=f"(tmp[ 0]), "=f"(tmp[ 1]), "=f"(tmp[ 2]), "=f"(tmp[ 3]), "=f"(tmp[ 4]), "=f"(tmp[ 5]), "=f"(tmp[ 6]), "=f"(tmp[ 7]),
                "=f"(tmp[ 8]), "=f"(tmp[ 9]), "=f"(tmp[10]), "=f"(tmp[11]), "=f"(tmp[12]), "=f"(tmp[13]), "=f"(tmp[14]), "=f"(tmp[15]),
                "=f"(tmp[16]), "=f"(tmp[17]), "=f"(tmp[18]), "=f"(tmp[19]), "=f"(tmp[20]), "=f"(tmp[21]), "=f"(tmp[22]), "=f"(tmp[23]),
                "=f"(tmp[24]), "=f"(tmp[25]), "=f"(tmp[26]), "=f"(tmp[27]), "=f"(tmp[28]), "=f"(tmp[29]), "=f"(tmp[30]), "=f"(tmp[31]),
                "=f"(tmp[32]), "=f"(tmp[33]), "=f"(tmp[34]), "=f"(tmp[35]), "=f"(tmp[36]), "=f"(tmp[37]), "=f"(tmp[38]), "=f"(tmp[39]),
                "=f"(tmp[40]), "=f"(tmp[41]), "=f"(tmp[42]), "=f"(tmp[43]), "=f"(tmp[44]), "=f"(tmp[45]), "=f"(tmp[46]), "=f"(tmp[47]),
                "=f"(tmp[48]), "=f"(tmp[49]), "=f"(tmp[50]), "=f"(tmp[51]), "=f"(tmp[52]), "=f"(tmp[53]), "=f"(tmp[54]), "=f"(tmp[55]),
                "=f"(tmp[56]), "=f"(tmp[57]), "=f"(tmp[58]), "=f"(tmp[59]), "=f"(tmp[60]), "=f"(tmp[61]), "=f"(tmp[62]), "=f"(tmp[63])
              : "r"(addr), "C"(SHAPE), "n"(NUM));
  if constexpr (NUM_REGS == 128)
  asm volatile("tcgen05.ld.sync.aligned%129.x%130.b32 "
              "{ %0,  %1,  %2,  %3,  %4,  %5,  %6,  %7, "
              "  %8,  %9, %10, %11, %12, %13, %14, %15, "
              " %16, %17, %18, %19, %20, %21, %22, %23, "
              " %24, %25, %26, %27, %28, %29, %30, %31, "
              " %32, %33, %34, %35, %36, %37, %38, %39, "
              " %40, %41, %42, %43, %44, %45, %46, %47, "
              " %48, %49, %50, %51, %52, %53, %54, %55, "
              " %56, %57, %58, %59, %60, %61, %62, %63, "
              " %64, %65, %66, %67, %68, %69, %70, %71, "
              " %72, %73, %74, %75, %76, %77, %78, %79, "
              " %80, %81, %82, %83, %84, %85, %86, %87, "
              " %88, %89, %90, %91, %92, %93, %94, %95, "
              " %96, %97, %98, %99,%100,%101,%102,%103, "
              "%104,%105,%106,%107,%108,%109,%110,%111, "
              "%112,%113,%114,%115,%116,%117,%118,%119, "
              "%120,%121,%122,%123,%124,%125,%126,%127}, [%128];"
              : "=f"(tmp[ 0]), "=f"(tmp[ 1]), "=f"(tmp[ 2]), "=f"(tmp[ 3]), "=f"(tmp[ 4]), "=f"(tmp[ 5]), "=f"(tmp[ 6]), "=f"(tmp[ 7]),
                "=f"(tmp[ 8]), "=f"(tmp[ 9]), "=f"(tmp[10]), "=f"(tmp[11]), "=f"(tmp[12]), "=f"(tmp[13]), "=f"(tmp[14]), "=f"(tmp[15]),
                "=f"(tmp[16]), "=f"(tmp[17]), "=f"(tmp[18]), "=f"(tmp[19]), "=f"(tmp[20]), "=f"(tmp[21]), "=f"(tmp[22]), "=f"(tmp[23]),
                "=f"(tmp[24]), "=f"(tmp[25]), "=f"(tmp[26]), "=f"(tmp[27]), "=f"(tmp[28]), "=f"(tmp[29]), "=f"(tmp[30]), "=f"(tmp[31]),
                "=f"(tmp[32]), "=f"(tmp[33]), "=f"(tmp[34]), "=f"(tmp[35]), "=f"(tmp[36]), "=f"(tmp[37]), "=f"(tmp[38]), "=f"(tmp[39]),
                "=f"(tmp[40]), "=f"(tmp[41]), "=f"(tmp[42]), "=f"(tmp[43]), "=f"(tmp[44]), "=f"(tmp[45]), "=f"(tmp[46]), "=f"(tmp[47]),
                "=f"(tmp[48]), "=f"(tmp[49]), "=f"(tmp[50]), "=f"(tmp[51]), "=f"(tmp[52]), "=f"(tmp[53]), "=f"(tmp[54]), "=f"(tmp[55]),
                "=f"(tmp[56]), "=f"(tmp[57]), "=f"(tmp[58]), "=f"(tmp[59]), "=f"(tmp[60]), "=f"(tmp[61]), "=f"(tmp[62]), "=f"(tmp[63]),
                "=f"(tmp[64]), "=f"(tmp[65]), "=f"(tmp[66]), "=f"(tmp[67]), "=f"(tmp[68]), "=f"(tmp[69]), "=f"(tmp[70]), "=f"(tmp[71]),
                "=f"(tmp[72]), "=f"(tmp[73]), "=f"(tmp[74]), "=f"(tmp[75]), "=f"(tmp[76]), "=f"(tmp[77]), "=f"(tmp[78]), "=f"(tmp[79]),
                "=f"(tmp[80]), "=f"(tmp[81]), "=f"(tmp[82]), "=f"(tmp[83]), "=f"(tmp[84]), "=f"(tmp[85]), "=f"(tmp[86]), "=f"(tmp[87]),
                "=f"(tmp[88]), "=f"(tmp[89]), "=f"(tmp[90]), "=f"(tmp[91]), "=f"(tmp[92]), "=f"(tmp[93]), "=f"(tmp[94]), "=f"(tmp[95]),
                "=f"(tmp[96]), "=f"(tmp[97]), "=f"(tmp[98]), "=f"(tmp[99]), "=f"(tmp[100]),"=f"(tmp[101]),"=f"(tmp[102]),"=f"(tmp[103]),
                "=f"(tmp[104]),"=f"(tmp[105]),"=f"(tmp[106]),"=f"(tmp[107]),"=f"(tmp[108]),"=f"(tmp[109]),"=f"(tmp[110]),"=f"(tmp[111]),
                "=f"(tmp[112]),"=f"(tmp[113]),"=f"(tmp[114]),"=f"(tmp[115]),"=f"(tmp[116]),"=f"(tmp[117]),"=f"(tmp[118]),"=f"(tmp[119]),
                "=f"(tmp[120]),"=f"(tmp[121]),"=f"(tmp[122]),"=f"(tmp[123]),"=f"(tmp[124]),"=f"(tmp[125]),"=f"(tmp[126]),"=f"(tmp[127])
              : "r"(addr), "C"(SHAPE), "n"(NUM));
}

template <int num>
__device__ inline void
tcgen05_ld_32x32b(float *tmp, int row, int col) {
  // each 32x32b tile uses 1 register per thread
  tcgen05_ld<num, SHAPE::_32x32b, num>(tmp, row, col);
}

template <int num>
__device__ inline
void tcgen05_ld_16x128b(float *tmp, int row, int col) {
  // each 16x128b tile uses 2 registers per thread
  tcgen05_ld<num * 2, SHAPE::_16x128b, num>(tmp, row, col);
}

template <int num>
__device__ inline
void tcgen05_ld_16x256b(float *tmp, int row, int col) {
  // each 16x256b tile uses 4 registers per thread
  tcgen05_ld<num * 4, SHAPE::_16x256b, num>(tmp, row, col);
}

constexpr int BLOCK_M = 128;
constexpr int BLOCK_N = 128;
constexpr int BLOCK_K = 256;
constexpr int NUM_WARPS = 6;
constexpr int TB_SIZE = NUM_WARPS * WARP_SIZE;

template <typename T>
__device__ __inline__
T warp_uniform(T x) { return __shfl_sync(0xFFFF'FFFF, x, 0); }

template <int NUM_GROUPS>
struct Arguments {
  CUtensorMap A_tmap_list[NUM_GROUPS];
  CUtensorMap B_tmap_list[NUM_GROUPS];
  char *SFA_ptr_list[NUM_GROUPS];
  char *SFB_ptr_list[NUM_GROUPS];
  half *C_ptr_list[NUM_GROUPS];
  int M_list[NUM_GROUPS];
  int N_list[NUM_GROUPS];
  int K_list[NUM_GROUPS];
};

template <int NUM_GROUPS, int NUM_STAGES>
__global__
__launch_bounds__(TB_SIZE)
void kernel_cutlass (const __grid_constant__ Arguments<NUM_GROUPS> args) {
  const int tid = threadIdx.x;
  const int bid_n = warp_uniform(blockIdx.x);
  const int bid_m = warp_uniform(blockIdx.y);
  const int group_id = warp_uniform(blockIdx.z);

  // is there a perf hit for indexing into argument list with a dynamic value?
  const int M = args.M_list[group_id];
  const int N = args.N_list[group_id];
  const int K = args.K_list[group_id];

  // early exit
  if (bid_m >= cdiv(M, BLOCK_M) || bid_n >= N / BLOCK_N)
    return;

  const int lane_id = tid % WARP_SIZE;
  const int warp_id = warp_uniform(tid / WARP_SIZE);

  const int off_m = bid_m * BLOCK_M;
  const int off_n = bid_n * BLOCK_N;

  // set up smem
  extern __shared__ __align__(1024) char smem_ptr[];
  const int smem = static_cast<int>(__cvta_generic_to_shared(smem_ptr));
  constexpr int A_size = BLOCK_M * BLOCK_K / 2;
  constexpr int B_size = BLOCK_N * BLOCK_K / 2;
  constexpr int SF_size = 128 * BLOCK_K / 16;
  constexpr int STAGE_SIZE = A_size + B_size + SF_size * 2;

  // set up mbarriers and tmem
  const int tma_mbar_addr = smem + NUM_STAGES * STAGE_SIZE;
  const int mma_mbar_addr = tma_mbar_addr + NUM_STAGES * 8;
  const int mainloop_mbar_addr = mma_mbar_addr + NUM_STAGES * 8;

  constexpr uint64_t cache_A = EVICT_NORMAL;
  constexpr uint64_t cache_B = EVICT_FIRST;

  constexpr int bar_epilogue = 2;
  const int rest_k = K / 16 / 4;

  auto A_tmap = args.A_tmap_list + group_id;
  auto B_tmap = args.B_tmap_list + group_id;

  if (warp_id == 0 && elect_sync()) {
    asm volatile("prefetch.tensormap [%0];" :: "l"(A_tmap) : "memory");
    asm volatile("prefetch.tensormap [%0];" :: "l"(B_tmap) : "memory");
  }
  else if (warp_id == 1 && elect_sync()) {
    // 1 thread init mbarrier
    for (int i = 0; i < NUM_STAGES; i++) {
      mbarrier_init(tma_mbar_addr + i * 8, 1);
      mbarrier_init(mma_mbar_addr + i * 8, 1);
    }
    mbarrier_init(mainloop_mbar_addr, 1);
    asm volatile("fence.mbarrier_init.release.cluster;");  // visible to async proxy
  }

  __syncthreads();

  const int num_iters = K / BLOCK_K;

  if (warp_id == NUM_WARPS - 2) {
    // TMA warp
    if (elect_sync()) {
      int stage_id = 0;
      int mma_phase = 1;

      const char *SFA_ptr = args.SFA_ptr_list[group_id];
      const char *SFB_ptr = args.SFB_ptr_list[group_id];

      for (int iter_k = 0; iter_k < num_iters; iter_k++) {
        // select tma mbar and smem
        const int mbar_addr = tma_mbar_addr + stage_id * 8;
        const int A_smem = smem + stage_id * STAGE_SIZE;
        const int B_smem = A_smem + A_size;
        const int SFA_smem = B_smem + B_size;
        const int SFB_smem = SFA_smem + SF_size;

        // wait MMA
        mbarrier_wait(mma_mbar_addr + stage_id * 8, mma_phase);

        // issue MMA
        tma_3d_gmem2smem(A_smem, A_tmap, 0, off_m, iter_k, mbar_addr, cache_A);
        tma_3d_gmem2smem(B_smem, B_tmap, 0, off_n, iter_k, mbar_addr, cache_B);

        const char *SFA_src = SFA_ptr + bid_m * rest_k * 512 + iter_k * 2048;
        const char *SFB_src = SFB_ptr + bid_n * rest_k * 512 + iter_k * 2048;
        tma_gmem2smem(SFA_smem, SFA_src, SF_size, mbar_addr, cache_A);
        tma_gmem2smem(SFB_smem, SFB_src, SF_size, mbar_addr, cache_B);

        mbarrier_arrive_expect_tx(mbar_addr, STAGE_SIZE);  // signal TMA done
        stage_id = (stage_id + 1) % NUM_STAGES;
        if (stage_id == 0)
          mma_phase ^= 1;
      }
    }
  }
  else if (warp_id == NUM_WARPS - 1) {
    // MMA warp
    tcgen05_alloc(mainloop_mbar_addr + 8, 512);  // allocate tmem

    // instruction desc
    constexpr uint32_t MMA_M = BLOCK_M;
    constexpr uint32_t MMA_N = BLOCK_N;
    constexpr uint32_t i_desc = (1U << 7U)   // atype=E2M1
                              | (1U << 10U)  // btype=E2M1
                              | (MMA_N >> 3U << 17U)
                              | (MMA_M >> 7U << 27U)
                              ;

    if (elect_sync()) {
      int stage_id = 0;
      int tma_phase = 0;

      for (int iter_k = 0; iter_k < num_iters; iter_k ++) {
        // select smem
        const int A_smem = smem + stage_id * STAGE_SIZE;
        const int B_smem = A_smem + A_size;
        const int SFA_smem = B_smem + B_size;
        const int SFB_smem = SFA_smem + SF_size;

        // set up smem desc
        // AB: 128-byte swizzling
        constexpr uint64_t AB_desc = (desc_encode(8 * 128) << 32ULL) | (1ULL << 46ULL) | (2ULL << 61ULL);
        uint64_t a_desc = AB_desc | (A_smem >> 4);
        uint64_t b_desc = AB_desc | (B_smem >> 4);

        // SF: no swizzling
        constexpr uint64_t SF_desc = (desc_encode(8 * 16) << 32ULL) | (1ULL << 46ULL);
        uint64_t sfa_desc = SF_desc | (SFA_smem >> 4);
        uint64_t sfb_desc = SF_desc | (SFB_smem >> 4);

        // first BLOCK_N columns are for acc
        // each SF consumes 16 columns per BLOCK_K=256
        int sfa_tmem = BLOCK_N;
        int sfb_tmem = BLOCK_N + 16;

        // wait TMA
        mbarrier_wait(tma_mbar_addr + stage_id * 8, tma_phase);

        // manual unroll 1st iteration
        tcgen05_cp_nvfp4(sfa_tmem, sfa_desc);
        tcgen05_cp_nvfp4(sfb_tmem, sfb_desc);
        tcgen05_mma_nvfp4(0, a_desc, b_desc, i_desc, sfa_tmem, sfb_tmem, iter_k);

        for (int k = 1; k < BLOCK_K / MMA_K; k++) {
          // next 4 columns
          sfa_tmem += 4;
          sfb_tmem += 4;

          // next 512-byte
          sfa_desc += (512 >> 4);
          sfb_desc += (512 >> 4);

          // next 32-byte
          a_desc += (32 >> 4);
          b_desc += (32 >> 4);

          tcgen05_cp_nvfp4(sfa_tmem, sfa_desc);
          tcgen05_cp_nvfp4(sfb_tmem, sfb_desc);
          tcgen05_mma_nvfp4(0, a_desc, b_desc, i_desc, sfa_tmem, sfb_tmem, 1);
        }

        tcgen05_commit(mma_mbar_addr + stage_id * 8);  // signal MMA done
        stage_id = (stage_id + 1) % NUM_STAGES;
        if (stage_id == 0)
          tma_phase ^= 1;
      }

      tcgen05_commit(mainloop_mbar_addr);  // signal mainloop done
    }
  }
  else {
    // epilogue warps
    if (warp_id == 0)
      mbarrier_wait(mainloop_mbar_addr, 0);
    asm volatile("bar.sync %0, %1;" :: "n"(bar_epilogue), "r"(4 * WARP_SIZE) : "memory");
    asm volatile("tcgen05.fence::after_thread_sync;");

    half *C_ptr = args.C_ptr_list[group_id];

    constexpr int WIDTH = std::min(BLOCK_N, 8);

    for (int m = 0; m < 32 / 16; m++)
      for (int n = 0; n < BLOCK_N / WIDTH; n++) {
        float out[WIDTH / 2];

        constexpr int num = WIDTH / 8;  // each 16x256b tile is 8-element wide
        const int tmem_row = warp_id * 32 + m * 16;
        tcgen05_ld_16x256b<num>(out, tmem_row, n * WIDTH);
        asm volatile("tcgen05.wait::ld.sync.aligned;");

        for (int i = 0; i < WIDTH / 8; i++) {
          const int row = off_m + tmem_row + lane_id / 4;
          const int col = off_n + n * WIDTH + i * 8 + (lane_id % 4) * 2;

          if (row < M)
            reinterpret_cast<half2 *>(C_ptr + (row + 0) * N + col)[0] = __float22half2_rn({out[i * 4 + 0], out[i * 4 + 1]});
          if (row + 8 < M)
            reinterpret_cast<half2 *>(C_ptr + (row + 8) * N + col)[0] = __float22half2_rn({out[i * 4 + 2], out[i * 4 + 3]});
        }
      }

    asm volatile("bar.sync %0, %1;" :: "n"(bar_epilogue), "r"(4 * WARP_SIZE) : "memory");

    if (warp_id == 0)
      tcgen05_dealloc(0, 512);
  }
}

void check_cu(CUresult err) {
  if (err == CUDA_SUCCESS) return;
  const char *error_msg_ptr;
  if (cuGetErrorString(err, &error_msg_ptr) != CUDA_SUCCESS)
    error_msg_ptr = "unable to get error string";
  TORCH_CHECK(false, "cuTensorMapEncodeTiled error: ", error_msg_ptr);
}

void check_cuda(cudaError_t err) {
  if (err == cudaSuccess) return;
  TORCH_CHECK(false, cudaGetErrorString(err));
}

void init_AB_tmap(
  CUtensorMap *tmap,
  void *ptr,
  uint64_t global_height, uint64_t global_width,
  uint32_t shared_height, uint32_t shared_width
) {
  constexpr uint32_t rank = 3;
  uint64_t globalDim[rank]       = {256, global_height, global_width / 256};
  uint64_t globalStrides[rank-1] = {global_width / 2, 128};  // in bytes
  uint32_t boxDim[rank]          = {256, shared_height, shared_width / 256};
  uint32_t elementStrides[rank]  = {1, 1, 1};

  auto err = cuTensorMapEncodeTiled(
    tmap,
    CUtensorMapDataType::CU_TENSOR_MAP_DATA_TYPE_16U4_ALIGN8B,
    rank,
    ptr,
    globalDim,
    globalStrides,
    boxDim,
    elementStrides,
    CUtensorMapInterleave::CU_TENSOR_MAP_INTERLEAVE_NONE,
    CUtensorMapSwizzle::CU_TENSOR_MAP_SWIZZLE_128B,
    CUtensorMapL2promotion::CU_TENSOR_MAP_L2_PROMOTION_NONE,
    CUtensorMapFloatOOBfill::CU_TENSOR_MAP_FLOAT_OOB_FILL_NONE
  );
  //check_cu(err);
}

template <int NUM_GROUPS>
void group_gemm_launch(
  at::TensorList A_list,
  at::TensorList B_list,
  at::TensorList SFA_list,
  at::TensorList SFB_list,
  at::TensorList C_list
) {
  Arguments<NUM_GROUPS> args;
  int max_M = 0;
  int max_N = 0;

  for (int i = 0; i < NUM_GROUPS; i++) {
    const int M = A_list[i].size(0);
    const int N = B_list[i].size(0);
    const int K = A_list[i].size(1) * 2;

    init_AB_tmap(args.A_tmap_list + i, A_list[i].data_ptr(), M, K, BLOCK_M, BLOCK_K);
    init_AB_tmap(args.B_tmap_list + i, B_list[i].data_ptr(), N, K, BLOCK_N, BLOCK_K);
    args.SFA_ptr_list[i] = reinterpret_cast<char *>(SFA_list[i].data_ptr());
    args.SFB_ptr_list[i] = reinterpret_cast<char *>(SFB_list[i].data_ptr());
    args.C_ptr_list[i] = reinterpret_cast<half *>(C_list[i].data_ptr());
    args.M_list[i] = M;
    args.N_list[i] = N;
    args.K_list[i] = K;

    max_M = std::max(max_M, M);
    max_N = std::max(max_N, N);
  }

  // launch max grid_m x grid_n to cover the largest group
  dim3 grid(max_N / BLOCK_N, cdiv(max_M, BLOCK_M), NUM_GROUPS);

  constexpr int AB_size = (BLOCK_M + BLOCK_N) * (BLOCK_K / 2);
  constexpr int SF_size = 128 * (BLOCK_K / 16) * 2;

  constexpr int sm100_size = 227'000;
  constexpr int dynamic_size = AB_size + SF_size + 2 * 8;  // 1 tma_mbar, 1 mma_mbar
  constexpr int static_size = 8 + 4;  // 1 mainloop_mbar, tmem_addr
  constexpr int NUM_STAGES = (sm100_size - static_size) / dynamic_size;

  constexpr int smem_size = dynamic_size * NUM_STAGES + static_size;

  // cutlass incantation (this affects ptxas)
  auto this_kernel = kernel_cutlass<NUM_GROUPS, NUM_STAGES>;
  cudaFuncSetAttribute(this_kernel, cudaFuncAttributeMaxDynamicSharedMemorySize, smem_size);
  this_kernel<<<grid, TB_SIZE, smem_size>>>(args);
}

void group_gemm(
  at::TensorList A_list,
  at::TensorList B_list,
  at::TensorList SFA_list,
  at::TensorList SFB_list,
  at::TensorList C_list
) {

#define LAUNCH(NUM_GROUPS) { \
  group_gemm_launch<NUM_GROUPS>(A_list, B_list, SFA_list, SFB_list, C_list); \
  }

  if (false) {}
  else if (A_list.size() == 2) LAUNCH(2)
  else if (A_list.size() == 3) LAUNCH(3)
  else if (A_list.size() == 4) LAUNCH(4)
  else if (A_list.size() == 8) LAUNCH(8)

#undef LAUNCH
}

TORCH_LIBRARY(my_module, m) {
  m.def("group_gemm(Tensor[] A_list, Tensor[] B_list, Tensor[] SFA_list, Tensor[] SFB_list, Tensor(a!)[] C_list) -> ()");
  m.impl("group_gemm", &group_gemm);
}
"""

load_inline(
    "group_gemm",
    cpp_sources="",
    cuda_sources=CUDA_SRC,
    verbose=True,
    is_python_module=False,
    no_implicit_headers=True,
    extra_cuda_cflags=[
        "-O3",
        "-gencode=arch=compute_100a,code=sm_100a",
        "--use_fast_math",
        "--expt-relaxed-constexpr",
        "--relocatable-device-code=false",
        "-lineinfo",
        "-Xptxas=-v",
        # "--keep",
        # "--keep-dir",
        # f"{Path(__file__).parent}/tmp",
    ],
    extra_ldflags=["-lcuda"],
)
group_gemm = torch.ops.my_module.group_gemm


def ref(A_list, B_list, SFA_list, SFB_list, C_list):
    for a, b, sfa, sfb, c in zip(A_list, B_list, SFA_list, SFB_list, C_list):
        torch._scaled_mm(
            a[..., 0],
            b[..., 0].T,
            sfa.permute(5, 2, 4, 0, 1, 3).view(-1),
            sfb.permute(5, 2, 4, 0, 1, 3).view(-1),
            out=c[..., 0],
        )


def custom_kernel(data: input_t) -> output_t:
    abc_list, _, sf_list, shape_list = data

    A_list, B_list, C_list = zip(*abc_list)
    SFA_list, SFB_list = zip(*sf_list)

    # don't handle K % 256 != 0
    if any(k % 256 for _, _, k, _ in shape_list):
        ref(A_list, B_list, SFA_list, SFB_list, C_list)
    else:
        group_gemm(A_list, B_list, SFA_list, SFB_list, C_list)

    # torch.cuda.synchronize()
    return C_list