Optimizing a CUDA Matrix Multiply Step by Step, From Naive to Tensor Cores

A fast CUDA matrix multiply is built in layers: make each warp read contiguous memory (coalescing), stage tiles in shared memory so each value is reused, keep a block of results in registers so shared memory stops being the bottleneck, then hand the inner product to tensor cores. Each step raises arithmetic intensity.

This post walks through that ladder one kernel at a time, at the level a GPU performance engineer would check line by line. For every stage you get the code, the arithmetic intensity maths, the bottleneck you should expect, and what Nsight Compute shows when you profile it. The published measurements quoted here come from two places: the NVIDIA CUDA C++ Best Practices Guide, which reports effective bandwidth for its shared memory examples on a Tesla V100, and Simon Boehm's December 2022 worklog, which measured ten SGEMM kernels against cuBLAS on an RTX A6000. Neither set was produced by us.

One plain statement before the code: we could not run any of this code. This post was written on a Mac with no NVIDIA GPU. The kernels follow the CUDA C++ Programming Guide, the PTX ISA and the cuBLAS documentation, and every number we computed ourselves (FLOP counts, intensities, occupancy, roofline ridge points) was checked with Python. No benchmark result in this post is ours.

If you want the hardware background first, our explainers on AI chip architecture, what an AI accelerator is and how a GPU is made, from wafer to HBM package cover it. For what these kernels cost to run at scale, see our GPU cost per million tokens estimate.

Key takeaways

• A square FP32 matrix multiply of size 4096 does about 137.4 billion floating point operations but only needs to move about 201 MB at minimum (268 MB if C is also read), so it is compute bound in principle. Every naive kernel turns it into a memory bound one by moving the same bytes many times.
• Coalescing is the first and cheapest win. The CUDA Programming Guide describes global memory as served in 32 byte transactions; a warp reading 32 consecutive floats needs 4 of them, a warp reading 32 floats with a large stride needs 32, and only 12.5% of the bytes fetched are used.
• Shared memory tiling with a tile of side T raises global arithmetic intensity to T/4 FLOP per byte for FP32, but it moves the bottleneck to shared memory bandwidth unless each thread computes more than one output.
• Register blocking (1D, then 2D) is what turns shared memory loads into many fused multiply adds per load. In Boehm's A6000 measurements, 1D blocktiling took the kernel from 2,980 to 8,475 GFLOP/s and 2D blocktiling to 15,972 GFLOP/s.
• Vectorized loads, double buffering with asynchronous copies, and warp tiling close most of the remaining gap. Boehm's warptiled kernel reached 21,779 GFLOP/s, 93.7% of cuBLAS on the same card, which itself ran at about 60% of the A6000's 38.7 TFLOPS FP32 peak.
• Tensor cores change the scale of the problem. The A6000's dense FP16 tensor peak with FP32 accumulate is 154.8 TFLOPS against 768 GB/s of memory bandwidth, a ridge point near 202 FLOP per byte. A tensor core kernel needs much larger tiles and deep pipelining to stay compute bound.
• In production, call cuBLAS or build on CUTLASS. Writing these kernels is how you learn what the libraries are doing and how to read a profile, not how you beat them.

The problem and the speed limits

What SGEMM computes

SGEMM is the single precision general matrix multiply from BLAS: C = alpha times A times B plus beta times C, where A is M by K, B is K by N and C is M by N. Throughout this post the matrices are row major and square, with M = N = K = 4096 unless a quoted measurement used a different size (Boehm used 4092).

Every element of C is a dot product of length K, and each step of the dot product is one fused multiply add, which counts as two floating point operations. So the work is:

That last row is the whole story of matrix multiply optimization. The algorithm reuses every input element N times, so a perfect implementation is compute bound by a huge margin. Every inefficiency in a real kernel shows up as moving the same bytes again, which drags the effective intensity down toward the memory bound region.

The two GPUs in the measurements

The published numbers below come from an RTX A6000 (GA102, compute capability 8.6) and, for some comparisons, an A100 (GA100, compute capability 8.0). Their headline specifications, all peak figures at boost clock from NVIDIA's own whitepapers, are:

Sources: the GA102 whitepaper (Table 3, RTX A6000 and A40 specs), the A100 whitepaper (Table 4), and the compute capability tables in the CUDA Programming Guide for shared memory and register limits. Sparse tensor figures (double the dense ones) are left out because a dense GEMM cannot use them.

You can rebuild the FP32 peaks from the table: 84 SMs x 128 lanes x 2 FLOP x 1.8 GHz is 38.7 TFLOPS, and 108 x 64 x 2 x 1.41 GHz is 19.5 TFLOPS. One caveat matters for SGEMM on GA102. The GA102 whitepaper explains that each SM partition has two datapaths and that GA10x added FP32 to the second one, which previous generations had limited to integer work. The second datapath is shared between FP32 and INT32, so a kernel that spends instructions on address arithmetic cannot use all 128 lanes for FP32 on every cycle. The 38.7 TFLOPS figure is a ceiling that assumes every issued instruction is an FP32 FMA.

Rooflines and ridge points

The roofline model divides kernels into two regions. Below the ridge point, performance is limited by memory bandwidth times arithmetic intensity. Above it, performance is limited by peak compute. The ridge point is peak FLOP/s divided by bandwidth:

All figures are peak, not sustained. Two things jump out. First, the A6000 needs four times the FP32 intensity of the A100 before it stops being memory bound, because it has twice the FP32 throughput and half the bandwidth. Second, moving to tensor cores multiplies the required intensity by two to four. The tiling techniques in this post exist to hit those numbers.

There is a second roofline inside the SM. Shared memory has 32 banks, and the CUDA Programming Guide states that each bank delivers 32 bits per clock, so one SM can serve 128 bytes of shared memory per clock when there are no bank conflicts. Across the A6000 that is 84 x 128 bytes x 1.8 GHz, about 19.4 TB/s in aggregate. That is 25 times the DRAM bandwidth, which sounds like plenty until you count how many shared memory loads a naive inner loop issues per FMA.

Lower bounds on time

At 38.7 TFLOPS the A6000 needs 3.54 ms for Boehm's 4092 problem (137.1 GFLOP including the beta term), and at 768 GB/s it needs 0.35 ms to move the 268 MB minimum. Compute dominates by a factor of about ten. Boehm makes the same estimate in his worklog and notes that cuBLAS loads about 500 MB from global memory for the whole calculation, which works out to an effective intensity of about 274 FLOP per byte. That is the kind of reuse the stages below build toward.

The harness: timing, checking and a reference

Before writing a single optimized kernel, set up three things: a timer, a correctness check against a reference, and the reference itself. The reference is cuBLAS. Its API is column major, so the usual trick for row major data is to ask for the transposed product, because the transpose of C equals the transpose of B times the transpose of A, and a row major matrix read as column major is its transpose.

// Host side harness (not executed here; follows the cuBLAS and CUDA runtime docs)
#include <cublas_v2.h>
#include <cuda_runtime.h>

void reference_sgemm(cublasHandle_t h, int M, int N, int K,
                     const float* dA, const float* dB, float* dC) {
  const float alpha = 1.0f, beta = 0.0f;
  // Row-major C = A * B  is column-major C^T = B^T * A^T
  cublasSgemm(h, CUBLAS_OP_N, CUBLAS_OP_N,
              N, M, K, &alpha,
              dB, N,     // B^T viewed column-major, leading dim N
              dA, K,     // A^T viewed column-major, leading dim K
              &beta, dC, N);
}

float time_kernel_ms(void (*launch)(), int reps) {
  cudaEvent_t start, stop;
  cudaEventCreate(&start); cudaEventCreate(&stop);
  launch();                          // warm-up launch
  cudaEventRecord(start);
  for (int i = 0; i < reps; ++i) launch();
  cudaEventRecord(stop);
  cudaEventSynchronize(stop);
  float ms = 0.f;
  cudaEventElapsedTime(&ms, start, stop);
  cudaEventDestroy(start); cudaEventDestroy(stop);
  return ms / reps;
}
// GFLOP/s = 2.0 * M * N * K / (ms * 1e6)

Three habits matter here. Warm up once before timing, so the first launch's module loading and clock ramp do not pollute the result. Average over several launches. And compare every kernel's output with the cuBLAS result using a tolerance that fits FP32 accumulation order differences; an exact comparison will fail even for a correct kernel, because the order of additions changes the rounding.

The kernels below also assume that M, N and K are multiples of the tile sizes. Real code needs bounds checks or padding for the edges. Boehm made the same simplification in his worklog and tested with sizes that satisfy it.

Stage 0: the naive kernel

The simplest correct kernel gives each thread one element of C.

// Stage 0: naive SGEMM. One thread per output element.
// Launch: block (32, 32), grid (M/32, N/32). Not executed here.
__global__ void sgemm_naive(int M, int N, int K, float alpha,
                            const float* A, const float* B,
                            float beta, float* C) {
  const int row = blockIdx.x * blockDim.x + threadIdx.x;  // x walks ROWS
  const int col = blockIdx.y * blockDim.y + threadIdx.y;
  if (row < M && col < N) {
    float acc = 0.0f;
    for (int k = 0; k < K; ++k)
      acc += A[row * K + k] * B[k * N + col];
    C[row * N + col] = alpha * acc + beta * C[row * N + col];
  }
}

This is the shape of Boehm's first kernel, and it is wrong in a way that is invisible to a correctness test. To see why, you have to think about warps rather than threads.

What a warp actually does

A warp is 32 threads that issue the same instruction together. In a 2D block, threadIdx.x varies fastest, so the 32 threads of a warp have consecutive threadIdx.x and the same threadIdx.y. In this kernel that means 32 consecutive rows of C and the same column.

Now trace one iteration of the k loop for one warp:

Warp lanes 0..31 -> rows r, r+1, ..., r+31 ; same column c

Load A[row*K + k]:  addresses  A + r*K + k, A + (r+1)*K + k, ...
                    stride = K floats = 16,384 bytes for K = 4096
                    -> 32 different 32-byte sectors, 4 useful bytes in each

Load B[k*N + col]:  same address for all 32 lanes
                    -> 1 sector, broadcast

Store C[row*N+col]: stride N floats again -> 32 sectors

The CUDA Programming Guide describes exactly this case. Global memory is accessed in 32 byte transactions. If 32 threads read 32 consecutive 4 byte words, the warp needs four transactions and uses all 128 bytes. If consecutive threads read words 32 bytes or more apart, the warp needs one transaction per thread, moves 1,024 bytes, and uses 128 of them: 12.5% efficiency.

Arithmetic intensity and the expected bottleneck

Ignore caches for a moment. Each FMA needs one value of A and one of B, so 8 bytes for 2 FLOP, an intensity of 0.25 FLOP per byte. At 768 GB/s that caps the kernel at 192 GFLOP/s, deep in the memory bound region of a roofline whose ridge sits at 50.4.

The real kernel does a little better because of caching. Boehm measured 309 GFLOP/s for this kernel on the A6000, 1.3% of cuBLAS, and about 0.5 s per 4092 multiply. That is above the no cache bound, which tells you L1 and L2 are catching some of the reuse. But the structure is still dominated by wasted sectors.

Expected Nsight Compute picture: Memory Workload Analysis shows many sectors per request for global loads (the profiling guide's own example of 32 sectors per request means each thread touched a different sector), and Warp State Statistics shows warps stalled mostly on long scoreboard, which the profiling guide defines as waiting on an L1TEX operation such as a global load.

Stage 1: coalesce global memory access

The fix is to change which index each lane walks. Map threadIdx.x to columns instead of rows:

// Stage 1: coalesced. Only the index mapping changes.
// Launch: block (32, 32), grid (N/32, M/32). Not executed here.
__global__ void sgemm_coalesced(int M, int N, int K, float alpha,
                                const float* A, const float* B,
                                float beta, float* C) {
  const int col = blockIdx.x * blockDim.x + threadIdx.x;  // x walks COLUMNS
  const int row = blockIdx.y * blockDim.y + threadIdx.y;
  if (row < M && col < N) {
    float acc = 0.0f;
    for (int k = 0; k < K; ++k)
      acc += A[row * K + k] * B[k * N + col];
    C[row * N + col] = alpha * acc + beta * C[row * N + col];
  }
}

Now a warp covers 32 consecutive columns of one row:

Warp lanes 0..31 -> row r ; columns c, c+1, ..., c+31

Load A[row*K + k]:  same address for all lanes -> 1 sector, broadcast
Load B[k*N + col]:  32 consecutive floats = 128 bytes -> 4 sectors, all used
Store C[...]:       32 consecutive floats -> 4 sectors, all used

Nothing about the arithmetic changed. The intensity in the no reuse model is still 0.25 FLOP per byte, because each lane still loads its own B value for every FMA. What changed is that every byte the memory system moves is a byte someone uses, and the A value is a single broadcast. Boehm measured 1,986.5 GFLOP/s for this kernel, 6.4 times the naive one and 8.5% of cuBLAS.

This is a good moment to be precise about what coalescing is. The CUDA Programming Guide says consecutive threads reading consecutive elements is the typical way to coalesce, not a requirement: what matters is that the lanes of a warp touch as few 32 byte segments as possible, in any order. A permutation within a segment is still coalesced.

Expected bottleneck: still memory, but now L1 and L2 bandwidth rather than wasted DRAM sectors. Many warps on the same SM read the same rows of A and overlapping columns of B, and the caches are doing all the reuse implicitly. The next stage makes that reuse explicit.

Stage 2: shared memory tiling

Shared memory is a programmer managed scratchpad inside each SM, shared by the threads of a block. The idea of tiling is to cooperatively load a square tile of A and a square tile of B into shared memory, synchronize, compute the partial products for that tile from shared memory, synchronize again, and slide the tiles along K.

// Stage 2: shared memory tiling, one output per thread.
// Launch: block (TILE, TILE), grid (N/TILE, M/TILE). Assumes M, N, K % TILE == 0.
// Not executed here.
template <int TILE>
__global__ void sgemm_smem(int M, int N, int K, float alpha,
                           const float* __restrict__ A,
                           const float* __restrict__ B,
                           float beta, float* __restrict__ C) {
  __shared__ float As[TILE][TILE];
  __shared__ float Bs[TILE][TILE];
  const int tx = threadIdx.x, ty = threadIdx.y;
  const int row = blockIdx.y * TILE + ty;
  const int col = blockIdx.x * TILE + tx;

  float acc = 0.0f;
  for (int t = 0; t < K; t += TILE) {
    As[ty][tx] = A[row * K + (t + tx)];        // coalesced: tx walks a row of A
    Bs[ty][tx] = B[(t + ty) * N + col];        // coalesced: tx walks a row of B
    __syncthreads();                           // tile fully loaded
    #pragma unroll
    for (int k = 0; k < TILE; ++k)
      acc += As[ty][k] * Bs[k][tx];
    __syncthreads();                           // everyone done before overwrite
  }
  C[row * N + col] = alpha * acc + beta * C[row * N + col];
}

The arithmetic intensity maths

Count global traffic per block per tile step. The block loads TILE x TILE floats of A and the same of B, 2 x TILE squared x 4 bytes. It then performs TILE cubed FMAs, 2 x TILE cubed FLOP. The ratio is:

global AI = (2 * T^3) / (2 * T^2 * 4) = T / 4   FLOP per byte (FP32)

T = 16  ->  4 FLOP/byte
T = 32  ->  8 FLOP/byte

A tile of 32 cuts global traffic by a factor of 32 compared with the no reuse model, from 0.25 to 8 FLOP per byte. That is still below the A6000's FP32 ridge of 50.4, so DRAM alone would cap this kernel near 8 x 768 GB/s, about 6.1 TFLOPS, if L2 gave no help at all. On the A100, whose ridge is 12.5, the same tile is much closer to compute bound territory, which is one reason tuned parameters differ between the two cards.

The new bottleneck: shared memory bandwidth

Now look at the inner loop from the SM's point of view. For each k, every lane issues two shared memory loads (one from As, one from Bs) and one FMA. With blockDim.x equal to 32, every lane in a warp has the same ty, so As[ty][k] is the same address for the whole warp and the hardware broadcasts it in one pass. Bs[k][tx] is 32 consecutive floats in 32 different banks, also one pass. So each warp level FMA instruction needs two shared memory wavefronts.

A simple throughput model makes the consequence clear. The Nsight Compute profiling guide describes L1TEX as processing one wavefront per cycle. An A6000 SM can issue FMAs for 128 lanes per cycle, which is four warp level FMA instructions. If each FMA instruction needs two wavefronts and the SM serves one wavefront per cycle, FMAs can only issue at half an instruction per cycle:

Model: 1 shared-memory wavefront per SM per clock
       2 wavefronts per warp FMA instruction  ->  0.5 warp FMA per clock
       0.5 x 32 lanes x 2 FLOP              =  32 FLOP per clock per SM
       peak (128 lanes x 2 FLOP)            = 256 FLOP per clock per SM
       ceiling = 32 / 256 = 12.5% of peak   ~ 4.8 TFLOPS on the A6000

That is a model, not a measurement, and it ignores the address arithmetic and loop overhead that the real kernel also pays. Boehm measured 2,980 GFLOP/s for his shared memory kernel, 12.8% of cuBLAS and 7.7% of the A6000's peak, comfortably under the model's ceiling. The measured gain over Stage 1 was only 1.5 times, which fits the model's message. Tiling alone does not make the GPU fast. It trades a DRAM bottleneck for a shared memory one.

Bank conflicts, and when padding helps

Shared memory has 32 banks, and successive 32 bit words map to successive banks. If two lanes of a warp access different words in the same bank, the access is split into multiple wavefronts. If they access the same word, it is a broadcast and costs nothing extra. The Nsight Compute profiling guide says the hardware splits a conflicting request into as many conflict free requests as necessary, reducing effective bandwidth by the degree of the conflict.

The kernel above is conflict free, because one operand is a broadcast and the other is unit stride. Conflicts appear when a warp walks a column of a square tile, because every row of a 32 wide float array starts in bank 0. The CUDA Best Practices Guide shows this with C = A times A transposed: copying a tile of A into shared memory so that it can be read transposed puts every lane of a warp in the same bank, a 32 way conflict. The fix is one extra column:

__shared__ float transposedTile[TILE_DIM][TILE_DIM + 1];  // stride 33 words

With a stride of 33, lane i lands in bank (33 x i) mod 32, which is i, so the column walk becomes conflict free. The guide reports the effect on a Tesla V100 as effective bandwidth:

And for the plain C = A times B case:

Note that these are effective bandwidth figures for the guide's teaching kernels, which compute a w by w tile with w = 32, not GFLOP/s for a full GEMM. Even the best of them is about 22% of the V100's 900 GB/s peak in the A100 whitepaper's comparison table. The guide's lesson, in its own words, is that the gain came "not due to improved coalescing" but from avoiding redundant transfers.

Hopper class hardware adds another tool. The Programming Guide's asynchronous copy chapter describes swizzled shared memory layouts for the tensor memory accelerator, which permute 16 byte chunks across groups of four banks so that both row and column walks avoid conflicts without padding.

Stage 3: 1D register blocking

The way out of the shared memory bottleneck is to make each value loaded from shared memory feed several FMAs. If one thread computes TM outputs in the same column of C, it can load one value of B into a register and reuse it TM times.

// Stage 3: 1D blocktiling. Each thread computes TM outputs in one column.
// BM = BN = 64, BK = 8, TM = 8  ->  512 threads, 4 KB of shared memory.
// Launch: block (BM*BN/TM), grid (N/BN, M/BM). Not executed here.
template <int BM, int BN, int BK, int TM>
__global__ void sgemm_1d(int M, int N, int K, float alpha,
                         const float* __restrict__ A,
                         const float* __restrict__ B,
                         float beta, float* __restrict__ C) {
  __shared__ float As[BM * BK];
  __shared__ float Bs[BK * BN];

  const int threadCol = threadIdx.x % BN;
  const int threadRow = threadIdx.x / BN;            // 0 .. BM/TM - 1

  A += blockIdx.y * BM * K;
  B += blockIdx.x * BN;
  C += blockIdx.y * BM * N + blockIdx.x * BN;

  // With these sizes BM*BK == BK*BN == blockDim.x, one element per thread each.
  const int innerColA = threadIdx.x % BK, innerRowA = threadIdx.x / BK;
  const int innerColB = threadIdx.x % BN, innerRowB = threadIdx.x / BN;

  float acc[TM] = {0.0f};
  for (int bk = 0; bk < K; bk += BK) {
    As[innerRowA * BK + innerColA] = A[innerRowA * K + innerColA];
    Bs[innerRowB * BN + innerColB] = B[innerRowB * N + innerColB];
    __syncthreads();
    A += BK;
    B += BK * N;

    for (int k = 0; k < BK; ++k) {
      const float b = Bs[k * BN + threadCol];        // one load, reused TM times
      #pragma unroll
      for (int i = 0; i < TM; ++i)
        acc[i] += As[(threadRow * TM + i) * BK + k] * b;
    }
    __syncthreads();
  }
  for (int i = 0; i < TM; ++i) {
    float* c = &C[(threadRow * TM + i) * N + threadCol];
    *c = alpha * acc[i] + beta * (*c);
  }
}

These are the parameters Boehm used for his fourth kernel (BM = BN = 64, BK = 8, TM = 8), with a shared memory footprint of 64 x 8 + 64 x 8 = 1,024 floats, 4 KB per block.

The maths

Global intensity. Per tile step the block loads (BM + BN) x BK floats and does BM x BN x BK FMAs:

global AI = (2 * BM * BN * BK) / ((BM + BN) * BK * 4) = BM*BN / (2*(BM + BN))

BM = BN = 64   ->  16 FLOP/byte

Shared memory intensity. Within the inner loop each lane issues 1 load of B plus TM loads of A for TM FMAs. Because a warp's 32 lanes share threadRow (the warp sits inside one 64 wide row of threads), the A loads are broadcasts and the B loads are unit stride. Using the same one wavefront per cycle model:

FMAs per shared-memory load = TM / (1 + TM) = 8 / 9
ceiling = (8/9) x 32 lanes x 2 FLOP = 56.9 FLOP per clock per SM
        = 56.9 / 256 = 22% of peak  ~ 8.6 TFLOPS on the A6000

Boehm measured 8,474.7 GFLOP/s for this kernel, 36.5% of cuBLAS and 2.8 times Stage 2. That sits just under the model's ceiling, which suggests shared memory bandwidth is still what limits it. The fix is to reuse both operands, not just B.

Expected Nsight Compute picture: warp stalls move from long scoreboard toward MIO throttle and short scoreboard. The profiling guide describes MIO throttle as warps waiting for the memory input/output instruction queue, which fills under heavy shared memory traffic, and short scoreboard as waiting on a non L1TEX memory operation, typically shared memory.

Stage 4: 2D register blocking

If each thread computes a TM by TN block of outputs, it can load TM values of A and TN values of B into registers and perform TM x TN FMAs from them. This is the outer product formulation, and it is the heart of every fast GEMM.

// Stage 4: 2D blocktiling. Each thread computes a TM x TN block of C.
// BM = BN = 128, BK = 8, TM = TN = 8 -> 256 threads, 8 KB shared memory.
// Launch: block (BM*BN/(TM*TN)), grid (N/BN, M/BM). Not executed here.
template <int BM, int BN, int BK, int TM, int TN>
__global__ void __launch_bounds__((BM * BN) / (TM * TN))
sgemm_2d(int M, int N, int K, float alpha,
         const float* __restrict__ A, const float* __restrict__ B,
         float beta, float* __restrict__ C) {
  constexpr int NT = (BM * BN) / (TM * TN);          // threads per block
  __shared__ float As[BM * BK];
  __shared__ float Bs[BK * BN];

  const int threadCol = threadIdx.x % (BN / TN);
  const int threadRow = threadIdx.x / (BN / TN);

  A += blockIdx.y * BM * K;
  B += blockIdx.x * BN;
  C += blockIdx.y * BM * N + blockIdx.x * BN;

  const int innerRowA = threadIdx.x / BK, innerColA = threadIdx.x % BK;
  constexpr int strideA = NT / BK;                   // rows of A per pass
  const int innerRowB = threadIdx.x / BN, innerColB = threadIdx.x % BN;
  constexpr int strideB = NT / BN;                   // rows of B per pass

  float acc[TM * TN] = {0.0f};
  float regM[TM], regN[TN];

  for (int bk = 0; bk < K; bk += BK) {
    for (int off = 0; off < BM; off += strideA)
      As[(innerRowA + off) * BK + innerColA] = A[(innerRowA + off) * K + innerColA];
    for (int off = 0; off < BK; off += strideB)
      Bs[(innerRowB + off) * BN + innerColB] = B[(innerRowB + off) * N + innerColB];
    __syncthreads();
    A += BK;
    B += BK * N;

    #pragma unroll
    for (int k = 0; k < BK; ++k) {
      #pragma unroll
      for (int i = 0; i < TM; ++i) regM[i] = As[(threadRow * TM + i) * BK + k];
      #pragma unroll
      for (int j = 0; j < TN; ++j) regN[j] = Bs[k * BN + threadCol * TN + j];
      #pragma unroll
      for (int i = 0; i < TM; ++i)
        #pragma unroll
        for (int j = 0; j < TN; ++j)
          acc[i * TN + j] += regM[i] * regN[j];      // TM*TN FMAs from TM+TN loads
    }
    __syncthreads();
  }

  #pragma unroll
  for (int i = 0; i < TM; ++i)
    #pragma unroll
    for (int j = 0; j < TN; ++j) {
      float* c = &C[(threadRow * TM + i) * N + threadCol * TN + j];
      *c = alpha * acc[i * TN + j] + beta * (*c);
    }
}

With NT = 256, strideA is 32, so each thread loads 4 values of A per tile step, and strideB is 2, so each thread loads 4 values of B. These are the parameters of Boehm's fifth kernel.

The maths

Per thread, per k:   TM + TN = 16 shared loads,  TM * TN = 64 FMAs
FMAs per shared load = 64 / 16 = 4

global AI = BM*BN / (2*(BM + BN)) = 128*128 / (2*256) = 32 FLOP/byte

Four FMAs per shared memory load means the simple wavefront model no longer caps the kernel below peak: 4 x 32 lanes x 2 FLOP is 256 FLOP per clock, exactly the SM's FP32 peak. The global intensity of 32 FLOP per byte is still under the A6000 ridge of 50.4, but L2 now catches a large share of the remaining reuse, because neighbouring blocks read the same rows of A and columns of B.

Boehm measured 15,971.7 GFLOP/s, 68.7% of cuBLAS, 1.9 times Stage 3. The kernel is now in the territory where instruction mix, register pressure and latency hiding matter more than raw bandwidth.

Register pressure and occupancy, worked

The accumulators alone are TM x TN = 64 registers per thread, plus 8 for regM, 8 for regN, and addresses and loop counters: comfortably over 80, and often over 100 after unrolling. You can see the real number at compile time with nvcc's -Xptxas -v (pass verbose to the assembler) or --resource-usage, both documented in the nvcc reference.

Occupancy is the ratio of resident warps to the maximum an SM supports. The CUDA Programming Guide's compute capability tables give 65,536 registers per SM and a 255 register per thread limit on all current parts, and a maximum of 1,536 resident threads per SM on compute capability 8.6 versus 2,048 on 8.0. Working through a few cases for a 256 thread block, ignoring allocation granularity:

Two lessons come out of this table. First, there is a cliff: going from 64 to 65 registers per thread halves the number of resident blocks. Second, low occupancy is normal for GEMM. The CUTLASS documentation says it plainly: accumulators typically take at least half of a thread's register budget, so occupancy is relatively low, and latency is hidden with software pipelining rather than with more warps. Boehm makes a related point, citing Vasily Volkov's thesis: at very high arithmetic intensity, high occupancy is not needed for peak throughput, because each warp has plenty of independent FMAs to issue while its loads are in flight.

Shared memory is not the limiter here. At 8 KB per block, the A6000's 100 KB per SM would fit 12 blocks, far more than registers allow.

Stage 5: vectorized loads

The 2D kernel issues 32 bit loads from shared memory for A because As is laid out with stride BK between the TM values a thread needs. Boehm's SASS listing in his 1D blocktiling section already shows the pattern: the compiler turned the contiguous Bs reads into 128 bit LDS.128 instructions, but the As reads stayed as 32 bit LDS loads, and he points out that this hints at transposing As. Stage 5 fixes that in two ways:

1. Transpose A while copying it into shared memory, so that the TM values a thread needs for one k are contiguous.
2. Use float4 for global loads, so each thread moves 16 bytes per instruction instead of 4.

// Stage 5 loader: float4 global loads, A transposed into shared memory.
// BM = BN = 128, BK = 8, 256 threads: one float4 of A and one of B per thread.
// Requires K and N to be multiples of 4 and 16-byte-aligned base pointers.
const int innerRowA = threadIdx.x / (BK / 4);        // 0..127
const int innerColA = threadIdx.x % (BK / 4);        // 0..1
const int innerRowB = threadIdx.x / (BN / 4);        // 0..7
const int innerColB = threadIdx.x % (BN / 4);        // 0..31

const float4 a4 = reinterpret_cast<const float4*>(&A[innerRowA * K + innerColA * 4])[0];
As[(innerColA * 4 + 0) * BM + innerRowA] = a4.x;     // As is now [BK][BM]
As[(innerColA * 4 + 1) * BM + innerRowA] = a4.y;
As[(innerColA * 4 + 2) * BM + innerRowA] = a4.z;
As[(innerColA * 4 + 3) * BM + innerRowA] = a4.w;

reinterpret_cast<float4*>(&Bs[innerRowB * BN + innerColB * 4])[0] =
    reinterpret_cast<const float4*>(&B[innerRowB * N + innerColB * 4])[0];
__syncthreads();

// Inner loop now reads As contiguously:
//   regM[i] = As[k * BM + threadRow * TM + i];   -> compiler can emit LDS.128

The alignment requirement is real. A float4 access must be 16 byte aligned. The Best Practices Guide states that cudaMalloc returns memory aligned to at least 256 bytes, so the base pointers are fine, but row offsets are only aligned if K and N are multiples of 4. Treat a misaligned vector access as a bug to prevent, not a slow path to tolerate.

What changes in the numbers. Global traffic in bytes is the same as Stage 4; the number of load instructions drops by four. Shared memory traffic in bytes is also the same; the number of shared memory instructions for As drops by up to four, and each LDS.128 moves 16 bytes per lane. The win is in instruction issue and in the MIO queue, not in bandwidth. Boehm measured 18,237.3 GFLOP/s, 78.4% of cuBLAS, 1.14 times Stage 4. He then autotuned the tile sizes (his kernel 9) and found that BM = BN = 128, BK = 16, TM = TN = 8 gained about 5% on the A6000 (his results table lists 19,721 GFLOP/s for the autotuned kernel), while the best setting on an A100 was different (BM = BN = 64, BK = 16, TM = TN = 4). Tile sizes are hardware specific, which is why libraries ship many kernels.

The transposed store into As has its own trap: four lanes writing four rows of a column can conflict. Boehm reports that his profile still showed shared memory bank conflicts at this stage, which cuBLAS avoids, and that his attempts to remove them with layout changes (his kernels 7 and 8) removed the conflicts but were slower overall. Swizzled layouts, where the column index is XORed with part of the row index, are the usual production answer.

Stage 6: double buffering with asynchronous copies

Every kernel so far has the same rhythm: load a tile, __syncthreads(), compute, __syncthreads(), repeat. While the tile loads, the FMA units idle; while the FMAs run, the memory system idles, except for whatever other resident blocks happen to be doing. With low occupancy, there are few other blocks to fill the gaps.

Double buffering allocates two shared memory tiles. While the block computes on one, it loads the next tile into the other. The CUTLASS documentation describes this at two scopes: two shared memory tiles per threadblock, and two register fragments per warp, so loads from shared memory into registers also overlap with math.

On compute capability 8.0 and later, there is hardware support for the first scope. The A100 whitepaper describes a new asynchronous copy instruction that loads from global memory directly into shared memory, bypassing the register file and optionally L1. In PTX it is cp.async, which the PTX ISA lists as requiring sm_80 or higher. In CUDA C++ the Programming Guide documents three ways to reach it: cuda::memcpy_async, cooperative_groups::memcpy_async, and the low level primitives __pipeline_memcpy_async, __pipeline_commit and __pipeline_wait_prior from cuda_pipeline.h. The primitives are the easiest to read in a GEMM main loop:

// Stage 6: double buffering with cp.async (compute capability 8.0+).
// BM = BN = 128, BK = 8, TM = TN = 8, 256 threads. A stays row-major in smem here,
// because cp.async copies bytes as-is and cannot transpose. Not executed here.
#include <cuda_pipeline.h>

template <int BM, int BN, int BK, int TM, int TN>
__global__ void __launch_bounds__((BM * BN) / (TM * TN))
sgemm_2d_db(int M, int N, int K, float alpha,
            const float* __restrict__ A, const float* __restrict__ B,
            float beta, float* __restrict__ C) {
  __shared__ __align__(16) float As[2][BM * BK];
  __shared__ __align__(16) float Bs[2][BK * BN];

  A += blockIdx.y * BM * K;
  B += blockIdx.x * BN;
  C += blockIdx.y * BM * N + blockIdx.x * BN;

  const int threadCol = threadIdx.x % (BN / TN);
  const int threadRow = threadIdx.x / (BN / TN);
  // One 16-byte copy of A and one of B per thread per stage:
  // BM*BK/4 == BK*BN/4 == 256 == blockDim.x for these sizes.
  const int rA = threadIdx.x / (BK / 4), cA = (threadIdx.x % (BK / 4)) * 4;
  const int rB = threadIdx.x / (BN / 4), cB = (threadIdx.x % (BN / 4)) * 4;

  auto load_stage = [&](int s, int k0) {
    __pipeline_memcpy_async(&As[s][rA * BK + cA], &A[rA * K + k0 + cA], 16);
    __pipeline_memcpy_async(&Bs[s][rB * BN + cB], &B[(k0 + rB) * N + cB], 16);
    __pipeline_commit();
  };

  float acc[TM * TN] = {0.0f};
  float regM[TM], regN[TN];

  load_stage(0, 0);                                  // prologue: fill stage 0
  int s = 0;
  for (int k0 = 0; k0 < K; k0 += BK) {
    if (k0 + BK < K) {
      load_stage(s ^ 1, k0 + BK);                    // start the NEXT tile
      __pipeline_wait_prior(1);                      // current tile has landed
    } else {
      __pipeline_wait_prior(0);                      // last tile: drain everything
    }
    __syncthreads();                                 // all threads' copies visible

    #pragma unroll
    for (int k = 0; k < BK; ++k) {
      #pragma unroll
      for (int i = 0; i < TM; ++i) regM[i] = As[s][(threadRow * TM + i) * BK + k];
      #pragma unroll
      for (int j = 0; j < TN; ++j) regN[j] = Bs[s][k * BN + threadCol * TN + j];
      #pragma unroll
      for (int i = 0; i < TM; ++i)
        #pragma unroll
        for (int j = 0; j < TN; ++j)
          acc[i * TN + j] += regM[i] * regN[j];
    }
    __syncthreads();                                 // done reading stage s
    s ^= 1;
  }
  // epilogue as in Stage 4
}

The ordering is the part people get wrong. At the top of iteration k0, the kernel starts filling stage s ^ 1. That buffer was last read in the previous iteration, which ended with __syncthreads(), so no thread can still be reading it. __pipeline_wait_prior(1) waits until at most one committed batch is still in flight for this thread, which means the batch for stage s is complete. The following __syncthreads() makes every thread's completed copies visible to every other thread.

The cost. Shared memory per block doubles, from 8 KB to 16 KB here, and with BK = 16 it would be 32 KB. That is still small next to the 100 KB available per SM on the A6000 or 164 KB on the A100, but tensor core kernels with larger tiles and three or four stages run into the limit quickly, which is why the A100 raised the maximum shared memory per SM to 164 KB and Hopper (compute capability 9.0) raised it to 228 KB, per the Programming Guide's tables.

Published measurements. Boehm's worklog lists double buffering as not implemented; his best kernel did not use it. We do not have a published, like for like measurement of this exact stage on the A6000, and we will not invent one. CUTLASS uses this technique in its pipelined main loop, and that is the strongest primary evidence that it matters.

Stage 7: warp tiling

The 2D kernel treats the 256 threads of a block as a flat grid of 16 by 16 threads, each owning an 8 by 8 output block. That ignores the fact that warps are the unit of scheduling and that shared memory bank behaviour is decided per warp. Warp tiling adds an explicit level to the hierarchy: the block tile is split into warp tiles, and each warp tile is split into per thread tiles.

Block tile  128 x 128  (8 warps, 256 threads)
+------------------+------------------+------------------+------------------+
| warp 0  64 x 32  | warp 1  64 x 32  | warp 2  64 x 32  | warp 3  64 x 32  |
|                  |                  |                  |                  |
+------------------+------------------+------------------+------------------+
| warp 4  64 x 32  | warp 5  64 x 32  | warp 6  64 x 32  | warp 7  64 x 32  |
|                  |                  |                  |                  |
+------------------+------------------+------------------+------------------+

Warp tile 64 x 32 = 2,048 outputs / 32 lanes = 64 outputs per lane (8 x 8)
Each lane's 64 outputs are split into sub-tiles so that, for each k, the
32 lanes read contiguous float4 runs of As and Bs with no bank conflicts.

// Stage 7 index mapping (sketch). The inner product loop is the Stage 5 loop,
// iterated over warp sub-tiles. Not executed here.
constexpr int WM = 64, WN = 32;                      // warp tile
const int warpId  = threadIdx.x / 32;
const int lane    = threadIdx.x % 32;
const int warpRow = warpId / (BN / WN);              // 0..1
const int warpCol = warpId % (BN / WN);              // 0..3
// Base of this warp's tile inside the block tile:
//   rows  warpRow * WM .. warpRow * WM + WM - 1
//   cols  warpCol * WN .. warpCol * WN + WN - 1
// Lane layout inside the warp tile is chosen so that each k-step of
// regM/regN loads is a broadcast or a contiguous, conflict-free LDS.128.

The benefits are indirect but real. Each warp's working set of As and Bs is smaller and more localized, which improves register cache behaviour. Lane layouts can be chosen to make every shared memory read a broadcast or a contiguous vector read. And the structure maps one to one onto how tensor cores work, because tensor core instructions are warp level: the warp tile is exactly where tensor core fragments live.

Boehm's warptiled kernel reached 21,779.3 GFLOP/s on the A6000, 93.7% of cuBLAS's 23,249.6 GFLOP/s. On an A100 he measured his warptiling kernel within a few percent of cuBLAS at sizes 2048 and 4096, and slower at small sizes, where cuBLAS switches to different kernels including split K.

Stage 8: tensor cores through WMMA

Everything up to now uses the FP32 CUDA cores. Tensor cores are separate matrix units that perform a small matrix multiply accumulate per instruction. The CUDA Programming Guide exposes them in C++ as warp matrix functions in the nvcuda::wmma namespace, available from compute capability 7.0. The basic contract:

• A fragment holds a piece of a matrix distributed across the 32 lanes of a warp. The mapping of elements to lanes is unspecified and can change between architectures.
• load_matrix_sync, mma_sync and store_matrix_sync are warp collective. All 32 lanes must call them, and calling them under a condition that differs between lanes can hang.
• The pointer passed to load and store must be 256 bit aligned, and the leading dimension must be a multiple of 8 for __half elements and of 4 for float.
• For FP16 inputs the common shape is m16n16k16 with FP32 accumulators. For TF32, the guide says the only supported shape is m16n16k8, and the inputs must be converted with __float_to_tf32 first.

A minimal, correct but slow WMMA kernel gives each warp one 16 by 16 tile of C and walks K:

// Stage 8a: one warp per 16x16 output tile, straight from global memory.
// A: M x K row-major half, B: K x N row-major half, C: M x N row-major float.
// M, N, K multiples of 16. Launch: 128 threads per block (4 warps).
// Not executed here.
#include <mma.h>
using namespace nvcuda;

__global__ void hgemm_wmma_naive(int M, int N, int K,
                                 const half* A, const half* B, float* C) {
  const int warpId  = (blockIdx.x * blockDim.x + threadIdx.x) / 32;
  const int tilesN  = N / 16;
  const int tileRow = warpId / tilesN;
  const int tileCol = warpId % tilesN;
  if (tileRow * 16 >= M) return;                     // warp-uniform condition

  wmma::fragment<wmma::matrix_a, 16, 16, 16, half, wmma::row_major> a;
  wmma::fragment<wmma::matrix_b, 16, 16, 16, half, wmma::row_major> b;
  wmma::fragment<wmma::accumulator, 16, 16, 16, float> acc;
  wmma::fill_fragment(acc, 0.0f);

  for (int k = 0; k < K; k += 16) {
    wmma::load_matrix_sync(a, A + tileRow * 16 * K + k, K);
    wmma::load_matrix_sync(b, B + k * N + tileCol * 16, N);
    wmma::mma_sync(acc, a, b, acc);                  // 16*16*16 = 4,096 FMAs
  }
  wmma::store_matrix_sync(C + tileRow * 16 * N + tileCol * 16, acc, N,
                          wmma::mem_row_major);
}

The early return is safe only because every lane of a warp computes the same warpId and so takes the same branch. Make the condition depend on the lane and the warp collective calls below it become undefined.

The maths, and why this kernel is slow

One mma_sync with m16n16k16 is 4,096 FMAs, 8,192 FLOP. To feed it, this kernel loads a 16 by 16 tile of A and one of B, 512 half values, 1,024 bytes, from global memory or L2. That is 8 FLOP per byte, the same intensity as the FP32 shared memory kernel in Stage 2, against a ridge point of about 202 FLOP per byte for FP16 tensor cores on the A6000. The tensor cores spend nearly all their time waiting.

The cure is everything from Stages 2 to 7, applied again at a larger scale:

Tensor core GEMM structure (FP16 in, FP32 accumulate)

Block tile    BM x BN, e.g. 128 x 128, staged through shared memory
              with 2 to 4 cp.async stages

Warp tile     e.g. 64 x 32 per warp = 4 x 2 WMMA 16x16 fragments
              -> each A fragment reused 2 times, each B fragment 4 times

Global AI per block tile (FP16 inputs, 2 bytes per element):
  2*BM*BN*BK / ((BM + BN) * BK * 2) = BM*BN / (BM + BN)
  128 x 128 ->  64 FLOP/byte
  256 x 128 ->  85.3 FLOP/byte

Even a 256 by 128 block tile is under the 202 FLOP per byte ridge if every block fetched its tiles from DRAM. It works in practice because neighbouring blocks share tiles through L2, which is much faster than DRAM. CUTLASS's documentation describes "threadblock rasterization", remapping block indices so that consecutively launched blocks cover a compact 2D region of C and hit the same A and B tiles in L2 at about the same time. That is a tensor core era optimization: on the FP32 path, the ridge was low enough that it mattered less.

The minimum possible intensity for the whole FP16 problem at size 4096 is N/3, about 1,365 FLOP per byte, so the problem itself is never the limit. The kernel's reuse is.

TF32 and precision

If your data is FP32 and you want tensor cores, the choices are converting to FP16 or BF16, or using TF32. The Programming Guide describes TF32 as having the range of FP32 with at least 10 bits of precision, and cuBLAS exposes it through the compute type CUBLAS_COMPUTE_32F_FAST_TF32 or the math mode CUBLAS_TF32_TENSOR_OP_MATH, both documented as rounding inputs to nearest even. On the A6000, the dense TF32 tensor peak is 77.4 TFLOPS, twice the FP32 CUDA core peak. Boehm notes that cuBLAS with TF32 or BF16 increased FLOP/s by 2.5 or 3.5 times in his setting. Whether TF32 is acceptable is a numerical decision for your workload, not a performance one: it changes results.

Beyond WMMA: mma.sync, ldmatrix, wgmma and tcgen05

WMMA is the portable C++ layer. Below it, the PTX ISA exposes the instructions that libraries actually use:

The direction of travel is clear from the table. Each generation moves more of the data movement and the matrix math into asynchronous, larger granularity operations, so that a few warps can keep the tensor cores busy while other warps only move data. CUTLASS describes this for Hopper as warp specialization: a producer warp group issues TMA copies into shared memory buffers, and one or two consumer warp groups run wgmma on them, with barriers signalling when buffers are full or empty. Its "ping pong" design gives each of two consumer warp groups its own output tile so one group's epilogue overlaps with the other's math. The CUTLASS docs note that Blackwell SM100 kernels use a substantially different warp specialization structure, though the producer and consumer idea carries over.

Hand writing kernels at this level is a project measured in months. The rest of the industry, including NVIDIA's own libraries, uses templates and generators for it.

cuBLAS and CUTLASS as the reference

cuBLAS

cuBLAS is the reference every kernel in this post is measured against. Two calls cover most cases:

// FP32 on CUDA cores (or tensor cores if the math mode allows it)
cublasSgemm(h, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K,
            &alpha, dB, N, dA, K, &beta, dC, N);

// FP32 data, TF32 tensor core compute, explicitly requested
cublasGemmEx(h, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K,
             &alpha, dB, CUDA_R_32F, N,
                     dA, CUDA_R_32F, K,
             &beta,  dC, CUDA_R_32F, N,
             CUBLAS_COMPUTE_32F_FAST_TF32, CUBLAS_GEMM_DEFAULT);

The cuBLAS documentation describes CUBLAS_DEFAULT_MATH as the default, highest performance mode that uses compute precision at least as high as requested, and CUBLAS_PEDANTIC_MATH as the mode intended for numerical studies, testing and debugging. If you are comparing your FP32 kernel against cuBLAS, make sure cuBLAS is not quietly using a reduced precision path you are not, or the comparison is meaningless in both directions.

A useful fact from Boehm's worklog: cuBLAS is not one kernel. He traced the calls for square sizes up to 4096 and found 16 different SGEMM kernels selected by problem size, including at size 256 a split K kernel followed by a reduction kernel. Split K divides the K dimension across blocks so that small output matrices still fill the GPU, at the cost of a second pass to sum the partial results. CUTLASS documents the same technique under parallelized reductions.

CUTLASS

CUTLASS is NVIDIA's open source C++ template library for GEMM and related operations. NVIDIA described its hierarchical structure in a technical blog post dated 5 December 2017, and the project's "Efficient GEMM in CUDA" document is the closest thing to an official explanation of everything in this post. Its loop nest is the whole ladder in one picture:

// Adapted from CUTLASS docs, "Efficient GEMM in CUDA" (structure only)
for cta_n, cta_m:                 // threadblock tiles of C
  for cta_k:                      // GEMM main loop, one iteration = one stage
    for warp_n, warp_m:           // warp tiles inside the block tile
      for warp_k:                 // fully unrolled across the block's K tile
        for mma_k, mma_n, mma_m:  // tensor core (or FMA) instructions
          mma_instruction(d, a, b, c);

Every stage above appears there. Block tiles are Stage 2. Per thread or per warp register accumulators are Stages 3, 4 and 7. Software pipelining with two shared memory tiles and two register fragments is Stage 6. The epilogue, where alpha, beta and any fused operation such as a bias or activation is applied while the results are still in registers, is a stage we only touched in passing, and it is where CUTLASS earns much of its value in real models.

The published results in one table

Here are all the measured numbers quoted in this post, with what was measured and where. We ran none of them.

The GFLOP/s and share of cuBLAS columns are Boehm's. The share of peak and time columns are our arithmetic on his figures (137.1 GFLOP per GEMM), checked in Python. Note that even cuBLAS sits at about 60% of the whitepaper FP32 peak on this card; part of that gap is the shared FP32 and INT32 datapath described earlier. On an A100, Boehm reports cuBLAS at 14.7 TFLOPS (75% of the 19.5 TFLOPS peak) and his autotuned kernel at 12.6 TFLOPS.

The stage to stage ratios tell you where effort pays. In Boehm's sequence the largest single jumps are coalescing (6.4 times), 1D register blocking (2.8 times) and 2D register blocking (1.9 times). Shared memory tiling on its own was only 1.5 times, and the last three stages together added about 36%. The total from naive to warptiled is about 70 times.

Profiling each stage with Nsight Compute

Nsight Compute is NVIDIA's kernel profiler. The command line tool is ncu. Build with line information so the profiler can map metrics back to source lines:

# Build (not executed here; flags from the nvcc reference)
nvcc -O3 -arch=sm_86 -lineinfo -Xptxas -v -o sgemm sgemm.cu -lcublas

# Full report for every kernel whose name matches, saved to a file
ncu --set full -k regex:sgemm -o sgemm_full ./sgemm

# Just the sections that answer "what is the bottleneck?"
ncu --section SpeedOfLight \
    --section MemoryWorkloadAnalysis_Tables \
    --section WarpStateStats \
    --section Occupancy --section LaunchStats \
    -k regex:sgemm_2d ./sgemm

# A few headline metrics across all kernels
ncu --metrics gpu__time_duration.sum,sm__throughput.avg.pct_of_peak_sustained_elapsed,launch__registers_per_thread,launch__occupancy_limit_registers \
    -k regex:sgemm ./sgemm

# Re-open a saved report on another machine
ncu --import sgemm_full.ncu-rep --page details

The Nsight Compute CLI documentation says that without --set, the basic set is collected, and that --set full collects every section. Full collection replays each kernel many times, so expect a large slowdown while profiling; time your kernels with events, not under the profiler.

What to look at, stage by stage

The stall reasons are defined in the profiling guide. Long scoreboard is a wait on an L1TEX operation (local, global, surface or texture memory). Short scoreboard is a wait on a memory input/output operation that is not L1TEX, most often shared memory. MIO throttle is a full MIO instruction queue, typically from heavy shared memory use. Barrier is a wait at __syncthreads() for sibling warps. Math pipe throttle is all eligible warps wanting the same, oversubscribed math pipe, which in a well tuned GEMM is the stall you want to see, because it means you are compute bound.

The Speed Of Light section includes hierarchical roofline charts for single, double and half precision. Put each of your stages on the single precision chart and you will see the story of this post as a path from bottom left (low intensity, low throughput) toward the ridge.

Common mistakes

• Mapping threadIdx.x to rows of a row major matrix. It passes every correctness test and costs you a factor of six or more. Check which index moves fastest within a warp before anything else.
• Comparing against cuBLAS in a different precision. If cuBLAS is allowed to use TF32 tensor cores and your kernel uses FP32 FMAs, a large gap is expected and tells you nothing about your kernel. Pin the compute type or math mode.
• Forgetting the second __syncthreads(). Without the barrier after the compute loop, fast warps overwrite a shared tile while slow warps are still reading it. The result is wrong only sometimes, which is worse than always.
• Exact equality checks. FP32 sums in a different order produce different rounding. Use a relative tolerance and report the maximum error.
• Misaligned vector loads. float4 and 16 byte cp.async need 16 byte alignment, which fails as soon as K or N is not a multiple of 4. Add a scalar fallback or pad the matrices.
• Chasing occupancy. Raising occupancy by shrinking per thread tiles lowers arithmetic intensity, which is the thing actually limiting a GEMM. Low occupancy with large register tiles is the normal, documented design point.
• Warp divergent WMMA calls. load_matrix_sync, mma_sync and store_matrix_sync must be called by all 32 lanes under the same condition. The Programming Guide warns that otherwise the code is likely to hang.
• Passing WMMA fragments across separately compiled code. The Programming Guide warns that fragment layout is architecture specific and passing fragments between functions compiled for different architectures can silently corrupt results. Store to memory at interfaces.
• Timing under the profiler. ncu replays kernels to collect metrics; timings it reports for full sets are not your production timings.
• Hard coding one tile size. Boehm's best configuration on the A6000 was not the best on the A100. Autotune per target, or let a library choose.

FAQ

Why is my naive CUDA matrix multiply so slow?

Usually because consecutive threads in a warp read memory with a large stride, so each 32 byte transaction delivers 4 useful bytes. Even when coalesced, a kernel that loads two values per FMA has an intensity of 0.25 FLOP per byte, far below the 50 FLOP per byte the A6000 needs on its FP32 path. You need tiling and register blocking to reuse data.

What is memory coalescing in CUDA?

It is the warp combining its 32 lanes' memory requests into as few 32 byte transactions as possible. The CUDA Programming Guide shows the best case (consecutive 4 byte words, 4 transactions for 128 bytes) and the worst case (stride of 32 bytes or more, 32 transactions, 12.5% of fetched bytes used).

How much does shared memory tiling speed up a matrix multiply?

On its own, less than people expect. In Boehm's A6000 measurements it went from 1,986.5 to 2,980.3 GFLOP/s, 1.5 times, because the bottleneck moved to shared memory bandwidth. It pays off fully once each thread computes several outputs from registers.

What is register blocking or thread coarsening in GEMM?

Each thread computes a small block of outputs, for example 8 by 8, holding the accumulators in registers. Per step it loads TM values of A and TN of B and performs TM x TN FMAs, so the ratio of math to shared memory loads rises from 1 to 2 up to 64 to 16 in our example.

What occupancy should a GEMM kernel target?

Lower than most kernels. CUTLASS's documentation notes that accumulators typically take at least half the register budget, so occupancy is relatively low, and latency is hidden with software pipelining. On compute capability 8.6, a 256 thread block at 128 registers per thread gives two blocks per SM, 33% occupancy, and that can be fine.

How do I use tensor cores from CUDA C++?

The simplest route is the WMMA API in mma.h, with warp wide fragments and load_matrix_sync, mma_sync and store_matrix_sync. Production kernels use PTX mma.sync and ldmatrix, or on Hopper wgmma with TMA, usually through CUTLASS rather than by hand.

Is TF32 accurate enough for FP32 matrix multiplication?

That depends on your workload. TF32 keeps FP32's range but has about 10 bits of mantissa precision, per the Programming Guide, and cuBLAS rounds inputs to nearest even when using it. Many training workloads accept it; numerical codes that need full FP32 precision should not enable it. Test against an FP32 reference.

Can a hand written kernel beat cuBLAS?

For a specific shape on a specific GPU, sometimes; in general, rarely. Boehm's best FP32 kernel reached 93.7% of cuBLAS at size 4092 on an A6000, and he found cuBLAS chooses among 16 kernels by size. Write your own when you need fusion that the library does not offer, and build on CUTLASS if you can.

Which Nsight Compute sections matter most for GEMM?

Speed Of Light for the overall bound and roofline position, Memory Workload Analysis for sectors per request and shared memory bank conflicts, Warp State Statistics for stall reasons, and Occupancy with Launch Statistics for register limits. Use ncu --set full once, then targeted sections while iterating.

Where Safeguard fits

Safeguard does not design GPUs, does not write or tune kernels, and does not make any of the code in this post faster. Nothing here changes because Safeguard exists.

Where we do fit is the software that surrounds these kernels in production. A GPU workload in a container ships a CUDA toolkit or runtime, cuBLAS, often CUTLASS headers compiled into your binary, cuDNN, NCCL, a Python framework, a model serving stack and dozens of open source dependencies, all on top of a host driver. Those are versioned components that receive security advisories like any other software. Safeguard's GA capabilities apply to them: SBOM generation and analysis in SPDX and CycloneDX, dependency scanning of repositories, container images and SBOM files from the CLI (safeguard scan --dir ., --image <name>, --sbom <file>), matching against vulnerability databases, reachability analysis, and policies and gates. What we do not cover is the host kernel driver installed outside your image, or GPU firmware; those updates come from the vendor and your operating system.

A tuned GEMM kernel is only one file in a GPU container. The libraries, runtimes and Python packages around it change every month, and a rebuilt image pulls whatever versions its lockfiles and base image allow, so keeping those dependencies patched is still ordinary supply chain work.

More in this series

Other posts in our GPUs, CPUs, QPUs and lithography series:

• AMD GPU Architecture Explained: GCN, RDNA, CDNA, Infinity Fabric and HIP
• CPU vs GPU vs QPU Explained: How Each Works and Why They Coexist
• CUDA Programming Model Explained: Grids, Blocks, Warps and Memory
• GPU Memory and Interconnects: HBM, GDDR, NVLink, PCIe, CXL
• CUDA vs HIP vs SYCL vs OpenCL vs Metal vs Vulkan vs Triton Compared
• GPU Software Supply Chain: Drivers, CUDA, cuDNN, NCCL and Your SBOM
• How a QPU Works: Qubit Chips, Fridges, Control Pulses and Readout
• How CPUs Work: Pipelines, Out-of-Order, Caches, Branches, SIMD
• How EUV Lithography Machines Work: Tin Plasma, Mirrors and High-NA
• How GPUs Work: Graphics Pipeline, Shader Cores, Ray Tracing and Compute
• Inside an NVIDIA SM: Warp Schedulers, Registers, Tensor Cores
• Intel GPU Architecture: Xe Cores, XMX, Ponte Vecchio and oneAPI SYCL
• NVIDIA GPU Architecture History: From Tesla and Fermi to Blackwell
• NVIDIA vs AMD vs Intel vs Apple GPUs: Architectures Compared
• Photolithography Explained: Contact Printing to DUV Immersion and Overlay
• Tile-Based GPUs Explained: Apple, Mali, Adreno, PowerVR and TBDR

Fix it without waiting for upstream

Whatever the vulnerability, most of the risk sits in two gaps. The first is the time between a fix being published and that fix actually running in your production builds. The second, and worse, is the time before any fix exists at all, while a maintainer works out a patch and a release schedule. Waiting through both gaps is a habit, not a requirement. Safeguard gives teams three ways to close them.

1. Autonomous remediation: pull requests, not tickets

When a fixed version exists, Safeguard's Griffin models open a tested, policy gated pull request that moves the affected dependency to the safe version. You can run it on one repository or in bulk across every project that carries the vulnerable version. Remediation covers eight categories, including dependencies, licenses, security posture and findings. Reachability analysis ranks the vulnerable code paths your application actually calls first, so the first pull requests you merge are the ones that matter most. Autonomous remediation ships switched off, with every category set to Manual, so your team decides how much runs on its own and what still needs a human approval.

2. Gold Open Source: the safe version in seconds, free

gold.safeguard.sh is free and needs no login. Search any CVE or package to see affected and fixed versions, CVSS, EPSS, CISA KEV status and a recommended safe version, or use the "Am I affected?" checker. The same data is available through a free read API, README status badges and RSS threat feeds. The safeguard-sh/gold-scan-action GitHub Action can fail a build on a known exploited vulnerability, on any known vulnerability, or never, so a version you have just removed cannot quietly come back in next week's merge.

3. Gold Registry: start from components that are already fixed

registry.safeguard.sh provides hardened, drop in replacements for upstream packages and container images. Every Gold artifact is held continuously to zero exploitable critical or high CVEs, zero malware and a completed license check, with SLSA L3 or higher provenance, a Sigstore signature and a public evidence page you can hand to an auditor. An artifact loses its Gold badge if it regresses. Teams can point their own package directory at the registry so hardened equivalents resolve without changing how developers install anything.

When upstream has no fix yet, or will not ship one, Safeguard first works to get the fix accepted upstream. If the issue is severe or upstream will not move, Safeguard builds a fork, fixes the vulnerability there and publishes a hardened Gold build, so you are not stuck on a maintainer's release schedule. You can also request a hardened, compatibility tested build of a specific package. This fork and patch flow is in early access today, so ask which packages it covers before you plan around it.

One honest limit: when the vulnerable code lives in a vendor appliance, a hosted service or an operating system runtime you do not build yourself, the fix still has to come from that vendor or distribution. The steps above apply to the open source dependencies and container images you build and ship, which is where most teams carry the largest and least tracked share of their exposure.

See how Safeguard would help

These demos show dependency, image and remediation workflows. None of them touches kernel code, GPU performance or the GPU driver. You can download Safeguard and follow the CLI docs to try them on your own projects.

This demo sets up the Safeguard runner on your own machine, approves it with its device code and scans from the terminal. Pointed at the container image you use for GPU work, the same scan lists the packages and versions inside it, including the Python and system packages layered on top of the CUDA base image, and their known vulnerabilities.

Click through it at your own pace, or open it full screen.

This demo is a real CLI session: ask it to fix the lodash CVE in a repository and watch it read the finding and the advisory, edit package.json to lodash 4.17.21, and run npm install and npm test, asking before any command that changes something and pushing nothing. The same flow applies when a dependency in the Python service that wraps your GEMM kernels needs a version bump.

Click through it at your own pace, or open it full screen.

This demo connects Claude to the Safeguard MCP server with OAuth and lets you choose which tools Claude can use. You can then ask which of your projects ship a given package version, for example a specific framework release used across several inference images, and get the answer from your own scan data.

Click through it at your own pace, or open it full screen.

Quick check

Sources

• NVIDIA CUDA Programming Guide: Writing CUDA kernels, coalesced global memory access and shared memory
• NVIDIA CUDA Programming Guide: C++ language extensions, warp matrix functions (WMMA)
• NVIDIA CUDA Programming Guide: Compute capabilities, register and shared memory limits
• NVIDIA CUDA Programming Guide: Asynchronous data copies and shared memory swizzling
• NVIDIA CUDA Programming Guide: Advanced kernel programming (SIMT)
• NVIDIA CUDA C++ Best Practices Guide: Shared memory in matrix multiplication
• NVIDIA PTX ISA: cp.async, ldmatrix, mma, wgmma and tcgen05
• NVIDIA cuBLAS documentation: math modes and compute types
• NVIDIA nvcc documentation: ptxas options and resource usage
• NVIDIA Nsight Compute Profiling Guide: sections, stall reasons, shared memory banks
• NVIDIA Nsight Compute CLI documentation
• NVIDIA A100 Tensor Core GPU Architecture whitepaper
• NVIDIA Ampere GA102 GPU Architecture whitepaper, v2 (RTX A6000 specs)
• NVIDIA CUTLASS: Efficient GEMM in CUDA
• NVIDIA Technical Blog: CUTLASS, Fast Linear Algebra in CUDA C++
• Simon Boehm: How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance, a Worklog (31 December 2022)

Last reviewed: 25 October 2026. Code samples follow the NVIDIA documentation listed above but were not compiled or executed for this post; measured results are the cited authors' own.