Real plans: mul!, allocation-free execution, dims path and a real 2D path - #132
Open
pankgeorg wants to merge 3 commits into
Open
Real plans: mul!, allocation-free execution, dims path and a real 2D path#132pankgeorg wants to merge 3 commits into
pankgeorg wants to merge 3 commits into
Conversation
With FFTW.jl loaded alongside FFTA, plan_rfft(::Vector{Float64}, ::Int) was
ambiguous between FFTW's StridedArray method and FFTA's method annotated
with region::RegionTypes, turning rfft(x) into a MethodError. Leave region
unannotated on the AbstractFFTs entry points (as plan_fft already does) and
normalise it in an internal function, so FFTW's methods are strictly more
specific and take over as AbstractFFTs intends.
A coexistence test runs in a subprocess (loading FFTW in the test process
would make every other test exercise FFTW).
Real-input/real-output plans only implemented *, so every rfft/irfft allocated, mul!(y, p, x) with a preallocated output was a MethodError, rfft along one dimension of an N-d array went through mapslices, and 2D real plans ran a full complex transform and discarded half of it. FFTAPlan_re now carries a scratch buffer and two pencil kernels (_rfft_pencil!/_brfft_pencil!) implement the even-length half-size trick and the odd-length transform on AbstractVectors, so views work as input and output. mul! is defined for 1D plans on 1D and N-d arrays (looping over pencils along the region) and for 2D plans on N-d arrays (real transform along the first region dimension, complex along the second), and * allocates the output and calls mul!. The 2D plan's first call graph is built for the half length like the 1D plan's.
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #132 +/- ##
==========================================
- Coverage 98.80% 97.52% -1.29%
==========================================
Files 5 4 -1
Lines 585 646 +61
==========================================
+ Hits 578 630 +52
- Misses 7 16 +9 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
…nels Replacing mapslices with strided views made real transforms along dims=2 of wide matrices (64 x N) 1.2-1.3x slower on x86-64: the mapslices copy had been an unlabelled copy-in that turned a stride-of-a-cache-line gather into one contiguous pass before the kernel. Pencils whose parent arrays are unit-stride along the transform dimension still go to the kernels directly (the dims=1 gain stays); any other pencil is copied to two plan-owned contiguous buffers first and copied back after, so execution stays allocation-free and the dims=2 result is now identical to the mapslices one.
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
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.
Item B of the plan in #130: real-input/real-output plans get a proper
mul!, allocation-free execution, an efficientdimspath, and a real 2D path. Stacked on #131 (sameplan_rfft/plan_brfftentry points); the diff includes that commit until it lands.Before
FFTAPlan_reonly implemented*, notmul!, so every real transform allocated (74 KiB–67 MiB per call in the benchmark sweep) and consumers that pass a preallocated output (mul!(Y, plan, X), which is how DSP.jl uses plans at 10 of its 13 call sites) got aMethodError;mul!into aSubArraywas not possible either.rfft/brfftalong one dimension of an N-d array went throughmapslices, allocating per pencil (≈10× FFTW).Now
FFTAPlan_recarries a scratch buffer (sized by_re_buflen:n/2for an even forward plan,nbackward,2nodd) and two pencil kernels,_rfft_pencil!/_brfft_pencil!, implement the even-length half-size trick (unchanged maths from the previous*) and the odd-length full transform onAbstractVectors, so views work as input and output.mul!(y, p, x)for 1D plans on 1D and N-d arrays (looping over pencils alongdims, type-stable via the same@nifpattern as the complex path) and for 2D plans on N-d arrays (real transform along the first region dimension, then the complexfft_along_dim!along the second; backward does the reverse).*allocates the output and callsmul!.test/argument_checking.jlpasses unchanged); new tests intest/real_mul.jlcovermul!for 1D/N-d/2D plans, views, zero allocations for 1D plans (Bluestein sizes still allocate their scratch until Performance roadmap: closing the gap to FFTW (proposed PR sequence) #130-A lands),Float32, and dimension errors.Base.complex(::FFTAPlan_re)is gone (nothing used it).Before/after (aarch64 Neoverse-N1, Julia 1.12.6,
benchmark/suite.jl --kinds rfft, planned execution, single thread;FFTA/FFTWvs FFTW 3.3.11ESTIMATE):Worst / best individual cases:
Strided pencils (found by the x86-64 companion run). The first version of this PR fed the real pencil kernels strided views directly. On an AVX2 x86-64 machine that made
rfftalongdims=2of wide64×Nmatrices 1.2–1.3× slower thanmain(whiledims=1got 1.7× faster), and a back-to-back five-branch probe there attributed it to this PR: the oldmapslicescopy had been an unlabelled copy-in optimisation — it turned a stride-of-a-cache-line gather into one contiguous pass before the kernel. The final version keeps the direct path only for pencils whose parent arrays are unit-stride along the transform dimension and copies every other pencil to two plan-owned contiguous buffers first (so execution stays allocation-free and thedims=2result is bit-identical tomapslices). With the fix this PR is at parity withmainon the strided shape and keeps the contiguous win: x86-64 back-to-back,rfftFloat6464×16384dims=237.1 ms (main) / 43.6 (pre-fix) / 38.3 (fixed), 64×65536 290 / 389 / 276 ms, whiledims=1stays 3.8–3.9× faster thanmain; aarch64 shows the same picture (29.8 → 27.1 ms and 157 → 148 ms againstmain's 28.2 / 157). The 1.3–1.7× improvement of that shape overmainarrives only with the rest of the stack (#134/#135/#137). Verified on x86-64 with 8 threads that the per-worker buffers of #137 keep the results bit-for-bit identical to the single-threaded ones at zero allocation.The 1D
rffttime itself is unchanged — it is dominated by the underlying half-length complex transform, which is what #130's items A and C address — so this PR is mostly about allocations, thedimspath and 2D.