Project 03 — First GPU Kernel

Write CUDA by hand, time it correctly, and find out why your first version loses to the CPU.

1. What you build#

Three kernels in raw CUDA C++, each with a correct timing harness:

1. SAXPY — the “hello world”; purely memory-bound.
2. Naive matmul — one thread per output element.
3. Tiled matmul — shared-memory tiles, the GPU version of Project 02’s blocking.

Then you compare all three against cuBLAS and against your Project 02 numbers.

2. Why it matters#

Until you have launched a kernel yourself, “the GPU” is a magic box that PyTorch talks to. After this project, terms like block, warp, shared memory, coalesced access, and launch overhead refer to code you wrote and bugs you hit.

3. Read first#

• VI.02 — Execution model
• VI.05 — Kernel launch, host and device
• VI.06 — Streams and synchronization
• VI.08 — CUDA programming
• VI.09 — Coalescing

4. Spec#

saxpy(n, a, x, y)                         n up to 2^28
matmul_naive(A, B, C, n)                  n ∈ {256 … 4096}
matmul_tiled(A, B, C, n)  TILE = 16, 32   shared memory per block
Reference: cublasSgemm
Timing:    cudaEvent around the kernel; separate timing for H2D/D2H copies

Report three times for every run: transfer in, kernel, transfer out. Never one blended number.

5. Milestones#

1. Toolchain. nvcc builds and runs; nvidia-smi shows your process.
2. SAXPY. Correct result, then sweep n. Compute effective GB/s (3 × 4 × n bytes: read x, read y, write y) and compare to the card’s rated bandwidth.
3. Small-n crossover. Find the n below which the CPU wins end to end. That number is launch overhead plus PCIe (VI.05, II.09).
4. Naive matmul. Works, and is disappointing. Profile it.
5. Tiled matmul. Load tiles into __shared__ arrays, __syncthreads(), accumulate.
6. Versus cuBLAS. Expect to lose by a wide margin. List what cuBLAS does that you don’t.
7. Pinned memory. Repeat transfers with cudaMallocHost; record the H2D bandwidth change.

6. Starter skeleton#

__global__ void saxpy(int n, float a, const float *x, float *y) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) y[i] = a * x[i] + y[i];
}

#define TILE 16
__global__ void matmul_tiled(const float *A, const float *B, float *C, int n) {
    __shared__ float As[TILE][TILE], Bs[TILE][TILE];
    int row = blockIdx.y * TILE + threadIdx.y;
    int col = blockIdx.x * TILE + threadIdx.x;
    float acc = 0.f;
    for (int t = 0; t < (n + TILE - 1) / TILE; ++t) {
        int ac = t * TILE + threadIdx.x, br = t * TILE + threadIdx.y;
        As[threadIdx.y][threadIdx.x] = (row < n && ac < n) ? A[row * n + ac] : 0.f;
        Bs[threadIdx.y][threadIdx.x] = (br < n && col < n) ? B[br * n + col] : 0.f;
        __syncthreads();
        for (int k = 0; k < TILE; ++k) acc += As[threadIdx.y][k] * Bs[k][threadIdx.x];
        __syncthreads();
    }
    if (row < n && col < n) C[row * n + col] = acc;
}

// Timing that measures the kernel, not the launch
cudaEvent_t t0, t1; cudaEventCreate(&t0); cudaEventCreate(&t1);
cudaEventRecord(t0);
saxpy<<<(n + 255) / 256, 256>>>(n, 2.0f, dx, dy);
cudaEventRecord(t1); cudaEventSynchronize(t1);
float ms; cudaEventElapsedTime(&ms, t0, t1);

7. What to measure#

Add the GPU’s measured bandwidth and peak GFLOP/s to numbers.md next to the CPU’s.

8. Done when#

• All kernels verified against a CPU reference.
• You can show, with numbers, a problem size where the GPU loses and explain it.
• Tiled beats naive, and you can explain it in terms of global-memory reads per output.
• You have opened one run in Nsight Systems or nvprof-equivalent and identified the kernel and the memcpy on the timeline.

9. Common pitfalls#

Timing without synchronizing. Kernel launches are asynchronous. A host timer around the launch measures microseconds no matter how big the job is.

Forgetting the bounds check when n is not a multiple of the block size — silent memory corruption.

Counting the first launch. It includes context creation and module load.

Column-major vs row-major when comparing to cuBLAS. cuBLAS is column-major.

Missing __syncthreads() — a race that only fails sometimes.

10. No GPU? Do this instead#

• Use a free Colab T4: !nvcc -o saxpy saxpy.cu && ./saxpy works in a notebook cell.
• Or run the logic on the Numba CUDA simulator (NUMBA_ENABLE_CUDASIM=1). You get correctness and the mental model, not the timings — borrow timings from a Colab session.

11. Stretch goals#

• Rewrite SAXPY in Triton-lang or as a torch.compiled function and compare the generated code path (preview of XIII.11-12).
• Add an FP16 matmul and look for tensor-core usage in the profiler (VI.11).
• Deliberately break coalescing (stride-n access) and measure the damage.
• Launch 10,000 tiny kernels vs one fused kernel — the motivation for CUDA graphs (VI.07).

12. Interview questions this project answers#

1. Why must you synchronize before reading a GPU timer?
2. Where does the time go for a small GPU workload?
3. What does shared memory buy you in a tiled matmul?
4. Why is SAXPY bandwidth-bound regardless of how many cores the GPU has?

13. Next#

Project 04 — Tiny transformer engine