Case Study 1: The Nondeterministic Sum

"The first rule of a data race is that it works perfectly on your machine — until the one run that matters."

Executive Summary

A colleague hands you a parallel Fortran routine that computes a vector dot product — the innermost kernel of half the numerical methods in this book. It ran fine in testing. In production, on a busier machine with more threads, it began returning a different answer every run, none of them right. This case study is a complete diagnosis: we reproduce the bug, locate the data race precisely in the read–add–write of a shared accumulator, fix it with a reduction, and prove the fix by showing the result is now identical on every thread count. Along the way we meet the one honest subtlety — floating-point reductions are not perfectly reproducible across thread counts for general real data — and we port a NumPy reduction to parallel Fortran to see where each language wins. By the end you can look at any parallel accumulation and know, before you run it, whether it is correct.

Skills applied: reading and reasoning about a parallel region (§33.1); recognizing a data race in an unscoped shared accumulator (§33.3); the reduction clause and why it is race-free (§33.3); default(none) as a diagnostic (§33.3); the determinism split — result vs. schedule (§33.3, Project Checkpoint); floating-point non-associativity (Chapter 20).

Background

A dot product $\mathbf{a} \cdot \mathbf{b} = \sum_i a_i b_i$ is the workhorse reduction of numerical computing: it is inside every matrix multiply, every residual norm, every projection. It is also the textbook shape of a parallel accumulation — many independent products, all summed into one scalar — and therefore the textbook place to get parallelism wrong. Here is the routine as your colleague wrote it, lightly cleaned up. Read it before reading on, and predict what it does.

function dot_broken(a, b) result(d)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  real(dp), intent(in) :: a(:), b(:)
  real(dp) :: d
  integer  :: i
  d = 0.0_dp
  !$omp parallel do shared(a, b, d) private(i)   ! d is SHARED -- the bug
  do i = 1, size(a)
    d = d + a(i)*b(i)
  end do
  !$omp end parallel do
end function dot_broken

It looks reasonable. The arrays are shared (correct — they are read-only), the loop index is private (correct). The accumulator d is shared so that "all threads add into it." That last decision is the bug, and it is worth seeing exactly why.

Phase 1 — Reproduce the Symptom

The report is "different answer every run." That phrase alone nearly diagnoses it: a deterministic bug (a wrong formula, an off-by-one) gives the same wrong answer every time. An answer that changes between runs of the same binary on the same input is the signature of a data race — a computation whose result depends on the nondeterministic timing of threads.

Set up an input we can check by hand. Take $a_i = i$ and $b_i = 1$ for $i = 1, \dots, 1000$, so the true dot product is $\sum_{i=1}^{1000} i = \frac{1000 \cdot 1001}{2} = 500500$, an integer exactly representable in double precision. The correct routine (built next phase) prints 500500.00 on every run. The broken routine, on several threads, prints a parade of wrong values — perhaps 498211.00, then 500073.00, then, infuriating everyone, 500500.00 by luck, then 499684.00. Those specific numbers are illustrative of the kind of garbage a race produces; the actual values depend on thread timing and are unrepeatable. The single reliable fact is that they disagree with each other and (almost always) with 500500.

The tell: if a number changes when you change nothing but the run, stop looking for a formula error. You have a race, and a race is a scoping mistake.

Phase 2 — Locate the Race

The statement d = d + a(i)*b(i) reads as one action but the processor performs three:

  1. read the current value of d from memory into a register,
  2. add a(i)*b(i) to it in the register,
  3. write the register back to d in memory.

With d shared and many threads executing this for different i at once, two threads can interleave those steps. Suppose d holds 100, and threads A and B both process their iterations at nearly the same moment:

Time Thread A Thread B d in memory
1 read d → 100 100
2 read d → 100 100
3 add 5 → 105 100
4 add 7 → 107 100
5 write 105 105
6 write 107 107

Thread A's contribution of 5 is gone — overwritten by B, which had read the stale 100. The final d is 107, not the correct 112. Every time two threads' three-step sequences overlap like this, one contribution is silently lost, and which ones are lost depends on the exact timing, so the deficit — and the answer — changes each run. That is the whole mechanism. It is not exotic; it is the default behaviour of a shared variable written by many threads.

Why default(none) would have caught it. Had your colleague written default(none), the compiler would still have accepted shared(a, b, d) — sharing d is legal, just wrong. default(none) prevents the accidental race (a variable you forgot to think about); this is a deliberate mis-scoping. The seat belt forces you to make a choice for d; it cannot stop you from choosing badly. The defence against this bug is knowing that a written-across-a-loop shared scalar is a race — which is exactly what this chapter drills.

Phase 3 — Fix It with a Reduction, and Prove the Fix

The correct tool is reduction(+:d). It gives each thread a private copy of d, initialized to 0 (the identity for +), lets each thread accumulate its share with no contention, and combines the private copies into the shared d once, at the end — no interleaving, no lost contributions.

function dot_omp(a, b) result(d)
  use, intrinsic :: iso_fortran_env, only: dp => real64
  real(dp), intent(in) :: a(:), b(:)
  real(dp) :: d
  integer  :: i
  d = 0.0_dp
  !$omp parallel do default(none) shared(a, b) private(i) reduction(+:d)
  do i = 1, size(a)
    d = d + a(i)*b(i)          ! private partial sums, combined at the end
  end do
  !$omp end parallel do
end function dot_omp

Here it is inside a complete, checkable program:

program dot_check
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer, parameter :: n = 1000
  real(dp) :: a(n), b(n), d
  integer  :: i
  do i = 1, n
    a(i) = real(i, dp)         ! a = [1, 2, ..., 1000]
    b(i) = 1.0_dp              ! b = [1, 1, ..., 1]
  end do
  d = dot_omp(a, b)
  print '(a, f12.2)', 'dot = ', d

contains
  function dot_omp(a, b) result(d)
    real(dp), intent(in) :: a(:), b(:)
    real(dp) :: d
    integer  :: i
    d = 0.0_dp
    !$omp parallel do default(none) shared(a, b) private(i) reduction(+:d)
    do i = 1, size(a)
      d = d + a(i)*b(i)
    end do
    !$omp end parallel do
  end function dot_omp
end program dot_check

$ gfortran -std=f2018 -fopenmp -Wall case-study-01-dot.f90 -o dotcheck
$ OMP_NUM_THREADS=1 ./dotcheck
dot =    500500.00
$ OMP_NUM_THREADS=8 ./dotcheck
dot =    500500.00
$ OMP_NUM_THREADS=64 ./dotcheck
dot =    500500.00

The complete program is code/case-study-01-dot.f90.

That is the proof. The value is 500500.00 on 1 thread, on 8, and on 64 — the same every run — because the reduction is correct regardless of how the iterations were divided or which thread finished first. We chose integer-valued inputs so every partial sum is exact and the result is bit-for-bit identical across thread counts. The schedule is nondeterministic; the result is not. This is the determinism discipline the whole chapter insists on, and it is the acceptance test for any parallel reduction: run it on several thread counts and confirm the answer does not move.

Phase 4 — The Honest Subtlety: Floating-Point Reductions

There is one caveat a careful engineer must know, and it is not the race. Floating-point addition is not perfectly associative: $(x + y) + z$ can differ from $x + (y + z)$ in the last bit or two, because each addition rounds (Chapter 20). A parallel reduction adds the elements in a different grouping than a serial loop — and a different grouping again for a different thread count — so for general real data the reduced sum can differ from the serial sum, and from itself across thread counts, in the last bit.

This is a rounding difference, not a race: it is tiny (relative error near machine epsilon, $\sim 10^{-16}$ for double precision), bounded, and not a bug in your code. Our integer-valued test hides it deliberately — those partial sums need no rounding. For real scientific data you should (a) not be surprised when an 8-thread sum differs from a 1-thread sum in the 15th digit, (b) never write a regression test that demands bit-for-bit agreement across thread counts, using a tolerance instead, and (c) reach for a compensated (Kahan) summation if you genuinely need reproducibility. The distinction to hold onto: a race gives you answers wrong in the first digit and different by large amounts; non-associativity gives you answers that agree to about fifteen digits and differ only in the last. One is a catastrophe; the other is the price of floating point.

Phase 5 — Port It: the Same Reduction from Python

Your colleague's original motivation was that Python was too slow. In NumPy the dot product is one call:

import numpy as np
a = np.arange(1, 1001, dtype=np.float64)   # [1, 2, ..., 1000]
b = np.ones(1000, dtype=np.float64)
print(np.dot(a, b))                        # 500500.0

and for this — a single, whole-array reduction — NumPy is genuinely fast, because np.dot drops into compiled, often OpenMP-threaded, BLAS. The Fortran dot_omp matches it and matches its answer. So why bother? Because the moment the reduction is inside an iterative computation — a dot product per timestep, ten thousand times, each depending on the last, as in a conjugate-gradient solver — the pure-Python loop around np.dot falls off the cliff Chapter 1 described, while the Fortran keeps the whole thing compiled and threaded. The lesson of §33.1's Python comparison, made concrete: NumPy wins the single call; Fortran wins the loop.

Discussion Questions

  1. The broken routine passed testing and failed in production "on a busier machine with more threads." Explain how a data race can hide on a lightly loaded machine (few threads, little contention) and surface under load. Why does this make races so dangerous?
  2. default(none) did not, by itself, prevent this bug. Argue both sides: is default(none) therefore worthless here, or does it still earn its place? What class of bug does it prevent?
  3. Phase 4 says never to write a regression test demanding bit-for-bit agreement across thread counts. What should the test assert instead, and how would you choose the tolerance?

Your Turn: Extensions

  • Option A. Instrument the broken dot_broken to reveal the race: have each thread also accumulate into a private reduction(+:) variable, and print the difference between the raced shared d and the correct reduced value. Run it on 2, 4, and 8 threads and tabulate the (varying) deficit.
  • Option B. Replace the + reduction with a max reduction to compute the infinity-norm $\lVert \mathbf{a} \rVert_\infty = \max_i |a_i|$, and verify it against a hand value. What is the identity element the private copies start from, and why must it be $-\infty$ rather than $0$?
  • Option C. Implement a Kahan (compensated) summation inside the reduction loop and compare its result, across thread counts, to the naive reduction on a deliberately ill-conditioned input (a few huge values among many tiny ones). Measure how many digits of reproducibility you buy.

Key Takeaways

  • An answer that changes from run to run on the same input is a data race, not a formula error — and a race is a data-scoping mistake.
  • A shared scalar written across a parallel loop (d = d + …) is the archetypal race: the read–add–write interleaves and contributions are silently lost. The fix is reduction(+:d), which gives each thread a private accumulator and combines them once.
  • Prove a reduction correct by running it on several thread counts and confirming the result does not move; the schedule may vary, the answer may not.
  • Floating-point reductions are not bit-for-bit reproducible across thread counts for general real data — a bounded, ~$10^{-16}$ rounding effect, not a bug. Test with tolerances, not exact equality.