Repository navigation
Performance roadmap: closing the gap to FFTW (proposed PR sequence) #130
Description
Activity
This output was eager instead of living in a fork 😅 .
@dannys4 @wheeheee the big question here is whether you'd want to take the time to review these performance improvements or whether we should take this as a performance-minded fork of FFTA.jl. The full set of PR is there and either approach is fine, but it would just be good to know what you'd prefer on the policy here as we'd be looking to complete a high-performance GPL-free FFT using the current FFTA as the base.
Apologies, had some spare compute and tokens and took a shot at this 🙈 sorry for the notification overload and thanks for the great work around here 🤗
I see. this is a lot to digest, as you alluded to. I mentioned this in another PR at some point, but I'm not the happiest about introducing excessive LLM-generated code into the actual operational-side of the repo (though I could care less about benchmark code and, to a lesser extent, documentation). I really prefer having a pretty tight codebase, so I'm going to review it as if it were hand-written etc etc.
That being said, a few of the things you mentioned are somewhat higher priorities:
- Implemented properly, the precalculated bit-twiddling would also address inaccurate results due to inaccurate trigonometric recurrences #118.
- I've been lazy and perhaps sloppy with the function signatures
- I know that
@generatedfunctions could solve some of the performance issues; I seem to recall, however, that there was some discussion with @wheeheee about having typestable generated functions, but I'm not sure...
A few points that I'll raise on this "roadmap" issue:
- I have no idea what you mean by "Smooth-length Bluestein padding." I don't have a SP background so it could be me, but I have just never heard of this term and have no idea a priori what this means. Having the ability to tune the Bluestein cutoff, though, seems reasonable enough.
- I'm not sure that ND-transforms belong in this repo. Because they're "just" a composition of embarrassingly-parallel one-dimensional FFTs, and such parallelization could be done any number of ways---distributed, MPI, threads, you name it---with varying degrees of how these communicate, I'd be skeptical this is something that we should maintain. I am, however, extremely biased toward the design of the multi-dimensional FFT package I've worked in before, which basically just wrapped one-dimensional FFTs from different C++ backends for higher-dimensional FFT applications.
IMO items A & B definitely should be done in some form and IIRC was discussed previously; it was just a matter of allocating time to do them up properly. Haven't taken a close look at them but even if the LLM code is a bit too verbose / spaghetti it probably can still be hammered into something nicer and more maintainable.
C falls along the same lines, although it being better might be extremely problem size and arch-specific. I seem to recall seeing many FFTW plans mostly not being improved by going beyond a 2/4 split-radix kernel. Could be wrong though, probably should just benchmark more.
The rest, more or less what @dannys4 said, although I did implement a chunk of the ND transforms, rather rudimentarily (only single-threaded though; fine IMO for light usage, but if people intend to use FFTA for big problems...). Mainly because multithreading/multiprocessing+buffers gives me pause. There are good arguments for that to be in a separate package...perhaps there is a good way to cut up FFTA so some parts can be reused even there, but I'm digressing.
There was also the whole issue about friction with the AbstractFFTs API, and with so many more changes it might get even worse. Plus as long as FFTA and FFTW are both using the AbstractFFTs API they shouldn't be loaded together so there's not much point in resolving that little ambiguity.
I continued the Fable... fable and it added another 20 "PRs" for which I didn't open a PR. To be honest, at this point I don't even know what it did (all of it sits on main, +138 commits). Won't pretend this is manageable. On the benchmark world (which was its feedback look), it looks great though ( 30% faster to 100% slower, only 5% slower on DSP.jl test suite).
I would propose the following:
- I'd close all of the PRs that are there now
- I'll try to make sense of the changes (with the help of magician 🪄 @oscardssmith)
- we'll try to come up with a better plan to upstream them, that whatever this is
when we have a better idea of what's there, it may also make sense to discuss putting optimizations somewhere else or naming this FFTC or something.
what do you think? Again, apologies for the spam.
Again, apologies for the spam.
no it's fine! I'm happy that things are happening, even if having them happen stresses me out 😅 .
I have a few thoughts:
Store per-node twiddle tables (and Bluestein chirp/scratch) in the CallGraph at plan time;
I think this is probably the single most important proposition; if you had $100, I'd spend $99 on this before thinking about anything else. It really solves a lot of performance and accuracy issues.
mul! for real plans and a zero-allocation rfft/irfft path
The
rfftproblems, I think, are probably the easiest to fix, at least partially? I would hope they don't require substantial effort, but I'm not sure.use dims path for real transforms reusing
fft_along_dim!instead of mapslicesI would prefer keeping the mapslices; see NDFFT discussion below.
C being better might be extremely problem size and arch-specific. I seem to recall seeing many FFTW plans mostly not being improved by going beyond a 2/4 split-radix kernel.
This is something I have 100% uncertainty about. I suspect that, by using more mature compilers, FFTW might get some boosts from, e.g., tail-call optimizations. Doubly so in the meso-scale base-cases (e.g., 8, 9, 16, etc).
On ND transforms, I think I strongly lean against adding tooling to support more advanced strategies within this package. These operations can be applied pretty broadly and higher-dimensional FFT isn't really specific to FFTA, but the optimal design can be extremely specific for different use-cases (e.g., which slabs/pencils to work with in which order, what architecture you're using, etc etc). I think it's probably the job of FFTA to have a competent two- and three-dimensional FFT. I don't think, however, that it's the job of FFTA to have the optimal NDFFT method---if that's your goal, you're definitely better off using something else.
That being said, I'm 100% willing to ensure that FFTA supports any such package by ensuring everything is done atomically (as needed) and is threadsafe, etc etc. I hope I'm not being inflexible, but I'm just worried about the maintenance that the parallel FFT problem incurs. I'm happy to provide any feedback or modification of the library to accomodate NDFFTs within reason.
This is something I have 100% uncertainty about. I suspect that, by using more mature compilers, FFTW might get some boosts from, e.g., tail-call optimizations. Doubly so in the meso-scale base-cases (e.g., 8, 9, 16, etc).
Planned some FFTs with
FFTW.MEASURE, and yeah, turns out it uses radix 8/16/32 quite often. So...probably confused radix with something like vector size, on my old CPU that loves AVX2...
Also I don't think FFTW uses TCO. Haven't seen that optimization in Julia's LLVM IR before either.Anyway I'll just poke around a bit in #134 and #135 (basically A & C), those two look fairly close to being in a merge-able state.
I ran a broad FFTA-vs-FFTW sweep (PR #128) and profiled where the time goes, with the aim of closing the gap on the sizes where FFTA is far behind. This issue is to agree on the plan before opening PRs, since a few of the items touch the plan/callgraph structure.
Where FFTA stands today (aarch64 Neoverse-N1, Julia 1.12.6, FFTW 3.3.11 single-threaded,
ComplexF64, planned execution, FFTA / FFTW; full tables and plots in the PR):(geomean and range of FFTA/FFTW execution time;
Float32is 1.3–1.5× worse than this because FFTA gets no SIMD benefit; FFTWMEASUREis a further 2.2–2.8× faster thanESTIMATEat n ≥ 2^19.)What profiling says (details in
benchmark/ANALYSIS.mdin #129)singleton_params(asincospi) runs once per output row of the O(n²)DFTleaf, once perj1infft_composite!, and per recursion level in the radix-4/3 kernels. For n = 5 that is the whole cost of the transform (172 ns vs 37 ns with a table vs 23 ns for FFTW's codelet); for n = 1000 it is ~1800sincospiper execution.@generatedstraight-line 16/32/64-point base case runs at 1.0–1.6× FFTW and compiles in < 1 s.*(nomul!), 2D/oddrfftruns a full complex transform,rfftalongdimsgoes throughmapslices(10× FFTW).plan_rfft(::Vector{Float64}, ::Int)is a method ambiguity (FFTA'sregion::RegionTypesvs FFTW'sStridedArray).Proposed PRs, in order (each with before/after numbers from the suite)
CallGraphat plan time; kernels read from them. Adds aVector{Vector{T}}-style field per node; no API change. Expected: 2–5× on sizes with factors ≥ 5 and on primes, ~10–20% on pow2/pow3.mul!for real plans and a zero-allocationrfft/irfftpath;dimspath for real transforms reusingfft_along_dim!instead ofmapslices.Val{16/32/64},@generated, gated onT <: Union{ComplexF32,ComplexF64}so generic element types keep the current path), and codelets / radix butterflies for 5 and 7 (ref PFA implementation #105).Threads.@threadsover columns with one workspace per task).DEFAULT_BLUESTEIN_CUTOFF; Rader for primes with smooth n−1 (later).::RegionTypesannotation on theAbstractFFTsentry points and convert inside).A is the one that changes the plan structure (twiddle storage per node) and is a prerequisite for C and E, so I'd like to hear whether you're happy with twiddle tables living in
CallGraph(vs. e.g. a per-nodeNamedTupleor a separateVectorparallel tonodes) before I open it. Happy to adjust scope or order.