Skip to content
Published on

Three Layers of Writing a Kernel — Comparing CUDA C++, Triton, and CUTLASS on the Same Problem

Share
Authors

Introduction — What Writing the Same Kernel Twice Reveals

Write a row-wise softmax kernel in CUDA C++ and, once you add the warp-shuffle reduction, the shared-memory partial sums, and two rounds of block synchronization, you end up with about 80 lines. Write the same kernel in Triton and it's 12 lines. Yet the performance is comparable, or Triton is slightly faster.

The first time you see this, you start to wonder, "what's the point of CUDA C++ then?" A few months later you run into a kernel that simply can't be expressed in Triton, and you open CUDA C++ again. Practical judgment lives somewhere between those two points.

This post measures the three layers of writing kernels with the same ruler. The versions checked against are Triton 3.7.1 (released June 2026), CUTLASS 4.6.1 (July 2026), and CUDA Toolkit 13.3 Update 1. All three layers are changing actively, so check API details against the documentation for your own version.

What Each Layer Handles for You

Writing a kernel ultimately comes down to making six decisions. The difference between the three layers is how many of those decisions you make yourself and how many the tool makes for you.

DecisionCUDA C++TritonCUTLASS / CuTe
Grid and block shapeYouYou (only the program grid)You (via tiling policy)
Thread layout within a blockYouCompilerDecided by the template
Shared-memory allocation and swing bufferingYouCompilerDecided by the template
Avoiding bank conflicts (swizzle)YouCompilerExpressed by the layout
Asynchronous copy and pipeliningYouCompilerDecided by the template
Tensor Core instruction selectionYou (or a library)CompilerMade explicit via atoms

Here's how to read it. CUDA C++ leaves all six decisions to us. That makes it the most flexible and the most time-consuming. Triton only decides the first two and hands the rest to the compiler. That makes it short, but it also means fewer levers to pull in patterns the compiler handles poorly. CUTLASS decides all six too, but instead of writing them by hand, you assemble them from already-validated parts. That's why it delivers the best performance in the GEMM family, at the cost of a steeper learning curve.

CUDA C++ — The Layer That Instructs the Hardware Directly

Let's start with the baseline. We write a row-wise softmax by hand. One block owns one row, and the max and the sum are each reduced across the whole block.

// softmax_cuda.cu
// Build: nvcc -O3 -arch=sm_80 softmax_cuda.cu -o softmax_cuda
#include <cstdio>
#include <cfloat>
#include <cuda_runtime.h>

static const int WARP = 32;
static const int BLOCK = 256;             // 8 warps per block

__inline__ __device__ float warpMax(float v) {
  for (int off = WARP / 2; off > 0; off >>= 1)
    v = fmaxf(v, __shfl_xor_sync(0xffffffffu, v, off));
  return v;
}
__inline__ __device__ float warpSum(float v) {
  for (int off = WARP / 2; off > 0; off >>= 1)
    v += __shfl_xor_sync(0xffffffffu, v, off);
  return v;
}

__global__ void softmaxKernel(float *out, const float *in, int nCols) {
  __shared__ float part[BLOCK / WARP];
  __shared__ float bcast;                 // Slot dedicated to broadcasting

  const long row = blockIdx.x;
  const float *src = in + row * nCols;
  float *dst = out + row * nCols;

  const int tid = threadIdx.x;
  const int lane = tid % WARP, wid = tid / WARP;
  const int nWarps = BLOCK / WARP;

  // Pass 1: max of the row
  float m = -FLT_MAX;
  for (int i = tid; i < nCols; i += BLOCK) m = fmaxf(m, src[i]);
  m = warpMax(m);
  if (lane == 0) part[wid] = m;
  __syncthreads();
  if (tid == 0) {
    float t = part[0];
    for (int w = 1; w < nWarps; w++) t = fmaxf(t, part[w]);
    bcast = t;
  }
  __syncthreads();
  m = bcast;
  __syncthreads();                        // Required before part[] is reused

  // Pass 2: sum of exp
  float s = 0.f;
  for (int i = tid; i < nCols; i += BLOCK) s += __expf(src[i] - m);
  s = warpSum(s);
  if (lane == 0) part[wid] = s;
  __syncthreads();
  if (tid == 0) {
    float t = 0.f;
    for (int w = 0; w < nWarps; w++) t += part[w];
    bcast = t;
  }
  __syncthreads();
  const float inv = 1.f / bcast;

  // Pass 3: normalize and write
  for (int i = tid; i < nCols; i += BLOCK) dst[i] = __expf(src[i] - m) * inv;
}

int main() {
  const int R = 8192, C = 4096;
  size_t bytes = (size_t)R * C * sizeof(float);
  float *h = (float *)malloc(bytes);
  for (long i = 0; i < (long)R * C; i++) h[i] = (float)((i * 37) % 100) * 0.01f;

  float *d_in, *d_out;
  cudaMalloc(&d_in, bytes); cudaMalloc(&d_out, bytes);
  cudaMemcpy(d_in, h, bytes, cudaMemcpyHostToDevice);

  for (int i = 0; i < 5; i++) softmaxKernel<<<R, BLOCK>>>(d_out, d_in, C);
  cudaDeviceSynchronize();

  cudaEvent_t a, b; cudaEventCreate(&a); cudaEventCreate(&b);
  cudaEventRecord(a);
  for (int i = 0; i < 50; i++) softmaxKernel<<<R, BLOCK>>>(d_out, d_in, C);
  cudaEventRecord(b); cudaEventSynchronize(b);

  float ms; cudaEventElapsedTime(&ms, a, b);
  double per = ms / 50.0;
  // Minimum traffic: one read, one write. (3 passes, but the cache absorbs it)
  printf("cuda   %7.3f ms  %7.1f GB/s\n", per,
         2.0 * bytes / (per * 1e-3) / 1e9);
  return 0;
}

Of those 80 lines, the actual math is three lines. Everything else is plumbing for "gathering values inside the block and broadcasting them back out." And that plumbing hides a quiet trap.

It's the __syncthreads() right after m = bcast; in the code above. Without it, a fast warp can overwrite pass 2's part[wid] before a slow warp has read pass 1's result out of it, wiping it out. The result comes out correct most of the time and wrong occasionally. On small inputs it almost never reproduces. This is the most expensive kind of bug in a hand-written kernel, and the odds go up the longer this plumbing gets.

Triton — Write Tile by Tile and Leave the Threads to the Compiler

Now let's write the same operation in Triton.

# softmax_triton.py
# Run: python softmax_triton.py   (based on triton 3.7.x, torch 2.13)
import torch
import triton
import triton.language as tl


@triton.jit
def softmax_kernel(out_ptr, in_ptr, in_stride, out_stride, n_cols,
                   BLOCK_SIZE: tl.constexpr):
    row = tl.program_id(0)
    cols = tl.arange(0, BLOCK_SIZE)
    mask = cols < n_cols

    # Load the whole row into SRAM at once. Padding slots are -inf, so they don't affect max/sum.
    x = tl.load(in_ptr + row * in_stride + cols, mask=mask, other=-float("inf"))
    x = x - tl.max(x, axis=0)
    num = tl.exp(x)
    y = num / tl.sum(num, axis=0)
    tl.store(out_ptr + row * out_stride + cols, y, mask=mask)


def softmax(x: torch.Tensor) -> torch.Tensor:
    n_rows, n_cols = x.shape
    BLOCK_SIZE = triton.next_power_of_2(n_cols)
    num_warps = 4 if BLOCK_SIZE < 2048 else (8 if BLOCK_SIZE < 8192 else 16)
    out = torch.empty_like(x)
    softmax_kernel[(n_rows,)](
        out, x, x.stride(0), out.stride(0), n_cols,
        BLOCK_SIZE=BLOCK_SIZE, num_warps=num_warps,
    )
    return out


if __name__ == "__main__":
    torch.manual_seed(0)
    x = torch.randn(8192, 4096, device="cuda", dtype=torch.float32)

    ours, ref = softmax(x), torch.softmax(x, axis=1)
    assert torch.allclose(ours, ref, atol=1e-5), "correctness check failed"

    gb = 2 * x.numel() * x.element_size() / 1e9
    for name, fn in [
        ("triton", lambda: softmax(x)),
        ("torch ", lambda: torch.softmax(x, axis=1)),
    ]:
        ms = triton.testing.do_bench(fn, warmup=25, rep=100)
        print(f"{name}  {ms:7.3f} ms  {gb / (ms * 1e-3):7.1f} GB/s")

The key difference is a single line: tl.max(x, axis=0). The whole-block reduction that took 40 lines in CUDA is a single array operation here. Whether to use warp shuffles, how many bytes of shared memory to allocate, where to place synchronization — the compiler decides all of it. And the synchronization-omission bug described above cannot happen, structurally. We never write a synchronization call in the first place.

The Compile Pipeline — What the Compiler Actually Does

Triton isn't magic. There are stages. Here's the path you can verify from the repository's structure.

Python function (@triton.jit)
   │  Walks the Python AST to generate IR
Triton IR (TTIR)          Tile operations. No hardware concepts yet
   │  Layout assignment, coalescing, pipelining passes
TritonGPU IR (TTGIR)      Warp/thread layout and shared memory are decided
LLVM IR
   │  NVPTX backend          │  AMDGPU backend
   ▼                        ▼
PTX → (ptxas) → SASS       AMDGCN → object

Two practically meaningful facts fall out of this.

First, the backends split below LLVM IR. The Triton repository's third_party directory holds the nvidia and amd backends side by side. That means the same @triton.jit kernel compiles for both vendors. This property is what the next two sections rest on.

Second, the step from TTIR down to TTGIR decides most of the performance. Get the layout assignment wrong and shared-memory accesses trigger bank conflicts; skip the pipelining pass and load latency stays fully exposed. The levers we can actually touch amount to BLOCK_SIZE, num_warps, and num_stages, so Triton tuning is mostly a matter of turning those three knobs.

Triton exposes that knob-turning as a decorator.

@triton.autotune(
    configs=[
        triton.Config({"BLOCK_M": 64,  "BLOCK_N": 64},  num_warps=4, num_stages=3),
        triton.Config({"BLOCK_M": 128, "BLOCK_N": 64},  num_warps=8, num_stages=3),
        triton.Config({"BLOCK_M": 128, "BLOCK_N": 128}, num_warps=8, num_stages=4),
        triton.Config({"BLOCK_M": 64,  "BLOCK_N": 128}, num_warps=4, num_stages=4),
    ],
    key=["M", "N", "K"],          # retunes whenever these values change
)
@triton.jit
def matmul_kernel(...):
    ...

Every time an argument listed in key changes, it compiles every candidate, measures them, and caches the fastest one. Doing the same thing in CUDA C++ means writing your own template instantiation and benchmark harness by hand.

When You Want to See What Got Built

There's a way to keep Triton from being a black box. You can pull the intermediate artifacts out directly.

import triton, torch

# Compile the kernel once to get a handle
k = softmax_kernel.warmup(
    torch.empty(1, 8, device="cuda"), torch.empty(1, 8, device="cuda"),
    8, 8, 8, BLOCK_SIZE=8, grid=(1,),
)
k._init_handles()

print(k.asm.keys())        # ttir, ttgir, llir, ptx, cubin
print(k.asm["ttgir"][:800])  # Eyeball the layout assignment result
print(k.n_regs, k.n_spills)  # Nonzero here means register spills occurred

If n_spills isn't zero, it means registers ran out and spilled to local memory, and that usually means BLOCK_SIZE is too large. It's the first value to check when a Triton kernel is inexplicably slow.

Where Triton Falls Short

This part deserves an honest accounting. There are clearly places you shouldn't reach for Triton.

  • Kernels that need precise, warp-level control. Algorithms that depend on shuffles between specific lanes, warp specialization (splitting producer warps from consumer warps), or __ballot-style voting operations don't fit into the tile abstraction. To close this gap, Triton has recently been experimenting with a lower-level language called triton.experimental.gluon, but as of 3.7.1 it's still experimental, so it's too early to treat it as a basis for production decisions.
  • Irregular indexing and data-dependent control flow. When the access pattern is only decided at runtime — graph traversal, sorting, dynamically sized sparse operations — the tile model doesn't fit.
  • The absolute-best-performance GEMM. Dense matrix multiplication is territory where CUTLASS and cuBLAS have accumulated years of tuning, and it's hard for Triton to close that last 10 to 20 percent. In practice there's no reason to try, either.
  • When compile time counts against response time. JIT compilation is tacked onto the first call. Turn on autotuning and that cost multiplies by the number of candidates. On a serving path, cache warming is a precondition.
  • When exact numerical reproducibility is required. The compiler decides the reduction order, so the last bit can change across versions.

Comparing the Same Kernel Written Twice

Here's how to measure the two implementations under the same conditions, and what shape the results take.

The measurement rules are the same as in the previous post: discard the warmup, run it many times and average, verify correctness first, and convert to effective bandwidth rather than raw time. On the Triton side, triton.testing.do_bench already handles warmup, repetition, and even L2 cache flushing, so there's no need to write it yourself. On the CUDA side, we match the same conditions with cudaEvent.

softmax is a memory-bound kernel with low arithmetic intensity, so both implementations ultimately run up against HBM bandwidth. So the results tend to come out looking roughly like this.

ImplementationCode size (kernel body)Effective bandwidthNotes
CUDA C++ (code above)~45 lines80% range of the ceiling3 passes. Mostly reduction plumbing
Triton (code above)12 lines85% range of the ceiling1 pass once the row fits in SRAM
torch.softmax0 lines85% range of the ceilingAlready a fused library kernel

The absolute numbers vary a lot with the GPU and the row length. Run the harness above on your own hardware and produce your own numbers. What matters in this table isn't the ranking, but two observations.

First, the code size differs by almost 4x, yet the performance is comparable. That's because, in a memory-bound kernel, both implementations run into the same wall. If the wall decides performance, then the shorter the code that gets you to the wall, the better.

Second, torch.softmax is already just as fast. That's the real lesson of this comparison. Don't hand-write standard operators. The value of writing a kernel yourself comes from operations the library doesn't have — fused combinations, or domain-specific masking. For example, if masking, scaling, and dropout are bolted onto both sides of a softmax and each one is running as a separate kernel, merging them into one drops the HBM round trips from four to one. The gain there is 3x or more, and that's the reason to reach for Triton.

Why Triton Became the De Facto Standard for Custom Attention

These days, when a new attention-variant kernel shows up, it's almost always Triton first. The reason isn't performance — it's structure.

  • Attention fits the tile model well. Taking one query block and sweeping keys and values block by block while accumulating with online softmax is, itself, Triton's programming model.
  • There are many variants and they're short-lived. Sliding window, ALiBi, soft-capping, and various sparse patterns keep appearing, and only some of them survive. You can't spend two weeks in CUDA C++ on a single experiment.
  • Write it once and it runs on both vendors. As we saw earlier, the backends split below LLVM IR, so the same kernel compiles for both NVIDIA and AMD. Written in CUDA C++, the AMD side is a separate effort.
  • The ecosystem already assumes Triton. PyTorch's TorchInductor generates Triton kernels as its GPU-targeted code. vLLM keeps TRITON_ATTN in its list of attention backends and on its priority table (vLLM attention backend docs). The simple fact that it's already installed and already running is a major advantage in itself.

CUTLASS and CuTe — Assembling GEMM from Templates

This is the third layer. CUTLASS is a CUDA C++ template library dating back to 2017 that hierarchically decomposes matrix multiplication and its surrounding operations into reassemblable parts.

The core idea is to stop treating GEMM as one monolithic block. You split the whole problem into thread-block tiles, split those again into warp tiles, and split those again down to the size a single Tensor Core instruction can handle. Each layer is an independent template parameter, and we pick the combination.

// Skeleton of a CUTLASS 3.x/4.x-style GEMM configuration
// See the examples directory in the repository for a full, buildable example.
using ElementA = cutlass::half_t;
using ElementB = cutlass::half_t;
using ElementC = float;

// What tile size, processed with how many pipeline stages
using TileShape   = cute::Shape<cute::_128, cute::_128, cute::_64>;
using ClusterShape = cute::Shape<cute::_1, cute::_1, cute::_1>;

// Assemble the main loop (data movement + Tensor Core accumulation)
// and the epilogue (post-processing applied to the result) separately
using CollectiveMainloop = /* CollectiveBuilder<...> */;
using CollectiveEpilogue = /* CollectiveBuilder<...> */;

using GemmKernel = cutlass::gemm::kernel::GemmUniversal<
    cute::Shape<int, int, int, int>, CollectiveMainloop, CollectiveEpilogue>;
using Gemm = cutlass::gemm::device::GemmUniversalAdapter<GemmKernel>;

The single most practically valuable concept here is the epilogue. Apply the bias addition and the activation function while the GEMM result is still sitting in the accumulator registers, and the entire round trip of writing the result to HBM and reading it back disappears. Call cuBLAS and apply the activation as a separate kernel, and that round trip stays right where it was. A good part of the reason to reach for CUTLASS directly lives here.

CuTe and Layout Algebra

The foundation from CUTLASS 3.x onward is CuTe. The idea is to make a tensor's layout — the mapping from logical coordinates to physical addresses — a first-class value, and to combine it algebraically.

In the previous post, we solved a shared-memory bank conflict by padding the array by one slot. In CuTe, that same kind of maneuver is expressed as a layout transform called a "swizzle," can be decomposed and composed into tiles, and gets checked at compile time. The kind of bug you get from writing index arithmetic by hand and getting it wrong is caught at the type level instead.

Honestly, the learning cost is high. Template error messages are long, and the concepts stack up in layers. That's why CUTLASS 4 added CuTe DSL, a Python interface. The CUTLASS 4.6.1 README describes it as "a low-level programming model that exposes CuTe's core concepts — layout, tensor, hardware atoms, full control over the thread and data hierarchy — without C++," and explicitly states that it is currently in public beta. It can be installed from PyPI as nvidia-cutlass-dsl, and 4.6.1 is the latest version.

It's worth simply taking the "beta" label at face value. It's usable for prototyping, but it's still too early to make it your only implementation of a production kernel.

Which Layer to Choose

Here's the judgment call compressed into a table.

SituationRecommendationReason
Standard GEMM, convolutioncuBLAS / cuDNNNo reason to hand-write it
Want to attach post-processing to a GEMM to remove a round tripCUTLASS epilogueThe layer designed for exactly this
Attention variants, fused elementwise chainsTritonFits the tile model, iterates fast
MoE routing, custom normalizationTritonSame
Algorithm is fundamentally warp specialization, lane-level shufflesCUDA C++Can't be expressed via tile abstraction
Irregular indexing, data-dependent control flowCUDA C++Same
Must use a new hardware instruction firstCUDA C++ or CUTLASSCompiler support arrives late
Must support NVIDIA and AMD from one codebaseTritonBackends split further down

Written out as a practical order of operations: check whether a library already does it. If not, write it in Triton. If Triton can't express it, or the performance isn't there, then drop down to CUDA C++, and only for the specific kernel the profiler points at. Starting from CUDA C++ from the very beginning only makes sense when you're already certain the kernel inherently demands warp-level control.

One more thing worth adding: mixing layers is normal. Open up an actual inference stack and you'll find GEMM handled by cuBLAS or CUTLASS, attention hand-written in CUDA C++ or in Triton, normalization and activations in Triton, and PyTorch gluing it all together. There's no real basis for the pressure to pick one layer and standardize on it everywhere.

Closing — The Layer of Abstraction Is Chosen for Iteration Speed, Not Performance

We've compared three layers, but the conclusion isn't "Triton won." CUDA C++ and Triton coming out with comparable performance on a memory-bound kernel is because both run into the same physical wall, and CUTLASS pulling ahead on GEMM is because that wall sits somewhere else. The layer itself doesn't decide the performance.

What the layer decides is something else: how long it takes to fix a kernel and re-measure it, the odds of a nondeterministic bug from a missed synchronization, how much code needs revisiting when a new GPU ships. When performance is comparable, these should be the deciding criteria — and in most real-world work, performance is comparable.

So the judgment call boils down to one line: how many more times am I likely to touch this kernel? If it's a kernel you'll write once and leave untouched for two years, dropping down to CUDA C++ to squeeze out the last 10 percent is reasonable. If it's a kernel the experiments will keep touching next month too, writing it in 12 lines and spending that saved time running one more experiment is almost always the better call.

References