Grids, blocks, threads, warps — how to map problem dimensions to launch configurations with worked index calculations.
This tutorial dives into the CUDA thread hierarchy — the fundamental mental model for mapping your data to GPU threads. We cover 1D and 2D indexing with concrete examples, launch configurations, and warp-level execution.
Tutorials 01 — GPU Architecture and 02 — Your First CUDA Kernel. You should be comfortable with the host/device model, __global__ functions, and nvcc compilation.
Every kernel launch creates a grid of blocks, and each block contains threads. This three-level hierarchy can be 1D, 2D, or 3D — you choose the dimensionality that best matches your data.
The simplest case — a flat array of blocks, each containing a flat array of threads. Ideal for 1D data like vectors.
dim3 TypeCUDA uses the dim3 type to specify grid and block dimensions. Unspecified dimensions default to 1.
dim3 block(256); // (256, 1, 1) — 1D block of 256 threads
dim3 block(16, 16); // (16, 16, 1) — 2D block of 256 threads
dim3 block(8, 8, 4); // (8, 8, 4) — 3D block of 256 threads
dim3 grid(10); // (10, 1, 1) — 1D grid of 10 blocks
dim3 grid(10, 10); // (10, 10, 1) — 2D grid of 100 blocks
kernel<<<grid, block>>>(...); // Launch with these dimensionsFor matrix operations, a 2D grid of 2D blocks maps naturally to rows and columns.
Max threads per block: 1024. Max block dimensions: (1024, 1024, 64). Max grid dimensions: (231-1, 65535, 65535). The product of block dimensions must not exceed 1024.
Every thread knows its position within its block (threadIdx.x), which block it belongs to (blockIdx.x), and the block size (blockDim.x). Combining these gives a unique global index.
Suppose we launch a kernel with 4 blocks of 8 threads each to process an array of 32 elements.
| blockIdx.x | threadIdx.x | blockDim.x | Global Index | Calculation |
|---|---|---|---|---|
| 0 | 0 | 8 | 0 | 0 * 8 + 0 = 0 |
| 0 | 5 | 8 | 5 | 0 * 8 + 5 = 5 |
| 1 | 0 | 8 | 8 | 1 * 8 + 0 = 8 |
| 2 | 3 | 8 | 19 | 2 * 8 + 3 = 19 |
| 3 | 7 | 8 | 31 | 3 * 8 + 7 = 31 |
#include <stdio.h>
#include <cuda_runtime.h>
__global__ void vectorAdd(const float *a, const float *b, float *c, int n) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < n) {
c[idx] = a[idx] + b[idx];
}
}
int main(void) {
const int N = 1000000;
size_t bytes = N * sizeof(float);
// Allocate host memory
float *h_a = (float *)malloc(bytes);
float *h_b = (float *)malloc(bytes);
float *h_c = (float *)malloc(bytes);
// Initialise host arrays
for (int i = 0; i < N; i++) {
h_a[i] = 1.0f;
h_b[i] = 2.0f;
}
// Allocate device memory
float *d_a, *d_b, *d_c;
cudaMalloc(&d_a, bytes);
cudaMalloc(&d_b, bytes);
cudaMalloc(&d_c, bytes);
// Copy data to device
cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);
cudaMemcpy(d_b, h_b, bytes, cudaMemcpyHostToDevice);
// Launch kernel: 256 threads per block, enough blocks to cover N
int blockSize = 256;
int gridSize = (N + blockSize - 1) / blockSize; // ceil(N / blockSize)
vectorAdd<<<gridSize, blockSize>>>(d_a, d_b, d_c, N);
// Copy result back to host
cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost);
// Verify result
int errors = 0;
for (int i = 0; i < N; i++) {
if (h_c[i] != 3.0f) errors++;
}
printf("Vector addition: %s (%d errors)\n",
errors == 0 ? "PASSED" : "FAILED", errors);
// Free memory
cudaFree(d_a);
cudaFree(d_b);
cudaFree(d_c);
free(h_a);
free(h_b);
free(h_c);
return 0;
}We launch ceil(N / blockSize) blocks, which may create more threads than elements. The if (idx < n) guard prevents out-of-bounds memory access. This pattern is used in virtually every CUDA kernel.
For matrices, images, and 2D data, use a 2D grid of 2D blocks. Each thread computes its row and column, then converts to a row-major linear index to access the flattened array in memory.
row = blockIdx.y * blockDim.y + threadIdx.y
col = blockIdx.x * blockDim.x + threadIdx.x
A 2×2 grid of 4×4 blocks covers an 8×8 matrix. Each cell shows its (col, row) global coordinate.
#include <stdio.h>
#include <cuda_runtime.h>
__global__ void matrixAdd(const float *a, const float *b, float *c,
int width, int height) {
int col = blockIdx.x * blockDim.x + threadIdx.x;
int row = blockIdx.y * blockDim.y + threadIdx.y;
if (col < width && row < height) {
int idx = row * width + col; // Row-major linear index
c[idx] = a[idx] + b[idx];
}
}
int main(void) {
const int WIDTH = 1024;
const int HEIGHT = 768;
size_t bytes = WIDTH * HEIGHT * sizeof(float);
// Allocate host memory
float *h_a = (float *)malloc(bytes);
float *h_b = (float *)malloc(bytes);
float *h_c = (float *)malloc(bytes);
// Initialise host arrays
for (int i = 0; i < WIDTH * HEIGHT; i++) {
h_a[i] = 1.0f;
h_b[i] = 2.0f;
}
// Allocate device memory
float *d_a, *d_b, *d_c;
cudaMalloc(&d_a, bytes);
cudaMalloc(&d_b, bytes);
cudaMalloc(&d_c, bytes);
// Copy data to device
cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);
cudaMemcpy(d_b, h_b, bytes, cudaMemcpyHostToDevice);
// 2D launch configuration
dim3 block(16, 16); // 256 threads per block
dim3 grid(
(WIDTH + block.x - 1) / block.x, // ceil(1024 / 16) = 64
(HEIGHT + block.y - 1) / block.y // ceil(768 / 16) = 48
);
matrixAdd<<<grid, block>>>(d_a, d_b, d_c, WIDTH, HEIGHT);
// Copy result back
cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost);
// Verify
int errors = 0;
for (int i = 0; i < WIDTH * HEIGHT; i++) {
if (h_c[i] != 3.0f) errors++;
}
printf("Matrix addition (%dx%d): %s (%d errors)\n",
WIDTH, HEIGHT, errors == 0 ? "PASSED" : "FAILED", errors);
// Free memory
cudaFree(d_a);
cudaFree(d_b);
cudaFree(d_c);
free(h_a);
free(h_b);
free(h_c);
return 0;
}The launch configuration (<<<grid, block>>>) determines how many threads run and how they're organised. Choosing well is critical for performance.
The grid must have enough blocks to cover all elements. The standard ceiling-division formula is:
| N (elements) | blockSize | gridSize (blocks) | Total threads | Excess threads |
|---|---|---|---|---|
| 1,000,000 | 256 | 3,907 | 999,936 + 256 = 1,000,192 | 192 |
| 1024 | 256 | 4 | 1,024 | 0 |
| 1000 | 256 | 4 | 1,024 | 24 |
| 100 | 128 | 1 | 128 | 28 |
Because the total thread count often exceeds the data size, every kernel must guard against out-of-bounds access.
#include <stdio.h>
#include <cuda_runtime.h>
// 1D boundary check
__global__ void kernel1D(float *data, int n) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < n) { // Guard against out-of-bounds
data[idx] *= 2.0f;
}
}
// 2D boundary check
__global__ void kernel2D(float *data, int width, int height) {
int col = blockIdx.x * blockDim.x + threadIdx.x;
int row = blockIdx.y * blockDim.y + threadIdx.y;
if (col < width && row < height) { // Guard both dimensions
int idx = row * width + col;
data[idx] *= 2.0f;
}
}
int main(void) {
// 1D example
const int N = 1000;
float *d_data;
cudaMalloc(&d_data, N * sizeof(float));
cudaMemset(d_data, 0, N * sizeof(float));
int blockSize = 256;
int gridSize = (N + blockSize - 1) / blockSize;
kernel1D<<<gridSize, blockSize>>>(d_data, N);
// 2D example
const int W = 1024, H = 768;
float *d_img;
cudaMalloc(&d_img, W * H * sizeof(float));
cudaMemset(d_img, 0, W * H * sizeof(float));
dim3 block2d(16, 16);
dim3 grid2d((W + 15) / 16, (H + 15) / 16);
kernel2D<<<grid2d, block2d>>>(d_img, W, H);
cudaFree(d_data);
cudaFree(d_img);
printf("Launch configurations executed successfully.\n");
return 0;
}Start with 256 threads per block. If you need more shared memory per block, reduce to 128. Only increase to 512 or 1024 if profiling shows it helps. Always use the CUDA occupancy calculator for production code.
The GPU doesn't execute individual threads — it executes warps of 32 threads in lockstep. Understanding warp behaviour is essential for writing fast CUDA code.
A block of 256 threads is divided into 8 warps. Threads are assigned to warps by their linear index within the block:
Warps never span block boundaries. If your block has 100 threads, the last warp (threads 96–127) will have 4 active threads and 28 idle lanes — wasted compute cycles.
When threads in the same warp take different branches of an if/else, the warp must execute both paths serially, masking inactive threads. This is called warp divergence and it halves throughput in the worst case.
#include <stdio.h>
#include <cuda_runtime.h>
// BAD: Divergent — threads within the same warp take different paths
__global__ void divergentKernel(float *data, int n) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < n) {
if (threadIdx.x % 2 == 0) { // Even threads in warp go here
data[idx] = data[idx] * 2.0f;
} else { // Odd threads in warp go here
data[idx] = data[idx] + 1.0f;
}
}
}
// GOOD: Non-divergent — entire warps take the same path
__global__ void nonDivergentKernel(float *data, int n) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < n) {
int warpId = threadIdx.x / 32;
if (warpId % 2 == 0) { // Even warps go here (all 32 threads)
data[idx] = data[idx] * 2.0f;
} else { // Odd warps go here (all 32 threads)
data[idx] = data[idx] + 1.0f;
}
}
}
int main(void) {
const int N = 1024;
size_t bytes = N * sizeof(float);
float *h_data = (float *)malloc(bytes);
for (int i = 0; i < N; i++) h_data[i] = 1.0f;
float *d_data;
cudaMalloc(&d_data, bytes);
cudaMemcpy(d_data, h_data, bytes, cudaMemcpyHostToDevice);
int blockSize = 256;
int gridSize = (N + blockSize - 1) / blockSize;
// Run divergent version
divergentKernel<<<gridSize, blockSize>>>(d_data, N);
cudaDeviceSynchronize();
// Reset and run non-divergent version
cudaMemcpy(d_data, h_data, bytes, cudaMemcpyHostToDevice);
nonDivergentKernel<<<gridSize, blockSize>>>(d_data, N);
cudaDeviceSynchronize();
printf("Warp divergence demo complete.\n");
cudaFree(d_data);
free(h_data);
return 0;
}In the divergent case, threads 0, 2, 4... execute the if branch while threads 1, 3, 5... are masked — then vice versa. Both passes happen; no thread is skipped.
Occupancy is the ratio of active warps to the maximum warps an SM can support. Higher occupancy helps hide memory latency — but maximum occupancy isn't always the goal.
Larger blocks = fewer blocks per SM. If blockSize=1024, only 1–2 blocks can fit on an SM, leaving less scheduling flexibility.
Each thread uses registers. If your kernel uses 64 registers/thread, a 256-thread block needs 16,384 registers — a quarter of a modern SM's 65,536-register file, so at most four such blocks fit.
Shared memory is partitioned among blocks. If each block uses 32 KB and the SM has 96 KB, only 3 blocks can be resident.
| Block Size | Warps/Block | Max Blocks/SM | Active Warps | Occupancy (of 48) |
|---|---|---|---|---|
| 64 | 2 | 16 | 32 | 67% |
| 128 | 4 | 12 | 48 | 100% |
| 256 | 8 | 6 | 48 | 100% |
| 512 | 16 | 3 | 48 | 100% |
| 1024 | 32 | 1 | 32 | 67% |
Example for a GPU with max 48 warps/SM and max 16 blocks/SM. Register and shared memory constraints not shown — they can further reduce these numbers.
The memory hierarchy — global, shared, local, constant, and texture memory — is the single biggest factor in CUDA performance. Occupancy matters, but memory access patterns matter more. We'll cover this in depth next.
Write a kernel that increases the brightness of each pixel by a given amount. Use a 2D grid to match the image dimensions. Clamp values to [0, 255].
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
__global__ void adjustBrightness(unsigned char *image,
int width, int height,
int adjustment) {
int col = blockIdx.x * blockDim.x + threadIdx.x;
int row = blockIdx.y * blockDim.y + threadIdx.y;
if (col < width && row < height) {
int idx = row * width + col;
int newVal = image[idx] + adjustment;
// Clamp to [0, 255]
if (newVal > 255) newVal = 255;
if (newVal < 0) newVal = 0;
image[idx] = (unsigned char)newVal;
}
}
int main(void) {
const int WIDTH = 1920;
const int HEIGHT = 1080;
const int BRIGHTNESS_INCREASE = 50;
size_t bytes = WIDTH * HEIGHT * sizeof(unsigned char);
// Create a test image (grayscale, all pixels at 100)
unsigned char *h_image = (unsigned char *)malloc(bytes);
for (int i = 0; i < WIDTH * HEIGHT; i++) {
h_image[i] = 100;
}
// Allocate and copy to device
unsigned char *d_image;
cudaMalloc(&d_image, bytes);
cudaMemcpy(d_image, h_image, bytes, cudaMemcpyHostToDevice);
// Launch with 2D configuration
dim3 block(16, 16); // 256 threads per block
dim3 grid(
(WIDTH + block.x - 1) / block.x,
(HEIGHT + block.y - 1) / block.y
);
adjustBrightness<<<grid, block>>>(d_image, WIDTH, HEIGHT,
BRIGHTNESS_INCREASE);
// Copy back and verify
cudaMemcpy(h_image, d_image, bytes, cudaMemcpyDeviceToHost);
int errors = 0;
for (int i = 0; i < WIDTH * HEIGHT; i++) {
if (h_image[i] != 150) errors++; // 100 + 50 = 150
}
printf("Brightness adjustment (%dx%d): %s (%d errors)\n",
WIDTH, HEIGHT, errors == 0 ? "PASSED" : "FAILED", errors);
cudaFree(d_image);
free(h_image);
return 0;
}Process the same 1920×1080 image using both a 1D and 2D launch configuration. Verify both produce identical results. Think about which is more natural for image data.
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cuda_runtime.h>
// 1D kernel — treats the image as a flat array
__global__ void brighten1D(unsigned char *image, int totalPixels, int adj) {
int idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < totalPixels) {
int val = image[idx] + adj;
if (val > 255) val = 255;
if (val < 0) val = 0;
image[idx] = (unsigned char)val;
}
}
// 2D kernel — uses row/col indexing
__global__ void brighten2D(unsigned char *image, int w, int h, int adj) {
int col = blockIdx.x * blockDim.x + threadIdx.x;
int row = blockIdx.y * blockDim.y + threadIdx.y;
if (col < w && row < h) {
int idx = row * w + col;
int val = image[idx] + adj;
if (val > 255) val = 255;
if (val < 0) val = 0;
image[idx] = (unsigned char)val;
}
}
int main(void) {
const int W = 1920, H = 1080;
const int TOTAL = W * H;
const int ADJ = 30;
size_t bytes = TOTAL * sizeof(unsigned char);
// Host arrays
unsigned char *h_orig = (unsigned char *)malloc(bytes);
unsigned char *h_result1 = (unsigned char *)malloc(bytes);
unsigned char *h_result2 = (unsigned char *)malloc(bytes);
for (int i = 0; i < TOTAL; i++) {
h_orig[i] = (unsigned char)(i % 256);
}
unsigned char *d_image;
cudaMalloc(&d_image, bytes);
// --- 1D Launch ---
cudaMemcpy(d_image, h_orig, bytes, cudaMemcpyHostToDevice);
int blockSize1D = 256;
int gridSize1D = (TOTAL + blockSize1D - 1) / blockSize1D;
brighten1D<<<gridSize1D, blockSize1D>>>(d_image, TOTAL, ADJ);
cudaMemcpy(h_result1, d_image, bytes, cudaMemcpyDeviceToHost);
// --- 2D Launch ---
cudaMemcpy(d_image, h_orig, bytes, cudaMemcpyHostToDevice);
dim3 block2D(16, 16);
dim3 grid2D((W + 15) / 16, (H + 15) / 16);
brighten2D<<<grid2D, block2D>>>(d_image, W, H, ADJ);
cudaMemcpy(h_result2, d_image, bytes, cudaMemcpyDeviceToHost);
// Compare results
int match = (memcmp(h_result1, h_result2, bytes) == 0);
printf("1D vs 2D results: %s\n", match ? "MATCH" : "DIFFER");
cudaFree(d_image);
free(h_orig);
free(h_result1);
free(h_result2);
return 0;
}Modify the brightness kernel to handle RGB images (3 channels per pixel). Each thread should process all 3 channels of a single pixel. How does this change your indexing?
dim3 type.blockIdx.x * blockDim.x + threadIdx.x gives a unique global index for each thread.This is the most important formula in CUDA. Every kernel starts with computing the global thread index from the built-in variables.
The warp is the true unit of execution. Keep block sizes as multiples of 32. Avoid divergent branches within a warp.
Memory Hierarchy — global, shared, local, constant, and texture memory. Learn how memory access patterns dominate CUDA performance, and how to use shared memory to achieve near-peak bandwidth.