Plan-owned per-thread workers, direct strided pencil reads, threading across pencils

[pankgeorg]

Item D of the plan in #130, stacked on #133 (the diff includes #132 and #133 until they land).

Before, every multidimensional mul! allocated two pencil buffers per call, copied each pencil into a contiguous buffer, transformed, and copied it back; a 1D plan applied along one dimension of a matrix transformed the pencils one after the other; and a plan's single call-graph workspace made it impossible to spread the work over threads.

Now

• A plan owns one Worker per thread it may use: call graphs that share the plan's nodes and twiddle tables but have their own workspace and Bluestein work arrays, plus the pencil buffers. The count is the new num_threads keyword of plan_fft/plan_rfft/… (same name as FFTW.jl's; default Threads.nthreads(); inv(p) keeps it).
• One driver, _foreach_pencil, runs every pencil loop: 1D plans on N-d arrays, N-d plans, and the real-input pencil kernels. Below THREAD_THRESHOLD (2^15 elements) or with one worker it runs serially and allocation-free; above it, pencils are split into one contiguous chunk per worker and the chunks run as Threads.@spawn tasks. Each chunk always uses its own worker, so results are bit-identical to the single-worker result and do not depend on the number of threads (tested with 1, 2, 3 and 5 workers on a 1-thread runtime and on 4 threads).
• Kernels read strided pencil views directly (they take strides anyway), so only the copy-out from the worker's output buffer remains; the per-call buffer allocations are gone (0 bytes/exec single-threaded; the threaded path allocates a few KiB per call for the tasks).
• The real 2D path and the rfft/brfft along dims path use the same machinery, so batched real transforms — DSP.jl's Welch/periodogram/fftfilt shape — are threaded too.

Numbers (aarch64 Neoverse-N1, Julia 1.12.6, ComplexF64 unless noted, planned execution; 1 vs 4 workers on julia -t 4):

Single worker vs main (nd/batched sections, FFTA 1 worker; the change here is the buffer/copy handling only):

class	kind	type	cases	speedup vs main geomean (min–max)	max bytes/exec main → D
2d	fft	Float32	9	0.99× (0.89–1.07)	32 KiB → 0
2d	fft	Float64	14	0.98× (0.79–1.05)	9 MiB → 9 MiB
2d	rfft	Float32	9	2.08× (1.59–2.56)	80 MiB → 0
2d	rfft	Float64	14	1.96× (1.12–3.01)	160 MiB → 6 MiB
3d	fft	Float32	5	1.01× (0.93–1.10)	2 KiB → 0
3d	fft	Float64	5	0.96× (0.89–1.01)	4 KiB → 0
batched_dim1	fft	Float32	6	1.01× (1.00–1.01)	0 → 0
batched_dim1	fft	Float64	6	1.00× (0.99–1.02)	0 → 0
batched_dim1	rfft	Float32	6	1.11× (1.04–1.21)	32 MiB → 0
batched_dim1	rfft	Float64	9	1.14× (1.07–1.23)	65 MiB → 0
batched_dim2	fft	Float32	6	0.99× (0.95–1.02)	0 → 0
batched_dim2	fft	Float64	6	1.00× (0.95–1.02)	0 → 0
batched_dim2	rfft	Float32	6	1.13× (1.07–1.23)	32 MiB → 0
batched_dim2	rfft	Float64	6	1.10× (0.96–1.27)	65 MiB → 0

Threading (8 threads, julia -t 8; FFTA planned with num_threads=8, FFTW with 8 threads):

size	dims	FFTA 1 worker	FFTA 8 workers	FFTA speedup	FFTW 1 thr	FFTW 8 thr	FFTA/FFTW (8 thr)	bytes/exec (8 workers)
256×256	1,2	2.92 ms	523.3 µs	5.59×	1.23 ms	207.9 µs	2.52×	8 KiB
512×512	1,2	19.14 ms	3.09 ms	6.19×	6.78 ms	930.7 µs	3.32×	8 KiB
1024×1024	1,2	66.06 ms	11.76 ms	5.62×	32.94 ms	4.30 ms	2.73×	8 KiB
2048×2048	1,2	416.14 ms	64.80 ms	6.42×	178.77 ms	26.41 ms	2.45×	8 KiB
1024×64	1	1.37 ms	224.0 µs	6.11×	434.8 µs	75.2 µs	2.98×	4 KiB
4096×64	1	6.86 ms	1.03 ms	6.68×	2.18 ms	356.2 µs	2.88×	4 KiB
16384×64	1	36.94 ms	4.81 ms	7.67×	16.58 ms	2.05 ms	2.34×	4 KiB
65536×64	1	185.59 ms	24.14 ms	7.69×	95.12 ms	14.51 ms	1.66×	4 KiB

With 8 workers FFTA's speedup over its single-worker time (5.6–7.7×) matches FFTW's own 8-thread speedup (5.9–6.8×), so the FFTA/FFTW ratio in the threaded case is back to the single-threaded 1.7–3.3× instead of the 7–20× measured before this PR. The single-worker rows are the reference: no change is expected for complex transforms (the buffer copies were a small fraction), and none is seen beyond noise; the real-transform rows carry #132's gains. Bluestein sizes still allocate in this branch (that is #134's job).

x86-64 note. An independent run of the fully merged branch on an AVX2 machine showed wide 64×N rfft along dims=2 slower than main; a five-branch probe there attributed it to #132's strided real pencil path, not to this PR (this PR's dims=1 control was neutral, 1.03–1.04×, and D/B was 0.95–1.04× on the affected rows), and #132 now copies strided pencils to contiguous buffers. With that fix in place the merged stack was re-verified on x86-64: the strided rows are 1.34–1.66× faster than main, the contiguous rows 3.8–3.9× faster, and the results at 8 threads are bit-for-bit identical to the single-threaded ones (no shared buffers) at zero allocation.

Not in this PR: threading within a single 1D transform (FFTW does this above ~2^16 elements); that needs a breadth-first pass structure and is the remaining item where FFTW's threading is ahead.

[codecov]

Codecov Report

❌ Patch coverage is 96.07843% with 10 lines in your changes missing coverage. Please review.
✅ Project coverage is 97.85%. Comparing base (7aeb327) to head (009c951).
⚠️ Report is 4 commits behind head on main.

Files with missing lines	Patch %	Lines
src/plan.jl	96.03%	10 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #137      +/-   ##
==========================================
- Coverage   98.80%   97.85%   -0.96%     
==========================================
  Files           5        4       -1     
  Lines         585      652      +67     
==========================================
+ Hits          578      638      +60     
- Misses          7       14       +7     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.