1 # bend-ml-tensor-gpu: tensor-array's products on the GPU, with tensor-array's bits. 2 # 3 # import ../tensor-array/main.bend as TA 4 # import bend-ml-tensor-gpu@0.1.2.1/main.bend as TG 5 # TG.matmul(n, k, m, par, a, b) # the same arguments, result type and bits as TA.Mat.matmul 6 # 7 # The six products (matmul, matmul_nt, matmul_tn and their accumulating forms) take tensor-array's Mat values and return 8 # its result types; par is ignored (the GPU decides the split). Each product is one `!` call (a bang): a balanced fork 9 # tree whose leaves each compute a 4 x 4 tile of C, every node sharing A, B and C through Base's Array.fork (two 10 # handles to one block, no copy), as upstream's one-block pattern does (bendlang/bend tests/run/stencil3d.bend). On a 11 # machine without a GPU, or with `--gpu off`, the same bang runs on the CPU threads. 12 # 13 # The bits: a cell of C is computed in tensor-array's order. Product p of a dot goes to partial sum p mod 4, in 14 # increasing p, each sum starting at +0.0, and the four sums are combined by TA.t8_comb with rot = k mod 4, as 15 # tensor-array's dotr and t8_dot do; an accumulating product adds the cell to C once with TA.cwr. The loop takes four 16 # steps of k a turn; the steps past k in the last turn add +0.0, which leaves a sum unchanged (a sum that starts at +0.0 17 # is never -0.0). Bend builds the device code with --fmad=false, so no multiply-add is fused. tests/bits.bend checks the 18 # bits against tensor-array on many shapes, on the GPU and the CPU. 19 # 20 # Not proved, and marked so: Array.fork and Array.join are @unsafe in Base, so the defs below that call them, or take 21 # a shared array, are @unsafe, and `bend main.bend` names every def that relies on them. Nothing checks that a leaf 22 # writes only its own tile of C, or never writes A or B; the tiles are disjoint by construction (leaf t writes tile t). 23 # All of it is Bend: no foreign code. Written by bend-ml's scripts/gen_tensor_gpu.py (4 x 4 tiles). 24 import Base 25 import bend-ml-tensor-array@0.1.14.0/main.bend as TA 26 27 # the sizes and strides of one product: C[i][j] = sum_p A[oa + i*arow + p*sa] * B[ob + p*sb + j*bcol], n x m, k steps 28 # (oa, ob: 0 except for matmul_precise's blocks) 29 type E is Data: 30 E{n: U32, m: U32, k: U32, tcols: U32, tiles: U32, arow: U32, sa: U32, sb: U32, bcol: U32, kg: Nat, rot: Nat, rem: U32, ac: Bool, oa: U32, ob: U32} 31 32 # a product that counts only when its step is inside k (the last turn's steps past k add +0.0) 33 def pk(ok: Bool, +x: F32) -> F32: 34 match ok: 35 case True{}: x 36 case False{}: 0.0 37 38 # the values of one turn: A's 4 rows and B's 4 columns at four consecutive steps, read through the arrays in order 39 type GR is Type: 40 GR{a: Array, b: Array, x0_0: F32, x0_1: F32, x0_2: F32, x0_3: F32, x1_0: F32, x1_1: F32, x1_2: F32, x1_3: F32, x2_0: F32, x2_1: F32, x2_2: F32, x2_3: F32, x3_0: F32, x3_1: F32, x3_2: F32, x3_3: F32, y0_0: F32, y0_1: F32, y0_2: F32, y0_3: F32, y1_0: F32, y1_1: F32, y1_2: F32, y1_3: F32, y2_0: F32, y2_1: F32, y2_2: F32, y2_3: F32, y3_0: F32, y3_1: F32, y3_2: F32, y3_3: F32} 41 42 def gr_31(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, +y0_0: F32, +y0_1: F32, +y0_2: F32, +y0_3: F32, +y1_0: F32, +y1_1: F32, +y1_2: F32, +y1_3: F32, +y2_0: F32, +y2_1: F32, +y2_2: F32, +y2_3: F32, +y3_0: F32, +y3_1: F32, +y3_2: F32, a: Array, rv: Array & F32) -> GR: 43 match rv: 44 case Tuple{b, +y3_3}: 45 GR{a, b, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, y0_1, y0_2, y0_3, y1_0, y1_1, y1_2, y1_3, y2_0, y2_1, y2_2, y2_3, y3_0, y3_1, y3_2, y3_3} 46 def gr_30(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, +y0_0: F32, +y0_1: F32, +y0_2: F32, +y0_3: F32, +y1_0: F32, +y1_1: F32, +y1_2: F32, +y1_3: F32, +y2_0: F32, +y2_1: F32, +y2_2: F32, +y2_3: F32, +y3_0: F32, +y3_1: F32, a: Array, rv: Array & F32) -> GR: 47 match rv: 48 case Tuple{b, +y3_2}: 49 gr_31(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, y0_1, y0_2, y0_3, y1_0, y1_1, y1_2, y1_3, y2_0, y2_1, y2_2, y2_3, y3_0, y3_1, y3_2, a, Array.get(F32, b, (ib + ((sb * 3 : U32) + (bcol * 3 : U32) : U32) : U32))) 50 def gr_29(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, +y0_0: F32, +y0_1: F32, +y0_2: F32, +y0_3: F32, +y1_0: F32, +y1_1: F32, +y1_2: F32, +y1_3: F32, +y2_0: F32, +y2_1: F32, +y2_2: F32, +y2_3: F32, +y3_0: F32, a: Array, rv: Array & F32) -> GR: 51 match rv: 52 case Tuple{b, +y3_1}: 53 gr_30(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, y0_1, y0_2, y0_3, y1_0, y1_1, y1_2, y1_3, y2_0, y2_1, y2_2, y2_3, y3_0, y3_1, a, Array.get(F32, b, (ib + ((sb * 3 : U32) + (bcol * 2 : U32) : U32) : U32))) 54 def gr_28(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, +y0_0: F32, +y0_1: F32, +y0_2: F32, +y0_3: F32, +y1_0: F32, +y1_1: F32, +y1_2: F32, +y1_3: F32, +y2_0: F32, +y2_1: F32, +y2_2: F32, +y2_3: F32, a: Array, rv: Array & F32) -> GR: 55 match rv: 56 case Tuple{b, +y3_0}: 57 gr_29(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, y0_1, y0_2, y0_3, y1_0, y1_1, y1_2, y1_3, y2_0, y2_1, y2_2, y2_3, y3_0, a, Array.get(F32, b, (ib + ((sb * 3 : U32) + (bcol * 1 : U32) : U32) : U32))) 58 def gr_27(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, +y0_0: F32, +y0_1: F32, +y0_2: F32, +y0_3: F32, +y1_0: F32, +y1_1: F32, +y1_2: F32, +y1_3: F32, +y2_0: F32, +y2_1: F32, +y2_2: F32, a: Array, rv: Array & F32) -> GR: 59 match rv: 60 case Tuple{b, +y2_3}: 61 gr_28(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, y0_1, y0_2, y0_3, y1_0, y1_1, y1_2, y1_3, y2_0, y2_1, y2_2, y2_3, a, Array.get(F32, b, (ib + ((sb * 3 : U32) + (bcol * 0 : U32) : U32) : U32))) 62 def gr_26(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, +y0_0: F32, +y0_1: F32, +y0_2: F32, +y0_3: F32, +y1_0: F32, +y1_1: F32, +y1_2: F32, +y1_3: F32, +y2_0: F32, +y2_1: F32, a: Array, rv: Array & F32) -> GR: 63 match rv: 64 case Tuple{b, +y2_2}: 65 gr_27(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, y0_1, y0_2, y0_3, y1_0, y1_1, y1_2, y1_3, y2_0, y2_1, y2_2, a, Array.get(F32, b, (ib + ((sb * 2 : U32) + (bcol * 3 : U32) : U32) : U32))) 66 def gr_25(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, +y0_0: F32, +y0_1: F32, +y0_2: F32, +y0_3: F32, +y1_0: F32, +y1_1: F32, +y1_2: F32, +y1_3: F32, +y2_0: F32, a: Array, rv: Array & F32) -> GR: 67 match rv: 68 case Tuple{b, +y2_1}: 69 gr_26(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, y0_1, y0_2, y0_3, y1_0, y1_1, y1_2, y1_3, y2_0, y2_1, a, Array.get(F32, b, (ib + ((sb * 2 : U32) + (bcol * 2 : U32) : U32) : U32))) 70 def gr_24(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, +y0_0: F32, +y0_1: F32, +y0_2: F32, +y0_3: F32, +y1_0: F32, +y1_1: F32, +y1_2: F32, +y1_3: F32, a: Array, rv: Array & F32) -> GR: 71 match rv: 72 case Tuple{b, +y2_0}: 73 gr_25(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, y0_1, y0_2, y0_3, y1_0, y1_1, y1_2, y1_3, y2_0, a, Array.get(F32, b, (ib + ((sb * 2 : U32) + (bcol * 1 : U32) : U32) : U32))) 74 def gr_23(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, +y0_0: F32, +y0_1: F32, +y0_2: F32, +y0_3: F32, +y1_0: F32, +y1_1: F32, +y1_2: F32, a: Array, rv: Array & F32) -> GR: 75 match rv: 76 case Tuple{b, +y1_3}: 77 gr_24(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, y0_1, y0_2, y0_3, y1_0, y1_1, y1_2, y1_3, a, Array.get(F32, b, (ib + ((sb * 2 : U32) + (bcol * 0 : U32) : U32) : U32))) 78 def gr_22(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, +y0_0: F32, +y0_1: F32, +y0_2: F32, +y0_3: F32, +y1_0: F32, +y1_1: F32, a: Array, rv: Array & F32) -> GR: 79 match rv: 80 case Tuple{b, +y1_2}: 81 gr_23(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, y0_1, y0_2, y0_3, y1_0, y1_1, y1_2, a, Array.get(F32, b, (ib + ((sb * 1 : U32) + (bcol * 3 : U32) : U32) : U32))) 82 def gr_21(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, +y0_0: F32, +y0_1: F32, +y0_2: F32, +y0_3: F32, +y1_0: F32, a: Array, rv: Array & F32) -> GR: 83 match rv: 84 case Tuple{b, +y1_1}: 85 gr_22(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, y0_1, y0_2, y0_3, y1_0, y1_1, a, Array.get(F32, b, (ib + ((sb * 1 : U32) + (bcol * 2 : U32) : U32) : U32))) 86 def gr_20(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, +y0_0: F32, +y0_1: F32, +y0_2: F32, +y0_3: F32, a: Array, rv: Array & F32) -> GR: 87 match rv: 88 case Tuple{b, +y1_0}: 89 gr_21(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, y0_1, y0_2, y0_3, y1_0, a, Array.get(F32, b, (ib + ((sb * 1 : U32) + (bcol * 1 : U32) : U32) : U32))) 90 def gr_19(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, +y0_0: F32, +y0_1: F32, +y0_2: F32, a: Array, rv: Array & F32) -> GR: 91 match rv: 92 case Tuple{b, +y0_3}: 93 gr_20(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, y0_1, y0_2, y0_3, a, Array.get(F32, b, (ib + ((sb * 1 : U32) + (bcol * 0 : U32) : U32) : U32))) 94 def gr_18(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, +y0_0: F32, +y0_1: F32, a: Array, rv: Array & F32) -> GR: 95 match rv: 96 case Tuple{b, +y0_2}: 97 gr_19(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, y0_1, y0_2, a, Array.get(F32, b, (ib + ((sb * 0 : U32) + (bcol * 3 : U32) : U32) : U32))) 98 def gr_17(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, +y0_0: F32, a: Array, rv: Array & F32) -> GR: 99 match rv: 100 case Tuple{b, +y0_1}: 101 gr_18(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, y0_1, a, Array.get(F32, b, (ib + ((sb * 0 : U32) + (bcol * 2 : U32) : U32) : U32))) 102 def gr_16(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, +x3_3: F32, a: Array, rv: Array & F32) -> GR: 103 match rv: 104 case Tuple{b, +y0_0}: 105 gr_17(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, y0_0, a, Array.get(F32, b, (ib + ((sb * 0 : U32) + (bcol * 1 : U32) : U32) : U32))) 106 def gr_15(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, +x3_2: F32, rv: Array & F32, b: Array) -> GR: 107 match rv: 108 case Tuple{a, +x3_3}: 109 gr_16(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, x3_3, a, Array.get(F32, b, (ib + ((sb * 0 : U32) + (bcol * 0 : U32) : U32) : U32))) 110 def gr_14(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, +x3_1: F32, rv: Array & F32, b: Array) -> GR: 111 match rv: 112 case Tuple{a, +x3_2}: 113 gr_15(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, x3_2, Array.get(F32, a, (ia + ((sa * 3 : U32) + (arow * 3 : U32) : U32) : U32)), b) 114 def gr_13(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, +x3_0: F32, rv: Array & F32, b: Array) -> GR: 115 match rv: 116 case Tuple{a, +x3_1}: 117 gr_14(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, x3_1, Array.get(F32, a, (ia + ((sa * 3 : U32) + (arow * 2 : U32) : U32) : U32)), b) 118 def gr_12(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, +x2_3: F32, rv: Array & F32, b: Array) -> GR: 119 match rv: 120 case Tuple{a, +x3_0}: 121 gr_13(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, x3_0, Array.get(F32, a, (ia + ((sa * 3 : U32) + (arow * 1 : U32) : U32) : U32)), b) 122 def gr_11(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, +x2_2: F32, rv: Array & F32, b: Array) -> GR: 123 match rv: 124 case Tuple{a, +x2_3}: 125 gr_12(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, x2_3, Array.get(F32, a, (ia + ((sa * 3 : U32) + (arow * 0 : U32) : U32) : U32)), b) 126 def gr_10(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, +x2_1: F32, rv: Array & F32, b: Array) -> GR: 127 match rv: 128 case Tuple{a, +x2_2}: 129 gr_11(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, x2_2, Array.get(F32, a, (ia + ((sa * 2 : U32) + (arow * 3 : U32) : U32) : U32)), b) 130 def gr_9(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, +x2_0: F32, rv: Array & F32, b: Array) -> GR: 131 match rv: 132 case Tuple{a, +x2_1}: 133 gr_10(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, x2_1, Array.get(F32, a, (ia + ((sa * 2 : U32) + (arow * 2 : U32) : U32) : U32)), b) 134 def gr_8(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, +x1_3: F32, rv: Array & F32, b: Array) -> GR: 135 match rv: 136 case Tuple{a, +x2_0}: 137 gr_9(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, x2_0, Array.get(F32, a, (ia + ((sa * 2 : U32) + (arow * 1 : U32) : U32) : U32)), b) 138 def gr_7(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, +x1_2: F32, rv: Array & F32, b: Array) -> GR: 139 match rv: 140 case Tuple{a, +x1_3}: 141 gr_8(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, x1_3, Array.get(F32, a, (ia + ((sa * 2 : U32) + (arow * 0 : U32) : U32) : U32)), b) 142 def gr_6(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, +x1_1: F32, rv: Array & F32, b: Array) -> GR: 143 match rv: 144 case Tuple{a, +x1_2}: 145 gr_7(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, x1_2, Array.get(F32, a, (ia + ((sa * 1 : U32) + (arow * 3 : U32) : U32) : U32)), b) 146 def gr_5(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, +x1_0: F32, rv: Array & F32, b: Array) -> GR: 147 match rv: 148 case Tuple{a, +x1_1}: 149 gr_6(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, x1_1, Array.get(F32, a, (ia + ((sa * 1 : U32) + (arow * 2 : U32) : U32) : U32)), b) 150 def gr_4(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, +x0_3: F32, rv: Array & F32, b: Array) -> GR: 151 match rv: 152 case Tuple{a, +x1_0}: 153 gr_5(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, x1_0, Array.get(F32, a, (ia + ((sa * 1 : U32) + (arow * 1 : U32) : U32) : U32)), b) 154 def gr_3(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, +x0_2: F32, rv: Array & F32, b: Array) -> GR: 155 match rv: 156 case Tuple{a, +x0_3}: 157 gr_4(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, x0_3, Array.get(F32, a, (ia + ((sa * 1 : U32) + (arow * 0 : U32) : U32) : U32)), b) 158 def gr_2(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, +x0_1: F32, rv: Array & F32, b: Array) -> GR: 159 match rv: 160 case Tuple{a, +x0_2}: 161 gr_3(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, x0_2, Array.get(F32, a, (ia + ((sa * 0 : U32) + (arow * 3 : U32) : U32) : U32)), b) 162 def gr_1(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +x0_0: F32, rv: Array & F32, b: Array) -> GR: 163 match rv: 164 case Tuple{a, +x0_1}: 165 gr_2(ia, ib, arow, sa, sb, bcol, x0_0, x0_1, Array.get(F32, a, (ia + ((sa * 0 : U32) + (arow * 2 : U32) : U32) : U32)), b) 166 def gr_0(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, rv: Array & F32, b: Array) -> GR: 167 match rv: 168 case Tuple{a, +x0_0}: 169 gr_1(ia, ib, arow, sa, sb, bcol, x0_0, Array.get(F32, a, (ia + ((sa * 0 : U32) + (arow * 1 : U32) : U32) : U32)), b) 170 # the turn whose A values start at ia and B values at ib 171 def turn(+ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, a: Array, b: Array) -> GR: 172 gr_0(ia, ib, arow, sa, sb, bcol, Array.get(F32, a, ia), b) 173 174 # the tile's cells after the four sums are combined, row by row 175 type T is Type: 176 T{a: Array, b: Array, c0_0: F32, c0_1: F32, c0_2: F32, c0_3: F32, c1_0: F32, c1_1: F32, c1_2: F32, c1_3: F32, c2_0: F32, c2_1: F32, c2_2: F32, c2_3: F32, c3_0: F32, c3_1: F32, c3_2: F32, c3_3: F32} 177 178 # kg more full turns; g holds the current turn's values. The last call (kg = 0) adds the turn that holds the steps 179 # past the last full turn (rem = k mod 4 of them), counting only those, and combines the sums. 180 @unsafe 181 def mm(kg: Nat, +ia: U32, +ib: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +rem: U32, +rot: Nat, +s0_0_0: F32, +s1_0_0: F32, +s2_0_0: F32, +s3_0_0: F32, +s0_0_1: F32, +s1_0_1: F32, +s2_0_1: F32, +s3_0_1: F32, +s0_0_2: F32, +s1_0_2: F32, +s2_0_2: F32, +s3_0_2: F32, +s0_0_3: F32, +s1_0_3: F32, +s2_0_3: F32, +s3_0_3: F32, +s0_1_0: F32, +s1_1_0: F32, +s2_1_0: F32, +s3_1_0: F32, +s0_1_1: F32, +s1_1_1: F32, +s2_1_1: F32, +s3_1_1: F32, +s0_1_2: F32, +s1_1_2: F32, +s2_1_2: F32, +s3_1_2: F32, +s0_1_3: F32, +s1_1_3: F32, +s2_1_3: F32, +s3_1_3: F32, +s0_2_0: F32, +s1_2_0: F32, +s2_2_0: F32, +s3_2_0: F32, +s0_2_1: F32, +s1_2_1: F32, +s2_2_1: F32, +s3_2_1: F32, +s0_2_2: F32, +s1_2_2: F32, +s2_2_2: F32, +s3_2_2: F32, +s0_2_3: F32, +s1_2_3: F32, +s2_2_3: F32, +s3_2_3: F32, +s0_3_0: F32, +s1_3_0: F32, +s2_3_0: F32, +s3_3_0: F32, +s0_3_1: F32, +s1_3_1: F32, +s2_3_1: F32, +s3_3_1: F32, +s0_3_2: F32, +s1_3_2: F32, +s2_3_2: F32, +s3_3_2: F32, +s0_3_3: F32, +s1_3_3: F32, +s2_3_3: F32, +s3_3_3: F32, g: GR) -> T: 182 match kg g: 183 case 0n GR{a, b, +x0_0, +x0_1, +x0_2, +x0_3, +x1_0, +x1_1, +x1_2, +x1_3, +x2_0, +x2_1, +x2_2, +x2_3, +x3_0, +x3_1, +x3_2, +x3_3, +y0_0, +y0_1, +y0_2, +y0_3, +y1_0, +y1_1, +y1_2, +y1_3, +y2_0, +y2_1, +y2_2, +y2_3, +y3_0, +y3_1, +y3_2, +y3_3}: 184 T{a, b, TA.t8_comb(rot, (s0_0_0 + pk(U32.is_lt(0, rem), (x0_0 * y0_0 : F32)) : F32), (s1_0_0 + pk(U32.is_lt(1, rem), (x1_0 * y1_0 : F32)) : F32), (s2_0_0 + pk(U32.is_lt(2, rem), (x2_0 * y2_0 : F32)) : F32), (s3_0_0 + pk(U32.is_lt(3, rem), (x3_0 * y3_0 : F32)) : F32)), TA.t8_comb(rot, (s0_0_1 + pk(U32.is_lt(0, rem), (x0_0 * y0_1 : F32)) : F32), (s1_0_1 + pk(U32.is_lt(1, rem), (x1_0 * y1_1 : F32)) : F32), (s2_0_1 + pk(U32.is_lt(2, rem), (x2_0 * y2_1 : F32)) : F32), (s3_0_1 + pk(U32.is_lt(3, rem), (x3_0 * y3_1 : F32)) : F32)), TA.t8_comb(rot, (s0_0_2 + pk(U32.is_lt(0, rem), (x0_0 * y0_2 : F32)) : F32), (s1_0_2 + pk(U32.is_lt(1, rem), (x1_0 * y1_2 : F32)) : F32), (s2_0_2 + pk(U32.is_lt(2, rem), (x2_0 * y2_2 : F32)) : F32), (s3_0_2 + pk(U32.is_lt(3, rem), (x3_0 * y3_2 : F32)) : F32)), TA.t8_comb(rot, (s0_0_3 + pk(U32.is_lt(0, rem), (x0_0 * y0_3 : F32)) : F32), (s1_0_3 + pk(U32.is_lt(1, rem), (x1_0 * y1_3 : F32)) : F32), (s2_0_3 + pk(U32.is_lt(2, rem), (x2_0 * y2_3 : F32)) : F32), (s3_0_3 + pk(U32.is_lt(3, rem), (x3_0 * y3_3 : F32)) : F32)), TA.t8_comb(rot, (s0_1_0 + pk(U32.is_lt(0, rem), (x0_1 * y0_0 : F32)) : F32), (s1_1_0 + pk(U32.is_lt(1, rem), (x1_1 * y1_0 : F32)) : F32), (s2_1_0 + pk(U32.is_lt(2, rem), (x2_1 * y2_0 : F32)) : F32), (s3_1_0 + pk(U32.is_lt(3, rem), (x3_1 * y3_0 : F32)) : F32)), TA.t8_comb(rot, (s0_1_1 + pk(U32.is_lt(0, rem), (x0_1 * y0_1 : F32)) : F32), (s1_1_1 + pk(U32.is_lt(1, rem), (x1_1 * y1_1 : F32)) : F32), (s2_1_1 + pk(U32.is_lt(2, rem), (x2_1 * y2_1 : F32)) : F32), (s3_1_1 + pk(U32.is_lt(3, rem), (x3_1 * y3_1 : F32)) : F32)), TA.t8_comb(rot, (s0_1_2 + pk(U32.is_lt(0, rem), (x0_1 * y0_2 : F32)) : F32), (s1_1_2 + pk(U32.is_lt(1, rem), (x1_1 * y1_2 : F32)) : F32), (s2_1_2 + pk(U32.is_lt(2, rem), (x2_1 * y2_2 : F32)) : F32), (s3_1_2 + pk(U32.is_lt(3, rem), (x3_1 * y3_2 : F32)) : F32)), TA.t8_comb(rot, (s0_1_3 + pk(U32.is_lt(0, rem), (x0_1 * y0_3 : F32)) : F32), (s1_1_3 + pk(U32.is_lt(1, rem), (x1_1 * y1_3 : F32)) : F32), (s2_1_3 + pk(U32.is_l 185 case 1n+q GR{a, b, +x0_0, +x0_1, +x0_2, +x0_3, +x1_0, +x1_1, +x1_2, +x1_3, +x2_0, +x2_1, +x2_2, +x2_3, +x3_0, +x3_1, +x3_2, +x3_3, +y0_0, +y0_1, +y0_2, +y0_3, +y1_0, +y1_1, +y1_2, +y1_3, +y2_0, +y2_1, +y2_2, +y2_3, +y3_0, +y3_1, +y3_2, +y3_3}: 186 mm(q, (ia + (sa * 4 : U32) : U32), (ib + (sb * 4 : U32) : U32), arow, sa, sb, bcol, rem, rot, (s0_0_0 + (x0_0 * y0_0 : F32) : F32), (s1_0_0 + (x1_0 * y1_0 : F32) : F32), (s2_0_0 + (x2_0 * y2_0 : F32) : F32), (s3_0_0 + (x3_0 * y3_0 : F32) : F32), (s0_0_1 + (x0_0 * y0_1 : F32) : F32), (s1_0_1 + (x1_0 * y1_1 : F32) : F32), (s2_0_1 + (x2_0 * y2_1 : F32) : F32), (s3_0_1 + (x3_0 * y3_1 : F32) : F32), (s0_0_2 + (x0_0 * y0_2 : F32) : F32), (s1_0_2 + (x1_0 * y1_2 : F32) : F32), (s2_0_2 + (x2_0 * y2_2 : F32) : F32), (s3_0_2 + (x3_0 * y3_2 : F32) : F32), (s0_0_3 + (x0_0 * y0_3 : F32) : F32), (s1_0_3 + (x1_0 * y1_3 : F32) : F32), (s2_0_3 + (x2_0 * y2_3 : F32) : F32), (s3_0_3 + (x3_0 * y3_3 : F32) : F32), (s0_1_0 + (x0_1 * y0_0 : F32) : F32), (s1_1_0 + (x1_1 * y1_0 : F32) : F32), (s2_1_0 + (x2_1 * y2_0 : F32) : F32), (s3_1_0 + (x3_1 * y3_0 : F32) : F32), (s0_1_1 + (x0_1 * y0_1 : F32) : F32), (s1_1_1 + (x1_1 * y1_1 : F32) : F32), (s2_1_1 + (x2_1 * y2_1 : F32) : F32), (s3_1_1 + (x3_1 * y3_1 : F32) : F32), (s0_1_2 + (x0_1 * y0_2 : F32) : F32), (s1_1_2 + (x1_1 * y1_2 : F32) : F32), (s2_1_2 + (x2_1 * y2_2 : F32) : F32), (s3_1_2 + (x3_1 * y3_2 : F32) : F32), (s0_1_3 + (x0_1 * y0_3 : F32) : F32), (s1_1_3 + (x1_1 * y1_3 : F32) : F32), (s2_1_3 + (x2_1 * y2_3 : F32) : F32), (s3_1_3 + (x3_1 * y3_3 : F32) : F32), (s0_2_0 + (x0_2 * y0_0 : F32) : F32), (s1_2_0 + (x1_2 * y1_0 : F32) : F32), (s2_2_0 + (x2_2 * y2_0 : F32) : F32), (s3_2_0 + (x3_2 * y3_0 : F32) : F32), (s0_2_1 + (x0_2 * y0_1 : F32) : F32), (s1_2_1 + (x1_2 * y1_1 : F32) : F32), (s2_2_1 + (x2_2 * y2_1 : F32) : F32), (s3_2_1 + (x3_2 * y3_1 : F32) : F32), (s0_2_2 + (x0_2 * y0_2 : F32) : F32), (s1_2_2 + (x1_2 * y1_2 : F32) : F32), (s2_2_2 + (x2_2 * y2_2 : F32) : F32), (s3_2_2 + (x3_2 * y3_2 : F32) : F32), (s0_2_3 + (x0_2 * y0_3 : F32) : F32), (s1_2_3 + (x1_2 * y1_3 : F32) : F32), (s2_2_3 + (x2_2 * y2_3 : F32) : F32), (s3_2_3 + (x3_2 * y3_3 : F32) : F32), (s0_3_0 + (x0_3 * y0_0 : F32) : F32), (s1_3_0 + (x1_3 * y1_0 : F32) : F32), 187 188 # writes cell (i, j) of C when it is inside n x m 189 def wr(ok: Bool, +ac: Bool, +idx: U32, +x: F32, c: Array) -> Array: 190 match ok: 191 case True{}: TA.cwr(ac, idx, x, c) 192 case False{}: c 193 194 # the tile at rows i0.., columns j0.. of C, from the loop's result; A and B are dropped (the caller holds them) 195 @unsafe 196 def put(+ac: Bool, +n: U32, +m: U32, +i0: U32, +j0: U32, c: Array, t: T) -> Array: 197 match t: 198 case T{a, b, +c0_0, +c0_1, +c0_2, +c0_3, +c1_0, +c1_1, +c1_2, +c1_3, +c2_0, +c2_1, +c2_2, +c2_3, +c3_0, +c3_1, +c3_2, +c3_3}: 199 wr(Bool.and(U32.is_lt((i0 + 3 : U32), n), U32.is_lt((j0 + 3 : U32), m)), ac, (((i0 + 3 : U32) * m : U32) + (j0 + 3 : U32) : U32), c3_3, wr(Bool.and(U32.is_lt((i0 + 3 : U32), n), U32.is_lt((j0 + 2 : U32), m)), ac, (((i0 + 3 : U32) * m : U32) + (j0 + 2 : U32) : U32), c3_2, wr(Bool.and(U32.is_lt((i0 + 3 : U32), n), U32.is_lt((j0 + 1 : U32), m)), ac, (((i0 + 3 : U32) * m : U32) + (j0 + 1 : U32) : U32), c3_1, wr(Bool.and(U32.is_lt((i0 + 3 : U32), n), U32.is_lt((j0 + 0 : U32), m)), ac, (((i0 + 3 : U32) * m : U32) + (j0 + 0 : U32) : U32), c3_0, wr(Bool.and(U32.is_lt((i0 + 2 : U32), n), U32.is_lt((j0 + 3 : U32), m)), ac, (((i0 + 2 : U32) * m : U32) + (j0 + 3 : U32) : U32), c2_3, wr(Bool.and(U32.is_lt((i0 + 2 : U32), n), U32.is_lt((j0 + 2 : U32), m)), ac, (((i0 + 2 : U32) * m : U32) + (j0 + 2 : U32) : U32), c2_2, wr(Bool.and(U32.is_lt((i0 + 2 : U32), n), U32.is_lt((j0 + 1 : U32), m)), ac, (((i0 + 2 : U32) * m : U32) + (j0 + 1 : U32) : U32), c2_1, wr(Bool.and(U32.is_lt((i0 + 2 : U32), n), U32.is_lt((j0 + 0 : U32), m)), ac, (((i0 + 2 : U32) * m : U32) + (j0 + 0 : U32) : U32), c2_0, wr(Bool.and(U32.is_lt((i0 + 1 : U32), n), U32.is_lt((j0 + 3 : U32), m)), ac, (((i0 + 1 : U32) * m : U32) + (j0 + 3 : U32) : U32), c1_3, wr(Bool.and(U32.is_lt((i0 + 1 : U32), n), U32.is_lt((j0 + 2 : U32), m)), ac, (((i0 + 1 : U32) * m : U32) + (j0 + 2 : U32) : U32), c1_2, wr(Bool.and(U32.is_lt((i0 + 1 : U32), n), U32.is_lt((j0 + 1 : U32), m)), ac, (((i0 + 1 : U32) * m : U32) + (j0 + 1 : U32) : U32), c1_1, wr(Bool.and(U32.is_lt((i0 + 1 : U32), n), U32.is_lt((j0 + 0 : U32), m)), ac, (((i0 + 1 : U32) * m : U32) + (j0 + 0 : U32) : U32), c1_0, wr(Bool.and(U32.is_lt((i0 + 0 : U32), n), U32.is_lt((j0 + 3 : U32), m)), ac, (((i0 + 0 : U32) * m : U32) + (j0 + 3 : U32) : U32), c0_3, wr(Bool.and(U32.is_lt((i0 + 0 : U32), n), U32.is_lt((j0 + 2 : U32), m)), ac, (((i0 + 0 : U32) * m : U32) + (j0 + 2 : U32) : U32), c0_2, wr(Bool.and(U32.is_lt((i0 + 0 : U32), n), U32.is_lt((j0 + 1 : U32), m)), ac, (((i0 + 0 : U 200 201 @unsafe 202 def tile(+e: E, +i0: U32, +j0: U32, +a: Array, +b: Array, c: Array) -> Array: 203 match e: 204 case E{+n, +m, +k, +tcols, +tiles, +arow, +sa, +sb, +bcol, +kg, +rot, +rem, +ac, +oa, +ob}: 205 put(ac, n, m, i0, j0, c, mm(kg, (oa + (i0 * arow : U32) : U32), (ob + (j0 * bcol : U32) : U32), arow, sa, sb, bcol, rem, rot, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, turn((oa + (i0 * arow : U32) : U32), (ob + (j0 * bcol : U32) : U32), arow, sa, sb, bcol, a, b))) 206 207 # leaf t computes tile t (row-major over the tiles); leaves past the last tile write nothing 208 @unsafe 209 def leaf_if(ok: Bool, +e: E, +t: U32, +tcols: U32, +a: Array, +b: Array, c: Array) -> Array: 210 match ok: 211 case True{}: tile(e, (U32.div(t, tcols) * 4 : U32), (U32.mod(t, tcols) * 4 : U32), a, b, c) 212 case False{}: c 213 214 @unsafe 215 def leaf(+e: E, +t: U32, +a: Array, +b: Array, c: Array) -> Array: 216 match e: 217 case E{+n, +m, +k, +tcols, +tiles, +arow, +sa, +sb, +bcol, +kg, +rot, +rem, +ac, +oa, +ob}: 218 leaf_if(U32.is_lt(t, tiles), e, t, tcols, a, b, c) 219 220 @unsafe 221 def tree(+e: E, d: Nat, +t: U32, +a: Array, +b: Array, +c: Array) -> Array: 222 match d: 223 case 0n: leaf(e, t, a, b, c) 224 case 1n++q: 225 l r = tree(e, q, (t * 2 : U32), a, b, c) tree(e, q, ((t * 2 : U32) + 1 : U32), a, b, c) 226 Array.join(F32, l, r) 227 228 # C made inside the bang: the root of the ! call allocates it on the device. Made on the host instead, a new C often 229 # reuses a block whose pages sit on the device, so zeroing it on the host pulls them back over PCIe and the bang moves 230 # them again (a 256 x 1536 product: 9.2 ms against 6.0, measured 2026-10-18). The leaves write every cell of C. 231 @unsafe 232 def tree_new(+e: E, d: Nat, +a: Array, +b: Array, +cd: Nat) -> Array: 233 tree(e, d, 0, a, b, Array.new(F32, cd, 0.0)) 234 235 @unsafe 236 def bang_new(+e: E, d: Nat, +a: Array, +b: Array, +cd: Nat) -> Array: 237 tree_new!(e, d, a, b, cd) 238 239 # one bang over the tiles; A and B are shared with every leaf and come back unchanged 240 @unsafe 241 def bang(+e: E, d: Nat, +a: Array, +b: Array, c: Array) -> Array: 242 tree!(e, d, 0, a, b, c) 243 244 type ABC is Type: 245 ABC{a: Array, b: Array, c: Array} 246 247 @unsafe 248 def run_e(+e: E, +tiles: U32, +a: Array, +b: Array, c: Array) -> ABC: 249 ABC{a, b, bang(e, TA.cap_depth(U32.to_nat(tiles)), a, b, c)} 250 251 # the product's sizes and strides 252 def e_ofo(+ac: Bool, +n: Nat, +m: Nat, +k: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +oa: U32, +ob: U32) -> E: 253 E{U32.from_nat(n), U32.from_nat(m), U32.from_nat(k), U32.from_nat(Nat.div(Nat.add(m, 3n), 4n)), U32.from_nat(Nat.mul(Nat.div(Nat.add(n, 3n), 4n), Nat.div(Nat.add(m, 3n), 4n))), arow, sa, sb, bcol, Nat.div(k, 4n), Nat.mod(k, 4n), U32.from_nat(Nat.mod(k, 4n)), ac, oa, ob} 254 255 def e_of(+ac: Bool, +n: Nat, +m: Nat, +k: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32) -> E: 256 e_ofo(ac, n, m, k, arow, sa, sb, bcol, 0, 0) 257 258 def tiles(+n: Nat, +m: Nat) -> Nat: 259 Nat.mul(Nat.div(Nat.add(n, 3n), 4n), Nat.div(Nat.add(m, 3n), 4n)) 260 261 # C (n x m) = A · B with k steps into a new C made on the device 262 @unsafe 263 def gemm_new(+n: Nat, +m: Nat, +k: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +a: Array, +b: Array) -> ABC: 264 ABC{a, b, bang_new(e_of(False{}, n, m, k, arow, sa, sb, bcol), TA.cap_depth(tiles(n, m)), a, b, TA.cap_depth(Nat.mul(n, m)))} 265 266 # C (n x m) = A · B with k steps, or C + A · B when ac, C[i][j] = sum_p A[i*arow + p*sa] * B[p*sb + j*bcol] 267 @unsafe 268 def gemm(+ac: Bool, +n: Nat, +m: Nat, +k: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, a: Array, b: Array, c: Array) -> ABC: 269 run_e(e_of(ac, n, m, k, arow, sa, sb, bcol), U32.from_nat(tiles(n, m)), a, b, c) 270 271 def mm_fin(-n: Nat, -k: Nat, -m: Nat, g: ABC) -> TA.MMul: 272 match g: 273 case ABC{a, b, c}: TA.MMul{TA.Mat{a}, TA.Mat{b}, TA.Mat{c}} 274 275 def mm_nt_fin(-n: Nat, -k: Nat, -m: Nat, g: ABC) -> TA.MMulNT: 276 match g: 277 case ABC{a, b, c}: TA.MMulNT{TA.Mat{a}, TA.Mat{b}, TA.Mat{c}} 278 279 def mm_tn_fin(-n: Nat, -k: Nat, -m: Nat, g: ABC) -> TA.MMulTN: 280 match g: 281 case ABC{a, b, c}: TA.MMulTN{TA.Mat{a}, TA.Mat{b}, TA.Mat{c}} 282 283 # C = A · B (n x k) · (k x m), tensor-array's TA.Mat.matmul bits 284 @unsafe 285 def matmul(+n: Nat, +k: Nat, +m: Nat, par: Nat, a: TA.Mat, b: TA.Mat) -> TA.MMul: 286 match a b: 287 case TA.Mat{x} TA.Mat{y}: 288 mm_fin(n, k, m, gemm_new(n, m, k, U32.from_nat(k), 1, U32.from_nat(m), 1, x, y)) 289 290 # C + A · B, into C (n x m) 291 @unsafe 292 def matmul_acc(+n: Nat, +k: Nat, +m: Nat, par: Nat, a: TA.Mat, b: TA.Mat, c: TA.Mat) -> TA.MMul: 293 match a b c: 294 case TA.Mat{x} TA.Mat{y} TA.Mat{z}: 295 mm_fin(n, k, m, gemm(True{}, n, m, k, U32.from_nat(k), 1, U32.from_nat(m), 1, x, y, z)) 296 297 # C = A · Bᵀ (n x k) · (m x k)ᵀ: for weights stored with one row per output 298 @unsafe 299 def matmul_nt(+n: Nat, +k: Nat, +m: Nat, par: Nat, a: TA.Mat, b: TA.Mat) -> TA.MMulNT: 300 match a b: 301 case TA.Mat{x} TA.Mat{y}: 302 mm_nt_fin(n, k, m, gemm_new(n, m, k, U32.from_nat(k), 1, 1, U32.from_nat(k), x, y)) 303 304 # C + A · Bᵀ, into C (n x m) 305 @unsafe 306 def matmul_nt_acc(+n: Nat, +k: Nat, +m: Nat, par: Nat, a: TA.Mat, b: TA.Mat, c: TA.Mat) -> TA.MMulNT: 307 match a b c: 308 case TA.Mat{x} TA.Mat{y} TA.Mat{z}: 309 mm_nt_fin(n, k, m, gemm(True{}, n, m, k, U32.from_nat(k), 1, 1, U32.from_nat(k), x, y, z)) 310 311 # C = Aᵀ · B (k x n)ᵀ · (k x m): weight gradient, Xᵀ · dY 312 @unsafe 313 def matmul_tn(+n: Nat, +k: Nat, +m: Nat, par: Nat, a: TA.Mat, b: TA.Mat) -> TA.MMulTN: 314 match a b: 315 case TA.Mat{x} TA.Mat{y}: 316 mm_tn_fin(n, k, m, gemm_new(n, m, k, 1, U32.from_nat(n), U32.from_nat(m), 1, x, y)) 317 318 # C + Aᵀ · B, into C (n x m): dW += Xᵀ · dY 319 @unsafe 320 def matmul_tn_acc(+n: Nat, +k: Nat, +m: Nat, par: Nat, a: TA.Mat, b: TA.Mat, c: TA.Mat) -> TA.MMulTN: 321 match a b c: 322 case TA.Mat{x} TA.Mat{y} TA.Mat{z}: 323 mm_tn_fin(n, k, m, gemm(True{}, n, m, k, 1, U32.from_nat(n), U32.from_nat(m), 1, x, y, z)) 324 325 # ===================================================================================================================== 326 # Bangs over any leaf code (0.1.1.0). runK and runK_new (K = 1..4 outputs) run one bang over 2^d leaves: leaf t of 327 # l = 2^d gets t, l, a record e of scalars (the same for every leaf), four input arrays and the K outputs, and returns 328 # the outputs. The leaf code is a template (~f), written by the caller in plain Bend: it takes every array as an ordinary 329 # argument (it reads them with Array.get and may write the outputs), so it needs no @unsafe; the sharing is here, where 330 # every node of the tree holds a handle to each array (Array.fork in Base, upstream's one-block pattern) and joins the 331 # outputs (Array.join). runK takes the outputs from the caller and writes them in place (an optimizer's weights and 332 # moments); runK_new makes them inside the bang with 2^c1.. cells, since an array made on the host pulls device pages 333 # back (see bang_new above). Both return the inputs unchanged with the outputs (RK). An unused input is any array, such 334 # as none(Unit{}). 335 # 336 # Not checked, as for the products: that a leaf writes only its own cells of the outputs (two leaves writing one cell 337 # race) and never an input. Leaves past the work must return their outputs unchanged; share(t, n, l) is how many of the 338 # cells t, t + l, t + 2l, .. lie below n (leaf t's cells when n cells are dealt round-robin, so neighbouring lanes read 339 # neighbouring cells), and depth(n) is a tree for n cells: 2^14 leaves once there are that many, else the fewest that 340 # cover them. The leaves run in a run-time order inside each group of 128 (bendlang/bend's CUDA lowering), so a leaf's 341 # result must not depend on another's. Every number a leaf computes is the caller's own code, so the bits are the 342 # caller's: the same on the GPU as on the CPU threads, except where CUDA's math functions differ from the host's (an exp 343 # on the device takes the float overload; bend-ml's docs/upstream/gpu-math-overloads.md). 344 # ===================================================================================================================== 345 346 def depth2(big: Bool, +n: Nat) -> Nat: 347 match big: 348 case True{}: 14n 349 case False{}: TA.cap_depth(n) 350 351 # leaves for n cells: 2^14 when there are that many, else the fewest powers of two that cover them 352 def depth(+n: Nat) -> Nat: 353 depth2(Nat.is_le(16384n, n), n) 354 355 # how many of the cells t, t + l, t + 2l, .. lie below n 356 def share(+t: U32, +n: U32, +l: U32) -> Nat: 357 U32.to_nat(Bool.pick(U32, U32.is_lt(t, n), U32.div(((n - t : U32) + (l - 1 : U32) : U32), l), 0)) 358 359 # an array of one cell, for an unused input 360 def none(u: Unit) -> Array: 361 Array.new(F32, 0n, 0.0) 362 363 type O2 is Type: 364 O2{a: Array, b: Array} 365 366 type O3 is Type: 367 O3{a: Array, b: Array, c: Array} 368 369 type O4 is Type: 370 O4{a: Array, b: Array, c: Array, d: Array} 371 372 type R1 is Type: 373 R1{x1: Array, x2: Array, x3: Array, x4: Array, a: Array} 374 375 type R2 is Type: 376 R2{x1: Array, x2: Array, x3: Array, x4: Array, o: O2} 377 378 type R3 is Type: 379 R3{x1: Array, x2: Array, x3: Array, x4: Array, o: O3} 380 381 type R4 is Type: 382 R4{x1: Array, x2: Array, x3: Array, x4: Array, o: O4} 383 384 @unsafe 385 def join2(x: O2, y: O2) -> O2: 386 match x y: 387 case O2{a, b} O2{a2, b2}: O2{Array.join(F32, a, a2), Array.join(F32, b, b2)} 388 389 @unsafe 390 def join3(x: O3, y: O3) -> O3: 391 match x y: 392 case O3{a, b, c} O3{a2, b2, c2}: O3{Array.join(F32, a, a2), Array.join(F32, b, b2), Array.join(F32, c, c2)} 393 394 @unsafe 395 def join4(x: O4, y: O4) -> O4: 396 match x y: 397 case O4{a, b, c, d} O4{a2, b2, c2, d2}: O4{Array.join(F32, a, a2), Array.join(F32, b, b2), Array.join(F32, c, c2), Array.join(F32, d, d2)} 398 399 @unsafe 400 def tree1(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> Array -> Array, +d: Nat, +t: U32, +l: U32, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, +oa: Array) -> Array: 401 match d: 402 case 0n: f(t, l, e, x1, x2, x3, x4, oa) 403 case 1n+q: 404 x y = tree1(~E, ~f, q, (t * 2 : U32), l, e, x1, x2, x3, x4, oa) tree1(~E, ~f, q, ((t * 2 : U32) + 1 : U32), l, e, x1, x2, x3, x4, oa) 405 Array.join(F32, x, y) 406 407 @unsafe 408 def root1(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> Array -> Array, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, oa: Array) -> Array: 409 tree1(~E, ~f, d, 0, U32.pow(2, d), e, x1, x2, x3, x4, oa) 410 411 @unsafe 412 def root1_new(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> Array -> Array, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, +ca: Nat) -> Array: 413 root1(~E, ~f, d, e, x1, x2, x3, x4, Array.new(F32, ca, 0.0)) 414 415 @unsafe 416 def go1(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> Array -> Array, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, oa: Array) -> R1: 417 R1{x1, x2, x3, x4, root1!(~E, ~f, d, e, x1, x2, x3, x4, oa)} 418 419 @unsafe 420 def go1_new(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> Array -> Array, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, +ca: Nat) -> R1: 421 R1{x1, x2, x3, x4, root1_new!(~E, ~f, d, e, x1, x2, x3, x4, ca)} 422 423 # one bang over 2^d leaves with 1 output written in place 424 @unsafe 425 def run1(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> Array -> Array, +d: Nat, +e: E, x1: Array, x2: Array, x3: Array, x4: Array, oa: Array) -> R1: 426 go1(~E, ~f, d, e, x1, x2, x3, x4, oa) 427 428 # one bang over 2^d leaves with 1 new output of 2^c.. cells, made on the device 429 @unsafe 430 def run1_new(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> Array -> Array, +d: Nat, +e: E, x1: Array, x2: Array, x3: Array, x4: Array, +ca: Nat) -> R1: 431 go1_new(~E, ~f, d, e, x1, x2, x3, x4, ca) 432 433 @unsafe 434 def tree2(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O2 -> O2, +d: Nat, +t: U32, +l: U32, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, +oa: Array, +ob: Array) -> O2: 435 match d: 436 case 0n: f(t, l, e, x1, x2, x3, x4, O2{oa, ob}) 437 case 1n+q: 438 x y = tree2(~E, ~f, q, (t * 2 : U32), l, e, x1, x2, x3, x4, oa, ob) tree2(~E, ~f, q, ((t * 2 : U32) + 1 : U32), l, e, x1, x2, x3, x4, oa, ob) 439 join2(x, y) 440 441 @unsafe 442 def root2(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O2 -> O2, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, oa: Array, ob: Array) -> O2: 443 tree2(~E, ~f, d, 0, U32.pow(2, d), e, x1, x2, x3, x4, oa, ob) 444 445 @unsafe 446 def root2_new(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O2 -> O2, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, +ca: Nat, +cb: Nat) -> O2: 447 root2(~E, ~f, d, e, x1, x2, x3, x4, Array.new(F32, ca, 0.0), Array.new(F32, cb, 0.0)) 448 449 @unsafe 450 def go2(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O2 -> O2, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, oa: Array, ob: Array) -> R2: 451 R2{x1, x2, x3, x4, root2!(~E, ~f, d, e, x1, x2, x3, x4, oa, ob)} 452 453 @unsafe 454 def go2_new(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O2 -> O2, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, +ca: Nat, +cb: Nat) -> R2: 455 R2{x1, x2, x3, x4, root2_new!(~E, ~f, d, e, x1, x2, x3, x4, ca, cb)} 456 457 # one bang over 2^d leaves with 2 outputs written in place 458 @unsafe 459 def run2(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O2 -> O2, +d: Nat, +e: E, x1: Array, x2: Array, x3: Array, x4: Array, oa: Array, ob: Array) -> R2: 460 go2(~E, ~f, d, e, x1, x2, x3, x4, oa, ob) 461 462 # one bang over 2^d leaves with 2 new outputs of 2^c.. cells, made on the device 463 @unsafe 464 def run2_new(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O2 -> O2, +d: Nat, +e: E, x1: Array, x2: Array, x3: Array, x4: Array, +ca: Nat, +cb: Nat) -> R2: 465 go2_new(~E, ~f, d, e, x1, x2, x3, x4, ca, cb) 466 467 @unsafe 468 def tree3(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O3 -> O3, +d: Nat, +t: U32, +l: U32, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, +oa: Array, +ob: Array, +oc: Array) -> O3: 469 match d: 470 case 0n: f(t, l, e, x1, x2, x3, x4, O3{oa, ob, oc}) 471 case 1n+q: 472 x y = tree3(~E, ~f, q, (t * 2 : U32), l, e, x1, x2, x3, x4, oa, ob, oc) tree3(~E, ~f, q, ((t * 2 : U32) + 1 : U32), l, e, x1, x2, x3, x4, oa, ob, oc) 473 join3(x, y) 474 475 @unsafe 476 def root3(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O3 -> O3, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, oa: Array, ob: Array, oc: Array) -> O3: 477 tree3(~E, ~f, d, 0, U32.pow(2, d), e, x1, x2, x3, x4, oa, ob, oc) 478 479 @unsafe 480 def root3_new(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O3 -> O3, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, +ca: Nat, +cb: Nat, +cc: Nat) -> O3: 481 root3(~E, ~f, d, e, x1, x2, x3, x4, Array.new(F32, ca, 0.0), Array.new(F32, cb, 0.0), Array.new(F32, cc, 0.0)) 482 483 @unsafe 484 def go3(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O3 -> O3, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, oa: Array, ob: Array, oc: Array) -> R3: 485 R3{x1, x2, x3, x4, root3!(~E, ~f, d, e, x1, x2, x3, x4, oa, ob, oc)} 486 487 @unsafe 488 def go3_new(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O3 -> O3, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, +ca: Nat, +cb: Nat, +cc: Nat) -> R3: 489 R3{x1, x2, x3, x4, root3_new!(~E, ~f, d, e, x1, x2, x3, x4, ca, cb, cc)} 490 491 # one bang over 2^d leaves with 3 outputs written in place 492 @unsafe 493 def run3(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O3 -> O3, +d: Nat, +e: E, x1: Array, x2: Array, x3: Array, x4: Array, oa: Array, ob: Array, oc: Array) -> R3: 494 go3(~E, ~f, d, e, x1, x2, x3, x4, oa, ob, oc) 495 496 # one bang over 2^d leaves with 3 new outputs of 2^c.. cells, made on the device 497 @unsafe 498 def run3_new(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O3 -> O3, +d: Nat, +e: E, x1: Array, x2: Array, x3: Array, x4: Array, +ca: Nat, +cb: Nat, +cc: Nat) -> R3: 499 go3_new(~E, ~f, d, e, x1, x2, x3, x4, ca, cb, cc) 500 501 @unsafe 502 def tree4(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O4 -> O4, +d: Nat, +t: U32, +l: U32, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, +oa: Array, +ob: Array, +oc: Array, +od: Array) -> O4: 503 match d: 504 case 0n: f(t, l, e, x1, x2, x3, x4, O4{oa, ob, oc, od}) 505 case 1n+q: 506 x y = tree4(~E, ~f, q, (t * 2 : U32), l, e, x1, x2, x3, x4, oa, ob, oc, od) tree4(~E, ~f, q, ((t * 2 : U32) + 1 : U32), l, e, x1, x2, x3, x4, oa, ob, oc, od) 507 join4(x, y) 508 509 @unsafe 510 def root4(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O4 -> O4, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, oa: Array, ob: Array, oc: Array, od: Array) -> O4: 511 tree4(~E, ~f, d, 0, U32.pow(2, d), e, x1, x2, x3, x4, oa, ob, oc, od) 512 513 @unsafe 514 def root4_new(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O4 -> O4, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, +ca: Nat, +cb: Nat, +cc: Nat, +cd: Nat) -> O4: 515 root4(~E, ~f, d, e, x1, x2, x3, x4, Array.new(F32, ca, 0.0), Array.new(F32, cb, 0.0), Array.new(F32, cc, 0.0), Array.new(F32, cd, 0.0)) 516 517 @unsafe 518 def go4(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O4 -> O4, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, oa: Array, ob: Array, oc: Array, od: Array) -> R4: 519 R4{x1, x2, x3, x4, root4!(~E, ~f, d, e, x1, x2, x3, x4, oa, ob, oc, od)} 520 521 @unsafe 522 def go4_new(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O4 -> O4, +d: Nat, +e: E, +x1: Array, +x2: Array, +x3: Array, +x4: Array, +ca: Nat, +cb: Nat, +cc: Nat, +cd: Nat) -> R4: 523 R4{x1, x2, x3, x4, root4_new!(~E, ~f, d, e, x1, x2, x3, x4, ca, cb, cc, cd)} 524 525 # one bang over 2^d leaves with 4 outputs written in place 526 @unsafe 527 def run4(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O4 -> O4, +d: Nat, +e: E, x1: Array, x2: Array, x3: Array, x4: Array, oa: Array, ob: Array, oc: Array, od: Array) -> R4: 528 go4(~E, ~f, d, e, x1, x2, x3, x4, oa, ob, oc, od) 529 530 # one bang over 2^d leaves with 4 new outputs of 2^c.. cells, made on the device 531 @unsafe 532 def run4_new(~E: Data, ~f: @+t: U32 -> @+l: U32 -> @+e: E -> Array -> Array -> Array -> Array -> O4 -> O4, +d: Nat, +e: E, x1: Array, x2: Array, x3: Array, x4: Array, +ca: Nat, +cb: Nat, +cc: Nat, +cd: Nat) -> R4: 533 go4_new(~E, ~f, d, e, x1, x2, x3, x4, ca, cb, cc, cd) 534 535 # ---- tensor-array 0.1.14.0's precise product on the GPU (opt-in) ---- 536 # TA.Mat.matmul_precise's sums, a bang a block of 1024 steps that starts at the block's place in A and B (oa, ob; no 537 # copy): the first block's bang makes C on the device, each next one adds its cells into C as the accumulating kernel 538 # does (once, as tensor-array's cwr), so every cell is ((x0 + x1) + x2) .., the same bits as tensor-array's 539 # (tests/precise.bend). Only plain (not accumulating) products, as tensor-array's matmul_precise. 540 541 def blk_len(+bi: U32, +k: U32) -> Nat: 542 U32.to_nat(Bool.pick(U32, U32.is_lt((k - (bi * 1024 : U32) : U32), 1024), (k - (bi * 1024 : U32) : U32), 1024)) 543 544 # the blocks after the first: nb more, from block bi, each added into C 545 @unsafe 546 def blocks_acc(nb: Nat, +bi: U32, +n: Nat, +m: Nat, +k: U32, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +a: Array, +b: Array, c: Array) -> Array: 547 match nb: 548 case 0n: c 549 case 1n++q: blocks_acc(q, (bi + 1 : U32), n, m, k, arow, sa, sb, bcol, a, b, bang(e_ofo(True{}, n, m, blk_len(bi, k), arow, sa, sb, bcol, ((bi * 1024 : U32) * sa : U32), ((bi * 1024 : U32) * sb : U32)), TA.cap_depth(tiles(n, m)), a, b, c)) 550 551 @unsafe 552 def precise_c(+n: Nat, +m: Nat, +k: Nat, +arow: U32, +sa: U32, +sb: U32, +bcol: U32, +a: Array, +b: Array) -> Array: 553 blocks_acc(Nat.sub(U32.to_nat(U32.div((U32.from_nat(k) + 1023 : U32), 1024)), 1n), 1, n, m, U32.from_nat(k), arow, sa, sb, bcol, a, b, bang_new(e_ofo(False{}, n, m, 1024n, arow, sa, sb, bcol, 0, 0), TA.cap_depth(tiles(n, m)), a, b, TA.cap_depth(Nat.mul(n, m)))) 554 555 @unsafe 556 def mm_precise3(+n: Nat, +k: Nat, +m: Nat, +a: Array, +b: Array) -> TA.MMul: 557 mm_fin(n, k, m, ABC{a, b, precise_c(n, m, k, U32.from_nat(k), 1, U32.from_nat(m), 1, a, b)}) 558 559 @unsafe 560 def mm_precise2(long: Bool, +n: Nat, +k: Nat, +m: Nat, +par: Nat, a: TA.Mat, b: TA.Mat) -> TA.MMul: 561 match long a b: 562 case True{} TA.Mat{x} TA.Mat{y}: mm_precise3(n, k, m, x, y) 563 case False{} a b: matmul(n, k, m, par, a, b) 564 565 # C = A · B with k summed in blocks of 1024: TA.Mat.matmul_precise's bits (opt-in) 566 @unsafe 567 def matmul_precise(+n: Nat, +k: Nat, +m: Nat, +par: Nat, a: TA.Mat, b: TA.Mat) -> TA.MMul: 568 mm_precise2(Nat.is_lt(1024n, k), n, k, m, par, a, b)