1D Block Tiling: Making Matrix Multiplication Fast in CUDA

Over the last few days, I have been learning about 1D block tiling and exploring how to push CUDA matrix multiplication beyond basic shared memory tiling. When you first start writing GPU kernels, naive matrix multiplication feels like an easy baseline, until you run a benchmark and realize it barely scratches 5% of what the hardware can do. The standard next step is shared memory tiling, which brings a big improvement by cutting down on slow DRAM reads, but still leaves a huge amount of performance on the table.

While experimenting in my gemm-kernels repository, I wanted to see how far a clean 1D block-tiling approach could go. By having each thread compute 8 outputs at once and reusing values stored directly in GPU registers, the kernel climbed to 1.71 TFLOPS, reaching over 42% of NVIDIA cuBLAS speed.

Here is the story of how 1D block tiling works, why registers are the real game changer, and how the code comes together line by line.

---

1. The Core Idea: Why 1 Output per Thread Leaves Performance on the Table

In a standard tiled matrix multiplication kernel, each thread is responsible for producing exactly one number in the output matrix C.

To compute that single output across several stages:

1. The thread block loads square tiles of input matrices A and B from global DRAM into shared memory.
2. Each thread reads values from shared memory, multiplies them, and adds them into its private accumulator.

Shared memory is much faster than global GPU memory, but it is not infinite in speed. Every multiply-accumulate step still has to issue read requests to shared memory banks.

Modern GPUs have an even faster storage layer that sits right next to the execution units: registers. Registers are essentially free to access in terms of latency and bandwidth compared to any memory buffer on the chip.

The key question is: what if each thread computes several numbers at once instead of just one?

If a thread computes a column strip of 8 numbers in matrix C simultaneously, it can load a single value from matrix B into a register once, and then multiply that single value against 8 different numbers from matrix A. That single register read replaces 8 separate reads from shared memory.

---

2. Kernel Design: Setting Up the Block and Thread Geometry

Let us look at how the tile sizes are chosen:

constexpr int BM = 64; // Rows of C computed by one block
constexpr int BN = 64; // Columns of C computed by one block
constexpr int BK = 8;  // Step size along the K dimension
constexpr int TM = 8;  // Number of output rows handled by each thread

Here is what these constants mean in plain English:

• Block output tile (BM x BN): Each thread block is responsible for calculating a 64 by 64 patch of the final output matrix C (4096 output elements in total).
• Work per thread (TM = 8): Instead of each thread computing 1 element, each thread calculates 8 output elements stacked vertically in a column.
• Threads in a block: Because each thread handles 8 elements, we do not need 4096 threads. We only need (64 / 8) * 64 = 8 * 64 = 512 threads.
• K step (BK = 8): We advance through the shared inner dimension in bites of 8 columns from A and 8 rows from B.

Why 512 Threads Is a Perfect Match

Notice how cleanly the math works out for loading data:

• tile_a has size BM * BK = 64 * 8 = 512 floats.
• tile_b has size BK * BN = 8 * 64 = 512 floats.
• Our thread block has exactly 512 threads.

Because the number of threads matches the number of elements in each tile, every thread loads exactly one float into tile_a and one float into tile_b during each step. There are no complicated loops or strided offsets required just to get the data into shared memory.

---

3. The Complete CUDA Kernel

Here is the full implementation of the 1D block-tiled GEMM kernel, including bounds checking for non-square or arbitrary matrix dimensions:

#include <cuda_runtime.h>

template <const int BM, const int BN, const int BK, const int TM>
__global__ void sgemm_blocktiling_1d_kernel(int num_rows_a, int num_cols_b,
                                            int num_cols_a, float alpha,
                                            const float *matrix_a,
                                            const float *matrix_b, float beta,
                                            float *matrix_c) {
  const uint block_row = blockIdx.x;
  const uint block_col = blockIdx.y;

  // Staging buffers in shared memory (2 KB each = 4 KB total)
  __shared__ float tile_a[BM * BK];
  __shared__ float tile_b[BK * BN];

  // Thread coordinates within the block
  const uint thread_row = threadIdx.x / BN;
  const uint thread_col = threadIdx.x % BN;

  const int global_row = block_row * BM + thread_row * TM;
  const int global_col = block_col * BN + thread_col;

  // Advance pointers to the start of this block's region
  matrix_a += block_row * BM * num_cols_a;
  matrix_b += block_col * BN;
  matrix_c += block_row * BM * num_cols_b + block_col * BN;

  // Private accumulators stored directly in registers
  float thread_results[TM] = {0.0f};

  // Outer loop: walk across the K dimension in chunks of BK
  for (int tile_idx = 0; tile_idx < num_cols_a; tile_idx += BK) {

    // 1. Cooperative load into tile_a (1 thread loads 1 float)
    const uint a_row = threadIdx.x / BK;
    const uint a_col = threadIdx.x % BK;
    if ((block_row * BM + a_row) < num_rows_a &&
        (tile_idx + a_col) < num_cols_a) {
      tile_a[a_row * BK + a_col] = matrix_a[a_row * num_cols_a + a_col];
    } else {
      tile_a[a_row * BK + a_col] = 0.0f;
    }

    // 2. Cooperative load into tile_b (1 thread loads 1 float)
    const uint b_row = threadIdx.x / BN;
    const uint b_col = threadIdx.x % BN;
    if ((tile_idx + b_row) < num_cols_a &&
        (block_col * BN + b_col) < num_cols_b) {
      tile_b[b_row * BN + b_col] = matrix_b[b_row * num_cols_b + b_col];
    } else {
      tile_b[b_row * BN + b_col] = 0.0f;
    }

    // Wait until all 512 threads finish loading their elements
    __syncthreads();

    matrix_a += BK;
    matrix_b += BK * num_cols_b;

    // 3. Compute loop: multiply and accumulate using register caching
    for (uint dot_idx = 0; dot_idx < BK; ++dot_idx) {
      // Load one value from tile_b into a register
      float b_tmp = tile_b[dot_idx * BN + thread_col];

      // Reuse that single value across all 8 output rows
      for (uint res_idx = 0; res_idx < TM; ++res_idx) {
        thread_results[res_idx] +=
            tile_a[(thread_row * TM + res_idx) * BK + dot_idx] * b_tmp;
      }
    }

    // Wait before the next iteration overwrites the shared memory tiles
    __syncthreads();
  }

  // 4. Write back the final results to global memory
  for (uint res_idx = 0; res_idx < TM; ++res_idx) {
    int row = global_row + res_idx;
    if (row < num_rows_a && global_col < num_cols_b) {
      matrix_c[(thread_row * TM + res_idx) * num_cols_b + thread_col] =
          alpha * thread_results[res_idx] +
          beta *
              matrix_c[(thread_row * TM + res_idx) * num_cols_b + thread_col];
    }
  }
}

---

4. Understanding How the Code Works

Let us break down the critical sections of this kernel to see why it runs smoothly.

Register Accumulation

The line float thread_results[TM] = {0.0f}; allocates 8 private floating point variables for each thread. Because TM is a constant known at compile time, the CUDA compiler places these 8 numbers directly into physical hardware registers. They never touch DRAM or cache during the entire computation loop.

The 8x Shared Memory Reuse Trick

Look closely at the inner compute loops:

for (uint dot_idx = 0; dot_idx < BK; ++dot_idx) {
  float b_tmp = tile_b[dot_idx * BN + thread_col];
  for (uint res_idx = 0; res_idx < TM; ++res_idx) {
    thread_results[res_idx] +=
        tile_a[(thread_row * TM + res_idx) * BK + dot_idx] * b_tmp;
  }
}

In every pass of dot_idx, the thread reads b_tmp once from tile_b. Then, that exact same b_tmp value is multiplied against 8 separate elements of tile_a.

If we had 8 separate threads doing this work independently, each one would have read b_tmp from shared memory on its own. By combining those 8 calculations into one thread, we reduce the shared memory traffic for matrix B by a factor of 8.

Tiny Memory Footprint Means Great Scheduling

Each block uses:

• tile_a: 512 floats (2048 bytes)
• tile_b: 512 floats (2048 bytes)

Together, that is only 4 KB of shared memory per block. Because modern streaming multiprocessors have tens to hundreds of kilobytes of shared memory, this small footprint allows the GPU warp scheduler to run multiple blocks concurrently on each multiprocessor. When one block is waiting at a __syncthreads() barrier, the hardware instantly switches to another block that is ready to compute.

---

5. Performance Benchmarks

To see the real-world impact of 1D block tiling, the kernel was tested against square matrices of increasing sizes using FP32 precision, with alpha = 1.0 and beta = 0.0. We compare it directly against NVIDIA cuBLAS:

Matrix Size (M = N = K)	Custom Kernel Time	Custom TFLOPS	cuBLAS Time	cuBLAS TFLOPS	% of cuBLAS
512	0.375 ms	0.72	0.132 ms	2.03	35.3%
1024	2.443 ms	0.88	0.752 ms	2.86	30.8%
2048	10.019 ms	1.71	4.254 ms	4.04	42.5%
4096	82.786 ms	1.66	34.473 ms	3.99	41.6%

What the Numbers Tell Us

• Small matrices (512 to 1024): The GPU does not have enough total work to completely saturate all execution pipelines, reaching roughly 0.7 to 0.9 TFLOPS.
• Large matrices (2048 to 4096): As matrix size grows, all multiprocessors are kept completely busy. Throughput stabilizes around 1.66 to 1.71 TFLOPS, which is 41.6% to 42.5% of cuBLAS.

Considering that naive matrix multiplication runs at less than 5% of peak capability, hitting over 42% of an industry standard library using clean, straightforward CUDA C++ proves the power of register tiling.