Working through a hello-world kernel to understand the async execution model, then a real-world SAXPY implementation for bandwidth measurement.

Compile and run the hello-world kernel from CUDA:

$ 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:

$ 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.

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.

// 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.