Case Study 2: Scaling the Stencil

"Correctness first, then speed — but a parallel program that is slower than the serial one is neither."

Executive Summary

The Project Checkpoint parallelized the heat solver's step correctly: fields shared, indices private, and the right answer on every thread count. Correct is not the same as fast. This case study takes that correct solver and optimizes it, the way you would to make a real scaling study honest. We find the overhead hiding in plain sight — a fresh team forked on every single time step — and hoist the parallel region outside the time loop so the team is forked once. We justify static scheduling for the uniform stencil, then add a parallel convergence check and use it to walk straight into the classic false-sharing trap and out the other side with a reduction. Finally we reason — honestly, without timing anything — about what strong scaling this solver can and cannot achieve, against the Amdahl ceiling from Chapter 31. The result is the same numbers as Chapter 24, computed in a way that could actually approach its speedup ceiling on a real machine.

Skills applied: parallel-region overhead and "make regions big" (§33.1); !$omp parallel with an inner !$omp do`, `!$omp single, and barriers (§33.2, §33.4); choosing a schedule for uniform work (§33.4); reduction(max:) for a convergence check (§33.3); diagnosing and fixing false sharing (§33.5); strong scaling vs. Amdahl (Chapter 31).

Background

The checkpoint's step wraps its interior sweep in !$omp parallel do, which forks a team and joins it every time it is called. A real run calls step once per time step — thousands, often millions, of times. So the checkpoint pays the fork/join cost on every step. On a large grid that cost is negligible against the sweep; on a modest grid, or at high step counts, it is exactly the overhead that makes measured speedup fall short of Amdahl's ideal. The fix is structural, and it is the single most important OpenMP optimization for an iterative solver: fork the team once, around the whole time loop, and work-share the sweep inside it.

Phase 1 — Find the Overhead

Recall the shape of the checkpoint's time loop (driver side) and step:

do s = 1, nsteps
  call step(plate, alpha, dt)      ! each call forks AND joins a team
end do

Each call step runs !$omp parallel do … !$omp end parallel do — one fork, one join. Forking a team is not free: the runtime must wake threads, distribute the loop, and synchronize them at the join, on the order of a microsecond of pure overhead per region. Put numbers to it, illustratively (Tier-3 round numbers, for the reasoning): if the stencil sweep on your grid takes 5 microseconds and the fork/join costs 1 microsecond, you are spending $1/6 \approx 17\%$ of every step on team management — a serial-fraction tax, invisible on paper, that caps your speedup well below the Amdahl estimate. Over a million steps that is a million needless forks.

The principle: parallel-region overhead is a fixed cost per region. Amortize it by making regions big and few. An iterative solver has one obvious big region: the whole time loop.

Phase 2 — Hoist the Parallel Region

We restructure so the team is forked once, outside the time loop, and the sweep is a work-sharing !$omp do inside it. The buffer commit — copying the new interior back — must happen exactly once per step by a single thread, which !$omp single provides; its implicit barrier guarantees every thread sees the committed field before the next step's sweep reads it.

subroutine run(field, alpha, dt, nsteps)
  use kinds,      only: dp
  use heat_types, only: field_t
  type(field_t), intent(inout) :: field
  real(dp),      intent(in)    :: alpha, dt
  integer,       intent(in)    :: nsteps
  real(dp) :: u_new(field%nx, field%ny)
  real(dp) :: rx, ry
  integer  :: i, j, s, nx, ny
  nx = field%nx;  ny = field%ny
  rx = alpha*dt / field%dx**2
  ry = alpha*dt / field%dy**2
  u_new = field%u                                  ! seed boundaries once

  !$omp parallel default(none)                            &
  !$omp   shared(field, u_new, nx, ny, rx, ry, nsteps)    &
  !$omp   private(i, j, s)
  do s = 1, nsteps                                 ! every thread runs the step loop, in lockstep
    !$omp do schedule(static)                      ! ...but the SWEEP is shared out
    do j = 2, ny-1
      do i = 2, nx-1
        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 do                                   ! implicit barrier: sweep complete
    !$omp single
    field%u(2:nx-1,2:ny-1) = u_new(2:nx-1,2:ny-1)  ! commit interior, one thread
    !$omp end single                               ! implicit barrier: field visible to all
  end do
  !$omp end parallel                               ! ONE join, after all steps
end subroutine run

The team is forked at !$omp parallel` and joined once at `!$omp end parallel, no matter how many steps run. Inside, each step's sweep is shared across the team by !$omp do; the two implicit barriers (end of do, end of single) keep the threads marching in step and the memory coherent. Note the scoping subtlety the hoist introduces: s, the step index, is now private — every thread runs the step loop for itself, kept synchronized by the barriers, so each needs its own s.

Run it on the same $5 \times 5$ plate as Chapter 24 for two steps and print the final field:

$ gfortran -std=f2018 -fopenmp -Wall case-study-02-hoisted.f90 -o hoisted
$ OMP_NUM_THREADS=4 ./hoisted
final field after 2 steps:
  100.00  100.00  100.00  100.00  100.00
    0.00   28.00   32.00   28.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

Identical to the checkpoint and to Chapter 24 — cell $(2,3)$ is $32$, cell $(2,2)$ is $28$, cell $(3,2)$ is $4$ — because we changed only when the team is forked, not the arithmetic. The optimization is invisible in the answer and decisive in the speed, which is exactly the property you want from a performance change: it must not move the result. The complete program is code/case-study-02-hoisted.f90.

⚠️ The hazard the hoist introduces. Forking once means the serial parts of each step — the buffer commit, and in a fuller solver the boundary re-imposition and any I/O — now live inside the parallel region and must be guarded. Forget the !$omp single around the commit and every thread would copy the whole interior, a redundant-write race on field%u. The checkpoint's per-step fork got this for free (the copy was serial, between calls); the hoisted version must name it single. Trading fork overhead for synchronization discipline is the essence of this optimization.

Phase 3 — Choose the Schedule

Which schedule belongs on that !$omp do? The stencil is uniform: every interior cell costs the same handful of operations, independent of its value or position. For uniform work, static is optimal — it assigns contiguous column blocks once, with essentially zero run-time bookkeeping, and (a bonus) reproducibly. dynamic would only add overhead here: its per-chunk hand-out coordination buys load balancing the uniform stencil does not need, so it would run slower for no benefit. The rule is worth stating crisply:

There is a case where the solver's work goes non-uniform — an adaptive or masked update that skips converged cells, or a domain with variable material properties triggering different branches — and there dynamic earns its keep. But the plain diffusion stencil is the textbook uniform loop, and static is its textbook schedule.

Phase 4 — A Convergence Check, and the False-Sharing Trap

A steady-state solver wants to stop when the plate stops changing. Add a per-step convergence measure: the largest change any interior cell underwent this step, $\max_{i,j}|u^{\text{new}}_{i,j} - u^{\text{old}}_{i,j}|$. When it drops below a tolerance, the plate has settled. This is a max reduction. On our two-step run the step-2 changes are $8$ at the corners of the warm band, $12$ at the centre $(2,3)$ (from $20$ to $32$), and $4$ along row 3, so the maximum change is $12$:

real(dp) :: dmax
dmax = 0.0_dp
!$omp parallel do default(none) shared(field, u_new, nx, ny) private(i, j) reduction(max:dmax)
do j = 2, ny-1
  do i = 2, nx-1
    dmax = max(dmax, abs(u_new(i,j) - field%u(i,j)))
  end do
end do
!$omp end parallel do
! step 2: dmax = 12.0  (the centre cell moved 20 -> 32)

Now the trap. A programmer worried about "reduction overhead" might try to hand-roll it: give each thread its own slot in an array and take the max at the end.

real(dp) :: slot(0:nthreads-1)          ! DANGER: false sharing
slot = 0.0_dp
!$omp parallel default(none) shared(slot, field, u_new, nx, ny) private(i, j, tid)
tid = omp_get_thread_num()
!$omp do
do j = 2, ny-1
  do i = 2, nx-1
    slot(tid) = max(slot(tid), abs(u_new(i,j) - field%u(i,j)))   ! each thread hammers its slot
  end do
end do
!$omp end do
!$omp end parallel
dmax = maxval(slot)                     ! correct answer... but slow

This is correct — it computes $12$ — and it can be dramatically slower than the reduction, because of false sharing (§33.5). Eight real(dp) slots are 64 bytes, exactly one cache line, so all threads' slots live on the same line. Every time any thread updates its slot, the cache-coherence hardware must invalidate that line in every other core's cache and shuttle it back, serializing what should be independent work. The threads are not sharing data logically — each touches only its own slot — but they share a cache line, and the hardware cannot tell the difference. The symptom is maddening: a correct parallel loop that scales worse than serial, with nothing wrong in the source to see.

The fix is the tool we already have. reduction(max:dmax) gives each thread a genuinely private accumulator — in a register or private storage, on no shared cache line — and combines them once at the end. It is simpler, faster, and race-free all at once. The reduction is not a convenience you pay for; it is the optimization. Hand-rolling per-thread arrays to "avoid" it is how you buy false sharing.

The diagnostic chain. When a correctly-scoped parallel loop scales badly: profile it → rule out a critical serializing the inner loop → suspect false sharing on any shared array threads write by index → replace the per-thread array with a reduction, or pad each thread's data onto its own cache line.

Phase 5 — Strong Scaling, Honestly

What speedup should this hoisted, correctly-scheduled solver achieve? We do not time code in this book, so this is expectation, grounded in Amdahl, not a measurement. From Chapter 31, the solver is ~98% parallel (the sweep) and ~2% serial (setup, the single commit, periodic output), giving an Amdahl ceiling of $1/(1-0.98) = 50\times$ and about $7\times$ on 8 threads for a large grid where the sweep dominates. Three honest caveats separate that ideal from what a real run would show:

1. Overhead we did not remove. Even hoisted, each step pays two barriers (end-do, end-single) and a serial commit. On a large grid these are tiny; on a small one they dominate — which is why the $5\times5$ plate would still run slower in parallel, its nine interior cells no match for even one barrier per step. Strong scaling needs a big enough grid.
2. Memory bandwidth. The stencil is memory-bound (Chapter 24): it does little arithmetic per value fetched. Several threads on one socket share one path to memory, so past a handful of threads you can saturate bandwidth and stop scaling regardless of core count — a ceiling Amdahl (which counts only the serial fraction) does not model. This is why a real stencil often scales well to one socket and then flattens.
3. False sharing and NUMA. The traps of Phase 4, plus the fact that a large plate's memory may be split across a node's memory banks, mean the measured curve trails the ideal. Closing that gap is the craft the Chapter 38 capstone documents.

So the honest expectation is: on a $2000\times2000$ plate, this solver should scale well across one socket's cores — a solid, high-efficiency single-node speedup — and then flatten as memory bandwidth bites, never quite reaching the overhead-free Amdahl ceiling. That is not a disappointment; it is the true shape of shared-memory stencil scaling, and reporting it accurately — with the thread count, and against the ideal — is the professional honesty the whole book insists on. To go past the single-node wall you must cross to distributed memory: Chapter 34's MPI, many nodes and a halo exchange.

Discussion Questions

1. The hoisted version is faster but harder to get right — the serial commit had to be wrapped in !$omp single. Under what circumstances is the checkpoint's simpler per-step fork the better engineering choice, despite its overhead?
2. Phase 4's hand-rolled per-thread array gives the correct answer yet is a bug worth fixing. In what sense is a correct-but-slow parallel loop a "bug," and how would you catch it (what would the profiler show)?
3. Phase 5 argues the solver will flatten before Amdahl's ceiling because the stencil is memory-bound. If you could not add memory bandwidth, what algorithmic change (hint: do more arithmetic per value fetched) might let the extra cores help — and what does that trade away?

Your Turn: Extensions

• Option A. Add the convergence check to the hoisted run (inside the same parallel region, with a reduction(max:dmax) and a !$omp single to test the tolerance and set a shared converged flag). Make the time loop exit when dmax falls below 1.0e-6_dp. What must be shared, and where is the barrier that makes the flag safe to read?
• Option B. Deliberately induce false sharing: implement the per-thread slot array of Phase 4, then "fix" it by padding each thread's datum to its own 64-byte cache line (e.g. a derived type with the value plus filler). Confirm both give $12$, and explain, from cache-line size, why the padded version would scale and the packed one would not.
• Option C. Write the parallel region so the boundary conditions are re-imposed each step inside the region (imagine time-varying edges). Which construct re-imposes them exactly once per step, and why would doing it in the plain !$omp do body be wrong?

Key Takeaways

• For an iterative solver, hoist the parallel region outside the time loop: fork the team once with !$omp parallel, work-share the sweep with !$omp do` inside, and guard the per-step serial commit with `!$omp single. Same answer, far less fork/join overhead.
• Match the schedule to the work: static for the uniform stencil (cheapest, reproducible); dynamic only when iterations are genuinely uneven.
• A correct parallel loop can still be a performance bug. Hand-rolled per-thread accumulator arrays invite false sharing — different threads' slots on one cache line, ping-ponged by the coherence hardware. The reduction clause is the fix and the optimization.
• Report scaling honestly: a memory-bound stencil scales across a socket and then flattens on bandwidth, trailing the overhead-free Amdahl ceiling. Naming the thread count and the ideal is the discipline; crossing to MPI (Chapter 34) is how you pass the single-node wall.