Let’s build a CUDA matrix multiplication kernel from scratch and optimize it step by step to get as close as possible to cuBLAS. This work-log walks through each kernel i will write.
All code is available in my repo: shehab-ashraf/matmul.cu.
In this post
What is a SGEMM?
SGEMM stands for Single-precision GEneral Matrix-Multiply (C = A × B).
It takes matrix A of shape (M, K) and matrix B of shape (K, N), producing C of shape (M, N). Each element C[i, j] is simply the dot product of row i of A and column j of B:
The L4
All kernels and benchmarks in this post are evaluated on a single NVIDIA L4 GPU.
The specs that matter for this work-log:
Device: NVIDIA L4 (AD104, Ada Lovelace, sm_89)
SMs: 58 | CUDA cores: 7,424 (128/SM) | Tensor cores: 240 (4/SM)
FP32 peak: 30.3 TFLOPS | Memory: 24 GB GDDR6, 300 GB/s
L2: 48 MB | Shared memory: 128 KB/SM (99 KB max per block)
Registers: 65,536 x 32-bit per SM | Max threads/SM: 1,536
cuBLAS implementation
Before writing any code, let’s see what NVIDIA’s cuBLAS can do. We run cublasSgemm on FP32 (no Tensor Cores, no TF32). This is our 100% baseline.
cublasHandle_t handle;
cublasCreate(&handle);
const float alpha = 1.0f;
const float beta = 0.0f;
cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
N, M, K,
&alpha,
matd, N,
Bd, K,
&beta,
Cd, N);
Benchmark result:
>> GPU allocation time: 0.993280 ms
>> Host to device transfer time: 32.254143 ms
------- cublas sgemm kernel ---------
>> Execution time: 9.759744 ms
>> Achieved (GFLOPS): 14082.229803
>> Theoretical max (GFLOPS): 30289.920000
>> Maximum memory bandwidth: 300.048000 GB/s
>> Achieves 46.491472 % of peak GFLOPS
>> Achieves 6.874989 % of peak Memory Bandwidth
---------------------------
>> Device to host transfer time: 14.312992 ms
cuBLAS hits 14,000 GFLOPS (about 46% of peak). This is our target. From here on, every kernel we write will be measured against this number.
Kernel 1: Naive Implementation
Let’s start with the most obvious way. We have an output matrix C. The simplest thing we can do is assign one thread to compute exactly one element of C. To get that single element, the thread walks across a row of A and a column of B, multiplying and adding.
Kernel:
__global__ void sgemm_naive(int M, int N, int K, float alpha, const float *A,
const float *B, float beta, float *C) {
const uint x = blockIdx.x * blockDim.x + threadIdx.x;
const uint y = blockIdx.y * blockDim.y + threadIdx.y;
if (x < M && y < N) {
float tmp = 0.0;
for (int i = 0; i < K; ++i) {
tmp += A[x * K + i] * B[i * N + y];
}
C[x * N + y] = alpha * tmp + beta * C[x * N + y];
}
}
Launch config: one 32x32 block of threads per 32x32 tile of C.
dim3 blockDim(32, 32, 1); // blockDim = (32, 32) -> 1024 threads per block
dim3 gridDim(CEIL_DIV(M, 32), CEIL_DIV(N, 32), 1); // one 32x32 tile of C per block
sgemm_naive<<<gridDim, blockDim>>>(M, N, K, alpha, A, B, beta, C);
One detail of this launch configuration deserves attention. A warp is 32 threads, and each row of a 32x32 block is exactly 32 threads wide, so every block row is exactly one warp. Threads in the same warp therefore have different x coordinates but share the same y coordinate. Keep that in mind: it decides how this kernel talks to memory.
Let’s see how this kernel works under the hood.
Running this kernel results in:
>> GPU allocation time: 0.678912 ms
>> Host to device transfer time: 28.547775 ms
------- v1_naive sgemm kernel ---------
>> Execution time: 605.906921 ms
>> Achieved (GFLOPS): 226.831793
>> Theoretical max (GFLOPS): 30289.920000
>> Maximum memory bandwidth: 300.048000 GB/s
>> Achieves 0.748869 % of peak GFLOPS
>> Achieves 0.110740 % of peak Memory Bandwidth
>> Correctness: PASS (0.0e+00)
---------------------------
>> Device to host transfer time: 14.557728 ms