From bbed2c103f839af7ede4b3d86b22b15138480372 Mon Sep 17 00:00:00 2001 From: eunhhu Date: Mon, 5 Oct 2026 21:41:34 +0100 Subject: 041 join ops: the recurrence interpreted and each round's rest as 11 operations apart, interpreted, 57.8 ms (C 32.6; fused by hand 48.4, the tier 2.46x); bfs: a node with more edges than the lanes runs its round alone (tested with lanes of 16, LANES_TEST); runs again --- .../041-batches/2026-10-05-batches.jsonl | 26 ++-- research/experiments/041-batches/bfs.c | 13 +- research/experiments/041-batches/join.c | 152 ++++++++++++++++++++- 3 files changed, 172 insertions(+), 19 deletions(-) (limited to 'research') diff --git a/research/experiments/041-batches/2026-10-05-batches.jsonl b/research/experiments/041-batches/2026-10-05-batches.jsonl index fd9eca6..bbdf939 100644 --- a/research/experiments/041-batches/2026-10-05-batches.jsonl +++ b/research/experiments/041-batches/2026-10-05-batches.jsonl @@ -1,15 +1,15 @@ -{"work":"matmul","c_natural_ms":7.65,"c_ms":7.09,"c_ikj_ms":8.49,"c_ikj_u32_ms":7.43,"c_ikj_vec_ms":4.89,"c_ikj_dvec_ms":7.89,"batch_k_ms":22.07,"batch_j_ms":17.75,"axpy_ms":6.68,"same":true} -{"work":"knapsack","c_natural_ms":6.25,"c_ms":6.37,"c_u32_ms":7.62,"c_vec_ms":9.50,"c_block_ms":1.97,"batch_ms":11.22,"u32_ms":2.71,"u32_ops_ms":4.55,"same":true} +{"work":"matmul","c_natural_ms":7.84,"c_ms":7.11,"c_ikj_ms":8.56,"c_ikj_u32_ms":7.42,"c_ikj_vec_ms":4.89,"c_ikj_dvec_ms":7.93,"batch_k_ms":22.11,"batch_j_ms":17.75,"axpy_ms":6.68,"same":true} +{"work":"knapsack","c_natural_ms":6.27,"c_ms":6.37,"c_u32_ms":7.63,"c_vec_ms":9.50,"c_block_ms":1.96,"batch_ms":11.23,"u32_ms":2.70,"u32_ops_ms":4.55,"same":true} {"work":"roots","c_ms":6.86,"batch_ms":6.79,"same":true} -{"work":"join","c_ms":33.09,"batch_ms":39.36,"interpreted_ms":52.81,"reciprocal_ms":47.64,"same":true} -{"work":"bfs","c_ms":10.01,"c_ahead_ms":7.05,"forms_ms":10.44,"batch_ms":10.28,"reached":99998,"same":true} -{"work":"matmul","c_natural_ms":8.06,"c_ms":7.10,"c_ikj_ms":8.59,"c_ikj_u32_ms":7.41,"c_ikj_vec_ms":4.90,"c_ikj_dvec_ms":7.98,"batch_k_ms":22.08,"batch_j_ms":17.74,"axpy_ms":6.68,"same":true} -{"work":"knapsack","c_natural_ms":6.24,"c_ms":6.38,"c_u32_ms":7.62,"c_vec_ms":9.50,"c_block_ms":1.94,"batch_ms":11.22,"u32_ms":2.70,"u32_ops_ms":4.53,"same":true} -{"work":"roots","c_ms":6.86,"batch_ms":6.79,"same":true} -{"work":"join","c_ms":28.45,"batch_ms":36.14,"interpreted_ms":50.29,"reciprocal_ms":45.72,"same":true} -{"work":"bfs","c_ms":9.96,"c_ahead_ms":6.99,"forms_ms":10.39,"batch_ms":10.34,"reached":99998,"same":true} -{"work":"matmul","c_natural_ms":7.69,"c_ms":7.09,"c_ikj_ms":8.51,"c_ikj_u32_ms":7.42,"c_ikj_vec_ms":4.89,"c_ikj_dvec_ms":7.88,"batch_k_ms":22.06,"batch_j_ms":17.74,"axpy_ms":6.68,"same":true} -{"work":"knapsack","c_natural_ms":6.25,"c_ms":6.38,"c_u32_ms":7.62,"c_vec_ms":9.50,"c_block_ms":1.97,"batch_ms":11.22,"u32_ms":2.71,"u32_ops_ms":4.54,"same":true} +{"work":"join","c_ms":32.59,"batch_ms":38.83,"interpreted_ms":52.10,"reciprocal_ms":47.93,"ops_ms":56.54,"same":true} +{"work":"bfs","c_ms":10.12,"c_ahead_ms":7.16,"forms_ms":10.53,"batch_ms":10.38,"reached":99998,"same":true} +{"work":"matmul","c_natural_ms":7.99,"c_ms":7.37,"c_ikj_ms":8.70,"c_ikj_u32_ms":7.42,"c_ikj_vec_ms":4.90,"c_ikj_dvec_ms":8.14,"batch_k_ms":22.11,"batch_j_ms":17.75,"axpy_ms":6.68,"same":true} +{"work":"knapsack","c_natural_ms":6.26,"c_ms":6.39,"c_u32_ms":7.63,"c_vec_ms":9.51,"c_block_ms":1.95,"batch_ms":11.23,"u32_ms":2.70,"u32_ops_ms":4.55,"same":true} +{"work":"roots","c_ms":6.87,"batch_ms":6.79,"same":true} +{"work":"join","c_ms":34.25,"batch_ms":40.98,"interpreted_ms":55.03,"reciprocal_ms":50.10,"ops_ms":58.01,"same":true} +{"work":"bfs","c_ms":9.97,"c_ahead_ms":7.04,"forms_ms":10.39,"batch_ms":10.27,"reached":99998,"same":true} +{"work":"matmul","c_natural_ms":7.66,"c_ms":7.10,"c_ikj_ms":8.62,"c_ikj_u32_ms":7.41,"c_ikj_vec_ms":4.90,"c_ikj_dvec_ms":8.01,"batch_k_ms":22.06,"batch_j_ms":17.74,"axpy_ms":6.68,"same":true} +{"work":"knapsack","c_natural_ms":6.23,"c_ms":6.37,"c_u32_ms":7.62,"c_vec_ms":9.50,"c_block_ms":1.97,"batch_ms":11.22,"u32_ms":2.70,"u32_ops_ms":4.55,"same":true} {"work":"roots","c_ms":6.87,"batch_ms":6.79,"same":true} -{"work":"join","c_ms":28.77,"batch_ms":36.50,"interpreted_ms":49.19,"reciprocal_ms":44.60,"same":true} -{"work":"bfs","c_ms":9.97,"c_ahead_ms":7.02,"forms_ms":10.34,"batch_ms":10.24,"reached":99998,"same":true} +{"work":"join","c_ms":31.77,"batch_ms":39.92,"interpreted_ms":52.69,"reciprocal_ms":48.44,"ops_ms":57.79,"same":true} +{"work":"bfs","c_ms":9.94,"c_ahead_ms":7.01,"forms_ms":10.31,"batch_ms":10.30,"reached":99998,"same":true} diff --git a/research/experiments/041-batches/bfs.c b/research/experiments/041-batches/bfs.c index 3c25a98..2c61226 100644 --- a/research/experiments/041-batches/bfs.c +++ b/research/experiments/041-batches/bfs.c @@ -12,7 +12,8 @@ * lane's bit read again first (an earlier lane may have added it) * Against natural.h's c_bfs (offsets per node: CSR) and c_bfs_ahead (that * read ahead), on the suite's graph and routes (100,000 nodes, 1,000,000 - * edges drawn, seed 121). Built cc -O2. */ + * edges drawn, seed 121). Built cc -O2; with -DLANES_TEST the edges' lanes + * are 16 (a node with more edges runs its round alone: the answers the same). */ #define _POSIX_C_SOURCE 200809L #include #include @@ -29,7 +30,11 @@ static double now(void) { return t.tv_sec * 1e3 + t.tv_nsec / 1e6; } +#ifdef LANES_TEST +enum { B = 256, LANES = 16 }; +#else enum { B = 256, LANES = 8192 }; +#endif /* seen: list of Seen by node, DENSE */ typedef struct { uint64_t *bits; uint32_t *depth; size_t count; } Seen; @@ -81,6 +86,12 @@ static void batch_bfs(uint32_t from, uint32_t target, uint32_t nodes, double *ou unsigned e = 0, rounds = 0; /* the edges into lanes, as many rounds as fit */ for (; rounds < m && e + (hi[rounds] - lo[rounds]) <= LANES; ++rounds) for (uint32_t k = lo[rounds]; k < hi[rounds]; ++k) edge_to[e] = to_col[k], edge_round[e] = rounds, ++e; + if (!rounds) { /* a node with more edges than the lanes: its round alone, round by round */ + for (uint32_t k = lo[0]; k < hi[0]; ++k) + if (!seen_has(&seen, to_col[k])) seen_add(&seen, to_col[k], depth[0]), queue[count++] = to_col[k]; + ++i; + continue; + } for (unsigned q = 0; q < e; ++q) has[q] = seen_has(&seen, edge_to[q]); /* seen.has(e.to), each lane */ for (unsigned q = 0; q < e; ++q) /* the adds and appends, in order */ if (!has[q] && !seen_has(&seen, edge_to[q])) seen_add(&seen, edge_to[q], depth[edge_round[q]]), queue[count++] = edge_to[q]; diff --git a/research/experiments/041-batches/join.c b/research/experiments/041-batches/join.c index 35b8462..7f53eaf 100644 --- a/research/experiments/041-batches/join.c +++ b/research/experiments/041-batches/join.c @@ -6,7 +6,8 @@ * setup: the appends, the adds in order; the orders: ids gathered, slots * probed, regions gathered, the sums added in order). Against natural.h's * c_join (its own hash table), both timed from nothing to the answer, the - * sizes read at run time. */ + * sizes read at run time. ops: the recurrence interpreted, the rest of each + * round a batch program of operations apart (11 for an order), interpreted. */ #define _POSIX_C_SOURCE 200809L #include #include @@ -117,19 +118,160 @@ static void pure_join(uint32_t customers, uint32_t orders, int interpreted, doub out[0] = bad ? -1 : (double)most, out[1] = (double)sum; free(ids), free(members), free(region); } +/* The same with the rest of each round as a batch program of operations apart, as a compiler would make it from + * the body (no operation fused but floor(a / k)), interpreted: an operation over the batch's lanes at a time. */ +enum { L_IOTA, L_FLOORDIVK, L_MULK, L_SUB, L_ADDK, L_GATHER32, L_PROBE, L_GATHERSLOT32, L_SCATTERADD, L_STORE32, L_INSERT }; +typedef struct { uint8_t op, d, a, b; int64_t k; uint64_t magic; unsigned shift; } LOp; +typedef struct { uint32_t *ids, *region; uint64_t *members; size_t cap; uint32_t customers; int64_t *totals; } Lists; +static void lanes_prepare(LOp *prog, unsigned ops) { + for (unsigned p = 0; p < ops; ++p) { + LOp *o = &prog[p]; + o->magic = 0; + if (o->op != L_FLOORDIVK || o->k <= 2 || !(o->k & (o->k - 1))) continue; + unsigned sh = 0; + while ((2ull << sh) < (uint64_t)o->k) ++sh; + o->magic = (uint64_t)((((unsigned __int128)1 << (64 + sh)) / (uint64_t)o->k) + 1), o->shift = sh; + } +} +__attribute__((optimize("O3"))) static void l_floordivk(int64_t *restrict d, const int64_t *restrict a, const LOp *o, unsigned n) { + for (unsigned l = 0; l < n; ++l) { + int64_t x = a[l], q; + if (x >= 0 && o->magic) q = (int64_t)(((unsigned __int128)(uint64_t)x * o->magic) >> 64 >> o->shift); + else q = x / o->k, q -= (x % o->k != 0) & ((x < 0) != (o->k < 0)); + d[l] = q; + } +} +__attribute__((optimize("O3"))) static void l_mulk(int64_t *restrict d, const int64_t *restrict a, int64_t k, unsigned n, int64_t *over) { + int64_t f = 0; + for (unsigned l = 0; l < n; ++l) { int64_t x = a[l] * k; f |= (x >> 53) ^ (x >> 63); d[l] = x; } + *over |= f; +} +__attribute__((optimize("O3"))) static void l_sub(int64_t *restrict d, const int64_t *restrict a, const int64_t *restrict b, unsigned n, int64_t *over) { + int64_t f = 0; + for (unsigned l = 0; l < n; ++l) { int64_t x = a[l] - b[l]; f |= (x >> 53) ^ (x >> 63); d[l] = x; } + *over |= f; +} +__attribute__((optimize("O3"))) static void l_addk(int64_t *restrict d, const int64_t *restrict a, int64_t k, unsigned n, int64_t *over) { + int64_t f = 0; + for (unsigned l = 0; l < n; ++l) { int64_t x = a[l] + k; f |= (x >> 53) ^ (x >> 63); d[l] = x; } + *over |= f; +} +static void lanes_run(const LOp *prog, unsigned ops, int64_t (*r)[B], unsigned n, Lists *t, int64_t *over) { + for (unsigned p = 0; p < ops; ++p) { + const LOp *o = &prog[p]; + int64_t *d = r[o->d], *a = r[o->a], *b = r[o->b]; + switch (o->op) { + case L_IOTA: for (unsigned l = 0; l < n; ++l) d[l] = o->k + l; break; + case L_FLOORDIVK: l_floordivk(d, a, o, n); break; + case L_MULK: l_mulk(d, a, o->k, n, over); break; + case L_SUB: l_sub(d, a, b, n, over); break; + case L_ADDK: l_addk(d, a, o->k, n, over); break; + case L_GATHER32: /* a list's field at a position (in range: checked) */ + for (unsigned l = 0; l < n; ++l) { uint64_t at = (uint64_t)a[l]; *over |= at >= t->customers; d[l] = t->ids[at < t->customers ? at : 0]; } + break; + case L_PROBE: /* a HASHED list's slot for a key (not there: a failed check) */ + for (unsigned l = 0; l < n; ++l) { + uint64_t bits = key_bits((double)a[l]); + size_t at = slot_of(bits, t->cap - 1); + while (t->members[at] != bits) { + if (t->members[at] == VACANT) { *over |= 1; break; } + at = (at + 1) & (t->cap - 1); + } + d[l] = (int64_t)at; + } + break; + case L_GATHERSLOT32: for (unsigned l = 0; l < n; ++l) d[l] = t->region[a[l]]; break; + case L_SCATTERADD: /* totals[a].amount += b, in order (DENSE by region: a key in range, checked) */ + if (*over) break; + for (unsigned l = 0; l < n; ++l) { + if ((uint64_t)a[l] >= 16) { *over |= 1; break; } + int64_t x = t->totals[a[l]] + b[l]; + *over |= (x >> 53) ^ (x >> 63); + t->totals[a[l]] = x; + } + break; + case L_STORE32: /* appended: the column at k on */ + if (*over) break; + for (unsigned l = 0; l < n; ++l) { *over |= a[l] >> 32; t->ids[o->k + l] = (uint32_t)a[l]; } + break; + case L_INSERT: /* added in order: key a, region b */ + if (*over) break; + for (unsigned l = 0; l < n; ++l) { + uint64_t bits = key_bits((double)a[l]); + size_t at = slot_of(bits, t->cap - 1); + while (t->members[at] != VACANT) at = (at + 1) & (t->cap - 1); + t->members[at] = bits, t->region[at] = (uint32_t)b[l]; + } + break; + } + } +} +static void pure_join_ops(uint32_t customers, uint32_t orders, double *out) { + rounds_prepare(minstd_ops, true); + Lists t = {.customers = customers, .cap = 16}; + while (5 * (size_t)customers > 4 * t.cap) t.cap *= 2; + t.ids = malloc(4 * (size_t)customers), t.members = malloc(8 * t.cap), t.region = malloc(4 * t.cap); + int64_t totals[16] = {0}; + t.totals = totals; + for (size_t i = 0; i < t.cap; ++i) t.members[i] = VACANT; + static int64_t r[8][B]; + LOp setup[] = { + {L_IOTA, 1, 0, 0, 0, 0, 0}, /* i */ + {L_FLOORDIVK, 2, 1, 0, 16, 0, 0}, /* floor(i / 16) */ + {L_MULK, 2, 2, 0, 16, 0, 0}, + {L_SUB, 3, 1, 2, 0, 0, 0}, /* region: i - that */ + {L_STORE32, 0, 0, 0, 0, 0, 0}, /* all.append(id: x) */ + {L_INSERT, 0, 0, 3, 0, 0, 0}, /* byid.add(id: x, region) */ + }; + LOp order[] = { + {L_FLOORDIVK, 1, 0, 0, customers, 0, 0}, /* floor(x / customers) */ + {L_MULK, 1, 1, 0, customers, 0, 0}, + {L_SUB, 1, 0, 1, 0, 0, 0}, /* x - that */ + {L_ADDK, 1, 1, 0, 0, 0, 0}, /* all.first + that: c */ + {L_GATHER32, 2, 1, 0, 0, 0, 0}, /* c.id */ + {L_PROBE, 3, 2, 0, 0, 0, 0}, /* byid[c.id] */ + {L_GATHERSLOT32, 4, 3, 0, 0, 0, 0}, /* .region */ + {L_FLOORDIVK, 5, 0, 0, 1000, 0, 0}, /* floor(x / 1000) */ + {L_MULK, 5, 5, 0, 1000, 0, 0}, + {L_SUB, 5, 0, 5, 0, 0, 0}, /* x - that */ + {L_SCATTERADD, 0, 4, 5, 0, 0, 0}, /* totals[region].amount += that */ + }; + lanes_prepare(setup, 6), lanes_prepare(order, 11); + int64_t over = 0, x = 1; + for (uint32_t i0 = 0; i0 < customers; i0 += B) { + unsigned m = customers - i0 < B ? customers - i0 : B; + int64_t reg[3] = {x}; + over |= rounds_run(minstd_ops, reg, (uint64_t *)r[0], m); + x = reg[0]; + setup[0].k = i0, setup[4].k = i0; + lanes_run(setup, 6, r, m, &t, &over); + } + for (uint32_t o0 = 0; o0 < orders; o0 += B) { + unsigned m = orders - o0 < B ? orders - o0 : B; + int64_t reg[3] = {x}; + over |= rounds_run(minstd_ops, reg, (uint64_t *)r[0], m); + x = reg[0]; + lanes_run(order, 11, r, m, &t, &over); + } + int64_t most = 0, sum = 0; + for (int k = 0; k < 16; ++k) most = totals[k] > most ? totals[k] : most, sum += totals[k]; + out[0] = over ? -1 : (double)most, out[1] = (double)sum; + free(t.ids), free(t.members), free(t.region); +} static int cmpd(const void *a, const void *b) { double x = *(const double *)a, y = *(const double *)b; return (x > y) - (x < y); } static volatile uint32_t sizes[2] = {100000, 1000000}; /* read at run time */ static void (*volatile natural)(const uint32_t *, double *) = c_join; int main(void) { - double tc[11], tp[11], ti[11], tm[11], oc[2], op[2], oi[2], om[2]; + double tc[11], tp[11], ti[11], tm[11], tu[11], oc[2], op[2], oi[2], om[2], ou[2]; uint32_t in[2] = {sizes[0], sizes[1]}; for (int r = 0; r < 11; ++r) { double t = now(); natural(in, oc); tc[r] = now() - t; t = now(); pure_join(in[0], in[1], 0, op); tp[r] = now() - t; t = now(); pure_join(in[0], in[1], 1, oi); ti[r] = now() - t; t = now(); pure_join(in[0], in[1], 2, om); tm[r] = now() - t; + t = now(); pure_join_ops(in[0], in[1], ou); tu[r] = now() - t; } - qsort(tc, 11, 8, cmpd), qsort(tp, 11, 8, cmpd), qsort(ti, 11, 8, cmpd), qsort(tm, 11, 8, cmpd); - printf("{\"work\":\"join\",\"c_ms\":%.2f,\"batch_ms\":%.2f,\"interpreted_ms\":%.2f,\"reciprocal_ms\":%.2f,\"same\":%s}\n", tc[5], tp[5], ti[5], tm[5], - oc[0] == op[0] && oc[1] == op[1] && oc[0] == oi[0] && oc[1] == oi[1] && oc[0] == om[0] && oc[1] == om[1] ? "true" : "false"); + qsort(tc, 11, 8, cmpd), qsort(tp, 11, 8, cmpd), qsort(ti, 11, 8, cmpd), qsort(tm, 11, 8, cmpd), qsort(tu, 11, 8, cmpd); + printf("{\"work\":\"join\",\"c_ms\":%.2f,\"batch_ms\":%.2f,\"interpreted_ms\":%.2f,\"reciprocal_ms\":%.2f,\"ops_ms\":%.2f,\"same\":%s}\n", tc[5], tp[5], ti[5], tm[5], tu[5], + oc[0] == op[0] && oc[1] == op[1] && oc[0] == oi[0] && oc[1] == oi[1] && oc[0] == om[0] && oc[1] == om[1] && oc[0] == ou[0] && oc[1] == ou[1] ? "true" : "false"); } -- cgit v1.2.3