It seemed to be ai assisted writing, but it was still interesting to me.

I wondered the same thing. _dp might mean double precision, but the rest is beyond me

Fortran

Fortran’s types have kinds, which are integer codes that further distinguish various flavors of each type. The ISO standard failed to define their meanings for the built-in types, but with any modern compiler worth using the kind values are the byte sizes. Literal constants can have a kind value as a suffix, and the suffix can be either an integer or a named constant. Fortran codes tend to establish named constants for the “working precision” kind and the next higher precision kind, so that one need only redefine those to compile various versions of a subroutine that work with different precisions.

Even when one isn’t defining named constants for various kinds, it’s still good practice to use them on literal floating constants, since those default to the default kind of REAL in Fortran, not to double precision as in most later languages. A very pernicious pitfall of the language is to (say) assign or initialize X=0.1 instead of X=0.1_8 when X is double precision.

Author here. It is Fortran indeed, but there is nothing specific to Fortran in what I discuss. Had you written the kernel in C or C++, you'd run in exactly the same problems.

Artificial intelligence is a lot like alternative medicine, if it was the thing it said it was it wouldn’t need a name.

Damn yall are too exhausted to even read now because you spent so much mental energy trying to pattern match llm-ims instead of reading i guess

Author here. I did use Claude to help me structure the architecture of the post a bit, but not the text. These are idioms I use all the time as an educator to keep easily-bored students engaged as much as possible. Really sorry that this triggers you.

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. 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.

Author here. Sure it does, and we'll get there! I just need to be done with my teaching duties a bit.

Author here. It is indeed specific to this 2D Poisson matrix, although loop-unrolling + substitution is a very general technique in performance engineering. Eventually, whenever I'm done teaching classes, we'll get to specialized sparse matrix solvers for this problem and we'll see they come with their own set of problems (and loop carried dependencies are definitely one of them).

In any case, the idea behind this series on Jacobi (this is like the 4th post I believe) is to help students (and people interested in scientific computing in general) get what actually matters computationally and lift the hood to get some understanding of the interplay between such algorithms and their actual implementation on hardware.