diff --git a/docs/SPEC.md b/docs/SPEC.md index 31cb4836..cdbf9f60 100644 --- a/docs/SPEC.md +++ b/docs/SPEC.md @@ -892,7 +892,19 @@ All in `PM/commonTest`, so they run wherever the interpreter loads. `min(a, b)` compares with `<` and keeps the first element in iteration order on ties (including `+0.0` versus `-0.0`), NaN behaviour depends on position because comparisons with NaN are false, and chunked combination must reproduce that chain bit for bit; if an implementation cannot, they - are left to a later step. **Float `+` and `*` reductions are excluded in v1**: chunked evaluation + are left to a later step. Step 2 implements `+` and `*` only: the compiled subset has no bitwise + operator (the front end skips a function using `&`, `|` or `^` on ints, there is no IR operation), + so `& | ^` are left until the subset has them. + A reduction loop is `acc = acc op e` statements only (no element-wise stores in the same loop, so + a rerun needs no restore), each `acc` a function-local i64 read nowhere else in the loop and + assigned before it, `e` built from i64 constants, locals, `b[i]` of a list or local int array + and i64 `+ - *`. Each chunk computes its partial from the identity with the sequential code's + overflow checks; an overflow in any chunk, in an element expression, or in the combination + discards the region (no accumulator written) and the sequential loop runs from its start state, + which deopts or promotes exactly as it does today, so the result is Python's unbounded int. The + exact integer result of a region that does not overflow equals the sequential one; a region may + complete where the sequential order would have overflowed in an intermediate prefix and the + interpreter's redo would have given the same value (the deopt counter is not part of the identity). **Float `+` and `*` reductions are excluded in v1**: chunked evaluation is a different expression tree and not bit-identical. This is not permanent: a fastmath opt-in (#153) awaits the owner's decision, and N-15 follows it. Partials combine in chunk order. - **Observability.** N-14 holds (a traced call runs the interpreted function, never a region); @@ -914,7 +926,7 @@ All in `PM/commonTest`, so they run wherever the interpreter loads. chunk, NaN/inf/-0.0), identical bytes for every thread count, failure at first/middle/last iteration, negative qualification cases, disabled-path test, reductions, race and GC tests, TSan/ASan, N-14, signals, fork and shutdown, and measurement tests for the threshold. - `Status: step 1 implemented (#277): element-wise F64 stores a[i] = into list copy-in or local arrays; no reductions, no failure path (eligible loops cannot fail or deopt); the rest planned (design)`. + `Status: steps 1 and 2 implemented (#277): element-wise F64 stores a[i] = into list copy-in or local arrays (no failure path: eligible loops cannot fail or deopt); i64 reductions acc += e and acc *= e with the overflow rerun above (stores and reductions never in one loop); & | ^ and min/max reductions, float reductions, the poll inside a region and the rest planned (design)`. --- diff --git a/python-multiplatform-ksp/src/main/python/typedpython/cgen.py b/python-multiplatform-ksp/src/main/python/typedpython/cgen.py index a10cc3b1..9a5a5124 100644 --- a/python-multiplatform-ksp/src/main/python/typedpython/cgen.py +++ b/python-multiplatform-ksp/src/main/python/typedpython/cgen.py @@ -170,6 +170,7 @@ """ from __future__ import annotations +import re from pathlib import Path from typedpython import ir @@ -2145,11 +2146,32 @@ def _emit_region(self, s: ir.ForRange, ts: str, tn: str, mag: int) -> None: "function, bad step or no ParInfo)") scalars: set[str] = set() arrays: set[str] = set() - for st in s.body: - if not isinstance(st, ir.StoreIndex): - raise CGenError(f"{f.name}: parallel loop body has a {type(st).__name__}") - arrays.add(st.array) - self._par_collect(st.value, scalars, arrays) + reds: list[tuple[str, str, ir.Expr]] = [] # (accumulator, op name, the other operand) + if info.reduce: + if info.stored or len(info.reduce) != len(s.body): + raise CGenError(f"{f.name}: parallel reduction loop with stores or a missing reduction") + for st, (acc, opn) in zip(sorted(s.body, key=lambda x: getattr(x, "target", "")), info.reduce): + v = st.value if isinstance(st, ir.Assign) else None + if (not isinstance(v, ir.BinOp) or st.target != acc or v.op.name != opn + or opn not in ("ADD", "MUL") or v.type is not Type.I64 + or self.types.get(acc) is not Type.I64 or acc == s.var): + raise CGenError(f"{f.name}: parallel reduction statement is not `acc = acc op e`") + if isinstance(v.left, ir.Local) and v.left.name == acc: + other = v.right + elif isinstance(v.right, ir.Local) and v.right.name == acc: + other = v.left + else: + raise CGenError(f"{f.name}: parallel reduction of {acc} does not read {acc}") + reds.append((acc, opn, other)) + self._par_collect(other, scalars, arrays) + if any(x in scalars for x, _, _ in reds): + raise CGenError(f"{f.name}: a parallel reduction reads its accumulator") + else: + for st in s.body: + if not isinstance(st, ir.StoreIndex): + raise CGenError(f"{f.name}: parallel loop body has a {type(st).__name__}") + arrays.add(st.array) + self._par_collect(st.value, scalars, arrays) scalars.discard(s.var) for a in arrays: if a not in self.arrays and a not in self.local_arrays: @@ -2180,6 +2202,8 @@ def _emit_region(self, s: ir.ForRange, ts: str, tn: str, mag: int) -> None: fields.append(f" {ct} {ident};") decls.append(f" {ct} {ident} = c->{ident};") sets.append(f"tp_pr.{ident} = {ident};") + if reds: + fields += [" unsigned char *tp_pst;"] + [f" int64_t *tp_part{j};" for j in range(len(reds))] saved = (self.lines, self.temps, self.depth, self.brk) self.lines, self.temps, self.depth, self.brk = [], [], 2, [] try: @@ -2188,30 +2212,73 @@ def _emit_region(self, s: ir.ForRange, ts: str, tn: str, mag: int) -> None: finally: self.lines, self.temps, self.depth, self.brk = saved text = "\n".join(body_lines) + if reds: + # the one failure path a reduction body has: the sequential overflow check, which here + # marks the chunk bad (the worker's own tp_rc / tp_exit); nothing else may remain + text = re.sub(r"^\s*tp_s = tp_(add|sub|mul)_i64\(.*\);\n?", "", text, flags=re.M) + text = re.sub(r"^\s*if \(tp_s != 0\) \{ tp_rc = tp_s; goto tp_exit; \}\n?", "", text, flags=re.M) for banned in ("goto", "tp_rc", "tp_s ", "return", "tp_st", "tp_dict", "PyObject", "tp_release"): if banned in text: raise CGenError(f"{f.name}: parallel loop body emitted `{banned}`") sign_var = _ident("l", s.var) out = [f"typedef struct {{", *fields, f"}} {ctx_t};", "", - f"static void {body_f}(const void *tp_c, uint64_t tp_k0, uint64_t tp_k1)", "{", + f"static void {body_f}(const void *tp_c, unsigned tp_ci, uint64_t tp_k0, uint64_t tp_k1)", "{", f" const {ctx_t} *c = (const {ctx_t} *)tp_c;", " uint64_t tp_k;", *decls] + if reds: + out += [" int tp_rc = 0;", " int tp_s = 0;"] + out += [f" int64_t {_ident('l', acc)} = {'1' if opn == 'MUL' else '0'};" for acc, opn, _ in reds] + out.append(" (void)tp_rc; (void)tp_s;") for ct, n in body_temps: out.append(f" {ct}{'' if ct.endswith('*') else ' '}{n} = 0;") out.append(" (void)c;" + "".join(f" (void){_ident('l', x)};" for x in sorted(scalars))) out.append(f" for (tp_k = tp_k0; tp_k < tp_k1; tp_k++) {{") out.append(f" {sign_var} = (int64_t)((uint64_t)c->tp_ts + tp_k * UINT64_C({mag}));") out += body_lines - out += [" }", "}", ""] + out.append(" }") + if reds: + out += [f" c->tp_part{j}[tp_ci] = {_ident('l', acc)};" for j, (acc, _, _) in enumerate(reds)] + out += [" c->tp_pst[tp_ci] = TP_PAR_OK;", " return;", "tp_exit:", + " c->tp_pst[tp_ci] = TP_PAR_BAD;"] + out += ["}", ""] self.par_defs.append("\n".join(out)) w = max(1, info.weight) + pok = self.tmp("int") if reds else None + if reds: + self.emit(f"{pok} = 0;") self.emit(f"if ({tn} >= (uint64_t)((TP_PAR_MIN_WORK + {w - 1}) / {w}) && tp_par_threads() > 1) {{") self.depth += 1 self.emit(f"{ctx_t} tp_pr;") + if reds: + self.emit("unsigned char tp_pst[TP_PAR_HARD_CAP] = {0};") + for j in range(len(reds)): + self.emit(f"int64_t tp_pp{j}[TP_PAR_HARD_CAP];") self.emit(f"tp_pr.tp_ts = {ts};") for line in sets: self.emit(line) + if reds: + self.emit("tp_pr.tp_pst = tp_pst;") + for j in range(len(reds)): + self.emit(f"tp_pr.tp_part{j} = tp_pp{j};") self.emit(f"(void)tp_par_run({body_f}, &tp_pr, {tn});") + if reds: + # Combine in chunk order with the sequential checks. All or none: an overflow anywhere + # (a chunk, or here) leaves every accumulator untouched and the sequential loop runs. + self.emit(f"{pok} = 1;") + for j, (acc, opn, _) in enumerate(reds): + self.emit(f"int64_t tp_rv{j} = {_ident('l', acc)};") + self.emit(f"{pok} = {pok} && tp_par_combine(TP_PAR_OP_{opn}, &tp_rv{j}, tp_pp{j}, tp_pst) == 0;") + self.emit(f"if ({pok}) {{") + for j, (acc, _, _) in enumerate(reds): + self.emit(f" {_ident('l', acc)} = tp_rv{j};") + self.emit(f" {_ident('l', s.var)} = (int64_t)((uint64_t){ts} + ({tn} - 1u) * UINT64_C({mag}));") + self.emit("} else {") + self.emit(" tp_par_rerun_note();") + self.emit("}") + self.depth -= 1 + self.emit("}") + self.emit(f"if (!{pok}) {{") + return self.emit(f"{_ident('l', s.var)} = (int64_t)((uint64_t){ts} + ({tn} - 1u) * UINT64_C({mag}));") self.depth -= 1 self.emit("} else {") @@ -2756,14 +2823,14 @@ def _state_functions(mg: _ModGen) -> str: def _par_stats_function() -> str: """`__typedpython_parallel__()`: (regions that used workers, workers started, failed thread - creations) since import, so a test can prove the parallel path ran (SPEC N-15).""" + creations, reduction regions rerun sequentially after an overflow) since import, so a test can prove the parallel path ran (SPEC N-15).""" return "\n".join([ "static PyObject *tp_par_stats_fn(PyObject *tp_module, PyObject *tp_unused)", "{", - " int64_t tp_v[3];", + " int64_t tp_v[4];", " (void)tp_module; (void)tp_unused;", " tp_par_stats(tp_v);", - " return Py_BuildValue(\"(LLL)\", (long long)tp_v[0], (long long)tp_v[1], (long long)tp_v[2]);", + " return Py_BuildValue(\"(LLLL)\", (long long)tp_v[0], (long long)tp_v[1], (long long)tp_v[2], (long long)tp_v[3]);", "}", "", "static PyMethodDef tp_par_stats_def = {\"__typedpython_parallel__\", tp_par_stats_fn, METH_NOARGS, NULL};", diff --git a/python-multiplatform-ksp/src/main/python/typedpython/ir.py b/python-multiplatform-ksp/src/main/python/typedpython/ir.py index ba31d3c9..723ed963 100644 --- a/python-multiplatform-ksp/src/main/python/typedpython/ir.py +++ b/python-multiplatform-ksp/src/main/python/typedpython/ir.py @@ -653,7 +653,7 @@ class While(Stmt): @dataclass(frozen=True) class ParInfo: """Why a `ForRange` may run its iterations on several threads (SPEC N-15, step 1): the body is - only `a[var] = ` stores, `stored` are the arrays stored into (each at the loop + only `a[var] = ` stores (or, step 2, only i64 `acc op= e` reductions: see `reduce`), `stored` are the arrays stored into (each at the loop variable, proven in bounds, so distinct iterations write distinct elements) and `read` the arrays read, all at the loop variable and proven in bounds. The expression nodes are Const, scalar Local (not assigned in the body), ToFloat of an I64, F64 + - * and negation, so no node can fail, @@ -664,6 +664,11 @@ class ParInfo: stored: tuple[str, ...] read: tuple[str, ...] weight: int = 1 + # Step 2 (#277): the i64 reductions `acc = acc op e` of a loop whose body is only such statements, + # as sorted (acc, BinOpKind name) pairs, op ADD or MUL. Then `stored` is empty. Each chunk computes + # its partial from the op's identity with the sequential code's overflow checks; an overflow in any + # chunk or in the combination discards the region and the sequential loop runs from its start. + reduce: tuple[tuple[str, str], ...] = () @dataclass(frozen=True) diff --git a/python-multiplatform-ksp/src/main/python/typedpython/runtime/API.md b/python-multiplatform-ksp/src/main/python/typedpython/runtime/API.md index 52dc203e..a3de3edb 100644 --- a/python-multiplatform-ksp/src/main/python/typedpython/runtime/API.md +++ b/python-multiplatform-ksp/src/main/python/typedpython/runtime/API.md @@ -492,7 +492,7 @@ for the epoch's increment. `tp_runtime.h` needs the GCC/Clang `__atomic` builtin ## Parallel regions (SPEC N-15, issues #151, #277) A `ForRange` the verifier marked `par` (`ir.ParInfo`; marked function, `a[i] = F64 expression` stores -only) runs its iterations on per-region C threads above a work threshold. Workers touch no Python +or i64 `+ *` reductions) runs its iterations on per-region C threads above a work threshold. Workers touch no Python object and no Python API, are not Python threads and never outlive the call. #define TP_PAR_MIN_WORK 65536 /* PLACEHOLDER, to be measured on #151: weighted iterations (trip count * ParInfo.weight) */ #define TP_PAR_MAX_THREADS 8 /* default cap; TYPEDPYTHON_THREADS overrides (1 disables), hard cap 64 */ @@ -501,21 +501,33 @@ object and no Python API, are not Python threads and never outlive the call. #define TP_PAR_ALIGN 8u /* chunk starts are multiples of 8 iterations (one cache line of doubles) */ #define TP_PAR_HAVE_THREADS /* defined with pthreads; -DTP_NO_THREADS compiles the threads out */ #define TP_PAR_TEST_FAIL_CREATE /* test hook: every pthread_create reports failure */ - typedef void (*tp_par_fn)(const void *ctx, uint64_t lo, uint64_t hi); + typedef void (*tp_par_fn)(const void *ctx, unsigned ci, uint64_t lo, uint64_t hi); /* ci: chunk index, 0.. */ int tp_par_threads(void); /* the cap, read once from TYPEDPYTHON_THREADS, else min(cores, 8) */ int tp_par_run(tp_par_fn fn, const void *ctx, uint64_t n); /* runs [0, n) in contiguous chunks, joined before return; returns workers started */ - void tp_par_stats(int64_t out[3]); /* regions with workers, workers started, failed creations */ + void tp_par_stats(int64_t out[4]); /* regions with workers, workers started, failed creations, reduction reruns */ + #define TP_PAR_NONE 0 / TP_PAR_OK 1 / TP_PAR_BAD 2 /* per-chunk state of a reduction chunk (zero-initialised array of TP_PAR_HARD_CAP) */ + #define TP_PAR_OP_ADD 0 / TP_PAR_OP_MUL 1 + int tp_par_combine(int op, int64_t *acc, const int64_t *part, const unsigned char *st); /* acc op= part[0], part[1], ... in chunk order, checked; 1 if a chunk was BAD or it overflowed (acc untouched) */ + void tp_par_rerun_note(void); /* counts a discarded reduction region */ `tp_par_run` always completes the whole range. The calling thread runs the first chunk; a chunk whose `pthread_create` failed runs on the calling thread after it (fewer workers, down to none: the sequential loop, no exception, not a semantic event). Workers start with all signals blocked, so Python signal handlers run on Python threads only. The chunk function must be safe for any chunk on any thread; generated code passes a context struct of raw array pointers and scalar copies. A module with a `par` loop also exports `__typedpython_parallel__()` returning `(regions, workers_started, -failed_creations)`, so tests can prove the parallel path ran. Free-threaded state: the four counters and +failed_creations, reduction_reruns)`, so tests can prove the parallel path ran. Free-threaded state: the four counters and the thread cap are atomics (relaxed); a region's chunk bounds, jobs and thread handles live on the caller's stack. The caller does not poll the eval breaker inside a region (array functions do not poll anyway, #152). +Reduction regions (step 2). The chunk function starts each accumulator at the op's identity (0 for `+`, +1 for `*`), runs the sequential statements `acc = acc op e` with the sequential overflow checks, and +writes `part[ci]` and `st[ci] = TP_PAR_OK`; a failed check sets `st[ci] = TP_PAR_BAD` and stops the chunk. +After `tp_par_run` the caller combines every accumulator with `tp_par_combine`; if any accumulator fails +(a BAD chunk or an overflow in the combination) no accumulator and no loop variable is written, the rerun +counter is incremented and the sequential loop runs from the loop's start state, which deopts or not +exactly as today. (`& | ^` and min/max are not reductions yet: the compiled subset has no bitwise operator.) + ## Development-mode import gate (Python API, not C; SPEC N-7b) typedpython.importgate.install(source_roots) # meta path finder; checks project modules before they run typedpython.importgate.uninstall() diff --git a/python-multiplatform-ksp/src/main/python/typedpython/runtime/tp_runtime.h b/python-multiplatform-ksp/src/main/python/typedpython/runtime/tp_runtime.h index 57e4334d..0d11f897 100644 --- a/python-multiplatform-ksp/src/main/python/typedpython/runtime/tp_runtime.h +++ b/python-multiplatform-ksp/src/main/python/typedpython/runtime/tp_runtime.h @@ -2943,12 +2943,13 @@ static inline int tp_traced(PyObject *module, PyThreadState *ts, PyCodeObject *c #include #endif -typedef void (*tp_par_fn)(const void *ctx, uint64_t lo, uint64_t hi); +typedef void (*tp_par_fn)(const void *ctx, unsigned ci, uint64_t lo, uint64_t hi); static int tp_par_cap_v = 0; /* 0: not read yet; then the thread cap (>= 1), set once */ static int64_t tp_par_regions_v = 0; /* regions that ran with at least one worker */ static int64_t tp_par_started_v = 0; /* worker threads started */ static int64_t tp_par_failed_v = 0; /* pthread_create calls that failed (their chunk ran on the caller) */ +static int64_t tp_par_reruns_v = 0; /* reduction regions discarded (overflow) and rerun sequentially */ /* The thread cap: TYPEDPYTHON_THREADS (a positive integer, 1 disables), else min(cores, 8); read once, * immutable afterwards. Two threads racing on the first call compute the same value. */ @@ -2987,21 +2988,23 @@ static inline int tp_par_threads(void) return c; } -/* out[0] regions that used workers, out[1] workers started, out[2] failed thread creations */ -static inline void tp_par_stats(int64_t out[3]) +/* out[0] regions that used workers, out[1] workers started, out[2] failed thread creations, + * out[3] reduction regions discarded and rerun sequentially */ +static inline void tp_par_stats(int64_t out[4]) { + out[3] = __atomic_load_n(&tp_par_reruns_v, __ATOMIC_RELAXED); out[0] = __atomic_load_n(&tp_par_regions_v, __ATOMIC_RELAXED); out[1] = __atomic_load_n(&tp_par_started_v, __ATOMIC_RELAXED); out[2] = __atomic_load_n(&tp_par_failed_v, __ATOMIC_RELAXED); } #ifdef TP_PAR_HAVE_THREADS -typedef struct { tp_par_fn fn; const void *ctx; uint64_t lo, hi; } tp_par_job; +typedef struct { tp_par_fn fn; const void *ctx; unsigned ci; uint64_t lo, hi; } tp_par_job; static void *tp_par_worker(void *p) { const tp_par_job *j = (const tp_par_job *)p; - j->fn(j->ctx, j->lo, j->hi); + j->fn(j->ctx, j->ci, j->lo, j->hi); return NULL; } #endif @@ -3033,6 +3036,7 @@ static inline int tp_par_run(tp_par_fn fn, const void *ctx, uint64_t n) } jobs[jn].fn = fn; jobs[jn].ctx = ctx; + jobs[jn].ci = (unsigned)jn; jobs[jn].lo = lo; jobs[jn].hi = lo + chunk < n ? lo + chunk : n; live[jn] = 0; @@ -3056,10 +3060,10 @@ static inline int tp_par_run(tp_par_fn fn, const void *ctx, uint64_t n) } } pthread_sigmask(SIG_SETMASK, &old, NULL); - fn(ctx, jobs[0].lo, jobs[0].hi); + fn(ctx, 0, jobs[0].lo, jobs[0].hi); for (k = 1; k < jn; k++) { if (!live[k]) { - fn(ctx, jobs[k].lo, jobs[k].hi); /* creation failed: this chunk runs here */ + fn(ctx, (unsigned)k, jobs[k].lo, jobs[k].hi); /* creation failed: this chunk runs here */ } } for (k = 1; k < jn; k++) { @@ -3074,8 +3078,47 @@ static inline int tp_par_run(tp_par_fn fn, const void *ctx, uint64_t n) return started; } #endif - fn(ctx, 0, n); + fn(ctx, 0, 0, n); return 0; } +/* Reductions (step 2, #277). A reduction chunk function writes its partial (the op applied from the + * op's identity, with the sequential code's overflow checks) into part[ci] and sets st[ci] to + * TP_PAR_OK, or sets st[ci] = TP_PAR_BAD when a check failed (the partial is then meaningless). + * st[] is zero-initialised by the caller, so chunks that never ran stay TP_PAR_NONE. Chunks are a + * prefix 0..m-1 of the 64 slots. */ +#define TP_PAR_NONE 0 +#define TP_PAR_OK 1 +#define TP_PAR_BAD 2 +#define TP_PAR_OP_ADD 0 +#define TP_PAR_OP_MUL 1 + +static inline int64_t tp_par_identity(int op) { return op == TP_PAR_OP_MUL ? 1 : 0; } + +/* acc = acc op part[0] op part[1] ... in chunk order, checked like the sequential code. Returns 0 and + * stores the result in *acc, or 1 (a chunk failed or the combination overflowed) leaving *acc alone. */ +static inline int tp_par_combine(int op, int64_t *acc, const int64_t *part, const unsigned char *st) +{ + int64_t v = *acc, w; + int k; + for (k = 0; k < TP_PAR_HARD_CAP; k++) { + if (st[k] == TP_PAR_BAD) { + return 1; + } + } + for (k = 0; k < TP_PAR_HARD_CAP && st[k] == TP_PAR_OK; k++) { + if (op == TP_PAR_OP_MUL ? tp_mul_i64(v, part[k], &w) : tp_add_i64(v, part[k], &w)) { + return 1; + } + v = w; + } + *acc = v; + return 0; +} + +static inline void tp_par_rerun_note(void) +{ + __atomic_fetch_add(&tp_par_reruns_v, 1, __ATOMIC_RELAXED); +} + #endif /* TP_RUNTIME_H */ diff --git a/python-multiplatform-ksp/src/main/python/typedpython/verify.py b/python-multiplatform-ksp/src/main/python/typedpython/verify.py index 1e34f73e..15312762 100644 --- a/python-multiplatform-ksp/src/main/python/typedpython/verify.py +++ b/python-multiplatform-ksp/src/main/python/typedpython/verify.py @@ -212,8 +212,16 @@ of `b[var]` where `b` is another array cannot alias a written element (array parameters are distinct lists by the entry guard, local arrays are distinct allocations), and a read of the stored array itself is the same element the iteration writes. -Anything else compiles sequentially, unchanged. This cannot find other eligible shapes (reductions, -stencils, calls, nested tiles); those are later steps. + * Step 2 (reductions): or the body is a non-empty sequence of `acc = acc op e` (`_par_reduce`): op + ADD or MUL, each `acc` a distinct i64 local, definitely assigned before the loop, not an array, + loop variable or entry global, and read nowhere else in the body (not in any `e`); `e` built from + i64 Const, i64 Local (the loop variable, or a scalar nothing in the body assigns), `b[var]` + proven on an i64 non-buffer list or local array, and i64 + - *. Stores and reductions are never + mixed in one loop. Those nodes can overflow, which is not a failure of the region: the runtime + discards it and reruns the loop sequentially. (& | ^ are not in the compiled subset: the front + end lowers no bitwise operator, so there is no such reduction to mark; min/max are later.) +Anything else compiles sequentially, unchanged. This cannot find other eligible shapes (stencils, +calls, nested tiles, min/max); those are later steps. Tuple. `Tuple(OBJ, elements)`: a tuple of OBJ expressions (a scalar needs Box), result OBJ, a reference creator in the ownership rule. Building one is not an effect. @@ -1634,14 +1642,15 @@ def s_ForRange(self, s: ForRange, state: _State): after: _State = state for br in loop.breaks: after = _meet(after, br) - par = self._par_info(s, step_ok, body, loop) + par = self._par_info(s, step_ok, body, loop, state) if s.par is not None and s.par != par: self.error(PAR_LOOP, f"ForRange over '{var}' declares par={s.par!r}, but the verifier " f"computes {par!r} (SPEC N-15)") return dataclasses.replace(s, start=start, stop=stop, body=body, par=par), after - def _par_info(self, s: ForRange, step_ok: bool, body: tuple, loop: _Loop) -> ParInfo | None: - """SPEC N-15 step 1: the ParInfo of this loop, or None (see "Parallel loops" above).""" + def _par_info(self, s: ForRange, step_ok: bool, body: tuple, loop: _Loop, + state: _State = None) -> ParInfo | None: + """SPEC N-15 steps 1 and 2: the ParInfo of this loop, or None (see "Parallel loops" above).""" if self.fn.parallel is not True or not step_ok or loop.bounds_array is None: return None if self.diags or not isinstance(body, tuple) or not body: @@ -1651,9 +1660,9 @@ def _par_info(self, s: ForRange, step_ok: bool, body: tuple, loop: _Loop) -> Par read: list[str] = [] weight = [0] - def array_ok(name: str) -> bool: + def array_ok(name: str, want: Type = Type.F64_ARRAY) -> bool: t = self.arrays.get(name) or self.larrays.get(name) - return t is Type.F64_ARRAY and name not in self.buffers + return t is want and name not in self.buffers def at_var(ix: object, proven: object) -> bool: return proven is True and type(ix) is Local and ix.name == var and ix.type is Type.I64 @@ -1685,6 +1694,8 @@ def ok(e: object) -> bool: return e.type is Type.F64 and e.op is UnaryOpKind.NEG and ok(e.operand) return False + if type(body[0]) is Assign: + return self._par_reduce(s, body, state, weight, array_ok, at_var) for st in body: if type(st) is not StoreIndex or not array_ok(st.array) or not at_var(st.index, st.proven): return None @@ -1694,6 +1705,64 @@ def ok(e: object) -> bool: stored.append(st.array) return ParInfo(tuple(sorted(set(stored))), tuple(sorted(set(read))), weight[0]) + def _par_reduce(self, s: ForRange, body: tuple, state: _State, weight: list, array_ok, + at_var) -> ParInfo | None: + """SPEC N-15 step 2: a body of only `acc = acc op e` (op ADD or MUL, either operand order), each + `acc` a distinct i64 local of this function (not a loop variable, array or entry global) that + is definitely assigned before the loop and appears nowhere else in the body, and `e` built + from i64 Const, i64 Local (the loop variable, or any scalar the body does not assign), `b[var]` + proven on an i64 non-buffer list or local array, and i64 + - * . Any of those nodes can overflow: + the runtime then discards the region and reruns the loop sequentially. No store is allowed in + the same loop (a rerun then needs no restore).""" + var = s.var + if state is None: + return None + accs: dict[str, str] = {} + for st in body: + if (type(st) is not Assign or type(st.value) is not BinOp or st.value.type is not Type.I64 + or st.value.op not in (BinOpKind.ADD, BinOpKind.MUL) or st.value.np + or st.target in accs): + return None + accs[st.target] = st.value.op.name + for acc in accs: + if (self.scalars.get(acc) is not Type.I64 or acc == var or acc in self.arrays + or acc in self.larrays or acc in self.entry or acc not in state): + return None + read: list[str] = [] + + def ok(e: object) -> bool: + weight[0] += 1 + t = type(e) + if t is Const: + return e.type is Type.I64 and type(e.value) is int + if t is Local: + return (e.type is Type.I64 and e.name not in accs and e.name not in self.arrays + and e.name not in self.larrays + and (e.name == var or e.name in self.scalars or e.name in self.entry)) + if t is Index: + if e.type is Type.I64 and array_ok(e.array, Type.I64_ARRAY) and at_var(e.index, e.proven): + read.append(e.array) + return True + return False + if t is BinOp: + return (e.type is Type.I64 and e.left.type is Type.I64 and e.right.type is Type.I64 + and e.op in (BinOpKind.ADD, BinOpKind.SUB, BinOpKind.MUL) and not e.np + and ok(e.left) and ok(e.right)) + return False + + for st in body: + weight[0] += 2 + b = st.value + if type(b.left) is Local and b.left.name == st.target and b.left.type is Type.I64: + other = b.right + elif type(b.right) is Local and b.right.name == st.target and b.right.type is Type.I64: + other = b.left + else: + return None + if not ok(other): + return None + return ParInfo((), tuple(sorted(set(read))), weight[0], tuple(sorted(accs.items()))) + def s_Return(self, s: Return, state: _State): returns = self.fn.returns if s.value is None: diff --git a/python-multiplatform-ksp/src/test/python/test_parallel_loops.py b/python-multiplatform-ksp/src/test/python/test_parallel_loops.py index 1aa28652..3b4f037d 100644 --- a/python-multiplatform-ksp/src/test/python/test_parallel_loops.py +++ b/python-multiplatform-ksp/src/test/python/test_parallel_loops.py @@ -417,3 +417,424 @@ def sh(): assert r["gil"] is False if (os.cpu_count() or 1) > 1: assert r["regions"] >= 6 * 6 + + +# ===================================================================================================== +# Step 2 (SPEC N-15, #277): i64 reductions `acc op= expr`, op in + and *. +# (& | ^ are not in the compiled subset: the front end skips `s &= x`, the IR has no bitwise BinOpKind.) +# +# Scope chosen: the body of a reduction loop is only reductions (no element-wise stores in the same +# loop), so a rerun needs no array restore. Each chunk computes its partial from the op's identity with +# the sequential code's own overflow checks; any chunk overflow or a combination overflow discards the +# whole region and the sequential loop runs from its start state (it deopts as today). +# +# Red record: written before the implementation. Before it, `ir.ParInfo` had no `reduce`, reduction +# loops were never `par` (so the region counter did not move and `__typedpython_parallel__()` had +# three fields, no rerun counter); see the report in #277. +# ===================================================================================================== + +RSOURCE = '''\ +# typedpython: compiled + +# typedpython: parallel +def r_sum(xs: list[int]) -> int: + s: int = 0 + for i in range(len(xs)): + s += xs[i] + return s + + +# typedpython: parallel +def r_sum_init(xs: list[int], s0: int) -> int: + s: int = s0 + for i in range(len(xs)): + s += xs[i] + return s + + +# typedpython: parallel +def r_prod(xs: list[int], c: int) -> int: + s: int = 1 + for i in range(len(xs)): + s *= xs[i] * c + return s + + +# typedpython: parallel +def r_sum_mul(xs: list[int], c: int) -> int: + s: int = 0 + for i in range(len(xs)): + s += xs[i] * c + return s + + +# typedpython: parallel +def r_two(xs: list[int]) -> int: + a: int = 0 + b: int = 7 + for i in range(len(xs)): + a += xs[i] + b += xs[i] * 2 + return a * 3 + b + + +# typedpython: parallel +def r_range_sum(xs: list[int]) -> int: + s: int = 5 + for i in range(len(xs)): + s += i - 2 + return s + + +# typedpython: parallel +def r_local(n: int) -> int: + t = [0] * n + for i in range(len(t)): + t[i] = i * 3 + 1 + s: int = 0 + for i in range(len(t)): + s += t[i] + return s + + +def r_unmarked(xs: list[int]) -> int: + s: int = 0 + for i in range(len(xs)): + s += xs[i] + return s + + +# typedpython: parallel +def n_acc_read(xs: list[int]) -> int: + s: int = 0 + for i in range(len(xs)): + s += xs[i] + s + return s + + +# typedpython: parallel +def n_two_ops(xs: list[int]) -> int: + s: int = 0 + for i in range(len(xs)): + s += xs[i] + s *= xs[i] + return s + + +# typedpython: parallel +def n_float(xs: list[float]) -> float: + s: float = 0.0 + for i in range(len(xs)): + s += xs[i] + return s + + +# typedpython: parallel +def n_copy(xs: list[int]) -> int: + s: int = 0 + w: int = 0 + for i in range(len(xs)): + s += xs[i] + w = s + return s + w + + +# typedpython: parallel +def n_mixed(xs: list[int]) -> int: + t = [0.0] * len(xs) + s: int = 0 + for i in range(len(xs)): + t[i] = t[i] + 1.5 + s += xs[i] + return s + + +# typedpython: parallel +def n_div(xs: list[int], d: int) -> int: + s: int = 0 + for i in range(len(xs)): + s += xs[i] // d + return s + + +# typedpython: parallel +def n_shift_step(xs: list[int]) -> int: + s: int = 0 + for i in range(0, len(xs), 2): + s += xs[i + 1] + return s +''' + +RSIZES = [0, 1, 2, 7, 1000, 200003, 1048576] +I64_MAX = 2**63 - 1 +I64_MIN = -(2**63) + + +@pytest.fixture(scope="module") +def rbuilt(tmp_path_factory): + return build(tmp_path_factory.mktemp("rpar"), RSOURCE, name="rparmod") + + +@pytest.fixture(scope="module") +def rmod(rbuilt): + return cbuild.load("rparmod", rbuilt[2]) + + +def rdata(n, seed=5, lo=-1000, hi=1000): + import random + r = random.Random(seed) + return [r.randint(lo, hi) for _ in range(n)] + + +def reruns(m): + return m.__typedpython_parallel__()[3] + + +def same(rmod, name, *args): + """compiled result == interpreter result (exact Python int), with fresh copies of list arguments""" + want = rmod.__typedpython_interpreted__[name](*[list(a) if isinstance(a, list) else a for a in args]) + got = getattr(rmod, name)(*[list(a) if isinstance(a, list) else a for a in args]) + assert got == want and type(got) is type(want), (name, got, want) + return got + + +RED = ["r_sum", "r_sum_init", "r_prod", "r_sum_mul", "r_two", "r_range_sum", "r_local"] +EXPECT_REDUCE = { + "r_sum": (("s", "ADD"),), "r_prod": (("s", "MUL"),), "r_sum_mul": (("s", "ADD"),), + "r_two": (("a", "ADD"), ("b", "ADD")), +} + + +@pytest.mark.parametrize("name", RED) +def test_reduction_loops_are_marked_par(rbuilt, name): + f = fn(rbuilt, name) + loops = [l for l in par_loops(f) if l.par is not None] + assert loops, name + for nm, want in EXPECT_REDUCE.items(): + if nm == name: + assert loops[-1].par.reduce == want and loops[-1].par.stored == () + + +@pytest.mark.parametrize("name", ["n_acc_read", "n_two_ops", "n_float", "n_copy", "n_mixed", "n_div", + "n_shift_step"]) +def test_ineligible_reductions_stay_sequential(rbuilt, name): + f = fn(rbuilt, name) + assert f.parallel is True + loops = par_loops(f) + assert loops and all(l.par is None for l in loops), name + + +def test_unmarked_reduction_never_parallel(rbuilt): + assert all(l.par is None for l in par_loops(fn(rbuilt, "r_unmarked"))) + + +def test_c_regions_only_for_eligible_reductions(rbuilt): + import re + owners = set(re.findall(r"static void tp_par_(\w+?)(?:_fast)?_\d+_body\(", rbuilt[1])) + assert owners == set(RED), owners + + +def _forge(rbuilt, name, par): + f = fn(rbuilt, name) + loop = par_loops(f)[-1] + fake = dataclasses.replace(loop, par=par) + bad = dataclasses.replace(f, body=tuple(fake if s is loop else s for s in f.body)) + return verify.verify(dataclasses.replace(rbuilt[3], functions=(bad,), skipped={})) + + +def test_verifier_rejects_a_forged_reduction(rbuilt): + good = par_loops(fn(rbuilt, "r_sum"))[-1].par + # forged on an ineligible loop (the accumulator is read in the body) + for name in ("n_acc_read", "n_two_ops", "n_copy", "n_mixed", "n_div"): + proved, diags = _forge(rbuilt, name, ir.ParInfo((), (), 3, (("s", "ADD"),))) + assert not proved.functions and any(d.rule == PAR_RULE for d in diags), name + # right loop, wrong reduction operator + proved, diags = _forge(rbuilt, "r_sum", dataclasses.replace(good, reduce=(("s", "MUL"),))) + assert not proved.functions and any(d.rule == PAR_RULE for d in diags) + # right loop, reduction dropped + proved, diags = _forge(rbuilt, "r_sum", dataclasses.replace(good, reduce=())) + assert not proved.functions and any(d.rule == PAR_RULE for d in diags) + # the honest value passes + proved, diags = _forge(rbuilt, "r_sum", good) + assert proved.functions and not diags + + +@pytest.mark.parametrize("n", RSIZES) +def test_reductions_match_interpreter_no_overflow(rmod, n): + xs, ys = rdata(n, 1), rdata(n, 2) + for name in ("r_sum", "r_range_sum"): + same(rmod, name, xs) + same(rmod, "r_sum_init", xs, 12345) + same(rmod, "r_sum_mul", xs, 977) + same(rmod, "r_two", xs) + same(rmod, "r_prod", [1, -1] * (n // 2) + [1] * (n % 2), 1) # +-1: never overflows + + +@pytest.mark.parametrize("n", RSIZES) +def test_sums_with_large_values_that_fit(rmod, n): + xs = rdata(n, 3, -(2**62) // max(n, 1), 2**62 // max(n, 1)) # partial sums stay inside i64 + same(rmod, "r_sum", xs) + same(rmod, "r_sum_mul", rdata(n, 4, -5, 5), 2**40) + + +def test_product_with_zero_and_small_values(rmod): + n = 300000 + xs = [1] * n + xs[10], xs[200000] = 7, -3 + same(rmod, "r_prod", xs, 1) + xs[150000] = 0 + same(rmod, "r_prod", xs, 1) + same(rmod, "r_prod", [2] * 62 + [1] * (n - 62), 1) # 2**62 fits + same(rmod, "r_prod", [-2] * 63 + [1] * (n - 63), 1) # -2**63 fits exactly + + +def test_local_array_reduction(rmod): + for n in (0, 1, 5, 300000): + same(rmod, "r_local", n) + + +# overflow partway: the whole region is discarded and the sequential rerun produces Python's int + +@pytest.mark.skipif((os.cpu_count() or 1) < 2, reason="needs a second core to start a region") +def test_sum_overflow_in_a_chunk_reruns_and_equals_python(rmod): + n = 400000 + base = rdata(n, 6, -3, 3) + xs = list(base) + xs[0] = xs[1] = 2**62 # chunk 0 overflows, and the prefix + xs[n - 1] = -(2**62) + r0, rr0 = regions(rmod), reruns(rmod) + got = same(rmod, "r_sum", xs) + assert got == sum(xs) + assert regions(rmod) == r0 + 1 and reruns(rmod) == rr0 + 1 + # not representable in i64 at the end: Python's big int + xs2 = list(base) + xs2[0] = xs2[n - 1] = I64_MAX + got = same(rmod, "r_sum", xs2) + assert got == sum(xs2) and got > I64_MAX + + +@pytest.mark.skipif((os.cpu_count() or 1) < 2, reason="needs a second core to start a region") +def test_sum_overflow_only_in_the_combination_reruns(rmod): + n = 400000 + xs = [0] * n + xs[0] = 2**62 # first chunk partial: 2**62 + xs[n - 1] = 2**62 # last chunk partial: 2**62; total 2**63 + rr0 = reruns(rmod) + got = same(rmod, "r_sum", xs) + assert got == 2**63 + assert reruns(rmod) == rr0 + 1 + + +@pytest.mark.skipif((os.cpu_count() or 1) < 2, reason="needs a second core to start a region") +def test_sum_that_overflows_and_comes_back_equals_python(rmod): + n = 400000 + xs = [0] * n + xs[1], xs[2], xs[3] = 2**62, 2**62, -(2**62) # prefix passes 2**63, ends at 2**62 + got = same(rmod, "r_sum", xs) + assert got == 2**62 + xs = rdata(n, 7, -5, 5) + xs[100], xs[101], xs[n - 5] = I64_MAX, I64_MAX, I64_MIN + got = same(rmod, "r_sum", xs) + assert got == sum(xs) + + +@pytest.mark.skipif((os.cpu_count() or 1) < 2, reason="needs a second core to start a region") +def test_product_overflow_reruns_and_equals_python(rmod): + n = 400000 + xs = [1] * n + for k in (3, 150000, 399999): + xs[k] = 2**22 + rr0 = reruns(rmod) + got = same(rmod, "r_prod", xs, 1) # 2**66 + assert got == 2**66 and reruns(rmod) == rr0 + 1 + got = same(rmod, "r_prod", [3] * 60 + [1] * (n - 60), 1) # 3**60 > 2**63 + assert got == 3**60 + xs = [1] * n + xs[0] = xs[n - 1] = 2**40 + xs[7] = 0 + got = same(rmod, "r_prod", xs, 1) # intermediate overflow, result 0 + assert got == 0 + + +@pytest.mark.skipif((os.cpu_count() or 1) < 2, reason="needs a second core to start a region") +def test_element_expression_overflow_reruns(rmod): + n = 400000 + xs = rdata(n, 8, -9, 9) + xs[250000] = 2**30 + rr0 = reruns(rmod) + got = same(rmod, "r_sum_mul", xs, 2**40) # xs[i] * c overflows at one element + assert reruns(rmod) == rr0 + 1 + assert got == __import__("functools").reduce(lambda a, b: a + (b * 2**40), xs, 0) + + +@pytest.mark.skipif((os.cpu_count() or 1) < 2, reason="needs a second core to start a region") +def test_region_and_rerun_counters(rmod): + r0, rr0 = regions(rmod), reruns(rmod) + rmod.r_sum([1] * 10) + assert regions(rmod) == r0 # below threshold: sequential + assert rmod.r_sum([1] * 1048576) == 1048576 + assert regions(rmod) == r0 + 1 and reruns(rmod) == rr0 # above: a region, no rerun + r1 = regions(rmod) + for name in ("r_unmarked",): + getattr(rmod, name)([1] * 1048576) + rmod.n_acc_read([0] * 1048576) + rmod.n_two_ops([1] * 1048576) + rmod.n_float([1.0] * 1048576) + rmod.n_copy([1] * 1048576) + rmod.n_mixed([1] * 1048576) + assert regions(rmod) == r1 + + +def test_negative_cases_correct(rmod): + xs = rdata(300000, 9, -5, 5) + for name, args in (("n_acc_read", (xs,)), ("n_two_ops", (xs,)), ("n_float", ([float(v) for v in xs],)), + ("n_copy", (xs,)), ("n_mixed", (xs,)), ("n_div", (xs, 3)), + ("n_shift_step", (xs,)), ("r_unmarked", (xs,))): + want = rmod.__typedpython_interpreted__[name](*[list(a) if isinstance(a, list) else a for a in args]) + got = getattr(rmod, name)(*[list(a) if isinstance(a, list) else a for a in args]) + assert got == want, name + + +RTHREAD_BODY = """ + import random + r = random.Random(11) + n = 400001 + small = [r.randint(-1000, 1000) for _ in range(n)] + ovf = list(small); ovf[0] = ovf[1] = 2**62; ovf[n - 1] = -(2**62) + ones = [1] * n; ones[5] = ones[200000] = ones[n - 1] = 2**22 + out = { + "sum": m.r_sum(small), "ovf": m.r_sum(ovf), "prod": m.r_prod(ones, 1), "two": m.r_two(small), + "mul": m.r_sum_mul(small, 2**40), + "regions": m.__typedpython_parallel__()[0], "reruns": m.__typedpython_parallel__()[3], + } + print(json.dumps(out)) +""" + + +def test_reduction_results_identical_for_every_thread_count(rbuilt): + import random + import functools + r = random.Random(11) + n = 400001 + small = [r.randint(-1000, 1000) for _ in range(n)] + ovf = list(small) + ovf[0] = ovf[1] = 2**62 + ovf[n - 1] = -(2**62) + ones = [1] * n + ones[5] = ones[200000] = ones[n - 1] = 2**22 + want = {"sum": sum(small), "ovf": sum(ovf), "prod": 2**66, "mul": sum(small) * 2**40, + "two": sum(small) * 3 + 7 + sum(small) * 2} + results = {t: run_child(rbuilt[2], RTHREAD_BODY, {"TYPEDPYTHON_THREADS": str(t)}) for t in (1, 2, 3, 7)} + for t, res in results.items(): + for k, v in want.items(): + assert res[k] == v, (t, k) + assert results[1]["regions"] == 0 and results[1]["reruns"] == 0 # 1 disables the feature + if (os.cpu_count() or 1) > 1: + assert results[2]["regions"] > 0 and results[7]["regions"] > 0 + assert results[2]["reruns"] == 2 and results[7]["reruns"] == 2 # ovf and prod rerun, once each + + +def test_no_reduction_region_without_a_marker(tmp_path): + src = RSOURCE.replace("# typedpython: parallel\n", "") + _, c_text, _, _ = build(tmp_path, src, name="rnomark") + assert "tp_par_run(" not in c_text and "__typedpython_parallel__" not in c_text