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.
The algebraic manipulation that you use to break the dependency seems to be specific to the matrix being Poisson on a finite-difference grid (is that right?). In that case wouldn’t it make more sense to compare against specialized Poisson solvers?
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.
I doubt outside of a few exceptional cases that one is going to do better on CPU-bound problems than a well-written LAPACK implementation built for the architecture you intend to run on.
Maybe pedagogy was the point? I didn’t really get that from the article, maybe some intention was lost by filtering it through AI.
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. 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. 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. 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.
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.
This article might be good but I always get scared when they use AI art for the header bc I wont know if I'm reading AI slop until i've invested a good amount of time into the article so I just avoid it all together.
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).
"... and this is the whole ballgame: ", "jump out immediately", "... and the idea is disarmingly simple", etc etc.
Written by Claude AI. Which is a problem because Claude text is riddled with engagement bait and filler, making it exhausting to read.
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.
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
Artificial intelligence is a lot like alternative medicine, if it was the thing it said it was it wouldn’t need a name.
The algebraic manipulation that you use to break the dependency seems to be specific to the matrix being Poisson on a finite-difference grid (is that right?). In that case wouldn’t it make more sense to compare against specialized Poisson solvers?
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.
It looks like this analysis extends to the successive over-relaxation (SOR) method. Does it extend to the symmetric successive over-relaxation method?
Author here. Sure it does, and we'll get there! I just need to be done with my teaching duties a bit.
I doubt outside of a few exceptional cases that one is going to do better on CPU-bound problems than a well-written LAPACK implementation built for the architecture you intend to run on.
Maybe pedagogy was the point? I didn’t really get that from the article, maybe some intention was lost by filtering it through AI.
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.
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.
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.
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.
What language is that, and what does `0.25_dp` mean?
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.
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.
I wondered the same thing. _dp might mean double precision, but the rest is beyond me
Fortran
This article might be good but I always get scared when they use AI art for the header bc I wont know if I'm reading AI slop until i've invested a good amount of time into the article so I just avoid it all together.
It seemed to be ai assisted writing, but it was still interesting to me.
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).