Remarks on some Implicit Mesh Code

I have not found a way to package code on the web to make it convenient both for study and execution. This Scheme code solves a set of linear equations defined on a mesh. The connectivity of the mesh is defined by input to the code. The code leaves hooks for topology savvy code to guide the order of computation. Great economy can result from this on some meshes.

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.

The Mesh

The soul of the mesh is a mutable array with an element per mesh point. (Arrays are called “vectors” in Scheme.) I call the elements of the array point states here. The initial state of the array is the formulation of the problem and our algorithm modifies the array until finally, point by point, the array explicitly provides the answers. There are three forms of point state and external logic is presumed to know which form a particular point state is in, although there is considerable internal logic to check that these states are not misused. Each point state form can be taken as an assertion, in the form of a linear equation, relating the unknown at the point to other unknowns. A mesh problem is posed by creating such an array with all point states relating that unknown as a linear combination of adjacent unknowns. Such a state is coded, for example, as (1.2 (2 . 1/4) (5 . 1/4) (6 . 1/8) (n . 0)). If this state is found in slot 3 of an array of length n then the meaning is that
T3 = 1.2 + T2/4 + T5/4 + T6/8. The last member of the Scheme value that codes this, (n . 0), is a hack that greatly simplifies a critical complex inner loop. Sorry. This list must be sorted by index. We refer to data in this form as an lc datum. The variable on the left must not be mentioned on the right. The debugging routine ln checks an array for compliance with these rules.

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.

Miscellaneous Code Details

The code exit and mer are taken from here where some commentary may be found. We use merge-assim and ms to sort things. They are taken from here. Do is the dumbest do loop construct that I could think of as I have not learned the standard elaborate Scheme do construct. (Do n proc) invokes the single parameter procedure proc once for each of the integers from 0 to n−1.

The routine solve invokes retire and bs in the most obvious manner, which is adequate for many problems including the demonstrations below.

Demonstrations

Most of the bugs were found by such trivially posed problems as #((2 (1 . 3) (2 . 0)) (5 (0 . 2)(2 . 0))) or #((1 (1 . 1/2) (3 . 0)) (2 (0 . 3/10) (2 . 3/10) (3 . 0)) (1 (1 . 1/2) (3 . 0))). Two bugs remained and the non-trivial problem below revealed them.

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.

Debugging the Code

Here is the code to debug the code. Debugging this messy code was greatly simplified by the presence of a strong theory regarding the meanings of the intermediate states. A debug jig could know the answers that the tested code was to find and check at each step of the solution that the meanings were indeed satisfied by these known (to the tester) solutions. The tester knows the solutions by running the retrodiction routine tfn which takes a posed problem and changes the constants to make the solution match an array argument to tfn.

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:

This leaves only termination of the procedure which is easy in this case as it involves retiring each element and then back substituting each element in the opposite order. The routine Solve solves the problem while applying intermediate validity testing.

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