Chapter 32: Coarrays — Fortran's Built-In Parallel Programming Model

"Co-Array Fortran is a small extension to Fortran… for parallel processing." — Robert W. Numrich and John K. Reid, "Co-Array Fortran for Parallel Programming" (1998)

Overview

Chapter 31 ended with a plan and a promise. You measured the solver's serial fraction, estimated its Amdahl ceiling, identified the stencil sweep as the data-parallel hot spot and the time loop as an unbreakable sequential dependency, and chose a model. Now the solver leaves the single core — and it leaves it, in this chapter, through a door that belongs to Fortran and almost no one else.

Every other language in mainstream scientific use reaches parallelism by calling out to a library. You link against MPI, or you sprinkle OpenMP pragmas that a separate runtime interprets, or you import multiprocessing. Fortran can do those things too — the next two chapters do exactly that — but it also has something none of the others has: a parallel programming model built into the language standard itself. Since Fortran 2008 you can take an ordinary array, add one set of square brackets to its declaration, and it becomes a coarray — a variable that exists as a separate copy on each of many parallel images, where any image can reach across and read or write any other image's copy by writing what looks like an array subscript. No #include, no MPI_Init, no handle objects. Parallelism becomes part of the type system.

This is the concrete cash value of two themes this book keeps returning to. Fortran is not dead: it standardized native parallelism in 2008, while other languages were still arguing about how to bolt threads on from the side. And modern Fortran is a modern language: parallel execution is a first-class feature you declare, not a library you import. When you write real :: temperature[*], you are doing something C, C++, and Python cannot do at all without external machinery.

We will build the model from the ground up — images and the SPMD idea, the codimension and remote access, the synchronization that keeps it correct, and the 2018 collectives that make common patterns one line — and then we will spend the payoff where the whole book spends it: on the heat solver. By the end of the chapter your plate is cut into vertical strips, one per image, each image marches the same Chapter 24 stencil over its own strip, and at every step the images exchange their shared edges — a halo exchange — so that the parallel answer is bit-for-bit the serial one.

In this chapter, you will learn to:

• Explain images and the SPMD model, and query them with this_image() and num_images().
• Declare a coarray with a codimension (real :: a[*]) and access a remote image's data with a coindexed reference a[q].
• Place sync all and sync images correctly by reasoning about segments, and use critical and lock/unlock for mutual exclusion.
• Reach for the 2018 collective subroutines — co_sum, co_max, co_broadcast — and know what teams are for.
• Build and run a coarray program two ways, and decide when coarrays beat OpenMP or MPI and when they do not.

Learning Paths

How to read this chapter by track. - ⚡ HPC ("I need parallel code") — this is a core chapter; read every section. §32.2 (remote access), §32.3 (synchronization), and the Project Checkpoint (halo exchange) are the mechanics you will use; §32.5 (coarrays vs MPI) is the judgement call you will make on Monday. - 🔬 Scientist ("my simulation is too slow") — read §32.1–32.3 and the Project Checkpoint closely; that is enough to parallelize a structured-grid code. Skim §32.4's teams and §32.5's Intel/tuning details. - 📖 Standard — read straight through. Coarrays are the part of the language most treatments skip; they are also the part that is uniquely Fortran, and worth meeting properly. - 🔧 Legacy ("I inherited old parallel code") — old Fortran parallel code is almost always MPI (Chapter 34); coarrays are the modern native alternative you may be asked to migrate toward. Read §32.1–32.2 for the model and §32.5 for the honest comparison.

32.1 Coarrays, Images, and the SPMD Model

Start with the execution model, because coarrays make no sense without it. When you run a coarray program, the system does not launch one process; it launches several identical copies of your program, running at the same time, each on its own core (or its own node). These copies are called images.

Definition (image). An image is one of the several concurrent instances of a coarray program that execute together. Each image runs the same executable from the same first line to the same last line, and each has its own complete, private set of the program's variables — its own copy of every array, every scalar, every allocation. Images are numbered $1$ to $N$, where $N$ is the number of images the program was launched with. The image is Fortran's unit of parallelism, the counterpart of an MPI "rank" or an OpenMP "thread" — but, uniquely, it is part of the language, not a library.

Because every image runs the same code, this is called the SPMD model — Single Program, Multiple Data. One program text; many data sets; all running at once. The single program is what you write; the "multiple data" is that each image works on its own share. If all the images ran the identical code on the identical data, they would merely do the same work $N$ times over, which is useless. The art of SPMD programming is having each image discover which image it is and, from that, which slice of the problem is its job. Two intrinsics are how it asks:

• this_image() returns the number of the image that calls it — an integer from $1$ to $N$. On image 3 it returns 3; on image 7 it returns 7. This is how an image learns its own identity.
• num_images() returns the total number of images $N$. Every image gets the same answer.

From those two numbers an image computes its share of the work. If there are 1000 grid columns and 4 images, image me = this_image() takes columns (me-1)*250 + 1 through me*250, and the four images sweep their 250 columns simultaneously. That single idea — identity in, work-share out — is the whole of SPMD.

Here is the smallest program that shows it. Every image announces itself and stores its own contribution in a coarray, which we will define properly in a moment:

program spmd_hello
  implicit none
  integer :: me, ni
  real    :: contribution[*]      ! a scalar COARRAY: one independent copy per image
  me = this_image()
  ni = num_images()
  contribution = real(me)         ! each image sets ITS OWN copy, to its own number
  print '(a,i0,a,i0,a,f4.1)', 'image ', me, ' of ', ni, &
        ': my contribution =', contribution
end program spmd_hello

Compile and run (single-image, for now — the flag and the parallel run come in §32.5):

$ gfortran -fcoarray=single -std=f2018 -Wall example-01-spmd-hello.f90 -o hello && ./hello
image 1 of 1: my contribution = 1.0

Compiled with -fcoarray=single, the program runs as exactly one image, so it prints one line. Run instead under OpenCoarrays on four images (cafrun -n 4 ./hello), and all four images execute the print, each with its own me, producing four lines:

image 1 of 4: my contribution = 1.0
image 2 of 4: my contribution = 2.0
image 3 of 4: my contribution = 3.0
image 4 of 4: my contribution = 4.0

with one crucial caveat you must internalize now.

⚠️ Common Pitfall — the order of output across images is not defined. The four lines above will appear, but in an unpredictable order, and sometimes even interleaved character-by-character, because the images run genuinely concurrently and share one terminal with no referee. Never write a coarray program that depends on which image prints first. When you need ordered or single-source output — as the Project Checkpoint does — you designate one image (conventionally image 1) to do all the printing, after a synchronization guarantees the data it prints is ready. "Let image 1 do the talking" is a habit worth forming immediately.

Now the declaration that made contribution special. The [*] is the whole story.

Definition (coarray). A coarray is a variable declared with a codimension — an extra bound written in square brackets, as in real :: a[*] — which tells the compiler that this variable exists as a separate copy on every image, and that any image is allowed to access any other image's copy. Without the brackets, real :: a is an ordinary local variable: each image still has its own, but no image can reach another's. Add [*] and the copies become collectively addressable: image 1 can read image 5's a. A coarray can be a scalar (real :: a[*]), an array (real :: a(100)[*] — a 100-element array on every image), or allocatable (real, allocatable :: a(:)[:]). The data in brackets () is the local shape, the same on every image; the data in brackets [] is the coshape, describing the grid of images.

The trailing * in [*] means "as many images as the program was launched with" — you do not hard-code the image count, and the same source runs on 1 or 1000 images. The bracket is deliberately chosen to rhyme with an ordinary subscript, because reaching another image's data will look exactly like an array access with one extra index — the subject of §32.2.

🚪 Threshold Concept — parallelism lives in the language, not in a library. Stop and feel how strange and powerful real :: a[*] is. In every other mainstream language, "this variable is distributed across parallel workers" is a fact expressed in library calls, buried in objects and handles that the compiler does not understand. In Fortran it is a fact expressed in the declaration, which the compiler understands completely — it type-checks your remote accesses, it knows the coshape, it can optimize the communication. Once you see parallelism as a property you declare on a variable, the same way you declare its kind or its dimension, the rest of coarray programming stops being exotic and becomes just… more Fortran. Remote data is an array access with an extra subscript. Synchronization is a statement. That is the doorway this chapter walks through, and it is one no other production language offers.

📜 From History: coarrays came from the Cray, and from two people. Coarrays were designed by Robert Numrich and John Reid and first appeared, under the name Co-Array Fortran, as an extension on Cray systems (the T3D and T3E massively parallel machines) in the 1990s, described in their 1998 paper. Their guiding principle was minimalism: add the smallest possible amount of new syntax — essentially just the square bracket — that lets a Fortran programmer write explicit parallel code, and let the compiler and runtime hide whether "another image's data" lives in the same memory or across a network. That extension was adopted, with refinements, into the Fortran 2008 standard, making Fortran the first (and still nearly the only) major language with standardized parallelism. The 2018 standard added the collectives and teams of §32.4. When you write a[q] you are using syntax hammered out on 1990s supercomputers and voted into an ISO standard — the opposite of a dead language.

🔄 Check Your Understanding. 1. What does it mean that a coarray program is "SPMD," and what two intrinsics let an image discover which slice of the work is its own? 2. What is the difference between real :: a and real :: a[*] in a coarray program — how many copies of each exists, and who can reach them? 3. Why must you never rely on the order in which different images print to the screen?

Answers

1. SPMD = Single Program, Multiple Data: every image runs the identical program text, but on its own share of the data. this_image() returns the calling image's number (1..N) and num_images() returns N; from those an image computes which part of the problem it owns. 2. Both give each image its own copy, so there are $N$ copies of each. But a (no brackets) is private — no image can access another's — while a[*] is a coarray, so any image may read or write any other image's copy via a[q]. 3. Because images execute concurrently and share the terminal with no ordering guarantee; the lines can appear in any order or interleave. For deterministic output, one image does the printing after a synchronization.

32.2 Codimensions and Accessing Remote Data

The square bracket in a[*] is a codimension, and it is worth naming carefully, because it is the one genuinely new idea in the syntax.

Definition (codimension). A codimension is a dimension of a coarray that ranges over images rather than over elements within one image. Where an ordinary dimension, written in ( ), indexes the elements of an array on a single image — a(3) is the third element — a codimension, written in [ ], indexes which image's copy you mean: a[3] is image 3's copy of a. A coarray therefore has two kinds of bounds: its ordinary shape (the ( ) part, identical on every image) and its coshape (the [ ] part, describing the arrangement of images). The number of codimensions is the corank; a[*] has corank 1. Just as an array can be multidimensional, a coarray can have several codimensions — grid(nx,ny)[np,*] lays the images out in a 2D np-by-something grid — which we return to below.

The payoff of the codimension is coindexed access — reaching another image's data by writing its image number in brackets:

• a with no brackets means this image's own copy — it is shorthand for a[this_image()]. Almost all of your code operates on the local copy, unbracketed, at full local speed.
• a[q] means image q's copy. Reading x = a[q] fetches image q's value (a "get"); writing a[q] = x stores into image q's copy (a "put"). For an array coarray, a(:,1)[q] is the first column of image q's array — a whole array section pulled from a remote image.

That last form is exactly what a halo exchange needs, and it is why this chapter's solver will read a neighbouring image's edge column with a single coindexed assignment. Let us see remote access in isolation first. Each image computes the square of its own number into a scalar coarray; then image 1 gathers and prints all of them:

program remote_access
  implicit none
  integer :: sq[*]              ! scalar coarray: each image holds one square
  integer :: q
  sq = this_image()**2          ! LOCAL write on every image (its own copy)
  sync all                      ! barrier: no image proceeds until ALL have written
  if (this_image() == 1) then
    do q = 1, num_images()
      print '(a,i0,a,i0)', 'square from image ', q, ' = ', sq[q]   ! coindexed READ of image q
    end do
  end if
end program remote_access

$ gfortran -fcoarray=single -std=f2018 -Wall example-02-remote-access.f90 -o remote && ./remote
square from image 1 = 1

On one image that is the whole output. On four images (cafrun -n 4 ./remote) image 1 reads all four copies and prints, deterministically because image 1 alone does the printing, in loop order:

square from image 1 = 1
square from image 2 = 4
square from image 3 = 9
square from image 4 = 16

The line sq[q] inside image 1's loop is the heart of it: image 1 reaches into image q's private memory and reads its sq. On a shared-memory machine that is a memory load; across a cluster it is a network message — and the source code is identical either way. You wrote a subscript; the runtime decided whether it was a load or a message.

That single line also hides the entire reason §32.3 exists. Look at the sync all between the write and the reads, and ask what would happen without it. Image 1 might reach the loop and read sq[3] before image 3 has executed sq = this_image()**2 — reading uninitialized garbage. The sync all is a barrier: no image moves past it until every image has arrived, which guarantees all the writes are done before any read begins. Remote access without synchronization is a race; the bracket gives you the reach, and the sync makes the reach safe. Hold that thought — it is the next section.

⚡ Performance Note — a coindexed access may be a network round-trip, so touch remote data in bulk. An unbracketed access to your own a is as fast as any local variable. A coindexed access a[q] to another image can be dramatically slower, because on a cluster it crosses the interconnect — latency measured in microseconds, versus nanoseconds for local memory. The practical rule that follows shapes every good coarray code: communicate rarely and in large chunks. Reading a neighbour's whole edge column in one coindexed section assignment, u(:,1) = u(:,n)[q], is one transfer; reading it element by element in a loop, u(i,1) = u(i,n)[q], may be hundreds of tiny messages for the same bytes. This is also why the solver decomposes the plate along columns: in Fortran's column-major layout (Chapter 5, Chapter 27) a column u(:,j) is contiguous in memory, so a whole-column halo is one contiguous block — one efficient message — while a row would be strided and gather-scatter. The memory-order lesson from Part I pays off again, now in communication cost.

🐍 Python Comparison — there is no native equivalent, which is the point. Python has no coarrays. To distribute an array across processes you import mpi4py and call comm.Send/comm.Recv or comm.bcast — explicit library calls, with the array marshalled through the MPI runtime, and a mental model of "processes exchanging messages" you must maintain by hand. There is nothing in the Python language that says "this variable lives on many workers." The nearest data-parallel tools (Dask, the multiprocessing shared arrays) are libraries with their own abstractions layered on top. Fortran's coarray is a language-level answer to the same need, and for structured numerical work it is markedly cleaner: u(:,1) = u(:,n)[q] versus a paired Send/Recv with matching tags and buffer types. This is one of the rare places where Fortran's syntax is not just competitive with Python's ecosystem but simpler than it.

Corank greater than one. When the natural decomposition is two-dimensional — a plate cut into a grid of tiles rather than strips — a corank-2 coarray lays the images out to match. Declare real :: u(nx,ny)[np,*] and the images form an np-by-$(N/np)$ grid; image me's place in that grid is this_image(u), which returns the cosubscripts (a small vector like [2,3]) instead of the plain image number, and image_index(u, [i,j]) converts a grid position back to an image number. We will keep this chapter's solver to a one-dimensional strip decomposition (corank 1) for clarity — it is the honest first step, and the 2D tiling is a natural exercise — but it is worth knowing that the codimension generalizes to match the shape of your problem exactly as an ordinary dimension does.

🔄 Check Your Understanding. 1. In real :: a(50)[*], what does a(10) refer to, and what does a(10)[4] refer to? 2. Why is reading a neighbour's edge as one array-section assignment u(:,1)=u(:,n)[q] usually far better than reading it element-by-element in a loop? 3. What does this_image(u) return for a corank-2 coarray u(...)[np,*], and how does it differ from this_image() with no argument?

Answers

1. a(10) is the 10th element of this image's copy of a (local, fast). a(10)[4] is the 10th element of image 4's copy — a coindexed access that may cross the network. 2. Because a coindexed access can be a network round-trip; one section assignment is a single (contiguous, if it is a column) transfer, while a loop of element accesses can be hundreds of tiny messages for the same data. Communicate rarely and in bulk. 3. this_image(u) returns the calling image's cosubscripts — its position in the coarray's image grid (e.g. [2,3] for a corank-2 coshape). this_image() with no argument returns the single scalar image number (1..N). image_index is the inverse of the cosubscript form.

32.3 Synchronization: Keeping Parallel Correct

Coarrays give you power — any image reaching any image's data — and that power is exactly enough rope to hang your program's correctness. If image A writes a coarray while image B reads it, and nothing orders the two, then B might see the old value, the new value, or (on some hardware) a torn mixture of both. The result depends on timing, so it changes from run to run: a race condition, the signature bug of all parallel programming. Synchronization is how you forbid it, and to place synchronization correctly you must think in segments.

The Fortran standard divides each image's execution into segments, bounded by image control statements (the synchronizations below, plus coarray allocation and a few others). The rule that governs correctness is this: within the ordering the synchronizations establish, if one image defines (writes) a coarray in one segment and another image references (reads) or defines it in an unordered segment, the program is undefined — a race. Synchronization orders segments, turning "who knows what happened first" into "this provably happened before that." You place a synchronization precisely at the boundary where one image's writes must be visible to another's reads.

The workhorse is the global barrier.

Definition (sync all). The statement sync all is a barrier across every image: an image that reaches sync all waits until all images have reached their own sync all, and only then does any image continue. It establishes an ordering point that all images share — everything each image did before its sync all is ordered before everything any image does after it. Consequently, a value written by any image before the barrier is safe to read by any image after the barrier. It is the simplest and most common coarray synchronization, and the right default when in doubt: exchange data, sync all, use it.

That is why the gather in §32.2 was correct: every image wrote sq in the segment before sync all, and image 1 read sq[q] in the segment after — an ordered pair, no race. Drop the sync all and the reads and writes fall in unordered segments, and the program is undefined even if it happens to print the right numbers on your machine today.

sync all is a sledgehammer — it makes every image wait for every other, even images that had nothing to exchange. When two specific images need to coordinate and the rest should keep working, use the scalpel:

integer :: x[*]
if (this_image() == 1) then
  x = 99                     ! image 1 produces a value
  sync images(*)             ! ...and signals every other image that it is ready
else
  sync images(1)             ! each other image waits specifically for image 1
  ! now it is safe to read x[1]  -- ordered after image 1's write
end if

sync images(image-set) synchronizes the executing image only with the listed images (pairwise), leaving all others free to run ahead. It is how you build point-to-point handshakes — a producer signalling consumers, a pipeline stage waiting for its predecessor — without the global cost of sync all. (* in the image set means "all other images.") The trade is that sync images is easier to get subtly wrong: the matching must line up, or images deadlock waiting for a signal that never comes.

Mutual exclusion: critical and lock/unlock. Sometimes the problem is not "wait for a value" but "only one image at a time may touch this." A shared counter incremented by many images, an output file appended to, an accumulator on image 1 — these need exclusion, not a barrier. The critical construct provides it:

critical
  counter[1] = counter[1] + 1     ! only one image executes this block at any instant
end critical

No two images ever execute a given critical block simultaneously, so the read-modify-write of counter[1] cannot be interleaved and corrupted; the increments are serialized and the count is exact. For finer control — several independent locks, or holding a lock across a longer region — Fortran offers lock/unlock with the lock_type from iso_fortran_env:

use, intrinsic :: iso_fortran_env, only: lock_type
type(lock_type) :: lck[*]         ! a lock is itself a coarray
...
lock(lck)                         ! acquire; other images block here until it is free
  ! ... exclusive access to the guarded resource ...
unlock(lck)                       ! release

critical is the simpler tool and the right one for most cases; lock/unlock is there when you need multiple distinct locks or a lock's lifetime to span more than one lexical block. (For the narrow case of a single scalar update, the atomic_add/atomic_ref/atomic_define intrinsics do a hardware-atomic operation even cheaper than a critical, but the construct is clearer to teach and correct here.)

🐛 Find the Bug. A colleague writes a "distributed maximum" by hand: every image puts its local maximum into a coarray, then image 1 reads them all and keeps the largest. It gives the right answer most of the time, but occasionally on a busy cluster it reports a value that is too small. Here is the core: fortran local_max = maxval(my_data) ! each image's own maximum best[1] = local_max ! PUT it into image 1's coarray slot ... but there is only one slot! if (this_image() == 1) result = best[1] Two things are wrong. Name them, and give the correct shape.

Answer

(1) A race with no synchronization, and (2) all images write the same single scalar best[1], clobbering each other. There is no sync all between the writes and image 1's read, so best[1] is read in a segment unordered with the writes — undefined. And even with a sync, every image stores into the one location best[1], so only one write survives (whichever landed last), not the maximum. The fix is to give each image its own slot and reduce: either make best an array indexed by image and take maxval after a sync all, or — far better — delete the hand-rolled code entirely and call the collective co_max(local_max) (§32.4), which does the whole reduction correctly and portably in one line. The bug is a small anthology of why §32.4's collectives exist.

⚠️ Common Pitfall — do not assume images share memory. The most damaging wrong mental model is to imagine all the images poking at one shared array, the way OpenMP threads (Chapter 33) share one address space. They do not. Each image has its own private copy of every coarray; the only way one image's changes reach another is through an explicit coindexed access, and the only way that access is safe is with synchronization. If you write a = a + 1 on every image expecting a single shared a to end up incremented $N$ times, you instead get $N$ separate as each incremented once — because a unbracketed is always the local copy. Coarrays are a distributed-memory model wearing convenient syntax, not a shared-memory one. Respect the privacy of each image's data, and reach across it only on purpose, only with brackets, and only around a synchronization.

🔄 Check Your Understanding. 1. In one sentence, what does sync all guarantee, and why does that make a "write, sync all, read" pattern race-free? 2. When would you reach for sync images instead of sync all? 3. Why does incrementing counter[1] inside a critical block give a correct total across images, while incrementing it without the construct may not?

Answers

1. sync all is a barrier: no image passes it until all images reach it, so everything before it (on any image) is ordered before everything after it (on any image). A value written before the barrier is therefore safe to read after it — the writes and reads fall in ordered segments, which is exactly what forbids the race. 2. When only a specific subset of images needs to coordinate and the others should keep running — a point-to-point handshake or a pipeline — so you avoid the global cost (and over-synchronization) of sync all. 3. critical guarantees only one image executes the block at a time, so the read-modify-write of counter[1] cannot be interleaved; each increment sees the previous one's result. Without it, two images can read the same old value, both add one, and both write back — losing an increment (a lost-update race).

32.4 Collectives and Teams

Some communication patterns are so common that hand-writing them with coindexed accesses and syncs — as the buggy "distributed maximum" above tried to — is both tedious and error-prone. Summing a value across all images, finding a global maximum, broadcasting one image's data to everyone: Fortran 2018 added intrinsic subroutines that do each of these correctly, portably, and in one line.

Definition (collective subroutine). A collective subroutine is an intrinsic that performs a coordinated operation across all images at once, and therefore must be called by every image (with matching arguments) to work. The 2018 collectives are co_sum, co_max, co_min (reductions that combine a value across images), co_broadcast (copy one image's value to all), and co_reduce (a reduction with a user-supplied operation). Because a collective inherently coordinates the images, it carries its own synchronization for the data it touches — you do not wrap a co_sum in sync all. Collectives are the high-level, safe counterpart to raw coindexed access, and they map directly onto MPI's MPI_Allreduce, MPI_Bcast, and friends (Chapter 34).

The reductions overwrite their argument with the combined result on every image:

• call co_sum(x) — after the call, every image's x holds the sum of all images' pre-call x values.
• call co_max(x) / call co_min(x) — every image's x becomes the maximum / minimum across images.
• call co_broadcast(x, source_image=k) — every image's x is set to the value x had on image k.

Optionally, co_sum(x, result_image=1) puts the result on image 1 only (leaving the others' x unchanged) — the exact analogue of MPI's MPI_Reduce versus MPI_Allreduce. Here they are, together with a critical accumulator for contrast:

program collectives
  implicit none
  integer :: me, ni
  integer :: total, biggest, tag
  integer :: counter[*]                    ! coarray used as a shared accumulator on image 1
  me = this_image();  ni = num_images()

  total   = me;  call co_sum(total)        ! -> 1 + 2 + ... + ni   on every image
  biggest = me;  call co_max(biggest)      ! -> ni                 on every image
  tag = 0
  if (me == 1) tag = 42
  call co_broadcast(tag, source_image=1)   ! -> 42 on every image (from image 1)

  counter = 0
  sync all                                 ! establish counter[1] = 0 before anyone increments
  critical
    counter[1] = counter[1] + 1            ! one image at a time; exact count of images
  end critical
  sync all                                 ! all increments done before image 1 reads

  if (me == 1) then
    print '(a,i0)', 'co_sum 1..ni      = ', total
    print '(a,i0)', 'co_max            = ', biggest
    print '(a,i0)', 'co_broadcast tag  = ', tag
    print '(a,i0)', 'critical counter  = ', counter[1]
  end if
end program collectives

$ gfortran -fcoarray=single -std=f2018 -Wall example-03-collectives.f90 -o coll && ./coll
co_sum 1..ni      = 1
co_max            = 1
co_broadcast tag  = 42
critical counter  = 1

On one image the sums are trivial (total = 1, biggest = 1, counter = 1); the broadcast still delivers 42. Run on four images and the collectives combine all four: total = 1+2+3+4 = 10, biggest = 4, counter = 4, and tag is 42 everywhere:

co_sum 1..ni      = 10
co_max            = 4
co_broadcast tag  = 42
critical counter  = 4

Every one of those numbers except tag depends on the image count — which is exactly why the field output of the Project Checkpoint, by contrast, will be the same on any image count: there, the collective work reconstructs the identical physics regardless of how the plate was cut.

Teams: partitioning the images. By default all images work together as one group. Fortran 2018 lets you split them.

Definition (team). A team is a named subset of images that, for the duration of a change team block, behaves as if it were the whole set of images: inside the block, this_image() and num_images() report the image's rank and count within its team, and coindexing a[i] refers to team members. You create teams by having each image choose a team number and calling form team, then execute per-team code inside change team ... end team. Teams let you run different parallel computations on disjoint groups of images at once — the classic use is a coupled climate model where an "atmosphere" team and an "ocean" team run concurrently and exchange data at their interface — or to recurse a decomposition hierarchically.

The shape of it:

use, intrinsic :: iso_fortran_env, only: team_type
type(team_type) :: half
integer :: which
which = merge(1, 2, this_image() <= num_images()/2)   ! lower half -> team 1, upper -> team 2
form team(which, half)
change team(half)
  ! inside here, this_image()/num_images() are RELATIVE to this image's team
  ! ... the two teams run independent parallel work ...
end team

⚠️ A compiler caveat you must not ignore — teams are the least-supported corner of coarray Fortran. Teams are fully standard (Fortran 2018), but they are, in practice, where compiler and runtime support has lagged the standard the most, exactly as parameterized derived types were in Chapter 9. As of the gfortran and OpenCoarrays versions current at this writing, form team / change team support is incomplete or absent — a program using them may fail to compile or link, or behave incorrectly, through no fault of your code. Treat teams as a feature to recognize and understand the model of, not one to build a deadline-critical program on with gfortran today; if you need them now, check the Intel ifx compiler, which has led on 2018 coarray features, and verify on your exact toolchain. The collectives (co_sum and friends), by contrast, are well supported by gfortran+OpenCoarrays and safe to use. We flag the difference plainly rather than pretend the ground is even.

🔗 Connection — the same reduction, three times over. co_sum is not an isolated trick; it is the coarray spelling of a pattern that recurs across this whole part. In Chapter 33 you will write !$omp ... reduction(+:total) to sum across OpenMP threads; in Chapter 34 you will call MPI_Allreduce to sum across MPI ranks. Three models, three syntaxes, one idea: combine a value across all the parallel workers. Learning to see co_sum, reduction(+:), and MPI_Allreduce as the same operation in three dialects is a large part of what makes the next two chapters feel like variations rather than fresh mountains. The Appendix G parallel reference tabulates the correspondences side by side.

🔄 Check Your Understanding. 1. Why must a collective subroutine like co_sum be called by every image, and what happens to the result — where does it end up? 2. What is the difference between co_sum(x) and co_sum(x, result_image=1)? 3. Inside a change team block, what do this_image() and num_images() report — and what is the honest status of team support in gfortran today?

Answers

1. Because it coordinates all images into one combined operation; if some image failed to call it, the others would wait forever (or the result would be ill-defined). By default the combined result is written back to the argument on every image (like MPI_Allreduce). 2. co_sum(x) leaves the sum on every image; co_sum(x, result_image=1) leaves it only on image 1 (the others' x is undefined/unchanged) — the analogue of MPI_Reduce versus MPI_Allreduce. 3. They report the image's rank and count within its team, not globally — the team behaves like a self-contained image set. Honest status: teams are standard but poorly supported by current gfortran/OpenCoarrays; recognize the model, but do not depend on change team compiling on gfortran today.

32.5 Building, Running, and When Coarrays Are the Right Choice

Coarray syntax compiles with gfortran, but how it runs depends on one flag, and there is an important honesty here about what gfortran can and cannot do alone.

Single-image builds — for developing and testing the logic. Compile with -fcoarray=single and gfortran produces an ordinary serial executable in which num_images() is always 1. Every coarray still compiles, every sync all is a no-op, every coindexed access a[1] reads the local copy — the program runs correctly as a degenerate one-image parallel program. This is invaluable: you can write, compile, and debug the correctness of your coarray code on a laptop with nothing installed but gfortran, then run it in parallel unchanged.

$ gfortran -fcoarray=single -std=f2018 -Wall program.f90 -o program && ./program

Multi-image builds — the real thing, via OpenCoarrays. Here is the caveat: gfortran does not, by itself, run a program on more than one image. Multi-image execution requires an external runtime, and the standard one for gfortran is OpenCoarrays, an open-source library that implements the coarray operations (usually on top of MPI). With it installed you compile with the caf wrapper and launch with cafrun:

$ caf -std=f2018 -Wall program.f90 -o program     # caf = gfortran + -fcoarray=lib + OpenCoarrays
$ cafrun -n 4 ./program                            # run on 4 images

caf is a thin wrapper that invokes gfortran with -fcoarray=lib and links the OpenCoarrays library; cafrun -n N launches N images (it wraps mpirun underneath). So the practical situation for gfortran users is: coarray syntax is standard and portable, but multi-image execution on gfortran carries a real dependency on OpenCoarrays. That is a genuine friction compared to OpenMP (built into gfortran with just -fopenmp) and worth weighing. The Intel compilers (ifx/ifort) take a different path — coarrays are built in, enabled with -coarray, and the image count is set by an environment variable (FOR_COARRAY_NUM_IMAGES) rather than a launcher — and Intel has historically led on the newest coarray features (notably teams). Appendix C and Appendix G collect the flags and the setup for each compiler.

A summary of what gfortran supports honestly:

Feature	gfortran status
Coarray syntax, this_image/num_images, coindexing	Supported (any recent gfortran)
Single-image run (-fcoarray=single)	Supported, built in — no external library
Multi-image run (-fcoarray=lib)	Requires OpenCoarrays (caf/cafrun) — a separate install
sync all, sync images, critical, lock/unlock	Supported (with OpenCoarrays for multi-image)
Collectives co_sum/co_max/co_broadcast/co_reduce	Supported with OpenCoarrays
Teams (form team/change team)	Incomplete/absent — verify on your version; prefer ifx

When are coarrays the right choice? Honestly weighed:

• Reach for coarrays when you want native, standardized parallelism with the communication in the language rather than in library calls; when your problem is a structured grid with regular neighbour exchange (the halo pattern of the Project Checkpoint fits coarrays like a glove); when you value source that reads as Fortran rather than as a sequence of MPI_ calls; and when you want one notation that runs on both a multicore laptop and a cluster (the PGAS promise of Chapter 31 §31.3).
• Prefer MPI (Chapter 34) when you need the most mature, most universally installed, most finely tunable option; when you are joining an existing code or community that is already MPI (which is most large HPC codes); or when you need capabilities — sophisticated non-blocking patterns, MPI-IO, the vast tuned-collective ecosystem — that MPI has spent thirty years perfecting.
• Prefer OpenMP (Chapter 33) when your grid fits in one node's memory and you simply want to use its cores with the least ceremony — a directive on a loop, no decomposition, no halos.

The pragmatic truth, which this book will not hide, is that MPI still dominates production HPC and coarrays remain the elegant, standardized, less-universally-deployed alternative. But coarrays are genuinely Fortran's own, they are genuinely cleaner for structured-grid halo exchange, and the gap in tooling narrows every year. For the heat solver — a structured grid whose only communication is a neighbour halo — they are an excellent and instructive fit, which is exactly why we build the parallel solver with them first.

Project Checkpoint

Time to cash it all in. The solver's step has, since Chapter 24, swept the five-point stencil over the whole plate on one core. We now decompose the plate across images — each image owns a vertical strip of columns — and at every timestep the images exchange their shared edge columns as halos with coindexed access, so that each image can apply the identical Chapter 24 stencil to its own strip. The physics does not change one bit; only who computes which columns changes.

Two design choices carry the whole checkpoint. First, we cut the plate along columns, not rows, because a column u(:,j) is contiguous in Fortran's column-major memory (Chapter 5) — so a halo is one contiguous block and one efficient message, per the §32.2 performance note. Second, each image's local array carries two extra halo columns, one on each side, holding copies of its neighbours' edge columns; the interior update then reads those halos exactly where the global stencil would have reached into the neighbouring strip. The neighbour columns arrive by coindexed reads, bracketed by sync all so no image reads a strip that another is still updating:

do step = 1, nsteps
  sync all                                        ! neighbours' u is committed & readable
  if (me > 1)  u(:, 1)      = u(:, nloc+1)[me-1]  ! left halo  <- left neighbour's last owned column
  if (me < ni) u(:, nloc+2) = u(:, 2)[me+1]       ! right halo <- right neighbour's first owned column
  sync all                                        ! all remote reads finished before any overwrite

  u_new(2:nx-1, 2:nloc+1) = u(2:nx-1, 2:nloc+1) + r*(          &   ! IDENTICAL stencil to Ch.24:
        u(1:nx-2, 2:nloc+1) + u(3:nx,   2:nloc+1)             &   ! up + down (all rows are local)
      + u(2:nx-1, 1:nloc)   + u(2:nx-1, 3:nloc+2)             &   ! left + right (halos supply the
      - 4.0_dp*u(2:nx-1, 2:nloc+1) )                              ! cross-image neighbours)
  u(2:nx-1, 2:nloc+1) = u_new(2:nx-1, 2:nloc+1)   ! commit new interior; Dirichlet edges untouched
end do

The two sync all per step are the correctness of the whole thing, and they are worth stating as a rule. The first separates the previous step's commit (a write to u) from this step's halo reads (a remote read of a neighbour's u) — so the value you pull from image me+1 is its committed value from last step, never a half-written one. The second separates the halo reads from this step's commit — so no image overwrites its u while a neighbour is still reading it as a halo. Remove either barrier and neighbouring images fall into unordered segments: a race, and a solver that gives different answers on different runs. This is the segment reasoning of §32.3 made physical.

The full self-contained program — a 5-row plate with 4 interior columns (5 × 6 including the cold Dirichlet edges), the top row held hot at 100, $\alpha = 1$, $\Delta x = \Delta y = 1$, and a CFL-safe $\Delta t = 0.2$ so $r = 0.2 \le \tfrac14$ — is code/project-checkpoint.f90. It decomposes the four interior columns evenly across the images, steps twice, then gathers every strip onto image 1 (more coindexed reads) and prints the whole plate. Because the halo exchange reconstructs exactly the neighbours the serial stencil would have seen, the assembled field is identical for any image count that divides the four interior columns — 1 image (-fcoarray=single), 2 images, or 4 — only the banner changes. Hand-computed, after two steps:

coarray heat solver: 2 image(s), 2 owned column(s) each
  100.0  100.0  100.0  100.0  100.0  100.0
    0.0   28.0   32.0   32.0   28.0    0.0
    0.0    4.0    4.0    4.0    4.0    0.0
    0.0    0.0    0.0    0.0    0.0    0.0
    0.0    0.0    0.0    0.0    0.0    0.0

Trace the cell that proves the halo works. With two images, image 1 owns global columns 2–3 and image 2 owns 4–5. Consider global cell $(2,3)$ — image 1's rightmost owned interior cell — at step 2. Its right neighbour is global column 4, which lives on image 2. After step 1 the whole second row of the interior is 20, so at step 2 cell $(2,3)$ computes $20 + 0.2\,(100 + 0 + 20 + 20 - 4\cdot 20) = 20 + 0.2(60) = 32$ — and the crucial 20 from its right came across the strip boundary, from image 2's column, through the halo u(:, nloc+2) = u(:, 2)[me+1]. Delete the halo exchange and that neighbour reads 0 instead, giving $20 + 0.2(40) = 28$ — the wrong answer. The 32 is the halo doing its job. (The two outer interior cells, global columns 2 and 5, need no halo — their stencils are fully local — which is why a real decomposition communicates only the strip edges, a tiny fraction of the work.)

This is the coarray branch of the parallel solver the whole part is building. Chapter 33 will parallelize the same update with OpenMP threads on one node; Chapter 34 will do this very halo exchange with MPI messages (and formalize domain decomposition and ghost cells); and the Chapter 38 capstone will run the parallel solver, measure its scaling against the Chapter 31 Amdahl estimate, and present it as a paper. The stencil has not changed since Chapter 24, and it never will — that stable step interface is what lets the same physics wear four different parallel clothes.

Summary

This chapter introduced Fortran's native parallel model — coarrays — and used it to make the heat solver run across images with halo exchange.

The two things to memorize. First: a coarray is an ordinary variable with [*] added, giving one copy per image, and a[q] reaches image $q$'s copy — parallelism you declare, not import. Second, and above all: a remote read must be ordered against remote writes by a synchronization — write, sync all, read — or the program has a race and is undefined. Everything correct in this chapter, from the gather to the halo exchange, is an instance of that one discipline.

Spaced Review

Two chapters underlie this one: modules (Chapter 8), which organize any real coarray code, and the parallel groundwork of Chapter 31. Answer before expanding.

1. (Ch. 8.) When you grow the coarray solver into a real program, you will put its field_t type and its step/exchange procedures in a module rather than leaving them loose. Beyond tidiness, what does a module give you for free that a coarray dummy argument specifically requires?
   

   Answer

   An **explicit interface**. Passing coarray dummy arguments (and assumed-shape or polymorphic ones) requires the caller to have an explicit interface for the procedure so the compiler knows the full shape and coshape of the dummies. A module supplies that interface automatically to every `use`-ing unit — exactly the point made for `class`/allocatable arguments in [Chapter 9](../../part-02-modern-fortran-features/chapter-09-derived-types/index.md). Loose external procedures (implicit interface only) cannot portably take coarray arguments.

2. (Ch. 8.) The solver module would declare private and then public :: field_t, step. From Chapter 8, what does that pair achieve, and why is hiding the halo-exchange helper procedure a good idea?
   

   Answer

   `private` makes everything in the module inaccessible from outside by default; `public :: field_t, step` re-exposes only those. The halo-exchange helper stays private — callers use `step` and never touch the internal communication, so you can change *how* halos are exchanged (coarrays now, MPI later) without breaking any caller. Clean interface, hidden implementation — the Chapter 8 discipline, and precisely what lets the same `step` interface wear four parallel clothes.

3. (Ch. 31.) Chapter 31 placed coarrays in both the shared-memory and distributed-memory rows of the parallelism taxonomy. What is the three-word name for the model that lets one notation span both, and how does this chapter's a[q] embody it?
   

   Answer

   The **partitioned global address space** (PGAS) model. `a[q]` is a single notation for "image $q$'s copy of `a`"; the compiler and runtime decide *per access* whether reaching it is a local memory read (shared memory, same node) or a network message (distributed memory, another node). The same coarray source therefore runs on both kinds of hardware — which is why coarrays sit in both rows.

4. (Ch. 31.) From Chapter 31's plan, the solver's stencil sweep is data-parallel but its time loop is a hard sequential dependency. Where do you see both facts in this chapter's Project Checkpoint code?
   

   Answer

   The **data parallelism** is the stencil update itself: each image sweeps its own strip's interior simultaneously with the others — the same operation on different data. The **sequential time-loop dependency** is the `do step` loop with its two `sync all` barriers *inside* each step: the images synchronize and exchange every step and cannot run step $n+1$ before step $n$ completes everywhere. We parallelize *within* a step (across strips) and march the steps in order — exactly Chapter 31's prescription.

What's Next

You have parallelized the solver the way only Fortran can — with the language itself, no external library in the source, images exchanging halos through square brackets. Chapter 33 takes the same stencil update and parallelizes it a second way, with OpenMP: not distributed images with private memory, but shared-memory threads that all see one copy of the plate, unleashed on the update loop by a single directive you write above it. Where coarrays made you think about who owns which columns and when to synchronize, OpenMP asks a different question — which variables are shared and which are private — and answers it with far less ceremony for the common case of one node's cores. Two models, one solver, one stencil; the contrast between them is where the real understanding of parallelism lives.