Programming Models & Tools

Complete revision notes. Everything in this unit, organised by topic — definitions, diagrams, formulas, comparison tables and worked examples. Revise from this alone and you have the whole unit.

Part A — OpenMP in Depth

A1. The model and the boilerplate

OpenMP is an API for shared-memory parallelism in C, C++ and Fortran, consisting of compiler directives, a runtime library and environment variables. It uses the fork-join model: the master thread forks a team at a parallel region, and the threads join at an implicit barrier at the end.

#include <omp.h>

#pragma omp parallel              // fork a team; all threads run this block
{
    int id = omp_get_thread_num();      // 0 .. n-1
    int n  = omp_get_num_threads();     // team size
}                                 // implicit barrier — threads join here

omp_set_num_threads(4);           // request a team size
double t = omp_get_wtime();       // wall-clock timer for measurement

compile:  gcc -fopenmp -O2 prog.c -o prog
control:  export OMP_NUM_THREADS=8

A key property: compile without -fopenmp and the pragmas are ignored, leaving a valid serial program. This is why OpenMP is described as incremental parallelisation — you annotate existing code rather than rewriting it.

A2. Work-sharing constructs

A parallel region alone makes every thread execute the same code. Work-sharing constructs divide the work between them.

// The combined shortcut you will use most: parallel + for in one line
#pragma omp parallel for
for (int i = 0; i < N; i++)
    a[i] = b[i] + c[i];        // each thread gets a slice of the iterations

// Task parallelism with sections
#pragma omp parallel sections
{
    #pragma omp section
    load_file();
    #pragma omp section
    compress();
}

A3. Scheduling — who gets which iterations

schedule(kind [, chunk]) decides how loop iterations map onto threads. This is Unit III's load balancing, made concrete.

Worked example — one loop, four schedules. 12 iterations, 4 threads; iteration i costs (i+1) ms, so 1, 2, 3 … 12 ms, total 78 ms.

A4. Data scoping — the number one bug source

SHARED BY ACCIDENT — the classic bug

int i, j, temp;                       int i, j, temp;
#pragma omp parallel for              #pragma omp parallel for private(j, temp)
for (i = 0; i < N; i++)               for (i = 0; i < N; i++)
  for (j = 0; j < N; j++) {             for (j = 0; j < N; j++) {
    temp = a[i][j] * 2;                   temp = a[i][j] * 2;
    b[i][j] = temp;                       b[i][j] = temp;
  }                                     }

  i is private (it's the loop           j and temp are private now.
  counter). j and temp are SHARED —     Correct, and deterministic.
  threads overwrite each other.
  Garbage output, different each run.

A5. Synchronisation, and the barrier you did not write

Worked example — the barrier you did not write. Every omp for, sections and single ends with an implicit barrier you never typed.

#pragma omp parallel
{
    #pragma omp for
    for (i=0;i<N;i++)  a[i] = f(i);      // <-- invisible barrier here
    #pragma omp for
    for (i=0;i<N;i++)  b[i] = g(i);      // <-- and here
}

Each loop averages 100 ms per thread, but the slowest thread takes 120 ms.
Four such loops:

  With the implicit barriers  -> every thread waits for the slowest,
                                 four times            = 480 ms
  With nowait on the first 3  -> waits only once, at
                                 the end of the region = ~420 ms

#pragma omp for nowait
for (i=0;i<N;i++)  a[i] = f(i);
#pragma omp for
for (i=0;i<N;i++)  b[i] = a[i] * 2;   // may read a[i] before it was written

A6. Pitfalls and good habits

Part B — MPI: Message Passing

B1. The MPI mental model

MPI (Message Passing Interface) is a standard library for distributed-memory parallel programming. It runs processes, not threads: mpirun -np 4 launches four full copies of your program, each with its own private memory. They cannot see each other's variables, ever. The only way to share data is to send a message.

MPI programs are SPMD (Unit I): one source file, run by every process, each behaving differently based on its rank. if (rank == 0) { ... } is how one process plays "master".

B2. The four calls every MPI program has

#include <mpi.h>
int main(int argc, char **argv) {
    MPI_Init(&argc, &argv);                    // 1. start MPI (first call)
    int rank, size;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);      // 2. which process am I?
    MPI_Comm_size(MPI_COMM_WORLD, &size);      // 3. how many of us?
    printf("Hello from rank %d of %d\n", rank, size);
    MPI_Finalize();                            // 4. shut MPI down (last call)
    return 0;
}

compile:  mpicc hello.c -o hello
run:      mpirun -np 4 ./hello

Output (order varies between runs!):
    Hello from rank 2 of 4
    Hello from rank 0 of 4
    Hello from rank 3 of 4
    Hello from rank 1 of 4

The print order changes between runs because the processes genuinely run at once. That non-determinism is your first sign of real parallelism — the same lesson as OpenMP thread ordering in Unit I.

B3. Point-to-point communication

MPI_Send(&value, count, MPI_INT, dest,   tag, MPI_COMM_WORLD);
MPI_Recv(&value, count, MPI_INT, source, tag, MPI_COMM_WORLD, &status);
//        buffer  how-many  type   who    label  group

B4. The deadlock everyone writes once

int partner = 1 - rank;                    // 0 <-> 1
MPI_Send(sendbuf, n, MPI_DOUBLE, partner, 0, MPI_COMM_WORLD);
MPI_Recv(recvbuf, n, MPI_DOUBLE, partner, 0, MPI_COMM_WORLD, &status);

Neither rank reaches its MPI_Recv, because each is stuck inside MPI_Send waiting for the other to reach a MPI_Recv it will never reach. Both processes sit at 100% CPU forever with no error message.

Three fixes, worst to best:

1. ORDER THEM BY RANK — somebody sends before everybody receives
   if (rank % 2 == 0) { MPI_Send(...); MPI_Recv(...); }
   else               { MPI_Recv(...); MPI_Send(...); }

2. USE MPI_Sendrecv — one call, the library handles the ordering
   MPI_Sendrecv(sendbuf, n, MPI_DOUBLE, partner, 0,
                recvbuf, n, MPI_DOUBLE, partner, 0,
                MPI_COMM_WORLD, &status);

3. USE NON-BLOCKING — post both, then wait; also overlaps comms with work
   MPI_Irecv(recvbuf, n, MPI_DOUBLE, partner, 0, comm, &req[0]);
   MPI_Isend(sendbuf, n, MPI_DOUBLE, partner, 0, comm, &req[1]);
   /* do useful computation here */
   MPI_Waitall(2, req, MPI_STATUSES_IGNORE);

B5. Collectives — the workhorses

You could build everything from Send and Recv, but you should not. Collectives move data among all processes at once, and the library implements them far more efficiently — typically in log(P) steps rather than P. Every process in the communicator must call the collective.

   BCAST                SCATTER              GATHER / REDUCE
   root: [X]            root: [A B C D]      ranks: [a][b][c][d]
      |  |  |  |              |  |  |  |          \  |  |  /
   [X][X][X][X]          [A][B][C][D]              [a+b+c+d]
   all get the same     each gets a chunk       combined onto root

   PARALLEL SORT (Lab 2) uses exactly this shape:
     Scatter chunks -> each rank sorts locally -> Gather -> final merge

Scatter is fan-out and Reduce is fan-in. Together they are the map-reduce pattern — the same shape as the large-scale data jobs in Unit V.

B6. OpenMP vs MPI

Hybrid MPI + OpenMP is the real supercomputer recipe: MPI between nodes, OpenMP within each node, and CUDA on the GPU inside each node. Messages across the cluster, shared memory inside each box.

Part C — CUDA and the GPU Stack

C1. Host and device

CUDA is C with a few extensions that let the CPU hand work to the GPU. Two worlds, connected by a bus:

• Host — the CPU. Runs main(), owns normal RAM, orchestrates: allocates GPU memory, copies data, launches kernels.
• Device — the GPU. Runs kernels across thousands of threads at once, owns its own separate memory.

A KERNEL is a function marked __global__ that runs on the GPU.

__global__ void saxpy(int n, float a, float *x, float *y) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) y[i] = a * x[i] + y[i];      // bounds check is essential
}

saxpy<<< blocks, threadsPerBlock >>>(n, 2.0f, d_x, d_y);
//     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^  the launch configuration
//     "run this kernel on (blocks x threadsPerBlock) threads"

Function qualifiers:  __global__  called from host, runs on device
                      __device__  called from device, runs on device
                      __host__    ordinary CPU function (the default)

C2. Thread → block → grid

GRID  (all the threads a kernel launches)
 +-----------------------------------------------+
 |  BLOCK 0            BLOCK 1            BLOCK 2 |
 | [t0 t1 t2 t3]      [t0 t1 t2 t3]      [t0 ...] |   threadIdx.x = 0..3
 +-----------------------------------------------+
   blockIdx.x=0        blockIdx.x=1       blockIdx.x=2    blockDim.x = 4

 GLOBAL INDEX OF A THREAD — memorise this line:

   i = blockIdx.x * blockDim.x + threadIdx.x

   e.g. block 2, thread 1, blockDim 4  ->  2 * 4 + 1 = 9

 LAUNCH CONFIGURATION for n elements:
   int threads = 256;
   int blocks  = (n + threads - 1) / threads;   // round UP
   kernel<<<blocks, threads>>>(n, ...);
   // rounding up means the last block has spare threads —
   // hence the  if (i < n)  bounds check inside the kernel

Threads within a block are further grouped by the hardware into warps of 32 (Unit II), which is why block sizes are conventionally multiples of 32 — 128 or 256 are typical. Threads in the same block can cooperate via shared memory and __syncthreads(); threads in different blocks cannot.

C3. The five-step pattern

The GPU cannot read your ordinary C arrays. Host memory and device memory are separate, and moving data across the PCIe bus is the real cost of GPU programming — so you do it as little as possible.

EVERY CUDA PROGRAM FOLLOWS THESE FIVE STEPS

1. cudaMalloc(&d_x, bytes);                    // allocate on the DEVICE
2. cudaMemcpy(d_x, h_x, bytes, H2D);           // copy inputs  HOST -> DEVICE
3. kernel<<<blocks, threads>>>(d_x, ...);      // launch
4. cudaMemcpy(h_y, d_y, bytes, D2H);           // copy result  DEVICE -> HOST
5. cudaFree(d_x);                              // free device memory

H2D = cudaMemcpyHostToDevice,  D2H = cudaMemcpyDeviceToHost

Kernel launches are ASYNCHRONOUS — control returns to the host
immediately. cudaDeviceSynchronize() waits for the GPU to finish.

Learn these five steps once and every kernel is just step 3 getting cleverer. Steps 2 and 4 are also why a GPU can lose to a CPU on small jobs (Unit II) — the transfer may cost more than the computation saves.

C4. GPU memory spaces

Getting speed out of a GPU is mostly about which memory you use: closer means faster and smaller. The standard trick is to load data once from slow global memory into fast shared memory, do all the back-and-forth work there, and write back once.

C5. Coalescing — worth 30×

A warp's 32 threads issue their memory request at the same instant. If those 32 addresses are consecutive, the hardware satisfies them in one transaction. If they are scattered, it needs 32.

COALESCED                          STRIDED
c[i] = a[i] * 2;                   c[i] = a[i*32] * 2;
warp reads a[0..31]                warp reads a[0], a[32], a[64] ... a[992]
= 128 consecutive bytes            = 32 separate transactions
= ONE transaction                  = 3% of each one actually used
100% of it used

   ~2,800 GB/s                        ~90 GB/s        -> 31x slower

Same instruction count. Same number of elements. 31x the runtime.

// thread per ROW — feels natural,      // thread per COLUMN — feels wrong,
// SLOW on a GPU                        // FAST on a GPU
for (j=0;j<N;j++)                       for (j=0;j<N;j++)
    sum += a[tid*N + j];                    sum += a[j*N + tid];

At any instant thread 0 reads a[0],     At any instant threads 0-31 read
thread 1 reads a[N], thread 2 reads     a[0..31] — consecutive.
a[2N] — scattered by N.                 ONE transaction.

C6. Parallel reduction on a GPU

Summing an array on a GPU is the canonical non-trivial kernel, because it needs threads to cooperate rather than work independently. The pattern is a tree — the same tree whose work and span you computed in Unit III.

TREE REDUCTION INSIDE ONE BLOCK (8 threads shown)

  data   [ 3  1  7  0  4  1  6  3 ]
                                        stride = 4
  step1  [ 7  2 13  3  |  .  .  .  ]    t0..t3 add data[i+4]
                                        stride = 2
  step2  [20  5  |  .  .  .  .  .  ]    t0,t1 add data[i+2]
                                        stride = 1
  step3  [25  |  .  .  .  .  .  .  ]    t0 adds data[i+1]

  log2(8) = 3 steps instead of 7 sequential adds.
  __syncthreads() between every step — or threads read stale values.

__global__ void reduce(float *in, float *out) {
    __shared__ float s[256];
    int tid = threadIdx.x;
    int i   = blockIdx.x * blockDim.x + threadIdx.x;
    s[tid] = in[i];                       // global -> shared, once
    __syncthreads();
    for (int stride = blockDim.x/2; stride > 0; stride /= 2) {
        if (tid < stride) s[tid] += s[tid + stride];
        __syncthreads();                  // ESSENTIAL
    }
    if (tid == 0) out[blockIdx.x] = s[0]; // one partial per block
}

Each block produces one partial sum; a second pass (or an atomic add) combines the partials. Note the structure: accumulate locally, combine once — the identical principle as an OpenMP reduction (Unit I), the fix for false sharing (Unit II) and the combiner in MapReduce (Unit V).

Quick recall

OPENMP
  #pragma omp parallel for reduction(+:sum) schedule(dynamic) private(j)
  omp_get_thread_num()  omp_get_num_threads()  omp_get_wtime()
  gcc -fopenmp -O2 prog.c

MPI
  MPI_Init / MPI_Comm_rank / MPI_Comm_size / MPI_Finalize
  MPI_Send(buf, count, type, dest, tag, comm)
  MPI_Recv(buf, count, type, source, tag, comm, &status)
  Bcast · Scatter · Gather · Reduce · Allreduce · Barrier
  mpicc prog.c -o prog   &&   mpirun -np 4 ./prog

CUDA
  __global__ void kernel(...)          kernel<<<blocks, threads>>>(...)
  i = blockIdx.x * blockDim.x + threadIdx.x        <-- the one to memorise
  blocks = (n + threads - 1) / threads             <-- round up
  cudaMalloc · cudaMemcpy(H2D) · launch · cudaMemcpy(D2H) · cudaFree
  __shared__   __syncthreads()   warp = 32 threads

Self-check

Answer each out loud before opening it.