Repository navigation
Oversample instead of Bluestein for non-nice sizes - #29
Open
oscardssmith wants to merge 7 commits into
Open
oscardssmith wants to merge 7 commits into
oscardssmith wants to merge 7 commits into
Conversation
Member
|
Benchmark of head Clear gains (both machines):
Regressions:
Questions:
|
Lengths with a prime factor above MAX_RADIX_PRIME take a type-2 NUFFT at the
points k·m/n (FFTJ's resampler): deconvolve, zero-pad to a smooth m ≥ 2n, one
transform, interpolate with a w-tap exponential-of-semicircle kernel. One
transform of length 2n where Bluestein takes two. FFTJ's settings: error
10 eps relative to ‖x‖ (w = 16 for F64, 8 for F32), σ = 2, γ = 0.97, weight
table up to a size (stored twice, for the re and im lanes of the interleaved
grid) and polynomial weights above it; the table limit is at measured parity
(F64 1 MiB, F32 4 MiB). Blocked interpolation: one load of the w interleaved
grid values per output, FMA folds, one tree reduction per vector of outputs.
- whole-length kernel (Resample), run-time plans (Resample over RChain; w
from the element type and the polynomial degree a field, so RResample is a
union of two concrete types), one long line on threads (ResampleMT: four-step
grid transform, deconvolution and output ranges over the tasks)
- composite kernel: the length-q DFTs of all s lanes in one batched pass
(RsBatch, as FFTJ's composite): lane-fastest grid rows in cache-sized chunks,
one batched chain, output rows as w FMAs of grid rows with broadcast weights
- own_buffers covers Resample (its buffers are per kernel)
Tiger Lake, resample / Bluestein (alternating rounds, min):
serial F64 1009 0.60, 10007 0.82, 100003 0.54
F32 1009 0.51, 10007 0.65, 100003 0.58
run-time F64 0.64, 0.94, 0.62 F32 0.50, 0.55, 0.63
4 threads F64 100003 0.60, 1000003 0.53 F32 0.63, 0.48
composite F64 1024·37 0.90, 64·1009 0.76, 128·37 0.79
F32 1024·37 0.64, 64·1009 0.81
Error vs a 128-bit reference: F64 2.5–4.7e-15 (Bluestein ~5e-16), F32
1.2–2.2e-7 (Bluestein 2.0–2.7e-7).
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
CompositeMT resampled lane by lane while the serial Composite runs the batched pass, so the two no longer agreed bitwise. The threaded kernel now runs the same pass, its tasks taking chunks of lanes (rs_chunk!) with copies of the RsBatch, at least one chunk per task where the lanes allow. The chain's vector width follows the lane count of a chunk, and vectors of one complex value round differently from wider ones, so a chunk runs its multiple of the vector's complex count and its remainder as separate groups; chunk starts are multiples of that count. Only the last s % CU lanes then take the narrow path, whatever the chunk size, and serial and threaded plans agree bitwise (all 46 cases of the threaded testset, 2–4 threads; F32 3·193·1021 on 4 threads is planned as a whole-length ResampleMT, as it was with Bluestein). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The reference was a direct O(n²) sum; the composite cases of 130462 points took about 11 minutes per call and four calls in the quick run, past a 30-minute timeout. FFTW.jl (already a test dependency) is an independent implementation at every size. The test of the AbstractFFTs generics without FFTW moves ahead of loading it; FFTW is imported, not used, so its exports do not meet FFTL's in the files included before fftw_compat.jl. Quick run: 9m15s (Tiger Lake), all passing. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Resampling has replaced it everywhere on the Float32/Float64 path (whole lengths, composite leftovers, run-time plans, threads), so Bluestein, BluesteinMT, bluestein_tables and the per-lane composite leftovers go. The helpers the resampler keeps are renamed after the grid: grid_length, grid_cost, GRID_RADIX, GRID_W, GRID_FOURSTEP_BYTES; use_resample_mt, resample_mt_dim and TUNING.resample_mt_min (same value); the fields Composite.rs, CompositeMT.rss and RPlan.rs. The composite split now sizes the leftover's grid as the resampler does (grid_length(T, 2q), was 2q - 1). Generic float types keep their exact Bluestein (generic.jl), and the KernelAbstractions extension its own. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The algorithms page describes the resampled DFT in place of Bluestein's algorithm (method, widths, weight tables, the measured comparison) and the batched pass of the composite kernel; the platforms table says per column which method a path uses (generic types and GPUs keep Bluestein). The Zen 2 timings of the composite kernel and the Rader paragraph, both measured or argued against Bluestein, are dropped. README and the other pages follow. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
The s resampled transforms of length q of a composite kernel (n = s·q) run as one batched pass (RsBatch) only from 32 lanes on. Below that every lane is one resampled transform (Resample), vectorised along the line. The batched pass vectorises across the lanes, so few lanes, or a lane count that is not a multiple of the vector width, leave most of the vector unused. Batched / lane by lane on an EPYC 9354 (AVX-512, one thread, q = 37, 1009, 10007), Float64 (Float32): s = 2 1.7–3.8 (2.1–7.3) s = 7 1.1–1.6 (1.6–2.9) s = 16 0.73–1.04 (0.84–1.31) s = 32 0.56–1.04 (0.87–0.96) s = 128 0.50–0.97 (0.66–1.13) The rule (RS_BATCH_LANES, leftover_kernel) depends on the length only, not on the number of threads, so Composite and CompositeMT take the same path and stay bitwise equal. The run-time pass plans take the same rule; RComposite is a union of four concrete types. Float64, one thread, same machine: 74 = 2·37 takes 0.26 of the time of the batched pass, 259 = 7·37 0.65, 3·193·1021 0.65; 1024·37, 64·1009 and 128·37 keep the batched pass and are unchanged. Not measured on AVX2 or aarch64. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
The grid of a single resampled transform (Resample, ResampleMT) lies between two margins. After the grid transform the margins hold copies of the other end of the grid, so the window of every output is one contiguous load and all outputs take the vector interpolation path. The outputs whose window wraps around the grid ends (7 per Float64 transform, 3 per Float32) otherwise need a scalar loop at about ten times the cost of the others: 250 of 559 ns at n = 37 and 290 of 1104 ns at n = 101. Time with margins over time without, Float64, AMD EPYC 9354, one thread, minima of interleaved runs (the same code against itself within 6 %): n = 37: 0.55, 101: 0.75, 127: 0.79, 1009: 0.96, 100003: 0.98, 101^3: 0.78. Against main: 0.81, 0.81, 0.77, 0.64, 0.80, 0.83, and 0.85 for 101^3 on 8 threads (1.17 without margins). The margin is a whole page, so the grid keeps the page offset of the work buffers: with a 128-byte margin the prime length 100003 took 1.03-1.05 of the time without margins. Measured on x86 only (512-bit and 256-bit vectors), not on aarch64. Co-Authored-By: Claude Fable 5.2 <noreply@anthropic.com>
oscardssmith
force-pushed
the
os/resample
branch
from
November 12, 2026 22:19
18582fa to
ea3b7f5
Compare
Member
Author
|
timings for short primes improved. |
Member
|
Follow-up on head Regressions that reproduce on the M1 (second, quieter run):
Please disregard the 8-thread Open on the EPYC:
How we measure (so you can reproduce on x64):
Quick check on x64: |
This branch has not been deployed
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.
Cells are main → PR (PR / main). Timing from AMD EPYC 9354 (32 cores) pinned to 1 numa node.
The rows with a short prime part (37 or 127:
74x259x127at 1 to 8 threads,fft 37,127x127and the small composite table) are fromea3b7f5againstmainatd63b21b. All other cells are from18582faagainstmainat103b9b6; that includes the 16- and 32-thread cells of74x259x127, which are not re-measured.Thread scaling
fft, 1D, Float64
fft, multidimensional, Float64
rfft, Float64
fft, 1D, Float32
fft, multidimensional, Float32
rfft, Float32
Small lengths and other paths, 1 thread
Float64
Float32
Small composite lengths (fft, Float64, 1 thread)