CUDA Programming Model

CUDA is the bridge between your algorithm and the GPU hardware. This chapter covers the fundamental execution model — how millions of threads are organized, launched, and scheduled — which is essential for understanding why kernels are fast or slow.

---

1. The Heterogeneous Model: CPU + GPU

CUDA programs run on a heterogeneous system: the CPU (Host) handles control logic and serial code, the GPU (Device) handles massively parallel computation. Think of it as a factory — the CPU is the manager who issues instructions, the GPU is a workshop with thousands of workers who execute in parallel.

Execution Flow

CPU: Initialize data
CPU → GPU: cudaMemcpy (transfer data to GPU)
CPU → GPU: kernel<<<grid, block>>>() (launch computation)
GPU: Thousands of threads execute in parallel
GPU → CPU: cudaMemcpy (transfer results back)
CPU: Post-processing

Every CUDA program follows this pattern: allocate → transfer → compute → transfer back → free.

Function Qualifiers

Qualifier	Runs on	Called from	Purpose
__global__	GPU	CPU (or GPU dynamic parallelism)	Kernel entry point
__device__	GPU	GPU	Helper function called from kernels
__host__	CPU	CPU	Regular CPU function (default)

// Can combine: compile for both CPU and GPU
__host__ __device__ float square(float x) {
    return x * x;
}

---

2. Thread Hierarchy: Grid → Block → Thread

CUDA organizes threads in a three-level hierarchy that maps directly to hardware:

Kernel Launch
└── Grid (all threads for this kernel)
    ├── Block 0
    │   ├── Thread 0, 1, 2, ... 255
    │   └── (scheduled onto one SM)
    ├── Block 1
    │   ├── Thread 0, 1, 2, ... 255
    │   └── (scheduled onto one SM)
    └── Block N
        └── ...

Level	Hardware mapping	Size limit	Can synchronize?
Grid	Entire GPU	Up to 2³¹ - 1 blocks	No (across blocks)
Block	One SM	Up to 1024 threads	Yes (__syncthreads())
Thread	One CUDA Core	—	—
Warp (implicit)	32 consecutive threads in a block	Always 32	Lockstep (SIMT)

Key constraint: Threads within a block can synchronize and share memory. Threads in different blocks cannot communicate (except through global memory with atomics).

Dimensions

Grids and blocks can be 1D, 2D, or 3D:

// 1D: process a vector
dim3 grid(numBlocks);
dim3 block(256);

// 2D: process a matrix
dim3 grid(width / 16, height / 16);
dim3 block(16, 16);  // 256 threads total

// 3D: process a volume
dim3 grid(X/8, Y/8, Z/8);
dim3 block(8, 8, 8);  // 512 threads total

---

3. Thread Indexing

Every thread knows its position through built-in variables:

Variable	Meaning
threadIdx.x/y/z	Thread index within its block
blockIdx.x/y/z	Block index within the grid
blockDim.x/y/z	Number of threads per block (in each dimension)
gridDim.x/y/z	Number of blocks in the grid (in each dimension)

1D Global Index

int idx = blockIdx.x * blockDim.x + threadIdx.x;

For grid of 4 blocks, 256 threads each:

• Block 0, Thread 0 → idx = 0
• Block 0, Thread 255 → idx = 255
• Block 1, Thread 0 → idx = 256
• Block 3, Thread 255 → idx = 1023

2D Global Index (for matrices)

int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
int idx = row * width + col;  // linear index into row-major matrix

Grid-Stride Loop

When there are more data elements than threads, use a grid-stride loop:

__global__ void add(float* a, float* b, float* c, int N) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = blockDim.x * gridDim.x;  // total number of threads

    for (int i = idx; i < N; i += stride) {
        c[i] = a[i] + b[i];
    }
}

This pattern handles any data size with any grid configuration — essential for production kernels.

---

4. Kernel Launch Syntax

// Triple angle brackets: <<<grid_size, block_size, shared_mem, stream>>>
kernel<<<numBlocks, threadsPerBlock>>>(args...);

// With shared memory and stream
kernel<<<grid, block, sharedMemBytes, stream>>>(args...);

Complete Example: Vector Addition

__global__ void vecAdd(float* A, float* B, float* C, int N) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N) {
        C[i] = A[i] + B[i];
    }
}

int main() {
    int N = 1000000;
    float *d_A, *d_B, *d_C;

    // Allocate GPU memory
    cudaMalloc(&d_A, N * sizeof(float));
    cudaMalloc(&d_B, N * sizeof(float));
    cudaMalloc(&d_C, N * sizeof(float));

    // Copy input to GPU
    cudaMemcpy(d_A, h_A, N * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, N * sizeof(float), cudaMemcpyHostToDevice);

    // Launch kernel
    int threadsPerBlock = 256;
    int numBlocks = (N + threadsPerBlock - 1) / threadsPerBlock;
    vecAdd<<<numBlocks, threadsPerBlock>>>(d_A, d_B, d_C, N);

    // Copy result back
    cudaMemcpy(h_C, d_C, N * sizeof(float), cudaMemcpyDeviceToHost);

    // Free GPU memory
    cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
}

Block Size Selection

Block size	Pros	Cons
32 (1 warp)	Simple, no intra-block sync needed	Low occupancy on most GPUs
128	Good balance, 4 warps	—
256	Most common default, 8 warps	—
512	Higher occupancy potential	More register pressure
1024 (max)	Maximum threads per block	Severe register/shared memory pressure

Rule of thumb: Start with 256 threads per block. Profile and adjust only if occupancy analysis shows a better option.

---

5. Warps and Synchronization

Warps: 32 Threads in Lockstep

The GPU groups every 32 consecutive threads in a block into a warp — the actual unit of execution. All 32 threads execute the same instruction each cycle (SIMT). If threads in a warp take different branches, both paths execute sequentially with masked threads — warp divergence.

Warps, divergence, occupancy, and scheduling are covered in depth in Warp Execution & Scheduling.

Block-Level Synchronization

Threads within a block can synchronize using __syncthreads():

__shared__ float tile[256];
tile[threadIdx.x] = globalData[idx];  // all threads write
__syncthreads();                       // WAIT: all threads must reach here
float val = tile[255 - threadIdx.x];  // safe to read others' data

Rules:

• All threads in the block must reach the same __syncthreads() call (no conditional sync)
• Only works within a block — there is no synchronization between blocks (except kernel completion)
• Needed whenever one thread reads data written by another thread in the same block

---

6. Error Handling

CUDA operations are asynchronous — errors may not appear immediately. Always check.

The CUDA_CHECK Macro (Use This Everywhere)

#define CUDA_CHECK(call) do { \
    cudaError_t err = call; \
    if (err != cudaSuccess) { \
        fprintf(stderr, "CUDA error at %s:%d: %s\n", \
                __FILE__, __LINE__, cudaGetErrorString(err)); \
        exit(EXIT_FAILURE); \
    } \
} while(0)

// Usage:
CUDA_CHECK(cudaMalloc(&d_A, size));
CUDA_CHECK(cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice));

// For kernel launches (two checks needed):
kernel<<<grid, block>>>(args);
CUDA_CHECK(cudaGetLastError());        // catches launch config errors
CUDA_CHECK(cudaDeviceSynchronize());   // catches execution errors

Always use this macro in every CUDA project. Silent CUDA errors are debugging nightmares — they produce wrong results without crashing.

---

Self-Check

• Draw the Grid → Block → Thread hierarchy and explain what maps to hardware
• Calculate the global thread index for a 2D grid
• Explain why block size should typically be a multiple of 32
• Implement a grid-stride loop and explain why it’s needed
• Explain warp divergence and its performance cost
• Use __shfl_down_sync to perform a warp-level reduction
• Explain why CUDA operations are asynchronous and how to check for errors

---

References

• CUDA C++ Programming Guide — Thread Hierarchy
• NVIDIA CUDA Samples
• CUDA by Example (Sanders & Kandrot)