Repository navigation
Conversation
Every SIMD planner wrapped the scalar BluesteinsAlgorithm around a SIMD inner FFT. SimdBluesteins is written against SimdVector instead, so there is one implementation for all four backends rather than four, and the scalar BluesteinsAlgorithm becomes an alias for the single-complex-per-vector case of it. The three pairwise multiply loops are the only SIMD part. They go through new SimdVector methods so that they land inside the backend's target feature, and use a new mul_complex_conjugated primitive: (a * b).conj() is a.conj() * b.conj(), so storing the precomputed side pre-conjugated turns every conjugated product into one operation. FCMA reaches it in the same instruction count as a plain complex multiply, by swapping vcmla_rot90 for vcmla_rot270. With one complex number per vector, interleaved data and a shuffle-based multiply, those loops lose to the plain scalar loop LLVM autovectorizes into deinterleaving loads: NEON f64 came out 3 to 5% slower at the larger sizes. The new FUSED_COMPLEX_MULTIPLY says whether a backend has native complex multiply instructions, and the loops take the scalar path when it is false and a vector holds a single complex number. That covers NEON, SSE and wasm f64, while FCMA f64 keeps the vector path it wins on. Measured on an Apple M1, full planned FFT against the scalar Bluestein's around the same inner FFT, at lengths 107 to 65579: NEON f32 1.6 to 5.2% faster, FCMA f32 2.0 to 12.5%, FCMA f64 1.8 to 13.0%, NEON f64 and the scalar planner unchanged.
The same treatment as SimdBluesteins, reusing its machinery: SimdRaders is written against SimdVector, so the four SIMD planners share one implementation instead of wrapping the scalar RadersAlgorithm, which becomes an alias for the single-complex-per-vector case. No new per-backend code was needed, only the pairwise multiply loops Bluestein's already added. Only the multiply between the two inner FFTs is SIMD. The input and output reorderings stay scalar because both are permuted copies, and without a gather instruction there is nothing for a vector to do. That caps the win: measured on an Apple M1 at lengths 101 to 64153, FCMA is 1 to 5% faster and NEON is unchanged, where Bluestein's reached 13%. The gather and scatter dominate the work outside the inner FFTs. The immutable and in-place paths differ only in where the inner FFT's scratch comes from, so they share one body now rather than repeating it as the scalar version did.
| // strength-reduced modular multiplies out of the hot loops. avx_raders.rs precomputes its | ||
| // output mapping for the same reason. This costs 4 * (len - 1) bytes, the same order as the | ||
| // twiddles avx_raders.rs and Bluestein's already store at this size. | ||
| permutation: Box<[u32]>, |
There was a problem hiding this comment.
The use of u32 here is nice because it cuts the amount of data in half, but it means users who plan a FFT larger than a u32::max don't get access to it. I did't check every planner, but I notice that SseSplanner for example doesn't account for this in its planning. So as-is, a user planning a prime FFT above u32::max will just crash or get data corruption.
We should add an assert to the constructor of SimdRaders that inner_fft.len() is u32::max or less, and add a step to the planner design_prime that switches to bluestein's if the fft size is larger than u32::max.
| /// This is the shared `SimdBluesteins` with a single complex number per vector, which is the scalar | ||
| /// case. Its pairwise multiply loops then compile to the same plain loops this algorithm had before | ||
| /// they were shared, so there is nothing left here but the name. | ||
| pub type BluesteinsAlgorithm<T> = SimdBluesteins<Complex<T>, T>; |
There was a problem hiding this comment.
A disadvantage of making this a type alias is that we lose the doc comment for the new() methods for this and raders. Run cargo doc --open and you'll see that the comments aren't visible.
I imagine it's rare these days that anyone would use anything but the planner, but in the interest of backwards compatibility I think we should prefer what's there. We can do that by making it a newtype instead that wraps SimdBluesteins etc.
| /// These are raw pointers so that the in-place case can pass the same one for `src` and `dest`. All | ||
| /// three must be valid for `count` elements, which the callers below check. | ||
| #[inline(always)] | ||
| unsafe fn mul_complex_vectors<V: SimdVector, const CONJUGATE: bool>( |
There was a problem hiding this comment.
We took this pointer approach for butterflies in the past, and got a lot of complaints, bug reports, even a CVE filed against it, because it's UB to have a mut and non-mut pointer pointing to the same memory.
Instead, I suggest having a single "srcdest" parameter which is an impl SimdComplexArrayMut<x>, and a multiplier slice. And it always reads from srcdest and writes back to it. If the methods on the trait do the same, we might be able to get rid of the inplace distiction, because inplace passes a mut slice directly, and out of place passes a DoubleBuf.
| /// those loops lose to the plain scalar loop LLVM autovectorizes into deinterleaving loads, so | ||
| /// they take the scalar path instead. FCMA's two-instruction multiply wins, and so does every | ||
| /// vector holding two complex numbers. | ||
| const FUSED_COMPLEX_MULTIPLY: bool; |
There was a problem hiding this comment.
Long shot: Could we benefit from hand-writing the deinterleave loop?
There was a problem hiding this comment.
I tried this and it showed promise for some situations, but it definitely isn't worth pursuing for this PR.
| /// See `SimdVector::FUSED_COMPLEX_MULTIPLY`. | ||
| #[inline(always)] | ||
| fn vectors_worthwhile<V: SimdVector>() -> bool { | ||
| V::COMPLEX_PER_VECTOR > 1 || V::FUSED_COMPLEX_MULTIPLY |
There was a problem hiding this comment.
While I was benchmarking my proposed handwritten deinterleave loop, my results suggest that this criteria for selecting algorithms doesn't capture what we need. Wasm simd f32 is faster to just write scalar code and let it autovectorize (and it generates almost identical assembly to my handwrriten deinterleave loop), but wasm 64, and sse 32 and 64 were all faster with the interleaved multiply loop.
I suggest another strategy: Create an enum for which algorithm to choose:
pub enum MulSliceAlgorithm {
Autovectorize,
Interleaved,
}
Then add an associated constant to SimdVector, and default it to interleaved:
pub trait SimdVector {
// ...
const MUL_SLICE_ALGORITHM : MulSliceAlgorithm = MulSliceAlgorithm::Interleaved;
// ...
}
Then in the mul slice code we can branch on that associated constant, instead of trying to infer it from other properties.
There was a problem hiding this comment.
Check out the branch i uploaded if you want to reproduce it for yourself: https://github.com/ejmahler/RustFFT/tree/mul-complex-slice-experiment
2 participants
Implements #190: one Bluestein's and one Rader's shared by all four SIMD backends, written against
SimdVector.SimdBluesteinsandSimdRadersinsrc/simd, which the SIMD planners now build instead of wrapping the scalar versions. The scalarBluesteinsAlgorithmandRadersAlgorithmbecome aliases for the single-complex-per-vector case, the wayRadixNalready is.mul_complex_conjugatedonSimdVector:(a * b).conj()isa.conj() * b.conj(), so storing the precomputed side pre-conjugated turns every conjugated product into one operation. On FCMA it costs the same as a plain multiply,vcmla_rot270instead ofvcmla_rot90.FUSED_COMPLEX_MULTIPLY, because with a single complex per vector and a shuffle-based multiply these loops lose to the scalar loop LLVM autovectorizes into deinterleaving loads. NEON f64 was 3 to 5% slower without it. NEON, SSE and wasm f64 take the scalar path, FCMA f64 keeps the vector path.boilerplate_simd_fft!besideboilerplate_fft!, for theFftimpl of theSimdVector-generic algorithms.Apple M1, planned FFT at lengths 107 to 65579 against the old path: Bluestein's is 1.6 to 5.2% faster on NEON f32 and 1.8 to 13.0% on FCMA, Rader's 1 to 5% on FCMA. NEON f64 and the scalar planner are unchanged. Rader's gains less because its reorderings are permuted copies with nothing for a vector to do.