# CUDA Examples

Working through a **[hello-world kernel](https://docs.nvidia.com/cuda/archive/12.3.1/cuda-c-programming-guide/index.html#simple-code)** to understand the async execution model, then a real-world SAXPY implementation for bandwidth measurement.

## Hello world walkthrough

Compile and run the hello-world kernel from [[cuda]]:

```bash
$ nvcc -o hello hello.cu
$ ./hello
block 1 thread 0
block 1 thread 1
block 1 thread 2
block 1 thread 3
block 0 thread 0
block 0 thread 1
block 0 thread 2
block 0 thread 3
```

Blocks are not ordered. Nothing guarantees block 0 runs before block 1, and nothing guarantees they run at the same time either. The GPU schedules blocks onto its streaming multiprocessors in whatever order it likes, which is exactly why blocks cannot communicate with each other during a kernel.

Now delete the `cudaDeviceSynchronize()` line and rebuild:

```bash
$ nvcc -o hello hello.cu
$ ./hello
$
```

No output, and exit status 0. The process ended before the GPU finished. This is the CUDA equivalent of dropping `-fopenmp`, except that it fails silently instead of falling back to something correct. When a CUDA program produces nothing, check synchronisation and check `cudaGetLastError()` before assuming the kernel logic is wrong.

## SAXPY: a bandwidth-bound kernel

A real example: [[saxpy]], the level-1 BLAS operation $y = \alpha x + y$, over 16 million elements. This kernel shows how to measure performance and identify what limits a real computation.

```c
// compile: nvcc -O2 -o saxpy saxpy.cu
// run: ./saxpy
// description: y = a*x + y on the GPU, timed with CUDA events

#include <stdio.h>

__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];
}

int main(void) {
    int n = 1 << 24;
    size_t bytes = n * sizeof(float);

    float *hx = (float *)malloc(bytes);
    float *hy = (float *)malloc(bytes);
    for (int i = 0; i < n; i++) { hx[i] = 1.0f; hy[i] = 2.0f; }

    float *dx, *dy;
    cudaMalloc(&dx, bytes);
    cudaMalloc(&dy, bytes);
    cudaMemcpy(dx, hx, bytes, cudaMemcpyHostToDevice);
    cudaMemcpy(dy, hy, bytes, cudaMemcpyHostToDevice);

    int threads = 256;
    int blocks = (n + threads - 1) / threads;

    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    cudaEventRecord(start);
    saxpy<<<blocks, threads>>>(n, 2.0f, dx, dy);
    cudaEventRecord(stop);

    cudaMemcpy(hy, dy, bytes, cudaMemcpyDeviceToHost);
    cudaEventSynchronize(stop);

    float ms = 0.0f;
    cudaEventElapsedTime(&ms, start, stop);
    printf("y[0]=%.1f  kernel=%.3f ms  %.0f GB/s\n",
           hy[0], ms, 3.0 * bytes / ms / 1e6);

    cudaFree(dx);
    cudaFree(dy);
    free(hx);
    free(hy);
    return 0;
}
```

The `blocks` calculation rounds up, so the last block is partly idle and the `if (i < n)` guard stops those threads from writing past the end. This guard is boilerplate in almost every CUDA kernel.

The kernel touches three arrays, reading `x`, reading `y`, and writing `y`, which is where the factor of 3 in the bandwidth figure comes from. Compare the number you get against your card's spec sheet. SAXPY does one multiply and one add per 12 bytes moved, so it is entirely memory-bound and should land close to peak bandwidth. If it does not, the access pattern is the first thing to look at ([[cuda-coalescing]]).

Two things worth trying. Vary `threads` from 32 up to 1024 and watch how little it matters once you are past a couple of hundred, since the kernel is bandwidth-bound rather than occupancy-bound. Then move the `cudaEventRecord` calls outside the `cudaMemcpy` calls so the transfers are included in the timing. For a kernel this cheap, the transfers dominate completely, and the GPU version loses to a plain CPU loop.