// HACKER NEWS — CYBERSECURITY
Measuring Gauss-Seidel loop-carried dependency and fixing it via loop unrolling
Last time, we explored the Gauss-Seidel method. We arrived at the following conclusion: the math says Gauss-Seidel takes half as many iterations to converge, and this is backed up by our numerical test case. Yet, despite requiring only half the number of iterations, the Gauss-Seidel solver takes 4 to 5 times as long as the Jacobi one. A seemingly puzzling fact if you think about it: how can a method which is better on paper perform worse when implemented numerically? As we’ve seen, the explanation lies in that, unlike Jacobi, the Gauss-Seidel update rule has loop-carried dependencies which prevent the compiler from vectorizing the code. The question we’ll try to answer today thus is: is this the end of the road or can we somehow recover Gauss-Seidel’s convergence advantage without sacrificing the hardware efficiency of Jacobi?
And it turns out that, yes you can! We have two different routes though: one fairly general where you help the compiler as much as you can (but you need to understand a bit how compilers and CPUs work), and another, quite specific to our particular 2D Poisson test case, where you let the math guide you. None of them will be a smooth ride, but we’ll learn a lot along the way. Eventually, we’ll explore both, but for now we’ll ride along the first one and discover what loop unrolling is.
Let’s rewind for a second. Last time, we pitted two update rules against each other on the same 2D Poisson problem: the humble Jacobi iteration, and its slightly more sophisticated cousin, Gauss-Seidel.
Jacobi always reads from the old iterate v and writes into a fresh array u. Gauss-Seidel, in contrast, updates in place:
That is a one-line difference that buys you a faster-converging method, at exactly half the number of sweeps to reach a given tolerance. On paper, Gauss-Seidel should win outright. And yet, when we actually ran the numbers, Gauss-Seidel took 4 to 5 times longer than Jacobi, wall-clock, despite needing only half as many iterations. Something in that one-line difference is costing us dearly, and it isn’t showing up anywhere in the convergence theory.
The culprit, as we found, is that u(i-1,j) in the Gauss-Seidel update has already been overwritten by the time you read it. It’s not the old value from the previous sweep, it’s this sweep’s freshly computed neighbor. That’s precisely what makes Gauss-Seidel converge faster: you’re propagating information within a single sweep instead of waiting a full iteration for it. But it also means iteration i cannot start until iteration i-1 has fully finished writing u(i-1,j). Jacobi has no such constraint: every read in a Jacobi sweep comes from v, which nothing in that sweep ever touches.
This is the textbook definition of a loop-carried dependency, and it’s exactly the kind of thing a compiler needs to reason about before deciding whether it can vectorize a loop, reorder instructions, or overlap iterations on an out-of-order core. Jacobi: no such dependency, compiler’s hands are free. Gauss-Seidel: hard dependency, compiler’s hands are tied.
That’s a nice story, but it’s still just a story until you can measure it. So this time, rather than taking the compiler’s word for it, we reached for OSACA — the Open Source Architecture Code Analyzer — and pointed it directly at the compiled assembly of both kernels.
OSACA reads a marked assembly loop and, assuming an idealized out-of-order core, reports two numbers per loop body: the critical path (CP), the length of the longest dependency chain running through one iteration in isolation, and the loop-carried dependency (LCD), the part of that chain that must cross from one iteration to the next. LCD is the number that matters here. It’s a hard floor on cycles-per-iteration, no amount of clever scheduling can get under it.
Here’s what came back for our -O3 -march=native kernels, compiled with plain do loops (no do concurrent, nothing fancy):