Case Study 1: The Energy Diagnostic That Depended on Loop Order

"Two runs, same data, same machine, different answers. One of them has to be lying — and it turns out both are telling the truth about floating-point arithmetic."

Executive Summary

A research group's 2D heat solver prints a "total thermal energy" diagnostic every few steps by summing the temperature field. After a refactor that changed the field from row-major to column-major traversal (for the cache-friendliness reasons of a later performance chapter), the diagnostic began reporting a different total for the same field — and worse, a conservation check that should have held to machine precision started failing by a wide margin. No cell value changed; only the order of a summation did.

This case study traces that bug to its root — absorption and the non-associativity of floating-point addition (§20.3) — reproduces it in a few lines you can compile, shows that the identical behaviour appears in Python (because it is a property of IEEE 754, not of the language), quantifies when it matters for a real field, and lays out the fix that Case Study 2 then builds.

Skills applied - Reading the anatomy of a real(dp) and its ULP to predict absorption (§20.1, §20.2). - Diagnosing catastrophic cancellation and non-associativity in a summation (§20.3). - Using spacing to reason quantitatively about which contributions are lost (§20.2). - Recognizing that floating-point behaviour is a hardware-standard property shared with Python (§20.1).

Background

The heat solver stores temperature as a 2D real(dp), allocatable :: u(:,:) array (Chapter 5). A common, cheap diagnostic is the sum of all cell temperatures, a proxy for total thermal energy that should change only through the boundaries. The original code summed the field in one order; the refactored code summed it in another. Mathematically, addition is associative and commutative, so the order cannot matter. In floating point, it can — and the plate under study had a few very hot cells (temperatures near $10^{8}$ in the code's non-dimensional units, from a concentrated source term) among many cells near $1$.

The essential ingredient is the one from §20.2: the gap between representable doubles — one ULP — grows with magnitude. Near $10^{17}$ it exceeds $1$; even near $10^{8}$ it is far larger than the rounding you would expect. When you add a small number to a large running sum, the small number can fall into the gap and vanish entirely. That is absorption, and it makes the sum depend on the order in which you add.

Phase 1 — Reproduce the Bug

Strip the diagnostic to its essence: a left-to-right sum of four numbers, whose multiset is fixed but whose order we vary. Two of the numbers are large and cancel exactly; two are small.

module kinds
  implicit none
  private
  public :: dp
  integer, parameter :: dp = selected_real_kind(15, 307)
end module kinds

module summation
  use kinds, only: dp
  implicit none
contains
  pure function sum_lr(v) result(s)   ! sum left-to-right, exactly as a do-loop does
    real(dp), intent(in) :: v(:)
    real(dp)             :: s
    integer              :: i
    s = 0.0_dp
    do i = 1, size(v)
       s = s + v(i)
    end do
  end function sum_lr
end module summation

program energy_order
  use kinds, only: dp
  use summation, only: sum_lr
  implicit none
  real(dp) :: order_a(4) = [ 1.0e17_dp, -1.0e17_dp,  1.0_dp,  1.0_dp ]
  real(dp) :: order_b(4) = [ 1.0e17_dp,  1.0_dp,  1.0_dp, -1.0e17_dp ]

  print '(a, f6.1)', 'sum, order A = ', sum_lr(order_a)
  print '(a, f6.1)', 'sum, order B = ', sum_lr(order_b)
end program energy_order

$ gfortran -std=f2018 -Wall case1.f90 -o case1 && ./case1
sum, order A =    2.0
sum, order B =    0.0

The same four numbers, summed in two orders, give 2.0 and 0.0. Neither is a bug in the code; both are the correct IEEE result for the order given.

Phase 2 — Diagnose It

Trace each order by hand, using the fact from §20.2 that at $10^{17}$ one ULP is $2^{57-53} = 2^{4} = 16$, so any value below half a ULP — below $8$ — added to $10^{17}$ is absorbed and disappears.

Order A [1e17, -1e17, 1.0, 1.0]:

$$ (10^{17} + (-10^{17})) + 1 + 1 = 0 + 1 + 1 = 2 $$

The two large values cancel first, exactly (they are equal in magnitude), and the small values then accumulate cleanly on a running sum of $0$. Result: 2.0, the true total.

Order B [1e17, 1.0, 1.0, -1e17]:

$$ ((10^{17} + 1) + 1) - 10^{17} = (10^{17} + 1) - 10^{17} = 10^{17} - 10^{17} = 0 $$

Here each 1.0 is added to $10^{17}$, falls into the ULP gap ($1 < 8$), and is absorbed — the running sum stays exactly $10^{17}$. When the large value is finally subtracted, nothing of the small values remains. Result: 0.0, off by the entire true total.

The core fact. Floating-point addition is commutative but not associative: $(a + b) + c$ need not equal $a + (b + c)$. Order A groups the cancellation early; Order B lets the large value swallow the small ones first. This is the same non-associativity you met in §20.3, now wearing the costume of a physics diagnostic.

Phase 3 — Confirm It Is the Hardware, Not Fortran

A tempting hypothesis is that Fortran's summation is somehow flawed. It is not — the behaviour is a property of IEEE 754, and the identical thing happens in Python, whose float is the same binary64:

>>> sum([1e17, -1e17, 1.0, 1.0])     # order A
2.0
>>> sum([1e17, 1.0, 1.0, -1e17])     # order B
0.0

🐍 Python Comparison. Python's built-in sum also adds left-to-right, so it reproduces Order A and Order B exactly. NumPy's np.sum, interestingly, often does pairwise summation for large arrays, so it can give yet a third answer — closer to the truth, but still order-sensitive at the bit level. The lesson for a mixed Fortran/Python workflow: do not expect a Fortran sum, a Python sum, and a NumPy np.sum to be bit-identical, because they may add in different orders. Reproducibility across the two requires agreeing on the summation algorithm, not just the data.

Phase 4 — Quantify: When Does It Bite the Solver?

For the real field, the question is: how many small contributions can a running sum silently lose? The answer follows directly from spacing. Once the running sum reaches magnitude $S$, every added value below $\tfrac{1}{2}\,\text{spacing}(S)$ is absorbed. The table shows the threshold:

For a plate with a few cells near $10^{8}$, the running sum reaches $\sim 10^{8}$ quickly, and thereafter every cool cell (temperature $\sim 1$) contributes only if it survives a $10^{-8}$ ULP — which it does, barely, so the error is a slow drift rather than a total wipeout. But over a large grid the drift accumulates, and it changes with traversal order, which is why the row-major-to-column-major refactor made a previously-hidden error suddenly visible. The diagnostic was never trustworthy; the refactor just stopped it from being wrong in a consistent way.

⚡ Performance Note. The fix must not cost much: the energy diagnostic runs every few steps over the whole field, so an $O(n)$ sum is fine but an $O(n \log n)$ sort every time is not. This rules out "sort ascending then add," elegant as it is, for the hot path — and points toward compensated summation (Case Study 2), which is $O(n)$ with a tiny constant factor.

Phase 5 — The Fix, and a Sanity Check

Three cures exist, in increasing order of robustness:

1. Add small-to-large. Sorting the values ascending before summing keeps the running sum small for as long as possible, minimizing absorption. Correct, but $O(n \log n)$ — too slow for a per-step diagnostic.
2. Pairwise (tree) summation. Recursively sum halves and add the two partial sums. $O(n)$ time, dramatically smaller error growth than the naive sweep, and what NumPy does. A good default.
3. Compensated (Kahan) summation. Carry a running correction term that captures the low-order bits lost at each addition. $O(n)$ with a small constant, and accurate almost to the last bit regardless of order. This is what Case Study 2 builds and folds into the solver.

As an immediate sanity check, the group re-ran the conservation test using the intrinsic sum (which the compiler may implement more carefully than a naive loop) on a scaled field — dividing every temperature by the maximum before summing, then multiplying back — so that the running sum never leaves the well-conditioned neighbourhood of $1$. The conservation error dropped back to a few ULPs, confirming the diagnosis: the problem was never the physics, only the arithmetic.

Discussion Questions

1. Order A gave the exact answer 2.0 by cancelling the large values first. Is "cancel large values first" a reliable general strategy, or did this example just get lucky? What could go wrong if the two large values were not exactly equal?
2. The intrinsic sum is permitted by the standard to add in any order (and to use extra-precision accumulators). Is that freedom a feature or a hazard for reproducibility? When would you want to forbid it?
3. spacing(1.0e16_dp) is 2.0. Explain to a colleague, without jargon, why that means "you cannot store every whole number near ten quadrillion."
4. The refactor that exposed the bug was made for cache performance (Chapter 27). Argue that the refactor was still the right thing to do, even though it "broke" the diagnostic.

Your Turn: Extensions

• Option A (analyze). Modify sum_lr to also return the number of additions in which the running sum did not change (i.e., the added value was fully absorbed). Run it on a field of one $10^{17}$ value followed by one thousand $1.0$s and report how many were lost.
• Option B (measure). Sum the array [0.1_dp, 0.1_dp, ..., 0.1_dp] (ten million copies) left-to-right and compare with 1.0e6_dp (the exact mathematical total). How large is the drift, and does it match the $\sim N u$ order-of-magnitude estimate from §20.2?
• Option C (port). Write the same ten-million-0.1 sum in Python (a plain loop, then math.fsum, then numpy.sum) and compare all three with the Fortran result. Explain the differences using this chapter, not language folklore.

Key Takeaways

• Floating-point addition is not associative: the order of a summation changes the result, sometimes catastrophically, because of absorption — small values falling into the ULP gap of a large running sum.
• The effect is a property of IEEE 754, shared identically by Fortran, Python, and NumPy (though their default summation orders differ, so their results can differ).
• Use spacing(S) to predict the absorption threshold: values below $\tfrac{1}{2}\,\text{spacing}(S)$ are lost once the running sum reaches $S$.
• A diagnostic whose value depends on loop order was never trustworthy. The cure is a better summation algorithm — pairwise or compensated — not a rearrangement of the same naive sweep. Case Study 2 builds it.