Implementing a Naive CUDA GEMM Kernel: Establishing a Baseline

17 Oct, 2026

Building a CUDA matrix multiplication kernel from scratch to understand GPU memory access, arithmetic intensity, and performance optimization.

Introduction: Why CUDA Matrix Multiplication Matters

Matrix multiplication is one of the most important operations in modern machine learning, deep learning, and scientific computing. Neural network layers, attention mechanisms, and many numerical workloads rely on efficient matrix multiplication.

On GPUs, this operation is commonly implemented as General Matrix Multiplication (GEMM).

For matrices $(A\in {ℝ}^{M\times K})and(B\in {ℝ}^{K\times N})$, GEMM computes:

Each output element is calculated as:

Although the mathematics is simple, implementing a fast CUDA GEMM kernel requires careful attention to GPU architecture, memory access patterns, data reuse, and thread organization.

I started Opti-GEMM as a hands-on project to explore these challenges by implementing CUDA matrix multiplication kernels and measuring how different optimization strategies affect performance.

In this first article, I focus on the simplest implementation: a naive CUDA GEMM kernel. The goal is to establish a baseline for correctness and performance before exploring more advanced GPU optimization techniques.

1. Implementing a Naive CUDA GEMM Kernel

The simplest approach to parallelizing matrix multiplication on a GPU is to assign one CUDA thread to each output element of matrix ©.

Each thread calculates a dot product between one row of (A) and one column of (B). It iterates through the reduction dimension (K), accumulates the products, and writes the result to global memory.

The core computation looks like this:

float sum = 0.0f;

for (int k = 0; k < K; k++) {
    sum += A[row * K + k] * B[k * N + col];
}

C[row * N + col] = sum;

Here, row and col identify the output element assigned to the current thread.

Because the matrices use row-major storage, the indexing expressions calculate the memory offsets for the corresponding elements.

A CUDA kernel maps thread and block indices to matrix coordinates:

const int row = blockIdx.y * blockDim.y + threadIdx.y;
const int col = blockIdx.x * blockDim.x + threadIdx.x;

if (row >= M || col >= N) return;

float sum = 0.0f;

for (int k = 0; k < K; k++) {
    sum += A[row * K + k] * B[k * N + col];
}

C[row * N + col] = sum;

The boundary check prevents threads outside the output matrix from accessing invalid memory. This matters because matrix dimensions do not always divide evenly into CUDA block dimensions.

This implementation is straightforward for several reasons:

• Each thread owns one output element.

• The computation directly follows the mathematical definition of matrix multiplication.

• There is no explicit shared-memory allocation.

• Threads do not need block-level synchronization for the calculation.

These properties make naive GEMM an effective starting point for learning CUDA programming and building a reference implementation.

However, simplicity does not necessarily translate into high GPU performance.

2. Understanding CUDA Thread Mapping and Memory Access

To understand the limitations of naive GEMM, we need to examine how CUDA threads access global memory.

CUDA threads execute in groups of 32 called warps. Memory accesses can be more efficient when threads in a warp access nearby memory addresses in a suitable pattern, allowing memory operations to be coalesced.

In this kernel, adjacent threads along the x-dimension calculate adjacent columns of the output matrix.

For a fixed value of (k), those threads access adjacent elements of matrix (B):

Since (B) is stored in row-major order, these elements are contiguous in memory. Output writes to © are also contiguous across adjacent columns.

Therefore, naive GEMM does have a useful memory-access property: its accesses are not inherently uncoalesced across neighboring threads.

The main limitation is how the kernel reuses input data.

3. Why Data Reuse Matters in GPU Matrix Multiplication

Consider the expression for an output element:

When different threads calculate elements of the same output row, they repeatedly need values from the same row of (A). Threads calculating different output rows also need overlapping portions of (B).

The naive CUDA kernel does not explicitly organize these shared input values for reuse. Each thread independently issues its loads and relies on the GPU's memory hierarchy to capture reuse through caches.

Caching can reduce the cost of repeated global-memory requests, but the kernel still performs repeated load instructions and address calculations.

This creates an opportunity for optimization: instead of having each thread independently request all its inputs, threads can cooperate to load data and reuse it across multiple output calculations.

That is the central idea behind tiled matrix multiplication.

4. Arithmetic Intensity and CUDA Kernel Performance

One useful metric for analyzing GPU kernels is arithmetic intensity, which measures the number of floating-point operations performed per byte of data transferred.

For a conventional matrix multiplication, the total computation is approximately:

The factor of two counts one multiplication and one addition for each multiply-accumulate operation.

For the naive kernel, each output thread performs approximately (2K) floating-point operations. It requests (K) elements from each input matrix and writes one output element.

Assuming 32-bit floating-point values and ignoring cache effects, the ratio of arithmetic operations to requested bytes is approximately:

For large (K), this approaches (0.25) FLOP/byte.

This is a simplified estimate of per-thread requested data traffic, not a measurement of actual DRAM traffic. Cache hits and reuse between threads can reduce external-memory transfers.

Nevertheless, the estimate illustrates the opportunity for optimization: a better kernel should perform more useful computation for each unit of data movement.

Shared-memory tiling addresses this problem by loading tiles of the input matrices and reusing those values across multiple output calculations.

5. Benchmarking a CUDA GEMM Kernel

GPU optimization requires reliable measurements. Before comparing different CUDA kernels, it is important to establish a repeatable benchmarking procedure and validate correctness.

For Opti-GEMM, the baseline experiments use an NVIDIA Tesla T4 GPU, based on the Turing architecture with compute capability 7.5.

The recorded benchmark configuration includes:

The kernel is checked against the CPU reference before its performance is evaluated.

Throughput is calculated using:

where (t) is the measured kernel execution time in seconds.

For square matrices of dimension (N), the operation count becomes (2N^3). Reporting both execution time and throughput is useful because the same latency can represent different computational throughput for different matrix sizes.

Benchmark results should also be interpreted in context. GPU clock speeds, thermal conditions, competing workloads, compiler settings, and measurement overhead can affect observed performance. A meaningful comparison requires consistent inputs, precision, correctness criteria, and timing procedures.

6. Naive CUDA GEMM Benchmark Results

The following table summarizes the recorded FP32 benchmark results for Opti-GEMM on the Tesla T4.

Two observations are worth investigating.

First, the naive kernel's measured throughput decreases as the matrix size increases in these results. Larger matrices change the working set and memory behavior, so performance does not necessarily improve simply because the GPU has more work to process.

Second, the register-blocked implementation reaches approximately 1,499 GFLOP/s for a (4096 \times 4096) multiplication, compared with approximately 330 GFLOP/s for naive GEMM.

That is approximately a 4.5× throughput improvement in this recorded experiment.

This comparison demonstrates why data reuse and thread-level work organization deserve attention. It does not establish that the register-blocked implementation will be faster for every workload, or how close it is to a highly optimized library such as NVIDIA cuBLAS. Those questions require further controlled experiments.

7. What Comes Next: Shared-Memory Tiling

The naive CUDA GEMM kernel provides a useful baseline, but its approach to data reuse leaves room for improvement.

The next step in Opti-GEMM is to investigate shared-memory tiling. Threads within a CUDA block will cooperatively load portions of the input matrices into shared memory, allowing those values to be reused across multiple output calculations.

This introduces additional considerations, including synchronization overhead, shared-memory capacity, tile dimensions, and memory access patterns.

Shared-memory tiling is not automatically faster in every situation. Its benefits depend on whether the reduction in data movement outweighs the added overhead for the workload and GPU architecture being tested.

That makes it a useful next experiment: the implementation provides a concrete way to investigate the relationship between data reuse and measured performance.

Conclusion

Implementing a naive CUDA GEMM kernel is a straightforward way to begin exploring GPU programming, but its real value is as a correctness and performance baseline.

This implementation illustrates how CUDA threads map to output elements, why coalesced memory access matters, and how limited explicit data reuse can create opportunities for optimization.

The next stage of Opti-GEMM will investigate shared-memory tiling, followed by more advanced approaches to register-level data reuse and thread organization.

The objective is not simply to make the kernel more complicated. It is to understand which changes improve performance, why they work, and when their overhead outweighs their benefits.

Explore the code: Opti-GEMM on GitHub

This is the first article in a series on CUDA GEMM optimization, GPU memory hierarchies, and performance engineering.

#c++ #cuda #gemm #gpu #performance-engineering