Lesson 12 — A complete iterative solver

Lesson 11 ran a fixed number of Jacobi sweeps. A solver instead runs until the answer is good enough, and every PE must agree on when that is. This lesson adds a sum over the whole grid after every sweep, and a loop that stops when the sum is small enough or after a maximum number of sweeps. It puts together what the main sequence has taught: a kernel with functions, a halo exchange, a reduction and a broadcast, and a layout written as Paint functions.

What you will be able to do

• Combine a halo exchange, a grid-wide sum and a convergence test into a program that stops by itself.
• Add one value from every PE of a grid and give the total to every PE: along the rows, up one column, and back by multicast.
• Give each PE its role at run time from its row and column ID registers, so that one code() covers the whole grid.
• Pass run-time parameters in a struct, and write the struct's input file from its layout in memory.
• Measure the three phases of a sweep, and tell which one grows with the size of the grid.

The program

The program solves the Poisson equation −∇²u = f on a field of H·B rows and W·B columns, with u = 0 outside the field. One Jacobi sweep replaces every point with

u_new[i][j] = 0.25 * (u[i-1][j] + u[i+1][j] + u[i][j-1] + u[i][j+1] + f[i][j])

f is the source term; the test sets it to 1 at every point. The sweeps start from random integers between 0 and 63 and approach the solution: a field that is highest in the middle and 0 at the edges, like a plate heated evenly and cooled at its rim.

After each sweep, every PE sums the squares of the changes of its block. The sum of these over all PEs is the residual. The kernel's loop is:

while residual > tol2 and sweeps < maxit:
    exchange the block edges
    update the block
    residual = the sum, over all PEs, of the squared changes

Every PE records the final residual, the number of sweeps and why it stopped. The test checks that:

• every PE's u equals a NumPy computation in float32, bit for bit;
• all PEs report the same residual, number of sweeps and stop reason;
• the number of sweeps and the stop reason are the expected ones;
• the residual is within 1e-4, relative, of the reference's.

The layout is lesson 11's, with five more colors for the sum. Partial sums move right along every row on rr, then up the last column on cc. The PE at the top right, marked *, then holds the total. It sends it on g, down the last column and left along every row. The figure shows a 3 × 3 grid:

  The sums, on rr and cc              The total, on g

y=3  [ ]-rr->[ ]-rr->[*]              [ ]<-g-[ ]<-g-[*]
                      ^                               |
                      cc                              g
                      |                               v
y=2  [ ]-rr->[ ]-rr->[ ]              [ ]<-g-[ ]<-g-[ ]
                      ^                               |
                      cc                              g
                      |                               v
y=1  [ ]-rr->[ ]-rr->[ ]              [ ]<-g-[ ]<-g-[ ]
     x=1     x=2     x=3              x=1    x=2    x=3

Row	B	W × H	K	maxit	Expected
tutorial	4	3 × 3	8	50	converged after 8 sweeps
w4h3	6	4 × 3	3	50	converged after 3 sweeps
w4h2	6	4 × 2	4	50	converged after 4 sweeps
w2h2	4	2 × 2	3	50	converged after 3 sweeps
maxit	4	3 × 3	8	5	stopped at maxit after 5 sweeps

How the test knows the number of sweeps

ref.py runs the same sweeps in NumPy and keeps the residual of each in the list res, so res[k - 1] is the residual of sweep K. It then sets the tolerance between the residuals of sweeps K − 1 and K, at their geometric mean:

tol2 = np.float32(np.sqrt(res[k - 2] * res[k - 1]))

The reference converges after exactly K sweeps. It adds the squares exactly, with Python's math.fsum; the kernel adds them in sp, in an order Carbide chooses. So the two residuals differ in the last bits. Both residuals sit far from tol2: in the tutorial row, sweep 7 gives 2886.59 and sweep 8 gives 2574.73, about 6 % above and below tol2 = 2726.2. So the kernel also stops after K sweeps, and u can be compared bit for bit. ref.py stops with an error if either residual lies within 1 % of tol2.

The code

The parameter struct

struct params_t {
    sp tol2;
    xp maxit;
};
struct params_t params;    // written by miset

A struct groups named fields, as in C, and params.tol2 reads one of them. The map shows how params lies in memory:

grid.params                  [ y<3 x<3 n<1 o<4         ] ( 0x2+o              1+x 1+y )

Each PE has one struct (n<1) of four words (o<4). tol2 takes words 0 and 1, maxit word 2, and word 3 is unused. An sp field starts at an even word: with a third field, sp f;, the struct takes six words, and f takes words 4 and 5.

The Makefile loads the struct like any array:

	miset image=$@ out=$@ map=bench.map symbol=grid.params order="y x n o" data=params.bin

ref.py writes params.bin with a NumPy type that copies this layout, the unused word included, one struct per PE:

params = np.zeros(w * h, dtype=[("tol2", "<f4"), ("maxit", "<i2"), ("unused", "<i2")])

Roles from the ID registers

xp row_id[1] address(0x7c45);
xp col_id[1] address(0x7c46);

    xp x = col_id[0];
    xp y = row_id[0];

Paint writes every PE's position into two registers of that PE: Paint's x into the column ID register and Paint's y into the row ID register. The positions are those of the whole layout. The grid starts at x = 1 and y = 1, because the shore takes x = 0 and y = 0, so its PEs read x from 1 to W and y from 1 to H. The map's x and y count the copies from 0 instead: the map's x = 0 is the PE whose col_id is 1.

Lesson 9 gave each role along a row its own kernel and its own code(). Here one kernel compares x and y with 1, W and H, and takes its role at run time. One code() covers every PE, name('grid) names all of them, and miget symbol=grid.resid reads every PE's value.

The allreduce

function allreduce(sp s, xp x, xp y) -> sp {
    if (x > 1) {
        s = s + rr[];
    }
    if (x < W) {
        rr[] <- s;
    } else {
        if (y > 1) {
            s = s + cc[];
        }
        if (y < H) {
            cc[] <- s;
        } else {
            g[] <- s;
        }
    }
    if (x < W || y < H) {
        s = g[];
    }
    return s;
}

s comes in as the PE's own sum and leaves as the grid's total. A scalar argument is a copy (lesson 2), so the function can change s freely. What each PE does:

PE	Receives	Sends
first column, x = 1		its sum on rr
1 < x < W	the sum from the left on rr, adds its own	on rr
x = W, y = 1	from the left on rr	the row's sum on cc
x = W, 1 < y < H	from the left on rr and from below on cc	on cc
x = W, y = H	the same	the total on g
every PE but (W, H)	the total on g

Each receive and send here runs on the main thread, so a PE waits at rr[] until its left neighbour has sent. The PE at (W, H) keeps the total it sent. || evaluates both operands, as in lesson 2; both are plain comparisons, so that is harmless.

The exchange and the update

exchange() is lesson 11's parallel block, moved into a function. It returns nothing, so its definition, function exchange() { ... }, has no -> type, and a call to it is a statement of its own. update() is lesson 11's update with two additions: the source term, and the sum of the squared changes, which it returns.

    ∀p ∀q v[p][q] <+- f[p][q];
    ∀p ∀q du[p][q] <- u[p+1][q+1];
    ∀p ∀q u[p+1][q+1] <- 0.25:sp * v[p][q];
    ∀p ∀q du[p][q] <- u[p+1][q+1] - du[p][q];
    sp s = 0.0:sp;
    ∀p ∀q s <+- du[p][q] * du[p][q];
    return s;

du keeps the old block, then the change. s <+- du[p][q] * du[p][q] adds B² products into one scalar. Carbide uses the dot-product instruction (fdps in the listing) and may add the terms in any order, which is why the residual check allows a tolerance. f is an array loaded by miset, like u.

The update takes 154 cycles at B = 4, against 41 in lesson 11, and the sum is most of the difference. Each of the other statements is one instruction over the whole block, as in lesson 11. Over a 2-D array, Carbide makes the sum a loop with one pass per row of the block: an fdps over the row's B values, then an fstdpas and an fadds that add the row's result to s. In a test kernel with the update's statements and B = 4, this loop took 86 cycles, and the other three statements that lesson 11 does not have, v += f and the two du lines, took 30 together. In a copy of this lesson with du declared with one dimension, sp du[i=B*B];, written as du[p*B + q], and summed with let i in [0, B*B); and ∀i s <+- du[i] * du[i];, the sum is a single fdps, and the update took 104 cycles at B = 4.

Neither function takes the arrays or sockets as arguments: they use the kernel's own u, f, v, du and sockets.

The loop

    sp res = 1.0e30;
    xp it = 0;
    tsc_control[0] <- tsc_control[0] | 4;
    while (res > params.tol2 && it < params.maxit) {
        stamp[0] <- tsc[0];
        exchange();
        stamp[1] <- tsc[0];
        sp s = update();
        stamp[2] <- tsc[0];
        res = allreduce(s, x, y);
        stamp[3] <- tsc[0];
        it = it + 1;
    }

• res starts at 10³⁰. In every test row tol2 is below that and maxit is positive, so the loop runs at least once. A literal with an exponent is always sp; adding :sp to it is a syntax error.
• && evaluates both comparisons, which is harmless here too.
• The test uses res, the total from the allreduce. Every PE receives the same 32 bits, so every PE leaves the loop after the same sweep. See Pitfalls for a test on the PE's own sum.
• The four timestamps bracket the three phases of each sweep, as in lesson 11. cycles.py prints the last sweep's phases.

After the loop, the kernel records why it stopped:

    if (res > params.tol2) {
        stop[0] = 2;
    } else {
        stop[0] = 1;
    }

The loop ends when either test fails. If the residual is still above tol2, the count reached maxit.

The layout: bench.paint

beach, shore and blocks are lesson 11's functions; blocks now places solver.w and passes it W and H. A fourth function routes the allreduce:

define reduction(w h) {
    let rows = rect(1 1 w h)
    rows : paint(rrp [CXR|R])
    rows : paint(rr  [L|LXC])
    rows : colorpair(rr rrp)

    let last = rect(w 1 1 h)
    last : paint(ccp [CXU|U])
    last : paint(cc  [D|DXC])
    last : colorpair(cc ccp)

    let top = rect(w h 1 1)
    top : paint(g [C>DL])
    let side = rect(w 2 1 h-2)
    side : paint(g [U>CDL])
    let bottom = rect(w 1 1 1)
    bottom : paint(g [U>CL])
    let middle = rect(2 1 w-2 h)
    middle : paint(g [R>CL])
    let first = rect(1 1 1 h)
    first : paint(g [R>C])

    union(rows last top side bottom middle first)
}

• rr and rrp are lesson 11's rightward colorpair, on the grid only. The |R term stops the last column from sending further right.
• cc and ccp are its upward colorpair, on the last column only.
• g is a plain color with one rule for each part of the grid. The total starts in the core of the PE at (W, H) and goes down and left. The rest of the last column takes it from above, gives it to its core and passes it down and left; its bottom PE does not pass it down. The middle columns take it from the right, give it to their cores and pass it left. The first column only gives it to its cores. One wavelet reaches every PE, because the routers copy it.

Each colorpair with | terms is painted on its own rect, as lesson 11 requires. With W = 2, middle is an empty rectangle, and with H = 2 so is side; Paint accepts both. W = 1 or H = 1 would put two g rules on one PE, so the layout needs W ≥ 2 and H ≥ 2.

The last line joins the four functions:

union(beach(W+2 H+2) shore(W+2 H+2) reduction(W H) blocks(W H))

What the hardware does

Three phases, measured

Cycles of the last sweep. A range is the smallest and largest value over the PEs.

Row	B	W × H	Exchange	Update	Allreduce
w2h2	4	2 × 2	94-100	154	83-89
tutorial	4	3 × 3	91-100	154	125-134
w4h2	6	4 × 2	115-116	322	120-121
w4h3	6	4 × 3	115-116	322	143-144

• The exchange and the update depend only on B. Each PE does the same work however large the grid is.
• The allreduce grows with the grid. The partial sum passes through W + H − 1 PEs one after another, and each one receives, adds and sends on its main thread. Between rows with the same B, each PE added to the path adds about 20 cycles: from w2h2 to tutorial (B = 4) the path grows by two PEs and the allreduce by 42-45 cycles, and from w4h2 to w4h3 (B = 6) by one PE and 23 cycles.
• The total comes back as one wavelet that the routers copy to every PE; no core forwards it.

So the exchange and the update scale weakly, as in lesson 11, and the allreduce does not. A global decision costs a trip across the grid.

The budget

Resource	solver.w uses	A PE has
Colors	13: eight for the halo, rr and rrp, cc and ccp, g	24
Input queues	7, none doubled	8
Output queues	7	8
General microthreads	8, all in the exchange	8
Memory	2(B+2)² + 6B² words for u, f, v and du	24,576 words (48 KB), shared with the code

The listing (test-tutorial/solver_artefacts/solver2.yeti) shows the queues Carbide chose: the halo sockets a, b, c and d on queues 0 to 3, and rr, cc and g on queues 4 to 6, all with .doubling=0. Queues 0 and 1 hold 8 wavelets, the others 4. The six allreduce transfers have .upr=NONE: they run on the main thread, not as microthreads.

Run it

In your copy of the tutorial (see Lesson 1):

cd tutorial/main-4-grids/12-iterative-solver
make TEST=tutorial
cat test-tutorial/make-log

make TEST='*' runs all five rows. The make-log of a passing row starts with the reference's residuals and ends with the kernel's result and the phases:

residual of each sweep: 56990.3 13767.2 7451.68 5114.77 3977.16 3318.9 2886.59 2574.73
tol2 2726.2, maxit 50: expect 8 sweeps, stop 1
...
9 PEs: residual 2574.73 (reference 2574.73), 8 sweeps, stop 1
last sweep, exchange : 91-100 cycles over the PEs
last sweep, update   : 154-154 cycles over the PEs
last sweep, allreduce: 125-134 cycles over the PEs

The test does not check the cycle counts. A failing row prints FAIL. The make-log then shows either cmp's first difference in u, or one of check.py's messages: the PEs disagree, wrong number of sweeps or wrong stop reason, or residual differs from the reference.

After you edit a source file, run make clean before the next test.

Pitfalls

• Test the total, not the PE's own sum. With res = s; in place of res = allreduce(s, x, y); (and the allreduce still called), each PE compares its own share with a tolerance meant for the whole grid. In the tutorial row all nine PEs stopped after 2 sweeps instead of 8, and cmp failed. With other data, PEs could leave the loop after different sweeps. A PE that has left stops sending, so its neighbours would wait forever in exchange().
• Read a struct field into a local before a ∀ statement. In a test kernel whose struct had an extra field sp f, ∀p ∀q v[p][q] <+- params.f took 365 cycles on an 8 × 8 block: Carbide computed one element at a time and loaded the field again for each. With sp c = params.f; first, ∀p ∀q v[p][q] <+- c took 18.
• The roles assume the layout. The kernel takes x = 1 as the first column because the grid starts at (1, 1). If you move the grid, the ID registers change with it, and the comparisons in allreduce must change too, or take the corner as parameters.
• Every PE must call allreduce equally often. A PE that skips one leaves its neighbours waiting for a partial sum that never comes.
• W and H must be at least 2, because of the g rules (see The layout).

Exercise

The allreduce is the one phase that grows with the grid. Test convergence only every fourth sweep, and measure what that saves:

1. Time the whole loop first, before you change it. In solver.w, declare u16 whole[k=2];, and write tsc[0] into whole[0] before the while loop and into whole[1] after it. Add this line to the Makefile's stamp-%.bin rule:

   	miget image=$< map=bench.map symbol=grid.whole order="y x k" data=whole-$*.bin

   Run make clean and make TEST=tutorial. Then this prints the time of the loop on each PE, whole[1] - whole[0] modulo 65536:

   python3 -c "import numpy as np; w = np.fromfile('test-tutorial/whole-sim.bin', dtype=np.uint16).astype(int).reshape(-1, 2); print((w[:, 1] - w[:, 0]) % 65536)"

2. In solver.w, declare let j in [0, 4); before the loop. In the while loop, run the sweeps as a loop, ∀j { exchange(); s = update(); }, with sp s = 0.0:sp; before it. Call allreduce once after it, and count it up by 4. Write the four sweeps as a loop: four written-out copies of the two calls made each exchange and update slower, and hid most of the saving.

3. In ref.py, set tol2 between the residuals of sweeps K − 4 and K: use res[k - 5] in place of res[k - 2], in the guard too.

4. In test.rc, make every K a multiple of 4 and at least 8, and give the maxit row a MAXIT that is a multiple of 4 below K, such as K=8 MAXIT=4. The loop can stop only after a multiple of 4 sweeps. The tutorial row keeps K = 8.

5. Run make clean and make TEST='*', then print the loop times of the tutorial row again, and compare them with step 1's.

Then estimate, from the table above, the share of each sweep that the allreduce takes on a grid of 100 × 100 PEs with B = 6.

Next

This completes the main sequence. The extension topics build on it; each names the lessons it needs. For grids, see Physics on a grid of points for other boundary conditions, and Reductions in the router: the fabric adder for a sum that the routers compute. The tutorial's front page lists them all.

References

• Structs
• Functions
• While loops
• Hardware register access (address())
• Dot product
• Typed literals
• Function definitions (define)
• Static color routing