LUNAROPS · OPERATIONAL UPLINK 100% UPTIME 1,247d POSTS 893 JEFF.MOON@LUNAROPS.DEV UTC --:--:--

MPI Programming Essentials: Collectives, Non-Blocking, and the Traps Nobody Warns You About

mpihpcparallel-computingopenmpimpichscientific-computingslurmcfortranpython

If you want to run a program across more than one machine in an HPC cluster, you are almost certainly going to write MPI. Three decades after its introduction, the Message Passing Interface is still the dominant model for distributed-memory parallelism — it runs every major weather model, every serious CFD code, every computational chemistry package, and the communication layer underneath most large-scale ML frameworks.

MPI is not hard to learn. MPI is hard to learn well. The API surface is enormous (600+ functions), the most common tutorials teach patterns that perform terribly at scale, and the error messages when you mess up are historically unhelpful. This post is the working-engineer companion: enough theory to understand what the runtime is doing, the ~20 functions you’ll actually use, and the failure modes nobody warns you about until you’ve spent a week debugging them.


The Mental Model

An MPI program is a single source file (usually) that is launched as N identical processes. Every process runs the same executable from main to exit. They are not threads — each has its own address space, possibly on a different physical machine. They coordinate by explicitly sending messages.

Three concepts sit at the center:

Rank. Every process gets an integer ID from 0 to N-1, called its rank, inside a given communicator. The rank is how you identify who you are sending to or receiving from.

Size. The total number of processes in a communicator.

Communicator. A named group of processes that can talk to each other. Every program gets MPI_COMM_WORLD by default, containing all launched processes. You can create smaller communicators by splitting or grouping.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
#include <mpi.h>
#include <stdio.h>

int main(int argc, char **argv) {
    MPI_Init(&argc, &argv);

    int rank, size;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);

    printf("Hello from rank %d of %d\n", rank, size);

    MPI_Finalize();
    return 0;
}

Compile and run:

1
2
3
4
5
6
mpicc hello.c -o hello
mpirun -np 4 ./hello
# Hello from rank 0 of 4
# Hello from rank 1 of 4
# Hello from rank 3 of 4
# Hello from rank 2 of 4    (order is not deterministic)

That’s the whole model. Everything else is variations on “rank A sends bytes to rank B, possibly via collective operations that coordinate many ranks at once.”


Picking an Implementation

MPI is a specification; you need a concrete implementation. The three that matter in 2026:

Implementation Strengths Typical use
OpenMPI Permissive license, wide platform support, flexible runtime Most university clusters, general HPC
MPICH Reference implementation, rock-solid correctness, base for vendor MPIs Government labs, embedded in Cray/Intel/HPE stacks
Intel MPI MPICH-derived, highly optimized for Intel fabrics and CPUs Commercial HPC on Intel hardware

Other notable derivatives: MVAPICH2 (MPICH-based, tuned for InfiniBand), HPE Cray MPI (closed, for Slingshot fabrics), NVIDIA HPC-X (OpenMPI + UCX tuned for NVIDIA NICs).

Practical advice: on a cluster, load whichever your site has installed via module load openmpi/4.1 or similar. Never mix implementations — a binary compiled against OpenMPI cannot run with MPICH’s mpirun. If you’re building your own code, use mpicc/mpicxx/mpifort (wrapper compilers) — they invoke your system C/C++/Fortran compiler with the right include and library paths.


Point-to-Point Communication

The two functions that started it all:

1
2
3
4
5
int MPI_Send(const void *buf, int count, MPI_Datatype dt,
             int dest, int tag, MPI_Comm comm);

int MPI_Recv(void *buf, int count, MPI_Datatype dt,
             int source, int tag, MPI_Comm comm, MPI_Status *status);

A trivial exchange: rank 0 sends an array of doubles to rank 1:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
int rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);

double data[100];
if (rank == 0) {
    for (int i = 0; i < 100; i++) data[i] = i * 3.14;
    MPI_Send(data, 100, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
} else if (rank == 1) {
    MPI_Recv(data, 100, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
    printf("Rank 1 got data[50] = %f\n", data[50]);
}

Things worth noticing:

  • The datatype is named, not inferred. MPI_DOUBLE, MPI_INT, MPI_CHAR, etc. For user structs, you build an MPI_Datatype (see later).
  • Tags are integer labels letting you distinguish different messages between the same pair of ranks. MPI_ANY_TAG matches any tag; MPI_ANY_SOURCE matches any sender.
  • MPI_Send is blocking but not synchronous. It returns when the buffer is safe to reuse — which might be because the message was copied to an internal buffer, not because rank 1 actually received it. For guaranteed synchronous behavior, use MPI_Ssend.

The Eager/Rendezvous Threshold

This is the first non-obvious thing every MPI programmer trips on. Internally, MPI_Send has two modes based on message size:

  • Eager (small messages): the data is sent immediately to a pre-allocated buffer on the receiver; MPI_Send returns as soon as the data leaves the sender.
  • Rendezvous (large messages): the sender first sends a small “ready to send N bytes” handshake, the receiver acknowledges once it has posted a matching MPI_Recv, then the data transfers directly.

The threshold defaults to somewhere around 64 KB and is configurable. Symptom of hitting this line unexpectedly: a program that runs fine on 10 KB messages deadlocks on 100 KB messages because every rank is stuck in MPI_Send waiting for its peer to post MPI_Recv, but the peer is also stuck in MPI_Send. The fix is either to reorder your sends/receives or to use non-blocking primitives.

The Classic Deadlock

1
2
3
4
5
6
7
8
// WRONG — deadlocks for large messages
if (rank == 0) {
    MPI_Send(buf, N, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD);
    MPI_Recv(buf, N, MPI_DOUBLE, 1, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
} else if (rank == 1) {
    MPI_Send(buf, N, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);   // Both ranks wait here
    MPI_Recv(buf, N, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
}

Three common fixes: swap the order on one side (ugly), use MPI_Sendrecv (clean), or switch to non-blocking (general):

1
2
3
4
// Using MPI_Sendrecv — symmetric, deadlock-free
MPI_Sendrecv(sendbuf, N, MPI_DOUBLE, peer, 0,
             recvbuf, N, MPI_DOUBLE, peer, 0,
             MPI_COMM_WORLD, MPI_STATUS_IGNORE);

Collective Operations

Collectives are the real workhorses. They express communication patterns involving all ranks in a communicator at once, and MPI implementations spend enormous engineering effort making them fast. If your algorithm can be expressed with collectives, use them — never hand-roll with point-to-point.

Every collective requires that every rank in the communicator call it with matching arguments. Forget one rank and the whole program hangs.

The Essential Eight

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
// Broadcast: rank 'root' sends buf to everyone
MPI_Bcast(buf, count, dt, root, comm);

// Reduce: combine values from all ranks into root
MPI_Reduce(sendbuf, recvbuf, count, dt, MPI_SUM, root, comm);

// Allreduce: like reduce, but every rank gets the result
MPI_Allreduce(sendbuf, recvbuf, count, dt, MPI_SUM, comm);

// Scatter: root splits sendbuf into N equal pieces, one per rank
MPI_Scatter(sendbuf, n, dt, recvbuf, n, dt, root, comm);

// Gather: inverse of scatter; each rank contributes n elements to root
MPI_Gather(sendbuf, n, dt, recvbuf, n, dt, root, comm);

// Allgather: gather, but every rank gets the full result
MPI_Allgather(sendbuf, n, dt, recvbuf, n, dt, comm);

// Alltoall: each rank sends a distinct piece to every other rank
MPI_Alltoall(sendbuf, n, dt, recvbuf, n, dt, comm);

// Barrier: synchronize all ranks (rarely needed — collectives synchronize)
MPI_Barrier(comm);

The most useful reduction operations: MPI_SUM, MPI_PROD, MPI_MIN, MPI_MAX, MPI_LAND/MPI_LOR (logical), MPI_BAND/MPI_BOR (bitwise), MPI_MINLOC/MPI_MAXLOC (which rank had the min/max). You can also define custom reductions with MPI_Op_create.

Worked Example: Parallel Sum

1
2
3
4
5
6
// Every rank computes a local sum of its chunk; we want the global sum.
double local_sum = 0;
for (int i = start; i < end; i++) local_sum += array[i];

double global_sum;
MPI_Allreduce(&local_sum, &global_sum, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);

After this call, every rank holds the same global_sum. Internally, the implementation uses a recursive doubling or ring algorithm — O(log N) or O(N) rounds depending on message size — that is dramatically faster than the naïve “everyone sends to rank 0, rank 0 sums, rank 0 broadcasts” approach most beginners write.

The Vector Variants

When ranks contribute different amounts of data, the v-suffix collectives accept per-rank counts and displacements:

1
2
3
4
int counts[size], displs[size];
// ... compute per-rank counts and displacements ...
MPI_Gatherv(sendbuf, my_count, MPI_DOUBLE,
            recvbuf, counts, displs, MPI_DOUBLE, 0, MPI_COMM_WORLD);

MPI_Scatterv, MPI_Gatherv, MPI_Allgatherv, MPI_Alltoallv exist for the same reason. Any time your data is unevenly distributed, reach for these.


Non-Blocking Communication

Every MPI_X has a non-blocking counterpart MPI_Ix. They return an MPI_Request handle immediately and let computation proceed while the network does its work. You call MPI_Wait or MPI_Test later to confirm completion.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
MPI_Request request;
MPI_Status status;

// Post the receive first (best practice)
MPI_Irecv(recvbuf, N, MPI_DOUBLE, peer, 0, MPI_COMM_WORLD, &request);

// Do useful work while data is in flight
compute_something(local_data);

// Now block until the receive completes
MPI_Wait(&request, &status);

This is the single biggest lever for MPI performance at scale. Any communication step that can be overlapped with independent computation should be non-blocking.

Multi-Request Wait

For a halo exchange (common in stencil codes), you typically post many sends and receives:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
MPI_Request reqs[8];
MPI_Irecv(north_buf, N, MPI_DOUBLE, north_rank, 0, comm, &reqs[0]);
MPI_Irecv(south_buf, N, MPI_DOUBLE, south_rank, 0, comm, &reqs[1]);
MPI_Irecv(east_buf,  N, MPI_DOUBLE, east_rank,  0, comm, &reqs[2]);
MPI_Irecv(west_buf,  N, MPI_DOUBLE, west_rank,  0, comm, &reqs[3]);
MPI_Isend(my_north, N, MPI_DOUBLE, north_rank, 0, comm, &reqs[4]);
MPI_Isend(my_south, N, MPI_DOUBLE, south_rank, 0, comm, &reqs[5]);
MPI_Isend(my_east,  N, MPI_DOUBLE, east_rank,  0, comm, &reqs[6]);
MPI_Isend(my_west,  N, MPI_DOUBLE, west_rank,  0, comm, &reqs[7]);

// Update the interior (doesn't need halo data) while network runs
update_interior(grid);

// Wait for all 8 to finish, then update boundary
MPI_Waitall(8, reqs, MPI_STATUSES_IGNORE);
update_boundary(grid);

MPI_Waitany, MPI_Waitsome, MPI_Testany, MPI_Testsome round out the family for “as soon as one completes” patterns.

Non-Blocking Collectives (MPI-3)

MPI-3 added non-blocking collectives: MPI_Ibcast, MPI_Iallreduce, MPI_Ibarrier, etc. These are genuinely useful — you can overlap a global all-reduce with independent computation. Few codes use them aggressively because the semantics are subtle (you cannot modify the input buffer until MPI_Wait returns), but the payoff at scale is substantial.

1
2
3
4
5
6
7
8
9
MPI_Request req;
MPI_Iallreduce(&local_norm, &global_norm, 1, MPI_DOUBLE, MPI_SUM,
               MPI_COMM_WORLD, &req);

// Start the next iteration's computation while the all-reduce flies
prepare_next_iteration(state);

MPI_Wait(&req, MPI_STATUS_IGNORE);
// global_norm is now valid; use it for convergence check

Derived Datatypes

If your data is not a contiguous array of primitives — a column of a matrix stored row-major, a struct with padding, a halo boundary — you have three options:

  1. Copy into a contiguous buffer before sending (simple, extra memory traffic).
  2. Make N small sends (correct, awful performance).
  3. Build a derived datatype that describes the in-memory layout, then send once.

Option 3 is the right answer for performance.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
// Send every 10th element of a large array: stride of 10 doubles
MPI_Datatype stride_type;
MPI_Type_vector(100,      // count of blocks
                1,        // elements per block
                10,       // stride between blocks (elements)
                MPI_DOUBLE,
                &stride_type);
MPI_Type_commit(&stride_type);

MPI_Send(array, 1, stride_type, dest, 0, MPI_COMM_WORLD);

MPI_Type_free(&stride_type);

Useful constructors:

Function Describes
MPI_Type_contiguous N contiguous copies of an existing type
MPI_Type_vector Regularly spaced blocks (constant stride)
MPI_Type_indexed Irregularly spaced blocks
MPI_Type_create_struct Heterogeneous fields, like a C struct
MPI_Type_create_subarray A subarray of a multi-dimensional array

Every derived type must be committed before use and freed when done. Failing to free leaks type handles.

Struct Example

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
typedef struct {
    int id;
    double x, y, z;
    float mass;
} Particle;

MPI_Datatype particle_type;
int blocklens[]    = {1, 3, 1};
MPI_Aint displs[]  = {offsetof(Particle, id),
                      offsetof(Particle, x),
                      offsetof(Particle, mass)};
MPI_Datatype types[] = {MPI_INT, MPI_DOUBLE, MPI_FLOAT};

MPI_Type_create_struct(3, blocklens, displs, types, &particle_type);
MPI_Type_commit(&particle_type);

// Now MPI treats Particle as a first-class type
MPI_Bcast(particles, 1000, particle_type, 0, MPI_COMM_WORLD);

This is the right pattern for sending arrays of structs — it respects alignment and padding, which hand-serialization frequently gets wrong.


Communicators, Groups, and Cartesian Topologies

MPI_COMM_WORLD is the default, but real codes benefit enormously from creating sub-communicators matched to the algorithm’s structure.

Splitting a Communicator

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
// Split into "even" and "odd" ranks
int color = rank % 2;
MPI_Comm sub_comm;
MPI_Comm_split(MPI_COMM_WORLD, color, rank, &sub_comm);

// Now sub_comm contains only ranks with the same color
int sub_rank, sub_size;
MPI_Comm_rank(sub_comm, &sub_rank);
MPI_Comm_size(sub_comm, &sub_size);

// ... work within sub_comm ...

MPI_Comm_free(&sub_comm);

Use cases: splitting a cluster into row/column groups for 2D matrix multiplication, separating I/O ranks from compute ranks, isolating subdomains in a multi-physics simulation.

Cartesian Topologies

For structured-grid codes, let MPI help you compute neighbors:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
int dims[2] = {0, 0};      // 0 = let MPI choose
MPI_Dims_create(size, 2, dims);

int periods[2] = {1, 1};   // Periodic boundaries on both dims
MPI_Comm cart_comm;
MPI_Cart_create(MPI_COMM_WORLD, 2, dims, periods, 1, &cart_comm);

int my_coords[2];
MPI_Cart_coords(cart_comm, rank, 2, my_coords);

int north, south, east, west;
MPI_Cart_shift(cart_comm, 0, 1, &north, &south);
MPI_Cart_shift(cart_comm, 1, 1, &west, &east);

// north/south/east/west now hold rank IDs for halo exchange
// (or MPI_PROC_NULL for non-periodic boundaries)

MPI_Cart_shift handles the periodic wrap-around and returns MPI_PROC_NULL for ranks at a non-periodic edge — sending to MPI_PROC_NULL is a no-op, which elegantly avoids special cases at boundaries.


One-Sided Communication (RMA)

MPI-2 added Remote Memory Access (RMA), where a process can read or write another process’s memory without explicit cooperation. This maps directly to hardware RDMA on modern fabrics and is often the fastest option for irregular communication patterns.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
double *data;
MPI_Win window;

// Allocate a window of shared memory
MPI_Win_allocate(N * sizeof(double), sizeof(double),
                 MPI_INFO_NULL, MPI_COMM_WORLD, &data, &window);

// Initialize local data
for (int i = 0; i < N; i++) data[i] = rank * N + i;

// Synchronize, then start the access epoch
MPI_Win_fence(0, window);

// Read 10 doubles from rank 0 into local buffer
double buf[10];
if (rank != 0) {
    MPI_Get(buf, 10, MPI_DOUBLE, 0, 5, 10, MPI_DOUBLE, window);
}

MPI_Win_fence(0, window);
// buf now contains rank 0's data[5..14]

MPI_Win_free(&window);

Three synchronization styles:

  • Fence (MPI_Win_fence) — simple, collective, barrier-like.
  • Post/start/complete/wait — fine-grained, scalable, complex.
  • Lock/unlock (passive target) — truly one-sided, no coordination from the target process.

RMA is powerful but tricky — memory ordering and consistency rules are stricter than most programmers expect. Use it when you’ve profiled and know point-to-point is the bottleneck.


MPI + Threading

By default, MPI makes no promises about thread safety. You must initialize with a specific thread level:

1
2
3
4
5
6
int provided;
MPI_Init_thread(&argc, &argv, MPI_THREAD_MULTIPLE, &provided);
if (provided < MPI_THREAD_MULTIPLE) {
    fprintf(stderr, "MPI does not support MPI_THREAD_MULTIPLE\n");
    MPI_Abort(MPI_COMM_WORLD, 1);
}

Four levels:

Level Meaning
MPI_THREAD_SINGLE Only one thread per process, period
MPI_THREAD_FUNNELED Multiple threads, only main thread calls MPI
MPI_THREAD_SERIALIZED Multiple threads call MPI, but only one at a time
MPI_THREAD_MULTIPLE Any thread can call MPI concurrently

MPI_THREAD_FUNNELED is the sweet spot for hybrid MPI + OpenMP codes — OpenMP parallel regions do their thing, and MPI calls happen outside them from the master thread. Requesting MPI_THREAD_MULTIPLE when you don’t need it costs performance because the implementation must serialize internally.

1
2
3
4
5
6
7
8
9
#pragma omp parallel
{
    // Many threads working on local data
    for (int i = omp_get_thread_num(); i < N; i += omp_get_num_threads()) {
        compute(data[i]);
    }
}
// Back in serial region — main thread calls MPI
MPI_Allreduce(...);

MPI-IO

Parallel I/O is not a solved problem, but MPI-IO is the most portable solution. It lets every rank write to a shared file with guaranteed consistency:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
MPI_File fh;
MPI_File_open(MPI_COMM_WORLD, "results.dat",
              MPI_MODE_CREATE | MPI_MODE_WRONLY,
              MPI_INFO_NULL, &fh);

// Each rank writes N doubles at offset rank*N*sizeof(double)
MPI_Offset offset = rank * N * sizeof(double);
MPI_File_write_at_all(fh, offset, local_data, N, MPI_DOUBLE,
                      MPI_STATUS_IGNORE);

MPI_File_close(&fh);

The _all suffix means the call is collective — the implementation can coordinate writes for drastically better performance on parallel filesystems (Lustre, GPFS, BeeGFS). Non-collective variants exist but perform worse at scale.

For structured scientific data, prefer higher-level libraries built on MPI-IO:

  • HDF5 parallel — the de facto standard in scientific computing
  • NetCDF-4 — HDF5-backed, widely used in climate/ocean models
  • ADIOS2 — newer, designed for exascale with staging/compression

Performance Traps

1. Too many small messages. Latency dominates below ~10 KB. Batch small updates into fewer larger messages whenever the algorithm allows it.

2. Load imbalance. If one rank takes 2× as long as the others, a collective’s time is dominated by that rank. Profile with MPI_Barrier bracketing to find imbalance.

3. The mythical MPI_Barrier. Most reasons to call MPI_Barrier are wrong. Collectives already synchronize their participants. Barrier is almost never needed for correctness; it’s occasionally useful for timing or debugging.

4. Oversubscription. Launching 64 MPI ranks on a 16-core node makes everything fight for CPUs and caches. Rule of thumb: 1 rank per physical core, or 1 rank per NUMA node for hybrid MPI + OpenMP codes.

5. Process binding matters. Without --bind-to core (OpenMPI) or equivalent, the OS scheduler can migrate ranks between cores, destroying cache locality. Always pin ranks.

6. Naïve broadcasts and reductions. Code like if (rank == 0) for (i=1; i<size; i++) MPI_Send(buf, ..., i, ...); is O(N). MPI_Bcast is O(log N). The difference at 10,000 ranks is 14 rounds versus 10,000.

7. Not using MPI_IN_PLACE for all-reduce. MPI_Allreduce(MPI_IN_PLACE, buf, ...) saves a buffer copy versus separate send/recv buffers. Use it when the input can be overwritten.

8. Blocking where non-blocking would work. Any step with ≥1 ms of computation that’s independent of in-flight data should overlap with MPI_Isend/MPI_Irecv.

9. Synchronous-mode MPI_Ssend by accident. Some users reach for MPI_Ssend thinking “more synchronous is safer.” It’s strictly slower for almost every real use case. Use MPI_Send and think about deadlock explicitly.

10. Not scaling communicators. MPI_COMM_WORLD collectives across 100,000 ranks are slower than the same operation on a tree of sub-communicators. For truly massive jobs, split.


Debugging MPI

MPI debuggers are their own world because you need to attach to N processes on M machines. Options:

  • printf with rank prefix. Always the first tool. printf("[rank %d] ...\n", rank); fflush(stdout); — fflush is critical; stdout buffering lies to you.
  • mpirun -np N xterm -e gdb ./program — launches a gdb per rank in its own terminal. Primitive but works over X forwarding.
  • TotalView, DDT (ARM Forge), Linaro MAP — commercial, powerful, expensive. Standard at most large centers.
  • Valgrind memcheck — works, but requires a memcheck-enabled MPI build. OpenMPI ships one; link against it for clean leak reports.
  • MPI profiling interface (PMPI) — every MPI function has a PMPI_ variant you can wrap. Every tracing tool uses this.
  • Score-P / Scalasca / TAU — the free tracing/profiling stack. Instruments PMPI, produces OTF2 traces viewable in Vampir or TraceCompass.

The One Trick That Fixes 40% of Crashes

When an MPI program crashes mid-run, the implementation kills all ranks, including any that were in the middle of I/O. Output is often lost. Before debugging a hard crash, add:

1
2
3
// At startup
setvbuf(stdout, NULL, _IOLBF, 0);
setvbuf(stderr, NULL, _IOLBF, 0);

Line-buffered output flushes on every newline — now a crash won’t eat your last 10 prints.


Running Under SLURM

On any cluster that uses SLURM, use srun instead of mpirun. It integrates with SLURM’s process tracking, cgroups, and PMIx for binding and spawning, so you get the resource guarantees your job requested.

1
2
3
4
5
6
7
8
9
#!/bin/bash
#SBATCH --nodes=4
#SBATCH --ntasks-per-node=32         # 128 MPI ranks total
#SBATCH --cpus-per-task=1
#SBATCH --mem-per-cpu=2G
#SBATCH --time=02:00:00

module load openmpi/4.1
srun ./my_mpi_program

For hybrid MPI + OpenMP:

1
2
3
4
5
6
7
8
#SBATCH --nodes=8
#SBATCH --ntasks-per-node=4          # 4 MPI ranks per node
#SBATCH --cpus-per-task=16           # 16 OpenMP threads per rank

export OMP_NUM_THREADS=$SLURM_CPUS_PER_TASK
export OMP_PLACES=cores
export OMP_PROC_BIND=close
srun --cpu-bind=cores ./hybrid_program

srun passes rank count, node list, and binding hints to the MPI runtime automatically — mpirun inside a SLURM job often re-queries the resource manager in ways that duplicate work or drop affinity. The one case to prefer mpirun: if you need OpenMPI-specific flags (--mca tuning parameters) that srun doesn’t translate.


mpi4py: MPI in Python

For rapid prototyping, mpi4py wraps MPI with a Pythonic API. It’s not a toy — it’s widely used in production ML and scientific codes.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
from mpi4py import MPI
import numpy as np

comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()

# Send a numpy array — zero-copy for contiguous arrays
if rank == 0:
    data = np.arange(1_000_000, dtype=np.float64)
    comm.Send(data, dest=1, tag=0)
else:
    data = np.empty(1_000_000, dtype=np.float64)
    comm.Recv(data, source=0, tag=0)

# Pickled Python objects — convenient but slower
comm.bcast({"config": "values"}, root=0)

# Collective all-reduce on a numpy array
local = np.array([rank], dtype=np.int32)
total = np.zeros(1, dtype=np.int32)
comm.Allreduce(local, total, op=MPI.SUM)
if rank == 0:
    print(f"Sum of ranks: {total[0]}")

Two conventions to remember:

  • Uppercase methods (Send, Recv, Bcast) work on buffer-protocol objects (numpy arrays) with no pickling — fast.
  • Lowercase methods (send, recv, bcast) pickle arbitrary Python objects — convenient but slow.

Launch the same way:

1
2
3
srun python my_script.py
# or
mpirun -np 4 python my_script.py

Putting It Together: A Real Iterative Solver

A compact, realistic pattern — Jacobi iteration on a 2D grid with halo exchange and global convergence check:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#include <mpi.h>
#include <math.h>
#include <stdlib.h>

int main(int argc, char **argv) {
    MPI_Init(&argc, &argv);
    int rank, size;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);

    // Set up 2D Cartesian topology
    int dims[2] = {0, 0};
    MPI_Dims_create(size, 2, dims);
    int periods[2] = {0, 0};
    MPI_Comm cart;
    MPI_Cart_create(MPI_COMM_WORLD, 2, dims, periods, 1, &cart);

    int coords[2];
    MPI_Cart_coords(cart, rank, 2, coords);

    int north, south, east, west;
    MPI_Cart_shift(cart, 0, 1, &north, &south);
    MPI_Cart_shift(cart, 1, 1, &west, &east);

    // Local grid with one-cell halo on each side
    int lnx = 256, lny = 256;
    double *u = calloc((lnx + 2) * (lny + 2), sizeof(double));
    double *un = calloc((lnx + 2) * (lny + 2), sizeof(double));

    // Column datatype for east/west halo (non-contiguous in row-major)
    MPI_Datatype col_type;
    MPI_Type_vector(lny, 1, lnx + 2, MPI_DOUBLE, &col_type);
    MPI_Type_commit(&col_type);

    double global_delta = 1.0;
    int iter = 0;
    while (global_delta > 1e-6 && iter < 10000) {
        MPI_Request reqs[8];

        // Post halo exchange (non-blocking)
        MPI_Irecv(&u[(lnx+2)*0 + 1], lnx, MPI_DOUBLE, north, 0, cart, &reqs[0]);
        MPI_Irecv(&u[(lnx+2)*(lny+1) + 1], lnx, MPI_DOUBLE, south, 0, cart, &reqs[1]);
        MPI_Irecv(&u[(lnx+2)*1 + 0], 1, col_type, west, 0, cart, &reqs[2]);
        MPI_Irecv(&u[(lnx+2)*1 + (lnx+1)], 1, col_type, east, 0, cart, &reqs[3]);

        MPI_Isend(&u[(lnx+2)*1 + 1], lnx, MPI_DOUBLE, north, 0, cart, &reqs[4]);
        MPI_Isend(&u[(lnx+2)*lny + 1], lnx, MPI_DOUBLE, south, 0, cart, &reqs[5]);
        MPI_Isend(&u[(lnx+2)*1 + 1], 1, col_type, west, 0, cart, &reqs[6]);
        MPI_Isend(&u[(lnx+2)*1 + lnx], 1, col_type, east, 0, cart, &reqs[7]);

        // Update interior cells that don't depend on halo
        double delta = 0.0;
        for (int j = 2; j <= lny - 1; j++) {
            for (int i = 2; i <= lnx - 1; i++) {
                int idx = j * (lnx + 2) + i;
                un[idx] = 0.25 * (u[idx-1] + u[idx+1] +
                                   u[idx-(lnx+2)] + u[idx+(lnx+2)]);
                double d = un[idx] - u[idx];
                delta += d * d;
            }
        }

        MPI_Waitall(8, reqs, MPI_STATUSES_IGNORE);
        // ... update boundary cells that needed halo data ...

        // Non-blocking global convergence check — overlaps with next swap
        MPI_Request conv_req;
        MPI_Iallreduce(&delta, &global_delta, 1, MPI_DOUBLE, MPI_SUM,
                       cart, &conv_req);

        double *tmp = u; u = un; un = tmp;   // swap buffers
        MPI_Wait(&conv_req, MPI_STATUS_IGNORE);
        global_delta = sqrt(global_delta);
        iter++;
    }

    if (rank == 0) printf("Converged in %d iterations\n", iter);

    MPI_Type_free(&col_type);
    free(u); free(un);
    MPI_Finalize();
    return 0;
}

Every pattern in this post appears: Cartesian topology for neighbor discovery, derived datatype for strided columns, non-blocking halo exchange overlapping with interior computation, non-blocking all-reduce overlapping with the buffer swap, collective convergence check.


Conclusion

MPI earns its place at the center of HPC not by being elegant, but by being predictable and portable. Code written to the standard runs on a laptop, a 16-node lab cluster, and the world’s largest supercomputers with identical semantics. The runtime can use TCP, InfiniBand, Slingshot, or shared memory — your source doesn’t change.

Ninety percent of real MPI code uses a tiny fraction of the API: Init/Finalize, Comm_rank/Comm_size, Send/Recv, the core collectives (Bcast, Reduce, Allreduce, Scatter, Gather), their non-blocking variants, a derived datatype or two, and Cartesian topology for structured grids. Master that subset and you can write most HPC applications well.

The pitfalls that burn the most engineering time — eager/rendezvous transitions, load imbalance, over-eager barriers, unbound processes, naïve fan-out patterns instead of collectives — all come from not understanding what the runtime does under the hood. Read your implementation’s tuning guide once a year. Profile with mpiP or Score-P before optimizing. And when things don’t make sense, the answer is almost always “someone’s MPI_Send is blocking waiting for a matching receive that will never be posted.”

Comments