Skip to content
Published on

What It Really Means to Hand-Tune a GPU Kernel — Making One Transpose Kernel 5x Faster

Share
Authors

Introduction — When the Kernel the Profiler Points At Is Not in Any Library

You hit a training step that takes 420ms, and 170ms of it sits inside one kernel with no recognizable name — and suddenly your options shrink. If it's a cuBLAS call, there's nothing to touch. If it's a PyTorch operator, you can try swapping in another operator. But if that kernel is something someone hastily wrote to handle your model's unusual masking logic, you're the only one who can fix it.

This post starts from that point. It covers what "hand-tuning" a GPU kernel actually involves — what you need to know, and what you need to measure. The short version: eight-tenths of kernel optimization is not about cutting compute, it's about changing the order in which memory moves. The kernel we fix in this post's hands-on section doesn't contain a single arithmetic operation. And yet it still gets more than 5x faster.

The reference environment is CUDA Toolkit 13.3 Update 1 and the Nsight Compute 2026.2 line. The concepts don't age with hardware generations, but tool flag names do change, so check your command lines against your own version's Nsight Compute CLI documentation.

What a Kernel Is, and Why You End Up Modifying One

A kernel is a single function that runs on the GPU. The difference from a CPU function is that one call launches tens of thousands of instances simultaneously. Each instance is a thread, and a thread figures out its own index from a built-in variable and processes only the data slice that corresponds to it.

Situations where hand-modifying a kernel is actually worthwhile are narrower in practice than you'd think. Checked in order, it looks like this.

SituationWhat to do firstIs it a reason to modify the kernel?
Standard GEMM or convolution is slowCheck cuBLAS/cuDNN version, data type, tensor-core pathAlmost never
Several small operators run back to backCoax fusion out of torch.compileUsually not
You need an attention variantCheck whether a FlashAttention-family library already has that variantIf not, yes
One kernel uses only 20 percent of theoretical bandwidthFind the cause with Nsight ComputeYes
Indexing, masking, or sparsity patterns unique to your domainNo library to substituteYes

The key is the last two rows. If a library can replace it, replacing it is always better. Library kernels carry far more tuning than you could ever spend the time on, and they get updated for you when a new architecture ships. A kernel you wrote by hand is one you have to maintain for the rest of its life.

The Execution Model — Threads, Warps, Blocks, Grids

To modify a kernel you need to know how the hardware groups threads. There are four levels of hierarchy.

Grid          The entire set of threads for one kernel launch
  └ Block     Placed together on the same SM, shares shared memory and __syncthreads()
      └ Warp  32 threads. The real unit of scheduling and instruction issue
          └ Thread   A single execution flow with its own registers

The single most practically important fact here is that the warp is the real unit. Programmers write code per-thread, but the hardware issues 32 of them as one bundle. Two consequences follow from this.

First, branch divergence. If threads within one warp take different branches, the hardware executes both paths sequentially, disabling whichever threads don't apply to the current path. A conditional that splits inside a warp adds to execution time. A conditional that splits cleanly at warp boundaries is free.

Second, memory access is also coalesced at the warp level. That's coalescing, the subject of the next section.

The default when choosing a block size is to make it a multiple of 32. If the block size is 33, the hardware issues two warps and 31 lanes of the second one sit idle.

The Memory Hierarchy — Coalescing and Bank Conflicts

GPU memory comes in layers, and latency and bandwidth differ by orders of magnitude between layers.

LayerScopeRough character
RegistersOne threadFastest. Count per thread governs occupancy
Shared memoryOne blockSRAM inside the SM. A programmer-managed cache
L1 / texture cacheOne SMPhysically shares the same storage as shared memory
L2 cacheWhole GPUShared by every SM. The last line of defense before HBM
HBM (global memory)Whole GPULarge capacity, but latency runs into the hundreds of cycles

Coalescing

Global memory accesses are handled as 32-byte transactions. If the 32 threads of one warp read 32 consecutive floats — 128 bytes — it finishes in four transactions. Conversely, if the same 32 threads read 32 floats spaced 4096 apart, it needs 32 transactions, and each one fetches 32 bytes just to use 4 and discard the rest. That's using only one-eighth of the bandwidth.

This is the single largest cause of kernel performance differences. In the hands-on section you'll see exactly this 8x waste.

Bank conflicts

Shared memory is interleaved across 32 banks. On a 4-byte-word basis, the bank number is the address divided by 4, then taken modulo 32. If the threads of one warp touch different banks, it's handled in one cycle; if they touch different addresses within the same bank, accesses get serialized by that many.

The classic accident is column-wise access into a square shared array. If you sweep tile[i][0] over i in a tile[32][32] array, every element lands on bank 0 — a 32-way conflict. The fix is to widen the array by one column. Declaring tile[32][33] shifts each row's starting bank by one, so column access spreads evenly across all 32 banks. It's a trade: 32 * 4 bytes, i.e. 128 extra bytes per block of shared memory, in exchange for removing the serialization.

Occupancy and the Roofline — What to Look At and What to Ignore

Occupancy is a symptom, not a goal

Occupancy is the number of warps actually resident on an SM divided by the maximum it can host. A common mistake for beginners is treating this as a target to maximize. It isn't.

Occupancy does exactly one job: hiding latency. If one warp is waiting on an HBM response and another warp runs in the meantime, the SM never sits idle. So occupancy is a proxy for the question "do we have enough warps to hide latency" — it is not performance in itself.

There are two cases where lower occupancy is actually faster.

First, when a thread uses many registers to extract instruction-level parallelism. If a single thread has four independent loads in flight at once, the total volume of in-flight memory requests stays the same even with a quarter as many warps. Vasily Volkov's Better Performance at Lower Occupancy addresses this thesis head-on. It's from 2010, but the argument still holds today.

Second, when the kernel is already pinned to the bandwidth ceiling. If the memory pipe is saturated, adding more warps has nowhere to go.

So the practical rule is: check occupancy only when it's low. If it's dropped below 25 percent and the kernel is latency-bound, suspect register usage or shared-memory allocation. If it reads 60 percent and the kernel is still slow, occupancy isn't the culprit — look elsewhere.

Roofline — most kernels are bandwidth-bound

The tool for deciding where to look is the roofline. Its axis is arithmetic intensity — operations performed per byte moved.

Achieved performance (FLOP/s)
   ^
   |            ______________  compute ceiling
   |           /
   |          /   slope = memory bandwidth
   |         /
   +--------+-------------------> arithmetic intensity (FLOP/Byte)
          ridge point

ridge point = (compute ceiling FLOP/s) / (memory bandwidth Byte/s)

The ridge point is a hardware property. On modern datacenter GPUs this value falls somewhere in the tens to hundreds of FLOP/Byte. But when you compute the arithmetic intensity of the kernels we actually use, most of them are single digits.

OperationRough arithmetic intensityWhere it lands
Elementwise addition1 FLOP / 12 ByteExtremely memory-bound
Activation functiona few FLOP / 8 ByteMemory-bound
Matrix transpose0 FLOP / 8 BytePure memory
LayerNormlow single-digit FLOP / ByteMemory-bound
GEMM (large matrices)up to hundreds, proportional to tile sizeCompute-bound
LLM decode stepunder 2 when batch is smallMemory-bound

Reading it is simple. If arithmetic intensity sits far below the ridge point, that kernel's performance ceiling is already fixed. Cutting compute instructions further does nothing; the only improvement that matters is reducing the bytes moved, or fixing how they move.

So the first question in kernel optimization is always: "what is this kernel's theoretical minimum traffic, and how many bytes is it actually moving right now?"

Hands-On — Fixing a Transpose Kernel in Four Stages

Now let's actually fix one. The target is a 4096 x 4096 float matrix transpose. Compute count is zero, so only the memory story remains — which makes it an ideal teaching example.

The theoretical minimum traffic is clear. One read and one write means 2 * 4096 * 4096 * 4 bytes, about 134MB. There's no way to move less than that. So the performance ceiling is "a kernel that just copies without transposing," and we measure that first as our baseline.

Full code

// transpose.cu
// build: nvcc -O3 -arch=sm_80 transpose.cu -o transpose
#include <cstdio>
#include <cstdlib>
#include <cuda_runtime.h>

static const int TILE = 32;
static const int BLOCK_ROWS = 8;   // 32x8 = 256 threads per block
static const int N = 4096;

#define CHECK(x) do { cudaError_t e_ = (x); if (e_ != cudaSuccess) { \
    printf("CUDA error: %s (line %d)\n", cudaGetErrorString(e_), __LINE__); \
    exit(1); } } while (0)

// Stage 0. Upper bound: copy only, no transpose.
__global__ void copyKernel(float *out, const float *in) {
  int x = blockIdx.x * TILE + threadIdx.x;
  int y = blockIdx.y * TILE + threadIdx.y;
  for (int j = 0; j < TILE; j += BLOCK_ROWS)
    out[(y + j) * N + x] = in[(y + j) * N + x];
}

// Stage 1. naive: the read is coalesced, but the write strides by N.
__global__ void transposeNaive(float *out, const float *in) {
  int x = blockIdx.x * TILE + threadIdx.x;
  int y = blockIdx.y * TILE + threadIdx.y;
  for (int j = 0; j < TILE; j += BLOCK_ROWS)
    out[x * N + (y + j)] = in[(y + j) * N + x];
}

// Stage 2. shared-memory tile: finish the transpose inside SRAM,
//          so both the global read and the global write are coalesced.
__global__ void transposeShared(float *out, const float *in) {
  __shared__ float tile[TILE][TILE];

  int x = blockIdx.x * TILE + threadIdx.x;
  int y = blockIdx.y * TILE + threadIdx.y;
  for (int j = 0; j < TILE; j += BLOCK_ROWS)
    tile[threadIdx.y + j][threadIdx.x] = in[(y + j) * N + x];

  __syncthreads();

  // Swap the block coordinates so the write also lands on consecutive addresses.
  x = blockIdx.y * TILE + threadIdx.x;
  y = blockIdx.x * TILE + threadIdx.y;
  for (int j = 0; j < TILE; j += BLOCK_ROWS)
    out[(y + j) * N + x] = tile[threadIdx.x][threadIdx.y + j];
}

// Stage 3. remove the shared-memory bank conflict with one column of padding.
__global__ void transposePadded(float *out, const float *in) {
  __shared__ float tile[TILE][TILE + 1];   // the only difference

  int x = blockIdx.x * TILE + threadIdx.x;
  int y = blockIdx.y * TILE + threadIdx.y;
  for (int j = 0; j < TILE; j += BLOCK_ROWS)
    tile[threadIdx.y + j][threadIdx.x] = in[(y + j) * N + x];

  __syncthreads();

  x = blockIdx.y * TILE + threadIdx.x;
  y = blockIdx.x * TILE + threadIdx.y;
  for (int j = 0; j < TILE; j += BLOCK_ROWS)
    out[(y + j) * N + x] = tile[threadIdx.x][threadIdx.y + j];
}

typedef void (*Kern)(float *, const float *);

static void bench(const char *name, Kern k, float *d_out, const float *d_in,
                  const float *h_ref, float *h_out, bool checkTranspose) {
  dim3 grid(N / TILE, N / TILE), block(TILE, BLOCK_ROWS);
  const int WARMUP = 5, ITERS = 50;
  const double bytes = 2.0 * N * N * sizeof(float);

  for (int i = 0; i < WARMUP; i++) k<<<grid, block>>>(d_out, d_in);
  CHECK(cudaDeviceSynchronize());

  cudaEvent_t t0, t1;
  CHECK(cudaEventCreate(&t0));
  CHECK(cudaEventCreate(&t1));
  CHECK(cudaEventRecord(t0));
  for (int i = 0; i < ITERS; i++) k<<<grid, block>>>(d_out, d_in);
  CHECK(cudaEventRecord(t1));
  CHECK(cudaEventSynchronize(t1));

  float ms = 0.f;
  CHECK(cudaEventElapsedTime(&ms, t0, t1));
  double perIter = ms / ITERS;
  double gbs = bytes / (perIter * 1.0e-3) / 1.0e9;

  // A performance number with no correctness check is meaningless.
  CHECK(cudaMemcpy(h_out, d_out, (size_t)N * N * sizeof(float),
                   cudaMemcpyDeviceToHost));
  long bad = 0;
  for (long r = 0; r < N && bad == 0; r++)
    for (long c = 0; c < N; c++) {
      float want = checkTranspose ? h_ref[c * N + r] : h_ref[r * N + c];
      if (h_out[r * N + c] != want) { bad++; break; }
    }

  printf("%-18s %8.3f ms   %8.1f GB/s   %s\n", name, perIter, gbs,
         bad ? "FAIL" : "ok");
  CHECK(cudaEventDestroy(t0));
  CHECK(cudaEventDestroy(t1));
}

int main() {
  size_t bytes = (size_t)N * N * sizeof(float);
  float *h_in = (float *)malloc(bytes), *h_out = (float *)malloc(bytes);
  for (long i = 0; i < (long)N * N; i++) h_in[i] = (float)(i % 1000);

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

  cudaDeviceProp p;
  CHECK(cudaGetDeviceProperties(&p, 0));
  printf("%s  peak HBM = %.1f GB/s\n\n", p.name,
         2.0 * p.memoryClockRate * (p.memoryBusWidth / 8) / 1.0e6);

  bench("copy (upper bound)", copyKernel,      d_out, d_in, h_in, h_out, false);
  bench("naive",              transposeNaive,  d_out, d_in, h_in, h_out, true);
  bench("shared tile",        transposeShared, d_out, d_in, h_in, h_out, true);
  bench("shared + padding",   transposePadded, d_out, d_in, h_in, h_out, true);

  cudaFree(d_in); cudaFree(d_out); free(h_in); free(h_out);
  return 0;
}

What matters about the measurement method

For the numbers to be trustworthy, the harness has to be honest first. The code above follows five rules.

  • Discard the warmup. The first call has context creation and module loading mixed into it.
  • Average over repeated measurements. If a single kernel runs around 1ms, clock variation alone can swing the result by 10 percent.
  • Use cudaEvent, not a CPU timer. Kernel launches are asynchronous, so a CPU timer only measures how long it took to submit.
  • Verify correctness. A broken index is usually faster. A GB/s number with no verification is just a number game.
  • Convert to effective bandwidth, not raw time. Absolute time depends on size and hardware, but the percentage of theoretical bandwidth is comparable anywhere.

The effective-bandwidth formula is simple: divide the minimum bytes that must move by the time it actually took. The "must" is the important word here — use the bytes the algorithm actually requires, not however many bytes were wasted in transit. That's what makes the waste show up in the number.

The shape of the results

Below is the typical shape you get running this harness on an A100 80GB (sm_80)-class machine. Absolute values vary a lot by device, driver, and clock state, so don't quote these numbers as-is — run the code above on your own GPU and build your own baseline. What matters is the relative ratio between stages.

StageEffective bandwidthVs. copyBottleneck
copy (ceiling)baseline 100100 percentNone. HBM saturated
naiveabout 18about 18 percentNon-coalesced write. Only 4 bytes used per transaction
shared tileabout 63about 63 percent32-way shared-memory bank conflict
shared + paddingabout 93about 93 percentEssentially none. Only tile-boundary effects remain

Three things are worth reading here.

First, naive sits at about a fifth of the ceiling. Zero compute, the same amount of data moved, and yet it's 5x slower. The only difference is order. The write is scattered at N-wide strides, so out of every 32-byte transaction only 4 bytes get used and 28 are discarded. An 8x waste shows up mixed with other effects as a 5x difference.

Second, the shared tile recovers most of the gap but doesn't get all the way there. Both global accesses are now fixed, but the bottleneck moved inside SRAM. tile[threadIdx.x][threadIdx.y + j] is a column-wise access, and in a 32x32 square array every column access lands on the same bank.

Third, the code difference in the final stage is one character in an array declaration. Changing [TILE] to [TILE + 1] is the whole fix. Kernel optimization often looks exactly like this — not an algorithmic change, just shifting the data layout by one slot.

Common failure modes

Landmines actually stepped on often in this exercise.

  • Missing __syncthreads(). Without a sync between filling the shared tile and reading it back, results go wrong non-deterministically. Small inputs happen to come out right often enough to be more dangerous.
  • Putting __syncthreads() inside a branch. Placed where only some threads in a block reach it, this is undefined behavior.
  • Forgetting to swap block coordinates in the index math. If you add the shared tile but leave the output index alone, the write becomes strided again and stage 2's gain disappears. The result is correct but not faster — the hardest form of this bug to notice.
  • Measuring without -O3. Without host-code optimization, the verification loop dominates the measured time and flips the conclusion.
  • N too small. Kernel launch overhead is a few microseconds; if the total time is tens of microseconds, you're measuring overhead.

What Nsight Compute Shows You

The profiler is what tells you the cause when a number comes back bad. Let's start with the command line.

# Collect every section. One kernel can take hundreds of ms, so narrow the target.
ncu --set full \
    --kernel-name regex:transpose \
    --launch-skip 5 --launch-count 1 \
    -o transpose_report \
    ./transpose

# Pull only the metrics you need to pin down the cause (much faster)
ncu --metrics \
  sm__throughput.avg.pct_of_peak_sustained_elapsed,\
gpu__dram_throughput.avg.pct_of_peak_sustained_elapsed,\
l1tex__data_bank_conflicts_pipe_lsu_mem_shared.sum,\
l1tex__average_t_sectors_per_request_pipe_lsu_mem_global_op_ld.ratio \
  --kernel-name regex:transpose --launch-count 1 ./transpose

# Open in the GUI
ncu-ui transpose_report.ncu-rep

Skipping the warmup runs with --launch-skip matters. The profiler shows you the cold-cache state of the first run exactly as it is, so without skipping you end up analyzing something that isn't steady state.

Opening the report gives you several sections, and in practice there's a fixed order to go through them.

1. Speed of Light. Shows compute throughput and memory throughput as a percentage of hardware theoretical maximum. This is where direction gets decided. If memory is above 80 percent, you're pinned to bandwidth; if both are under 30 percent, it's a latency or occupancy problem. Our naive transpose kernel comes back low on both in this view — not because the pipe is saturated, but because it's wasting.

2. Memory Workload Analysis. This is the crux. It shows how many sectors were fetched per request; a perfectly coalesced 32-thread float load is 4 sectors per request. The naive kernel's write comes back at 32 sectors per request. The 8x waste shows up in this one line exactly as it is. This single metric answers "is it a coalescing problem" directly.

3. Shared Memory metrics. Shows the bank-conflict count. Stage 2's kernel spikes here, and stage 3 drops it to near zero. This is where you confirm the padding actually worked.

4. Warp State Statistics. Shows, by category, the reasons warps were stalled. If Stall Long Scoreboard dominates, that's waiting on a global memory response; if it's Stall MIO Throttle, that's congestion on shared memory or special-function units. It splits the cause into memory versus instructions.

5. Occupancy. Look at this last. When the first four are clean and it's still slow, that's when you check whether resident warps are too few. Skip this order and look at occupancy first, and you usually end up optimizing in the wrong direction.

Nsight Compute also has a roofline section, which shows in a chart whether our kernel sits on the slope or on the flat part. The transpose kernel has zero arithmetic intensity, so it plots at the far left edge — and that alone tells you "there is no compute to optimize here."

When to stop this work

Once you've pushed the transpose kernel to 93 percent of the ceiling, there's little reason to chase the remaining 7 percent. It's worth having a stopping rule set in advance.

  • Stop once you've crossed 90 percent of theoretical minimum traffic. In a memory-bound kernel, anything above that is tile-boundary and TLB effects, where the return on effort drops off sharply.
  • Re-measure this kernel's share of total runtime. If you took it from 170ms to 40ms, the bottleneck is now somewhere else. Amdahl's law holds just as much for kernel optimization.
  • Factor in maintenance cost. A hand-written kernel is up for revalidation every time a new architecture ships. A kernel that's 20 percent faster than the library today can be 30 percent slower two years from now.
  • Check first whether it's solved one layer up. If the next post's Triton version of the same kernel gets written in 20 lines and lands at similar performance, that's much less reason to keep maintaining the CUDA C++ version.

Closing — Fixing a Kernel Means Fixing the Order Data Moves In

Not a single floating-point operation appeared in this post's hands-on section. And yet there was a 5x gap between the first version and the last. What changed was only this: in what order the same data was read, where it briefly sat, and in what order it was written back.

That fact defines the character of GPU kernel work. Algorithmic improvements that cut compute are usually already done by a library, or simply not available to change in our problem. What's left as a lever is data placement and movement order within the memory hierarchy — and fortunately, that's the much bigger lever.

The work order boils down to one line: compute the theoretical minimum traffic, measure with a harness what percentage you're currently using, pinpoint the location of the waste with Nsight Compute, fix the layout, and measure again. Fixing by feel and calling it faster skips two measurements out of this sequence, and a conclusion reached that way flips on the next machine.

References