Case Study 2: Reducing on the GPU — the Plate's Total Energy, Two Ways

"A sum looks trivial until a thousand threads try to write it at once."

Executive Summary

Every simulation needs diagnostics: a total, a norm, a maximum, computed by reducing a whole field to a single number. On a CPU that is a one-line sum. On a GPU it is the subtlest common operation there is, because a reduction asks thousands of threads to combine their results into one place — and doing that correctly, let alone quickly, is a genuine design problem. In this study you will build the same reduction — the plate's total thermal energy and its $L_2$ norm — three ways: the naive version that looks right and silently corrupts its answer, the OpenACC reduction clause that gets it right in one line, and the explicit CUDA Fortran kernel that shows you what the clause is doing under the hood (and what it costs to do by hand). By the end you will know why, for reductions above all, the directive path earns its keep — and you will have met an honesty this book insists on: a floating-point reduction on the GPU may not be bit-for-bit identical to the serial one, and why that is a feature to manage, not a bug to fear.

Skills applied: OpenACC offload and the reduction clause (§35.2); CUDA Fortran kernels, device variables, one-based indexing (§35.3); host–device transfer and residency (§35.4); the data-race hazard (shared-memory parallelism, Chapter 33); floating-point associativity and reproducibility (Chapter 20, Chapter 37).

Background

We use the $5 \times 5$ plate exactly as the solver left it after two steps in the Chapter 24 / Chapter 35 Project Checkpoint — hot top edge, warmth crept inward — because its numbers are small enough to check every reduction by hand:

  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

We want two scalar diagnostics of this field: the total $\sum u_{ij}$ (a proxy for total thermal energy) and the sum of squares $\sum u_{ij}^2$, whose square root is the $L_2$ norm the capstone's validation study will track. By hand: the top row contributes $5 \times 100 = 500$; row 2 contributes $28 + 32 + 28 = 88$; row 3 contributes $4 + 4 + 4 = 12$; the rest are zero. So the total is $600$. The squares: $5 \times 100^2 = 50{,}000$, plus $28^2 + 32^2 + 28^2 = 784 + 1024 + 784 = 2592$, plus $3 \times 4^2 = 48$, giving sum of squares $= 52{,}640$ and $L_2 = \sqrt{52640} \approx 229.43$. Hold those three numbers; every version must reproduce them.

Here is a plain-Fortran reference (no GPU) that computes them, so the target is concrete and hand-checked:

program plate_energy
  use, intrinsic :: iso_fortran_env, only: dp => real64
  implicit none
  integer, parameter :: nx = 5, ny = 5
  real(dp) :: u(nx,ny), total, sumsq
  integer  :: i, j

  u = 0.0_dp
  u(1,:) = 100.0_dp                          ! hot top edge
  u(2,2) = 28.0_dp;  u(2,3) = 32.0_dp;  u(2,4) = 28.0_dp
  u(3,2) =  4.0_dp;  u(3,3) =  4.0_dp;  u(3,4) =  4.0_dp

  total = 0.0_dp;  sumsq = 0.0_dp
  do j = 1, ny
    do i = 1, nx
      total = total + u(i,j)
      sumsq = sumsq + u(i,j)**2
    end do
  end do

  print '(a, f8.2)',  'total thermal energy (sum) = ', total
  print '(a, f10.2)', 'sum of squares             = ', sumsq
  print '(a, f8.2)',  'L2 norm = sqrt(sumsq)      = ', sqrt(sumsq)
end program plate_energy

$ gfortran -std=f2018 -Wall plate_energy.f90 -o energy && ./energy
total thermal energy (sum) =   600.00
sum of squares             =  52640.00
L2 norm = sqrt(sumsq)      =   229.43

Phase 1 — The Naive Version That Corrupts Its Answer

The obvious first attempt at a GPU reduction is to offload the loop and let every thread add its cell into a shared accumulator:

total = 0.0_dp
!$acc parallel loop collapse(2) copyin(u) copy(total)   ! WRONG: a data race on `total`
do j = 1, ny
  do i = 1, nx
    total = total + u(i,j)          ! every thread reads-modifies-writes the SAME variable
  end do
end do

This is a data race, the shared-memory hazard from Chapter 33, now amplified across thousands of threads instead of a handful. Each thread reads total, adds its cell, and writes it back; when two threads do that concurrently, one overwrites the other's contribution and cells are silently lost. The result is not just wrong — it is nondeterministically wrong, changing from run to run, which is the worst kind of bug in scientific code. The compiler may even refuse it, but if it compiles, it lies. A reduction cannot be an ordinary parallel loop; the combining step needs special handling.

Phase 2 — The OpenACC Way: One Clause

OpenACC's reduction clause is that special handling, and it is a single word. It gives each thread its own private partial sum, then combines the partials safely at the end — exactly as OpenMP's reduction did on the CPU in Chapter 33, now on the device across thousands of threads:

total = 0.0_dp;  sumsq = 0.0_dp
!$acc parallel loop collapse(2) reduction(+:total,sumsq) copyin(u)
do j = 1, ny
  do i = 1, nx
    total = total + u(i,j)
    sumsq = sumsq + u(i,j)**2
  end do
end do
! total -> 600.00,  sumsq -> 52640.00,  L2 = sqrt(52640) ~ 229.43  (matches the reference)

That is the whole build. The GPU performs a tree of partial sums — threads combine within a warp, warps within a block, blocks across the grid — and hands you the final scalars, race-free, in total and sumsq. The result reproduces the reference exactly: $600$ and $52{,}640$. For the overwhelming majority of scientific reductions, this is where you should start and stop.

Phase 3 — The CUDA Fortran Way: What the Clause Hides

To see what reduction is doing for you — and what it costs to do by hand — build the reduction explicitly in CUDA Fortran. The simplest correct approach has every thread atomically add its cell to a device accumulator:

module energy_mod
  use, intrinsic :: iso_fortran_env, only: dp => real64
  use cudafor
  implicit none
contains
  attributes(global) subroutine energy_kernel(u, total, n)
    real(dp) :: u(n)                  ! device array (flattened field)
    real(dp) :: total                 ! device scalar accumulator (set to 0 on the host first)
    integer, value :: n
    integer  :: i
    real(dp) :: old
    i = (blockIdx%x - 1)*blockDim%x + threadIdx%x     ! one-based global index (§35.3)
    if (i <= n) old = atomicadd(total, u(i))          ! atomic: no lost updates
  end subroutine energy_kernel
end module energy_mod

atomicadd(total, u(i)) performs the read-modify-write as one indivisible operation, so no update is ever lost — it is the race of Phase 1, made safe by hardware. Launched over the flattened $25$-element field, it too returns $600$. But two honest caveats show why this is the hard road:

• Double-precision atomicadd needs compute capability ≥ 6.0 (NVIDIA's Pascal generation, 2016 and later). On older devices you would have to emulate it — real, fiddly code. The directive hides that portability concern.
• Every thread atomically updating one location serializes those updates. With 25 cells it is nothing; with ten million, thousands of threads queue for the same memory address and the reduction crawls. An efficient hand-written CUDA reduction avoids that with a shared-memory tree inside each block, then combines one partial per block — a well-known but genuinely intricate kernel, easy to get subtly wrong. The reduction clause generates essentially that tree for you.

So for a reduction, CUDA Fortran's "more control" is mostly more burden: you inherit a hardware caveat and a performance trap that the OpenACC clause quietly handles. This is the clearest case in the chapter where the directive path is not just easier but better for most people.

Phase 4 — An Honesty About Floating Point

There is one more thing the parallel reduction changes, and a scientific programmer must know it. The GPU adds the cells in a different order than the serial loop — it combines partial sums in a tree, not left to right — and floating-point addition is not associative (Chapter 20): $(a + b) + c$ can differ, in the last bits, from $a + (b + c)$. For our plate the values ($100$, $28$, $32$, $4$) are small integers, exactly representable, and the sums ($600$, $52{,}640$) are exact — so every version agrees to the last digit. But for a general real field, expect the GPU reduction to differ from the serial sum in the final bit or two, and to possibly differ from run to run if the thread scheduling varies.

That is not a bug; it is the nature of parallel floating-point summation, and it matters for reproducibility (Chapter 37): a regression test that demands bit-for-bit agreement between the CPU and GPU reductions will fail for a correct program. The right discipline is to compare within a tolerance, not bit-for-bit, and to record that the GPU reduction is non-deterministic at the last-bit level. Knowing this before it bites you in a test suite is the difference between trusting your GPU results and chasing a phantom bug.

Phase 5 — The Verdict

Three builds, one number ($600$), and a clear lesson about the trade at the heart of the chapter:

For reductions, the OpenACC clause is not merely the easy path (§35.2) — it is, for most scientific code, the right one, because it hides a hardware caveat, a serialization trap, and a tree of index arithmetic behind a single, correct word. Reach for explicit CUDA only when profiling proves the directive's reduction is your bottleneck and you have control to gain. And whichever you choose, compare reductions by tolerance, not bits.

Discussion Questions

1. Phase 1's naive version sometimes prints the correct $600$ anyway — especially on a tiny grid or a slow GPU. Why is "it worked on my test" the most dangerous possible outcome for a data race, and how does that argue for the reduction clause even when the naive version "seems fine"?
2. The atomicadd version is correct but can be slow for a large field. Explain the serialization, and sketch in words how a per-block shared-memory tree reduction avoids it. Why does OpenACC not make you write that?
3. A regression test asserts the GPU and CPU total energies are exactly equal and fails on a real (non-integer) field. Is the code wrong? What should the test assert instead, and which two chapters' ideas justify the change?

Your Turn: Extensions

• Option A. Extend the reference program to also compute maxval(u) and its location, then write the OpenACC version using reduction(max:peak). Predict the maximum for the plate (it is on the hot edge) before you run it. Which reductions does OpenACC support besides +?
• Option B. Fold the energy diagnostic into the resident solver of the Project Checkpoint: compute the total energy every 100 steps without copying the whole field back — the reduction runs on the device and only the scalar result comes to the host. What is the transfer cost of a scalar versus the whole field, and why is an on-device reduction the right way to monitor a long run?
• Option C. Construct a small real (non-integer) field — say u(i,j) = 1.0_dp/real(i+j, dp) — and argue, without a GPU, why the serial sum and a tree-ordered sum could differ in the last bit. What tolerance would you set for a CPU-vs-GPU regression test, and how would you justify it?

Key Takeaways

• A reduction (sum, norm, max) cannot be an ordinary parallel loop: thousands of threads writing one accumulator is a data race. It needs special combining.
• OpenACC's reduction(+:var) clause does that combining correctly in one word — private partials, a safe tree combine — and is the right first (and usually last) tool for a GPU reduction.
• The CUDA Fortran equivalent (atomicadd, or a hand-written tree) shows what the clause hides, and inherits real burdens: a double-precision compute-capability floor and a serialization trap the directive avoids.
• A parallel floating-point reduction is not bit-for-bit identical to the serial one (addition is not associative). Compare diagnostics by tolerance, not bits, and treat GPU reductions as last-bit non-deterministic — a reproducibility fact to manage, not a bug to chase.