MPI from Zero: Ranks, Messages and a 1-D Heat Equation with Halo Exchange

Most of the earlier posts in this series ran inside one address space: threads shared an array and the cache-coherence protocol kept their copies of it consistent. A cluster has no shared array. Each node has its own memory, and the only way data moves is that one process explicitly sends it and another explicitly receives it. MPI, the Message Passing Interface, is the standard API for that model, and it has stayed the standard since its first version in 1994, largely because it maps so directly onto the hardware: a message is a copy, and the programmer decides when copies happen.

WARNING

I have no cluster. The code below was compiled against Open MPI 5.0.9 and run with 1, 2 and 4 ranks on an Apple-silicon laptop, where all ranks live on one machine and messages travel through shared memory. That is enough to check correctness, deadlock-freedom and the decomposition logic; it says nothing about performance on a real interconnect, so there are no bandwidth or scaling claims in this post.

The programming model

An MPI program is one executable started $p$ times by a launcher (mpirun -np 4 ./heat1d). Each copy is a rank, numbered $0 \dots p - 1$ inside a communicator, and MPI_COMM_WORLD is the communicator that contains everyone. The program is the same for all ranks (SPMD: single program, multiple data); behaviour differs only because each rank reads its own number and acts on its own slice of the problem.

#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("rank %d of %d\n", rank, size);
    MPI_Finalize();
    return 0;
}

The output order of the four printf lines is undefined; nothing in MPI orders the ranks' side effects, and getting used to that is the first mental shift. Communicators matter beyond MPI_COMM_WORLD because libraries create their own (so their messages cannot be confused with yours) and because a sub-communicator is how you say "the ranks in this row of the process grid".

Point-to-point: send, receive, match

A message is matched by three things: the source rank, the tag (an integer the program chooses) and the communicator. The count and datatype describe the buffer; MPI converts representations between heterogeneous machines, which is why you say MPI_DOUBLE instead of sizeof(double).

Call	Blocking?	Completes when
`MPI_Send`	yes	the buffer may be reused; the library may or may not have delivered it
`MPI_Recv`	yes	the message is in the buffer
`MPI_Sendrecv`	yes	both the send and the receive have completed; safe for exchanges
`MPI_Isend` / `MPI_Irecv`	no, returns a request	never by itself; `MPI_Wait` or `MPI_Test` completes it
`MPI_Bcast`, `MPI_Reduce`, `MPI_Allreduce`, `MPI_Scatter`, `MPI_Gather`	collective	this rank's part is done and its buffers may be reused; others may not have started
`MPI_Barrier`	collective	every rank has arrived; rarely needed for correctness, useful for timing

The subtle entry is MPI_Send. The standard allows it to return as soon as the data has been copied into an internal buffer (the "eager" protocol used for small messages) or to block until the matching receive is posted (the "rendezvous" protocol for large messages). Both are legal, and which one you get depends on the library, the message size and the transport. Code that works only because the eager path happened to be taken is a deadlock waiting for a larger input. Collectives have the same kind of latitude: a rank may leave MPI_Bcast or MPI_Reduce before other ranks have even entered it, so none of them can stand in for a barrier. Only MPI_Barrier promises that every rank has arrived.

The exchange everyone writes first (left) relies on buffering the standard does not promise. MPI_Sendrecv (middle) or non-blocking calls (right) are the two idiomatic fixes.

Collectives: say what you mean

When every rank participates in the same operation, a collective is both clearer and faster than a hand-written loop of sends. MPI_Bcast distributes one buffer from a root to all; MPI_Reduce combines one value per rank with an operation (MPI_SUM, MPI_MAX, …) onto the root; MPI_Allreduce does the same and gives everyone the answer. Libraries implement these with tree or butterfly algorithms, so the cost of a reduction over $p$ ranks grows like $\log_{2} p$ message latencies, not $p$ :

T_{reduce} (p) \approx ⌈ \log_{2} p ⌉ (α + β m),

where $α$ is the per-message latency, $β$ the inverse bandwidth and $m$ the message size. The heat solver below uses two reductions at the end: MPI_MAX of the per-rank error and MPI_MAX of the per-rank elapsed time.

The 1-D heat equation, decomposed

The model problem is heat conduction along a rod with its ends held at fixed temperatures. Discretising $u_{t} = α u_{x x}$ with forward differences in time and central differences in space gives the explicit update

u_{i}^{t + 1} = u_{i}^{t} + r (u_{i - 1}^{t} - 2 u_{i}^{t} + u_{i + 1}^{t}), r = \frac{α Δ t}{Δ x^{2}} \leq \frac{1}{2},

where the bound on $r$ is the stability condition of the scheme. Each new value needs its two neighbours' old values, and that is the whole communication pattern: a rank that owns points $1 \dots n$ needs point $0$ from its left neighbour and point $n + 1$ from its right neighbour at every step. Those two extra slots are ghost cells (or halo cells), and refreshing them each step is a halo exchange.

One-dimensional domain decomposition. Each rank's local array has n interior points plus two ghost cells; the ghosts are the neighbours' edge values and are refreshed every time step.

The interior points $2 \dots n - 1$ depend only on local data, so they can be updated while the ghost messages are in flight; only $u_{1}$ and $u_{n}$ have to wait. That is what the non-blocking version buys.

/* heat1d.c : explicit 1-D heat equation with a halo exchange between MPI ranks.
 *   mpicc -O2 -std=c11 -o heat1d heat1d.c
 *   mpirun -np 4 ./heat1d 512 1000000
 */
#include <mpi.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.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);

    int N     = argc > 1 ? atoi(argv[1]) : 1024;    /* interior points, global */
    int steps = argc > 2 ? atoi(argv[2]) : 20000;   /* time steps              */
    if (N % size != 0) {
        if (rank == 0) fprintf(stderr, "N=%d must be divisible by %d ranks\n", N, size);
        MPI_Abort(MPI_COMM_WORLD, 1);
    }
    int n = N / size;
    double dx = 1.0 / (N + 1), dt = 0.4 * dx * dx, r = dt / (dx * dx);   /* r = 0.4 <= 1/2 */

    /* neighbours; MPI_PROC_NULL turns the outermost sends and receives into no-ops */
    int left  = rank == 0        ? MPI_PROC_NULL : rank - 1;
    int right = rank == size - 1 ? MPI_PROC_NULL : rank + 1;

    double *u = calloc(n + 2, sizeof *u), *unew = calloc(n + 2, sizeof *unew);
    if (rank == 0)        u[0]     = 1.0;   /* T_L: a physical boundary, never overwritten */
    if (rank == size - 1) u[n + 1] = 0.0;   /* T_R */

    double t0 = MPI_Wtime();
    for (int s = 0; s < steps; s++) {
        MPI_Request req[4];
        /* 1. post receives first so arriving data has somewhere to land */
        MPI_Irecv(&u[0],     1, MPI_DOUBLE, left,  0, MPI_COMM_WORLD, &req[0]);
        MPI_Irecv(&u[n + 1], 1, MPI_DOUBLE, right, 1, MPI_COMM_WORLD, &req[1]);
        /* 2. send my edge values into the neighbours' ghost cells */
        MPI_Isend(&u[1],     1, MPI_DOUBLE, left,  1, MPI_COMM_WORLD, &req[2]);
        MPI_Isend(&u[n],     1, MPI_DOUBLE, right, 0, MPI_COMM_WORLD, &req[3]);
        /* 3. interior points need only local data: overlap with the messages */
        for (int i = 2; i <= n - 1; i++)
            unew[i] = u[i] + r * (u[i - 1] - 2.0 * u[i] + u[i + 1]);
        /* 4. the two edge points need the ghosts */
        MPI_Waitall(4, req, MPI_STATUSES_IGNORE);
        unew[1] = u[1] + r * (u[0]     - 2.0 * u[1] + u[2]);
        unew[n] = u[n] + r * (u[n - 1] - 2.0 * u[n] + u[n + 1]);
        unew[0] = u[0]; unew[n + 1] = u[n + 1];     /* keep boundaries; ghosts are refreshed anyway */
        double *tmp = u; u = unew; unew = tmp;
    }
    double elapsed = MPI_Wtime() - t0;

    /* check against the steady state u(x) = 1 - x */
    double local_err = 0.0;
    for (int i = 1; i <= n; i++) {
        double x = (rank * n + i) * dx;
        local_err = fmax(local_err, fabs(u[i] - (1.0 - x)));
    }
    double max_err, max_time;
    MPI_Reduce(&local_err, &max_err,  1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
    MPI_Reduce(&elapsed,   &max_time, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
    if (rank == 0)
        printf("ranks=%d N=%d steps=%d  max |u - (1-x)| = %.3e  time %.3f s\n",
               size, N, steps, max_err, max_time);

    free(u); free(unew);
    MPI_Finalize();
    return 0;
}

Four details carry the correctness of this loop:

Tags label the directions. A value sent to the right with tag 0 is received from the left with tag 0; the leftward message uses tag 1. In this chain they document intent more than they prevent errors: every receive names its source rank, and an interior rank's left and right neighbours are different ranks, so even a single tag could not put the left neighbour's value into the right ghost cell. Tags become necessary when two messages from the same source must be told apart, as in a periodic ring of two ranks, where each rank's left and right neighbour are the same process.
MPI_PROC_NULL removes the edge cases. Rank 0 has no left neighbour, so its left send and receive are addressed to MPI_PROC_NULL, which completes immediately and touches nothing. The code has no if (rank == 0) inside the time loop.
Receives are posted before sends. A message that arrives before its receive exists is an "unexpected message". A small one sent eagerly has to be parked in a temporary buffer and copied again when the receive appears; a large one under the rendezvous protocol simply waits until the receive is posted. Posting receives first cannot guarantee direct delivery, since a faster neighbour may send before this rank reaches its receives, but it makes landing straight in u[0] or u[n+1] the common case.
No buffer is touched between Isend/Irecv and Waitall. The interior loop reads u[1] and u[n] (being sent) but never writes them, and it writes unew, not the ghost cells (being received). Writing into an in-flight buffer is undefined behaviour in MPI and produces results that are wrong only sometimes.

The same exchange with MPI_Sendrecv, and why the naive version deadlocks

If overlap is not needed, one blocking call per direction is simplest and provably safe:

MPI_Sendrecv(&u[n], 1, MPI_DOUBLE, right, 0,          /* send my right edge rightwards   */
             &u[0], 1, MPI_DOUBLE, left,  0,          /* receive left ghost from the left */
             MPI_COMM_WORLD, MPI_STATUS_IGNORE);
MPI_Sendrecv(&u[1],     1, MPI_DOUBLE, left,  1,
             &u[n + 1], 1, MPI_DOUBLE, right, 1,
             MPI_COMM_WORLD, MPI_STATUS_IGNORE);

The version to avoid is MPI_Send(right); MPI_Send(left); MPI_Recv(left); MPI_Recv(right); on every rank. With one double per message it will almost certainly work, because the library buffers it eagerly. Change the halo to a few thousand elements (a 2-D or 3-D stencil, or a halo of width 2 for a higher-order scheme) and every rank blocks in its first MPI_Send, waiting for a receive nobody has reached. The other classic fix is to break the symmetry: even ranks send first and odd ranks receive first.

What the laptop run showed

With $N = 512$ interior points and one million steps, the solution is within $2 \times 10^{- 7}$ of the steady state $u (x) = 1 - x$ (what remains is the slowest mode still decaying, not rounding), and the three runs agree in every printed digit because every grid point sees exactly the same arithmetic regardless of which rank owns it. Built with Open MPI 5.0.9 and Apple clang 17.0.0 (-O2), run with mpirun --oversubscribe -np 1|2|4:

ranks=1 N=512 steps=1000000  max |u - (1-x)| = 1.945e-07  time 0.142 s
ranks=2 N=512 steps=1000000  max |u - (1-x)| = 1.945e-07  time 0.201 s
ranks=4 N=512 steps=1000000  max |u - (1-x)| = 1.945e-07  time 0.367 s

The times are included because the program prints them, not because they mean anything about MPI performance. With 128 points per rank the halo exchange is two messages per step against a few hundred flops, a communication-to-computation ratio of $2 / n$ ; on a laptop where ranks compete for cores and the heavy ranks are scheduled unpredictably, adding ranks makes this tiny problem slower, which is the expected outcome and the reason real stencil codes use large subdomains and 2-D or 3-D decompositions whose surface-to-volume ratio shrinks as the subdomain grows.

TIP

Build the habit of checking an MPI code with a single rank first (mpirun -np 1): with MPI_PROC_NULL on both sides the halo exchange disappears and any difference from the serial answer is a bug in the decomposition, not in the physics. Then run with a rank count that does not divide the problem evenly, or with a count larger than the number of cores, to flush out assumptions. The identical error printed for 1, 2 and 4 ranks above is this check passing.

The limits of non-blocking overlap

The loop above overlaps the interior update with the ghost messages, and in a textbook the picture is that communication becomes free. Two caveats:

Progress is not guaranteed in the background. Many MPI implementations only move data when the application is inside an MPI call. If the interior loop runs for ten milliseconds, a large message may sit untouched until MPI_Waitall. Implementations offer progress threads or asynchronous-progress options, and the only way to know is to measure the time inside MPI_Waitall with and without the overlap.
The overlap needs enough interior. With $n = 128$ there are 126 interior points, a few hundred nanoseconds of work, less than the latency of a message even over shared memory. Overlap pays when the local computation is at least as long as the communication it hides, which in stencil codes means subdomains of thousands of points per dimension.

The structure is still the right one to write, because it costs nothing when it does not help and it is the structure production stencil codes use: post receives, post sends, compute what you can, wait, finish the edges.

Where to go next

The two-dimensional version of this solver is the natural follow-up. MPI_Cart_create builds a process grid and MPI_Cart_shift returns the neighbour ranks (including MPI_PROC_NULL at the edges) so that the decomposition code stays as short as the 1-D case; column halos are non-contiguous in memory and are described with MPI_Type_vector so that the library gathers them without a manual pack loop. The same pattern, with the stencil replaced by a sparse matrix, is a parallel conjugate-gradient solver; with the stencil replaced by a convolution, it is the spatial partitioning used to train networks on images too large for one device.

MPI from Zero: Ranks, Messages and a 1-D Heat Equation with Halo Exchange

WARNING