Unified threading - #44
Open
hakkelt wants to merge 48 commits into
Open
Conversation
…t merge Adds AcceleratedDCTs test dependency, disables Aqua's persistent_tasks check, and carries forward the test coverage expansion, test refactors, and small dead-code/API cleanups developed on gpu-refactor after the GPU-support work was squash-merged upstream as 861af22.
It requires Julia >= 1.11 and is already added dynamically at runtime for that case in test/utils.jl; listing it statically forced dependency resolution to fail on the Julia 1.10 LTS CI job.
The private thread_count_functions registry and the @restrict_threading / @enable_full_threading macros in src/utils.jl are replaced by the new standalone NestedThreading.jl package. Beyond relocating the mechanism this fixes a latent correctness bug: the old save/restore was per scope, so two concurrent batched mul! calls captured each other's already-restricted thread counts and left BLAS permanently pinned at 1 with no scope active (reproducible before this commit, fixed after). NestedThreading refcounts budgets process-wide, snapshotting once and restoring once, and applies the minimum over all active scopes so nesting can only narrow a budget, never widen it. - src/utils.jl: delete the threading machinery. - SpreadingBatchOp/SimpleBatchOp: @restrict_threading @threads -> the @budgeted_threads macro; the LOCKING @sync/@Spawn sites use the function form with_restricted_threads. - FFTWOperators: drop the __init__ that reached into AbstractOperators' private registry; the FFTW pool now registers via NestedThreading's own package extension. - NFFTOperators: delete set_nfft_threading_expr and its two macros in favour of a with_nfft_threading helper over NestedThreading's scoped API, at all six call sites. - Inner libraries now get nthreads() / trip_count threads instead of always 1, so a 2-element batch on 8 threads runs BLAS at 4. - Add a regression test for the concurrency bug and for budget nesting. NestedThreading is wired in via [sources] pending registration in General, so these changes are not ready to push yet. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
BroadCast.jl's threaded adjoint mul! runs an outer @threads loop whose body calls mul! on arbitrary sub-operators, which may themselves use BLAS/FFTW. It never went through the old @restrict_threading, so it was a genuine oversubscription source that the migration had left in place. This was the last raw @threads in src/, so @threads is dropped from the imports; @Spawn and nthreads are still used. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
`TestUtils` is the setup module of nearly every testitem, and it unconditionally ran `GPUEnv.activate(; persist=true)` plus the FFTW-specific `AcceleratedDCTs` install. Any filtered run paid full GPU environment activation even when it selected no GPU test. Move both into a new `GpuEnvSetup` @testmodule (test/gpu_env_setup.jl). TestItemRunner evaluates setup modules lazily, only for testitems that survive the filter, so activation now happens iff a `:gpu` testitem is selected. All 47 `:gpu` testitems gain `GpuEnvSetup` in their setup list. `runtests.jl` publishes `ABSTRACTOPERATORS_TEST_GPU` from the filter closure, set from the tags of the items it actually accepts. Deriving it from the filter *string* would be wrong: GPU testitems also carry their category tags, so `:linearoperator` selects `Eye (GPU)` too. `GpuEnvSetup` errors when the flag is false, which can only mean a testitem missing its `:gpu` tag pulled the module in. Audit results: - `:gpu` testitems keep their inline `using GPUEnv`; the import is a name binding and costs nothing, only `activate` was expensive. - test_gpu_quality.jl's two items are tagged `:gpu` but exercise GpuExt through JLArrays alone, a direct test dep, so they deliberately do not use `GpuEnvSetup`. Measured: single non-GPU testitem 3m20s -> 8.7s. Pass counts unchanged -- `:jet` 224/224 and `:misc,:quality` 185/185 with zero GPU activation, `:calculus` 1557/1557 (still activates, correctly: it selects the `*(GPU)` items by their `:calculus` tag), `"Eye (GPU)"` 27/27 through the new setup module. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
Phases 1-4 of THREADING_PLAN.md. Threading was previously expressed three different ways with four scattered thresholds, and `copy_operator`'s deepcopy fallback silently ignored both of its keywords. Infrastructure (src/threading_policy.jl): - `is_threaded` trait, `supports_threading`, `default_threaded`/`threading_threshold` per-operator policy, and `adapt_operator`. - `copy_operator` now always copies; the "share when thread-safe" short-circuit moved to `adapt_operator`. The deepcopy fallback refuses `storage_type`/`threaded` instead of ignoring them, naming the type that needs a `_copy_operator_impl`. - Thresholds are transcriptions of benchmark/threading_sweep.jl (new), not guesses; each carries a PROVENANCE line. The sweep contradicted the plan's hypothesis that FastBroadcast wins or ties everywhere: for a pure copy `@.. thread=true` is time-identical to serial at every swept size, while `@batch` crosses over at 2^15. Nesting safety (the bug the plan called out): - create_BatchOp, create_threaded_SpreadingBatchOp and OperatorBroadCast adapted only the per-thread *copies*, leaving instance 1 threaded. All three now go through `_per_thread_operators`, which also branches on `is_thread_safe` to decide share vs copy -- `require_thread_safe` cannot manufacture thread safety, so it must not pretend to. - `opType` had to move after the adaptation: threading is a type parameter. Elementwise threading: all 10 nonlinear operators (SoftMax excluded, with the reduction reason recorded), FiniteDiff, and forwarding traits for the calculus operators. FiniteDiff `mul!` no longer allocates: `b[idx_1] .- b[idx_2]` materialised two temporaries per call. With `@views` it is allocation-free and 45x faster at n=1e6 (2581us -> 57us). Two problems the tests surfaced, both fixed in source rather than in the assertion: - `adapt_operator(op; threaded=false)` silently no-opped on forwarders, which lacked `is_threaded` and so inherited the `false` default while their children stayed threaded. - JET flagged runtime dispatch in `FiniteDiff(T, dims)`, which fell through to the `dir::Int` method and its `Val(dir)`. Split out a D=1 literal method. The two batching perf tests needed wider margins, not weaker ones: the FiniteDiff allocation fix sped the serial baseline up ~3x, shrinking an honest 2.5x speedup to 1.7x and letting runner noise flip a bare `<`. Both now take a minimum across repetitions at a workload size measured to favour batch parallelism (smaller-items/more-items was measured to hurt). Tests: :batching 283, :jet 245 (incl. new @test_opt/@test_call threading coverage), :quality 61, :calculus+:linearoperator+:nonlinearoperator+:misc 2698, :Threading+:Syntax 276 -- all passing. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
Phase 5 of THREADING_PLAN.md, for DCAT: its blocks are fully independent in both directions, with no accumulation and no shared output. The policy needed a second condition the plan did not anticipate. Aggregate work alone does not predict a win: 2 blocks x 2^18 elements is a 2^19 aggregate but measured only 1.03x, because two blocks cap the gain at 2x and the threading overhead eats it, while the same 2^19 aggregate across 8 blocks measured 1.9x. So `default_block_threaded` requires both `THRESHOLD_BLOCK_PARALLEL` (2^18, measured) and `MIN_BLOCKS_FOR_PARALLEL` (4). That pair admits every measured win and excludes every measured loss: nb=8 bs=2^16 1.90x threaded nb=8 bs=2^12 0.26x not threaded nb=4 bs=2^16 1.18x threaded nb=2 bs=2^18 1.03x not threaded nb=8 bs=2^18 9.67x threaded The threaded `mul!` is deliberately not `@generated`: a generated body must be pure and the threading macros expand to closures. The parallel loop is therefore a plain function and only the single-block body is generated, selected by `Val(i)`. That costs one dynamic dispatch per block -- nanoseconds against a whole child `mul!` -- and keeps the per-block type stability a closure over a heterogeneous tuple would have lost. `@budgeted_threads` rather than `@batch` because each body may itself reach BLAS/FFTW. Nesting safety applies here too: when the block loop threads, the constructor adapts the blocks to `threaded = false`. VCAT-forward and HCAT-adjoint from the plan are not included; they need their own measurement before shipping, per the plan's own rule that unmeasured variants are dropped. Tests: :Threading,:DCAT,:calculus 1769; :jet,:quality,:misc 462 -- all passing. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
…ression
Phase 6 of THREADING_PLAN.md, plus a regression this phase uncovered.
The regression: Phase 1's strict `copy_operator` fallback made batching *any* FFTW/DSP
operator raise. Those operators own scratch buffers, so a threaded batch needs one private
copy per thread; but they have no `_copy_operator_impl`, and the fallback refused the
`threaded = false` request the batch layer makes. The fix belongs in the fallback's rule,
not in per-operator patches: an operator with no threaded path cannot have a `threaded`
request violated, so deepcopy is a correct answer for it. That is the same rule
`_satisfies_constraints` already used, now applied consistently. `storage_type` is still
refused, since deepcopy genuinely cannot honour it.
Per subpackage:
- NFFTOp already carried a `threaded` field; it gains the traits and a `_copy_operator_impl`
that shares the immutable plan and dcf while giving the copy its own scratch buffer.
Requests that would need a different plan (thread count, storage backend) raise rather
than being silently ignored.
- DFT/IDFT gain `threaded` as the package-wide spelling of the existing `num_threads`, and
record the plan-time count so `is_threaded` can report it. FFTW is a counted pool, so this
is fixed at construction; switching it replans.
- DSP and Wavelet operators declare `supports_threading = false` explicitly rather than
inheriting the default, which is what lets a threaded batch wrap them.
Two bugs caught while writing this, both mine:
- A "cleanup" of a trailing `'` deleted the adjoint postfix from five IDFT constructors,
turning IDFT into a plain DFT and silently disabling the DFT/IDFT combination rules. The
FFTW combination tests caught it; there is now a test asserting IDFT stays adjoint-wrapped.
- The DFT copy replanned from the codomain element type, which is wrong for a real-input
DFT (`C == Complex{D}`) and produced an operator with a complex domain.
DCT/RDFT/IRDFT are not wired for plan-time thread counts -- they have no `num_threads`
plumbing at all today -- so they are declared unthreaded rather than half-wired. Recorded in
the plan's "Not done" section along with VCAT/HCAT block parallelism and the four legacy
thresholds that still gate the pre-existing operators.
THREADING_PLAN.md gains an execution record: measured-vs-guessed thresholds, the three plan
hypotheses that measurement contradicted, and what was deliberately left out.
Tests: full suite 3856 passed, 0 failed, 0 errored (21m47s).
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
…all operators
Follow-up round covering VCAT/HCAT block threading, the remaining FFTW transforms, the
legacy thresholds, per-operator thresholds, and one semantics change that ties them
together.
`threaded` is now `Bool` only, with one meaning everywhere, resolved in one place
(`_resolve_threaded`):
false veto -- never threads. The only hard directive, and it must be: nesting
safety depends on a block/batch loop switching its children off and
being able to rely on it.
true (default) permission -- threads only if the policy also agrees.
Three families previously disagreed. Elementwise and block operators treated `true` as an
absolute override, batch operators as a permission, and the FFTW transforms bypassed the
policy entirely -- so a 256-point threaded DFT reported is_threaded == true while measuring
0.02x. `is_threaded` now reports what the operator will actually do, never what was asked
for; consequently a single-threaded session or GPU storage reports false without any
special-casing. `copy_operator`/`adapt_operator` keep accepting `nothing`, since there
`threaded` is a *constraint* ("no constraint, preserve") rather than a setting -- without it
a plain copy could not preserve an explicitly serial operator.
VCAT-forward and HCAT-adjoint now thread their block loops, and measuring them forced two
corrections to the Phase 5 policy:
- per-block, not aggregate: at a fixed 2^18 aggregate VCAT measures 1.32x with 4 blocks of
2^16 but 0.85x with 16 blocks of 2^14;
- per-operator: at 4 blocks of 2^16 VCAT measures 1.32x while HCAT measures 0.76x, because
HCAT's adjoint carries per-block ArrayPartition indexing. One shared constant would have
shipped that regression.
FFTW does thread r2r and r2c (measured 3.08x r2c, 2.25x DCT, 1.91x IDCT at n=2^22), so
DCT/IDCT/RDFT/IRDFT gained plan-time thread counts, and all six transforms are now gated on
measured crossovers (c2c 2^13, r2r 2^15, r2c 2^15). `num_threads` stays an explicit command
rather than a permission -- the escape hatch for callers who know their workload.
The four legacy thresholds are gone, which fixed two live defects:
- Variation's two constructors thresholded in different *units* (bytes vs elements), so the
same 10000-element input produced opposite threading;
- `_should_thread(::AbstractOperator)` had no size component at all, so batch operators
threaded four-element work.
Scale's legacy 1e4 cutoff was ~400x below its measured 2^22 crossover.
Thresholds are now per-operator transcriptions from benchmark/operator_thresholds.jl (new),
which measures each operator's real `mul!`. The cost classes proved too coarse: within
"transcendental" the crossovers span 2^8 (SoftPlus) to 2^11 (integer Pow), and Pow splits
internally by exponent kind since `x^0.5` lowers to `exp(p*log(x))`.
Also unified: SignAlternation had its own `threaded && nthreads() > 1` with no size policy
and no `is_threaded` method, so its trait always reported false regardless of what it did.
JET caught a regression that no functional test could: reassigning the captured `A` inside
the DCAT/VCAT/HCAT constructors boxed it to `Any`, cascading into 11 runtime-dispatch sites
in a constructor that should infer fully. Results were correct throughout; only @test_opt
saw it.
Tests: full suite 3971 passed, 0 failed. Three batching tests needed updating because they
encoded the old "threaded=true forces threading" assumption -- one of them was reading
`batch_op.operator[1]`, which on the single-threaded struct silently *slices* the wrapped
operator instead of erroring; it now goes through the `_wrapped_operator` accessor.
Known, unfixed, pre-existing: Variation's adjoint throws BoundsError when the trailing
dimension is exactly 2. Present on master; recorded in THREADING_PLAN.md.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
The adjoint's per-dimension branch chain tested `i == 2` before `i == size(y, d)`, so the last column of a dimension of length exactly 2 took the interior formula and read one element past the end. This hit *any* dimension of size 2, not only the trailing one -- (2,8) and (8,2,3) failed alongside (8,2) -- and the result was mathematically wrong there, not merely out of bounds: column n of the transpose is `y_1 + y_n` when n == 2, while the interior formula gives `y_1 + y_2 - y_3`. Rewrite the contribution as three independent terms (`_variation_adjoint_term`) rather than mutually exclusive branches, so the n == 2 case correctly takes both the `i == 2` and `i == n` terms. Both mul! methods now share one body so the threaded and serial paths cannot drift apart. `N` is passed to the shared body as a plain Int, deliberately not a Val: inside Polyester's `@batch` the loop body becomes a closure that does not carry the enclosing method's static parameters, so a `Val(N)` argument is built from a runtime value and every element dispatches dynamically. That allocated 293 MB at n = 2^22 and turned a 4.5x speedup into 0.25x, while leaving every correctness test green -- hence the new allocation testitem, timings being too noisy to assert on. Also reject singleton dimensions at construction. They have no finite difference to take and already failed with a BoundsError inside the *forward* kernel, so an ArgumentError naming the cause is strictly better; the array-based constructor now delegates to the tuple one so the guard and the threading policy each exist in one place. Verified against the dense forward matrix's exact transpose over 16 shapes x 2 element types x 2 threading modes. Threaded adjoint is back to baseline (12527 us vs 12668 us at n = 2^22, zero allocations) and the serial path is ~1.75x faster at mid sizes (3528 vs 6195 us at 2^18), the three-term form replacing unpredictable branches with straight-line arithmetic. Measured crossover is unchanged at 2^10. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
The plan is fully executed and its content is now recorded where it is actually used: threshold provenance in `threading_threshold` docstrings, policy rationale in `src/threading_policy.jl`, and the measurement methodology in `benchmark/operator_thresholds.jl`. Nothing references the file. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
Threaded coverage went from 15 entries to 45. Every operator whose policy
can grant threading now has a matched `-single` / `-threaded` pair:
FiniteDiff, Scale, HCAT/VCAT/DCAT, DFT, NFFTOp and the nine threading-
capable non-linear operators join the five that were already covered.
Three defects surfaced while doing it:
1. FiniteDiff and all nine non-linear operators were *implicitly*
threaded. Their constructors default to `threaded = true` and every
benchmark size clears the operator's threshold, so a single entry
measured the threaded path on a dev box and the serial path on a
one-core runner under the same name. Base entries are now pinned
`threaded = false` and renamed to `-single`, matching the convention
DiagOp/Variation/BroadCast already used.
2. `threaded = true` is a permission, not a command, so a mis-sized
benchmark silently measures the serial path -- and would keep
reporting a healthy "threaded" number straight through a threading
regression. Every threaded entry now routes through `check_threaded`
(or `check_block_threaded` for the *CATs, where `is_threaded` is also
true when merely the blocks thread). Sizes that the policy rejects are
given their own constants: Scale needs 2^22, and the *CATs need four
blocks rather than the two the serial states use.
3. At `evals = 1` the threaded entries were dominated by first-touch page
faults on the several MB `setup` allocates per sample. DiagOp's
threaded forward measured a median 4.15x its own minimum and came out
*slower* than the serial path it beats by 8x once warm. Threaded
entries now use `evals = 50`, which drops the spread to 1.5x and makes
the minimum agree with a standalone hot-loop measurement. Serial
entries keep `evals = 1` so their baselines stay comparable.
Multithreaded entries are skipped whenever `CI` is set, and when Julia has
one thread. `ABSTRACTOPERATORS_BENCH_THREADED=true|false` overrides both,
for a dedicated runner with a pinned core count. On CI the suite is 87
deterministic entries; locally with threads it is 134.
benchmarks.jl is now an entry point over `bench_common.jl` (probes, gate,
constants, shared fixtures) plus one file per suite in `suites/`, each
runnable on its own:
julia --project=benchmark benchmark/suites/linearoperators.jl
SUITE keys and grouping are unchanged by the split: 134 entries across 9
groups, and the per-suite counts sum to 134.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
Co-Authored-By: Claude Haiku 4.5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
hakkelt
force-pushed
the
unified-threading
branch
from
August 19, 2026 15:01
0287558 to
a8a806a
Compare
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #44 +/- ##
==========================================
- Coverage 90.13% 89.83% -0.30%
==========================================
Files 51 52 +1
Lines 3678 4111 +433
==========================================
+ Hits 3315 3693 +378
- Misses 363 418 +55 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Benchmark Results (Julia v1.12.7)Run once per thread count ( 1 thread🚀 6 benchmarks improved in time · 🚀 2 benchmarks use less memory Time benchmarks
Memory benchmarks
2 threads🚀 25 benchmarks improved in time · 🚀 3 benchmarks use less memory · 🐢 3 memory regressions detected Time benchmarks
Memory benchmarks
|
_fbthread(b::Bool) dispatched through Val(b) across two singleton methods, which Julia LTS's inference could not resolve at compile time (unlike current stable, which constant-folds the whole chain) -- JET flagged it as runtime dispatch in DiagOp/Scale @test_opt. Replaced with a direct ternary, which is a single method returning a small Union instead of a multi-method dispatch. Documenter's checkdocs=:exports also failed: is_threaded, adapt_operator, and supports_threading are exported but had no @docs entry. Added a Threading section to docs/src/properties.md. That exposed two @ref links (_resolve_threaded, used twice) pointing at an internal, non-exported helper with no rendered docstring; unlinked them to plain code spans since they're implementation details, not public API. Verified locally: JET/threading test suite (787 tests) and `docs/make.jl` both pass.
…-loss bug Codecov flagged this PR's patch coverage at 68.65% (211 lines missing). Most of the gap was the new `is_threaded`/`supports_threading`/`_children` one-liners and `_copy_operator_impl` paths added across AffineAdd, Ax_mul_Bx, Ax_mul_Bxt, Axt_mul_Bx, HadamardProd, BroadCast, SimpleBatchOp, and SpreadingBatchOp, none of which any existing test called directly -- plus HCAT's threaded adjoint block loop, whose only exercising test used a block size below HCAT's own (higher) threshold and so never actually took the threaded path. Writing the DiagOp/SpreadingBatchOp coverage surfaced a real bug: DiagOp's `_copy_operator_impl` ran `op.d` (the operator's actual diagonal data, not a scratch buffer) through `_convert_buffer`, which is `similar`-based and intentionally does not copy values -- correct for `mul!` scratch buffers that get overwritten every call, silently wrong for the one field that isn't. Any `copy_operator(diag_op; storage_type = ...)` was replacing the diagonal with uninitialized memory, and a scalar `d` MethodError'd outright since `_convert_buffer` has no `Number` method. Fixed with a `copyto!`-based `_copy_diag` mirroring AffineAdd's existing `_copy_displacement`, and strengthened the existing DiagOp copy_operator test (which asserted `!==` but never checked the copied values were still correct) to catch this class of regression going forward. Also fixed two test-writing bugs surfaced while adding this coverage: the new SpreadingBatchOp/HCAT threading tests initially used per-item/per-block sizes below the relevant threshold, so `threaded = true` was correctly declined by the policy and the assertions were exercising the serial path while believing they were exercising the threaded one. Verified locally (JULIA_NUM_THREADS=2, matching CI): full :batching/:Threading tag run (1375 tests) and the full :jet suite (245 tests), both passing. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
…dingBatchOp gap Two more codecov-flagged gaps from the threading migration: - Scale and AdjointOperator gained the same is_threaded/supports_threading/ _children/_copy_operator_impl additions as the other operators, with zero test coverage. Scale also needed its ArrayPartition-codomain branch exercised (codomain_type_for_policy/codomain_array_type_for_policy), covered here via a Scale wrapping a DCAT. - Bigger issue: nearly every existing SpreadingBatchOp test that names a specific threading strategy (COPYING, LOCKING, FIXED_OPERATOR, AUTO) used a per-item domain size of 5-15 elements. This PR's threading policy gates `threaded = true` on MIN_BATCH_WORK_FOR_PARALLEL (2^10 elements) before it ever reaches strategy selection, so every one of those tests was silently falling back to the single-threaded branch regardless of which strategy it named -- `test_failing_nonthreadsafe_spreading_batch_op` already carried a comment explaining this exact trap for its own workload. Bumped the shared `test_nonthreadsafe_spreading_batch_op` helper and the four named-strategy testitems to sizes that actually cross the gate, and added `isa` assertions so a future regression back to the wrong strategy fails loudly instead of silently. One power-iteration opnorm assertion needed its tolerance loosened to match the larger workload (same relaxation already used one line below it for the same estimate). Verified locally (JULIA_NUM_THREADS=2): 1082 tests passing across :batching/:Threading/:Scale/:AdjointOperator/:jet tags. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
…nd BroadCast's threaded construction
Third round closing the codecov gap on the threading-migration diff, plus a
cleanup found along the way:
- OperatorWrapper: same is_threaded/supports_threading gap as the other
operators covered in the previous two commits.
- HCAT: every mul!/adjoint mul!/_hcat_block_adj! branch that handles a
multi-element ("tuple-idxs") sub-operator was untested -- all existing
coverage used single-index HCATs, so the `Pi <: Integer` false branch never
ran anywhere, in either the natural or the permuted (indexed) form, serial
or threaded. Added one operator built from a HCAT-of-DiagOp wrapped in
Compose(FiniteDiff, ...) to keep it non-diagonal (diagonal wrappers get
simplified straight back into a flat HCAT by the package's own combination
rules) and cheap enough to push past the block-threading size threshold.
Verified with `test_op`'s forward/adjoint dot-product identity rather than
hand-derived expected values, which is what caught two dimension bugs in
earlier drafts of this test before they were committed.
- BroadCast: OperatorBroadCast's threaded-construction branch (per-thread
domain buffers, per-thread operator copies) needs a broadcast sized above
THRESHOLD_MEMORY_BOUND; the existing "non-compact threaded" test used ~60
elements, so `threaded = true` was always declined by the policy.
- Removed a stray `>>>>>>> d84c552 (...)` conflict marker left at the end of
test/calculus/test_hcat.jl from the branch's original rebase, before any of
this session's work. Harmless to test execution (outside any @testitem) but
should not have been committed.
Verified locally (JULIA_NUM_THREADS=2): instrumented coverage now shows zero
missed executable lines in HCAT.jl for the :HCAT tag scope (was 8), and
BroadCast.jl's only remaining misses are two pre-existing lines outside this
PR's diff. Full :HCAT/:BroadCast/:OperatorWrapper/:jet run: 618 tests passing.
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
Same is_threaded/supports_threading/_children gap as the operators covered in the previous three commits: SimpleBatchOp.jl and VCAT's threaded-forward mul! path were already fully exercised (SimpleBatchOp by existing tests, VCAT by the "VCAT-forward and HCAT-adjoint block loops" threading-contract test at the right block size) but neither test asserted is_threaded/supports_threading themselves, and Sum had no threading-trait coverage at all. Verified locally (JULIA_NUM_THREADS=2): instrumented coverage now shows zero missed executable lines in SimpleBatchOp.jl and VCAT.jl for the :batching/ :Threading tag scope. Sum.jl's remaining local misses are pre-existing @generated-function lines outside this PR's diff. Full :Sum/:VCAT/:Threading run: 496 tests passing. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
Both gained a _copy_operator_impl in this PR (the whole reason for it is honouring storage_type, since neither has a threaded path) with no test exercising it. Verified locally: instrumented coverage now shows zero missed executable lines in both files. Full :LMatrixOp/:MyLinOp/:jet run: 353 tests passing. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
A review comment on this branch claimed "None of the DSP operators has a Julia-level threaded path: the work happens inside DSP.jl/FFTW, which manage their own parallelism." That's wrong on two counts: DSPOperators does not depend on DSP.jl at all (only AbstractFFTs/FFTW), and Conv/Xcorr plan their own FFTW transforms directly -- they were simply never wired into FFTW's thread count the way FFTWOperators' DFT/DCT/RDFT already are, so they were hard-coded to `supports_threading = false` alongside Filt/MIMOFilt, which genuinely have no threaded path (IIR filtering is a sequential recursion, not an FFT). Gave Conv and Xcorr the same `num_threads`/`threaded` constructor keywords, `is_threaded`/`supports_threading`, and `_copy_operator_impl` (replan on a threaded-flag change, share plans and just re-copy scratch buffers otherwise) as FFTWOperators' DFT, via a small `_dsp_fftw_num_threads`/ `_dsp_with_fftw_threads` pair mirroring `FFTWOperators._fftw_num_threads` (DSPOperators does not depend on FFTWOperators, so this isn't reused directly; the threshold is marked provisional, borrowed from FFTWOperators' measured :c2c class since both plan a c2c/r2c FFT of the padded convolution length). Two bugs surfaced while getting this working: - The first draft wrapped each constructor's *whole* plan-building block (including buffer allocation and `if domain_type <: Real` branching) in the `_dsp_with_fftw_threads` closure. That runs correctly but JET's `@test_call` couldn't infer the closure's return type through the branching and reported spurious "local variable not defined" errors. Every existing FFTW-based operator in this codebase (RDFT, DCT) only wraps the bare `plan_*` calls themselves -- matching that convention fixed the JET report, not just the symptom. - Xcorr's `_copy_operator_impl` replan path called the `(dim_in, h)` constructor Xcorr never actually defines (unlike Conv, which does); fixed to call the `(domain_type, dim_in, h)` form it does have. Also updated the pre-existing "subpackage operators declare their threading" and "batching an operator that has no threaded path" tests, which had used Conv/Xcorr as their canonical no-threaded-path example -- replaced with Filt, which is the operator that actually fits that description. Verified locally (JULIA_NUM_THREADS=2): full :dsp/:batching/:Threading/:jet/ :quality run (1057 tests) and a full Documenter doctest build, both passing. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
…eturn-before-loop _copy_operator_impl's purpose is to produce a new operator satisfying either a storage_type or threaded request (or both) -- refusing storage_type outright, which the previous commit did by copying FFTWOperators' DFT.jl pattern verbatim, defeats that. Unlike DFT (whose plans are the only per-instance state), Conv/Xcorr can genuinely support it: `h` is copied to the requested array type with `copyto!` (it's the operator's actual filter data, not a scratch buffer, so uninitialized allocation would silently corrupt it -- the same class of bug fixed for DiagOp earlier in this branch) and the operator is rebuilt through the public constructor, which replans FFTW against whatever backend the copied `h` now lives on. Also fixed a Runic formatting artifact: `_xcorr_fir_adj!`'s tiled FIR adjoint loop ended in `return @inbounds while ... end`, wrapping a loop (whose value is always `nothing`) in `return` rather than following it with a bare `return`. Grepped the full branch diff against upstream master for the same `return for`/`return while` pattern; this was the only occurrence. Verified locally (JULIA_NUM_THREADS=2): storage_type round-trips correctly for both operators (copied array is a distinct object with equal values, `copy * x` matches the original), combined with threaded=true, and the full :dsp/:jet/:Threading run (433 tests) passes. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
…sed it Same issue as the previous Conv/Xcorr fix, found across the rest of the branch: DFT, RDFT, IRDFT, DCT/IDCT (FFTWOperators) and NFFTOp (NFFTOperators) all threw ArgumentError on a storage_type request instead of honouring it. _copy_operator_impl's contract is to produce a new operator satisfying either kind of request; refusing one outright isn't a valid implementation of that contract regardless of how the operator is built internally. FFTWOperators (DFT/RDFT/IRDFT/DCT/IDCT): none of these hold persistent input data, only plans (and scratch, for RDFT/DCT) built from a prototype array's shape and element type -- unlike Conv's `h`, there's nothing to carry values for. A storage-type change just means replanning against an uninitialized prototype on the requested backend, exactly like the existing threaded-change path already did for `Array`; the fix widens that path to any storage_type and folds the two conditions (unchanged storage AND unchanged threading) into one guard for the "just share the existing plan" fast path. NFFTOp: the harder case, since it plans against a `trajectory` array the operator doesn't retain as its own field, and rebuilding needs it back. It turns out the NFFT plan already carries it (`plan.k`, in the flattened 2D form `create_plan` reshapes every trajectory into before planning) -- this codebase's own `NFFTPlan` override already depends on that exact internal field layout to *construct* NFFT.jl's plan type in the first place, so reading `.k` back out is no more fragile than what was already there. The image size comes from the public `NFFT.size_in`, and `dcf`'s own shape is what determines how to un-flatten `plan.k` back into a paddable trajectory (reshaping to `plan`'s own flattened shape would have been wrong whenever the original trajectory had more than one non-node dimension). This also incidentally fixes the same wrong-refusal bug for `threaded` on NFFTOp, which carried an identical "cannot change after construction, rebuild instead" argument even though the same trajectory-recovery path resolves it too. Updated the two tests that asserted the old refusals (RDFT/NFFTOp in "Threading contract: batching an operator that has no threaded path" and "...NFFTOp reports its plan-time threading") to instead verify the request is now honoured and numerically correct. Verified locally (JULIA_NUM_THREADS=2): full :fftw/:nfft/:Threading/:jet run (906 tests) and a full Documenter doctest build, both passing. Manually verified NFFTOp's storage_type copy reproduces identical forward/adjoint results from a distinct plan object. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
…bumps - Move the _dsp_fftw_num_threads docstring to sit above the function itself instead of above the unrelated THRESHOLD_C2C constant, and correct its claim that threaded accepts nothing (it is a plain Bool). - Fix the same false "or the default nothing" claim in NFFTOp's constructor and _nfft_threaded docstrings. - Document every parameter of Conv/Xcorr's constructors, not just the new num_threads/threaded ones. - Revert the package version bump (0.5.0 -> 0.4.0) and the matching subpackage AbstractOperators compat bounds (0.5 -> 0.4), which were bumped without being part of an actual release. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
…g one Comments across the threading code narrated the change itself (what the previous/legacy/bespoke version did) instead of documenting current behavior, which rots as the code moves on. Rewrote each to state the current invariant directly: FFTWOperators/Shift.jl, properties.jl, batching/BatchOp.jl, batching/SimpleBatchOp.jl, calculus/BroadCast.jl, calculus/Scale.jl, linearoperators/Variation.jl, linearoperators/FiniteDiff.jl, threading_policy.jl (MIN_BATCH_WORK_FOR_PARALLEL docstring), and the corresponding test comments/testitem name in test_threading_policy.jl and test_threading_operators.jl. test_hcat.jl's comment was outright wrong, not just stale: it claimed a wrong-length permutation vector "silently corrupts invpermute!"; tested directly, it raises a clear DimensionMismatch instead. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
…text - benchmark/compare.jl: Runic had inserted a `return` before `@info`, making it look like the function returns the logged message instead of `nothing`. - benchmark/operator_thresholds.jl: use `@belapsed` instead of `median(@benchmark(...))` for the per-operator crossover measurements. - benchmark/threading_sweep.jl / operator_thresholds.jl: cross-reference each other's purpose (class-level defaults vs. per-operator overrides) now that it wasn't obvious which one was current. - benchmark/bench_common.jl: replace the BENCH_NONLIN_N::Dict with one named constant per operator, matching every other constant's naming convention in the file. - threading_policy.jl: record the Julia version the measured thresholds were swept under, alongside the existing machine/thread-count/date provenance. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
Previously, when the block-parallel direction (HCAT's adjoint, VCAT's forward) threaded, the constructor forced *every* wrapped block to threaded=false and stored only that forced-serial copy -- so the other, always-serial direction (HCAT's forward, VCAT's adjoint) lost each block's own internal parallelism for no reason, since that direction never nests under the block-parallel loop. Store both: `A` keeps each block's natural (possibly individually threaded) state and is used by the always-serial direction and by every introspection method (_children, is_threaded, domain_type, size, ...); `A_par` is the forced-serial copy, used only by the threaded block loop to avoid nesting a block's own parallelism inside it. `A_par` aliases `A` (no extra allocation) whenever the block loop itself isn't threaded. The extra type parameter is appended last so essentially none of the existing partial-parameter dispatch patterns needed touching. _copy_operator_impl no longer force-serializes children before handing them to the constructor, since the constructor now derives the forced-serial copy itself. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
supports_threading(::LBFGS) is false (it has no threaded path at all), so per the documented copy_operator/adapt_operator contract a threaded request -- true or false -- is vacuous and should be a no-op, exactly like the generic fallback treats any operator without a threaded path. _copy_operator_impl instead threw an ArgumentError for threaded=true specifically, while silently accepting threaded=false -- an asymmetric, contract-violating special case. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
BLAS's mul!(gemv!/gemm!) does thread for large enough operands, but
neither operator declared supports_threading or did anything to control
it -- BLAS just ran with whatever the ambient global thread count
happened to be. Add a threaded::Bool=true constructor keyword, resolved
through the shared size policy (_blas_threaded, using the matrix's
element count as the FLOP proxy), and scope BLAS's thread count around
each mul! call via a new _with_blas_threading helper.
That scoping goes through NestedThreading's refcounted budget
(with_full_threads/with_restricted_threads), not a raw
BLAS.set_num_threads save/restore: a mul! can run concurrently with
other mul!s (e.g. several MatrixOp blocks inside a threaded HCAT/VCAT
block loop), and a naive save/restore of a process-global from
concurrent callers is exactly the interleaving bug NestedThreading
exists to prevent. Unlike FFTW/NFFT, BLAS has no plan to bake a thread
count into, so this is a per-call scope rather than a construction-time
choice.
threading_threshold(::Type{<:MatrixOp}) has no measured override yet
(falls back to the conservative default) -- a benchmark run (square
gemv, 8 threads) shows a real win starting somewhere between 65,536 and
1,048,576 elements (3.36x at 1024x1024) and a regression below that
(0.72x at 64x64), consistent with but not pinning down the current
default; a proper operator_thresholds.jl sweep is a follow-up.
Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
…ssover mul! is two reductions (a stabilizing max, a normalizing sum) around an elementwise exp, and the Jacobian adjoint has a third (a dot product) -- not a plain elementwise map, so it couldn't use the @.. thread=Th split every other nonlinear operator uses. It was previously marked deliberately unthreaded with a comment claiming a threaded rewrite wasn't applicable "without a multi-pass rewrite with its own reductions". Polyester's `@batch reduction=((op, var), ...)` (+, max here) is exactly that rewrite: each reduction becomes its own @Batch loop over a plain local accumulator. threading_threshold(::Type{<:SoftMax}) is now measured rather than falling back to the conservative default: with the size policy bypassed to find the true crossover (sweeping through the exported threaded=true keyword alone only ever reproduces whatever the threshold already was), forward crosses at 2^9 (Float64) / 2^11 (Float32); the Jacobian-adjoint, with one more reduction pass, crosses later, at 2^12 (Float64) / 2^11 (Float32). Set to 2^12, the latest of all four, matching how every other multi-direction/multi-dtype threshold in this package is chosen. is_thread_safe stays false: buf is shared, mutable state written by both mul! paths, which is an orthogonal concern to whether a single call threads internally. Also fixes the same throw-on-threaded=true bug as the LBFGS commit: _copy_operator_impl now treats the request as vacuous rather than erroring, though it's now moot for the threaded=true case since SoftMax genuinely supports threading. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
…h singleton Atan, Cos, Exp, Pow, Sech, Sigmoid, Sin, SoftPlus, and Tanh each had two near-identical mul! methods (Th=false/Th=true dispatch, differing only in the @.. thread=false/true literal). Collapse each pair into one method parameterized by Th, using @.. thread = Th -- the same pattern DiagOp/Scale already use. This requires switching Th's encoding from a plain Bool to FastBroadcast.True()/False() (via the existing _fbthread helper): `@..`'s thread= argument only wraps a *literal* true/false into the singleton type it dispatches on; passed a variable, it forwards the variable's runtime value unchanged, and fast_materialize! has no method for a plain Bool. is_threaded/_copy_operator_impl convert back via the existing _fbbool helper. SoftMax is deliberately not unified: its threaded path is a structurally different reduction-based rewrite (see the previous commit), not a thread= toggle on the same expression. Also switched the serial branch from plain broadcasting (y .= f.(x)) to @.. thread=false, for one broadcast mechanism across both branches. Benchmarked serial-mode FastBroadcast against plain broadcasting first (sin and cheap-arithmetic kernels, 2^4-2^22, Float64/Float32): no measurable difference (ratios ~1.00x +/- 10% noise, no consistent direction), so this is purely a consistency simplification, not a performance change, and needed no threshold remeasurement. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
- "leaf operators without a threaded path" listed MatrixOp, which now genuinely supports threading (a prior commit added real BLAS threading support), so `supports_threading(op) == false` no longer holds for it. Remove it from that list and add a dedicated threading-contract test for MatrixOp/LMatrixOp covering is_threaded/supports_threading, copy_operator/adapt_operator round-tripping, and numerical correctness -- forward (gemv) is bit-identical between serial and threaded since BLAS splits it by output row, but the adjoint (gemv on A') is not, so that check uses ≈ rather than ==, both confirmed by direct measurement. - "VCAT-forward and HCAT-adjoint block loops" asserted `all(!is_threaded, nested.A)`, which was the correct check before HCAT/VCAT split their block storage into `A` (natural, used by the always-serial direction) and `A_par` (forced serial, used only by the threaded block loop). It's `A_par` that must be all-serial now; `A` legitimately stays threaded since nothing nests under the always- serial direction. Updated the assertion and its comment accordingly. Fixes: https://github.com/kul-optec/AbstractOperators.jl/actions/runs/32706205442/job/97367749529 Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
…atrixOp/LMatrixOp The convenience constructors forwarded the new threaded keyword via a generic `kwargs...` splat instead of naming it explicitly. That is exactly the kwarg-to-kwarg forwarding pattern this repo's JET guidelines warn about: it breaks static resolution of the downstream `array_type`/`threaded` keywords, cascading into runtime dispatch in _normalize_array_type and _blas_threaded that JET's @test_opt correctly flagged (`@test_opt constructors: MatrixOp(M)`, 4 possible errors). Name `threaded::Bool = true` explicitly in every MatrixOp/LMatrixOp constructor overload instead, matching every other operator's convention in this codebase. Confirmed locally: test_opt_constructors.jl 33/33 (was 32/33), full jet/ suite and MatrixOp/LMatrixOp tests clean. Fixes: https://github.com/kul-optec/AbstractOperators.jl/actions/runs/32709361574/job/97377268735 Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
…tform The new MatrixOp/LMatrixOp threading test asserted exact equality (==) between serial and threaded results for MatrixOp's forward direction and for LMatrixOp's gemm. CI (Julia LTS, ubuntu) showed LMatrixOp's gemm reassociates its reduction across threads -- not bit-identical, unlike what local measurement showed for MatrixOp's gemv. Since MatrixOp's own "bit-identical" claim was also only ever verified on one local machine/BLAS build, weaken both to ≈ rather than risk the same surprise on a different platform/BLAS build in the CI matrix (windows/macOS runners, or a different OpenBLAS version). Fixes: https://github.com/kul-optec/AbstractOperators.jl/actions/runs/32713440255/job/97389566230 Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01JHYpJyvm5XwhtMXLE1N8i9
Addresses 11 correctness/perf findings from the unified-threading review: - SignAlternation had no _copy_operator_impl, so copy_operator/BatchOp/HCAT/etc. threw for any threaded instance. - MatrixOp's _copy_operator_impl dropped the NC (n_colons) batching parameter. - RDFT re-derived `dims` by diffing dim_out/dim_in, which silently picked the wrong axis when the transformed dimension has length 1 or 2; now stored as a type parameter like IRDFT already does. - SoftPlus's (DomainDim::NTuple) constructor was missing the `threaded` kwarg. - SimpleBatchOp/SpreadingBatchOp copies shared the wrapped operator(s) (and their scratch buffers) with the original instead of always copying. - NoOperatorBroadCast's copy spliced a bare storage wrapper into a slot that needs the full parameterized type, losing the element type. - Scale's remove_displacement/permute/get_normal_op silently re-enabled threading by not forwarding the current Th veto. - Pow's integer-exponent threading threshold was dead code: the elementwise threading resolver was called with the bare `Pow` UnionAll, which never dispatches to the Integer-specific threshold method. - SimpleBatchOpMultiThreaded's get_normal_op parameterized opT from the pre-adaptation operator instead of the stored (threading-forced-off) ones. - SpreadingBatchOp's copy dropped threading_strategy, always falling back to AUTO. - DFT's replanning path (storage/thread-count change) dropped flags/timelimit, silently degrading a MEASURE/PATIENT plan back to ESTIMATE. Verified each fix with targeted reproductions of the review's failure scenarios, then ran the filtered suite: julia --project=test test/runtests.jl ":SignAlternation,:MatrixOp,:RDFT,:SoftPlus,:SimpleBatchOp,:SpreadingBatchOp,:BroadCast,:Scale,:Pow,:DFT" 1842 passed, 1 failed (failure isolated to a stale .temp/base snapshot file, not the live test suite). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01KP9PMHZeEvckXJo45NPveg
Dedup Scale's threaded-rebuild logic into _rethread_scale, derive MatrixOp's copy domain shape from size(op,2) instead of re-deriving it, and make SpreadingBatchOp's _threading_strategy_of exhaustive so a future leaf type errors instead of silently falling back to AUTO. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
hakkelt
marked this pull request as ready for review
August 24, 2026 19:44
…nstructors SimpleBatchOp: no CI test ever crossed MIN_BATCH_WORK_FOR_PARALLEL through the real BatchOp(...) constructor path, so create_BatchOp's threaded branch and the actual @budgeted_threads mul! loop went unexercised (the one test large enough, benchmark_threading, is skipped when CI=="true"). Add a correctness-only test that crosses the threshold. Nonlinear operators (Atan/Cos/Exp/Sech/Tanh/Sin): the tuple-only constructor (no domain_type) and the construct-from-array constructor were never called by any test. Exercise both, and extend the shared test_NLop helper to check show/fun_name, is_thread_safe, and supports_threading for every nonlinear operator. Verified via local --code-coverage=user + Coverage.jl runs (before/after) and :batching + :nonlinearoperator regression (569/569 passing, GPU excluded). Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MGe9CxEyPmji64JYLtN8DN
WaveletOp only defined is_threaded/supports_threading, so copy_operator requests with a storage_type would hit the generic fallback and throw. Add a _copy_operator_impl that rebuilds the storage-tracking type parameter (threaded stays a no-op, matching supports_threading = false), plus a copy_operator test. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MGe9CxEyPmji64JYLtN8DN
hakkelt
force-pushed
the
unified-threading
branch
from
August 26, 2026 23:05
0787719 to
525f485
Compare
…default CI's benchmark job (JULIA_NUM_THREADS=2) silently exercised each operator's own size-based threading policy for whichever benchmarks lacked a manual -single/-threaded split, while BENCH_THREADED (gating the manually split entries) was unconditionally disabled under CI. So a "default" CI benchmark run measured neither reliably-serial nor reliably-threaded execution -- it measured whatever each operator's policy happened to decide at 2 threads, invisibly. compare.jl now runs the whole suite once per thread count (`-t 1`, `-t 2` by default, via --threads), pinning the count explicitly per subprocess launch rather than relying on JULIA_NUM_THREADS (which -t always overrides anyway). BENCH_THREADED simply follows Threads.nthreads() > 1, so the threaded-suffixed entries are exercised for real under -t 2 instead of being skipped. Each thread count gets its own summary + time/memory tables in the PR comment body, so a "1 thread" comparison and a "2 threads" comparison are both visible and neither is contaminated by the other. Simplified DiagOp/FiniteDiff/Variation/SimpleBatchOp accordingly: they used the same input size for their -single and -threaded variants, so under this per-thread-count design the plain default-threaded entry already gives a genuinely serial measurement at -t 1 and a genuinely threaded one above that -- the explicit veto/check_threaded pair was pure duplication once thread count itself is the controlled axis. Left the split in place wherever it does real work: a dedicated larger size for the threaded path (Scale, HCAT/VCAT/DCAT, DFT, NFFTOp), multiple threading strategies to compare (SpreadingBatchOp), or an explicit veto's overhead under otherwise-available parallelism. Running the suite twice roughly doubles CI cost, so the per-leaf `seconds` budget is picked from Threads.nthreads() (1.0s at one thread, 0.5s above that) instead of a flat 5s, and the benchmark manifest is instantiated once per revision instead of once per (revision, thread count) subprocess. Measured end to end through compare.jl: ~19 minutes, confirmed by two full runs -- a few minutes over the 10-15 min target, accepted because a correct, unambiguous single/multithreaded comparison is worth more than shaving the estimate further. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MGe9CxEyPmji64JYLtN8DN
GPUEnv 0.2.1 (registered today, a caching/idempotency rework of activate() in hakkelt/GPUEnv.jl#5) breaks GpuEnvSetup's `GPUEnv.activate(; persist = true)` (test/gpu_env_setup.jl): the active project ends up on neither the GPU overlay nor test/Project.toml afterward, so every testitem that runs later in the same process fails with "Package AbstractOperators ... does not seem to be installed". Reproduced twice in CI (run 33018812036) against the same commit with no code changes in between, ruling out a one-off flake. Pinned to "0.1 - 0.2.0" rather than "0.1, 0.2.0": a bare compat entry is caret-bounded (0.2.0 alone means >=0.2.0 <0.3.0 and would still admit 0.2.1), so only an explicit range or exact pin actually excludes it -- verified with Pkg.Types.semver_spec before picking this form. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MGe9CxEyPmji64JYLtN8DN
hakkelt
force-pushed
the
unified-threading
branch
from
August 27, 2026 05:11
525f485 to
c9c9fda
Compare
Reproduced the original bug locally to confirm both the diagnosis and the fix, on Julia 1.10 (the "lts" CI matrix leg, the only one that actually hits it -- the same scenario on Julia 1.12 doesn't reproduce it): pinning GPUEnv exactly to 0.2.1 and running the same develop+activate(;persist=true)+`using DSPOperators` sequence tests.yml's LTS job runs fails identically to the CI log; pinning to 0.2.2 from hakkelt/GPUEnv.jl#6 (branch fix-instantiate-skip-after-develop) instead succeeds. test/Project.toml now points GPUEnv's `[sources]` at that branch and widens compat to admit 0.2.2 alongside the existing "not 0.2.1" exclusion. That alone isn't enough, though: `[sources]` requires Pkg >= 1.11 and is silently ignored on Julia 1.10 -- confirmed by testing both the inline-table form and the `[sources.GPUEnv]` dotted form Pkg itself rewrites it to, neither took effect there. Since 1.10 is the one job that needs the pin, tests.yml's "Prepare LTS test environment" step now also does an imperative `Pkg.add(url=, rev=)` for GPUEnv, which is old, version-independent API and is what actually exercises the fix in CI. Once #6 merges and 0.2.2 is registered: drop the `[sources]` entry and the imperative Pkg.add in tests.yml, and relax the compat bound to admit registered 0.2.2+ normally. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MGe9CxEyPmji64JYLtN8DN
CI testing of hakkelt/GPUEnv.jl#6 (pinned via Pkg.add(url=, rev=) in the LTS job, previous commit) found the branch's first fix incomplete: the real failure needed a second round from the maintainer (two Julia-1.10-specific defects from the earlier PR #3 compat work, not just the Pkg.instantiate() skip from PR #5), released together as the registered 0.2.2. Verified: a fresh `Pkg.resolve()` against a plain `GPUEnv = "0.2.2"` compat bound resolves cleanly from the registry. Drops the [sources] git-branch pin and its now-inaccurate comments, the imperative Pkg.add workaround in tests.yml's "Prepare LTS test environment" step, and the "0.1 - 0.2.0, =0.2.2" compat gymnastics -- back to a single ordinary compat entry. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MGe9CxEyPmji64JYLtN8DN
…ntwise ops) Same pattern as the earlier DiagOp/FiniteDiff/Variation/SimpleBatchOp simplification, just missed on the first pass: DFT, NFFTOp, and every elementwise nonlinear operator (Pow, Exp, Sin, Cos, Atan, Tanh, Sech, Sigmoid, SoftMax, SoftPlus) used the identical input size for their -single and -threaded variants, so the split was pure duplication once the per-thread-count run (see BENCH_THREADED) is what actually controls serial vs threaded execution. Collapsed each to a single default-threaded entry; dropped the now-unused `threaded`/`threadable` parameters and check_threaded gates along with them. Scale/HCAT/VCAT/DCAT (dedicated larger threaded-only size) and SpreadingBatchOp (three named threading strategies, no single default to fall back on) are unaffected -- their split does real work, unlike the ones collapsed here. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MGe9CxEyPmji64JYLtN8DN
The PR benchmark comparison flagged five time regressions and seventeen memory
regressions. Four independent causes, each measured rather than guessed.
BLAS pinned to one thread on every MatrixOp/LMatrixOp mul!
`_blas_threaded` resolved through `_default_threaded`, which requires
`Threads.nthreads() > 1` and `length(A) >= 2^18`. Both are the wrong question to
ask of BLAS: it owns a thread pool separate from Julia's, and it already applies
its own per-call size heuristic. LinearAlgebra initialises it to
`max(1, jl_effective_threads() / 2)`, driven by the machine and CPU affinity and
never by `-t`, so a `-t 1` session still has a fully threaded BLAS that the
`nthreads()` gate was pinning to one. The 192x192 benchmark matrix (36864
elements) also fell far below the element threshold.
Worse, the *permitted* branch called `with_full_threads`, which costs ~5.5us per
call in lock and setter round trips and can raise BLAS past a count the caller
deliberately set.
`threaded = true` now opens no scope at all; only the `false` veto opens
`with_restricted_threads`, which is the direction nesting safety depends on.
The policy's only remaining say is that BLAS threading is meaningless for
non-CPU storage.
MatrixOp 192^2, BLAS=4: 354.8us / 64 B -> 188.8us / 0 B
(raw gemm for reference: 186.5us)
threaded = false still pins to one: 362.9us (raw BLAS=1: 356.4us)
LMatrixOp 1024: 0 allocs; BLAS=1 128.5us, BLAS=4 55.1us
This also covers LMatrixOp/forward and all six Ax_mul_B* entries, whose children
are MatrixOps, and the one- and two-allocation memory regressions on those rows.
Two situations leave BLAS at one thread through no fault of this package: fewer
than four effective CPUs, and a Distributed worker, where Distributed calls
`Base.disable_library_threading()` whose hook is `BLAS.set_num_threads(1)`. Both
are deliberate policy from outside, so `threaded = true` leaves them alone rather
than overriding them. Documented on `_with_blas_threading`.
Sigmoid Jacobian-adjoint evaluated exp twice per element
Fusing the multi-pass kernel into a single broadcast duplicated
`exp(-gamma * x)`, and a broadcast evaluates each occurrence separately. `exp`
dominates this kernel, so the fused form cost roughly double. `_sigmoid_jac_adj`
binds it to a local, keeping the single fused pass over the arrays.
inline double-exp 419.6us -> helper 239.8us (Base broadcast: 282.6us)
threaded, 4 threads: 69.1us
Output is bit-identical to the pre-PR formulation. The helper stays on
FastBroadcast's fast path: a function call does not change the broadcast style,
which is what its fallback keys on.
THRESHOLD_C2C was four to six powers of two too low for Conv/Xcorr
The value was borrowed from FFTWOperators' bare-transform crossover and flagged
provisional in its own docstring. A Conv/Xcorr mul! is two transforms plus a
pointwise product, so FFTW's threaded plan synchronises twice; its crossover sits
well above a single plan_fft of the same length. Swept over the padded length at
one versus four FFTW threads:
fftlen 2^15: 0.89x 2^17: 0.96x 2^19: 1.12x 2^21: 1.30x
Raised to 2^19, the first swept length that actually wins, and the PROVENANCE
line rewritten as measured.
Building an above-threshold Xcorr takes about 55 seconds under FFTW.MEASURE, far
too slow for a contract test, so the above-threshold half of that contract now
asserts on `_dsp_fftw_num_threads` directly and the operator-level checks use
`num_threads` to get a threaded plan cheaply.
NFFT threading never pays with two workers
Its gate was `nthreads() > 1` with no size component, which its own comment
called out. Sweeping one process per thread count shows the deciding variable is
the thread count, not the workload (ratios are serial/threaded):
threads | 48^2 fwd/normal | 96^2 fwd/normal | 192^2 fwd/normal
2 | 0.59 / 0.85 | 0.78 / 0.99 | 0.91 / 0.90
3 | 0.87 / 0.81 | 1.03 / 1.07 | 2.29 / 1.41
4 | 0.87 / 0.96 | 1.93 / 1.56 | 1.43 / 2.21
At two workers threading loses at every size measured, up to a 5ms mul!, so
growing the workload does not rescue it. Added `MIN_THREADS_FOR_NFFT = 3`. The
threaded path also allocates 4-12 KiB per mul! against 112 B for the serial one.
No size gate was added for the 48^2 column despite it sitting at or below 1.0
everywhere: a repeat of the four-thread run put it at 1.20 / 1.18, so on a shared
machine it is inside the noise and a threshold fitted to it would be curve
fitting rather than measurement.
Benchmark workflow comment
Its stated rationale for omitting OPENBLAS_NUM_THREADS was that the operators
clamp BLAS to `Threads.nthreads()` per mul!. That clamped to `threadpoolsize()`,
and after this change the permitted branch clamps nothing. Replaced with the
actual mechanism, plus the reason not to set the variable there: doing so makes
LinearAlgebra skip its own default entirely, which would pin BLAS for the whole
run and stop these operators measuring what a user gets.
Not changed, deliberately: the ZeroPad and Eye deltas (mul! untouched by this PR,
and the one- and two-thread rows disagree in direction); the two-allocation rows
on the elementwise Jacobian-adjoints, which are intrinsic to `@.. thread = true`
and bought two- to fivefold time wins; and the `-threaded` benchmark entries,
which have no base counterpart and whose per-element cost is explained by cache
residency at 128x and 4x the serial workload.
Tests: 718 passed, 0 failed, filtered on
:Threading,:MatrixOp,:LMatrixOp,:Sigmoid,:dsp,:nfft at -t 2.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01MGe9CxEyPmji64JYLtN8DN
These changes were already present, uncommitted, in the unified-threading
worktree before the benchmark regression work; they are committed separately so
they stay reviewable on their own.
`domain_array_type`/`codomain_array_type` on NFFTOp returned the raw buffer types
(`typeof(op.plan.tmpVec)` and the `K` parameter) rather than a wrapper type
parameterised by the operator's own element type. Routing both through
`_array_wrapper_type` and re-parameterising with `domain_type`/`codomain_type`
brings them in line with the convention the rest of the package follows.
Two testitems cover paths nothing reached before:
- `BatchOp(operators, mask::NTuple{M,Symbol})`, the two-positional-argument
overload with a bare non-Pair mask and no batch size, which forwards to
`BatchOp(operators, (), mask => mask)`.
- `_policy_storage`/`_storage_eltype_or_float`, including the real call path
where a multi-domain operator's `domain_array_type` is an `ArrayPartition` and
`_should_thread` must still resolve to a Bool rather than error on it.
Tests: 346 passed, 0 failed on :batching; the threading-policy and NFFT items are
covered by the 718-pass run recorded on the preceding commit. Both at -t 2.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01MGe9CxEyPmji64JYLtN8DN
Finishes the -single/-threaded cleanup that 02794e7 applied to DiagOp, Variation, DFT and the elementwise operators. BroadCast was missed, and it turned out to have two different problems rather than one. `identity-single` and `identity-threaded` were the plain redundancy: same operator, same size, differing only by the keyword. `BENCH_CALC_N * 8` is exactly THRESHOLD_MEMORY_BOUND, so a single default-constructed entry is serial at `-t 1` and threaded at `-t 2` and reproduces the pair exactly. Measured on the merged entry: `is_threaded` false at one thread and true at two, 70.9us against 43.3us. `operator-threaded-forward` and `operator-threaded-adjoint` were the worse case. At 256 * 8 = 2048 elements they sit three orders of magnitude below THRESHOLD_MEMORY_BOUND, so the policy declines `threaded = true` at every thread count and both entries were measuring the serial path under a threaded name. The run that prompted this shows it plainly: 261 ns for `operator-threaded-forward` against 270 ns for `operator-single-forward`, and 1.31us against 1.29us on the adjoint. Neither carried a `check_threaded` guard, which is exactly the failure that guard exists to make loud rather than silent. They are dropped rather than resized. The pair measures per-child dispatch overhead in a broadcast over a small wrapped operator, which is a different cost shape from `identity`; sizing it up to clear the threshold would measure the same memory-bound kernel twice instead of keeping the small-operator case. Not touched: the `-threaded` entries on Scale, HCAT, VCAT and DCAT. Those are not redundant with their plain counterparts, because they deliberately run at a different workload chosen to clear a gate the shared size cannot -- Scale at 2^22 against its own crossover, the CATs at four blocks against MIN_BLOCKS_FOR_PARALLEL. The plain entries can never reach those paths, and all four are check_threaded/check_block_threaded guarded so a mis-size fails loudly. `SpreadingBatchOp`'s `-single` entries stay too: there `-single` is the serial baseline of a three-way threading-strategy comparison, not a thread-count split. This renames keys, so these rows will show no base counterpart in this PR's comparison table until the base script catches up -- the same trade 02794e7 already accepted. Verified by including bench_common.jl and suites/calculus.jl at `-t 1` and `-t 2` and running each BroadCast entry, which also exercises the setup blocks. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01MGe9CxEyPmji64JYLtN8DN
Collaborator
Author
|
@lostella, this PR is also ready for review. Sorry, this is quite a large PR again, but the goal is straightforward: add multithreading support for all operators where possible via a unified API (and fix minor issues encountered along the way while implementing the feature). The proposed changes don't break the current API, and the default value of the parameter governing multithreading aims to maximize out-of-the-box performance. Based on the operator type and input size, it enables multithreading only when preliminary benchmarks indicate it is beneficial. The goal of adding the parameter named This is the suggested unified API: op1 = Variation(Float64, (64,64))
x, y = rand(64, 64), (4096, 2)
mul!(y, op1, x) # this runs multithreaded
op2 = Variation(Float64, (64,64); threaded=false)
x, y = rand(64, 64), (4096, 2)
mul!(y, op2, x) # this runs single-threaded
op3 = Variation(Float64, (10,10))
x, y = rand(10, 10), (100, 2)
mul!(y, op3, x) # this runs single-threaded because the input is too small to benefit from multithreading (2^10 elements is the threshold for Variation)A side-effect of this work is that the |
lostella
self-requested a review
August 29, 2026 05:17
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
This PR unifies the interface for enabling/disabling threading across operators and adds multithreading support to all possible operators, replacing scattered per-operator implementations with a consistent policy infrastructure.
Important Changes
Threading Policy Infrastructure
src/threading_policy.jlmodule introducing a unified threading framework:is_threadedtrait for operators with optional threading supportsupports_threading,default_threaded, andthreading_thresholdper-operator policiesadapt_operatorandcopy_operatorfunctions for consistent threading controlbenchmark/threading_sweep.jl) replace previous guessesOperator Threading Coverage
Elementwise nonlinear operators (9 operators: Pow, Exp, Sin, Cos, Atan, Tanh, Sech, Sigmoid, SoftPlus)
FiniteDiff optimization: eliminated allocations from
b[idx_1] .- b[idx_2]materialization using@viewsCalculus operators (Sum, HCAT, VCAT, DCAT, Compose, Scale, BroadCast, etc.) add forwarding traits
Batching operators (BatchOp, SimpleBatchOp, SpreadingBatchOp)
Bug Fixes
_per_thread_operators(branch onis_thread_safe)Test Expansion
@test_optand@test_callfor threading APIsInfrastructure & Cleanup
Performance Baselines
benchmark/threading_sweep.jl) captures 243 measurements for future regression detection