GPU Programming

GPUs (Graphics Processing Units) are throughput-oriented processors designed to apply similar operations across large numbers of data elements in parallel. They were originally built for graphics workloads, where millions of pixels or vertices must be processed, but the same architecture is also effective for scientific simulation, image processing, machine learning, numerical computing, and other data-parallel workloads.

What makes GPU programming different?

• GPUs contain many execution lanes organized for massive parallelism.
• They use a specialized memory hierarchy that includes global memory, shared memory, registers, caches, constant memory, texture memory, and banked memory structures.
• Programming models such as CUDA, HIP, OpenCL, and SYCL let developers launch large numbers of lightweight threads or work-items.
• High performance depends on memory access patterns, thread organization, occupancy, synchronization, and limiting unnecessary branch divergence.

A GPU is fast when many threads perform similar work on different data with predictable memory access.

Time-Space Duality

In parallel architecture, time-space duality describes two complementary ways to achieve high throughput:

1. Use many replicated processing elements at the same time (spatial parallelism).
2. Reuse fewer processing elements over time through a pipeline (temporal parallelism).

Modern GPUs combine both ideas: they contain many replicated execution resources, and those resources are also pipelined.

====================================================
ARRAY PROCESSOR: Spatial Parallelism
====================================================
         |   T0   |   T1   |   T2   |   T3   |
----------------------------------------------------
 PE0     |  LD0   |  ADD0  |  MUL0  |  ST0   |
 PE1     |  LD1   |  ADD1  |  MUL1  |  ST1   |
 PE2     |  LD2   |  ADD2  |  MUL2  |  ST2   |
 PE3     |  LD3   |  ADD3  |  MUL3  |  ST3   |
----------------------------------------------------
Many processing elements work at the same time.

====================================================
VECTOR PIPELINE: Temporal Parallelism
====================================================
Time -->    T0     T1     T2     T3     T4     T5     T6
----------------------------------------------------------------
 Elem0:     LD0    ADD0   MUL0   ST0
 Elem1:            LD1    ADD1   MUL1   ST1
 Elem2:                   LD2    ADD2   MUL2   ST2
 Elem3:                          LD3    ADD3   MUL3   ST3
----------------------------------------------------------------
Each element moves through pipeline stages over time.

Analogy: Imagine a bakery that must make thousands of identical cookies before sunrise. One strategy is to give many bakers identical workstations so they can perform the same step on different batches at the same time. This resembles an array processor, where many processing elements work side by side. Another strategy is to use a pipelined cookie-making machine: while one batch is being shaped, the next can already be mixed and another can be loaded. This resembles a vector pipeline, where hardware stages are reused over time. A modern GPU combines both ideas by providing many parallel execution resources, each backed by deeply pipelined hardware.

I. Array Processors: Spatial Parallelism

An array processor contains multiple processing elements that work on different data elements, often under a common instruction stream.

• Many elements are processed at the same time.
• Throughput is high when all processing elements have useful work.
• Hardware cost is higher because more execution resources are replicated.
• This model works well for image filters, vector arithmetic, matrix operations, and stencil computations.

II. Vector Processors: Temporal Parallelism

A vector processor applies one vector instruction to many elements. The hardware may pipeline the operation so that a new element enters the pipeline every cycle.

• The same functional units are reused across elements over time.
• Loop overhead is reduced because one instruction represents many scalar operations.
• Performance is best when memory access is regular and predictable.
• Startup cost matters: short vectors may not fully utilize the pipeline.

How This Relates to GPUs

A modern GPU is not just an array processor or just a vector processor. It combines:

• Spatial parallelism: many streaming multiprocessors / compute units.
• Temporal parallelism: pipelined execution inside each unit.
• Thread-level parallelism: many resident warps or wavefronts hide latency.
• SIMT execution: programmers write scalar thread code, while hardware executes groups of threads in SIMD-like fashion.

Vector Processor

A vector processor executes operations on vectors, which are ordered collections of elements. Instead of issuing one scalar instruction per element, a vector instruction expresses the same operation across many elements.

For example, a vector add instruction can represent:

$$ C[i] = A[i] + B[i] \quad \text{for many values of } i $$

Important vector-processor concepts include:

Benefits

• Fewer loop-control instructions.
• High arithmetic throughput on regular data.
• Efficient memory access for contiguous or predictable strides.
• Less scheduling complexity when elements are independent.

Limitations

• Poor fit for irregular control flow.
• Poor fit for pointer-heavy or dependency-heavy code.
• Performance drops when memory access is scattered or badly aligned.

Strides

A stride is the distance in memory between consecutive elements of a logical vector.

• Stride 1: consecutive elements are adjacent in memory.
• Stride k: consecutive elements are separated by k elements in memory.

Stride matters because regular, contiguous memory access is usually more efficient than scattered access. On GPUs, adjacent threads should ideally access nearby addresses so the hardware can combine their requests into a small number of memory transactions.

Row-Major Layout Example

Consider a $4 \times 4$ matrix $A$ stored in row-major order:

$$ A = \begin{bmatrix} A(0,0) & A(0,1) & A(0,2) & A(0,3) \\ A(1,0) & A(1,1) & A(1,2) & A(1,3) \\ A(2,0) & A(2,1) & A(2,2) & A(2,3) \\ A(3,0) & A(3,1) & A(3,2) & A(3,3) \end{bmatrix} $$

In row-major memory, the layout is:

A(0,0), A(0,1), A(0,2), A(0,3),
A(1,0), A(1,1), A(1,2), A(1,3),
A(2,0), A(2,1), A(2,2), A(2,3),
A(3,0), A(3,1), A(3,2), A(3,3)

So in a row-major $N \times M$ matrix, moving across a row has stride 1, while moving down a column has stride M.

Matrix Multiplication Example

For matrix multiplication,

$$ C = A \times B $$

where $A$ is $m \times n$, $B$ is $n \times p$, and $C$ is $m \times p$:

$$ C(i,j) = \sum_{k=0}^{n-1} A(i,k) \cdot B(k,j) $$

To compute one element $C(i,j)$:

I. Load row $i$ of $A$:

$$ [A(i,0), A(i,1), \dots, A(i,n-1)] $$

In row-major order, this is contiguous, so the stride is 1.

II. Load column $j$ of $B$:

$$ [B(0,j), B(1,j), \dots, B(n-1,j)] $$

In row-major order, this is not contiguous. Since $B$ has $p$ columns, the stride is p.

III. Multiply corresponding elements and accumulate the sum.

GPU performance note: A single thread walking down a column of a row-major matrix performs strided accesses. Whether those accesses are coalesced at the GPU level also depends on how neighboring threads map to matrix rows and columns. Optimized matrix multiplication usually uses tiling so that blocks load contiguous tiles of $A$ and $B$ into shared memory and reuse those values many times.

Loading and Storing Vectors from Memory

Vector loads and stores move multiple elements between memory and registers. A closely related GPU concern is coalesced memory access, where threads in the same warp access nearby addresses that the hardware can serve with a small number of memory transactions.

Why Memory Layout Matters

Good pattern: consecutive threads read consecutive elements
Thread:   T0   T1   T2   T3   T4   T5   T6   T7
Address:  A0   A1   A2   A3   A4   A5   A6   A7

Poor pattern: consecutive threads read far-apart elements
Thread:   T0   T1   T2   T3   T4   T5   T6   T7
Address:  A0   A8   A16  A24  A32  A40  A48  A56

Efficient GPU programs often organize data so adjacent threads access adjacent elements. This is one reason an array of structures is sometimes converted into a structure of arrays.

// Array of Structures: convenient, but fields may be interleaved
typedef struct {
  float x, y, z;
} Particle;
Particle particles[N];

// Structure of Arrays: often better for coalesced GPU access
float x[N];
float y[N];
float z[N];

Memory Banking

Memory banking divides storage into multiple banks so several accesses can be serviced in parallel. The idea appears in vector machines, GPU shared memory, register files, caches, and memory controllers.

A simplified banked-memory system looks like this:

                 +----------------+
                 |  Processor /   |
                 | Memory Control |
                 +--------+-------+
                          |
        +-----------------+-----------------+
        |                 |                 |
        v                 v                 v
   +---------+       +---------+       +---------+      ...
   | Bank 0  |       | Bank 1  |       | Bank 2  |
   +---------+       +---------+       +---------+

If consecutive addresses map to different banks, the system can serve multiple accesses efficiently:

Address:  0   1   2   3   4   5   6   7
Bank:     0   1   2   3   0   1   2   3

If many threads access different addresses that map to the same bank, a bank conflict may occur. The hardware may need to serialize some of those accesses, reducing throughput.

Example: Shared Memory Bank Conflict

Good: each thread hits a different bank
T0 -> Bank 0
T1 -> Bank 1
T2 -> Bank 2
T3 -> Bank 3

Bad: several threads hit the same bank
T0 -> Bank 0
T1 -> Bank 0
T2 -> Bank 0
T3 -> Bank 0

Banked memory increases bandwidth only when access patterns spread requests across banks. Poor strides can create bank conflicts.

Vectorizable Loops

A loop is vectorizable when its iterations are sufficiently independent to execute in parallel without changing the program's result.

Vectorizable Loop

for (int i = 0; i < n; ++i) {
  C[i] = A[i] + B[i];
}

Each iteration reads A[i] and B[i] and writes C[i]. No iteration depends on another, so the loop can be vectorized or mapped to GPU threads.

Not Safely Vectorizable Without Transformation

for (int i = 1; i < n; ++i) {
  C[i] = C[i - 1] + A[i];
}

Here, C[i] depends on C[i - 1], which was computed by the previous iteration. This loop has a loop-carried dependency.

Average Example

for (int i = 0; i < 50; ++i) {
  C[i] = (A[i] + B[i]) / 2;
}

This is vectorizable because each output element is independent. For integer arrays, a shift such as (A[i] + B[i]) >> 1 may look like an alternative to division by 2, but modern compilers often optimize division by powers of two automatically, and the two expressions are not always equivalent:

• Signed right shifts and integer division can differ for negative values.
• A[i] + B[i] can overflow before either the shift or division.
• The shift form is not appropriate for floating-point averages.

A safer integer average may require a wider intermediate type or a formula that avoids overflow.

GPUs Are SIMT Engines

GPUs are often described as SIMD-like, but many GPU programming models are better described as SIMT: Single Instruction, Multiple Threads.

In SIMT:

• The programmer writes code as if each thread is scalar and independent.
• The hardware groups threads into fixed execution groups: warps on NVIDIA GPUs and wavefronts/waves on AMD GPUs.
• Threads in the same group normally issue instructions together across active lanes.
• If threads diverge into different branches, the hardware tracks active lanes and executes the required paths separately, reducing efficiency.

Good GPU Workload

C[i] = A[i] + B[i];

All threads perform the same operation on different elements.

Less Efficient GPU Workload

if (A[i] > threshold) {
  C[i] = expensive_path(A[i]);
} else {
  C[i] = cheap_path(A[i]);
}

If neighboring threads choose different branches, the warp may have to run both paths with different active lanes.

SIMD vs. SIMT

SIMD

===============================================
SIMD: Single Instruction, Multiple Data
===============================================
          +------------------------+
          |  Instruction Fetch     |
          |  one instruction       |
          +-----------+------------+
                      |
                      v
    +-----------+-----------+-----------+-----------+
    |  Lane 0   |  Lane 1   |  Lane 2   |  Lane 3   |
    +-----------+-----------+-----------+-----------+

A single vector instruction drives several data lanes.

SIMT

===============================================
SIMT: Single Instruction, Multiple Threads
===============================================
              +----------------------+
              |   Warp Scheduler     |
              +----------+-----------+
                         |
                         v
    +-----------+-----------+-----------+-----------+
    | Thread 0  | Thread 1  | Thread 2  | Thread 3  |
    | scalar PC | scalar PC | scalar PC | scalar PC |
    +-----------+-----------+-----------+-----------+

Each thread has its own logical execution state, while the hardware schedules and issues work in groups of threads.

Note: Warps are fixed groups of neighboring threads, not arbitrary groups that are dynamically re-formed each cycle. Divergence is handled through per-lane execution state and reconvergence mechanisms rather than by regrouping matching threads across the GPU.

Fine-Grained Multithreading

Fine-grained multithreading interleaves instructions from many resident hardware thread groups to hide latency. GPUs use this idea heavily.

When one warp stalls on a long-latency operation, the scheduler can issue instructions from another ready warp without performing a heavyweight software context switch.

Cycle:       0      1      2      3      4      5
Issued:   Warp0  Warp1  Warp2  Warp3  Warp0  Warp1

This improves throughput because execution units spend less time idle.

Benefits

• Hides memory latency.
• Keeps pipelines busy.
• Reduces the need for complex out-of-order execution.
• Works well when many independent warps are resident.

Costs

• Each resident thread needs registers and state.
• Single-thread latency may be poor compared with CPUs.
• Too many registers per thread can reduce occupancy.
• If all resident warps stall, the GPU still cannot issue work.

High-Level GPU Architecture

A modern GPU is organized around many compute units, commonly called Streaming Multiprocessors (SMs) on NVIDIA GPUs and Compute Units (CUs) on AMD GPUs.

+------------------------------------------------------------------+
|                              GPU                                 |
|                                                                  |
|  +----------------+   +----------------+   +----------------+    |
|  | SM / CU        |   | SM / CU        |   | SM / CU        |    |
|  | - schedulers   |   | - schedulers   |   | - schedulers   |    |
|  | - registers    |   | - registers    |   | - registers    |    |
|  | - shared mem   |   | - shared mem   |   | - shared mem   |    |
|  | - ALUs / lanes |   | - ALUs / lanes |   | - ALUs / lanes |    |
|  +-------+--------+   +-------+--------+   +-------+--------+    |
|          |                    |                    |             |
|          +--------------------+--------------------+             |
|                               |                                  |
|                       +-------v-------+                          |
|                       |      L2       |                          |
|                       |     Cache     |                          |
|                       +-------+-------+                          |
|                               |                                  |
|                  +------------v-------------+                    |
|                  | Memory Controllers       |                    |
|                  | GDDR / HBM Device Memory |                    |
|                  +--------------------------+                    |
+------------------------------------------------------------------+

Main Components

General-Purpose Programming on GPUs

General-purpose GPU programming (GPGPU) uses GPUs for non-graphics computation. It works best when a workload exposes substantial data parallelism.

Good GPU workloads usually have:

• Many independent elements.
• High arithmetic intensity or coalesced memory access.
• Limited synchronization.
• Regular control flow.
• Enough work to amortize launch overhead and, on discrete GPUs, host-device transfer costs.

Poor GPU workloads often have:

• Heavy or highly irregular branching.
• Small input sizes.
• Irregular pointer chasing.
• Strong sequential dependencies.
• Frequent CPU-GPU transfers on systems with separate host and device memory.

Typical GPU Execution Flow

+-------------------------+       1. Copy input        +------------------------+
|      CPU / Host         | -------------------------> |      GPU / Device      |
|                         |                            |                        |
| Runs host code          |       3. Copy output       | Executes kernels       |
| Allocates device memory | <------------------------- | Stores device data     |
+-------------------------+                            +-----------+------------+
                                                               ^
                                                               |
                                                       2. Kernel launch

I. Copy input data to the GPU

• Allocate device memory.
• On a discrete-memory system, transfer input arrays from host memory to device memory.

II. Launch a GPU kernel

• The CPU launches a GPU function called a kernel.
• The kernel runs across many blocks and threads.

III. Copy results back to the CPU

• Transfer output data from device memory back to host memory.
• Free device memory when finished.

Performance note: On discrete GPUs, transfers across PCIe or other host-device links can be expensive. Reusing data on the GPU across multiple kernels is often faster than repeatedly moving it between host and device memory.

Warps and Blocks

CUDA organizes GPU work into a hierarchy:

Grid
└── Block 0
    ├── Warp 0: threads 0-31
    ├── Warp 1: threads 32-63
    └── ...
└── Block 1
    ├── Warp 0
    ├── Warp 1
    └── ...

A kernel launch specifies a grid of blocks, and each block contains threads.

kernel<<<numBlocks, threadsPerBlock>>>(args);

Grid

The grid is the full set of blocks launched for one kernel. It represents the whole problem.

Block

A block is a group of threads that can cooperate.

• Threads in the same block can share __shared__ memory.
• Threads in the same block can synchronize with __syncthreads().
• Blocks are normally scheduled independently, so a conventional kernel should not rely on execution order or direct synchronization between blocks.

Thread

A thread is the basic programming unit. Each thread has:

• Its own registers.
• Its own local variables.
• Built-in indices such as threadIdx, blockIdx, and blockDim.

Warp / Wavefront

A warp is a hardware execution group. On NVIDIA GPUs, a warp contains 32 threads. On AMD GPUs, the corresponding group is usually called a wavefront or wave, and its width depends on the architecture and execution mode.

+-------------------------------------------------------------+
|                         GPU Grid                            |
|                                                             |
|  +------------------+    +------------------+               |
|  |    Block 0       |    |    Block 1       |      ...      |
|  |                  |    |                  |               |
|  | +--------------+ |    | +--------------+ |               |
|  | | Warp 0       | |    | | Warp 0       | |               |
|  | | th0..th31    | |    | | th0..th31    | |               |
|  | +--------------+ |    | +--------------+ |               |
|  | +--------------+ |    | +--------------+ |               |
|  | | Warp 1       | |    | | Warp 1       | |               |
|  | | th32..th63   | |    | | th32..th63   | |               |
|  | +--------------+ |    | +--------------+ |               |
|  +------------------+    +------------------+               |
+-------------------------------------------------------------+

Warp-Based SIMD

A warp executes one instruction across many threads. When all lanes are active and follow the same path, the warp behaves like a SIMD vector operation.

Efficient Warp Execution

int i = blockIdx.x * blockDim.x + threadIdx.x;
C[i] = A[i] + B[i];

Every thread executes the same instruction sequence.

Divergent Warp Execution

if (threadIdx.x % 2 == 0) {
  C[i] = A[i] + B[i];
} else {
  C[i] = A[i] - B[i];
}

Threads in the same warp take different branches. The hardware may execute both paths while masking off inactive lanes, reducing efficiency.

Traditional SIMD vs. Warp-Based SIMT

Program Structure in CUDA

A CUDA program usually has two parts:

1. Host code running on the CPU.
2. Device code running on the GPU.

A Simple Kernel

A CUDA kernel is marked with __global__. It runs on the GPU and is launched from host code.

__global__ void scale(float* d, float s, int n) {
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i < n) {
    d[i] *= s;
  }
}

Allocate Device Memory

float* d_data = nullptr;
size_t bytes = n * sizeof(float);

cudaError_t err = cudaMalloc(&d_data, bytes);
if (err != cudaSuccess) {
  // Handle allocation error
}

The CPU variable d_data is a device pointer that refers to memory allocated in the GPU's address space.

Copy Data to the GPU

cudaMemcpy(d_data, h_data, bytes, cudaMemcpyHostToDevice);

The copy direction must be specified explicitly.

Launch the Kernel

int block = 256;
int grid  = (n + block - 1) / block;

scale<<<grid, block>>>(d_data, 2.0f, n);
cudaDeviceSynchronize();

The grid size is rounded up so every element is covered. The bounds check inside the kernel prevents out-of-range memory access.

Copy Results Back and Free Memory

cudaMemcpy(h_data, d_data, bytes, cudaMemcpyDeviceToHost);
cudaFree(d_data);

Shared Memory Example: Block Reduction

Shared memory lets threads in a block cooperate. The example below computes one partial sum per block.

__global__ void blockSum(const float* in, float* out, int n) {
  extern __shared__ float buf[];

  int t = threadIdx.x;
  int i = blockIdx.x * blockDim.x + threadIdx.x;

  buf[t] = (i < n) ? in[i] : 0.0f;
  __syncthreads();

  for (int stride = blockDim.x / 2; stride > 0; stride >>= 1) {
    if (t < stride) {
      buf[t] += buf[t + stride];
    }
    __syncthreads();
  }

  if (t == 0) {
    out[blockIdx.x] = buf[0];
  }
}

Launch example:

int block = 256;
int grid = (n + block - 1) / block;
size_t sharedBytes = block * sizeof(float);

blockSum<<<grid, block, sharedBytes>>>(d_in, d_partial, n);

This version addresses two common issues:

• It checks i < n before reading input.
• It uses dynamically sized shared memory so the buffer matches the block size.

The halving reduction shown here assumes a power-of-two block size, such as the 256 threads used in the launch example.

Common CUDA Concepts

Performance Checklist

Use this checklist when improving a GPU kernel:

Common Mistakes