CUDA

CUDA is NVIDIA's platform for running general-purpose code on GPUs. You write a function called a kernel, and the GPU runs it in thousands of threads at once. It consists of a language extension to C and C++ (compiled with nvcc), a runtime library, and a driver. Where OpenMP parallelises across CPU cores that already share your address space, and MPI runs separate processes with separate memory, CUDA targets a separate device with its own memory. You have to move data there and back explicitly.

The execution model is SIMT (single instruction, multiple threads). A kernel launch creates a grid of thread blocks, and each block contains up to 1024 threads. Threads are scheduled in groups of 32 called warps, and every thread in a warp executes the same instruction at the same time. Unlike OpenMP, where the runtime picks the thread count for you, the launch geometry is yours to choose on every call.

__global__ void hello(void) {
    printf("block %d thread %d\n", blockIdx.x, threadIdx.x);
}
 
int main(void) {
    hello<<<2, 4>>>();       // 2 blocks of 4 threads = 8 threads total
    cudaDeviceSynchronize(); // wait for the GPU before exiting
    return 0;
}

The <<<blocks, threads>>> syntax is the kernel launch. Inside the kernel, blockIdx, blockDim, and threadIdx tell each thread which piece of the work it owns. The usual pattern is int i = blockIdx.x * blockDim.x + threadIdx.x, which flattens the hierarchy back into a single global index. Compile with nvcc -o prog prog.cu.

A kernel launch is asynchronous. Control returns to the CPU immediately, before the GPU has run anything, which is why the cudaDeviceSynchronize() above is needed for the printf output to appear. This is also the source of the most common beginner bug: a kernel that fails to launch reports no error at the launch site, so a broken program can exit with status 0 and no output at all. See Error handling (CUDA).

A GPU is not faster at any individual operation. It wins by having enormous latency-hiding capacity, so it needs enough parallel work to keep thousands of threads resident. Below that threshold the CPU is usually faster, and for small problems the PCIe transfer alone can cost more than the whole computation. Whether a kernel is limited by arithmetic or by memory bandwidth is the question the Roofline model exists to answer, and on GPUs the answer is almost always memory.

Compile the hello-world above and run it:

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

Here is a real example: SAXPY, the level-1 BLAS operation $y = \alpha x + y$, over 16 million elements.

// 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 (Coalescing (CUDA)).

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.

1. Kernels
2. Thread hierarchy
3. Warps
4. Warp divergence
5. Memory model
6. Global memory
7. Shared memory
8. Constant memory
9. Unified memory
10. Memory transfer
11. Coalescing
12. Bank conflicts
13. Synchronization
14. Atomics
15. Reduction
16. Streams
17. Events
18. Occupancy
19. Error handling
20. Tensor cores

__global__ void f(...)          // kernel: called from host, runs on device
__device__ int  g(...)          // device function: callable only from device code
__host__   int  h(...)          // host function (the default)
__shared__ float tile[256];     // per-block scratchpad, shared by threads in a block
__constant__ float coeff[16];   // read-only, broadcast-optimised device memory
 
f<<<blocks, threads>>>(...)                  // basic launch
f<<<blocks, threads, shmem_bytes>>>(...)     // with dynamic shared memory
f<<<blocks, threads, shmem_bytes, stream>>>  // on a specific stream
 
threadIdx.x   // thread index within its block
blockIdx.x    // block index within the grid
blockDim.x    // threads per block
gridDim.x     // blocks per grid
 
__syncthreads()   // barrier across all threads in a block

cudaMalloc(&ptr, bytes)              // allocate device memory
cudaMallocManaged(&ptr, bytes)       // allocate unified memory, migrated on demand
cudaFree(ptr)                        // release device memory
cudaMemcpy(dst, src, bytes, kind)    // copy; kind is cudaMemcpyHostToDevice etc.
cudaMemcpyAsync(...)                 // same, queued on a stream
 
cudaDeviceSynchronize()              // block the host until the device is idle
cudaStreamCreate(&stream)            // create a stream for concurrent work
cudaStreamSynchronize(stream)        // block until one stream drains
 
cudaEventCreate(&ev)                 // create a timing event
cudaEventRecord(ev, stream)          // mark a point in the stream
cudaEventElapsedTime(&ms, a, b)      // milliseconds between two recorded events
 
cudaGetLastError()                   // fetch and clear the last error
cudaGetErrorString(err)              // human-readable message for an error code
cudaGetDeviceProperties(&prop, dev)  // query core count, bandwidth, limits