loiseaujc

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.

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.

Author here. Glad you liked it. Claude has been used to help me structure the post, but not the writing itself (except maybe proof-reading it). All the techniques I use are pretty standard tricks to retain someone's attention (in particular easily-bored students in class).

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.

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.

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.

OP here.

> This is a big red flag. If they have little background in programming, teach good standards of Python.

That is a big issue, I'll admit it. And actually, in my Uni, I am a very strong proponent of making an intro to programming class mandatory before even being able to enroll in any other engineering classes. Unfortunately, I don't see this happening anytime soon mostly because the higher-ups have a distorted view of what programming actually is. And I believe it is the same in many French universities. May-be different elsewhere, I don't know, but in the mean time I have to make do with what I have.

> Otherwise they will get worst patterns, and think that GOTO is the preferred way of coding.

People keep on referring to `goto`. But that is a construct that has been considered bad since the Fortran 1990 standard, 35 years ago. Fair enough, there are plenty of legacy codes using it. Just like there are plenty of C code written 35 years which are just terrible by today standard. Modern Fortran (and by modern I mean anything following 2003 standard and more recent) does not use `goto`. And academics writing Fortran code today do not use `goto`, nor do they teach it.

> If you think why a lot of academic code looks like mess, one of the reason is using archaic tools and techniques as a standard.

Again, I think this is a very distorted view. When you need to run your code on 1000+ CPU, which is what I and my colleagues do, you better make sure your code makes use of all the most recent good practices and use industry-standard tooling not to waste CPU hours because of a stupid bug or what not.

OP here. 0-based vs 1-based indexing is one these hot debates I don't want get in as I've mentioned in the post. Yet, point is, if you take pretty much any math textbook, indexing starts at 1, whether we like it or not. And for many students, understanding the mechanics of a given algorithm is already enough of an effort that they do not need, on top of that, to translate every indices in the book's pseudo-code from 1-based to 0-based. That's all I'm saying. No other value judgement.