Why I decided to stop treating the GPU as a black box — and what I found when I opened it.
I'm building an autonomous drone that learns to fly using reinforcement learning. The policy is trained in simulation using a custom SAC implementation, and the eventual goal is to deploy that trained policy onto a physical drone I'm building myself.
At some point I asked myself: once that policy runs on hardware, what does inference actually look like at the GPU level? What is PyTorch actually doing when I call torch.matmul? How much performance am I leaving on the table by treating everything as a black box?
That question led me down a rabbit hole — raw CUDA C kernel development. This post documents where I've got to so far.
My home machine runs Ubuntu 24.04 with an NVIDIA GTX 1080 — Pascal architecture, Compute Capability 6.1, 8GB VRAM, 2560 CUDA cores. Not a new card, but a real one. Everything here runs on consumer hardware with no cloud involved.
nvcc verified workingcu121 wheel to match the driver ceiling~/cuda_matmulOne thing worth knowing upfront: Pascal (GTX 10-series) does not have Tensor Cores. Those arrived with Turing (RTX 20-series). Everything here is on raw CUDA cores — no hardware-accelerated FP16 matrix ops. That actually makes it a cleaner learning environment because there's no magic happening underneath. Whatever performance I get, I earned.
Before writing a single line of CUDA C, the mental model you need is the thread/block/grid hierarchy. It sounds simple in diagrams but it only really clicks when you're writing the indexing arithmetic yourself.
When you write a kernel launch like matmul_naive<<<dim3(32,32), dim3(32,32)>>>, you're saying: launch a grid of 32×32 blocks, each containing 32×32 threads. For a 1024×1024 matrix that's 1,048,576 threads — one per output element.
The qualifiers matter too: __global__ runs on GPU called from CPU, __device__ runs on GPU called from GPU, and __shared__ is memory that lives on-chip shared across a block.
The thing nobody tells you upfront is that writing a fast CUDA kernel is mostly a memory access problem, not a compute problem. The GTX 1080 has theoretical compute throughput of ~11 TFLOPS. The bottleneck is almost never arithmetic — it's getting data to and from the right memory tier at the right time.
Every time a thread reads from global memory, it pays ~600ns. The goal of kernel optimisation is to minimise global memory reads by exploiting the faster tiers.
The first kernel is the straightforward approach — each thread computes one element of the output matrix C by walking across a row of A and a column of B:
__global__ void matmul_naive(float *A, float *B, float *C, int n) {
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
if (row < n && col < n) {
float sum = 0.0f;
for (int k = 0; k < n; k++) {
sum += A[row * n + k] * B[k * n + col];
}
C[row * n + col] = sum;
}
}The indexing arithmetic is the first thing to fully understand. blockIdx.y * blockDim.y + threadIdx.y translates a thread's position in the grid into a row index in the matrix. Every thread computes its own row and col, and they all run in parallel.
| Kernel | Time |
|---|---|
matmul_naive | 54.257ms |
| PyTorch / cuBLAS | 0.393ms |
| Gap | 138× |
The fix is to load tiles of A and B into shared memory before doing arithmetic. Instead of each thread independently hammering global VRAM, threads in a block cooperate to load a tile into fast shared memory, then everyone computes from that local copy.
The concept: divide the output matrix into TILE_SIZE × TILE_SIZE tiles. For each tile, load the corresponding slice of A (rows) and B (columns) into shared memory, synchronise all threads in the block, then compute partial sums from shared memory. Repeat across all tiles and accumulate.
__shared__ float tileA[TILE_SIZE][TILE_SIZE];
__shared__ float tileB[TILE_SIZE][TILE_SIZE];
for (int t = 0; t < n / TILE_SIZE; t++) {
tileA[ty][tx] = A[row * n + (t * TILE_SIZE + tx)];
tileB[ty][tx] = B[(t * TILE_SIZE + ty) * n + col];
__syncthreads();
// ... partial sum from shared memory
__syncthreads();
}| Kernel | Time |
|---|---|
matmul_naive | 54.257ms |
matmul_tiled | 33.594ms |
| PyTorch / cuBLAS | 0.393ms |
A 38% improvement. Still 85× off PyTorch. But the gap is closing — and more importantly, the mechanism is now understood. The remaining gap is memory coalescing, register pressure, warp scheduling, and the occupancy tuning that NVIDIA's engineers have spent years optimising in cuBLAS.
When you call torch.matmul on a CUDA tensor, PyTorch calls cuBLAS. The hierarchy looks like this:
cuBLAS is general purpose and highly optimised for each GPU architecture. The reason to write raw kernels isn't to beat it on matmul — it's to understand what it's doing, and to eventually write fused operations that cuBLAS can't do generically. A custom kernel that does matmul + bias add + activation in one pass will beat calling cuBLAS three times separately, because you eliminate intermediate global memory writes between operations.
I'm working through Simon Boehm's "How to Optimize a CUDA Matmul Kernel for cuBLAS-like Performance". It's one of the best technical resources I've found — goes from a naive kernel all the way to near-cuBLAS performance, with benchmark numbers and clear explanations at every step. The progression maps exactly to what I'm doing here.
I'm also recording the learning process for a YouTube channel — the idea being that there's very little content that bridges raw CUDA C and practical deep learning in an applied context. Most engineers use PyTorch as a black box. The channel is about opening it.
This isn't a detour from the drone project — it's a prerequisite. The SAC policy running on the physical drone will run inference on a Raspberry Pi, not a desktop GPU. Understanding the compute primitives at this level means I can make informed decisions about kernel fusion, quantisation, and inference optimisation when it comes to deploying the policy on constrained hardware.
The longer-term goal is to write custom CUDA kernels for the SAC inference loop, then benchmark them against the PyTorch baseline — the same way I'm benchmarking matmul here.
matmul_tiled.cunvprof / Nsight