Measuring Gauss-Seidel Loop-Carried Dependency and Fixing it via Loop Unrolling

Author: Jean-Christophe Loiseau
Published: September 14, 2026

In our previous discussion, we analyzed the Gauss-Seidel method. The mathematical conclusion was clear: Gauss-Seidel requires only half the iterations to reach convergence compared to the Jacobi method—a fact confirmed by our numerical tests.

However, a paradox emerged: despite the reduction in iterations, the Gauss-Seidel solver is $\approx 4$ to $5$ times slower in terms of wall-clock time than the Jacobi solver.

The Puzzle: How can a method that is mathematically superior perform so poorly when implemented on actual hardware?

The answer lies in the update rule. Unlike Jacobi, Gauss-Seidel introduces loop-carried dependencies, which effectively "handcuff" the compiler and prevent it from utilizing vectorization.

The Path Forward

Is this a fundamental limitation, or can we reclaim the convergence benefits of Gauss-Seidel without losing the hardware efficiency of Jacobi? We have two potential strategies:

• The General Route: Assist the compiler by understanding the interplay between the CPU and the generated machine code.
• The Specific Route: Leverage the mathematical properties of our 2D Poisson test case to rewrite the approach.

Neither path is easy, but both are educational. For now, we will pursue the first route and dive into the concept of loop unrolling.

---

Comparing the Update Rules

To understand the bottleneck, let's revisit the two methods applied to the 2D Poisson problem.

1. The Jacobi Iteration

Jacobi is straightforward: it reads from a "frozen" previous state v and writes the results into a new array u.

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

2. The Gauss-Seidel Iteration

Gauss-Seidel updates the values in place. This one-line change is the secret to its faster convergence.

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

Mathematically, the update for a cell $(i, j)$ can be represented as: $u_{i,j}^{(k+1)} = \frac{1}{4} \left( \Delta x^2 b_{i,j} + u_{i+1,j}^{(k)} + u_{i-1,j}^{(k+1)} + u_{i,j+1}^{(k)} + u_{i,j-1}^{(k+1)} \right)$

While this allows information to propagate across the grid within a single sweep (leading to $\text{half the iterations}$), it creates a strict sequential requirement.

---

Loop-Carried Dependencies: The Culprit

The performance hit comes from the term u(i-1, j). In Gauss-Seidel, by the time the loop reaches index $i$, the value at $i-1$ has already been updated in the current sweep.

The Dependency Chain: Iteration $i$ $\rightarrow$ must wait for $\rightarrow$ Iteration $i-1$ $\rightarrow$ must wait for $\rightarrow$ Iteration $i-2$.

In contrast, Jacobi reads from v, which remains untouched during the sweep. There is no dependency between iterations, allowing the CPU to process multiple elements simultaneously.

---

Quantifying the Bottleneck with OSACA

To move beyond theory, we used OSACA (Open Source Architecture Code Analyzer) to examine the compiled assembly. OSACA analyzes loops on an idealized out-of-order core and provides two key metrics:

1. Critical Path (CP): The longest dependency chain within a single iteration.
2. Loop-Carried Dependency (LCD): The portion of that chain that must cross from one iteration to the next.

The LCD represents a "hard floor" on the cycles-per-iteration; no amount of scheduling can bypass it.

Analysis Results

Compiled with gfortran 15.3 using -O3 -mtune=native -march=native.

Interpreting the Data

1. The Critical Path is identical (24 cycles). This is expected because the underlying math is the same. Both kernels perform:

• Three additions: tmp1 = v(i+1, j) + v(i-1, j) + v(i, j+1) + v(i, j-1)
• One multiplication: tmp2 = b(i, j) * dx2
• One Fused Multiply-Add (FMA): tmp3 = 0.25_dp * tmp1 + 0.25_dp * tmp2
• Plus associated loads and stores.

2. The LCD is the "Ballgame." Jacobi's LCD is only 1 cycle (simply incrementing the loop counter). Gauss-Seidel's LCD is 12 cycles. OSACA pinpointed the exact instruction:

46 | 12.0 | vaddsd %xmm1, %xmm0, %xmm0 | [46, 47, 48]

This is the "west-neighbor" term. The value u(i, j) written in one trip is read as u(i-1, j) in the next. The 12-cycle latency comes from three chained instructions (west-neighbor add $\rightarrow$ north-neighbor add $\rightarrow$ final multiply), each costing 4 cycles.

The compiler is being lazy. The compiler is being honest. The algorithm itself demands that $u(i, j)$ cannot be computed until $u(i-1, j)$ is finalized.

The difference is also visible in the assembly mnemonics: Jacobi uses packed double-precision instructions (SIMD), while Gauss-Seidel is forced into scalar operations.