Author here. Performance-wise, it will never ever be a serious solver. No matter how hard you try. And for this particular matrix, (2D Poisson, block tridiagonal with tridiagonal blocks, quarter millions unknows), LAPACK isn't gonna do it either unless you take special care.

Pedagogy is indeed the point. This is the 4th (3rd?) post in this Jacobi series. Gauss-Seidel and Jacobi are sufficiently simple that very little math is needed to understand how they work. But they also offer a very nice opportunity to illustrate the interplay of mathematical algorithms and actual hardware implementation as well as gradually introducing some elements of performance engineering.

To give you idea of where we're going: after discussing SOR and make connections with other fields of applied math, we'll eventually get to multigrid solvers. And Jacobi and Gauss-Seidel play a really important role in this. Multigrid are also some of the best iterative solvers possible for such discretized elliptic operators. But I can only blog now and then, and I want to take one step at a time so that upper-level undergradute in either math or computer science can follow along.