Technical blog: ill-Conditioned Opinions

Uncategorized
1

Hej,

I’ve just published (yet) another post in the Jacobi vs Gauss-Seidel series. You can find it here. Moreover, rather than creating a new thread every time I publish something, I figured I might as well create just one and update the first message each time. At least, discussions and news will be somewhat centralized.

ill-Conditioned Opinions

Notes from the space between fields — iterative solvers, Krylov subspaces, reduced-order models, and the occasional strongly-held opinion on numerical methods. Written for people who’d rather understand why something breaks than just get it to work. Fortran included, apologies not.

Current series: From Jacobi to multigrid

  • Is Fortran better than Python for teaching the basics of numerical linear algebra?
    Date: September, 9th 2025
    TL-DR Not about speed or elegance. It’s about teaching. After 10+ years running an intro scientific computing class in Python/numpy, I argue that modern Fortran’s strong typing, 1-based indexing, and explicit do/end do blocks may actually make it a better first language for numerical linear algebra than Python, removing incidental cognitive load (indentation errors, off-by-one bugs, numpy’s many dot-product syntaxes) so students can focus on the algorithms. Walks through two classic examples (Jacobi solver, QR least-squares) comparing typical student code in both languages. Not a call to ditch Python, just a case for giving Fortran more credit in the classroom.

  • Jacobi method: From a naïve implementation to a modern Fortran multithreaded one
    Date: September, 23rd 2025
    TL-DR Covers the theory (convergence conditions, spectral radius) behind the Jacobi method for the 2D Poisson equation, then optimizes a standard-compliant Fortran implementation step by step: eliminating a redundant array copy, switching to do concurrent, and reducing unnecessary residual-norm checks, taking the baseline from 58s down to 16s. Finishes by multithreading via -ftree-parallelize-loops=n alone (no OpenMP/MPI), reaching a 27x speedup on a quarter-million-unknown system in under 3 seconds with the exact same source code.

  • Gauss-Seidel is faster than Jacobi. Or is it?
    Date: September, 8th 2026
    TL-DR Gauss-Seidel needs half the iterations of Jacobi to converge on the 2D Poisson test problem. But implemented in Fortran, it ends up 5x slower than the optimized Jacobi solver. The culprit isn’t the compiler but a loop-carried dependence: lexicographic Gauss-Seidel updates points sequentially in-place, blocking instruction-level-parallelism and vectorization, while Jacobi’s two-buffer scheme has no such dependency and lets do concurrent fully vectorize. Wraps up pointing toward a way to recover Gauss-Seidel’s convergence advantage without losing SIMD-friendliness (covered in the next post).

  • Make Gauss-Seidel great again!
    Date: September, 14th 2026
    TL-DR Follow-up to the previous post’s puzzle: Gauss-Seidel is slow not because of a lack of vectorization but a 12-cycle loop-carried dependency (confirmed via OSACA static analysis of the compiled assembly). The fix requires algebraically rewriting the recursion (solving two steps ahead) so consecutive updates become independent functions of prior data. Implementing this degree-2 “unrolled” Gauss-Seidel kernel in Fortran drops the LCD from 12 to 4 cycles, nearly matching Jacobi’s per-sweep speed while keeping Gauss-Seidel’s 2x faster convergence, a genuine single-core win. The catch: the compiler achieves this via independent scalar instruction chains, not true SIMD vectorization, so the kernel still can’t be multithreaded.

  • :tada: :tada: :tada: Red or Black?
    Date: September, 21st 2026
    TL-DR Last post’s “unrolled” Gauss-Seidel kernel was a serial-only win. It still can’t be vectorized or multithreaded because the compiler achieves speed via scalar instruction-level parallelism, not SIMD. This post switches strategy entirely, using red-black coloring to eliminate loop-carried dependencies at the algorithmic level rather than patching around them. This is shown mathematically to be equivalent to block Gauss-Seidel on a permuted, still-SPD system, so convergence is preserved. The resulting Fortran kernel gets true SIMD vectorization (confirmed via OSACA) and scales with OpenMP threading, beating Jacobi at every thread count up to 8 cores.

8 Likes