GPU architecture and CUDA programming

Imagine drawing a textured triangle on screen. The renderer starts with three vertices (points in 3D), a camera, lights, and a texture that supplies surface colors. It performs two broad kinds of work:

• Transform the triangle's vertices to determine where its shape lies on the 2D screen.
• For each covered image sample, calculate a candidate color. A covered sample from one triangle is called a fragment; several triangles may produce fragments for the same pixel.

A fragment shader is a small program that computes one such candidate color. This simplified GLSL example reads a texture at coordinates supplied for the fragment, then darkens the color according to the surface's direction relative to a light:

#version 330 core
uniform sampler2D myTexture;  // texture available to every invocation
uniform vec3 lightDir;      // one light direction for this draw
in vec3 norm;                // interpolated normal for this fragment
in vec2 uv;                  // interpolated texture coordinates
out vec4 fragColor;

void main() {
  vec3 kd = texture(myTexture, uv).rgb;
  float brightness = max(dot(normalize(lightDir), normalize(norm)), 0.0);
  fragColor = vec4(kd * brightness, 1.0);  // red, green, blue, alpha
}

The uniform inputs have one value for this draw; norm and uv can differ across fragments. The example is map-like: many invocations perform the same instructions on different inputs. A fragment shader is not required to be a pure function, and one pixel can have several candidate fragments. Depth testing and other graphics operations decide which results contribute to the final image. The abundance of similar, independent candidate-color calculations is what made wide parallel hardware useful.

The shader pattern also describes non-graphics work. If an array has 512 × 512 independent elements, a function could run once for each element just as a fragment shader runs for each image location. Early researchers used this data-parallel match even before GPUs offered a direct compute interface.

The hack

On the graphics interfaces of that period, launching the computation meant posing it as a drawing problem:

1. Put input values in textures and create a 512 × 512 output image to hold the result.
2. Draw two triangles covering that output image, so the graphics pipeline generates a fragment for every location.
3. Write the array operation as a fragment shader. Its invocation for location (i, j) reads the corresponding inputs and writes the corresponding result.

Stanford's Brook language (2004) let programmers write the array view instead. A stream held many values, and a kernel described the operation on one value. In the example below, element k of output_stream becomes scale_amount × input_stream[k]; Brook translated that work into graphics commands and shaders for the available GPUs.

kernel void scale(float amount, float a<>, out float b<>) {
  b = amount * a;
}
float input_stream<1000>;    // stream declarations
float output_stream<1000>;
scale(scale_amount, input_stream, output_stream);  // map the kernel over the streams

Compute mode

The detour through triangles was useful, but it forced a numerical program to speak in graphics terms. NVIDIA's Tesla generation introduced a compute interface in 2007: a program could allocate GPU buffers, copy input data, and launch a function directly over a chosen collection of work items.

CUDA is NVIDIA's programming platform for this compute route. Its C++ extensions let the host CPU name the kernel, choose how many GPU threads to launch, and decide what data they can access. The programmer no longer needs to draw triangles, but does need to understand how thousands of logical threads are grouped and scheduled. That is the next part of the lecture.

Adding two matrices is a good first kernel: output C[j][i] needs only A[j][i] and B[j][i], so each output element can have its own worker. CUDA makes those workers explicit. A kernel launch creates a grid of thread blocks, each containing many CUDA threads. A block and a thread can each have one-, two-, or three-dimensional IDs; two dimensions fit the rows and columns here.

The host CPU chooses a 4 × 3 block, so each block covers twelve output cells. For a 12 × 6 matrix, it needs 3 blocks across and 2 down. The arrays below are flattened in row-major order: element (i, j) is at offset j * Nx + i.

// host code (runs on the CPU)
const int Nx = 12, Ny = 6;
dim3 threadsPerBlock(4, 3);
dim3 numBlocks((Nx + threadsPerBlock.x - 1) / threadsPerBlock.x,
               (Ny + threadsPerBlock.y - 1) / threadsPerBlock.y);
// devA, devB, devC are device allocations; launch 6 blocks, 72 threads
matrixAdd<<<numBlocks, threadsPerBlock>>>(devA, devB, devC, Nx, Ny);

// device code (runs on the GPU)
__global__ void matrixAdd(const float* A, const float* B,
                          float* C, int Nx, int Ny) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  int j = blockIdx.y * blockDim.y + threadIdx.y;
  if (i < Nx && j < Ny) {
    int k = j * Nx + i;
    C[k] = A[k] + B[k];
  }
}

For example, thread (1, 2) in block (2, 1) computes i = 2 × 4 + 1 = 9 and j = 1 × 3 + 2 = 5, so it writes C[5][9]. blockIdx names the block, threadIdx names the thread inside it, and blockDim is the block's size. __global__ marks code the host can launch on the device. The launch queues GPU work and normally returns to the CPU before the threads finish; the host must wait or use an ordered copy before consuming their results.

The programmer also marks where functions can run. __device__ marks a GPU helper callable by device code; an ordinary C++ function runs on the host unless annotated for device use. Next we need to put the matrices somewhere the kernel can read them.

In the explicit-memory version used here, host arrays and device arrays are separate allocations. The CPU allocates GPU storage with cudaMalloc and copies input bytes into it with cudaMemcpy. A pointer returned for ordinary device memory is passed to a kernel; the CPU does not dereference it as if it were a host array:

float* A = new float[N];              // host memory
for (int i = 0; i < N; i++) A[i] = (float)i;
float* deviceA;
cudaMalloc(&deviceA, sizeof(float) * N);    // device memory
cudaMemcpy(deviceA, A, sizeof(float) * N, cudaMemcpyHostToDevice);
// deviceA[i] = ... here would be invalid: deviceA is not a pointer
// into the host's address space

For this workflow, cudaMemcpy plays a role similar to sending a message: bytes move explicitly between host and device storage. CUDA also supports managed memory, where the runtime can migrate data rather than requiring these explicit copies. Inside the GPU, it helps to distinguish three scopes of storage:

Consider a short sliding-window calculation, often called a stencil: output[i] = (input[i] + input[i+1] + input[i+2]) / 3. Output 0 needs inputs 0, 1, and 2; output 1 needs 1, 2, and 3. The two threads therefore ask for some of the same values. First, give each thread one output and let it read all three values from device global memory:

#define THREADS_PER_BLK 128
__global__ void convolve(int N, float* input, float* output) {
  int index = blockIdx.x * blockDim.x + threadIdx.x;  // thread-local
  float result = 0.0f;                                 // thread-local
  for (int i = 0; i < 3; i++)
    result += input[index + i];
  output[index] = result / 3.f;
}
// host: N = 1024 * 1024 outputs (divisible by 128);
// device input has N + 2 elements, device output has N
convolve<<<N/THREADS_PER_BLK, THREADS_PER_BLK>>>(N, devInput, devOutput);

For a block producing outputs 0 through 127, the union of required inputs is only 0 through 129: 130 values. Yet version 1 writes three load operations per thread, or 384 across the block. Version 2 first has the block cooperatively copy that 130-value support region into shared memory. Each thread loads one input, and threads 0 and 1 load the two extra values at the right edge. Then all threads read their three inputs from the shared copy:

__global__ void convolve(int N, float* input, float* output) {
  __shared__ float support[THREADS_PER_BLK+2];     // one copy per block
  int index = blockIdx.x * blockDim.x + threadIdx.x;
  support[threadIdx.x] = input[index];                 // each thread loads one
  if (threadIdx.x < 2)                                  // two threads load the extra 2
    support[THREADS_PER_BLK + threadIdx.x] = input[index + THREADS_PER_BLK];
  __syncthreads();                                    // barrier: all loads done
  float result = 0.0f;
  for (int i = 0; i < 3; i++)
    result += support[threadIdx.x + i];
  output[index] = result / 3.f;
}

The __syncthreads() call is essential to version 2's correctness. Thread 1 loads support[1], which thread 0 needs for output 0; without coordination, thread 0 might read that slot before thread 1 writes it. A barrier makes every thread in the block wait until all have finished the loading phase, then makes those writes visible to the block. The barrier sits outside the if, so all 128 threads reach it. In this example N is divisible by 128; a general kernel with a partial final block would need carefully guarded loads without letting some threads skip the barrier. Unit 5, Lesson 23 revisits the pattern for tiled matrix multiplication.

CUDA's synchronization constructs

• __syncthreads(): participating threads in one block wait for one another, and their prior memory writes become visible at the barrier. It does not synchronize separate blocks.
• Atomic operations, such as atomicAdd, perform one indivisible update to an address in supported shared or global memory.
• Host/device synchronization: a kernel launch normally returns to the CPU before completion. The host can wait with cudaDeviceSynchronize() or cudaStreamSynchronize(); an ordered blocking device-to-host copy can also wait for the result it copies.

Shared memory has a cost beyond the barrier: every resident block reserves space for its own 130-float support array. The next part shows how that resource requirement affects how many blocks a GPU can run at once.

The convolution launch names 1,048,576 CUDA threads in 8,192 blocks. That does not mean a million operating-system threads or a million full stacks must exist at once. It describes a million logical work items. Only blocks currently resident on the GPU consume per-thread register slots and per-block shared-memory space; the rest wait for a turn.

This separation lets the same launch run on GPUs with different numbers of streaming multiprocessors (SMs). The code specifies what each block does, not which SM receives it. A large CUDA grid resembles a data-parallel forall more than a loop calling pthread_create().

What the compiler records

The compiled kernel contains instructions and resource information such as registers used per thread and its static shared-memory allocation. In version 2, support contributes 130 floats, or 520 bytes, of shared memory per resident block. The launch, not the binary, supplies 128 threads per block and the number of blocks. Scheduling must fit both the launch shape and the kernel's resource use.

Think of blocks as jobs waiting for a limited set of workers. An SM takes more blocks only while their combined resource needs fit. As blocks finish, their slots can be reused. CPU task runtimes, including those used with ISPC, can also maintain a worker pool, but CUDA's workers are GPU SMs rather than host threads. Next we look at the hardware inside one of them.

The V100 has 80 streaming multiprocessors (SMs). Each SM is divided into four processing partitions, called sub-cores in this diagram. A partition has a scheduler that selects work from resident warps, arithmetic and memory units, and a share of the SM's registers:

Warps

CUDA groups consecutive threads of a block into warps of 32: threads 0–31 form one warp, 32–63 the next, so a 256-thread block has eight warps. Each V100 partition can keep up to 16 warps resident, subject to the block's register and shared-memory needs. When one warp waits for data, the scheduler can issue ready work from another without saving and restoring an OS thread's stack.

SIMT means single instruction, multiple threads. Threads retain separate values and logical control flow, while a warp can execute a common instruction for several active threads together. If a branch sends threads down different paths, the hardware may execute each path with only its participating threads active; the unused lanes reduce throughput. This is divergence. Volta can schedule threads within a warp more independently than a simple lockstep drawing suggests, but divergent work still has a cost.

The whole SM and the whole chip

These resident counts are upper bounds, not the number a kernel necessarily achieves. The convolution block's registers and shared array also take space. A million-thread launch is still possible: blocks beyond the current resident capacity remain pending until earlier blocks finish. We can now trace that process on a smaller GPU where every allocation is visible.

Return to the convolution block: 128 threads and a 520-byte shared array must be placed together on one SM. To make placement visible, the lecture uses a fictitious two-core GPU. Each core has room for 384 resident threads (12 warps) and 1,536 bytes of shared memory. The real launch has thousands of blocks, but this explorer follows only the first few.

Try version 1: with no shared array, three 128-thread blocks fit on a core. The 256-thread preset needs 1,032 shared bytes, so only one of those larger blocks fits. The limiting resource changes with the kernel. In every case, the source program names blocks rather than SMs; the GPU places pending blocks as space opens. This is why a single launch can adapt to machines of different sizes.

Why must a whole block be resident at once?

Suppose a block has 256 threads but an SM could hold only 128 of them. Running the first half to completion before admitting the second would break a block-wide barrier: the first half could wait at __syncthreads() for threads that have not started. CUDA therefore places a whole block on one SM with resources for all its threads before it runs. The warps then take turns issuing instructions; “resident together” does not mean every thread performs an instruction in the same clock.

Blocks: any order, but not necessarily independent

In an ordinary kernel without an explicit grid-wide synchronization mechanism, blocks may run in any order or at the same time, so one block cannot require another particular block to finish first. That does not forbid every shared update. For example, if a histogram's global counters start at zero and each input is counted once with atomicAdd(&counts[A[i]], 1), the final integer counts do not depend on block order. A different pattern is dangerous:

// thread block 0
// do stuff here
atomicAdd(&myFlag, 1);

// thread block 1
while (atomicAdd(&myFlag, 0) == 0)
  { }            // wait for block 0
// do stuff here

Persistent threads: taking scheduling into your own hands

An advanced persistent-kernel style launches a limited number of long-lived blocks. The GPU scheduler still places those blocks initially; each resident block then loops, claiming new chunks of work from a counter. This moves chunk assignment into software. The V100's 80 SMs × (2,048 maximum threads per SM ÷ 128 threads per block) = 1,280 blocks is only a thread-capacity ceiling: registers, shared memory, and other block limits can lower the number that fits. Persistent kernels must size and coordinate their blocks carefully; residency is not the same as every thread executing simultaneously.

CUDA combines ideas from several models. A grid describes many similar work items. Threads inside a block are SPMD workers that can share data and synchronize. Host-to-device copies in the explicit-memory example resemble message passing. These descriptions apply at different scopes; none alone describes the whole program.

The GPU still offers a graphics pipeline, but a CUDA kernel launch does not pass through its graphics-specific fixed-function stages. It runs on the programmable SMs and uses memory and caches shared with other GPU workloads. The SM also has tensor cores for selected matrix operations; NVIDIA rates the V100 at up to 125 Tensor TFLOP/s for supported mixed-precision work. Unit 6 returns to matrix units. The next time you read a GPU kernel, start with its thread and block mapping, then ask which resources and synchronization rules make that mapping safe and fast.

Shader

A programmable graphics stage. A fragment shader computes candidate results for covered image samples; several fragments may contribute to one pixel.

GPGPU

General-purpose computation on GPUs, originally by disguising it as rendering.

Compute interface

A way to launch GPU kernels directly on data without posing the work as drawing.

Kernel

A __global__ function the host launches over a grid of threads.

Grid, block, thread

The launch's hierarchy: a grid of thread blocks, each containing threads, with 1–3D IDs.

Host, device

The CPU side and GPU side of a CUDA program. Explicit-memory code allocates and copies data between them; managed memory is another option.

Shared memory

Fast on-chip memory shared by the threads of one block.

__syncthreads()

A barrier across all threads of a block.

SM, sub-core

A V100 streaming multiprocessor, divided into four processing partitions with separate warp schedulers and execution units.

Warp, SIMT

A warp groups 32 neighboring CUDA threads; SIMT execution can apply a common instruction to its active threads.

Block scheduler

GPU hardware that assigns thread blocks to SMs as resources allow.

Persistent kernel

A limited set of long-lived blocks that repeatedly claim new chunks of work after hardware schedules them.

Based on Stanford CS149 (Fall 2023), Lecture 7, "GPU Architecture & CUDA Programming", by Kayvon Fatahalian and Kunle Olukotun. The V100 and programming examples draw on the lecture; the CUDA snippets are adapted for this guide. The explanations, diagrams and interactive models are original study material. For CUDA's exact semantics see NVIDIA's CUDA C++ Programming Guide. For the course's slides and videos, see the CS149 course site.