Chapter 34: MPI — Distributed-Memory Parallelism for Cluster Computing

"MPI is the assembly language of parallel computing." — a common saying among HPC programmers

Overview

Everything your solver has done in parallel so far has happened inside a single computer. Coarrays in Chapter 32 split the plate across images, and OpenMP in Chapter 33 turned the update loop loose on the threads of one node — both models where, ultimately, there is one shared pool of memory that the parallel workers can reach. That model has a hard ceiling: a single node, however fat, holds only so many cores and so much RAM. When the weather service runs a global forecast, or the astrophysicist evolves a galaxy, the computation does not fit on one machine and never will. It runs across a cluster — hundreds or thousands of separate computers, each with its own private memory, wired together by a fast network. To program that machine you need a way for those separate computers to cooperate, and for thirty years the answer, across essentially every supercomputer on Earth, has been one thing: MPI, the Message Passing Interface.

This is the chapter where your solver leaves the single node behind and becomes a genuine cluster code. It is also, not coincidentally, the model that runs the production Fortran named in Chapter 1 — WRF, CESM, and the quantum-chemistry and astrophysics codes that consume most of the cycles on the TOP500 machines are, overwhelmingly, Fortran plus MPI. Learning MPI is learning the lingua franca of real high-performance computing. It is more work than OpenMP — there is no sprinkling a directive on a loop; you must decide who owns which data and write every exchange by hand — but that explicitness is exactly what lets it scale to a million cores when shared-memory models cannot. Fortran is not dead; it is running, right now, on the largest machines humanity has built, and MPI is how it talks to itself across them.

The intellectual core of the chapter is a single pattern — domain decomposition with halo exchange — that you have been quietly prepared for since Chapter 24. The five-point stencil there reaches only to a cell's four nearest neighbours. Cut the plate into strips, give each process one strip, and the only thing a process needs from its neighbours is the single row of cells along each shared edge. Exchange those edge rows each step and every process can march its own strip forward as if it owned the whole plate. That is the whole idea, and it is how a stencil code — the workhorse of computational physics — is distributed across a cluster.

In this chapter, you will learn to:

• Describe the MPI model — many processes, each with private memory, coordinated by messages — and write the mpi_init … mpi_finalize skeleton that every MPI program shares.
• Send and receive data between processes with point-to-point calls, get the argument order exactly right, and diagnose and cure the notorious send/recv deadlock.
• Reach for the right collective — broadcast, reduce, all-reduce, gather, scatter — instead of building those patterns by hand.
• Decompose the heat plate across processes, wrap each piece in ghost cells, and exchange halos each step with mpi_sendrecv, deadlock-free — the Project Checkpoint that makes your solver a cluster code.
• Compile with mpif90, launch with mpirun -np N, and place MPI-IO and hybrid MPI+OpenMP on the map.

Learning Paths

How to read this chapter by track. - ⚡ HPC ("I need cluster code") — this is the chapter you came for; read every section. §34.4 (domain decomposition and halo exchange) and the Project Checkpoint are the pattern behind essentially every distributed stencil code in production. Hand-trace the small example; that is where the understanding lives. - 🔬 Scientist ("my simulation outgrew one machine") — read §34.1 for the model, §34.4 for how a grid is split and stitched, and the Checkpoint. Skim §34.2–34.3 for the vocabulary you will meet in others' code. - 📖 Standard — read straight through. MPI is a library, not part of the Fortran standard, so this chapter is about a binding — how the Fortran interface expresses the message-passing model. - 🔧 Legacy ("I inherited an MPI code") — §34.1–34.3 decode the mpi_* calls you are staring at, and the ⚠️ deadlock pitfall in §34.2 explains the hang you may have hit. §34.5 covers the use mpi vs mpi_f08 interface question you will face when modernising it.

34.1 The MPI Model: Processes, Communicators, and Ranks

Start with what is different about this model, because it is the source of everything else. In OpenMP, several threads shared one address space: a variable written by one thread could be read by another, and the whole difficulty was keeping them from colliding. MPI is the opposite world. Here the parallel workers are not threads but processes — separate running programs, each with its own private memory — and no process can touch another's variables at all. If process 0 has a value that process 1 needs, process 0 must package it up and send it, and process 1 must receive it. There is no shared state; there is only communication. This is the distributed-memory model from Chapter 31, and MPI is its dominant expression.

Definition (MPI). MPI — the Message Passing Interface — is a standardised library for distributed-memory parallel programming, in which many independent processes, each with private memory, coordinate by explicitly sending and receiving messages. It is not part of any language; it is a specification (first standardised in 1994) with bindings for Fortran, C, and C++, implemented by libraries such as Open MPI and MPICH. A Fortran program uses it by use mpi and by calling mpi_* procedures. MPI is the substrate of nearly all large-scale scientific computing: when a code runs across the thousands of nodes of a cluster, MPI is almost always what carries the data between them.

Every MPI program follows the same shape. You start the MPI runtime, ask it who you are and how many processes there are, do your work with that knowledge, and shut MPI down. Here is that skeleton — the "hello, world" of message passing — and it already introduces the four calls you will write in every MPI program you ever build.

program mpi_hello
  use mpi
  implicit none
  integer :: ierr, rank, nprocs

  call mpi_init(ierr)                              ! start the MPI runtime
  call mpi_comm_rank(MPI_COMM_WORLD, rank,   ierr) ! which process am I? (0-based)
  call mpi_comm_size(MPI_COMM_WORLD, nprocs, ierr) ! how many of us are there?

  print '(a,i0,a,i0)', 'Hello from rank ', rank, ' of ', nprocs

  call mpi_finalize(ierr)                          ! shut the MPI runtime down
end program mpi_hello

$ mpif90 -std=f2018 -Wall -O2 example-01-hello.f90 -o hello && mpirun -np 4 ./hello
Hello from rank 0 of 4
Hello from rank 1 of 4
Hello from rank 2 of 4
Hello from rank 3 of 4

Several things in those eleven lines deserve unpacking, because they define the whole model.

One program, many processes. You wrote one program, but mpirun -np 4 launched four copies of it, running at the same time, each in its own memory. Every copy executes the identical source — the same print statement runs four times — yet they are not clones, because each one gets a different answer from mpi_comm_rank. This is the SPMD model: Single Program, Multiple Data. One program text, run by many processes, each specialising its behaviour based on which process it is. Almost all MPI programs are SPMD; the if (rank == 0) branch that does something special on one process is the model's signature move.

Definition (rank). A rank is a process's integer identity within a group of MPI processes, numbered from 0 to nprocs - 1. It is how a process knows which of the running copies it is, and how other processes address it: to send data "to rank 2" is to name its rank. mpi_comm_rank returns the calling process's own rank; mpi_comm_size returns how many processes are in the group. Between them, those two numbers are the entire basis on which an SPMD program divides its work — rank 0 takes the first slice, rank 1 the next, and so on.

Communicators name the group. Notice MPI_COMM_WORLD in both queries. That is a communicator, and it is MPI's way of naming a set of processes that can talk to each other.

Definition (communicator). A communicator is an MPI object identifying a group of processes and providing a private communication context for them. Every message is sent within a communicator, and a rank is only meaningful relative to one — "rank 2" means rank 2 in this communicator. The predefined MPI_COMM_WORLD contains every process the program was launched with, and it is the only communicator most programs ever need. You can create sub-communicators to split the processes into cooperating teams (say, one per subdomain), but that is an advanced move; for this chapter, communicator means MPI_COMM_WORLD. In the use mpi interface a communicator is just a default integer handle.

The ierr last-argument convention. Every mpi_* call ends with an extra integer argument, ierr, into which MPI writes a status code (MPI_SUCCESS on success). This is the Fortran binding's universal convention and it is not optional in the use mpi interface — leave it off and the call is wrong. It looks alien beside modern Fortran's intent(out) results, and it is: MPI is a C-flavoured library from 1994, and the Fortran binding mirrors the C functions' integer return codes as a trailing argument. In practice almost nobody checks ierr (a failed MPI call usually aborts the whole job anyway), but you must always pass it. We will meet the more modern mpi_f08 interface, which makes ierr genuinely optional, in §34.5.

⚠️ Common Pitfall — forgetting ierr, or checking the wrong thing. The single most common beginner error in Fortran MPI is omitting the trailing ierr argument — call mpi_init() will not compile against the use mpi interface, and call mpi_comm_rank(MPI_COMM_WORLD, rank) is missing its last argument. Every mpi_* subroutine in this chapter takes ierr last, without exception. A subtler trap: the value rank is written into the second-to-last argument and ierr into the last, so mpi_comm_rank(comm, rank, ierr) — get those two in the wrong order and you will "successfully" read the error code into your rank variable and wonder why every process thinks it is rank 0.

📜 From History: how MPI came to be. Before 1994 every parallel-computer vendor shipped its own incompatible message-passing library, so a code written for one machine had to be rewritten for the next — a catastrophe for scientific software meant to outlive the hardware. The MPI Forum, a consortium of vendors, national labs, and academics, met through 1993–94 to agree on a single portable interface, and MPI-1 was the result. The bet paid off spectacularly: three decades later, the same MPI calls compile and run on a laptop and on the largest supercomputer in the world. It is the same standardisation instinct that keeps fifty-year-old Fortran running (Chapter 1's ⭐ From History) — the unglamorous act of agreeing on an interface is what lets scientific software accumulate instead of evaporate.

🔄 Check Your Understanding. 1. In the MPI model, can process 1 read a variable that process 0 declared and set? How does it obtain the value if it needs it? 2. What do mpi_comm_rank and mpi_comm_size return, and why does an SPMD program need both? 3. Why does call mpi_init() fail to compile against the use mpi interface?

Answers

(1) No — each process has private memory; there is no shared state. Process 0 must mpi_send the value and process 1 must mpi_recv it. (2) mpi_comm_rank returns the calling process's own rank (its identity, 0-based); mpi_comm_size returns the total number of processes. An SPMD program runs one identical text on all processes, so each must learn which it is (rank) and how many there are (size) to take its own slice of the work. (3) Every mpi_* call requires the trailing ierr status argument; mpi_init() is missing it. The correct call is call mpi_init(ierr).

34.2 Point-to-Point Communication: Send, Receive, and Deadlock

With processes launched and each knowing its rank, the first real question is how two of them exchange data. The most basic operation is point-to-point communication: one named process sends a message, one named process receives it.

Definition (point-to-point communication). Point-to-point communication is the exchange of a message between exactly two processes — one calls a send operation naming the receiver, the other calls a matching receive operation naming the sender. It is the fundamental building block of MPI; every more elaborate pattern can, in principle, be built from it. The core pair is mpi_send and mpi_recv.

Here is the archetype: rank 0 packages an array and sends it to rank 1, which receives it and works with it.

program point_to_point
  use mpi
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer  :: ierr, rank, nprocs
  integer  :: status(MPI_STATUS_SIZE)              ! receive status (use mpi interface)
  real(dp) :: buf(3)

  call mpi_init(ierr)
  call mpi_comm_rank(MPI_COMM_WORLD, rank,   ierr)
  call mpi_comm_size(MPI_COMM_WORLD, nprocs, ierr)

  if (rank == 0) then
    buf = [10.0_dp, 20.0_dp, 30.0_dp]
    call mpi_send(buf, 3, MPI_DOUBLE_PRECISION, 1, 0, MPI_COMM_WORLD, ierr)
    !              buf  count  datatype        dest tag  comm         ierr
  else if (rank == 1) then
    call mpi_recv(buf, 3, MPI_DOUBLE_PRECISION, 0, 0, MPI_COMM_WORLD, status, ierr)
    !              buf  count  datatype        src tag  comm         status  ierr
    print '(a,f6.1)', 'rank 1 received, sum = ', sum(buf)
  end if

  call mpi_finalize(ierr)
end program point_to_point

$ mpif90 -std=f2018 -Wall -O2 example-02-send-recv.f90 -o pp && mpirun -np 2 ./pp
rank 1 received, sum =   60.0

The two calls are the heart of MPI, so learn their arguments cold. mpi_send takes, in order: the buffer (the data), the count (how many elements), the datatype, the destination rank, a tag, the communicator, and ierr. mpi_recv takes the same first three, then the source rank, the tag, the communicator, a status, and ierr. Six things must line up for a message to be delivered, and getting any of them wrong is the source of most MPI bugs:

The datatype deserves a word, because it is where Fortran meets MPI's C heritage. MPI does not know your Fortran kinds; you must tell it, by hand, what kind of element you are sending, using one of its named constants. Our project's reals are real(dp) with dp = selected_real_kind(15, 307), i.e. IEEE double precision, so the matching MPI datatype is MPI_DOUBLE_PRECISION. Send an array of real(dp) and declare it MPI_INTEGER, and MPI will faithfully copy the wrong bytes and hand the receiver garbage — with no error, because MPI trusts you. This is the first of several places where MPI is powerful and unforgiving in equal measure.

💡 Intuition: think of mpi_send/mpi_recv as posting a parcel. The buffer is what is in the box; the count and datatype are the customs declaration ("three double-precision reals"); the destination is the address; the tag is a reference number written on the outside so the recipient can tell this parcel apart from others you sent them; and the communicator is the postal system you are using. The parcel is delivered only when all of these agree between sender and receiver. The status the receiver gets back is the delivery slip — it records who really sent it and under what tag, useful when you accepted a parcel from MPI_ANY_SOURCE.

The deadlock trap

Now the mistake that every MPI programmer makes exactly once, and never forgets. Suppose two neighbouring processes each need to give the other a row of data — precisely the halo exchange your solver will need. The obvious code has each process send its row first, then receive the neighbour's:

! BOTH rank 0 and rank 1 run this -- and it can HANG FOREVER
call mpi_send(mine,  n, MPI_DOUBLE_PRECISION, other, 0, MPI_COMM_WORLD, ierr)
call mpi_recv(yours, n, MPI_DOUBLE_PRECISION, other, 0, MPI_COMM_WORLD, status, ierr)

This looks symmetric and correct, and it is a deadlock waiting to happen.

⚠️ Common Pitfall — the send/recv deadlock (and why it "works" in testing). When both processes call mpi_send first, each blocks inside its send waiting for the other to receive — but neither can reach its mpi_recv, because both are stuck in mpi_send. Each waits for the other forever; the program hangs with no error. The vicious part is that it often runs fine in testing. mpi_send is "standard mode," which means MPI may copy small messages into an internal buffer and return immediately (the "eager" protocol), so for a few elements both sends complete, both recvs run, and everything works. Cross the implementation's buffer threshold — a large enough row, which real problems always have — and mpi_send switches to the "rendezvous" protocol, blocking until the matching recv is posted, and the code deadlocks in production after passing every small test. Correctness must never depend on MPI buffering it for you. There are three clean fixes:

```fortran ! Fix 1 -- order the calls: one side sends-then-receives, the other receives-then-sends. if (mod(rank,2) == 0) then call mpi_send(mine, n, MPI_DOUBLE_PRECISION, other, 0, MPI_COMM_WORLD, ierr) call mpi_recv(yours, n, MPI_DOUBLE_PRECISION, other, 0, MPI_COMM_WORLD, status, ierr) else call mpi_recv(yours, n, MPI_DOUBLE_PRECISION, other, 0, MPI_COMM_WORLD, status, ierr) call mpi_send(mine, n, MPI_DOUBLE_PRECISION, other, 0, MPI_COMM_WORLD, ierr) end if

! Fix 2 -- ONE call that does both, which MPI guarantees will not deadlock: call mpi_sendrecv(mine, n, MPI_DOUBLE_PRECISION, other, 0, & ! send half yours, n, MPI_DOUBLE_PRECISION, other, 0, & ! recv half MPI_COMM_WORLD, status, ierr) ```

mpi_sendrecv exists for exactly this situation, so prefer it: it is one call, it cannot deadlock, and it reads as what it is — an exchange. (Fix 3, non-blocking mpi_isend/mpi_irecv, comes in §34.4.) The project's halo exchange uses mpi_sendrecv.

mpi_sendrecv bundles a send and a receive into a single deadlock-free call. Its argument list is just mpi_send's followed by mpi_recv's, sharing one communicator and status at the end: send buffer, count, type, dest, sendtag; then recv buffer, count, type, source, recvtag; then comm, status, ierr. It is the right tool for any structured pairwise exchange, and it is the backbone of stencil halo swapping.

🐛 Find the Bug. A colleague's two-process exchange runs perfectly on their laptop with a test grid of 100 cells per row, then hangs on the cluster with 4000 cells per row. They swear "nothing changed but the problem size." What is happening, and what is the one-line change that fixes it for good?

Answer

It is the send-first/recv-second deadlock above. At 100 elements the message fit under the MPI implementation's eager-buffer threshold, so mpi_send returned immediately and the code worked; at 4000 elements the message exceeds the threshold, mpi_send switches to rendezvous mode and blocks until a matching recv is posted, and since both ranks are stuck in mpi_send, they wait forever. The fix is to replace the separate send-then-recv with a single mpi_sendrecv (or to order the calls so one side receives first). The bug was always there; the larger problem merely exposed it — which is why you never rely on buffering for correctness.

🔄 Check Your Understanding. 1. List the seven arguments of mpi_send, in order. 2. Why can a send-then-receive exchange pass every small test and then deadlock on a large production run? 3. What single call replaces a paired send-and-receive and is guaranteed not to deadlock?

Answers

(1) buffer, count, datatype, destination rank, tag, communicator, ierr. (2) mpi_send is standard mode: for small messages MPI may buffer the data and return immediately, so both sends complete and the recvs run; for large messages it blocks until the matching recv is posted, and if both ranks send first, neither ever posts a recv — deadlock. Correctness must not depend on buffering. (3) mpi_sendrecv.

34.3 Collective Communication: Broadcast, Reduce, Gather, Scatter

Point-to-point is the foundation, but a great deal of what parallel programs actually do is not one process talking to one other — it is all the processes cooperating in a pattern: everyone needs the same configuration value; everyone contributes a partial sum that must be totalled; one process holds an array that must be dealt out to all. You could build these from loops of mpi_send/mpi_recv, and it would be tedious, slow, and buggy. MPI provides them directly, as collective operations.

Definition (collective communication). A collective operation is one that all processes in a communicator call together, cooperating in a single communication pattern — broadcasting, combining, gathering, or scattering data across the whole group. Every process must make the matching call (a collective that some processes skip will hang), and in return MPI implements the pattern efficiently, often with clever tree algorithms that no hand-rolled loop of point-to-point calls would match. Collectives are both clearer (one call states the whole intent) and faster than doing it yourself.

The five you will use constantly:

Here is a program using the two most common — a broadcast to hand everyone the same starting value, and an all-reduce to combine a per-process contribution into a global total that everyone then knows:

program collectives
  use mpi
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer  :: ierr, rank, nprocs
  real(dp) :: base, mine, total, biggest

  call mpi_init(ierr)
  call mpi_comm_rank(MPI_COMM_WORLD, rank,   ierr)
  call mpi_comm_size(MPI_COMM_WORLD, nprocs, ierr)

  if (rank == 0) base = 100.0_dp                  ! only the root has a meaningful value...
  call mpi_bcast(base, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)   ! ...now everyone does

  mine = real(rank + 1, dp)                        ! each rank contributes rank+1
  call mpi_allreduce(mine, total,   1, MPI_DOUBLE_PRECISION, MPI_SUM, MPI_COMM_WORLD, ierr)
  call mpi_allreduce(mine, biggest, 1, MPI_DOUBLE_PRECISION, MPI_MAX, MPI_COMM_WORLD, ierr)

  if (rank == 0) then
    print '(a,f6.1)', 'broadcast base = ', base
    print '(a,f6.1)', 'sum of (rank+1) over all ranks = ', total
    print '(a,f6.1)', 'max of (rank+1) over all ranks = ', biggest
  end if

  call mpi_finalize(ierr)
end program collectives

$ mpif90 -std=f2018 -Wall -O2 example-03-collectives.f90 -o coll && mpirun -np 4 ./coll
broadcast base = 100.0
sum of (rank+1) over all ranks = 10.0
max of (rank+1) over all ranks = 4.0

Work the output by hand to trust it, and note which parts depend on the process count. mpi_bcast copies rank 0's base = 100.0 to every process, so all four now hold 100.0. Each rank sets mine = rank + 1, so the four contributions are 1, 2, 3, 4; mpi_allreduce with MPI_SUM totals them to 1+2+3+4 = 10, and with MPI_MAX returns the largest, 4. With -np N the sum would be N(N+1)/2 and the max would be N — those two numbers are process-count-dependent; the broadcast value 100.0 is not. The crucial word is all-reduce: after the call, every process holds total and biggest, not just rank 0. Had we used mpi_reduce (which takes an extra root argument), only the named root would hold the answer and the others would have garbage — the right choice when only one process needs the result, e.g. to print it.

The argument shapes are worth fixing in memory, because they differ in small but important ways:

• mpi_bcast(buffer, count, datatype, root, comm, ierr) — one buffer, played from root to all. No separate send/recv buffers: the root's buffer is the source, everyone else's buffer is the destination.
• mpi_reduce(sendbuf, recvbuf, count, datatype, op, root, comm, ierr) — separate in/out buffers, an operation (MPI_SUM, MPI_MAX, MPI_MIN, MPI_PROD, …), and a root; the result lands only on root.
• mpi_allreduce(sendbuf, recvbuf, count, datatype, op, comm, ierr) — identical to reduce but with no root, because everyone gets the result.

⚡ Performance Note — collectives are not just convenient, they are fast. A naive "sum everyone's value onto rank 0" written by hand sends N - 1 messages one after another to rank 0, taking time proportional to N. A good MPI library implements mpi_reduce as a tree: pairs combine, then pairs of pairs, and so on, finishing in time proportional to $\log_2 N$ — for 1024 processes, 10 steps instead of 1023. The library also knows the machine's network topology and can route the combination along it. This is the same lesson as Chapter 21's "don't hand-roll what a tuned library does better": performance is not accidental, and a collective hands the implementation the whole pattern to optimise, exactly as a whole-array operation hands the compiler the whole loop. Always prefer a collective to a hand-built equivalent.

🐍 Python Comparison. The same patterns exist in mpi4py, the Python binding: comm.bcast, comm.reduce, comm.allreduce, comm.Scatter, comm.Gather. mpi4py is genuinely useful for orchestration and for the parts of a workflow that are not performance-critical — and, true to this book's refrain, the productive pattern is Fortran and Python better together: run the heavy stencil in Fortran+MPI and script the launch, the parameter sweep, or the post-processing in Python. But for the inner loop of a real simulation, the Fortran binding wins decisively, because there is no interpreter between your code and the network, and because the data is already in the contiguous, column-major arrays MPI wants to send. Python is the conductor; Fortran+MPI is the orchestra.

🔄 Check Your Understanding. 1. What is the difference between mpi_reduce and mpi_allreduce, and when would you choose each? 2. In the program above, if you launch with -np 8, what will the sum and max lines print? 3. Why must every process in the communicator call a collective, even one that seems to only "receive"?

Answers

(1) mpi_reduce delivers the combined result to a single root process only; mpi_allreduce delivers it to every process. Use reduce when one process needs the answer (say, to print or write it); use allreduce when every process needs it to proceed (say, a global convergence test each step). (2) Sum $= 1+2+\dots+8 = N(N+1)/2 = 36$; max $= N = 8$ — both scale with the process count. (3) A collective is a single cooperative operation across the whole group; MPI's algorithm (e.g. a tree) needs every process to participate at its point in the pattern. If one process skips its call, the others wait for it forever — the collective hangs.

34.4 Domain Decomposition and Halo Exchange

Now we assemble the pieces into the pattern that distributes a stencil code — and it is the same pattern whether the code is your heat solver, a global weather model, or a galaxy simulation. The idea is domain decomposition: cut the physical domain into pieces, give one piece to each process, and let each process compute its own piece while exchanging just enough boundary data to keep the stencil fed.

Definition (domain decomposition). Domain decomposition is the strategy of dividing a simulation's spatial domain into subdomains, assigning one subdomain to each process, and having each process compute the update for its own subdomain. It is the dominant way to parallelise grid-based physics on distributed memory. Because a finite-difference stencil couples only nearby points, a subdomain's update needs data from other subdomains only along its shared boundaries — so the communication is small (boundaries) relative to the computation (interiors), which is exactly why the approach scales.

Recall the geometry. Your plate is a 2D array; its five-point stencil (Chapter 24) computes each interior cell from its four nearest neighbours. Cut the plate into horizontal strips, one per process — rank 0 takes the top strip, rank 1 the next, and so on down. Each process now owns a block of rows and can update almost all of them from data it already holds. The only problem is the cells along each cut: to update its topmost owned row, a process needs the row just above it, which lives on the neighbour above; to update its bottom row, it needs the row just below, on the neighbour below. Those two borrowed rows are the whole of the communication.

The standard device for holding borrowed boundary data is the ghost cell.

Definition (ghost cell). A ghost cell (or halo cell) is a grid cell a process stores but does not own — a copy of a neighbouring subdomain's boundary data, kept so the stencil can reach "across" the subdomain edge without special-casing it. A process surrounds its owned block with a one-cell-thick layer of ghost cells (a halo), fills that layer each step with the neighbour's current edge values, and then runs the ordinary interior update over its owned cells — the stencil reads the ghosts exactly as if they were real neighbours. The act of refreshing the halo each step is the halo exchange (or halo swap).

Here is the layout for one process's strip. We store it as u(nx, 0:nloc+1): the first index runs across the full plate width (nx columns), and the second index is the process's local row number, with owned rows 1..nloc bracketed by a ghost row at 0 (above) and nloc+1 (below):

   one process's strip, stored as u(nx, 0:nloc+1)
   +---------------------------------------------+
   | u(:,0)      GHOST row (copy of neighbour ↑) |  <- filled by halo exchange
   +---------------------------------------------+
   | u(:,1)      owned row 1        \            |
   | u(:,2)      owned row 2         |  the rows  |
   |   ...                          |  I update   |
   | u(:,nloc)   owned row nloc     /            |
   +---------------------------------------------+
   | u(:,nloc+1) GHOST row (copy of neighbour ↓) |  <- filled by halo exchange
   +---------------------------------------------+

🚪 Threshold Concept — a process computes as if it owned the whole grid. The power of the ghost-cell idea is that it localises the parallelism. Once the halo is filled, the interior update on a process is byte-for-byte the same serial stencil from Chapter 24 — it reads u(i,k-1) and u(i,k+1) without ever knowing or caring that some of those rows are ghosts copied from another machine. All the distributed-memory complexity is quarantined into one routine, the halo exchange, that runs before each update; the physics kernel never changes. This is why the same solver survives every parallel model in this part: coarrays, OpenMP, and MPI all reduce to "get the neighbour's edge data into place, then run the ordinary update." Once you see a parallel stencil code as serial-update-plus-halo-exchange, every such code — including the ones running the world's weather — reads as a variation on this one structure.

Why store the strip with the split (row) direction as the second index, and the full width as the first? Because of column-major order (Chapter 5). In Fortran the first index varies fastest through memory, so a whole row in this layout — u(:, k), the thing we exchange — is a block of nx contiguous reals. That means each halo message is a plain contiguous buffer we can hand straight to mpi_sendrecv with count = nx, no derived MPI datatypes, no packing. This is arrays as Fortran's superpower in a distributed setting: laying the field out to respect the memory order makes the communication trivial. (In C, which is row-major, you would split the other way, into column strips, for the same reason — the memory layout decides the decomposition.)

Now the exchange itself. Each process has an "up" neighbour (rank − 1) and a "down" neighbour (rank + 1), except the top and bottom strips, which run into the physical edge of the plate. MPI has a beautiful device for the edges: MPI_PROC_NULL, a null rank. A send to or receive from MPI_PROC_NULL is a no-op that returns immediately — so if we set the missing neighbour to MPI_PROC_NULL, the same exchange code runs on every process with no if for the boundaries, and the physical-edge ghost rows (holding the fixed Dirichlet boundary) are simply never overwritten:

subroutine exchange_halos(u, nx, nloc, up, down, comm)
  use mpi
  use, intrinsic :: iso_fortran_env, only: dp => real64
  real(dp), intent(inout) :: u(nx, 0:nloc+1)       ! owned rows 1..nloc, ghosts 0 and nloc+1
  integer,  intent(in)    :: nx, nloc, up, down, comm
  integer :: ierr, status(MPI_STATUS_SIZE)

  ! send my top owned row up; receive my top ghost from up   (one tag: source rank disambiguates)
  call mpi_sendrecv(u(:,1),      nx, MPI_DOUBLE_PRECISION, up,   0,   &
                    u(:,0),      nx, MPI_DOUBLE_PRECISION, up,   0,   &
                    comm, status, ierr)
  ! send my bottom owned row down; receive my bottom ghost from down
  call mpi_sendrecv(u(:,nloc),   nx, MPI_DOUBLE_PRECISION, down, 0,   &
                    u(:,nloc+1), nx, MPI_DOUBLE_PRECISION, down, 0,   &
                    comm, status, ierr)
end subroutine exchange_halos

Two mpi_sendrecv calls — one swapping with the neighbour above, one with the neighbour below — and every process's halo is current. We use a single tag (0) for all four messages: because each call names both a specific destination and a specific source rank, the rank already tells each message apart, so no separate tags are needed here. (Tags earn their keep when several messages of different kinds flow between the same pair of ranks — then a distinct tag on each keeps them from being mismatched.) Because mpi_sendrecv cannot deadlock, this is safe for any strip size; because MPI_PROC_NULL neutralises the edges, there is not a single boundary special-case. The neighbours are computed once at setup:

up   = rank - 1;  if (rank == 0)          up   = MPI_PROC_NULL
down = rank + 1;  if (rank == nprocs - 1) down = MPI_PROC_NULL

💡 Intuition — surface versus volume, and why this scales. A strip's computation is proportional to its area (all the cells it owns); its communication is proportional to its perimeter (the edge rows it swaps). As you use more processes on a fixed plate, each strip gets thinner, and the ratio of perimeter to area grows — you spend relatively more time talking and less computing. This surface-to-volume ratio is the fundamental economics of domain decomposition: it scales well as long as each subdomain stays fat enough that its interior work dwarfs its boundary exchange. Cut the plate into too many too-thin strips and communication overtakes computation — precisely the "rising overhead" Case Study 1 of Chapter 31 diagnosed with the Karp–Flatt metric. The cure is 2D decomposition (square tiles have less perimeter per unit area than thin strips), but 1D strips are the right place to learn the pattern, and they scale well to the modest process counts most runs use.

Non-blocking exchange: overlapping communication with computation

mpi_sendrecv is blocking: the process waits inside it until the exchange is done, doing nothing useful meanwhile. On a large run, the network latency of that wait is pure overhead. The remedy is non-blocking communication: start the sends and receives, go compute something that does not depend on the halo (the strip's deep interior), and only wait for the exchange to finish when you finally need the ghost rows (for the edge rows of the update).

Definition (non-blocking communication). A non-blocking MPI operation — mpi_isend, mpi_irecv (the i is for immediate) — initiates a send or receive and returns at once, handing you a request handle instead of waiting for completion. The data is not safe to use until you complete the operation with mpi_wait (or mpi_waitall for several requests). The point is overlap: between the i-call and the wait, the process is free to do independent work, so the communication happens underneath the computation instead of stalling it — often hiding the network cost entirely.

The pattern for the halo, sketched:

integer :: reqs(4), stats(MPI_STATUS_SIZE, 4), ierr

call mpi_irecv(u(:,0),      nx, MPI_DOUBLE_PRECISION, up,   1, comm, reqs(1), ierr) ! ghost from up
call mpi_irecv(u(:,nloc+1), nx, MPI_DOUBLE_PRECISION, down, 0, comm, reqs(2), ierr) ! ghost from down
call mpi_isend(u(:,1),      nx, MPI_DOUBLE_PRECISION, up,   0, comm, reqs(3), ierr) ! my top row up
call mpi_isend(u(:,nloc),   nx, MPI_DOUBLE_PRECISION, down, 1, comm, reqs(4), ierr) ! my bottom row down

call update_deep_interior(u)                 ! rows 2..nloc-1 need no ghosts -- compute NOW, overlapping

call mpi_waitall(4, reqs, stats, ierr)       ! now the ghosts have arrived...
call update_edge_rows(u)                     ! ...update rows 1 and nloc, which needed them

Post the receives first (so the incoming data has somewhere to land the instant it arrives), post the sends, compute the part of the update that needs no ghosts, then mpi_waitall and finish the edges. On a big enough strip the interior work fully hides the exchange, and the halo becomes almost free. The Project Checkpoint uses the simpler mpi_sendrecv; the exercises build this non-blocking version, which is the form production codes use.

⚠️ Common Pitfall — using a non-blocking buffer before mpi_wait. The buffer you hand to mpi_isend or mpi_irecv is off-limits until the matching mpi_wait returns: reading a receive buffer early gives you stale data, and modifying a send buffer early can corrupt the message in flight. The whole point of non-blocking is that the transfer happens later, so touching the data before you have waited is a race. Do independent work between the i-call and the wait; never touch the halo buffers themselves until after mpi_waitall.

🔗 Connection — this is how real models run. The halo exchange you just wrote is, in structure, identical to what happens inside WRF, CESM, and every other production grid code from Chapter 1 — decompose the globe into tiles, give each MPI rank a tile, exchange halos of atmospheric state each timestep. The coarray version you built in Chapter 32 expresses the same decomposition with the same halo, but reaches the neighbour's data through a coindexed reference u(:, :)[q] instead of an explicit message; the physics is identical, only the delivery mechanism differs. Understanding one is most of understanding the other — and understanding this pattern is most of understanding how distributed scientific computing works at all.

🔄 Check Your Understanding. 1. What is a ghost cell, and why does adding a halo let a process run the unmodified serial stencil? 2. Why is MPI_PROC_NULL useful at the top and bottom strips? 3. What does non-blocking communication (mpi_isend/mpi_irecv + mpi_wait) buy you over mpi_sendrecv, and what must you not do before the wait?

Answers

(1) A ghost cell is a stored copy of a neighbour subdomain's boundary data that the process does not own. With a one-cell halo filled each step, the interior update reads u(i,k±1) uniformly — some of those are ghosts, but the stencil neither knows nor cares — so the exact serial Chapter 24 update runs unchanged. (2) A send/recv to MPI_PROC_NULL is a no-op, so setting the missing neighbour of the top and bottom strips to it lets the identical exchange code run everywhere with no boundary if, and leaves the physical-edge ghost rows (holding the fixed Dirichlet values) untouched. (3) Overlap: you can compute the ghost-independent interior while the exchange is in flight, hiding the network cost. You must not read or write the halo buffers passed to mpi_isend/mpi_irecv until the matching mpi_wait/mpi_waitall returns.

34.5 Running MPI Programs: mpif90, mpirun, and a Note on MPI-IO

MPI is a library, so building an MPI program means compiling and linking against it. You could invoke gfortran with the right -I include path and -l link flags, but every MPI installation ships a compiler wrapper that does it for you.

Definition (mpif90). mpif90 (also mpifort) is the MPI Fortran compiler wrapper: it calls your underlying Fortran compiler — usually gfortran — with all the MPI include and library flags added automatically. You use it exactly as you would gfortran, passing the same flags: mpif90 -std=f2018 -Wall -O2 mycode.f90 -o mycode. It is not a different compiler, just gfortran with MPI's paperwork filled in, so everything you know about optimisation flags (Chapter 30) still applies.

To run it, you do not launch the executable directly — you launch it through mpirun (or mpiexec), which starts the requested number of process copies and wires up their communication:

$ mpif90 -std=f2018 -Wall -O2 project-checkpoint.f90 -o heat_mpi
$ mpirun -np 4 ./heat_mpi

-np 4 launches four processes. On a real cluster you would not type mpirun at a shell at all; you would submit a batch script to a scheduler (Slurm's srun, PBS, LSF), which places the processes across the allocated nodes for you — but the program is identical, and mpirun -np N on your laptop is the right way to develop and test it. A subtlety worth stating: -np can exceed your physical core count (the processes time-share), which is how you test an 8-process decomposition on a 4-core laptop — useful, though of course you will not see a speedup that way.

⚡ Performance Note — MPI ranks and physical cores, and the hybrid model. For real performance you want roughly one MPI process per physical core (or per node, in the hybrid model below), and you want ranks that exchange halos placed near each other on the network. The scheduler handles placement, but the principle matters: MPI's cost is the network, so a decomposition that keeps neighbours physically close and messages large-and-few beats one that scatters them and chatters. This is the distributed-memory face of the same truth as cache locality on one core — performance is not accidental; it comes from respecting where the data physically lives, whether "far" means a different cache line or a different rack.

Hybrid MPI + OpenMP. The two models you have learned compose. A modern cluster node has many cores sharing memory; a cluster has many such nodes. The scalable pattern is hybrid: run one MPI process per node (which handles the between-node halo exchange), and inside each process use OpenMP (Chapter 33) to spread the strip's update loop across that node's cores. MPI crosses the nodes; OpenMP fills each node. Your solver is already most of the way there — the OpenMP step from Chapter 33 and the MPI halo exchange from this chapter drop into the same program, MPI around the outside and OpenMP inside. The capstone in Chapter 38 assembles exactly this.

A note on MPI-IO. One serial bottleneck remains: output. If every process sends its strip to rank 0 to write one file, rank 0 becomes a funnel and the write is serial — the 2% that Chapter 31's Amdahl estimate warned would cap your speedup.

Definition (MPI-IO). MPI-IO is the part of the MPI standard for parallel file I/O: it lets all the processes write to (or read from) a single shared file concurrently, each into its own region, so the output does not funnel through one rank. Combined with a parallel file system, it keeps I/O from becoming the serial bottleneck on large runs. In practice you rarely call raw MPI-IO; you use it through a self-describing library — parallel HDF5 or parallel NetCDF (Chapter 25) — which is built on MPI-IO and gives you portable, metadata-rich files as well as parallel writes. For the project, gathering to rank 0 is fine at small scale; MPI-IO (via parallel HDF5) is the answer when it is not.

🔗 Connection — where MPI comes from. MPI is not something you install per project; it is a system library — Open MPI or MPICH — that a cluster provides, and that you can install on a workstation from your package manager. It is one of the pillars of the Fortran-adjacent numerical ecosystem named in Chapter 16, alongside LAPACK, FFTW, and NetCDF/HDF5. The full MPI reference for this book — the calls, datatypes, and the use mpi versus mpi_f08 question — is collected in Appendix G, the parallel-programming reference.

🔧 Modern vs Legacy — use mpi versus mpi_f08. You have three ways to access MPI from Fortran, in ascending modernity. The ancient include 'mpif.h' is a raw text include with no compile-time checking — you will see it in old codes; avoid it in new ones. The use mpi module (which this chapter uses) is the F90-style module: handles are plain integers, ierr is mandatory, and the compiler checks argument counts. The modern use mpi_f08 module (from MPI-3, mirroring Fortran 2008) is the one to prefer in new code: handles become real derived types — type(MPI_Comm), type(MPI_Datatype), type(MPI_Status), type(MPI_Request) — so the compiler catches a communicator used where a datatype was meant, and ierr becomes a genuinely optional final argument. The same skeleton in mpi_f08:

fortran program hello_f08 use mpi_f08 ! modern interface: typed handles, optional ierr implicit none integer :: rank, nprocs call mpi_init() ! ierr is OPTIONAL now -- omit it call mpi_comm_rank(MPI_COMM_WORLD, rank) call mpi_comm_size(MPI_COMM_WORLD, nprocs) print '(a,i0,a,i0)', 'Hello from rank ', rank, ' of ', nprocs call mpi_finalize() end program hello_f08

We teach with use mpi because it is what you will most often encounter — the overwhelming majority of existing Fortran MPI code, and most tutorials, use it — and because the ierr convention it forces on you is worth seeing plainly once. For code you write fresh, reach for mpi_f08. Modern Fortran is a modern language, and its MPI binding has kept pace.

Project Checkpoint

This is the checkpoint the whole part has been aiming at: your solver becomes a distributed program that runs across MPI processes, each owning a strip of the plate and exchanging halos each step. The physics does not change — the five-point stencil and the FTCS update are exactly Chapter 24's — and, crucially, neither does the answer: a distributed run computes the identical numbers a single-process run would, because the halo exchange simply reconstructs, by message, the neighbour rows a serial program has for free.

We build a self-contained MPI heat solver. Each process stores its strip as u(nx, 0:nloc+1) with ghost rows at 0 and nloc+1; the global top edge (rank 0's upper ghost) is held at 100, the bottom edge (rank nprocs−1's lower ghost) and the two side columns at 0 — Dirichlet everywhere, matching the plate of Chapter 1. Each step is: exchange halos, then run the ordinary interior update. Here is the core (the full program, with setup and a gather-to-rank-0 for output, is in code/project-checkpoint.f90):

do step = 1, nsteps
  call exchange_halos(u, nx, nloc, up, down, MPI_COMM_WORLD)     ! refresh ghost rows (§34.4)
  u_new = u
  do k = 1, nloc                          ! owned rows
    do i = 2, nx-1                        ! interior columns (column-major: i inner)
      u_new(i,k) = u(i,k) + r*( u(i-1,k) + u(i+1,k)   &          ! east + west
                              + u(i,k-1) + u(i,k+1)   &          ! north + south (may be ghosts)
                              - 4.0_dp*u(i,k) )                  ! - 4*centre  (identical to Ch.24)
    end do
  end do
  u(:,1:nloc) = u_new(:,1:nloc)           ! commit owned rows; ghosts refreshed next step
end do

The update loop is Chapter 24's stencil verbatim; the only addition is the exchange_halos call that precedes it. That is the entire distributed-memory cost, and it is quarantined in one routine.

A hand-computed run. Take a plate 5 columns wide and 6 rows tall (top row 100, all other edges 0), $\alpha = 1$, $\Delta x = \Delta y = 1$, so $r = \alpha\,\Delta t/h^2 = 0.2$ with $\Delta t = 0.2$ (safely under the CFL limit $r \le 1/4$). Run on 2 processes, so each owns 2 interior rows: rank 0 owns global rows 2–3, rank 1 owns global rows 4–5. Marching three steps, the assembled interior (gathered to rank 0 for printing) is:

after step 1        after step 2        after step 3
 100 100 100 100 100  100 100 100 100 100  100  100  100  100  100
   0  20  20  20   0    0  28  32  28   0    0 32.8 38.4 32.8   0     <- rank 0
   0   0   0   0   0    0   4   4   4   0    0  7.2  8.8  7.2   0     <- rank 0
   0   0   0   0   0    0   0   0   0   0    0  0.8  0.8  0.8   0     <- rank 1  (halo!)
   0   0   0   0   0    0   0   0   0   0    0   0    0    0    0     <- rank 1
   0   0   0   0   0    0   0   0   0   0    0   0    0    0    0

Trace the moment the halo earns its keep. Through steps 1 and 2, heat has spread only within rank 0's rows; rank 1's rows are still all zero, and — note — rank 0's step-2 interior, 28 32 28 and 4 4 4, is the same pattern Chapter 24 computed for its plate, now produced by two cooperating processes. At the start of step 3, exchange_halos sends rank 0's bottom owned row (global row 3, now 0 4 4 4 0) down into rank 1's upper ghost. Rank 1's update of global row 4 then reads that ghost: for its centre column, $0 + 0.2\,(4 + 0 + 0 + 0 - 0) = 0.8$. That 0.8 exists only because the halo carried rank 0's row across the process boundary — without the exchange, rank 1 would see a zero ghost and its row would stay cold, silently computing the wrong physics. The full program also calls mpi_allreduce with MPI_MAX each step on the maximum cell change, giving every process a global convergence measure — a second, honest use of a collective.

Honesty note. The temperatures above are exact and deterministic: run on 1, 2, or 3 processes, this solver prints the identical field, because domain decomposition changes how the work is divided, not what is computed. What is not deterministic, and what this book never presents as measured, is timing and the interleaving of process output — those depend on the machine and the run. The result is physics; the schedule is illustrative.

This is the distributed heat solver. Chapter 35 offers the optional GPU variant, and the Chapter 38 capstone assembles the fully hybrid MPI+OpenMP solver, validates it, and presents it as a paper. Your solver now runs on a cluster.

Summary

This chapter took the running solver off a single node and onto a distributed-memory cluster with MPI — the model that runs essentially all large-scale scientific computing.

The two things to memorize. First, the mpi_send argument order — buffer, count, datatype, destination, tag, communicator, ierr — and that mpi_recv inserts a status before ierr; getting this wrong is the most common MPI bug, and the deadlock cure is mpi_sendrecv. Second, the shape of a distributed stencil code: decompose into strips, wrap each in ghost cells, exchange the halo each step, then run the ordinary serial update — the pattern that scales a finite-difference solver from your laptop to a supercomputer, unchanged in its physics.

Spaced Review

Retrieval practice on the two parallel models this chapter joins — coarrays (Chapter 32) and OpenMP (Chapter 33). Answer before expanding.

1. (Chapter 32.) Coarrays and MPI both distribute the plate with the same 1D decomposition and the same halo. In one sentence each, how does a process reach its neighbour's edge row in each model, and which is part of the Fortran language itself?
   

   Answer

   In **coarrays**, a process (image) reaches the neighbour's data directly through a *coindexed reference* — e.g. `u(:, :)[q]` — and the compiler/runtime turns that into a local read or a network transfer; coarrays are part of the Fortran *standard* (since Fortran 2008). In **MPI**, the process sends and receives the edge row with an explicit *message* (`mpi_sendrecv`); MPI is an external *library*, not part of the language. Same decomposition and same halo, different delivery mechanism — implicit coindexing versus explicit message.

2. (Chapter 33.) OpenMP parallelised the update loop with a directive and no change to the data layout, whereas MPI required you to split the array across processes and add ghost cells. What is the one property of shared versus distributed memory that forces that difference?
   

   Answer

   **Shared memory.** Under OpenMP all threads see the *same* single copy of the field in one address space, so parallelising the loop needs no data partitioning — every thread already has every cell. MPI processes have *private* memory and cannot see each other's arrays at all, so the field must be physically split, and any cell a process needs but does not own must be *copied in* as a ghost and refreshed by message. Distributed memory is what forces decomposition and halos; shared memory does not.

3. (Chapters 32 & 33.) A colleague asks why they should not just use OpenMP for everything, since it is so much less code than MPI. Give the one-sentence hardware reason MPI exists, and name the hybrid scheme that uses both.
   

   Answer

   OpenMP is confined to the shared memory of a *single node* — it cannot use more cores or more RAM than one machine has — so when a problem needs many nodes of a cluster, you must cross to distributed memory, which is MPI's domain. The scalable combination is the **hybrid MPI+OpenMP** model: one MPI process per node handles the between-node halo exchange, and OpenMP threads fill each node's cores with the update loop.

4. (Chapter 33.) In OpenMP you feared a data race when two threads wrote the same cell; the MPI solver in this chapter has no such fear inside a process. Why does the distributed-memory model make data races on the field impossible — and what replaces that worry?
   

   Answer

   Because each process has *private* memory and updates only its *own* owned rows, no two processes ever write the same cell — there is nothing shared to race over. What replaces the data-race worry is the *communication* correctness worry: you must exchange halos at the right time (before each update) and match every send with a receive, or a process computes from stale or missing neighbour data. Shared memory trades communication for race conditions; distributed memory trades race conditions for communication.

What's Next

Your solver now spans a cluster, but there is one more class of hardware it has not touched: the GPU, a processor with thousands of tiny cores built to do the same operation on enormous quantities of data at once — which is exactly what your stencil sweep is. Chapter 35 closes Part VIII by offloading the update to a graphics processor with OpenACC directives and, for the adventurous, CUDA Fortran — the last parallel model in the book, and the one that packs the most arithmetic into the least silicon. Same solver, same stencil, a different kind of chip. Let's put a GPU to work.