The problem is to solve a set of linear equations with one unknown for each point in a provided mesh. There is a symmetric adjacency relation among the points. For each point there is one equation expressing the unknown for that point as a constant plus a linear combination of the unknowns for the adjacent points.
(When the mesh is a line of points there is a particular solution that is especially easy, simple and fast. You should not study the solution here unless you are curious or really need it. Even then the 1D case is a good intro.)
The points are numbered consecutively and the numbering has nothing to do with logic of most of the code.
For me this code is an unhappy mixture of patterns from the world of mutable arrays and for (do) loops, and the world of Lisp and immutable lists. It could be well transcribed to C but better debugged here in the Scheme world.
Several of the routines take or return a block which is a Scheme pair of a list of dimensions and an array of lc’s. ((5 8) . vector-of-lcs) would be a 5 by 8 rectangular mesh with an array of 40 elements representing 40 points. Most of the current code is either oblivious to the dimensions or assumes that there are exactly two dimensions.
Of the procedures named in the top level Scheme name space just retire and bs modify the lc array within their block argument. One must keep this in mind to follow the program logic. The routine retire takes an array and an index into the array and eliminates that unknown from the set of equations. These concepts are taken directly from the first lecture on simultaneous equations from a high school algebra course. The routine replaces the point state by another in the form ("d" . lc) where lc is in the original form. The rules are also relaxed to allow the new expression to refer to non-adjacent points. The meaning, however, remains the same as above. Point states of unretired points will remain in the first form but will be replaced by states that do not mention the retired point.
After all of the points have been retired, the routine bs (back substitution) is invoked for each element of the array in the inverse order that the elements were retired. These invocations modify only the one array element, replacing it by one of the form ("k" . num) where the number within is the final answer that particular unknown.
The routine solve invokes retire and bs in the most obvious manner, which is adequate for many problems including the demonstrations below.
The code vg generates an array for a rectangular mesh where only nearest neighbors are adjacent. At the top and bottom edges there are imaginary neighbors with values 1 whereas at the right and left edges there are such neighbors with value –1. Within the mesh proper points are required to be the average of their four neighbors, real and imaginary. This is the meaning of the Poisson equation. It is curious and profligate to solve this problem with Scheme’s rational arithmetic. A 10 by 10 mesh is solved in about 30 seconds on a 400 MHz PowerPC. Intermediate values are rational fractions with numerators and denominators of 25 digits and more. The n by n mesh problem requires n4 floating point operations but requires n6 steps with rational arithmetic.
A strong test of the results is to note that the solution has the same symmetry as that of the posed problem. Routine flipH flips the array horizontally and flipV vertically. flipD is a diagonal flip that also flips the sign. cmpv reports the maximum absolute difference between two arrays. The line of code (cmpv s t) compares the exact solution with the one produced with double precision floating point arithmetic. All output is less than 2–50 which is remarkable given that IEEE floating point provides 53 bits and that many of the value magnitudes are near 1.
The routine verify takes the cheat sheet, (known answers) and an array that is part way thru the sequence of steps towards those answers, and verifies that the developing theory of the solving code is still in conformance with the truth.
Bugs were found in the following hierarchy:
Here is some speculation on how to do this in multiprocessor system.