Step-by-step optimisation from a naive O(n³) global-memory kernel to a shared-memory tiled implementation, with complete compilable code and performance analysis.
This tutorial takes you from a straightforward but slow matrix multiplication to a high-performance tiled implementation using shared memory. We write complete, compilable CUDA programs at each step.
Tutorials 01–04 of this series — you should be comfortable with kernel launches, thread indexing, memory allocation (cudaMalloc/cudaMemcpy), and the basics of shared memory.
Matrix multiplication is the single most important kernel in GPU computing. It sits at the heart of deep learning, scientific simulation, computer graphics, and signal processing.
Every fully-connected layer, every attention head, every convolution (via im2col) is a matrix multiply. Training and inference are dominated by GEMM.
Finite element methods, molecular dynamics, fluid simulation — all rely on dense and sparse linear algebra with MatMul at the core.
Transformations, projections, lighting calculations — every vertex passes through matrix-vector and matrix-matrix products.
For matrices of size N×N, matrix multiplication requires O(N³) arithmetic operations on O(N²) data. This high compute-to-data ratio means there is massive parallelism: each of the N² output elements can be computed independently.
MatMul is embarrassingly parallel at the output level — but memory-bound if you're not careful about data reuse. The gap between a naive and optimised implementation can be 10× or more.
The simplest approach: one thread per output element. Each thread reads an entire row of A and an entire column of B from global memory, computes the dot product, and writes the result to C.
Every element of A is read N times (once per column of C). Every element of B is read N times (once per row of C). That is 2N³ global memory reads for N³ multiply-add operations — an arithmetic intensity of just 0.5 FLOP/byte.
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#include <math.h>
// Error-checking macro
#define CUDA_CHECK(call) do { \
cudaError_t err = call; \
if (err != cudaSuccess) { \
fprintf(stderr, "CUDA error at %s:%d: %s\n", \
__FILE__, __LINE__, cudaGetErrorString(err)); \
exit(EXIT_FAILURE); \
} \
} while (0)
// ── Naive matrix multiply kernel ──
// C = A * B, where A is (M x K), B is (K x N), C is (M x N)
__global__ void matmulNaive(const float *A, const float *B, float *C,
int M, int K, int N) {
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
if (row < M && col < N) {
float sum = 0.0f;
for (int i = 0; i < K; i++) {
sum += A[row * K + i] * B[i * N + col];
}
C[row * N + col] = sum;
}
}
// ── Host helper: initialise matrix with random values ──
void initMatrix(float *mat, int rows, int cols) {
for (int i = 0; i < rows * cols; i++)
mat[i] = (float)(rand()) / RAND_MAX;
}
// ── Host reference multiply for verification ──
void matmulCPU(const float *A, const float *B, float *C,
int M, int K, int N) {
for (int i = 0; i < M; i++)
for (int j = 0; j < N; j++) {
float sum = 0.0f;
for (int k = 0; k < K; k++)
sum += A[i * K + k] * B[k * N + j];
C[i * N + j] = sum;
}
}
int main(void) {
// Matrix dimensions: A(M x K) * B(K x N) = C(M x N)
const int M = 1024, K = 1024, N = 1024;
size_t sizeA = M * K * sizeof(float);
size_t sizeB = K * N * sizeof(float);
size_t sizeC = M * N * sizeof(float);
// Allocate host memory
float *h_A = (float *)malloc(sizeA);
float *h_B = (float *)malloc(sizeB);
float *h_C = (float *)malloc(sizeC);
float *h_ref = (float *)malloc(sizeC);
srand(42);
initMatrix(h_A, M, K);
initMatrix(h_B, K, N);
// Allocate device memory
float *d_A, *d_B, *d_C;
CUDA_CHECK(cudaMalloc(&d_A, sizeA));
CUDA_CHECK(cudaMalloc(&d_B, sizeB));
CUDA_CHECK(cudaMalloc(&d_C, sizeC));
// Copy inputs to device
CUDA_CHECK(cudaMemcpy(d_A, h_A, sizeA, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_B, h_B, sizeB, cudaMemcpyHostToDevice));
// Launch kernel: 16x16 threads per block
dim3 blockDim(16, 16);
dim3 gridDim((N + blockDim.x - 1) / blockDim.x,
(M + blockDim.y - 1) / blockDim.y);
// Warm up
matmulNaive<<<gridDim, blockDim>>>(d_A, d_B, d_C, M, K, N);
CUDA_CHECK(cudaDeviceSynchronize());
// Timed run
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
CUDA_CHECK(cudaEventRecord(start));
matmulNaive<<<gridDim, blockDim>>>(d_A, d_B, d_C, M, K, N);
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
// Copy result back
CUDA_CHECK(cudaMemcpy(h_C, d_C, sizeC, cudaMemcpyDeviceToHost));
// Verify against CPU
matmulCPU(h_A, h_B, h_ref, M, K, N);
float maxErr = 0.0f;
for (int i = 0; i < M * N; i++)
maxErr = fmaxf(maxErr, fabsf(h_C[i] - h_ref[i]));
printf("Naive MatMul (%dx%d x %dx%d)\n", M, K, K, N);
printf(" Time: %.3f ms\n", ms);
printf(" GFLOPS: %.1f\n",
2.0 * M * N * K / (ms * 1e6));
printf(" Max error: %e\n", maxErr);
// Cleanup
cudaEventDestroy(start);
cudaEventDestroy(stop);
cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
free(h_A); free(h_B); free(h_C); free(h_ref);
return 0;
}nvcc -o naive_matmul naive_matmul.cu && ./naive_matmul — the program times the kernel, computes GFLOPS, and verifies against a CPU reference.
Let's quantify why the naive kernel is slow. The bottleneck is not compute — it's memory bandwidth.
For an N×N matrix multiply, each output element requires N multiply-adds (2N FLOP) and reads 2N floats (8N bytes) from global memory.
// Per output element:
FLOPs = 2 * N // N multiplies + N adds
Bytes loaded = 2 * N * 4 // N floats from A row + N floats from B col
Arith intensity = 2N / (8N) = 0.25 FLOP/byte
// For N = 1024:
Total FLOPs = 2 * 1024^3 = 2.15 billion
Total bytes = 2 * 1024^3 * 4 = 8.59 GB // with redundant reads!A modern GPU (e.g., RTX 3090) has ~936 GB/s memory bandwidth and ~36 TFLOPS FP32. The crossover point (ridge point) is at ~38 FLOP/byte. Our naive kernel at 0.25 FLOP/byte is firmly memory-bound — operating at less than 1% of the compute roofline.
| Metric | Value | Notes |
|---|---|---|
| Kernel time | ~5–15 ms | Varies by GPU generation |
| Effective GFLOPS | ~150–400 | Far below peak capability |
| Global loads | 2N³ = 2.1 billion | Massive redundancy |
| Arithmetic intensity | 0.25 FLOP/byte | Memory-bound regime |
| Data reuse | None | Every thread reads independently |
We need to increase data reuse. If a block of threads cooperatively loads a chunk of A and B into fast shared memory, each element is fetched from global memory once and reused many times. This is the core idea behind tiling.
Instead of each thread independently streaming through global memory, we divide the computation into tiles. A tile is a small square sub-matrix (typically 16×16 or 32×32) that fits in shared memory.
To compute one TILE_SIZE×TILE_SIZE block of C, we sweep across K in steps of TILE_SIZE, loading one tile from A's row-band and one from B's column-band each step.
With a tile size of T×T, each element loaded into shared memory is used by T threads. Global memory traffic drops from 2N³ to 2N³/T reads. For T=16, that's a 16× reduction in global memory bandwidth demand.
| Tile Size | Shared Mem per Block | Global Load Reduction | Arith. Intensity |
|---|---|---|---|
| No tiling | 0 bytes | 1× (baseline) | 0.25 FLOP/byte |
| 16 × 16 | 2 KB | 16× | 4.0 FLOP/byte |
| 32 × 32 | 8 KB | 32× | 8.0 FLOP/byte |
The tiled kernel uses __shared__ memory arrays to cache sub-matrices, __syncthreads() barriers to coordinate loading and computation, and boundary checks to handle matrices that aren't multiples of the tile size.
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#include <math.h>
#define TILE_SIZE 16
// Error-checking macro
#define CUDA_CHECK(call) do { \
cudaError_t err = call; \
if (err != cudaSuccess) { \
fprintf(stderr, "CUDA error at %s:%d: %s\n", \
__FILE__, __LINE__, cudaGetErrorString(err)); \
exit(EXIT_FAILURE); \
} \
} while (0)
// ── Tiled matrix multiply kernel ──
// C = A * B, where A is (M x K), B is (K x N), C is (M x N)
__global__ void matmulTiled(const float *A, const float *B, float *C,
int M, int K, int N) {
// Shared memory tiles for A and B
__shared__ float As[TILE_SIZE][TILE_SIZE];
__shared__ float Bs[TILE_SIZE][TILE_SIZE];
int tx = threadIdx.x, ty = threadIdx.y;
int row = blockIdx.y * TILE_SIZE + ty;
int col = blockIdx.x * TILE_SIZE + tx;
float sum = 0.0f;
// Sweep across tiles along the K dimension
int numTiles = (K + TILE_SIZE - 1) / TILE_SIZE;
for (int t = 0; t < numTiles; t++) {
// ── Load tile from A into shared memory ──
int aCol = t * TILE_SIZE + tx;
if (row < M && aCol < K)
As[ty][tx] = A[row * K + aCol];
else
As[ty][tx] = 0.0f;
// ── Load tile from B into shared memory ──
int bRow = t * TILE_SIZE + ty;
if (bRow < K && col < N)
Bs[ty][tx] = B[bRow * N + col];
else
Bs[ty][tx] = 0.0f;
// Wait for all threads to finish loading
__syncthreads();
// ── Compute partial dot product from this tile ──
for (int i = 0; i < TILE_SIZE; i++)
sum += As[ty][i] * Bs[i][tx];
// Wait before loading next tile
__syncthreads();
}
// ── Write result ──
if (row < M && col < N)
C[row * N + col] = sum;
}
// ── Host helper: initialise matrix with random values ──
void initMatrix(float *mat, int rows, int cols) {
for (int i = 0; i < rows * cols; i++)
mat[i] = (float)(rand()) / RAND_MAX;
}
// ── Host reference multiply for verification ──
void matmulCPU(const float *A, const float *B, float *C,
int M, int K, int N) {
for (int i = 0; i < M; i++)
for (int j = 0; j < N; j++) {
float sum = 0.0f;
for (int k = 0; k < K; k++)
sum += A[i * K + k] * B[k * N + j];
C[i * N + j] = sum;
}
}
int main(void) {
// Matrix dimensions: A(M x K) * B(K x N) = C(M x N)
const int M = 1024, K = 1024, N = 1024;
size_t sizeA = M * K * sizeof(float);
size_t sizeB = K * N * sizeof(float);
size_t sizeC = M * N * sizeof(float);
// Allocate host memory
float *h_A = (float *)malloc(sizeA);
float *h_B = (float *)malloc(sizeB);
float *h_C = (float *)malloc(sizeC);
float *h_ref = (float *)malloc(sizeC);
srand(42);
initMatrix(h_A, M, K);
initMatrix(h_B, K, N);
// Allocate device memory
float *d_A, *d_B, *d_C;
CUDA_CHECK(cudaMalloc(&d_A, sizeA));
CUDA_CHECK(cudaMalloc(&d_B, sizeB));
CUDA_CHECK(cudaMalloc(&d_C, sizeC));
// Copy inputs to device
CUDA_CHECK(cudaMemcpy(d_A, h_A, sizeA, cudaMemcpyHostToDevice));
CUDA_CHECK(cudaMemcpy(d_B, h_B, sizeB, cudaMemcpyHostToDevice));
// Launch kernel: TILE_SIZE x TILE_SIZE threads per block
dim3 blockDim(TILE_SIZE, TILE_SIZE);
dim3 gridDim((N + TILE_SIZE - 1) / TILE_SIZE,
(M + TILE_SIZE - 1) / TILE_SIZE);
// Warm up
matmulTiled<<<gridDim, blockDim>>>(d_A, d_B, d_C, M, K, N);
CUDA_CHECK(cudaDeviceSynchronize());
// Timed run
cudaEvent_t start, stop;
CUDA_CHECK(cudaEventCreate(&start));
CUDA_CHECK(cudaEventCreate(&stop));
CUDA_CHECK(cudaEventRecord(start));
matmulTiled<<<gridDim, blockDim>>>(d_A, d_B, d_C, M, K, N);
CUDA_CHECK(cudaEventRecord(stop));
CUDA_CHECK(cudaEventSynchronize(stop));
float ms = 0.0f;
CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
// Copy result back
CUDA_CHECK(cudaMemcpy(h_C, d_C, sizeC, cudaMemcpyDeviceToHost));
// Verify against CPU
matmulCPU(h_A, h_B, h_ref, M, K, N);
float maxErr = 0.0f;
for (int i = 0; i < M * N; i++)
maxErr = fmaxf(maxErr, fabsf(h_C[i] - h_ref[i]));
printf("Tiled MatMul (%dx%d x %dx%d), TILE_SIZE=%d\n",
M, K, K, N, TILE_SIZE);
printf(" Time: %.3f ms\n", ms);
printf(" GFLOPS: %.1f\n",
2.0 * M * N * K / (ms * 1e6));
printf(" Max error: %e\n", maxErr);
// Cleanup
cudaEventDestroy(start);
cudaEventDestroy(stop);
cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
free(h_A); free(h_B); free(h_C); free(h_ref);
return 0;
}__shared__ float As[TILE_SIZE][TILE_SIZE] — declares shared memory visible to all threads in the block.__syncthreads() per tile iteration — one after loading (so all data is ready) and one after computing (so we don't overwrite data another thread still needs).nvcc -o tiled_matmul tiled_matmul.cu && ./tiled_matmul — compare the GFLOPS to the naive version. On most GPUs you'll see a significant speedup.
The following benchmarks were collected on an NVIDIA RTX 3090 (GA102, Ampere). Your results will vary but the relative speedup should be consistent across GPU generations.
| Matrix Size | Naive (ms) | Tiled (ms) | Speedup | Naive GFLOPS | Tiled GFLOPS |
|---|---|---|---|---|---|
| 256 × 256 | 0.12 | 0.06 | 2.0× | 280 | 560 |
| 512 × 512 | 0.85 | 0.22 | 3.9× | 316 | 1,217 |
| 1024 × 1024 | 6.2 | 1.1 | 5.6× | 347 | 1,954 |
| 2048 × 2048 | 48.5 | 7.8 | 6.2× | 355 | 2,203 |
| 4096 × 4096 | 385.0 | 58.2 | 6.6× | 357 | 2,362 |
Larger matrices have more tiles, which means better GPU occupancy and more opportunity for the memory hierarchy to hide latency. The tiled kernel's advantage widens as the compute-to-overhead ratio improves.
| Kernel | Global Mem Reads (1024×1024) | Effective BW | % of Peak (936 GB/s) |
|---|---|---|---|
| Naive | ~8.6 GB | ~700 GB/s | ~75% |
| Tiled (T=16) | ~0.54 GB | ~490 GB/s | ~52% |
The naive kernel saturates memory bandwidth but wastes it on redundant loads. The tiled kernel uses less total bandwidth by reusing data, freeing the compute units to actually do useful work.
Our tiled kernel is a major improvement, but production GEMM libraries like cuBLAS go much further. Here are the main techniques they employ.
Each thread computes a small sub-tile (e.g., 4×4 elements) rather than a single element. This increases register-level reuse: each shared memory load feeds multiple FMA operations.
Use two shared memory buffers: compute on buffer 0 while loading the next tile into buffer 1. This hides global memory latency behind computation.
float4 to load 128 bits at once, improving memory throughput.__shfl_sync() to share data within a warp without shared memory.wmma API for mixed-precision matrix multiply-accumulate (FP16 inputs, FP32 accumulators).#include <cublas_v2.h>
cublasHandle_t handle;
cublasCreate(&handle);
float alpha = 1.0f, beta = 0.0f;
// C = alpha * A * B + beta * C
cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
N, M, K, &alpha,
d_B, N, d_A, K, &beta, d_C, N);
cublasDestroy(handle);cuBLAS on an RTX 3090 achieves ~25+ TFLOPS for large FP32 GEMM — roughly 10× faster than our tiled kernel. Writing your own kernel is great for learning; for production code, use the vendor library.
Hands-on practice to solidify your understanding of tiled matrix multiplication and GPU performance analysis.
Modify the tiled kernel to multiply non-square matrices: A(512 × 768) × B(768 × 1024) = C(512 × 1024). Verify the result against the CPU reference. The boundary handling in our tiled kernel already supports this — but make sure you understand why by tracing through the boundary conditions.
Run the tiled kernel with TILE_SIZE = 8, 16, and 32 on a 2048×2048 matrix. Record the execution time and GFLOPS for each. Which tile size is fastest on your GPU? What limits you from going to TILE_SIZE = 64?
Profile both the naive and tiled kernels using NVIDIA Nsight Compute:
# Profile naive kernel
ncu --set full -o naive_profile ./naive_matmul
# Profile tiled kernel
ncu --set full -o tiled_profile ./tiled_matmul
# Compare specific metrics
ncu --metrics sm__throughput_pct,\
dram__bytes_read.sum,\
l1tex__t_bytes_pipe_lsu_mem_global_op_ld.sum \
./naive_matmulCompare these metrics between the two kernels:
Modify the tiled kernel so each thread computes a 2×2 sub-tile of C instead of a single element. This means each block computes a (2×TILE_SIZE) × (2×TILE_SIZE) output region. Measure the improvement. Hint: you'll need to load the same shared memory tile but accumulate into 4 register variables.
__syncthreads() per tile step are essential — one after loading, one after computing — to prevent data races.| Concept | Description |
|---|---|
__shared__ |
Declares shared memory visible to all threads in a block |
__syncthreads() |
Block-level barrier — all threads must reach it before any proceed |
| Tiling | Decomposing a problem into sub-problems that fit in fast on-chip memory |
| Arithmetic Intensity | FLOP per byte of memory traffic — determines memory-bound vs compute-bound |
| Roofline Model | Framework for understanding performance limits based on compute and memory bandwidth |
Synchronisation & Atomic Operations — explore thread synchronisation beyond __syncthreads(): atomic operations, memory fences, cooperative groups, and patterns for reductions and histograms.