Foundational kernel exercises

Before anyone optimizes a matrix multiply, they work through a shorter list of exercises: transpose, reduction, scan, softmax. None of them is interesting as a computation, which is the point. Each one isolates a single mechanism that the big kernels combine. Transpose is the cleanest example because it does no arithmetic at all, so its performance is purely a question of coalescing. NVIDIA’s classic walkthrough makes the cost concrete: a naive transpose reads its input coalesced but writes its output with a stride of 1024 elements, 4096 bytes, between neighboring threads on a 1024×1024 matrix, and on a Tesla M2050 it reaches 18.8 GB/s of effective bandwidth where a plain copy of the same data reaches 105.2 GB/s.

The fix is the same staging idea the next section applies to matmul: a warp reads a 32×32 tile row by row into shared memory, the block synchronizes at a barrier, and then warps write columns of the tile back out so that the global-memory writes become contiguous again. The barrier matters because threads now consume data that other threads staged, the cooperation discipline that reduction and scan exercises then make the entire kernel. Shared memory brings its own lesson: in a 32×32 tile every element of a column lands in the same memory bank, so reading a column is a worst-case 32-way bank conflict, and the cure is almost comically small: declare the tile 33 elements wide instead of 32 so columns spread across banks. With both fixes the transpose reaches about 95% of copy throughput.

The kernel itself is short. Every kernel in the walkthrough launches blocks of 32×8 threads to move a 32×32 tile, so each thread handles four elements and the index arithmetic is amortized across them. The tiled version reads rows of the input, waits at the barrier, then writes columns of the tile out as rows of the output:

__global__ void transposeCoalesced(float *odata, const float *idata)
{
  __shared__ float tile[TILE_DIM][TILE_DIM];

  int x = blockIdx.x * TILE_DIM + threadIdx.x;
  int y = blockIdx.y * TILE_DIM + threadIdx.y;
  int width = gridDim.x * TILE_DIM;

  for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)
     tile[threadIdx.y+j][threadIdx.x] = idata[(y+j)*width + x];

  __syncthreads();

  x = blockIdx.y * TILE_DIM + threadIdx.x;  // transpose block offset
  y = blockIdx.x * TILE_DIM + threadIdx.y;

  for (int j = 0; j < TILE_DIM; j += BLOCK_ROWS)
     odata[(y+j)*width + x] = tile[threadIdx.x][threadIdx.y + j];
}

On the Tesla M2050 this kernel reaches 51.3 GB/s, up from the naive 18.8 but still half of copy throughput, and the post rules out the obvious suspect with a control experiment: a copy kernel routed through the same shared-memory tile and barrier runs at 104.6 GB/s, essentially full speed. The staging is not the cost. What remains is the bank conflict described above, and the one-line padding fix, tile[TILE_DIM][TILE_DIM+1], takes the transpose to 99.5 GB/s.

Reduction, collapsing an array to a single sum, is the exercise where that cooperation discipline gets optimized end to end. Mark Harris’s NVIDIA walkthrough takes one kernel through seven versions on a G80 GPU whose theoretical bandwidth is 86.4 GB/s, and since a reduction performs one flop per element loaded, bandwidth is the only score that matters. Each block builds a tree in shared memory, halving the number of active threads each step, and a second kernel launch reduces the per-block results, because a kernel launch is CUDA’s only global synchronization point across blocks. The first version reads:

__global__ void reduce0(int *g_idata, int *g_odata) {
  extern __shared__ int sdata[];

  // each thread loads one element from global to shared mem
  unsigned int tid = threadIdx.x;
  unsigned int i = blockIdx.x*blockDim.x + threadIdx.x;
  sdata[tid] = g_idata[i];
  __syncthreads();

  // do reduction in shared mem
  for (unsigned int s=1; s < blockDim.x; s *= 2) {
    if (tid % (2*s) == 0) {
      sdata[tid] += sdata[tid + s];
    }
    __syncthreads();
  }

  // write result for this block to global mem
  if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}

The modulo test looks innocent and is the whole problem: within every warp, which threads pass tid % (2*s) == 0 alternates, so the warps are highly divergent, and the kernel manages 2.083 GB/s on 4M elements. The walkthrough then removes one bottleneck at a time. A strided index makes the branch non-divergent (2.33× faster, but now bank-conflicted), sequential addressing makes shared-memory access conflict-free (4.68× cumulative), doing a first add while loading from global memory stops half the threads idling on the first pass (8.34×), unrolling the last warp drops the barrier and the branch once only 32 threads remain (15.01×), templating the block size unrolls the rest (21.16×), and giving each thread many elements in a grid-strided loop, what Harris calls algorithm cascading, lands at 62.671 GB/s, a 30× cumulative speedup, 73 GB/s on 32M elements.

The deck’s closing arithmetic is the part worth memorizing: of that 30×, the algorithmic changes, addressing and cascading, contributed 11.84×, and the code-level unrolling contributed 2.54×. Fixing how threads cooperate bought almost five times more than fixing how instructions are emitted.

Softmax teaches the streaming trick. The numerically safe version every major framework uses subtracts the vector’s maximum before exponentiating, which costs three passes over the input, the max, the normalizer, then the outputs, four memory accesses per element. The online softmax of Milakov and Gimelshein folds the first two passes into one: carry a running maximum and a running sum together, and whenever a new maximum appears, multiply the sum by e raised to the old max minus the new max before adding the next term. That cuts memory accesses from four to three per element, measured at up to a 1.3× speedup alone and up to 5× fused with top-k. The deeper payoff is structural: a normalizer that can absorb one new element at a time can absorb one new block at a time, which is exactly what the attention section below needs when it walks the score matrix tile by tile without ever holding it whole. The reading list collects worked versions of all four exercises.