Author here. The second-order accurate finite difference approximation of the 2D Poisson operator is a block tridiagonal matrix with tridiagonal blocks. There are ways to handle it with LAPACK, but you'll still need special care. And 512x512 grid points (quarter millions unknowns) really is a toy problem. The linear systems I work with typically have in the billions of unknows, and LAPACK is not gonna cut it no matter how smart you are.

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.