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.

T

Several of the routines take or return a ` block` which is a Scheme pair of a list of dimensions and an array of

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

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

The routine ` solve` invokes

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 n^{4} floating point operations but requires n^{6} 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:

- Scheme reveals that the program is not a valid Scheme program.
`ln`reveals that the problem is ill-posed, or the information becomes malformed.`verify`reveals that some step, including the originally posed problem, is not in conformance with the known truth.

Here is some speculation on how to do this in multiprocessor system.