Atlas Cloud

Posted on May 26

Writing High-Performance Kernels in TileLang, from GEMM to MLA

#deeplearning #llm #performance #python

If you write GPU kernels, you live somewhere on a spectrum. At one end is Triton: quick to write, but the compiler makes most of the layout and shared-memory decisions for you. At the other end is CUTLASS / CuTe: total control, at the cost of a lot of template machinery. TileLang sits in the middle. You write Python, but you say explicitly what lives in shared memory, how the pipeline is staged, and how warps split the work — and a layout inference pass fills in the rest.

In this post we'll cover the mental model, write a GEMM, and then build up to a real production kernel: DeepSeek's MLA decode, where the interesting decisions actually show up. The goal is not to be exhaustive. It's to show what you think about tiles, and where TileLang quietly does the hard parts for you. We'll finish with a more typical story from production — a kernel where the win wasn't speed at all.

The mental model

Here's the whole idea in three points.

A tile is a first-class object. A shaped chunk of data (block_M × block_K, say) is owned and operated on by a thread block, a warp, or a thread. You stop thinking purely at the thread-block level the way you do in Triton, and you stop hand-managing individual threads the way you do in CUDA.
You place buffers in the memory hierarchy yourself. You declare what goes to shared memory (T.alloc_shared), what goes to registers (T.alloc_fragment), and what's thread-local. This is the biggest difference from Triton, which hides shared-memory allocation and staging inside the compiler.
The compiler infers the thread mapping. Once you've said where a tile lives and what operation runs on it (a copy, a gemm, a reduce), a layout inference pass parallelizes it across threads and works out the register and shared-memory layouts. You can override it when you need to, but most of the time you don't. This pass is the load-bearing feature — by the time we get to MLA you'll see why.

If you're coming from Triton, here's the rough mapping.

	Triton	TileLang
Granularity	thread block + implicit vectorization	tile (block / warp / thread)
Shared memory	managed by the compiler	explicit alloc_shared + copy
Layout	the compiler decides	inferred, but you can annotate
Pipelining	tl.range + compiler	explicit T.Pipelined(num_stages=)
Tensor Core	tl.dot	T.gemm with a selectable warp policy
Backends	NVIDIA (mainly) / AMD	NVIDIA / AMD / CPU / WebGPU / CuTeDSL, plus Ascend & MUSA forks

The short version: if you want fine control over blocking, pipeline depth, and warp partitioning without writing CUTLASS, TileLang is the sweet spot. For simple elementwise or light fusion, Triton is still quicker to reach for.

Getting set up

conda create -n tilelang python=3.10 -y
conda activate tilelang
pip install tilelang                 # prebuilt wheel, easiest path

If you're going to touch the compiler passes, build from source instead (you'll need a local LLVM/CUDA toolchain):

git clone --recursive https://github.com/tile-ai/tilelang.git
cd tilelang && pip install -r requirements-dev.txt
pip install -e . -v --no-build-isolation

Let's write a GEMM

We'll start with the kernel everyone starts with: C = ReLU(A @ B). It's small, but it touches every primitive that matters — explicit buffers, parallel copy, software pipelining, the Tensor Core call, and an L2 swizzle.

import tilelang
import tilelang.language as T
import torch

@tilelang.jit
def matmul(M, N, K, block_M, block_N, block_K,
           dtype="float16", accum_dtype="float"):

    @T.prim_func
    def matmul_relu_kernel(
        A: T.Tensor((M, K), dtype),
        B: T.Tensor((K, N), dtype),
        C: T.Tensor((M, N), dtype),
    ):
        # grid dims: (#blocks along N, #blocks along M); 128 threads per block
        with T.Kernel(T.ceildiv(N, block_N), T.ceildiv(M, block_M),
                      threads=128) as (bx, by):

            # Say where each tile lives, explicitly.
            A_shared = T.alloc_shared((block_M, block_K), dtype)         # shared memory
            B_shared = T.alloc_shared((block_K, block_N), dtype)
            C_local  = T.alloc_fragment((block_M, block_N), accum_dtype) # register accumulator

            T.use_swizzle(panel_size=4, order="col")   # optional: better L2 reuse
            T.clear(C_local)                           # zero the accumulator

            for ko in T.Pipelined(T.ceildiv(K, block_K), num_stages=3):
                T.copy(A[by * block_M, ko * block_K], A_shared)   # global -> shared
                T.copy(B[ko * block_K, bx * block_N], B_shared)
                T.gemm(A_shared, B_shared, C_local)               # tile-level MMA

            for i, j in T.Parallel(block_M, block_N):             # fused ReLU
                C_local[i, j] = T.max(C_local[i, j], 0)

            T.copy(C_local, C[by * block_M, bx * block_N])        # write back

    return matmul_relu_kernel


M = N = K = 1024
kernel = matmul(M, N, K, block_M=128, block_N=128, block_K=64)
a = torch.randn(M, K, device="cuda", dtype=torch.float16)
b = torch.randn(K, N, device="cuda", dtype=torch.float16)
c = torch.empty(M, N, device="cuda", dtype=torch.float16)
kernel(a, b, c)

torch.testing.assert_close(c, torch.relu(a @ b), rtol=1e-2, atol=1e-2)
print("gemm ok")

Here is what each piece is doing.

Three buffers, three levels. A_shared and B_shared live in shared memory; C_local lives in registers. Accumulator in registers, operands staged through shared memory — that's the standard GEMM recipe, except here you write it down. That's the whole difference from Triton in one line.
T.copy is sugar for a parallel copy. It expands into a T.Parallel-style move, and the compiler derives a vectorized, coalesced global→shared transfer from it. When the copy sits inside T.Pipelined, it becomes cp.async automatically.
T.Pipelined(extent, num_stages=N) is a software pipeline. num_stages=3 means triple buffering — while you compute K-tile ko, the loads for ko+1 and ko+2 are already in flight. In Triton, this is a compile flag; here it's just the loop, which is easier to reason about.
T.gemm(A, B, C) is the tile-level matmul. It lowers to CuTe/MMA on NVIDIA and the matching intrinsically on AMD. It also takes transpose_A / transpose_B and a policy=T.GemmWarpPolicy.* that controls how warps split the output tile. Hold onto that policy argument — it's the whole story when we get to MLA.
T.use_swizzle reorders how thread blocks are scheduled so that blocks adjacent in L2 run close together in time. Usually a few percent of free bandwidth.

The figure below maps all of this onto the hardware. It's worth reading against the code, because the labeled spots are exactly where TileLang hands you control that Triton keeps for itself.

A few primitives you'll reach for

You can write most kernels with a small vocabulary.

Allocate: T.alloc_shared, T.alloc_fragment (registers), T.alloc_local.
Move and init: T.copy(src, dst) between any two levels; T.clear, T.fill.
Compute: T.gemm(...); T.Parallel(d0, d1, ...) for elementwise loops (this is the entry point for layout inference); T.reduce_max / T.reduce_sum; scalar math like T.exp, T.exp2, T.max, T.infinity.
Schedule: T.Pipelined(extent, num_stages=), T.use_swizzle(...), T.annotate_layout(...) when you need a specific layout (bank-conflict avoidance, custom swizzle).
Dynamic shape: M = T.dynamic("m") so you don't recompile per shape (it's called T.symbolic in some versions).

Checking your work

Two things you'll want often. To see what the compiler actually emitted:

print(kernel.get_kernel_source())     # generated CUDA / HIP

And to time it:

profiler = kernel.get_profiler(tensor_supply_type=tilelang.TensorSupplyType.Normal)
print(f"latency: {profiler.do_bench()} ms")

T.print(buf) prints a tile from inside the kernel, and the repo's examples/plot_layout draws the memory layout, which is handy when you're chasing a bank conflict or checking a swizzle.

Now a real one: MLA decode

The GEMM shows the mechanics. This next one shows why they matter. We'll walk through DeepSeek's MLA (Multi-Head Latent Attention) decode kernel, because it's the cleanest example of TileLang earning its keep. The TileLang reference lands at roughly FlashMLA's H100 performance (benchmarked at batch 64/128 in fp16, comfortably ahead of Triton and FlashInfer) in about 80 lines of Python. The interesting question is how, because the hard part of MLA isn't the math — it's register pressure.

Let's review the loop everyone knows. Every FlashAttention-family kernel has the same shape. Per query block, you stream over key/value blocks and keep a running max and denominator, so the full score matrix never lands in memory:

# acc_s : [block_M, block_N]  scores for this KV block
# acc_o : [block_M, dim]      output accumulator
for i in range(num_kv_blocks):
    acc_s = Q @ K[i].T
    m_prev       = scores_max
    scores_max   = max(m_prev, rowmax(acc_s))
    scores_scale = exp(m_prev - scores_max)
    acc_o *= scores_scale                       # rescale prior output
    acc_s  = exp(acc_s - scores_max)            # probabilities
    acc_o += acc_s @ V[i]

Both acc_s and acc_o want to stay in registers. For MHA or GQA, that's fine. For MLA, it isn't.

Here's where it gets hard. MLA's head dimensions are big: query and key are 576 wide (a 512-wide "nope" part with no positional encoding, plus a 64-wide "rope" part), and value is 512. So acc_o = [block_M, 512], and it has to stay resident in registers across the whole KV loop.

Now bring in the hardware. On Hopper, the fast path is wgmma.mma_async, which ties 4 warps (128 threads) into one warpgroup and requires a minimum M of 64. So the smallest M one warpgroup can own is 64, which means one warpgroup would be holding a 64 × 512 accumulator. That's too big for a single warpgroup's register file. It spills, and performance falls off a cliff.

The fix is to split the output across two warpgroups. You can't shrink M below 64, so the only axis left is dim. Use two warpgroups: WG0 owns acc_o[:, :256], WG1 owns acc_o[:, 256:]. Now each one holds a 64 × 256 accumulator, which fits. That creates a second problem, though: the P @ V step (with policy=FullCol, each warpgroup producing one column slab of the output) needs the complete acc_s, but in Q @ K each warpgroup only naturally computed half of it. The resolution is a shared-memory swap. During Q @ K, each warpgroup writes its half of acc_s to shared memory and reads back the other warpgroup's half, so afterward both hold the full acc_s and can each compute their slab of acc_o. The diagram above is exactly that: split the scores, swap through S_shared, split the output.

In CuTe you'd hand-write the layouts, the swizzles, the Tensor Core alignment, and the producer/consumer sync to pull this off. The reason it collapses to ~80 lines here is layout inference.

Let's break down what layout inference does. You annotate intent on the T.gemm calls, and it propagates the constraints through the program for you:

policy=FullCol on P @ V means each warpgroup needs the full acc_s, so acc_s = [block_M, block_N].
That propagates back to the staging buffer, so S_shared in T.copy(S_shared, acc_s) is also [block_M, block_N].
And forward into Q @ K: with FullCol, each warpgroup's score slab is [block_M, block_N/2].

The key insight is that you never write any of those shapes. You pick the warp policy and write the math; the shapes, the swizzled layouts, and the warp-specialized producer/consumer code all come out of inference.

The kernel skeleton. In MLA decode the query splits into a "nope" part (Q, dim 512) and a "rope" part (Q_pe, dim 64), and the compressed latent serves as both K and V. So the score is a sum of two GEMMs, and the output is one more. The inner loop looks like this (a representative skeleton, not line-exact — see example_mla_decode.py):

# acc_s = Q_nope @ KV^T + Q_rope @ K_pe^T
T.gemm(Q_shared,    KV_shared,   acc_s, transpose_B=True,
       policy=T.GemmWarpPolicy.FullCol, clear_accum=True)
T.gemm(Q_pe_shared, K_pe_shared, acc_s, transpose_B=True,
       policy=T.GemmWarpPolicy.FullCol)

# online softmax
T.copy(scores_max, scores_max_prev)
T.fill(scores_max, -T.infinity(accum_dtype))
T.reduce_max(acc_s, scores_max, dim=1, clear=False)
# ... exp, rescale acc_o by scores_scale, reduce_sum into logsum ...

# acc_o += P @ V  (V is the same latent KV)
T.copy(acc_s, acc_s_cast)
T.gemm(acc_s_cast, KV_shared, acc_o, policy=T.GemmWarpPolicy.FullCol)

The S_shared exchange between the two warpgroups is the part inference inserts for you, once the FullCol policies force acc_s to be full per warpgroup.

The nice part: the optimizations are one line each. This is where TileLang pays off — the whole performance toolkit is one-liners, and the messy lowering is handled for you.

Threadblock swizzling for L2 reuse: T.use_swizzle(panel_size, order="row").
Shared-memory swizzling for bank conflicts: T.annotate_layout({S_shared: T.layout.make_swizzled_layout(S_shared)}) — XOR-style address remapping so concurrent accesses spread across banks instead of serializing.
Warp specialization: you write a plain script, and it's lowered into a producer warpgroup (TMA loads) plus consumer warpgroups, with all the mbarrier sync generated. None of it shows up in your code.
Pipelining: T.Pipelined(range, num_stages) overlaps loads with compute — more stages, more overlap, but more shared memory, so it's a knob.
Split-KV (FlashDecoding-style): when the batch is small and the SMs are idle, split the KV context across SMs and merge. It's a num_split parameter plus a combine kernel.

So the genuinely hard reasoning — register budget against the M≥64 floor, who owns what across warpgroups, the shared-memory swap — you express by choosing a policy and writing the math. Everything that would be hundreds of fragile lines in CuTe is inference and codegen. That's the pitch, and MLA is where it's most convincing.

One of our own: a drop-in RMSNorm at AtlasCloud

The last example is one of our own production kernels at AtlasCloud, from the Wan video-generation VAE on H100/H200. It's a great illustration of the other thing TileLang is excellent at: covering a config a hand-tuned kernel can't reach, with a clean drop-in.

The setup. We already ship a hand-tuned fused RMSNorm + SiLU kernel. It's fast, and it's compiled for the hidden dims D ∈ {96, 192, 384} that one model config uses. A newer config needs channel widths like {160, 256, 320, 512, 640, 1024}, so on that config the hand-tuned fast path can't run. We wrote a TileLang drop-in to cover exactly that gap.

The TileLang kernel. A drop-in with the same interface (BTHWC in/out, same math, same eps) that supports any C that's a multiple of 32. Two passes, fully coalesced, FP32 accumulator:

@T.prim_func
def main(X:     T.Tensor((M, C), dtype),      # M = B*T*H*W rows
         gamma: T.Tensor((C,),  dtype),
         Y:     T.Tensor((M, C), dtype)):
    with T.Kernel(T.ceildiv(M, BLOCK_M), threads=128) as bm:
        X_chunk = T.alloc_shared((BLOCK_M, BLOCK_C), dtype)
        ss      = T.alloc_fragment((BLOCK_M,), accum_dtype)   # FP32 sum-of-squares
        # pass 1: loop over C in BLOCK_C chunks, accumulate sum of squares in FP32
        # rinv = rsqrt(ss / C + 1e-5)
        # pass 2: re-load X, y = silu(x * gamma * rinv), write back

BLOCK_C is 128/64/32 depending on C, to respect the TMA boxDim ≤ 256 limit, and the FP32 accumulator keeps the sum of squares from overflowing in FP16. Dispatch keeps the hand-tuned path where it works and only falls back when it has to:

_ATLAS_SUPPORTED_D = (96, 192, 384)

def rms_silu_dispatch(x, gamma, out):
    if x.shape[-1] in _ATLAS_SUPPORTED_D:
        atlas_rms_norm_silu(x, gamma, out=out)        # keep the hand-tuned path
    else:
        tilelang_rms_silu_bthwc(x, gamma, out=out)    # cover the gap

What it gained us. All upside, and it's a true drop-in — same interface, same math, same eps, so it slots in behind the existing dispatch with no call-site changes.

What	Gain
Previously-unsupported config	0 → 1 — it runs now (the headline win)
Attention-block RMSNorm vs the eager PyTorch norm it replaced	42 μs → 20 μs (~2× faster)
End-to-end VAE at production resolution (720×1280, 21 frames)	~1.79× encode, ~1.78× decode

The first row is the real point: TileLang let us serve a model config that previously had no fast path at all, without touching the hand-tuned kernel that already works for the other config. One drop-in, written in Python, and a whole model path went from "throws" to "ships."

Where TileLang shines

Fine control over blocking, pipeline stages, and warp partitioning, without writing CUTLASS/CuTe.
Structurally complex, layout-sensitive kernels: GEMM variants, the FlashAttention family, MLA, linear attention, dequant-fused quant GEMM, MoE routing.
Covering an op or config your hand-tuned kernels don't reach (a hidden dim outside the instantiated set, an unusual layout) — and beating an eager fallback while you're there.
One kernel body across backends (NVIDIA / AMD / vendor forks).
The whole optimization toolkit is one call at a time — T.use_swizzle, T.annotate_layout, T.Pipelined, warp specialization, split-KV — with the lowering handled for you.

Wrapping up

The cool part of TileLang is that the hard reasoning stays in your head, not in boilerplate. You decide how to split work across warps, where buffers live, and how deep the pipeline runs — and then layout inference and warp specialization turn that into the register layouts, the swizzles, and the producer/consumer dance that would otherwise be hundreds of lines of CuTe. You pick a policy and write the math. That's the whole pitch, and it's why an 80-line MLA kernel can sit next to a hand-tuned CUTLASS one.

DEV Community