LAPACK doesn't do sparse systems.

(On the other hand, I skimmed the article and this might be a banded system, which LAPACK can handle).

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.

▲

sheafification 2 days ago | parent | prev [-]

True, 2D Poisson is only approximately banded. Still, I don’t think I’d bother rolling my own sparse solver. It’s a well-trod problem.

	▲	loiseaujc a day ago \| parent [-]
		Author here. It actually is a block tridiagonal matrix with tridiagonal blocks. The test problem only has a quarter millions of unknowns so, for a production run, I wouldn't bother write my own sparse solver either. Here, it is done mainly for the sake of pedagogy to help students understand the interplay between a mathematical algorithm that looks good on paper and its hardware implementation which is not as promising as one would expect. Jacobi and Gauss-Seidel wouldn't be solvers I'd even consider for a real problem. But they are simple enough that anyone with a basic understanding of linear algebra and programming can follow along. But much research requires me to run simulations on thousands (if not hundreds of thousands) of cores, and there, off-the-shelf solver implementations will often not cut it.