Skip to content

Real plans: mul!, allocation-free execution, dims path and a real 2D path - #132

Open
pankgeorg wants to merge 3 commits into
JuliaMath:mainfrom
JuliaComputing:feat/real-plan-mul
Open

Real plans: mul!, allocation-free execution, dims path and a real 2D path#132
pankgeorg wants to merge 3 commits into
JuliaMath:mainfrom
JuliaComputing:feat/real-plan-mul

Conversation

@pankgeorg

@pankgeorg pankgeorg commented Aug 29, 2026

Copy link
Copy Markdown

Item B of the plan in #130: real-input/real-output plans get a proper mul!, allocation-free execution, an efficient dims path, and a real 2D path. Stacked on #131 (same plan_rfft/plan_brfft entry points); the diff includes that commit until it lands.

Before

  • FFTAPlan_re only implemented *, not mul!, 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 a MethodError; mul! into a SubArray was not possible either.
  • rfft/brfft along one dimension of an N-d array went through mapslices, allocating per pencil (≈10× FFTW).
  • 2D real plans ran a full complex transform and discarded half of it.

Now

  • FFTAPlan_re carries a scratch buffer (sized by _re_buflen: n/2 for an even forward plan, n backward, 2n odd) 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 on AbstractVectors, so views work as input and output.
  • mul!(y, p, x) for 1D plans on 1D and N-d arrays (looping over pencils along dims, type-stable via the same @nif pattern as the complex path) and for 2D plans on N-d arrays (real transform along the first region dimension, then the complex fft_along_dim! along the second; backward does the reverse). * allocates the output and calls mul!.
  • Error messages and argument checks of the old methods are preserved (test/argument_checking.jl passes unchanged); new tests in test/real_mul.jl cover mul! 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.
  • The 2D real plan's first call graph is now built for the half length, like the 1D plan; 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/FFTW vs FFTW 3.3.11 ESTIMATE):

class type cases FFTA speedup geomean (min–max) FFTA/FFTW before → after max bytes/exec before → after
1d/awkward Float32 29 1.04× (0.99–1.30) 7.73× → 7.60× 33793 KiB → 24576 KiB
1d/awkward Float64 29 1.03× (0.99–1.18) 6.17× → 6.05× 67585 KiB → 49152 KiB
1d/pow2 Float32 20 1.06× (0.97–1.45) 4.47× → 4.30× 16384 KiB → 0 KiB
1d/pow2 Float64 20 1.07× (0.98–1.29) 3.29× → 3.14× 32768 KiB → 0 KiB
1d/prime Float32 19 1.03× (1.00–1.19) 6.57× → 6.52× 118556 KiB → 98308 KiB
1d/prime Float64 19 1.03× (0.98–1.21) 6.25× → 6.11× 237112 KiB → 196620 KiB
1d/smooth Float32 24 1.05× (1.01–1.17) 10.47× → 10.10× 16279 KiB → 0 KiB
1d/smooth Float64 24 1.08× (1.01–1.18) 8.66× → 8.45× 32558 KiB → 0 KiB
2d Float32 9 2.05× (1.54–2.52) 13.23× → 6.51× 81973 KiB → 32 KiB
2d Float64 14 1.93× (1.10–2.71) 11.83× → 6.24× 163948 KiB → 6232 KiB
batched_dim1 Float32 6 1.11× (1.05–1.19) 6.39× → 5.74× 33034 KiB → 0 KiB
batched_dim1 Float64 9 1.15× (1.06–1.23) 4.90× → 4.26× 66061 KiB → 0 KiB
batched_dim2 Float32 6 1.13× (1.07–1.22) 3.54× → 3.15× 33034 KiB → 0 KiB
batched_dim2 Float64 6 1.10× (0.95–1.27) 3.05× → 2.76× 66061 KiB → 0 KiB

Worst / best individual cases:

  • 0.95× — rfft Float64 64×16384 dims=(2,): 28248.7 µs → 29838.5 µs (FFTW 16671.0 µs)
  • 0.97× — rfft Float32 65536 dims=(1,): 1320.9 µs → 1361.2 µs (FFTW 251.7 µs)
  • 0.98× — rfft Float64 262144 dims=(1,): 6829.4 µs → 6985.8 µs (FFTW 2628.5 µs)
  • 2.52× — rfft Float64 512×512 dims=(1, 2): 17474.8 µs → 6931.3 µs (FFTW 1988.9 µs)
  • 2.52× — rfft Float32 512×512 dims=(1, 2): 16156.1 µs → 6405.4 µs (FFTW 1474.3 µs)
  • 2.71× — rfft Float64 2048×2048 dims=(1, 2): 402560.4 µs → 148299.9 µs (FFTW 56682.1 µs)

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 rfft along dims=2 of wide 64×N matrices 1.2–1.3× slower than main (while dims=1 got 1.7× faster), and a back-to-back five-branch probe there attributed it to this PR: the old mapslices copy 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 the dims=2 result is bit-identical to mapslices). With the fix this PR is at parity with main on the strided shape and keeps the contiguous win: x86-64 back-to-back, rfft Float64 64×16384 dims=2 37.1 ms (main) / 43.6 (pre-fix) / 38.3 (fixed), 64×65536 290 / 389 / 276 ms, while dims=1 stays 3.8–3.9× faster than main; aarch64 shows the same picture (29.8 → 27.1 ms and 157 → 148 ms against main's 28.2 / 157). The 1.3–1.7× improvement of that shape over main arrives 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 rfft time 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, the dims path and 2D.

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

codecov Bot commented Aug 29, 2026

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 95.93023% with 7 lines in your changes missing coverage. Please review.
✅ Project coverage is 97.52%. Comparing base (7aeb327) to head (600d633).
⚠️ Report is 4 commits behind head on main.

Files with missing lines Patch % Lines
src/plan.jl 95.93% 7 Missing ⚠️
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.
📢 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.

…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.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant