From 306 ms to 86 ms

I implemented C = A × B for square N × N single-precision matrices in CUDA C and timed it from N = 128 to 4096. At N = 4096, each of the 16.8 million elements of C is a sum of 4096 products.

The first kernel is the direct translation: one thread per element of C, in blocks of 16 × 16 threads. Each thread walks its row of A and its column of B in global memory, the GPU’s main memory.

How each kernel was timed

Timer

gettimeofday around one launch and cudaDeviceSynchronize, so copies from host to GPU are excluded

Inputs

the same seeded random matrices for every kernel; C set to zero before each

Runs

one Colab session, one timed launch per kernel and size

Check

the sum of all elements of C, in double precision

The session’s first launch, the naive kernel at N = 128, took 0.159 s, about 300 times its time at N = 256, and nothing like it happened again. I counted it as start-up cost and left N = 128 out of the comparisons.

I did not profile the kernels, so this step is reasoning from the code. The naive kernel loads two values from global memory for every multiply-add: 8,192 loads for each element of C at N = 4096.

Most of those loads repeat. The 16 threads in a row of a block read the same row of A, the 16 in a column the same column of B, and each fetches its own copy. What had to change was not the arithmetic but where its operands came from.

Two changes, each measured before the next.

Tiling

The 256 threads of a block copy a 16 × 16 tile of A and one of B into shared memory, the fast on-chip memory they share, then wait at a barrier. Each thread does 16 multiply-adds from the tiles before the next pair is loaded. Every value fetched now serves 16 threads: 512 global loads per element of C instead of 8,192.

Measured at N = 4096: 306.0 ms to 190.7 ms 1.6× speedup

Register blocking

Each thread computes a 4 × 4 block of C and keeps its 16 running sums in registers, the fastest storage on the GPU. At each step it reads four values of A and four of B from shared memory and multiplies every pair: 16 multiply-adds for 8 reads, where the tiled kernel did 1 for 2.

Measured at N = 4096: 190.7 ms to 86.1 ms 2.2× speedup

That kernel also doubled the tiles to 32 × 32 and cut blocks to 64 threads, which halves the global loads again, to 256 per element. The 2.2× belongs to that combination, not to register blocking alone. And 16 times fewer global loads bought 1.6×, not 16×: they were not the whole cost, and without a profiler I cannot say what the rest was.

The lesson: work out where every operand comes from before changing code, and measure after each change. Next time I would change one thing at a time, so register blocking gets a number of its own, repeat every timing after a warm-up, and profile, so the reasoning about memory becomes a measurement.