Chapter 33: OpenMP — Shared-Memory Parallelism with Compiler Directives

"The OpenMP API supports multi-platform shared-memory parallel programming in C/C++ and Fortran." — The OpenMP Architecture Review Board, openmp.org

Overview

In Chapter 31 you learned why the free lunch ended and how much parallelism can buy you; in Chapter 32 you met Fortran's own native parallel model, coarrays, and split the plate across images. This chapter is the workhorse of the whole part. OpenMP is how the overwhelming majority of scientific codes actually use the cores of a single machine, and it does so with a property that will feel almost too good: you parallelize an existing serial loop by adding a comment. A line that begins !$omp is invisible to a compiler that does not understand it — so the same source file compiles and runs correctly as a serial program, and compiles and runs in parallel the moment you add one flag, -fopenmp. No rewrite, no separate parallel version, no external library. You annotate the loop you already have, and the compiler splits its iterations across threads.

That convenience is real, and it is also a trap, which is why this is an advanced chapter. Shared-memory parallelism means every thread can see and touch every variable, and the instant two threads write the same location — or one reads while another writes — you have a data race: a bug whose answer changes from run to run, that vanishes when you add a print statement to hunt it, and that a code can carry silently for years until it corrupts a result nobody double-checks. The single skill this chapter teaches, above every piece of syntax, is data scoping: deciding, for every variable in a parallel region, whether the threads share one copy or each keep their own. Get the scoping right and OpenMP is the gentlest parallelism there is. Get it wrong and the compiler will not warn you; only a wrong number will, if you are lucky enough to notice.

In this chapter, you will learn to:

• Picture parallel execution as fork–join — a serial master thread that forks a team, runs work across it, and joins back to one — and create a team with !$omp parallel.
• Hand the iterations of a loop to that team with the !$omp do work-sharing construct, and control who gets which iterations with a schedule.
• Give every variable a correct data-sharing attribute (shared, private, firstprivate, reduction), and use default(none) so the compiler forces you to decide each one on purpose.
• Combine per-thread partial results safely with reduction, and recognize the race that a missing reduction creates.
• Synchronize with barrier, critical, and atomic; vectorize with !$omp simd; spawn irregular work with !$omp task; and diagnose the silent performance killer, false sharing.
• Parallelize the heat solver's stencil sweep — and get exactly the Chapter 24 answer, on any number of threads.

Learning Paths

How to read this chapter by track. - ⚡ HPC ("I need parallel code") — this is your bread and butter; read every section closely. §33.3 (data scoping) and the Project Checkpoint are the ones you will use daily; §33.4–33.5 (scheduling, false sharing) are where the performance is won or lost. - 🔬 Scientist ("my simulation is too slow") — OpenMP is the fastest path from your serial solver to a parallel one. Read §33.1–33.3 and the Checkpoint carefully; skim §33.5. The reduction in §33.3 is the one construct you cannot do without. - 📖 Standard — read straight through. This is the shared-memory counterpart to the coarrays of Chapter 32 and sets up the distributed MPI of Chapter 34. - 🔧 Legacy ("I inherited parallel code") — old Fortran is full of !$omp directives; §33.1–33.4 teach you to read them, and the ⚠️ pitfalls show you the races their authors may have left behind.

33.1 The Fork–Join Model and the Parallel Region

Everything in OpenMP is built on one picture, and if you hold it clearly the rest is detail. A program starts as a single thread of execution — call it the master (or initial) thread — running your code top to bottom, exactly as every program in this book has so far. When that lone thread reaches a parallel region, it forks: the runtime wakes up a team of additional threads, and every thread in the team — the master included — executes the code inside the region. When they all reach the end of the region, they join: the extra threads go dormant, and the master alone continues past the region, serial again. Fork, run as a team, join, back to serial. That is the whole model.

Definition (OpenMP). OpenMP (Open Multi-Processing) is a standard set of compiler directives, library routines, and environment variables for shared-memory parallel programming in Fortran, C, and C++. You express parallelism by annotating ordinary code with directives — in Fortran, specially formatted comments beginning !$omp — that tell the compiler how to split work across a team of threads that all share one memory. Because the directives are comments, a compiler built without OpenMP support (or invoked without the -fopenmp flag) ignores them entirely and produces a correct serial program from the same source. OpenMP is defined by the OpenMP Architecture Review Board and is supported by essentially every serious Fortran compiler, gfortran included.

Definition (fork–join). The fork–join model is OpenMP's execution model: a single master thread runs serially until it reaches a parallel region, where it forks a team of threads that execute the region concurrently; at the end of the region the threads join — synchronize and terminate — leaving the master to continue alone. A program is thus a sequence of serial stretches punctuated by parallel regions, and you add parallelism incrementally, one region at a time, without restructuring the whole program.

Here is the model as a diagram. Time runs downward; each vertical line is a thread.

      master
        |            <- serial: one thread
        |
   =====|=====  !$omp parallel   (FORK: a team of 4 is created)
    /   |   \  \
   |    |    |   |   <- parallel region: all 4 threads run the block
    \   |   /   /
   =====|=====  !$omp end parallel (JOIN: the team synchronizes and ends)
        |
        |            <- serial again: master alone
        v

The Fortran syntax for a parallel region is a matched pair of directives around the block that the team should run:

!$omp parallel
  ! every thread in the team executes this block
!$omp end parallel

The !$omp at the start of the line is the sentinel that marks an OpenMP directive. To a normal Fortran compiler it is just a comment (it begins with !), so the program is valid serial Fortran. Compile with -fopenmp and gfortran treats those lines as directives instead. This dual life is the source of OpenMP's famous gentleness — and the reason you always show the flag, because without it your "parallel" program is quietly serial.

To do anything useful a thread must know who it is and how many threads there are. Two library functions, made available by use omp_lib, answer exactly those questions:

• omp_get_thread_num() returns this thread's identifier, an integer from 0 to (team size − 1). The master is always thread 0.
• omp_get_num_threads() returns the number of threads in the current team.

Let us watch the fork and join happen. This program prints one line before the region (serial), one line per thread inside it (parallel), and one line after (serial again):

program fork_join
  use omp_lib
  implicit none

  print '(a)', 'before the region: one thread'          ! serial: master only

  !$omp parallel
  print '(a,i0,a,i0)', 'hello from thread ', omp_get_thread_num(), &
                       ' of ', omp_get_num_threads()
  !$omp end parallel

  print '(a)', 'after the region: one thread again'      ! serial: master only
end program fork_join

Compile and run. Note the new flag, -fopenmp; it both enables the directives and links the runtime that provides omp_lib. We set the team size for this run with the OMP_NUM_THREADS environment variable.

$ gfortran -std=f2018 -fopenmp -Wall example-01-fork-join.f90 -o forkjoin
$ OMP_NUM_THREADS=4 ./forkjoin
before the region: one thread
hello from thread 0 of 4
hello from thread 2 of 4
hello from thread 1 of 4
hello from thread 3 of 4
after the region: one thread again

Read that output carefully, because it teaches two things at once — one deterministic, one not. What is deterministic: there are exactly four "hello" lines, one per thread, each reporting a team size of 4 (because we asked for four), and the two serial lines bracket them, first and last. What is not deterministic: the order of the four middle lines. Thread 2 printed before thread 1 here; on the next run it might be 0, 1, 3, 2, or any other order, and two lines can even interleave mid-print, because the threads run genuinely at once and nothing coordinates their access to the screen. If you compile and run this and see a different order than the one printed above, nothing is wrong — you have simply observed the defining property of parallel execution. (The output above is one possible ordering, computed by hand; yours will vary.)

💡 Intuition: think of the master thread as a chef who, at the parallel region, calls in a team of line cooks. Inside the region every cook does the same station's work; omp_get_thread_num() is the cook's number on their apron, omp_get_num_threads() is how many cooks showed up. At end parallel the extra cooks clock out and the chef finishes the plating alone. You do not rebuild the kitchen to add cooks — you just call them in for the busy stretch and send them home after.

Who decides the team size? In order of increasing precedence: a compiled-in default (typically the number of hardware cores the machine reports), the OMP_NUM_THREADS environment variable, a call to the library routine omp_set_num_threads(n), and finally a num_threads(n) clause on the parallel directive itself, which wins for that one region. In practice you set OMP_NUM_THREADS at run time and leave the source alone — the same binary then uses 1 thread on a laptop and 64 on a compute node, no recompile.

⚡ Performance Note — forking is not free. Creating and tearing down a team costs real time: the runtime must wake threads, hand them work, and synchronize them at the join. That overhead is a fixed cost paid per parallel region, on the order of microseconds — negligible if the region does millions of operations, ruinous if you fork a fresh team for a region that does almost nothing. This is the theme performance is not accidental in a parallel dress — good scaling comes from deliberately amortizing overhead, never for free. The practical rule, which the Project Checkpoint and Case Study 2 both lean on: make parallel regions big. Put the region outside the hot loop when you can, so you fork once and reuse the team, rather than forking anew on every iteration. Overhead is why a real measured speedup always trails the Amdahl ceiling of Chapter 31 — the fork/join cost is part of the serial fraction you could not see on paper.

📜 From History. Before OpenMP, every vendor shipped its own incompatible set of parallel directives — Cray had one dialect, SGI another — so a parallel Fortran code was welded to one manufacturer's compiler. In 1997 a group of hardware and compiler vendors formed the OpenMP Architecture Review Board and published a single portable standard, and — tellingly for this book — the very first OpenMP specification was for Fortran; the C/C++ binding followed a year later. It is one more data point for two of this book's themes — Fortran is not dead and modern Fortran is a modern language: when the industry needed to standardize shared-memory parallelism, it standardized it on Fortran first, because Fortran is where the heavy numerical loops live.

🐍 Python Comparison. You might expect Python's threading module to do what OpenMP does — run a loop's iterations on several threads at once. It cannot, for CPU-bound work: CPython's Global Interpreter Lock (the GIL) permits only one thread to execute Python bytecode at a time, so threads take turns rather than running in parallel, and a numerical loop gets no speedup (sometimes a slowdown, from lock contention). Python programmers reach instead for multiprocessing — separate processes with separate memories, which is really the message-passing model of Chapter 34, not shared memory — or they drop the loop into NumPy or compiled Fortran, where genuine threading (often OpenMP underneath) happens out of Python's reach. Fortran has no GIL: !$omp parallel do gives you true, simultaneous, shared-memory parallelism directly, which is exactly why, in the theme this book keeps returning to — Fortran and Python are better together — the hot kernel belongs in Fortran and the orchestration in Python.

🔄 Check Your Understanding. 1. In the fork–join model, how many threads run the code between two parallel regions? 2. What does omp_get_thread_num() return for the master thread, and what is the range of its value on a team of 8? 3. You compile the fork_join program without -fopenmp and run it. How many "hello" lines print, and why?

Answers

1. One — the master thread alone. Parallel regions are islands of many threads in a sea of serial execution; outside them, only the master runs. 2. It returns 0 for the master. On a team of 8, values range over 0, 1, …, 7 (that is, 0 to omp_get_num_threads() - 1). 3. One. Without -fopenmp the !$omp lines are ordinary comments, so there is no team — the master executes the block once. This is the feature, not a bug: the same source is a correct serial program.

33.2 Work Sharing: Splitting a Loop Across the Team

The fork_join program has every thread do the same thing, which is rarely what you want. The point of a team is to divide work — to have thread 0 handle some iterations of a loop, thread 1 handle others, and so on, so the loop finishes in a fraction of the time. Directives that split work among the existing team, rather than replicating it, are called work-sharing constructs, and the one you will use ninety-nine times out of a hundred shares out the iterations of a do loop.

Definition (work sharing). A work-sharing construct divides the execution of a code block among the threads of the team that encounters it, so that the block's work is done once collectively rather than once per thread. OpenMP's work-sharing constructs for Fortran are !$omp do (partition a loop's iterations among the threads), !$omp sections` (give different code blocks to different threads), `!$omp single (exactly one thread runs the block), and !$omp workshare (partition whole-array Fortran statements). A work-sharing construct must appear inside a parallel region — it shares work among a team that already exists — and it carries an implicit barrier at its end: no thread proceeds past it until all have finished their share.

The !$omp do construct sits immediately before a do loop and tells the team to split its iterations:

!$omp parallel                 ! fork the team
!$omp do                       ! share the loop's iterations among the team
do i = 1, n
  ! iteration i is executed by exactly one thread
end do
!$omp end do
!$omp end parallel             ! join

The pattern "a parallel region whose entire body is one work-shared loop" is so common that OpenMP provides a combined directive, !$omp parallel do, which forks the team and shares the loop in one stroke (closed by !$omp end parallel do). It is exactly equivalent to the two nested directives above and is what you will write most often:

!$omp parallel do
do i = 1, n
  ! ...
end do
!$omp end parallel do

Crucially, !$omp do does not run the loop n times per thread — it runs each iteration once, on one thread. The iterations are partitioned. Which thread gets which iterations is decided by the schedule (the full story is §33.4), but the default, static, has a property worth relying on: it divides the iteration range into contiguous, nearly equal blocks, one per thread, and — importantly — it does so deterministically. The same loop, the same thread count, and static scheduling always produce the same iteration-to-thread assignment, run after run. That reproducibility lets us hand-predict the mapping. Here we record, for each iteration, which thread executed it:

program work_sharing
  use omp_lib
  implicit none
  integer, parameter :: n = 8
  integer :: who(n), i

  who = -1
  !$omp parallel do default(none) shared(who) private(i) schedule(static)
  do i = 1, n
    who(i) = omp_get_thread_num()      ! record who ran iteration i
  end do
  !$omp end parallel do

  do i = 1, n
    print '(a,i0,a,i0)', 'iteration ', i, ' ran on thread ', who(i)
  end do
end program work_sharing

$ gfortran -std=f2018 -fopenmp -Wall example-02-work-sharing.f90 -o worksharing
$ OMP_NUM_THREADS=4 ./worksharing
iteration 1 ran on thread 0
iteration 2 ran on thread 0
iteration 3 ran on thread 1
iteration 4 ran on thread 1
iteration 5 ran on thread 2
iteration 6 ran on thread 2
iteration 7 ran on thread 3
iteration 8 ran on thread 3

This output is fully deterministic, and you can compute it by hand: 8 iterations across 4 threads under static gives 2 contiguous iterations each — thread 0 takes {1, 2}, thread 1 takes {3, 4}, thread 2 takes {5, 6}, thread 3 takes {7, 8}. Notice why the final print loop can report a clean result: the parallel loop wrote the who array in parallel (each thread writing its own distinct elements — no conflict), then the region ended, and only afterward, serially, did we read it back in order. Separating the parallel write from the serial read is a habit worth forming early; it is exactly what the heat solver does with its two field buffers.

⚠️ Common Pitfall — a work-sharing directive with no team does nothing. Write !$omp do` (or `!$omp sections) outside any parallel region and the loop simply runs serially on the master — no error, no warning, just no parallelism. The construct shares work among a team, and if there is no team there is nothing to share it among. The combined !$omp parallel do avoids this by creating the team itself, which is one reason to prefer it until you are deliberately reusing one team across several loops.

The other work-sharing construct named in the outline, !$omp sections, is for a different shape of problem: not "many iterations of one loop" but "a few different tasks that can run at once." Each !$omp section block is given to one thread:

!$omp parallel
!$omp sections
  !$omp section
    call compute_energy(field)       ! one thread does this
  !$omp section
    call compute_flux(field)         ! another thread does this, concurrently
!$omp end sections
!$omp end parallel

This is task parallelism (Chapter 31) — different operations at once — and it scales only as far as you have distinct sections: two sections keep at most two threads busy, no matter how large the team. It is genuinely useful for overlapping a handful of independent phases, but the heavy lifting in scientific computing is data parallelism — the same operation over a huge array — which is the province of !$omp do. This is one more payoff of the theme arrays are Fortran's superpower: the whole-array operations of Chapter 5 already state the independent, data-parallel work that !$omp do then hands to a team. For the modern, more flexible form of task parallelism (recursive and irregular work), OpenMP later added the task construct, which we meet in §33.5.

🔗 Connection. !$omp workshare deserves a mention because it speaks Fortran's native dialect. Recall from Chapter 5 that Fortran expresses the stencil as a whole-array statement, lap(2:n-1,2:n-1) = u(1:n-2,…) + …. You cannot put !$omp do on that — there is no explicit loop — but !$omp workshare parallelizes exactly such array assignments:

fortran !$omp parallel !$omp workshare c = a + b ! the whole-array add is split across the team !$omp end workshare !$omp end parallel

It is elegant and matches Fortran's array syntax, but in practice workshare is unevenly optimized across compilers, and most performance-minded codes write the explicit !$omp do loop instead, for the fine control over scheduling that §33.4 describes. We will use the explicit loop in the solver.

33.3 Data Scoping: Shared, Private, Firstprivate, and Reduction

We now arrive at the heart of the chapter, and of shared-memory programming generally. When a team of threads runs a region, some variables should be held in common — all threads reading and writing one shared copy — and others should be duplicated, each thread getting its own private instance the others cannot see. Deciding which is which, for every variable in the region, is data scoping, and it is where correctness is won or lost.

Definition (data scoping). Data scoping is the assignment of a data-sharing attribute to each variable referenced in a parallel region, declaring whether the team shares one instance of it or each thread holds a private one. The principal attributes are shared (one instance, visible to all threads), private (each thread gets its own uninitialized instance, invisible to the others and gone at the region's end), firstprivate (private, but each thread's copy is initialized to the value the variable had before the region), and reduction (private per thread during the region, then combined into one shared result at the end). Getting an attribute wrong does not produce a compiler error — it produces a data race or silently wrong values, which is why scoping demands deliberate care.

The rules of thumb are short and you should memorize them:

• Read-only data is shared. If threads only read a variable — the grid dimensions, a coefficient, the input field — one shared copy is correct and efficient. There is no conflict in many threads reading the same value.
• A variable each thread writes as scratch is private. Loop indices, temporary accumulators computed fresh inside each iteration, a local real used and reused — each thread needs its own, or they will clobber one another. Loop iteration variables are the canonical case: they must be private, and the loop variable of the !$omp do loop is made private for you automatically. Inner loop variables and other temporaries are not automatic — you must scope them yourself.
• A per-thread result to be combined is a reduction. A running sum, a maximum, a count — see below.
• A private variable that must start from the pre-region value is firstprivate.

Because forgetting a variable is so easy and so dangerous, OpenMP gives you a seat belt: the default(none) clause. It removes all implicit scoping and requires you to name an attribute for every variable the region uses, or the code will not compile. Use it always. The compile error it produces is a gift — it is the compiler forcing you to think about a variable you would otherwise have raced on.

🚪 Threshold Concept — every variable is either shared or private, and you must know which. Serial programming lets you be careless about a variable's lifetime and ownership; it works either way because there is only one thread. Parallel programming makes ownership a correctness property. Once you internalize that every name in a parallel region has an attribute — that the compiler will, by default, guess shared for you, and that its guess is a data race waiting to happen for anything you write — you stop writing directives by pattern-matching and start writing them by reasoning. The reflex default(none) installs is the mark of someone who has crossed this threshold: they never let the compiler guess, because they have learned that its guess is exactly the bug. Every parallel bug in this chapter, and most in your career, is a scoping mistake in disguise.

The reduction: combining partial results correctly

Consider the most ordinary computation imaginable — summing an array:

s = 0.0_dp
do i = 1, n
  s = s + x(i)
end do

Every iteration reads s, adds to it, and writes it back. Now try to parallelize it naively, sharing s across the team so all threads accumulate into it. Disaster: s = s + x(i) is not one indivisible action but three — read s, add x(i), write s — and two threads can interleave those steps. Thread A reads s as 10; thread B reads s as 10 before A has written; A writes 13, B writes 12, overwriting A's contribution. x(i) from thread A is simply lost, and which contributions vanish depends on the exact timing, so the answer is both wrong and different every run. This is the data race made concrete, and it is the single most common OpenMP bug in existence.

The fix is not to share s and not to make it plainly private (a plain private s would give each thread its own sum that is then thrown away). The fix is a reduction.

Definition (reduction). A reduction(op:var) clause tells OpenMP that var accumulates a combination of per-thread contributions under the associative operator op (+, *, max, min, .and., .or., and others). The runtime gives each thread its own private copy of var, initialized to the operator's identity element (0 for +, 1 for *, -huge for max, …); each thread accumulates into its private copy with no contention; and at the end of the region the runtime combines all the private copies into the single shared var using op. You get the parallelism of private accumulation and the correctness of a single combined result, with no race and no explicit synchronization.

Here is the sum done right, with reduction(+:s), alongside a hand-checkable input — the integers 1 through 100 as reals, whose sum is $\frac{100 \cdot 101}{2} = 5050$:

program reduction_demo
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer,  parameter :: n = 100
  real(dp) :: x(n), s
  integer  :: i

  do i = 1, n
    x(i) = real(i, dp)                 ! x = [1.0, 2.0, ..., 100.0]
  end do

  s = 0.0_dp
  !$omp parallel do default(none) shared(x) private(i) reduction(+:s)
  do i = 1, n
    s = s + x(i)                       ! each thread sums its share into a private s;
  end do                               ! the runtime adds the private sums at the end
  !$omp end parallel do

  print '(a, f10.2)', 'sum = ', s
end program reduction_demo

$ gfortran -std=f2018 -fopenmp -Wall example-03-reduction.f90 -o reduction
$ OMP_NUM_THREADS=4 ./reduction
sum =    5050.00

The printed value is 5050.00, and — this is the whole point — it is 5050.00 on 1 thread, on 4, on 64, and every single run, because the reduction combines the partial sums correctly no matter how the iterations were divided or in what order the threads finished. The schedule is nondeterministic; the result is not. (We chose small whole-number inputs so that every partial sum is exact in floating point and the answer is bit-for-bit identical across thread counts; for general real data, different reduction orders can differ in the last bit or two of rounding — a real but usually negligible effect, and a reason Chapter 20's floating-point care matters here too.)

⚠️ Common Pitfall — the shared accumulator race (the classic wake-up call). Take the program above and make one "innocent" change: drop the reduction(+:s) clause and declare s as shared instead.

fortran !$omp parallel do default(none) shared(x, s) private(i) ! WRONG: s is raced do i = 1, n s = s + x(i) ! read-add-write, unsynchronized end do !$omp end parallel do

This compiles without complaint and runs without crashing. But because every thread reads-adds-writes the one shared s without coordination, contributions are lost to the interleaving described above, and you get a wrong answer that changes from run to run — perhaps 4128.00, then 4873.00, then, maddeningly, 5050.00 by luck, then 4519.00. (Those specific numbers are illustrative of the kind of garbage you get; the actual values depend on timing and are unrepeatable — that is the nature of the bug.) The lesson every OpenMP programmer learns once, and never forgets: an unsynchronized shared variable that is written in a parallel loop is a race, and a race is not a rare crash — it is a quietly wrong number. When you find yourself writing to one shared scalar across a loop, stop: you almost certainly want a reduction.

firstprivate: private, but with a running start

A plain private variable enters the region uninitialized — each thread's copy is undefined until the thread assigns it. Usually that is what you want for scratch. But sometimes each thread needs a private copy that begins at the value the variable held before the region — a coefficient computed once, say, that each thread will then modify locally. That is firstprivate:

base = 100.0_dp
!$omp parallel do default(none) shared(y) firstprivate(base) private(i)
do i = 1, size(y)
  y(i) = base + real(i, dp)        ! each thread's 'base' starts at 100.0, as intended
end do
!$omp end parallel do

Had base been plain private, each thread's copy would be undefined and y would fill with garbage; had it been shared, it would be fine here (nobody writes it) but fragile the moment a thread modified it. The symmetric clause lastprivate copies the value from the last iteration back out to the shared variable after the region — useful, occasionally, when you need the final iteration's scratch value to survive.

🐛 Find the Bug. A colleague parallelizes a nested loop that normalizes each row of a matrix, and gets different results every run:

fortran !$omp parallel do default(none) shared(a, n) private(i) ! only i is private do i = 1, n rowsum = 0.0_dp do j = 1, n rowsum = rowsum + a(i,j) end do do j = 1, n a(i,j) = a(i,j) / rowsum end do end do !$omp end parallel do

The code will not even compile with default(none) — which is the point. What is missing, and why is the missing piece a race rather than merely an omission?

Answer

Two variables are used in the region but not scoped: the inner loop index j and the per-row scratch rowsum. Under default(none) the compiler refuses to guess and errors until you scope them, which is precisely the seat belt working. Both must be private: each thread processes different rows i, but they would share one j and one rowsum if those were left shared/default, so thread A's rowsum for its row would be corrupted by thread B accumulating its own row into the same location — a classic race giving run-to-run garbage. The fix is private(i, j, rowsum). Note the deeper lesson: the outer index i is automatically private as the !$omp do variable, but the inner index j and every scratch temporary are your responsibility. Forgetting an inner-loop temporary is the most common scoping race there is.

🔄 Check Your Understanding. 1. You sum the elements of an array across a team. Which data-sharing attribute does the running total need, and what goes wrong if you make it shared instead? 2. What is the difference between private(x) and firstprivate(x)? 3. Why is default(none) recommended on every parallel region, even though it means more typing?

Answers

1. It needs reduction(+:total). Made shared, the unsynchronized read-add-write of total is a data race: contributions are lost to interleaving and the answer is wrong and different every run. 2. private(x) gives each thread its own uninitialized copy; firstprivate(x) gives each thread its own copy initialized to x's value from before the region. 3. Because it forces you to assign an attribute to every variable explicitly, turning a silent, dangerous default (shared, which races on anything you write) into a compile error you must resolve on purpose. It converts the most common class of parallel bug into a message from the compiler.

33.4 Synchronization and Scheduling

Data scoping keeps threads from corrupting each other's variables. Two more concerns remain: sometimes threads genuinely must coordinate — wait for one another, or take turns at a shared resource — and the runtime must decide how to hand loop iterations out. These are synchronization and scheduling.

Synchronization: barrier, critical, atomic

The bluntest coordination is the barrier: a line no thread passes until every thread has reached it. Work-sharing constructs (!$omp do`, `!$omp sections) and the end of a parallel region carry an implicit barrier, so you rarely write one by hand, but you can — !$omp barrier — when one phase of a region must complete across all threads before the next begins.

When threads must take turns at something that genuinely cannot be done concurrently — appending to a shared list, writing a log line, updating one shared counter in a way no reduction covers — you need mutual exclusion. OpenMP offers two grades of it:

• !$omp critical wraps a block that only one thread may execute at a time. Any thread reaching it waits until no other thread is inside, then enters; the others queue. It is fully general — any code can go in a critical section — but it is relatively expensive (it manages a lock) and it serializes that block, so overusing it throws away your parallelism.
• !$omp atomic applies to a single update statement of the form x = x op expr and guarantees just that one read-modify-write happens indivisibly. It is far cheaper than critical because it maps to a single hardware atomic instruction, but it is restricted to one simple update.

!$omp parallel do default(none) shared(hits, x) private(i)
do i = 1, n
  if (inside_circle(x(i))) then
    !$omp atomic
    hits = hits + 1              ! indivisible increment: no two threads collide
  end if
end do
!$omp end parallel do

That counter could also be a reduction(+:hits), and where a reduction applies it is faster still (no per-update synchronization at all — threads accumulate privately and combine once at the end). Reach for atomic or critical only when the update does not fit a reduction's fixed set of operators, or when the thing you are protecting is not a simple accumulation.

⚡ Performance Note — the ranking. For combining per-thread results, prefer, in order: reduction (no contention; combine once at the end) → atomic (one cheap hardware-atomic update) → critical (a general but comparatively slow lock). A critical section inside a hot loop is often a sign you have reintroduced serialization: if every iteration must pass single-file through a critical block, Amdahl's serial fraction (Chapter 31) has crept back in through the side door, and your speedup collapses. Measure, and push the exclusion out of the inner loop wherever you can.

Scheduling: who gets which iterations, and when

By default OpenMP hands out loop iterations in fixed contiguous blocks, decided before the loop runs. That is fine when every iteration costs the same. But when iterations vary wildly in cost — some grid cells trigger an expensive branch, some pixels of a fractal take a hundred times longer to converge than others — a fixed split leaves some threads idle while one grinds through the heavy iterations. Scheduling is how you control the hand-out.

Definition (scheduling). Scheduling is the policy by which a work-sharing loop assigns its iterations to threads, set with a schedule(kind[, chunk]) clause. The three principal kinds are static — the iteration range is divided into equal chunks once, before the loop runs, and dealt out round-robin to threads; low overhead and reproducible, ideal for uniform work. dynamic — threads grab a chunk of iterations, and when a thread finishes its chunk it comes back for another, at run time until the loop is exhausted; higher overhead but self-balancing, ideal for uneven work. guided — like dynamic, but the chunk size starts large and shrinks as the loop proceeds, trading a little balance for less hand-out overhead. The optional chunk sets the block size (schedule(dynamic,16) deals 16 iterations at a time).

The trade-off is balance versus overhead. static has essentially zero run-time bookkeeping — every thread knows its iterations in advance — so it is fastest when the work is even. dynamic pays for a bit of coordination on every chunk hand-out, but it keeps every thread busy right to the end when the work is lumpy, which more than repays the overhead. A picture:

   Uniform work (each iteration costs the same):
     static  -> [ thr0 ][ thr1 ][ thr2 ][ thr3 ]   all finish together. Best.

   Lumpy work (a few iterations are very expensive, '#'):
     static  -> [ thr0 ][ thr1# ][ thr2 ][ thr3 ]  thr1 stuck on '#',
                 thr0/2/3 idle, waiting at the barrier.  Bad.
     dynamic -> threads grab small chunks; whoever finishes early
                 grabs more, so the '#' work is shared out.  Good.

The heat-solver stencil is uniform — every interior cell costs exactly the same handful of operations — so its right schedule is static, and that is what the Project Checkpoint uses. A load-imbalanced problem, where dynamic earns its keep, is the subject of Case Study 2.

🔗 Connection — !$omp do vs do concurrent. You met another way to say "these iterations are independent" in Chapter 29: do concurrent, the standard Fortran construct that asserts a loop's iterations may run in any order. The two are complementary. do concurrent is portable, standard Fortran with no directives, and a modern compiler may parallelize or vectorize it — but the standard only lets you state independence, not control the team, the schedule, or the reductions; you are at the compiler's mercy for whether and how it goes parallel. !$omp parallel do is not standard Fortran (it is OpenMP), but it gives you explicit, portable-across-compilers control: the thread count, the schedule, the exact scoping, the reduction. For a production kernel where you must guarantee parallel execution and tune it, OpenMP remains the workhorse; do concurrent is the cleaner expression of intent, increasingly able to target GPUs (Chapter 35) as compilers mature.

33.5 SIMD, Tasks, and the Pitfall of False Sharing

Three final topics round out a working knowledge of OpenMP: a second, finer level of parallelism (SIMD), a construct for irregular work (tasks), and the performance trap that catches everyone once (false sharing).

SIMD: parallelism inside a single thread

Threading is not the only parallelism in a modern core. Each core has SIMD units — Single Instruction, Multiple Data — that apply one operation to several array elements at once, the vectorization you met in Chapters 27 and 29. Threading and SIMD are orthogonal and multiply: a loop can be split across 8 threads, each of which processes 4 elements per instruction, for 32-fold throughput. The !$omp simd directive asks the compiler to vectorize a loop, and !$omp do simd does both at once — share iterations across threads and vectorize each thread's chunk:

!$omp parallel do simd default(none) shared(y, a, x) private(i)
do i = 1, n
  y(i) = a * x(i) + y(i)         ! threaded across the team AND vectorized per thread
end do
!$omp end parallel do simd

Where the compiler would auto-vectorize anyway (Chapter 29), !$omp simd mostly confirms your intent and can push vectorization through cases the compiler would otherwise play safe on. It is the explicit lever for the SIMD parallelism Fortran's array style already courts.

Tasks: parallelism for irregular work

!$omp do needs a countable loop with a known trip count. Some problems are not shaped that way — walking a tree, following a linked list, a recursive divide-and-conquer — and for those OpenMP 3.0 added the task construct. A thread reaching !$omp task packages the block as a unit of work and hands it to the runtime, which farms queued tasks out to idle threads; !$omp taskwait waits for the spawned tasks to finish. Tasks are typically generated by one thread (inside !$omp single) and executed by the whole team:

!$omp parallel
!$omp single                      ! one thread generates the tasks...
  !$omp task
    call process(left_subtree)    ! ...and the team executes them
  !$omp end task
  !$omp task
    call process(right_subtree)
  !$omp end task
  !$omp taskwait
!$omp end single
!$omp end parallel

Tasks are how you parallelize the irregular, pointer-chasing computations that a data-parallel do loop cannot express. For the regular, array-shaped work of a stencil solver you will not need them — but it is good to know the tool exists when the computation stops being a rectangle.

False sharing: the invisible performance killer

Now the trap. You have scoped everything correctly — no races, correct answers — and yet your parallel loop is slower than serial, or barely faster on 8 cores. The culprit is very often false sharing, and it is invisible in the source because it is a fact about the hardware, not the code.

Definition (false sharing). False sharing occurs when threads on different cores write to different variables that happen to occupy the same cache line — the fixed-size block (commonly 64 bytes) that is the unit of transfer between memory and a core's cache. Although the threads share no data logically, the cache-coherence hardware sees writes to one line from multiple cores and must ping-pong that line between their caches to keep them consistent, serializing what should be independent work. The result is a severe, silent slowdown with no incorrect answer to signal the problem — the code is right, just mysteriously slow.

The textbook trigger is a per-thread accumulator array, one slot per thread, meant to avoid a reduction:

real(dp) :: partial(0:nthreads-1)          ! DANGER: adjacent slots share a cache line
partial = 0.0_dp
!$omp parallel default(none) shared(partial, x, n) private(i, tid)
tid = omp_get_thread_num()
!$omp do
do i = 1, n
  partial(tid) = partial(tid) + x(i)       ! each thread hammers its own slot...
end do
!$omp end do
!$omp end parallel
! ...but partial(0), partial(1), ... sit side by side in ONE cache line,
!    so the cores fight over that line on every update. Slow, though correct.

Eight real(dp) values are 64 bytes — one cache line — so all eight threads' slots live on the same line, and every update by any thread invalidates the line in all the others' caches. The computation is correct and crawls. The fixes are exactly the ones this chapter has been building toward: use a reduction (which gives each thread a truly private accumulator, not a shared-array slot, sidestepping the whole problem), or if you must keep the array, pad each thread's data onto its own cache line. Ninety-nine times in a hundred, the answer is "you should have written reduction(+:s)," which is why we taught it first.

⚡ Performance Note. False sharing is why the humble reduction clause is not just a convenience but usually the fastest way to combine per-thread results: it keeps each thread's partial in a private register or private storage with no shared cache line to contend over, then combines once. Hand-rolling a per-thread array to "avoid the reduction overhead" is a classic way to make your code slower and more fragile at the same time. When Chapter 28's profiler shows a parallel loop scaling far worse than Amdahl predicts and you have already ruled out a critical bottleneck, suspect false sharing on a shared array that threads write by index.

🔄 Check Your Understanding. 1. Threading and SIMD are said to "multiply." What does each contribute, and why are they independent? 2. When would you reach for !$omp task` instead of `!$omp do? 3. Your parallel loop gives the right answer but runs slower than the serial version. Name the likely hardware cause and the one-clause fix.

Answers

1. Threading splits iterations across cores; SIMD makes each core process several elements per instruction. They are independent because they exploit different hardware (multiple cores vs. a core's vector unit), so their speedups multiply — 8 threads × 4-wide SIMD ≈ 32×. 2. When the work is irregular — recursion, tree/list traversal, an unknown or dynamic trip count — that a countable do loop cannot express. !$omp task queues units of work for the team. 3. False sharing — threads writing different variables that share a cache line, forcing the coherence hardware to ping-pong the line between cores. The usual fix is to replace the shared per-thread array with a reduction clause (a genuinely private accumulator).

Project Checkpoint

The moment the whole part has been building to: your solver leaves the single core. In Chapter 24 the interior update became a real five-point stencil sweep; Chapter 31 established, on paper, that this sweep is the parallel fraction (the time loop around it is inherently serial — step $n+1$ needs step $n$) and that the sweep is embarrassingly parallel: within one step, every interior cell's new value depends only on the old values of its neighbours, so no cell's update depends on any other's. That independence is not luck — it is the two-buffer structure Chapter 24 built in — and it is exactly what makes OpenMP a two-line change.

We keep the frozen step(field, alpha, dt) signature from Chapter 6; only the body changes, gaining an !$omp parallel do over the interior. The scoping is the whole lesson: the fields are shared (all threads read the same old field%u and write distinct cells of u_new), the loop indices are private (each thread needs its own i and j), and the read-only coefficients are shared. We write it with an explicit loop rather than array sections precisely so the scoping is visible:

subroutine step(field, alpha, dt)                 ! frozen signature (Ch. 6); field_t (Ch. 9)
  type(field_t), intent(inout) :: field
  real(dp),      intent(in)    :: alpha, dt
  real(dp) :: u_new(field%nx, field%ny)
  real(dp) :: rx, ry
  integer  :: i, j, nx, ny
  nx = field%nx;  ny = field%ny
  rx = alpha*dt / field%dx**2                      ! diffusion numbers (Ch. 24)
  ry = alpha*dt / field%dy**2
  u_new = field%u                                  ! copy keeps Dirichlet edges; interior overwritten
  !$omp parallel do default(none)                        &
  !$omp   shared(field, u_new, nx, ny, rx, ry)           &
  !$omp   private(i, j) schedule(static)
  do j = 2, ny-1                                    ! outer index j: auto-private, listed for clarity
    do i = 2, nx-1                                  ! inner index i: MUST be private -- else a race
      u_new(i,j) = field%u(i,j)                                             &
                 + rx*(field%u(i-1,j) - 2.0_dp*field%u(i,j) + field%u(i+1,j)) &
                 + ry*(field%u(i,j-1) - 2.0_dp*field%u(i,j) + field%u(i,j+1))
    end do
  end do
  !$omp end parallel do
  field%u = u_new                                  ! commit the whole new field at once
end subroutine step

Three points of correctness earn their place. First, the two buffers are mandatory now, not merely tidy: if threads wrote back into field%u in place, a thread computing column j might read field%u(i, j-1) after the thread owning column j-1 had overwritten it — a race on the read. Reading the untouched old field%u and writing a separate u_new is what makes the sweep independent, so the two-buffer discipline of Chapter 24 is exactly the property that licenses the parallelism. Second, i must be private. It is the inner loop index, not the !$omp do variable, so it is not scoped automatically; leave it shared and every thread stamps on the same i, producing a race and nonsense. This is the single most common mistake in parallelizing a nested loop, and default(none) catches it at compile time. Third, static scheduling is right here because every cell costs the same — uniform work, so the cheapest, reproducible schedule wins.

The full self-contained program is code/project-checkpoint.f90: the same $5 \times 5$ plate as Chapter 24, $\Delta x = \Delta y = 1$, top edge held at $100°$, $\alpha = 1$, and a CFL-safe $\Delta t = 0.2$ (so $r = 0.2 \le \tfrac14$), marched two steps. Because the result is deterministic, its hand-computed output is identical to Chapter 24's, to the digit:

after step 1                            after step 2
  100.00  100.00  100.00  100.00 100.00   100.00  100.00  100.00  100.00 100.00
    0.00   20.00   20.00   20.00   0.00     0.00   28.00   32.00   28.00   0.00
    0.00    0.00    0.00    0.00   0.00     0.00    4.00    4.00    4.00   0.00
    0.00    0.00    0.00    0.00   0.00     0.00    0.00    0.00    0.00   0.00
    0.00    0.00    0.00    0.00   0.00     0.00    0.00    0.00    0.00   0.00

Trace the centre cell $(2,3)$ to confirm nothing changed but the speed. Step 1: neighbours are the initial field — $100$ above, three $0$s — so $0 + 0.2(100 - 0) + 0.2(0) = 20$. Step 2: it now sees $100$ above, $0$ below, and two warm $20$s beside it: $20 + 0.2(100 + 0 - 2\cdot20) + 0.2(20 + 20 - 2\cdot20) = 20 + 0.2(60) + 0.2(0) = 32$. Identical to Chapter 24. This is the honesty at the core of parallel programming: the answer is deterministic and must match the serial solver exactly — run it on 1 thread or 64 and you get this same grid — while the thread schedule (which core computes which column, in what order they finish) is nondeterministic and invisible in the result. If your parallel solver ever disagrees with the serial one, you have a scoping bug, not a faster answer.

An illustrative scaling note (not a measurement). We do not time code in this book, so treat this as expectation, not result. From Chapter 31's Amdahl estimate — the sweep at ~98% of run time — the ideal ceiling is about $1/(1-0.98) = 50\times$, with roughly $7\times$ on 8 threads. On a genuinely large grid (say $2000 \times 2000$), a correctly scoped OpenMP stencil should approach that on a single node; on the tiny $5\times5$ plate above it would run slower in parallel, because the fork/join overhead per step dwarfs the twenty-odd operations of real work — a concrete instance of the "make regions big" rule. To close the gap toward Amdahl on a real run you would hoist the parallel region outside the time loop (fork once, not once per step) and confirm the schedule and cache behaviour — precisely the optimization worked in Case Study 2. The same step interface now hides a serial body, an OpenMP body, and — in Chapters 32 and 34 — coarray and MPI bodies, all assembled in the Chapter 38 capstone.

Summary

This chapter turned the serial stencil sweep into a shared-memory parallel one by annotating it with OpenMP directives — and taught the data scoping that keeps that annotation correct.

The two things to memorize. First, default(none) and the scoping reflex: in a parallel region every variable is shared or private, read-only data is shared, and scratch and loop indices are private — let the compiler force you to say which. Second, the shared-accumulator race and its cure: s = s + x(i) across a shared s is the archetypal OpenMP bug (wrong, nondeterministic, silent), and reduction(+:s) is the fix — which also, by giving each thread a truly private accumulator, avoids the false sharing that a hand-rolled per-thread array would suffer.

Spaced Review

Retrieval practice on the two chapters this one stands on: arrays (Chapter 5), whose column-major layout and sections shape how we parallelize, and the parallel foundations of Chapter 31, whose Amdahl arithmetic tells us how much to expect. Answer before expanding.

1. (Ch. 5) The stencil sweep is do j …; do i …. Recalling column-major order, which index belongs on the inner loop for cache-friendly access — and note that this is also the index that must be private in the OpenMP version. Why both?
   

   Answer

   The **first** index, `i`, on the inner loop: Fortran stores arrays column-major, so `u(i,j)` and `u(i+1,j)` are adjacent in memory and the inner sweep over `i` walks contiguous addresses, using each cache line fully. Separately, `i` must be **`private`** in `!$omp parallel do` because it is the *inner* loop index (not the auto-private `!$omp do` variable `j`), so each thread needs its own copy; shared, it would race. Same variable, two different reasons for care — memory layout and thread ownership.

2. (Ch. 5) In Chapter 5 you could write the whole interior Laplacian as one array-section statement. Which OpenMP work-sharing construct parallelizes such a whole-array statement directly, and why do we nonetheless use an explicit !$omp do loop in the solver?
   

   Answer

   **`!$omp workshare`** parallelizes whole-array Fortran statements. We prefer the explicit `!$omp do` loop because `workshare` is unevenly optimized across compilers and, more importantly, gives no control over the **schedule** or fine scoping, whereas an explicit loop lets us choose `static` scheduling and see every variable's attribute — the control a production kernel wants.

3. (Ch. 5) Why can all interior cells of one time step be updated in parallel with no synchronization within the step, and what property of the update makes that true?
   

   Answer

   Because each new value depends only on the **old** neighbour values (the two-buffer FTCS structure): the sweep reads the unchanging old field and writes a separate new buffer, so no cell's update depends on another cell's *new* value. The interior updates are therefore independent — *embarrassingly parallel* — needing no locks or ordering, only correct shared/private scoping.

4. (Ch. 31) Chapter 31 estimated the solver at ~98% parallel. State Amdahl's ceiling that implies, and explain why the tiny $5\times5$ plate would actually run slower with OpenMP than serially.
   

   Answer

   Ceiling $= 1/(1-0.98) = 50\times$. The $5\times5$ plate slows down because its interior is nine cells — a few dozen operations per step — while forking and joining a team costs on the order of microseconds *per step*; the parallel overhead dwarfs the work. Amdahl's ideal assumes overhead-free parallelism; on a trivially small problem the overhead *is* the runtime, so real speedup needs a big enough grid (and, better, the parallel region hoisted outside the time loop).

5. (Ch. 31) Distinguish the data parallelism of the stencil sweep from task parallelism, and name the OpenMP construct that expresses each.
   

   Answer

   **Data parallelism** applies the *same* operation (the five-point update) to *many* independent data elements (every interior cell) — expressed by **`!$omp do`** (or `parallel do`). **Task parallelism** runs *different* operations concurrently (e.g. compute energy while computing flux) — expressed by **`!$omp sections`** or, for irregular work, **`!$omp task`**. The solver's heavy work is data parallelism, which is why it scales.

What's Next

OpenMP took the solver across the cores of one machine, sharing its memory. But one node's cores and RAM are a ceiling: to reach the thousands of cores of a cluster you must cross to distributed memory, where each process owns private memory and they cooperate by passing messages. Chapter 34 introduces MPI, the Message Passing Interface that runs the world's largest simulations: you will cut the plate into tiles, give each process a tile, and teach the tiles to exchange their shared edges — the halo exchange — each step. Where OpenMP shared one field across threads, MPI will have many processes each holding a piece of the field and trading boundaries at the seams. The directives become explicit calls; the shared memory becomes a network. Let's distribute the solver.