Make Gauss-Seidel great again! – Dr. Jean-Christophe Loiseau
Pangram verdict · v3.3
We believe that this text is a mix of AI and human-written content.
AI likelihood · overall
MixedArticle text · 855 words · 5 segments analyzed
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. Wait, what was the problem again? 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. do j = 2, n-1 do i = 2, n-1 u(i, j) = 0.25_dp*(dx2*b(i,j) + v(i+1,j) + v(i-1,j) & + v(i,j+1) + v(i,j-1)) end do end do Jacobi always reads from the old iterate v and writes into a fresh array u. Gauss-Seidel, in contrast, updates in place: do j = 2, n-1 do i = 2, n-1 u(i, j) = 0.25_dp*(dx2*b(i,j) + u(i+1,j) + u(i-1,j) & + u(i,j+1) + u(i,j-1)) end do end do 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. Loop-carried dependencies, made concrete 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. Putting a number on “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): Summary of the osaca analysis of the Jacobi and standard Gauss-Seidel kernels.
Both kernels have been compiled using gfortran 15.3 with options -O3 -mtune=native -march=native. The precise numbers reported may depend on the exact CPU and compiler you’re using.
CP (cycles) LCD (cycles) Port-pressure floor (cycles) Governing bound Jacobi (vectorized, 4 elements/it.) 24 1 3 throughput Gauss-Seidel (scalar, 1 element/it.) 24 12 2.5 latency Two things jump out immediately. First, the critical path is identical for both kernels — 24 cycles.
That makes sense: it’s essentially the same chain of floating-point operations in both cases: Three adds : tmp1 = v(i+1, j) + v(i-1, j) + v(i, j+1) + v(i, j-1) One multiply : tmp2 = b(i, j) * dx2 One FMA : tmp3 = 0.25_dp * tmp1 + 0.25_dp*tmp2 along with a a load and a store.