Lesson 3 — How a PE runs a loop: vector instructions and cycles

Lesson 1's loop compiled to one instruction for the whole array; lesson 2's ran one element at a time. This lesson measures the difference. One PE runs twelve computations over the same inputs and reads a cycle counter between them, and the compiler's listing shows what each one became.

What you will be able to do

• Read the PE's cycle counter from a kernel, and time a piece of code with it.
• Find a kernel's listing, and read a vector instruction and the descriptors of its operands.
• Tell from the source which loops become one instruction and which run element by element.
• Estimate a loop's cost from its type and its operands, and know why multiply-adds of three arrays, conditionals and divisions are slow.
• Write a reference that rounds the way the compiled code does.

The program

loops.w runs on one PE. Its inputs are two HP arrays a_h and b_h of N values between 1 and 4, which ref.py draws at random from a fixed seed, so every run gets the same values. miset writes them into the image. The kernel first copies them to two SP arrays, a_s and b_s. It then runs these loops, each over all N elements. k is the constant 0.3, and each c starts as a copy of b:

Output	Loop
scale_h	k * a in HP
add_h, mul_h	a + b and a * b in HP
axpy_h	c + k * a in HP
fma_h	c + a * b in HP
prod_h, split_h	the same multiply-add as two statements: p = a * b, then c + p
add_s, mul_s, axpy_s	a + b, a * b and c + k * a in SP
sel_h	the minimum, with the ternary a < b ? a : b
if_h	the minimum, with if/else
div_s	a / b in SP

The ternary condition ? x : y gives x when the condition is true and y otherwise. Unlike C, it evaluates both x and y; use if/else when only one of them may run.

Between the loops it stores the cycle counter into t. After the run, ref.py checks every output, and cycles.py prints how many cycles each loop took. The test passes when every output is right; the cycle counts never decide it.

ref.py --a_h.bin, b_h.bin--> miset --> +-----------------+ --> miget --> outputs --> ref.py check
                                       | PE: loops.w     |
                                       | cycle counter   |              t --> cycles.py (table)
                                       +-----------------+

Reading the cycle counter

u16 tsc[1] address(0x7b30);
u16 tsc_control[1] address(0x7b34);
u16 t[n=14];

address(0x7b30) places the array tsc on the PE's timestamp counter instead of in memory, so reading tsc[0] reads the counter. [1] gives tsc one element and no dimension name: a dimension name is optional, and the map names unnamed dimensions n0, n1 and so on. The counter is 48 bits wide, in three 16-bit words starting at 0x7b30; tsc[0] is the low word, which wraps to 0 every 65,536 cycles. tsc_control is the counter's control register, and its bit 2 makes the counter count every cycle:

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

Each timed loop then sits between two reads:

    t[2] <- tsc[0];
    ∀i add_h[i] <- a_h[i] + b_h[i];
    t[3] <- tsc[0];

cycles.py subtracts consecutive values of t modulo 65,536, which is right for any interval shorter than 65,536 cycles. The first two reads have nothing between them. Their difference, 2 cycles in this run, is in every other interval too. It can change when the code around it changes, so read it from each run.

The full kernel: loops.w

xp param N;

hp a_h[n=N];
hp b_h[n=N];
sp a_s[n=N];
sp b_s[n=N];

hp scale_h[n=N];
hp add_h[n=N];
hp mul_h[n=N];
hp axpy_h[n=N];
hp fma_h[n=N];
hp prod_h[n=N];
hp split_h[n=N];
sp add_s[n=N];
sp mul_s[n=N];
sp axpy_s[n=N];
hp sel_h[n=N];
hp if_h[n=N];
sp div_s[n=N];

u16 tsc[1] address(0x7b30);
u16 tsc_control[1] address(0x7b34);
u16 t[n=14];

main {
    let i in [0, N);
    hp k_h = 0.3:hp;
    sp k_s = 0.3:sp;

    tsc_control[0] <- tsc_control[0] | 4;
    ∀i a_s[i] <- (sp)a_h[i];
    ∀i b_s[i] <- (sp)b_h[i];
    ∀i axpy_h[i] <- b_h[i];
    ∀i fma_h[i] <- b_h[i];
    ∀i split_h[i] <- b_h[i];
    ∀i axpy_s[i] <- b_s[i];

    t[0] <- tsc[0];
    t[1] <- tsc[0];
    ∀i scale_h[i] <- k_h * a_h[i];
    t[2] <- tsc[0];
    ∀i add_h[i] <- a_h[i] + b_h[i];
    t[3] <- tsc[0];
    ∀i mul_h[i] <- a_h[i] * b_h[i];
    t[4] <- tsc[0];
    ∀i axpy_h[i] <- axpy_h[i] + k_h * a_h[i];
    t[5] <- tsc[0];
    ∀i fma_h[i] <- fma_h[i] + a_h[i] * b_h[i];
    t[6] <- tsc[0];
    ∀i prod_h[i] <- a_h[i] * b_h[i];
    ∀i split_h[i] <- split_h[i] + prod_h[i];
    t[7] <- tsc[0];
    ∀i add_s[i] <- a_s[i] + b_s[i];
    t[8] <- tsc[0];
    ∀i mul_s[i] <- a_s[i] * b_s[i];
    t[9] <- tsc[0];
    ∀i axpy_s[i] <- axpy_s[i] + k_s * a_s[i];
    t[10] <- tsc[0];
    ∀i sel_h[i] <- a_h[i] < b_h[i] ? a_h[i] : b_h[i];
    t[11] <- tsc[0];
    ∀i {
        if (a_h[i] < b_h[i]) if_h[i] <- a_h[i];
        else if_h[i] <- b_h[i];
    }
    t[12] <- tsc[0];
    ∀i div_s[i] <- a_s[i] / b_s[i];
    t[13] <- tsc[0];
}

The counter is enabled before the copy loops, not just before the first read: a read right after the enabling write can still return 0.

The listing

Carbide writes a listing of the code it generates for each kernel. For the tutorial row it is test-tutorial/loops_artefacts/loops1.yeti: the kernel's name, then the number Paint gives that compile. test-tutorial/bench-listing.json holds the listings of all kernels in one JSON file, for the debug tools.

The listing first declares the descriptors. A data structure descriptor (DSD) describes a vector operand: here, an array in memory by its base address, its length and its stride. These three describe add_h, a_h and b_h for the HP addition:

vector add_h_dst_14 register(T2b) = {.base=(add_h+0), .length=64, .stride=1, .save=0, .vp=0, .ua=0, .rchoasjt=0, };
vector a_h_src0_15 register(T01) = {.base=(a_h+0), .length=64, .stride=1, .save=0, .vp=0, .ua=0, .rchoasjt=0, };
vector b_h_src1_16 register(T1c) = {.base=(b_h+0), .length=64, .stride=1, .save=0, .vp=0, .ua=0, .rchoasjt=0, };

register(T2b) is the descriptor register that holds the descriptor while the program runs. The loop ∀i add_h[i] <- a_h[i] + b_h[i]; is then one instruction, whose operands are those three registers, destination first:

	faddh [T2b] = [T01], [T1c];  //@../loops.w:43

The comment gives the source line. A few instructions carry line 1 instead: the ternary's compare, flth, is marked loops.w:1. The descriptor names and register numbers can differ in your build. The other loops in the listing:

• c + k * a is one multiply-add, fmach [T2d] = [T03], [T1e], R0;. Register R0 holds k.

• c + a * b is a loop that runs once per element. Before the loop, setjt stores the address where the program goes on afterwards, main_bb3. The loop's first line reads T04, a descriptor that only counts: after N uses (64 in the tutorial row) its length is 0, and the CE then jumps to main_bb3, so jmp main_bb1 repeats the body N times, not forever. The next three lines move one element of each array into a register, fmach multiplies and adds on registers, and the last move stores the result:

  	setjt JTR0 = main_bb3;  //@../loops.w:49
  
  main_bb1: 
  	imov16 Ra = [T04];  //@../loops.w:49
  	fmov16 Rb = [T1f];  //@../loops.w:49
  	fmov16 Rc = [T05];  //@../loops.w:49
  	fmov16 Rd = [T10];  //@../loops.w:49
  	fmach Re = Rb, Rc, Rd;  //@../loops.w:49
  	fmov16 [T2e] = Re;  //@../loops.w:49
  	jmp main_bb1;  //@../loops.w:49

  Each of its descriptors moves to the next element every time it is used.

• The two-statement form is two instructions, fmulh and faddh.

• The ternary is a loop with a compare, flth, and a select, select16, per element.

• The if/else is a loop with a compare and a conditional jump per element.

• The division is a loop that runs the library routine __tungsten_lib_divsp, which Carbide copies into the loop body.

The results

Each count includes the 2 cycles of one timer read.

Loop	What Carbide made of it	N=64 cycles	per element	N=62 cycles	per element
two timer reads		2		2
HP c = k*a	one fmulh	12	0.19	11	0.18
HP c = a + b	one faddh	23	0.36	20	0.32
HP c = a * b	one fmulh	22	0.34	20	0.32
HP c = c + k*a	one fmach	21	0.33	20	0.32
HP c = c + a*b	a loop per element	518	8.09	502	8.10
HP p = a*b; then c = c + p	fmulh, faddh	41	0.64	39	0.63
SP c = a + b	one fadds	34	0.53	33	0.53
SP c = a * b	one fmuls	132	2.06	67	1.08
SP c = c + k*a	one fmacs	135	2.11	69	1.11
HP c = a < b ? a : b	a loop per element	393	6.14	257	4.15
HP if (a < b) c = a; else c = b;	a loop per element	1303	20.36	1196	19.29
SP c = a / b	a loop per element	6092	95.19	5902	95.19

The rows of test.rc are tutorial (N=64) and n62 (N=62).

What the hardware does

One instruction for a whole vector

Most CE instructions take DSDs as operands, so one instruction processes a whole vector. The datapath is 128 bits wide: one instruction does up to eight 16-bit or four 32-bit operations per cycle. An instruction also has a fixed cost of a few cycles to start. The HP c = k*a took 12 cycles for 64 elements: 8 cycles of work at eight per cycle, and 4 for the start and the timer read.

Three limits lower the rate:

• SP multiplication. The SP multiplier does one multiply, or one multiply-add, per cycle. HP multiplies run eight per cycle.
• Memory reads. Memory delivers at most four 32-bit words per cycle, one from each of four memory banks. One HP array read at eight elements per cycle uses all four, so an instruction with two HP array operands gets four elements per cycle (c = a + b, c = c + k*a), and one with two SP array operands gets two (SP c = a + b).
• Memory-bank conflicts. Consecutive 32-bit words lie in consecutive memory banks, so two 32-bit words whose addresses differ by a multiple of 8 (16-bit words) lie in the same memory bank. When two operands need the same memory bank in the same cycle, one waits. With N=64, each SP array is 128 words long, so a_s[i] and b_s[i] always share a memory bank, and the SP multiply takes two cycles per element. With N=62 the arrays are 124 words apart, and it takes one. The SP addition reads the same two arrays and is not slowed at N=64. It reads four elements of each array per step, one word from every memory bank, so its two operands take turns at the memory banks wherever the arrays sit: that is the limit of two elements per cycle above. The SP multiply reads one element of each array per cycle, so it waits only when those two elements share a memory bank. The ternary loop also depends on the placement: 6.14 cycles per element at N=64, 4.15 at N=62.

Loops that are not one instruction

A multiply-add of three arrays. This has two causes, one in the hardware and one in the compiler:

• In a multiply-add instruction (fmach, fmacs), the third operand is always a register. A multiply-add of two arrays and a scalar fits: the scalar k goes in the register.
• Carbide fuses c + a * b into one multiply-add on purpose. With three arrays no operand can be the register, so Carbide runs the fused form one element at a time, at 8 cycles per element.

Written as two statements, p = a * b and then c = c + p, the same work is two vector instructions, 13 times faster. The two forms round differently (next section).

Conditionals. The ternary becomes a loop that compares and selects, one element at a time. The if/else becomes a loop of branches, which is slower still. Both compute the same minimum.

Division. The CE has no instruction that computes a quotient. A division calls a library routine that builds the quotient from multiplies and adds, about 95 cycles per element. The routine is not correctly rounded: ref.py finds that 17 of the 64 quotients differ from the correctly rounded quotient in the last bit.

The rule of thumb: a loop whose body is one operation over arrays, with at most two array operands, becomes one vector instruction, and its rate follows from the instruction, the operand types and where the arrays sit, as above. A multiply-add of three arrays, a condition or a division runs element by element, at 4 to 100 cycles per element. The listing shows which case a loop is.

Rounding, and the reference

A multiply-add instruction rounds once: it adds the exact product and rounds the sum. The two-statement form rounds twice: the product to HP, then the sum. ref.py models each form. Its a64 and b64 are the inputs as 64-bit Python floats, which hold every product and sum here exactly, so each hp(...) is one rounding:

prod = hp(a64 * b64)
expect = {
    ...
    "axpy_h": hp(b64 + k_h * a64),
    "fma_h": hp(b64 + a64 * b64),
    "prod_h": prod,
    "split_h": hp(b64 + prod),
    ...
}

Each output matches its own reference bit for bit, and fma_h and split_h differ from each other in 17 of the 64 elements.

Carbide fuses a product and a sum wherever it can (see Floating-point evaluation). A reference that rounds after every operation, as NumPy does in float16, would fail on fma_h. A reference must round the way the compiled code does, or compare within a tolerance. ref.py checks the division within one unit in the last place (ulp).

The test

The Makefile writes the inputs into the image with miset, as in lesson 2, reads every output and t with miget, and runs the two scripts:

test/$(MODE): out-$(MODE).bin
	python3 ../ref.py check $(N)
	python3 ../cycles.py $(N) t.bin

python3 ../ref.py inputs N writes the inputs a_h.bin and b_h.bin. python3 ../ref.py check N compares each output with its reference, prints one line per output, and fails the test if any output differs, or if a quotient is more than one ulp off. cycles.py prints the table and never fails: cycle counts change with the compiler, so the test does not check them.

Run it

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

cd tutorial/main-1-one-pe/03-loops-and-cycles
make TEST='*'

Each row prints a PASS line. The checks and the table are in each row's log:

sed -n '/ref.py check/,$p' test-tutorial/make-log

For the tutorial row the check prints:

scale_h  64 of 64 exact
add_h    64 of 64 exact
mul_h    64 of 64 exact
axpy_h   64 of 64 exact
fma_h    64 of 64 exact
prod_h   64 of 64 exact
split_h  64 of 64 exact
add_s    64 of 64 exact
mul_s    64 of 64 exact
axpy_s   64 of 64 exact
sel_h    64 of 64 exact
if_h     64 of 64 exact
fma_h and split_h differ in 17 of 64 elements
div_s    47 of 64 exact, 17 off by one ulp

and then the cycle table. A failing row prints FAIL, and the log shows which output had fewer exact elements.

Pitfalls

• Time only intervals shorter than 65,536 cycles, or read all three words of the counter, low word first: the low word wraps.
• Enable the counter well before the first read, as loops.w does.
• Check in the listing that the code you time sits between its two reads. Carbide can move instructions across a timer read. Here it moved a descriptor load of line 58, lds1wds T15 = [a_s_src1_36];, into the interval that times lines 51 and 52.
• Split a multiply-add of three arrays into two statements when it is on a hot path, and check what the second rounding does to your results.
• Multiply by a reciprocal instead of dividing where the last bit allows it, and compare quotients within a tolerance.
• Look at the placement when a loop is twice as slow as expected. Two arrays whose elements always share a memory bank halve the rate of an SP multiply.

Exercise

Add two loops at the end of the kernel, each between its own pair of timer reads:

1. q1[i] <- a_s[i] / 3.0:sp;
2. q2[i] <- a_s[i] * r;, where sp r = 1.0:sp / 3.0:sp; is computed once before the loop.

Before you run them, predict their cycle counts from the table, and predict whether q1 and q2 are equal. Extend t, cycles.py, the Makefile's OUTPUTS and ref.py (compare both with the correctly rounded quotient, within one ulp), then run make clean and the test.

Next

Lesson 4 — Sockets, wavelets, colors and routes moves data between PEs.

References

• Hardware register access: address()
• Execution domains and hardware instructions
• Floating-point evaluation
• Ternary operator
• Vector operands (DSDs)
• Compiler usage: output files