Make Gauss-Seidel great again!
The art of 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.
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 doJacobi 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 doThat 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):
| 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. The arithmetic hasn’t changed; only what’s allowed to overlap has.
Second, and this is the whole ballgame: Jacobi’s loop-carried dependency is a throwaway 1 cycle (just the loop counter incrementing, nothing to do with the actual computation). Gauss-Seidel’s is 12 cycles, and OSACA even tells you exactly where it comes from:
46 | 12.0 | vaddsd %xmm1, %xmm0, %xmm0 | [46, 47, 48]Trace it back and it’s the west-neighbor term: %xmm1 holds this sweep’s freshly written u(i,j), and it gets read straight back in as u(i-1,j) on the very next loop trip. Three chained instructions — the west-neighbor add, the north-neighbor add, and the final multiply — each with 4-cycle latency, sum to exactly 12. That’s not a compiler being lazy; it’s the compiler correctly reporting that there is no way to start computing u(i,j) before u(i-1,j) exists, because the algorithm itself demands it.
The other giveaway, readable straight off the assembly without even needing OSACA, is the instruction mnemonics themselves. Jacobi’s inner loop is built entirely out of packed double-precision instructions — vaddpd, vmulpd, vfmadd213pd — each one crunching four grid points at a time in a single 256-bit AVX register. Gauss-Seidel’s inner loop uses the scalar forms of the exact same instructions — vaddsd, vmulsd, vfmadd213sd — one grid point at a time. The pd/sd suffix alone tells you, at a glance, which loop the compiler was able to vectorize and which one it wasn’t.
Does the number add up?
For each kernel, the CPU’s actual sustained cost per iteration is bounded below by both constraints at once, the throughput floor and the LCD, so the one that actually governs is whichever of the two is larger. For Gauss-Seidel, that’s the LCD (12 cycles, versus a 2.5-cycle throughput floor). For Jacobi, it’s the throughput floor (3 cycles, versus a negligible 1-cycle LCD). In other words, Gauss-Seidel is latency-bound and Jacobi is throughput-bound. Two different bottlenecks entirely, which is precisely why comparing them isn’t just a matter of reading off the same column for both.
Normalizing each kernel’s governing bound per grid point (dividing Jacobi’s by 4, since it processes four elements per vector iteration) gives a predicted per-sweep slowdown of roughly
\[\dfrac{12}{3 / 4} = 16\times\]
That’s higher than the 4–5× we measured wall-clock. But remember, wall-clock measures total time to solution, and Gauss-Seidel needs only half as many sweeps as Jacobi to converge. Once you back that out, the actual per-sweep slowdown implied by the measurements is
\[\left(\frac{T_{GS}}{T_J}\right) \Big/ \left(\frac{N_{GS}}{N_J}\right) \approx (4\text{–}5) \times 2 \approx 8\text{–}10\times\]
which lands within a factor of two of OSACA’s prediction, a healthy level of agreement for a static model that has zero visibility into the memory system. The mechanism OSACA identifies — Gauss-Seidel latency-bound at a hard 12-cycle wall, Jacobi comfortably throughput-bound — is exactly right, even if the precise number needs a pinch of salt.
So, the diagnosis is confirmed in triplicate: theory, static analysis, and measurement all agree. Gauss-Seidel’s faster convergence is bought entirely at the price of a serialized, 12-cycle-per-element dependency chain that the compiler simply cannot break on its own. Which raises the obvious question: is that dependency chain fundamental to the algorithm, or is it an artifact of how we happened to write the loop?
So what is loop unrolling then?
Loop unrolling is one of the oldest tricks in the optimizing-compiler playbook, and the idea is disarmingly simple. Instead of executing a loop body once per iteration and paying the overhead of a branch and a counter update every single time, you write out the body two, four, or eight times per iteration, advancing the counter by that many steps at once.
Take a trivial example, nothing to do with our Poisson problem yet:
Unrolled by a factor of two, this becomes:
modulo a bit of bookkeeping if n is odd, which we’ll conveniently ignore for now. Fewer branches, fewer counter increments, and — crucially — twice as much independent work exposed to the compiler within a single loop trip. That last point is really the whole appeal: with two statements sitting side by side and no dependency between them, the compiler is free to interleave them, pack them into a single vector instruction, or schedule them however the hardware likes. Unrolling doesn’t make the CPU do less work; it makes the scheduler’s job easier by giving it more to look at at once. So, naturally, the question we should ask is: does the same trick save Gauss-Seidel?
A tempting shortcut that doesn’t work
Let’s try the obvious thing and unroll our Gauss-Seidel inner loop by two, changing nothing else:
do i = 2, n-2, 2
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))
u(i+1, j) = 0.25_dp*(dx2*b(i+1,j) + u(i,j) + u(i+2,j) &
+ u(i+1,j-1) + u(i+1,j+1))
end doLook closely at the second statement: it reads u(i,j) which is exactly the value the first statement just finished computing, two lines above. We haven’t removed the loop-carried dependency at all. We’ve just moved a copy of it inside the loop body, where it’s now a same-iteration dependency instead of a cross-iteration one. The compiler still can’t touch these two statements independently, still can’t vectorize them together, and still has to execute them strictly in order. We’ve paid for the bookkeeping of unrolling (larger code, an odd-n remainder to handle) and bought ourselves precisely nothing.
This is worth sitting with for a second, because it’s a genuinely common trap: unrolling looks like it should help with dependencies, since it’s so often reached for as a vectorization enabler, but by itself it does nothing whatsoever to break a true dependency chain. It’s very good at exposing parallelism that already exists (as in our trivial a(i) = 2*a(i) example) but is completely powerless against parallelism that has to be created. Gauss-Seidel falls squarely in the second category: the dependency isn’t a scheduling accident the compiler failed to notice, it’s baked into the mathematics of the update rule.
So if unrolling alone is a dead end, what’s missing? The answer is that the dependency itself has to be rewritten away before unrolling can do any good. That means reaching for a bit of algebra, not just a bit of loop restructuring.
Fair enough, let’s do that for Gauss-Seidel!
To keep the algebra honest and readable, let’s forget about the full 2D grid for now, and work with the 1D Poisson equation, \(u'' = -f\), discretized on a uniform grid as
\[\frac{u_{i+1} - 2 u_i + u_{i-1}}{\Delta x^2} = -f_i.\]
Solving for \(u_i\) gives the update rule
\[u_i = \dfrac{1}{2}\left(u_{i-1} + u_{i+1} + \Delta x^2 f_i\right),\]
so, unsurprisingly, the constant playing the role of our earlier \(0.25\) is now \(c = \tfrac12\). Sweeping left to right and updating in place, Gauss-Seidel reads:
Exactly one neighbor here is a problem: u(i-1) was overwritten by the previous loop trip, while u(i+1) is still the untouched, old value (the sweep hasn’t reached it yet). That’s the whole dependency, laid completely bare with nothing else to distract from it.
Writing it as a recursion
Group everything that doesn’t depend on the sweep’s own progress into a single per-column quantity,
\[t_i = u_{i+1} + \Delta x^2 f_i,\]
which is fully known before the block starts (it’s old data, untouched by anything we’re computing right now). The update collapses to a clean, one-term affine recursion:
\[u_i = c\left(u_{i-1} + t_i\right), \qquad u_{i-1} \equiv \texttt{um},\]
where um is shorthand for the last already-updated value feeding into this block. This is precisely the shape that makes the dependency chain visible as just a recursion in \(u\) — nothing more, nothing less.
Solving the recursion two steps ahead
Rather than computing \(u_i\) and then \(u_{i+1}\) from it, we solve the recursion explicitly for both, purely in terms of quantities that existed before the block started:
\[u_i = c\,(\texttt{um} + t_i)\]
\[u_{i+1} = c\,(u_i + t_{i+1}) = c\Big(c\,(\texttt{um}+t_i) + t_{i+1}\Big) = c^2\,\texttt{um} + c^2\,t_i + c\,t_{i+1}\]
Look at what happened: \(u_{i+1}\) no longer references \(u_i\) at all. Both outputs are now pure functions of um, \(t_i\), and \(t_{i+1}\). All quantities that existed before this pair of updates began. The intra-pair dependency that plain unrolling couldn’t touch (the “unroll-only” kernel from the previous section) has been algebraically dissolved.
In code, defining \(c_2 \equiv c^2\), this gives the following kernel:
do i = 2, n-2, 2
um = u(i-1)
t0 = u(i+1) + dx2*b(i)
t1 = u(i+2) + dx2*b(i+1)
u(i) = c *(um + t0)
u(i+1) = c2*um + c2*t0 + c*t1 ! Does not depend on u(i) !
end doTwo independent expressions, no read of one output by the other’s computation. Exactly what the compiler needs to schedule them concurrently instead of serially.
Back to 2D: nothing new needed
This is where the 1D detour pays off. In the real 2D stencil, the north, south, and source contributions play no role whatsoever in creating the i-dependency — they’re either values from a row Gauss-Seidel has already finished, or values from a row it hasn’t touched yet, but never anything the i-loop itself is computing. So they simply fold into the same lumped \(t_i\) we just defined:
\[t_i \;\longrightarrow\; \texttt{tmp}_i = \Delta x^2 b(i,j) + u(i+1,j) + u(i,j+1) + u(i,j-1),\]
and the exact same two-line substitution applies, verbatim, giving you the unrolled kernel:
real(dp), parameter :: c = 0.25_dp
real(dp), parameter :: c2 = 0.0625_dp
do j = 2, n-1
do i = 2, n-2, 2
um = u(i-1, j)
tmp1 = dx2*b(i, j) + (u(i+1, j) + u(i, j+1) + u(i, j-1))
tmp2 = dx2*b(i+1, j) + (u(i+2, j) + u(i+1, j+1) + u(i+1, j-1))
u(i, j) = c *(um + tmp1)
u(i+1, j) = c2*um + (c2*tmp1 + c*tmp2)
end do
! Handles the case where an odd number of grid points is used.
if (mod(n, 2) == 1) then
tmp1 = dx2*b(n-1, j) + (u(n, j) + u(n-1, j+1) + u(n-1, j-1))
u(n-1, j) = c*(u(n-2, j) + tmp1)
end if
end doThe trailing if block just mops up a leftover column when n is odd, i.e. the same remainder bookkeeping any unrolled loop needs, nothing conceptually new.
Does it actually work?
Ok, so the new kernel looks a bit more complicated than our original one. But mathematically, it’s all the same thing, so is it truly more efficient? Time to check the theory against the machine. Same workflow as before: compile with -O3 -march=native -S, mark the inner loop, and let OSACA loose on it.
Port pressure in cycles
| 0 - 0DV | 1 | 2 - 2D | 3 - 3D | 4 | 5 | 6 | 7 || CP | LCD |
--------------------------------------------------------------------------------------------------
74 | | | | | | | | || | | .L5:
75 | | | 0.50 0.50 | 0.50 0.50 | | | | || | | vmovsd 16(%rax), %xmm1
76 | | | 0.50 0.50 | 0.50 0.50 | | | | || | | vmovsd 8(%rax), %xmm0
77 | 0.00 | 0.00 | | | | 0.50 | 0.50 | || | | addq $16, %rdx
78 | 0.50 | 0.50 | 0.50 0.50 | 0.50 0.50 | | | | || 4.0 | | vaddsd 8(%rax,%rsi,8), %xmm1, %xmm1
79 | 0.50 | 0.50 | 0.50 0.50 | 0.50 0.50 | | | | || | | vaddsd (%rax,%rsi,8), %xmm0, %xmm2
80 | 0.50 | 0.50 | 0.50 0.50 | 0.50 0.50 | | | | || 4.0 | | vaddsd 8(%rax,%rcx,8), %xmm1, %xmm1
81 | 0.50 | 0.50 | 0.50 0.50 | 0.50 0.50 | | | | || | | vaddsd (%rax,%rcx,8), %xmm2, %xmm2
82 | 0.00 | 0.00 | | | | 0.50 | 0.50 | || | | addq $16, %rax
83 | 0.50 | 0.50 | 0.50 0.50 | 0.50 0.50 | | | | || 4.0 | | vfmadd231sd -8(%rdx), %xmm5, %xmm1
84 | 0.50 | 0.50 | 0.50 0.50 | 0.50 0.50 | | | | || | | vfmadd231sd -16(%rdx), %xmm5, %xmm2
85 | 0.50 | 0.50 | | | | | | || 4.0 | | vmulsd %xmm4, %xmm1, %xmm1
86 | 0.50 | 0.50 | | | | | | || | | vaddsd %xmm3, %xmm2, %xmm0
87 | 0.50 | 0.50 | | | | | | || | | vmulsd %xmm4, %xmm0, %xmm0
88 | 0.50 | 0.50 | | | | | | || 4.0 | | vfmadd132sd %xmm6, %xmm1, %xmm2
89 | 0.50 | 0.50 | | | | | | || 4.0 | 4.0 | vfmadd132sd %xmm6, %xmm2, %xmm3
90 | | | | | | 1.00 | | || 1.0 | | vunpcklpd %xmm3, %xmm0, %xmm0
91 | | | 0.00 | 0.00 | 1.00 | | | 1.00 || 0.0 | | vmovupd %xmm0, -16(%rax)
92 | 0.00 | 0.00 | | | | 0.00 | 1.00 | || | | cmpq %rax, %rdi
93 | | | | | | | | || | | * jne .L5
5.50 5.50 4.00 4.00 4.00 4.00 1.00 2.00 2.00 1.00 29.0 4.0If you look at the last two columns, two things should stand out. The critical path CP got a bit longer, but the loop-carried dependencies (LCD) dropped from 12 to 4 cycles per loop trip — and since each trip now covers two elements, that’s roughly a 6× drop per grid point.
| LCD (per element) | Throughput floor (per element) | Governing bound | |
|---|---|---|---|
| Gauss-Seidel (scalar, degree 1) | 12.0 | 2.5 | latency: 12.0 |
| Kernel 2 (degree 2) | 2.0 | 2.75 | throughput: 2.75 |
That’s the static prediction. Time to see if the actual machine agrees — and, more interestingly, to see why it agrees, one compilation flag at a time.
Experiment 1 — baseline, no optimization
First, a sanity check. Compiled with -O0, so no scheduling cleverness of any kind gets to play a role:
| Solver | Iterations | Time / iteration | Total |
|---|---|---|---|
| Jacobi | 138 000 | 2 ms | 273 s |
| Textbook Gauss-Seidel | 74 000 | 2.4 ms | 178 s |
| Unrolled Gauss-Seidel | 74 000 | 2 ms | 147 s |
Two things worth checking here before moving on. First, all three solvers reach the solution in the same number of iterations as the earlier posts — confirming the unrolled kernel really is computing the same thing as the textbook one, algebra and all, not some subtly different (and subtly wrong) approximation. Second, at -O0, the unrolled kernel is already reliably faster than the textbook kernel, by almost exactly the ratio convergence theory predicts for Gauss-Seidel over Jacobi (2×) — a nice confirmation that even without any compiler scheduling to exploit, exposing two independent computations per loop trip buys you something, purely from halving the loop overhead. But nobody ships -O0 code, so this is a baseline, not the result.
Experiment 2 — the critical one
Now compile with real optimization, but deliberately keep vectorization switched off (-O3 -march=native -mtune=native -fno-tree-vectorize), so that any gains we see can’t be attributed to SIMD at all:
| Solver | Iterations | Time / iteration | Total |
|---|---|---|---|
| Jacobi | 138 000 | 166 µs | 23 s |
| Textbook Gauss-Seidel | 74 000 | 1054 µs | 78 s |
| Unrolled Gauss-Seidel | 74 000 | 175 µs | 13 s |
This is the number that matters. With vectorization explicitly disabled, the unrolled kernel is still six times faster per sweep than the textbook kernel, and lands within 5% of Jacobi’s per-sweep cost. Since neither Gauss-Seidel kernel can be vectorized under this flag anyway, this gain has nothing to do with packing multiple grid points into one instruction. It comes entirely from something more subtle: reducing the loop-carried dependency lets the CPU’s out-of-order scheduler overlap the scalar work of consecutive iterations, filling execution ports that the textbook kernel’s 12-cycle dependency chain left idle. Fewer cycles spent waiting, same number of instructions issued. This is the empirical confirmation of the LCD-versus-throughput story from the OSACA table above, isolated from vectorization entirely.
Experiment 3 — turning vectorization back on
Finally, the same three kernels, same flags, but with vectorization allowed:
| Solver | Iterations | Time / iteration | Total |
|---|---|---|---|
| Jacobi | 138 000 | 115 µs | 16 s |
| Textbook Gauss-Seidel | 74 000 | 1054 µs | 78 s |
| Unrolled Gauss-Seidel | 74 000 | 175 µs | 13 s |
Both Gauss-Seidel kernels are completely unchanged, down to the microsecond. That’s expected once you know the assembly: gfortran never vectorizes either GS kernel regardless of this flag, so there’s nothing for the flag to act on. Vectorization was never the mechanism behind experiment 2’s speedup. It’s a separate story that only applies to Jacobi.
And even for Jacobi, notice the gain is real but modest: 166 µs → 115 µs, about 1.4×, not the ~4× you might expect from packing four doubles into a single AVX register. This is the by-product the LCD story predicted but didn’t fully explain on its own: vectorization only pays off to the extent that a kernel’s compute is the bottleneck, and Jacobi’s inner loop, despite having no loop-carried dependency at all, is largely memory-bound at this problem size. It streams through the whole grid every sweep with essentially no data reuse, so DRAM bandwidth caps the achievable speedup well below what the instruction-level analysis alone would suggest.
Put the three experiments side by side and the conclusion is unambiguous: it’s the absence (or reduction) of a loop-carried dependency, not vectorization, that lets the scheduler do its job. Vectorization is a nice bonus when it’s available. It does still help Jacobi, just less than you’d naively guess. But it was never the reason Gauss-Seidel was slow, and reducing the LCD was already enough, on its own, to close almost the entire performance gap with Jacobi. Actually, this unrolled Gauss-Seidel kernel hands us the solution almost 3 seconds faster than the Jacobi one. Mission accomplished! Or so it appears…
But look closely at the mnemonics
Here’s the detail that’s easy to miss if you only look at the CP/LCD numbers and declare victory. If you look at the far right of the OSACA output, that is the assembly for inner-most loop (i) of our kernel. Every single arithmetic instruction in that table is still a scalar op — vaddsd, vmulsd, vfmadd231sd, vfmadd132sd, all operating on individual doubles in %xmm registers, not a single packed pd instruction among them. Compare that to Jacobi’s inner loop, built entirely out of vaddpd/vmulpd/vfmadd213pd on 256-bit %ymm registers, four grid points at a time. The compiler never vectorized this loop. It didn’t need to, to hit Jacobi’s throughput. It got there by exposing two independent scalar chains (one for u(i,j), one for u(i+1,j)) that fill different execution ports at the same time, rather than by packing four elements into one instruction the way Jacobi does. The one and only place the two lanes actually meet is right at the end, at line 90: a vunpcklpd stitching the two independently-computed scalars into a single 128-bit register just before the store. That’s not vectorized computation. It’s a vectorized write of two scalar results, a small, free bonus riding on top of the real work, not the source of the speedup itself.
This distinction matters more than it might look like at first glance, and it’s the seed of the next post’s cliffhanger: because the compiler is achieving Jacobi-level throughput through instruction-level parallelism across ports rather than true SIMD, this kernel still can’t be handed a wider vector width for free, and it certainly can’t be split across threads. That tiny 4-cycle pair-to-pair dependency, however cheap, is still a sequential handoff that no do concurrent or OpenMP directive could safely break apart. We’ve matched (and even slightly exceeded) Jacobi’s single-core speed, but not its parallelizability. Those turn out to be two different prizes entirely. And only one of them is in hand.
Can we go even further?
Could we push the substitution one step further? Unroll by 4 instead of 2, and squeeze out even more performance? Not really. We’d be at risk of marginally reducing an LCD that’s already a non-issue, while greatly increasing the pressure on the CPU ports.
Here’s why. Going back to the general recursion from before,
\[u_i = c^{i+1}\,u_{-1} + \sum_{m=0}^{i} c^{i-m}\,t_m,\]
each additional element you fold into the substitution costs you one more multiply-add. By the time you reach the degree-4 unrolling, you’re resolving a 4-term cascading sum, every term weighted by its own power of \(c\). Averaged across the four outputs, that’s roughly twice the multiply-add work per element compared to degree 2. And recall where we ended up last section: the throughput floor is already the governing bound at degree 2 (2.75 cycles/element, versus an LCD of only 2.0). Doubling the arithmetic per element doesn’t touch the LCD’s role in the story at all. It just raises the throughput floor further, from the bound that’s actually constraining us. You’d be spending real port-pressure budget to shave cycles off a dependency chain that was already comfortably hidden behind other work. A losing trade, and one you can see coming from the algebra alone, without even needing to compile anything.
So degree 2 isn’t an arbitrary stopping point. It’s exactly the amount of unrolling needed to push the LCD below the throughput floor, and not a bit more. Past that point, every extra unit of unrolling is pure cost with no corresponding benefit. At least for Gauss-Seidel.
So, where does that leave us?
We made it: a Gauss-Seidel kernel that keeps the mathematical convergence advantage over Jacobi. Same iteration count as the textbook version, while very nearly matching Jacobi’s raw computational efficiency per sweep. A serial win on both fronts at once, which is exactly the “have your cake and eat it too” outcome we set out chasing at the start of this post.
But, there’s always a but, it’s a serial win only. The kernel is still built entirely out of scalar instructions. The compiler (gfortran 15.3) never vectorized it for my particular CPU. It got to Jacobi-level throughput through a different route (independent scalar chains sharing execution ports, not packed SIMD lanes). Worse, that tiny residual pair-to-pair dependency, however cheap, is still a genuine sequential handoff. No do concurrent, no OpenMP directive, no compiler flag can safely split this loop across threads or hand it a wider vector register, because the math still insists that pair \(k\) cannot start until pair \(k-1\) has written its result.
For a single core on a laptop, that’s a complete success. But it’s a dead end for anything that needs to scale. And scaling is, after all, the entire point of numerical computing once your grids get big enough to matter. If we want a Gauss-Seidel kernel that vectorizes and multithreads, unrolling and clever algebra won’t get us there. We’ll need to change the order in which we visit grid points altogether, which is exactly where red-black ordering comes in, next time.
Acknowledgement – I have to give credit where credit is due. It is Ivan Pribec who showed me this particular unrolling trick for Gauss-Seidel. It is also he who put osaca under my radar, so none of what I presented here would have been possible without him. Thanks a lot!
If you want to read more of my stuff