Lesson 8 — Space and time: a wavefront along a row

So far every PE had its own code() call and its own route rules. This lesson runs one kernel on a whole row of PEs, and spreads a computation over the row in space and over steps in time. The program aligns two DNA sequences with the Smith-Waterman algorithm: each PE owns one column of the score matrix, and the rows of the matrix stream through the PEs as a wavefront.

What you will be able to do

• Place one kernel on a row of PEs with a single code(), and give each PE its own data with miset.
• Route a row with a colorpair, so that every PE sends and receives on one socket under the same two route rules.
• Close the two ends of a row with moats: one that feeds the row, one that absorbs what leaves it.
• Map a recurrence onto PEs and steps, and check the mapping by reading the cycle at which each value was computed.

The program

Smith-Waterman scores how well parts of a query sequence (Q bases) match parts of a reference sequence (W bases). It fills a Q × W matrix H, row i for query base i and column j for reference base j:

H[i][j] = max(0,  H[i-1][j-1] + score(q[i], ref[j]),  H[i-1][j] - 1,  H[i][j-1] - 1)

score is +2 for a match and -1 for a mismatch; the two - 1 terms are gaps. Cells outside the matrix count as 0. The largest H[i][j] is the score of the best local alignment. Bases are coded A=1, C=2, G=3, T=4, and 0 pads a short query and never matches.

The row has W + 2 PEs. PE j holds reference base j. A moat is a PE at the edge of the computing region that supplies or absorbs its boundary data: the left moat holds the query, and the right moat absorbs what the last PE sends.

  x=0           x=1        x=2              x=W          x=W+1
+-----------+  +-------+  +-------+       +---------+  +------------+
| left_moat |->| PE 0  |->| PE 1  |-> .. ->| PE W-1  |->| right_moat |
| query[]   |  | ref 0 |  | ref 1 |       | ref W-1 |  | absorbs    |
+-----------+  +-------+  +-------+       +---------+  +------------+
     every PE, moats included:  rp [CXR|R]  and  r [L|LXC]

Each arrow carries two hp values per step on one socket, r. The test compares every cell of H with a Python reference (ref.py). The three rows of test.rc:

Row	W	Q	Sequences	Largest H
tutorial	32	16	ACGT repeated: the query is part of the reference	32, at H[15][15]
snp	32	16	a 14-base query with one mismatched base	25, at H[13][15]
long	8	40	ACGT repeated: the reference is part of the query	16, at H[7][7]

For snp, the best alignment has 13 matches and one mismatch: 13 × 2 - 1 = 25.

The mapping: a column per PE, a row per step

H[i][j] needs three neighbours: the cell above, H[i-1][j]; the cell to the left, H[i][j-1]; and the diagonal, H[i-1][j-1]. Cells with the same t = i + j form an anti-diagonal. Every cell on anti-diagonal t depends only on anti-diagonals t - 1 and t - 2, so all of its cells can be computed at once.

PE j computes column j, one row per iteration: iteration i computes H[i][j]. PE j can start a row only when PE j-1 has sent it its value for that row, so the PEs start one step apart, and at step t PE j works on row t - j:

         PE 0     PE 1     PE 2     PE 3
t = 0:   H[0][0]
t = 1:   H[1][0]   H[0][1]
t = 2:   H[2][0]   H[1][1]   H[0][2]
t = 3:   H[3][0]   H[2][1]   H[1][2]   H[0][3]

This moving front is the wavefront. Each dependency has a place to live:

Value	Computed by	How PE j gets it
H[i-1][j], above	PE j, one iteration earlier	kept in h_above
H[i][j-1], left	PE j-1	received on r
H[i-1][j-1], diagonal	PE j-1, one iteration earlier	the previous left value, kept in h_diag

The query base q[i] travels with the left value: PE j-1 forwards it along with H[i][j-1], so only the left moat stores the whole query. The matrix has W + Q - 1 anti-diagonals, so the run takes W + Q - 1 steps, not W × Q.

The code

The kernel, sw.w

Every one of the W PEs runs this file:

xp param Q = 16;

socket r;

xp ref_base[n=1];
hp history[n=Q];
u16 when[n=Q];

u16 tsc[1] address(0x7b30);
u16 tsc_control[1] address(0x7b34);

• = 16 gives Q a default, which a Q: binding in code() overrides. The layout passes Q:Q, the value from test.rc.
• ref_base holds this PE's reference base. It has no initializer: miset writes a different value into each PE's copy (below).
• history[i] keeps H[i][j], so the test reads the whole matrix from the PEs afterwards.
• when[i] keeps the timestamp at which the PE finished row i, read from the timestamp counter of lesson 3.
• r is the only socket: the PE receives on it and sends on it.

The program starts the counter and the two values a PE keeps between rows:

    tsc_control[0] <- tsc_control[0] | 4;

    hp h_above;
    hp h_diag;
    h_above <- 0:hp;
    h_diag  <- 0:hp;

Both start at 0, the matrix's row -1.

The loop body is one cell:

    let i in [0, Q);
    ∀i {
        hp h_left;
        hp q_char;
        h_left <- r[];
        q_char <- r[];

        // Code 0 pads a short query and matches nothing.
        hp score;
        if (q_char == (hp)ref_base[0] && q_char != 0:hp) {
            score <- 2:hp;
        } else {
            score <- -1:hp;
        }

        hp from_diag;  from_diag  <- h_diag  + score;
        hp from_above; from_above <- h_above - 1:hp;
        hp from_left;  from_left  <- h_left  - 1:hp;

        hp h_new;
        h_new <- 0:hp;
        if (from_diag  > h_new) h_new <- from_diag;
        if (from_above > h_new) h_new <- from_above;
        if (from_left  > h_new) h_new <- from_left;

        history[i] <- h_new;

        // Row i+1 needs H[i][j-1] as its diagonal and H[i][j] as above.
        h_diag  <- h_left;
        h_above <- h_new;

        r[] <- h_new;
        r[] <- q_char;
        when[i] <- tsc[0];
    }

The loop variable i is the query row. It is not the anti-diagonal t: the PE never needs t, because each iteration waits for its left neighbour's values.

Two values travel per step on one socket: H[i][j] first, then q[i]. Nothing in a wavelet says which is which, so the receiver reads them in the order they were sent.

The scores are small integers, which hp holds exactly up to 2048. A match adds 2, so no score exceeds 2 × min(Q, W).

The moats

The left moat, moat_left.w, sends what PE 0's missing left neighbour would send: H[i][-1] = 0 and the query base.

    let i in [0, Q);
    ∀i {
        r[] <- 0:hp;
        r[] <- (hp)query[i];
    }

The right moat, moat_right.w, reads every wavelet the last PE sends and discards it:

    forever dispatch(r) {
        wavelet() { }
    }

It is the handler with no arguments from lesson 7: it runs once per full data wavelet and needs no count. The last PE sends two plain hp values per step, and they share one wavelet, so the handler runs once per step: with a counter in it, it ran 16 times for Q = 16. Lesson 7 warns that a handler counts wavelets, not values; that does no harm here, because the handler only discards.

wavelet is the reserved name for the handler of plain data wavelets. A handler with any other name, such as edge() { }, is a named message (lesson 7), and plain sends never run it: with a counter in it, it counted 0, and nothing reported a problem. This moat still let the row finish, because the wavelets are taken off its queue whether or not a handler runs. A moat that must answer each wavelet, as in lesson 11, would never answer.

The right moat has to have code. Without a reader at the end of the row, the last PE's wavelets wait in the buffers of the router and the input queue, and once those are full its sends block, as in lesson 6. With no code on the right moat, the last PE finished 12 of its 16 rows.

One kernel on a rectangle

The layout, bench.paint, places sw.w on W PEs with one code() call:

let interior = rect(1 0 W 1 'j)
interior : code('sw [r:r Q:Q])
interior : name('grid)

rect(x y w h) is a w × h rectangle of PEs: columns x to x+w-1, rows y to y+h-1. One code() on it puts the same kernel, with the same parameters, on every PE of the rectangle. Paint builds sw.w for the whole rectangle, not once per PE: every PE runs the same compiled code, and only the data differ.

The quoted names after the size name the rectangle's dimensions, rect(x y w h 'xdim 'ydim). Here the x dimension is called j, after the matrix column. The map records it with every symbol of grid:

grid.ref_base                [ j<32 n<1          ] ( 0x0          1+j 0 )
grid.history                 [ j<32 n<16         ] ( 0x1+n        1+j 0 )

grid.history has dimensions j (32 PEs) and n (16 values). PE j sits at fabric column 1 + j, and its history[n] at word 0x1 + n. A rectangle dimension counts from the rectangle's own corner, not from column 0. Without a name, Paint calls the dimensions x and y, and leaves out a dimension of size 1: this rectangle is one PE tall, so its y does not appear in these entries. A named dimension stays in the map at any size, even 1.

Per-PE data with miset

The Makefile loads one reference base into each PE:

	miset image=$@ out=$@ map=bench.map \
	    symbol=grid.ref_base order="j n" data=ref.bin
	miset image=$@ out=$@ map=bench.map \
	    symbol=left_moat.query order="n" data=query.bin

order lists the symbol's dimensions, slowest first. Lesson 2 used order="n" for one PE. order="j n" reads ref.bin as one block of n values per PE, PE 0 first: the file's W values go to the W PEs. Had the rectangle no names, its dimension would be x, and the order "x n". display.py reads the results with order="n j": history[0] of every PE, then history[1], which is H row by row.

The colorpair beach

A color takes its input from one port at a time (lesson 4). PE j must receive from the left and also send its own data right, which takes two inputs, so it needs two colors. With two plain colors the rules would alternate from PE to PE. A colorpair lets every PE use one socket and the same two rules:

let r  = color()
let rp = color()

let beach = rect(0 0 W+2 1)
beach : paint(rp [CXR|R])
beach : paint(r  [L|LXC])
beach : colorpair(r rp)

The beach is the rectangle that carries the routes: the W PEs and both moats. colorpair(r rp) makes Paint number the two colors as partners, one apart (r = rp ^ 1). In a rule, X takes the place of > and makes the color take its input from its partner. Term by term:

Term	in rp [CXR|R]	in r [L|LXC]
the color being routed	rp	r
a | before the input	none	L|: on the left edge, no input from the left
the input	C, the core	L, the left neighbour
X in place of >	read the partner: rp takes what the core sends on r	read the partner: r takes what arrives from the left on rp
the outputs	R, the right neighbour	C, the core
a | after the outputs	|R: on the right edge, no output to the right	none

Every PE reads and writes the one socket r, and the wavelets cross between PEs on rp.

In this program the two | terms change nothing: the right moat never sends, and nothing arrives from the left of the left moat. They keep the row's route inside the row when an edge PE does send, or when the row sits next to other PEs.

A | term acts at the edges of each rectangle the painted region is built from, not at the outline of the region. Paint the rules on one rect that covers the whole row, as beach does. See Pitfalls for what happens otherwise.

The union

union(beach interior left_moat right_moat)

union joins regions at their own positions. Here they overlap: beach covers all W + 2 PEs and carries only routes, and the three code regions cover parts of it and carry only code. Paint accepts the overlap because only one of the overlapping regions has code on any PE. It rejects a union in which two regions both put code on one PE, even the same kernel.

What the hardware does

The swap in the router

Each color in a router chooses which color it takes its input from. A color can read its own color, or a color one or two numbers away. The choice is made separately for the north/south inputs, the east/west inputs and the input from the core, and the core's input can read only the color itself or its partner c ^ 1. X sets all three to the partner. The router does the swap while routing, with no work from the core.

       PE j's router                             PE j+1's router
core --r--> rp takes it (CXR) --rp, east--> r takes it (LXC) --r--> core

The architecture notes draw the same pattern.

Space and time, measured

make-log holds the when table for the first 16 PEs: the cycle at which PE j finished row i. The first rows and PEs of the tutorial row:

         PE 0   PE 1   PE 2   PE 3   PE 4   PE 5
i = 0:     42     86    130    174    215    259
i = 1:     84    125    172    216    260    301
i = 2:    126    170    211    258    302    346
i = 3:    168    212    256    297    344    388

Each PE has its own counter, which starts when its program enables it. In the simulator, every PE starts its program, and so its counter, in the same cycle, so the table can compare counters across PEs. On the hardware, PEs do not all start together: there, compare only times read on the same PE. The cells of one anti-diagonal finish within a few cycles of each other: H[0][3], H[1][2], H[2][1] and H[3][0] at cycles 174, 172, 170 and 168. A step takes about 43 cycles: the scalar work of one cell on one PE.

The last cell, H[15][31], finishes at cycle 2024: 47 anti-diagonals of about 43 cycles. The long row also has W + Q - 1 = 47 anti-diagonals, and its last cell finishes at cycle 1964, although it computes 320 cells against 512. The run time grows with W + Q, the number of anti-diagonals, not with the number of cells. At most min(W, Q) PEs work at once: 16 of the 32 in tutorial, all 8 in long.

Run it

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

cd tutorial/main-3-rows/08-wavefront
make TEST=tutorial
cat test-tutorial/make-log

Success is a line ending in PASS tutorial [ W=32 Q=16 SEED=0 ARCH=sdr ]. make-log shows the simulated matrix, the best alignment, the when table and PASS: All 512 DP cells match the CPU reference. A failure prints FAIL, and make-log lists the first mismatched cells. make TEST='*' runs all three rows. After editing a file, run make clean first: make does not rebuild an existing test-tutorial/. make animate TEST=snp writes the figure above to test-snp/sw_animation.gif; it needs Python's matplotlib and Pillow.

Pitfalls

• Painting the row's rules on a joined region. With the two rules painted on hstack(left_moat interior right_moat) instead of beach, Paint gives no message, but every PE gets its route input and output turned off. Nothing moves: every cell of history stays 0, and the test fails. Paint the rules on one rect.
• Sends and receives out of step. Every sender must send exactly two values per step. A third send, or a missing one, would shift every later value into the wrong variable, and nothing checks this.
• A wrong order string. order must use the names in the map. order="x n" for grid.ref_base fails with 'x' is not a name or number, because the dimension is called j. order="n" leaves j out, which asks for one value shared by every PE, so miset stops with ref.bin: too long.

Exercise

Change the gap penalty from -1 to -2: the two - 1:hp gap terms in sw.w and GAP in ref.py and display.py. Before you run it, predict the largest score of the tutorial and snp rows, and which cells of the first query row change. Then run make clean and make TEST='*'.

Next

Lesson 9 — Broadcast and reduce along a row computes on the stream itself: every PE adds its own vector to a running sum as it passes.

References

• Regions and code attachment
• Color pairs, color swap modifiers and boundary operators
• What dimension names are for
• miset and miget
• Dispatch blocks