diff --git a/README.md b/README.md index e451eac0..3ce4bb75 100644 --- a/README.md +++ b/README.md @@ -84,8 +84,8 @@ yd = FFTL.fft(MtlArray(randn(ComplexF32, 1024, 1024))) **Fourier and discrete cosine transforms written in Julia, with wide support for number types and CPU and GPU execution.** - Complex, real-input and all eleven real-to-real kinds of FFTW, along any subset of - dimensions, for any length (the part of a length made of prime factors above 31 through - Bluestein's algorithm). + dimensions, for any length (the part of a length made of prime factors above 31 resampled onto + a smooth grid). - `Float32` and `Float64` on the fast path (SIMD codelets, threads). `BigFloat`, `Double64`, `Float128` and other floating-point types are transformed in their own precision by a generic path, accurate to a few ulps of the type. @@ -130,7 +130,8 @@ figures there are regenerated from the committed data on every documentation bui stride S is a chain started with stride S, the S neighbouring columns forming the vectorised index. No transposition. - **Lengths with prime factors above 31** split off their smooth part in one Cooley–Tukey step - and use Bluestein's algorithm only for the rest. + and resample only the rest onto a smooth grid (one transform of twice its length and a short + interpolation). - **Threads** split each dimension over slices or blocks of columns; long 1D transforms use the four-step algorithm. diff --git a/docs/AGENTS.md b/docs/AGENTS.md index 8f2fc521..2467f588 100644 --- a/docs/AGENTS.md +++ b/docs/AGENTS.md @@ -34,7 +34,7 @@ rules there apply to everything written here as well. Everything in this folder plans `RPlan`, `RRFFTPlan`, `RR2RPlan`), contribute. - **Explanation** (`src/explanation/`): the DFT and its conventions, FFTW and FFTL (history, papers, what FFTL takes from FFTW; written from papers and documentation only), algorithms (Stockham chains, - codelets, Bluestein, multidimensional without transposition, four-step, real and real-to-real + codelets, resampling, multidimensional without transposition, four-step, real and real-to-real reductions), implementation (type specialisation, SIMD, kernel cache, cache sizes, the two chain forms, why libfftl needs `RChain`), design decisions (what was decided and why, with the evidence), package size and load time, performance (measurement protocol). @@ -51,9 +51,9 @@ Design decisions page. These concern the manual itself: - **Interface in the text:** the manual writes the exported names without prefix; next to `using FFTL` it uses `import AbstractFFTs`, because `using AbstractFFTs` would make `fft` ambiguous. -- **Feature/platform table:** built from the README, the code and the libfftl header. The Bluestein - convolution-length rule is taken from the code (aarch64: smallest 5-smooth M ≥ 2n−1; x86: cost model, - 7-smooth M in [2n−1, 1.4·(2n−1)]); the FFTW-padded real layout of in-place r2c/c2r in libfftl comes +- **Feature/platform table:** built from the README, the code and the libfftl header. The grid-length + rule of resampled lengths is taken from the code (aarch64: smallest 5-smooth m ≥ 2n; x86: cost model, + 7-smooth m in [2n, 2.8·n]); the FFTW-padded real layout of in-place r2c/c2r in libfftl comes from the header. - **Figures:** CairoMakie in the tutorials. Literate runs with `execute = false` (its default for markdown), so the scripts execute as Documenter `@example` blocks during `makedocs`; a figure cell ends diff --git a/docs/src/explanation/algorithms.md b/docs/src/explanation/algorithms.md index ac848b16..67d8a0bf 100644 --- a/docs/src/explanation/algorithms.md +++ b/docs/src/explanation/algorithms.md @@ -19,8 +19,8 @@ The radices are chosen by a fixed rule (`src/factor.jl`). The rule is the result searches over the decompositions of many lengths on the development machines, not of a runtime search. It prefers radix 4 and radix 8 for powers of two. On aarch64 it merges a factor 3 with a 4 into a 12; the x86 rules differ and are in the same file. Every other prime up to 31 is a radix of -its own. Larger primes are not radices; the part of the length they make up goes to Bluestein's -algorithm (below). +its own. Larger primes are not radices; the part of the length they make up is resampled onto a +smooth grid (below). ## Codelets @@ -38,26 +38,37 @@ The generated code runs on complex SIMD vectors (`src/complexvec.jl`), with the the element type. On aarch64 the complex add-and-rotate and the complex multiply are single instructions (`FCADD`, `FCMLA`). -## Bluestein's algorithm - -Lengths with a prime factor above 31 go through Bluestein's identity. With ``jk = (j^2 + k^2 - (k-j)^2)/2``, - -```math -X_k = c_k \sum_{j} (x_j c_j)\, \overline{c_{k-j}}, \qquad c_j = e^{\sigma \pi i j^2/n}. -``` - -The sum is a cyclic convolution with a kernel that depends only on ``n``. Its length ``M`` -is a smooth number with ``M \ge 2n - 1``, so the convolution is two transforms of smooth lengths -and one pointwise product. The kernel's transform is computed once when the plan is made. - -Per line, Bluestein costs two transforms of length ``M`` and a pointwise product, so it is used -only for the part of a length that the smooth codelets cannot cover. On aarch64 the convolution -length is the smallest 5-smooth number with ``M \ge 2n - 1``. On x86 a cost model fitted to -measured chains picks the fastest 7-smooth length in ``[2n - 1,\, 1.4\,(2n - 1)]``: a length with -many radix 9, 5 or 3 stages transforms several times slower than a power of two a few percent -larger (`src/plan.jl`). - -On x86 a convolution of 8 MiB or more (``M \ge 2^{19}`` in `Float64`) runs as a +## Lengths with large prime factors: resampling + +Lengths with a prime factor above 31 are transformed by oversampled interpolation onto a smooth +grid (a type-2 nonuniform FFT evaluated at the points ``k\,m/n``; `src/resample.jl`): + +1. the input, divided by the Fourier transform of an interpolation kernel ``\psi``, is + zero-padded to a smooth length ``m \ge 2n``, +2. one transform of length ``m``, +3. output ``k`` is the sum of the ``w`` grid values around ``k\,m/n``, weighted by + ``\psi(z) = e^{\beta(\sqrt{1 - z^2} - 1)}`` (the exponential of a semicircle). + +The result is not exact: its error is about ``e^{-\pi w \sqrt{1 - n/m}}`` relative to ``\lVert x \rVert``. +FFTL takes the narrowest ``w`` (a multiple of 4) for an error of ``10\,\varepsilon``: 16 taps in +`Float64`, 8 in `Float32`. The weights of every output are tabulated up to a table size (1 MiB +in `Float64`, 4 MiB in `Float32`, where the two measured about even), above it they are +evaluated from piecewise polynomials. Generic number types (`BigFloat` and others, +`src/generic.jl`) keep Bluestein's algorithm, which is exact: there the width grows with the +precision and the interpolation would cost more than the transform it saves. + +Compared with Bluestein's algorithm, which needs two transforms of a length ``\ge 2n - 1``, one +transform is saved for ``n\,w`` multiply-adds. Measured against it on a Tiger Lake core (one +thread, before Bluestein was removed), the time ratio was 0.51–0.82 for primes from 1009 to +100003, at an error of 2–5·10⁻¹⁵ in `Float64` (Bluestein about 5·10⁻¹⁶) and 1–2·10⁻⁷ in +`Float32` (Bluestein about 2.5·10⁻⁷). + +The grid length is a smooth number. On aarch64 it is the smallest 5-smooth ``m \ge 2n``. On x86 +a cost model fitted to measured chains picks the fastest 7-smooth length in ``[2n,\, 2.8\,n]``: a +length with many radix 9, 5 or 3 stages transforms several times slower than a power of two a +few percent larger (`grid_length` in `src/plan.jl`). + +On x86 a grid of 8 MiB or more (``m \ge 2^{19}`` in `Float64`) is transformed by a [four-step](@ref "Long one-dimensional transforms: four-step") transform on the calling thread instead of a chain. A chain of that length no longer fits in cache and makes one pass over memory per stage; the four-step makes two. On a Zen 2 core this is 1.06× faster at ``2^{19}`` @@ -67,7 +78,7 @@ and 1.3× from ``2^{20}`` on. On the M1 the chain is faster at every length meas ### Lengths with a smooth part A length ``n = s\,q`` whose large prime factors make up ``q`` and whose small ones make up -``s > 1`` does not need a convolution of length ``2n``. Since ``s`` and ``q`` are coprime, one +``s > 1`` does not need a grid of length ``2n``. Since ``s`` and ``q`` are coprime, one Cooley–Tukey step with ``j = j_2 + q\,j_1`` and ``k = k_1 + s\,k_2`` splits the transform: ```math @@ -76,26 +87,30 @@ X_{k_1 + s k_2} = \sum_{j_2=0}^{q-1} e^{\sigma 2\pi i\, j_2 k_2/q}\; ``` 1. A chain of length ``s`` over the ``q`` interleaved lines ``j_2`` (the inner sum), -2. multiplication by the twiddle factors ``e^{\sigma 2\pi i\, j_2 k_1/n}``, folded into the first - chirp product of the next step, -3. ``s`` Bluestein transforms of length ``q`` (the outer sum), each written with stride ``s``. - -The arithmetic is about the same as one Bluestein transform of length ``n``: ``s`` convolutions -of length about ``2q`` are again about ``2n`` points. The gain is where the data lives. A -convolution of length ``2q`` fits in cache, one of length ``2n`` often does not. On a Zen 2 core -(one thread, `Float64`) ``262145 = 65 \cdot 4033`` takes 6.6 ms instead of 16.8 ms, and -``1048577 = 17 \cdot 61681`` takes 43 ms instead of 89 ms. +2. multiplication by the twiddle factors ``e^{\sigma 2\pi i\, j_2 k_1/n}``, folded into the + deconvolution of the next step, +3. ``s`` resampled transforms of length ``q`` (the outer sum), each written with stride ``s``. + +From 32 lanes on, the ``s`` transforms of the last step run as one batched pass (as in FFTJ): +the grid holds rows of lanes, lane fastest, so the transform of length ``m`` is one chain over +the lanes and every output row is ``w`` multiply-adds of whole grid rows. Lanes go in chunks +whose grid stays in cache. A grid of length about ``2q`` fits in cache where one of length +``2n`` often does not. + +With fewer than 32 lanes each lane is one resampled transform of its own, vectorised along the +line. A 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: on an EPYC 9354 (AVX-512, one +thread, ``q`` = 37, 1009, 10007) the batched pass took 1.7–3.8 times as long as the lanes one by +one at ``s = 2`` in `Float64` (2.1–7.3 in `Float32`), 1.1–1.6 (1.6–2.9) at ``s = 7``, 0.73–1.04 +(0.84–1.31) at ``s = 16`` and 0.56–1.04 (0.87–0.96) at ``s = 32`` (`RS_BATCH_LANES` in +`src/resample.jl`). The rule depends on the length only, not on the number of threads, so a +threaded plan computes the same values as a serial one. When ``q`` itself is large it is split again: ``q`` becomes the full power of its largest prime, -and the first step of length ``s = n/q`` is again a split transform (or a Bluestein transform -when ``s`` has no smooth part). ``3660001 = 73 \cdot 181 \cdot 277`` runs as Bluestein transforms of -length 13213 and 277. Each split costs one more pass over memory, so a leftover ``q`` stays one -Bluestein transform as long as its convolution takes at most 4 MiB. Below that, splitting -measured 1.15–1.5× slower; above it, 1.02–1.4× faster (`composite_split` in `src/plan.jl`). - -Rader's algorithm, which turns a prime length ``p`` into a cyclic convolution of length -``p - 1``, would halve the convolution where ``p - 1`` is smooth. Of 45 random lengths with a -large prime factor between ``2^{10}`` and ``2^{22}``, five qualify; FFTL does not implement it. +and the first step of length ``s = n/q`` is again a split transform (or a resampled transform +when ``s`` has no smooth part). Each split costs one more pass over memory, so a leftover ``q`` +stays one resampled transform as long as its grid takes at most 4 MiB (`composite_split` in +`src/plan.jl`; the threshold was measured with Bluestein's algorithm in place of resampling). ## Multidimensional transforms without transposition diff --git a/docs/src/explanation/design-decisions.md b/docs/src/explanation/design-decisions.md index 437bef8a..27f40a00 100644 --- a/docs/src/explanation/design-decisions.md +++ b/docs/src/explanation/design-decisions.md @@ -43,7 +43,7 @@ line, all buffers at fixed page offsets, time of the radix-16 sequence over the ## Rules for lengths are derived on some lengths and checked on others -**Decision (2026-11-04).** The radix choice and the Bluestein length are rules on the prime +**Decision (2026-11-04).** The radix choice and the grid length of a resampled transform are rules on the prime factorisation, not tables per length. A rule is derived from one part of the lengths (a brute-force search over all decompositions, or a cost model) and checked on lengths that were not used for the derivation, for example by parity of n or by a random split; only the checked result counts as @@ -56,23 +56,44 @@ not single points. the M1. A dense sweep of every length 2…1024 against FFTW MEASURE (one thread, Zen 2, 2046 cases) showed 100 `Float32` and 64 `Float64` lengths slower than 1.5×, the worst at 3.4×. -## Lengths with large prime factors: split off the smooth part, Bluestein on the rest +## Lengths with large prime factors: split off the smooth part, resample the rest **Decision (2026-11-07).** For n = s·q, with q the part made of primes above 31: one Cooley–Tukey -step, a chain of length s, twiddles, and s Bluestein transforms of length q (`Composite`). A leftover q -stays one Bluestein while its convolution is at most 4 MiB; above that, q is split into its prime -powers, recursively. Bluestein convolutions of 8 MiB and more run as a serial four-step algorithm on -x86 (1.3× from 2²⁰ on, Zen 2) but not on aarch64, where the chain is faster up to 2²² (M1 Max). With -threads, `CompositeMT` is used when s is at least the number of threads and is bitwise equal to the -serial kernel; otherwise the whole length is one `BluesteinMT`. - -**Why.** The convolutions of length about 2q stay in cache, and the arithmetic equals one Bluestein of -length n. On 45 non-smooth lengths (Zen 2) FFTL moved from ×0.64 to ×1.04 of the fastest other library -with one thread, and from ×0.77 to ×1.16 with 16 threads. +step, a chain of length s, twiddles, and s resampled transforms of length q (`Composite`): one batched +pass from 32 lanes on, lane by lane below (2026-11-12). A leftover q stays one resampled transform while its grid is at most 4 MiB; above that, +q is split into its prime powers, recursively. Grids of 8 MiB and more are transformed by a serial +four-step algorithm on x86 (1.3× from 2²⁰ on, Zen 2) but not on aarch64, where the chain is faster up +to 2²² (M1 Max). With threads, `CompositeMT` is used when s is at least the number of threads and is +bitwise equal to the serial kernel; otherwise the whole length is one `ResampleMT`. + +**Why.** The grids of length about 2q stay in cache. On 45 non-smooth lengths (Zen 2) FFTL moved from +×0.64 to ×1.04 of the fastest other library with one thread, and from ×0.77 to ×1.16 with 16 threads +(measured with Bluestein's algorithm on the leftover). Resampling needs one transform of a length of +at least 2q where Bluestein's algorithm needs two: it took 0.51–0.82 of Bluestein's time for primes +from 1009 to 100003 (Tiger Lake, one thread), at an error of 2–5·10⁻¹⁵ in `Float64` (Bluestein about +5·10⁻¹⁶). The batched pass vectorises across the lanes and loses with few of them: against the lanes one +by one it took 1.7–3.8 times as long at s = 2, 1.1–1.6 at s = 7, 0.73–1.04 at s = 16 and 0.56–1.04 at +s = 32 (`Float64`, q = 37, 1009, 10007, EPYC 9354, one thread; `Float32` 2.1–7.3, 1.6–2.9, 0.84–1.31, +0.87–0.96). The switch is on the lane count and not on the number of threads, which keeps serial and +threaded plans bitwise equal. + +The grid of a single resampled transform has a margin on each side that holds copies of the other end +(2026-11-12), so every output is interpolated from one contiguous load. Without it the outputs whose +window wraps around the grid ends (7 per `Float64` transform, 3 per `Float32`) took 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. With the +margin, time over the time without it: n = 37: 0.55, 101: 0.75, 127: 0.79, 1009: 0.96, 100003: 0.98, +101³: 0.78; against `main`: 0.81, 0.81, 0.77, 0.64, 0.80, 0.83, and 0.85 for 101³ on 8 threads (1.17 +before). `Float64`, EPYC 9354, minima of interleaved runs, the same code against itself within 6 %. +The margin is a whole page: a margin of 128 bytes moves the grid off the page offset of the work +buffers, and the prime length 100003 then took 1.03–1.05 of the time without margins. + +**Replaced.** Bluestein's algorithm on the leftover and `BluesteinMT` for a whole line on threads +(2026-11-07). Generic number types and the GPU passes keep Bluestein's algorithm. **Not done / not checked.** Rader's algorithm (5 of the 45 sweep lengths qualify). Both thresholds were measured on Zen 2 only, and the 4 MiB rule is checked on 12 lengths that also took part in its -derivation. +derivation. The 32-lane rule was measured on one AVX-512 machine only, not on AVX2 or aarch64, and the grid +margin on x86 only (512-bit and 256-bit vectors), not on aarch64. ## Two chain forms: compiled for short smooth lengths, run-time data for the rest @@ -144,7 +165,7 @@ sets it to empty, because `juliac` compiles nothing at run time. **Decision (2026-11-10).** A pass plan (`plan_pass`, `plan_r2c_pass`, `plan_c2r_pass`) built with `_runtime` (`plan_pass_runtime`, `plan_r2c_pass_runtime`, `plan_c2r_pass_runtime`) uses run-time kernels -only: RChain, Bluestein and composite transforms over RChain, the four-step algorithm over RChain. Its +only: RChain, resampled and composite transforms over RChain, the four-step algorithm over RChain. Its type depends on the element type and the number of dimensions, not on the length. The specialised plans are unchanged. @@ -158,8 +179,8 @@ and matches the direct DFT to 3·10⁻¹⁴ when called from C. Planning on 64 specialised and the run-time plans is within noise of each other on the same views. **Costs.** The run-time kernels dispatch per stage, which costs a little on short lines (the reason the -specialised plans stay the default). For real lengths with large prime factors, the run-time kernel is a -Bluestein transform, without the composite split of the specialised plan. +specialised plans stay the default). For real lengths with large prime factors, the run-time kernel is one +resampled transform, without the composite split of the specialised plan. ## The interface: FFTW.jl's names, FFTL's own functions, AbstractFFTs as a backend diff --git a/docs/src/explanation/implementation.md b/docs/src/explanation/implementation.md index 65e62587..ff008293 100644 --- a/docs/src/explanation/implementation.md +++ b/docs/src/explanation/implementation.md @@ -31,8 +31,8 @@ is how the KernelAbstractions kernels use the codelet generator without a second ## Kernel cache Kernels are cached by kind, element type, length, sign and radix sequence. The cache has 128 -entries and drops the least recently used one when it is full. A Bluestein kernel is two entries, -its chain and its tables. `clear_cache!` empties it. Within one plan, dimensions with the same +entries and drops the least recently used one when it is full. A resampling kernel is several +entries: its chain and its tables. `clear_cache!` empties it. Within one plan, dimensions with the same length and stride class share a kernel object. ## Cache sizes and blocks diff --git a/docs/src/explanation/performance.md b/docs/src/explanation/performance.md index 9b37eb9d..5ebeba40 100644 --- a/docs/src/explanation/performance.md +++ b/docs/src/explanation/performance.md @@ -50,6 +50,7 @@ application they happen at different times. - DCT-I and DST-I are computed through a real DFT of about twice the length, and DCT-IV of odd length likewise. The extra length shows in the measurements. - On the GPU, the real-to-real kinds are not tuned yet. -- Bluestein transforms a line through two transforms of length ``M ≥ 2q - 1``, where ``q`` is - at most the part of the length made of prime factors above 31 (the smooth part is split off, - see [Algorithms](@ref "Lengths with a smooth part")). A prime length pays this in full. +- A length with prime factors above 31 is resampled: one transform of a smooth length + ``m ≥ 2q`` and ``w`` multiply-adds per output, where ``q`` is at most the part of the length + made of those prime factors (the smooth part is split off, see + [Algorithms](@ref "Lengths with a smooth part")). A prime length pays this in full. diff --git a/docs/src/howto/contribute.md b/docs/src/howto/contribute.md index bea83aec..db89e284 100644 --- a/docs/src/howto/contribute.md +++ b/docs/src/howto/contribute.md @@ -72,7 +72,8 @@ The [Performance](@ref "Performance") page describes the protocol. | `src/codegen.jl` | the symbolic generator of the straight-line codelets | | `src/kernels.jl`, `src/complexvec.jl` | codelets on SIMD vectors, the complex arithmetic | | `src/factor.jl` | the radix sequence chosen for a length | -| `src/plan.jl`, `src/rchain.jl` | Stockham chains and Bluestein's algorithm; `rchain.jl` is the form that libfftl compiles | +| `src/plan.jl`, `src/rchain.jl` | Stockham chains and the composite kernel; `rchain.jl` is the form that libfftl compiles | +| `src/resample.jl` | lengths with prime factors above 31: resampling onto a smooth grid | | `src/kernel_methods.jl` | the stage kernels and short chains of the fixed codelet set as plain methods (see [Package size and load time](@ref "Package size and load time")) | | `src/threads.jl` | threads, lane blocks, the four-step algorithm | | `src/interface.jl` | the complex plans and the public functions | diff --git a/docs/src/howto/threads.md b/docs/src/howto/threads.md index 8a2a9707..bda19d33 100644 --- a/docs/src/howto/threads.md +++ b/docs/src/howto/threads.md @@ -31,9 +31,10 @@ number of values the transform touches: - A long line whose length has a large prime factor (at least ``2^{14}``) runs on threads as well. When the smooth part of the length is split off (``n = s\,q``, see [Algorithms](@ref "Lengths with a smooth part")) and ``s`` is at least the number of threads, - the ``s`` Bluestein transforms of length ``q`` are spread over the threads; the result is - bitwise the same as on one thread. Otherwise one Bluestein transform of the whole length - runs with its convolution transforms as four-step transforms on the threads. + the ``s`` resampled transforms of length ``q`` are spread over the threads in chunks of lanes; + the result is bitwise the same as on one thread. Otherwise one resampled transform of the whole + length runs with its grid transform as a four-step transform on the threads, and its + interpolation split over them. A 1D transform on a machine where some cores are busy is held up by the slowest piece, so fewer threads than cores can be the faster choice there. diff --git a/docs/src/howto/trim.md b/docs/src/howto/trim.md index 6cd9e1d5..a0bd35f4 100644 --- a/docs/src/howto/trim.md +++ b/docs/src/howto/trim.md @@ -17,7 +17,7 @@ For trimmed programs, FFTL has plans whose radix sequences and sizes are data: `T` is `Float32` or `Float64`; sizes, regions and kinds are vectors. All three take the `threads` keyword. They compute the same transforms with the same strategies as the normal plans -(Bluestein, four-step, fused planes); per stage they dispatch at run time, which costs a little +(resampling, four-step, fused planes); per stage they dispatch at run time, which costs a little speed. [libfftl](@ref "C library (libfftl)") is built from exactly these plans, so they are the tested route through `juliac --trim=safe`. diff --git a/docs/src/index.md b/docs/src/index.md index 0b1c4567..bb1a8fd8 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -32,7 +32,7 @@ shares. [Contribute](@ref "Contribute") says where help is most useful. - Complex, real-input and all eleven real-to-real kinds (DCT-I to DCT-IV, DST-I to DST-IV, Hartley, halfcomplex), along any set of dimensions of any rank. - Any length: smooth lengths through generated codelets, lengths with large prime factors - through Bluestein's algorithm, no padding to powers of two. + resampled onto a smooth grid, no padding to powers of two. - `BigFloat` and other floating-point types in their own precision, not converted to `Float64`. - GPUs through KernelAbstractions.jl, with the same functions. - Multithreading where the transform is large enough, and the four-step algorithm for one long line @@ -103,8 +103,8 @@ idct(C) ≈ r A plan is a Julia object that carries the radix sequence, the twiddle tables and the work buffers for one size and one set of dimensions. Radix sequences are chosen by a fixed rule, and the codelets they run are generated for the radix, element type and SIMD width actually used. -Long transforms run on several threads, and lengths with large prime factors go through -Bluestein's algorithm. The details are in [Algorithms](@ref "Algorithms"). +Long transforms run on several threads, and lengths with large prime factors are resampled onto a +smooth grid. The details are in [Algorithms](@ref "Algorithms"). ## Package status diff --git a/docs/src/reference/libfftl.md b/docs/src/reference/libfftl.md index e62cea93..c66640a4 100644 --- a/docs/src/reference/libfftl.md +++ b/docs/src/reference/libfftl.md @@ -25,7 +25,7 @@ library: All eleven real-to-real kinds (`FFTW_R2HC`, `FFTW_HC2R`, `FFTW_DHT`, `FFTW_REDFT00` … `FFTW_REDFT11`, `FFTW_RODFT00` … `FFTW_RODFT11`), any rank, any lengths (lengths with large -prime factors use Bluestein's algorithm), and every layout the manual describes: strides, +prime factors are resampled onto a smooth grid), and every layout the manual describes: strides, embedded arrays, vector loops, different input and output layouts, in place and out of place, interleaved and split complex arrays. diff --git a/docs/src/reference/platforms.md b/docs/src/reference/platforms.md index 25df9c92..00f5b4e0 100644 --- a/docs/src/reference/platforms.md +++ b/docs/src/reference/platforms.md @@ -17,7 +17,7 @@ available. | DHT, `R2HC`, `HC2R` | ✓ | ✓ | ◐ not tuned | ✓ | | batched 1D passes (`plan_pass`, `plan_r2r_pass`) | ✓ | – | – | – | | any rank and any subset of dimensions | ✓ | ✓ | ✓ | ✓ | -| any length (Bluestein above prime factor 31) | ✓ | ✓ | ✓ | ✓ | +| any length (prime factors above 31: resampling, or Bluestein where noted) | ✓ resampling | ✓ Bluestein | ✓ Bluestein | ✓ resampling | ## Platforms diff --git a/src/FFTL.jl b/src/FFTL.jl index 7bdd2ca0..11466e3d 100644 --- a/src/FFTL.jl +++ b/src/FFTL.jl @@ -30,6 +30,7 @@ include("codegen.jl") include("kernels.jl") include("factor.jl") include("plan.jl") +include("resample.jl") include("rchain.jl") include("kernel_methods.jl") include("threads.jl") diff --git a/src/factor.jl b/src/factor.jl index 869fd86f..9a6c1535 100644 --- a/src/factor.jl +++ b/src/factor.jl @@ -26,8 +26,8 @@ """ MAX_RADIX_PRIME -Largest prime factor handled by a codelet; lengths with larger prime factors use Bluestein's -algorithm. +Largest prime factor handled by a codelet; lengths with larger prime factors are resampled onto a +smooth grid (see `Resample`). """ const MAX_RADIX_PRIME = 31 diff --git a/src/interface.jl b/src/interface.jl index 7f68cb26..947969f9 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -72,7 +72,7 @@ function FFTLPlan{T}(sz::NTuple{N,Int}, region, σ::Int; threads::Integer = Thre function dim_kernel(d) fourstep_dim(sz, d, nt) && return FourStep{T}(sz[d], σ, nt) composite_mt_dim(T, sz, d, nt) && return CompositeMT{T}(sz[d], σ, nt) - bluestein_mt_dim(sz, d, nt) && return BluesteinMT{T}(sz[d], σ, nt) + resample_mt_dim(sz, d, nt) && return ResampleMT{T}(sz[d], σ, nt) S, _ = dim_extent(sz, d) lines_mt_dim(sz, d, nt) && return LinesMT{T}(sz[d], σ, S, nt) get!(() -> kernel(T, sz[d], σ, S), cache, (sz[d], S > LARGE_STRIDE || !smooth(sz[d]) ? S : 0)) @@ -99,11 +99,11 @@ function Base.show(io::IO, ::MIME"text/plain", p::FFTLPlan{T}) where {T} for (d, k) in zip(p.region, p.kernels) print(io, "\n dim $d: ") k isa RChain && print(io, "radices ", radix_tuple(k)) - k isa Bluestein && print(io, "Bluestein, convolution length ", k.fwd.n) - k isa Composite && print(io, "chain of length ", k.s, ", Bluestein of length ", k.q, ", convolution length ", k.blue.fwd.n) - k isa CompositeMT && print(io, "chain of length ", k.s, ", Bluestein of length ", k.q, " on threads, convolution length ", k.blues[1].fwd.n) - k isa BluesteinMT && print(io, "Bluestein on threads, convolution length ", k.fwd.n) + k isa Resample && print(io, "resampled, grid length ", k.fwd.n) + k isa Composite && print(io, "chain of length ", k.s, ", resampled DFT of length ", k.q, ", grid length ", grid_size(k.rs)) + k isa CompositeMT && print(io, "chain of length ", k.s, ", resampled DFT of length ", k.q, " on threads, grid length ", grid_size(k.rss[1])) k isa LinesMT && print(io, "one kernel per thread, blocks of ", k.B, " lines") + k isa ResampleMT && print(io, "resampled on threads, grid length ", k.fwd.n) k isa FourStep && print(io, "four-step $(k.n1) × $(k.n2), radices ", radix_tuple(k.k1), " and ", radix_tuple(k.k2)) end end @@ -228,7 +228,7 @@ end # These kernels represent a single long line and own their transform buffers. Reuse the plan's # serial workspace to gather/scatter only when that line is not contiguous. -function run_strided_dim!(py::Ptr{T}, px::Ptr{T}, p, k::Union{FourStep{T},BluesteinMT{T},CompositeMT{T}}, +function run_strided_dim!(py::Ptr{T}, px::Ptr{T}, p, k::Union{FourStep{T},ResampleMT{T},CompositeMT{T}}, d::Int, sy::NTuple{N,Int}, sx::NTuple{N,Int}) where {T,N} n = p.sz[d] if sx[d] == sy[d] == 1 diff --git a/src/pass.jl b/src/pass.jl index f7843c54..710c4fac 100644 --- a/src/pass.jl +++ b/src/pass.jl @@ -3,20 +3,23 @@ # Kernels of run-time pass plans (`plan_pass_runtime`). Their type depends on the element type only: # a smooth length runs on RChain, a leftover of large prime factors is a composite of an RChain and a -# Bluestein transform over RChain, and a single long line on several threads runs the four-step +# resampled transform over RChain, and a single long line on several threads runs the four-step # algorithm (RChain kernels). A library can therefore hold plans of any size in one container and # compile them ahead of time (juliac --trim). -const RComposite{T} = Composite{T,RChain{T},RBluestein{T}} -const RPassKern{T} = Union{RChain{T},RBluestein{T},RComposite{T},RFourStep{T}} +const RComposite{T} = Union{Composite{T,RChain{T},RsBatch{T,RChain{T},RTable{T}}}, + Composite{T,RChain{T},RsBatch{T,RChain{T},RPoly{T}}}, + Composite{T,RChain{T},Resample{T,RChain{T},RTable{T}}}, + Composite{T,RChain{T},Resample{T,RChain{T},RPoly{T}}}} +const RPassKern{T} = Union{RChain{T},RResample{T},RComposite{T},RFourStep{T}} function runtime_kernel(::Type{T}, n::Int, σ::Int, S::Int) where {T} smooth(n) && return RChain{T}(n, σ, radices(n, S; esize = sizeof(Complex{T}))) s, q = composite_split(T, n) - s == 1 && return Bluestein{T}(n, σ, RChain) + s == 1 && return Resample{T}(n, σ, RChain) chain = RChain{T}(s, σ, radices(s, S * q; esize = sizeof(Complex{T}))) - blue = Bluestein{T}(q, σ, RChain) + rs = leftover_kernel(T, q, σ, s, RChain) tw = cached(() -> Complex{T}.(unit_roots(σ, n)), CacheKey(:composite, T, n, σ, Int[]))::Vector{Complex{T}} - RComposite{T}(n, s, q, chain, blue, tw, Vector{Complex{T}}(undef, S * n)) + Composite{T,typeof(chain),typeof(rs)}(n, s, q, chain, rs, tw, Vector{Complex{T}}(undef, S * n)) end """ @@ -118,7 +121,7 @@ end """ plan_pass_runtime(v::StridedArray{<:Complex}, d; sign = -1, threads = 1) -The same transform as [`plan_pass`](@ref), with run-time kernels only (RChain, Bluestein and composite +The same transform as [`plan_pass`](@ref), with run-time kernels only (RChain, resampling and composite over RChain, four-step over RChain). The type of the plan is `PassPlan{T,N,RPassKern{T}}`, the same for every length, so that a library can hold plans of any size in one container and compile them ahead of time (`juliac --trim`). The kernels are the same as those of `plan_pass` for lengths with run-time diff --git a/src/plan.jl b/src/plan.jl index f5d1bcd2..30aadfd7 100644 --- a/src/plan.jl +++ b/src/plan.jl @@ -1,6 +1,6 @@ # Kernels for one dimension: a Stockham chain for lengths whose prime factors are -# ≤ MAX_RADIX_PRIME, Bluestein's algorithm on top of two chains otherwise. Both transform a batch -# of S interleaved sequences (element j of lane i at offset i + S·j), see kernels.jl. +# ≤ MAX_RADIX_PRIME, a resampled DFT on top of a chain otherwise (resample.jl). Both transform a +# batch of S interleaved sequences (element j of lane i at offset i + S·j), see kernels.jl. abstract type Kernel{T} end @@ -16,7 +16,7 @@ struct Chain{T,Rs,σ,K} <: Kernel{T} tw1::Vector{Complex{T}} # first stage: w^(kp) at (k-1)m + p end -# Read-only parts of kernels (chains, four-step twiddles, Bluestein chirp and filter) are shared +# Read-only parts of kernels (chains, four-step twiddles, resampling weights) are shared # between plans and threads: a second plan of a size costs a lookup. Bounded; when full, an # arbitrary entry is dropped. `clear_cache!()` frees everything. # Keys and values are concretely typed at every call site (kind, element type, sizes; values are @@ -297,24 +297,7 @@ end run_body(T, Rs, σ, K) end -""" - Bluestein{T} - -DFT of any length n as a cyclic convolution of length M ≥ 2n - 1 (M smooth): - X_k = c_k · Σ_j (x_j c_j) · conj(c_{k-j}), c_j = exp(σ·πi·j²/n). -""" -struct Bluestein{T,F,B} <: Kernel{T} - n::Int - chirp::Vector{Complex{T}} # c_j, j = 0…n-1 - kernel::Vector{Complex{T}} # DFT of the convolution kernel, scaled by 1/M - fwd::F - bwd::B - buf::Vector{Complex{T}} - work1::Vector{Complex{T}} - work2::Vector{Complex{T}} -end - -# 5-smooth numbers: the convolution lengths of the aarch64 rule. +# 5-smooth numbers: the grid lengths of the aarch64 rule. function next_smooth(n::Int) m = n while true @@ -324,28 +307,28 @@ function next_smooth(n::Int) end end -# Convolution length of Bluestein's algorithm: a 7-smooth M in [L, 1.4 L] with the lowest estimated -# time of a chain of length M, time = M·(c0 + Σ_stages w[radix]) in µs. The stage costs are fitted -# to measured chains (Zen 2, M = 70…300000, one thread); a length with many radix-9, -5 or -3 +# Grid length of a resampled DFT (see resample.jl): a 7-smooth M in [L, 1.4 L] with the lowest +# estimated time of a chain of length M, time = M·(c0 + Σ_stages w[radix]) in µs. The stage costs are +# fitted to measured chains (Zen 2, M = 70…300000, one thread); a length with many radix-9, -5 or -3 # stages transforms several times slower than a power of two a few percent larger (Float32: # 2025 = 3^4·5^2 takes 11 µs, 2048 takes 2.5 µs). The aarch64 stage costs have not been measured: # there the smallest 5-smooth M is kept. -const BLUESTEIN_RADIX = (2, 3, 4, 5, 6, 7, 8, 9, 10, 12) -# weights for BLUESTEIN_RADIX, then any other radix, then the constant per point -const BLUESTEIN_W = (Float64 = (0.0, 1.11e-3, 3.5e-4, 7.4e-4, 2.9e-4, 9.9e-4, 1.7e-4, 1.52e-3, 5.6e-4, 7.6e-4, 2.38e-3, 2.6e-4), - Float32 = (4.1e-4, 8.2e-4, 9.0e-5, 4.0e-4, 4.9e-4, 6.2e-4, 0.0, 1.08e-3, 3.4e-4, 1.6e-4, 1.6e-3, 5.1e-4)) +const GRID_RADIX = (2, 3, 4, 5, 6, 7, 8, 9, 10, 12) +# weights for GRID_RADIX, then any other radix, then the constant per point +const GRID_W = (Float64 = (0.0, 1.11e-3, 3.5e-4, 7.4e-4, 2.9e-4, 9.9e-4, 1.7e-4, 1.52e-3, 5.6e-4, 7.6e-4, 2.38e-3, 2.6e-4), + Float32 = (4.1e-4, 8.2e-4, 9.0e-5, 4.0e-4, 4.9e-4, 6.2e-4, 0.0, 1.08e-3, 3.4e-4, 1.6e-4, 1.6e-3, 5.1e-4)) -function bluestein_cost(::Type{T}, m::Int) where {T} - w = BLUESTEIN_W[T === Float32 ? :Float32 : :Float64] +function grid_cost(::Type{T}, m::Int) where {T} + w = GRID_W[T === Float32 ? :Float32 : :Float64] c = w[end] for r in radices(m) - i = findfirst(==(r), BLUESTEIN_RADIX) + i = findfirst(==(r), GRID_RADIX) c += i === nothing ? w[end-1] : w[i] end m * c end -function bluestein_length(::Type{T}, L::Int) where {T} +function grid_length(::Type{T}, L::Int) where {T} Sys.ARCH === :aarch64 && return next_smooth(L) hi = ceil(Int, 1.4L) best = 0; bestc = Inf @@ -358,7 +341,7 @@ function bluestein_length(::Type{T}, L::Int) where {T} p7 = p5 while p7 <= hi if p7 >= L - c = bluestein_cost(T, p7) + c = grid_cost(T, p7) c < bestc && (best = p7; bestc = c) end p7 *= 7 @@ -372,74 +355,11 @@ function bluestein_length(::Type{T}, L::Int) where {T} best end -# Convolution transforms of at least this many bytes run as a four-step on one thread: on Zen 2 a +# Grid transforms of at least this many bytes run as a four-step on one thread: on Zen 2 a # Stockham chain beyond the cache drops to 7 GFLOP/s, the four-step is 1.06× faster at 2^19 # (8 MiB, Float64) and 1.3× from 2^20 on, 0.8× at 3·2^17 and below. On the M1 Max the chain is # faster at every length up to 2^22. -const BLUESTEIN_FOURSTEP_BYTES = Sys.ARCH === :aarch64 ? typemax(Int) : 1 << 23 - -# C: the chain type of the two convolution transforms (Chain, or RChain for run-time plans) -function Bluestein{T}(n::Int, σ::Int, ::Type{C} = Chain) where {T,C} - M = bluestein_length(T, 2n - 1) - # one assignment each: the closure below captures fwd, and a variable assigned in two - # branches would be boxed (untyped, which juliac --trim cannot resolve) - fwd, bwd = convolution_transforms(T, M, C) - w1 = Vector{Complex{T}}(undef, M); w2 = Vector{Complex{T}}(undef, M) - chirp, b = cached(() -> bluestein_tables(T, n, σ, fwd, w1, w2), CacheKey(:bluestein, T, n, σ, Int[]) - )::Tuple{Vector{Complex{T}},Vector{Complex{T}}} - Bluestein{T,typeof(fwd),typeof(bwd)}(n, chirp, b, fwd, bwd, zeros(Complex{T}, M), w1, w2) -end - -# The forward and backward transforms of length M of a Bluestein kernel. The run-time plans keep -# their chains (one concrete type for juliac --trim). -convolution_transforms(::Type{T}, M::Int, ::Type{C}) where {T,C} = (C{T}(M, -1), C{T}(M, 1)) -function convolution_transforms(::Type{T}, M::Int, ::Type{Chain}) where {T} - if M * sizeof(Complex{T}) >= BLUESTEIN_FOURSTEP_BYTES && fourstep_split(M) > 0 - return (FourStep{T}(M, -1, 1), FourStep{T}(M, 1, 1)) - end - (RChain{T}(M, -1), RChain{T}(M, 1)) -end - -# chirp c_j = exp(σ iπ j²/n) = W_2n^(j² mod 2n) and the filter FFT_M(conj(c) wrapped) / M -function bluestein_tables(::Type{T}, n::Int, σ::Int, fwd, w1::Vector{Complex{T}}, w2::Vector{Complex{T}}) where {T} - M = length(w1) - W = unit_roots(σ, 2n) - c = Complex{T}[Complex{T}(W[mod(j * j, 2n)+1]) for j in 0:n-1] - b = zeros(Complex{T}, M) - b[1] = conj(c[1]) - for j in 1:n-1 - b[j+1] = b[M-j+1] = conj(c[j+1]) - end - GC.@preserve b w1 w2 run!(cptr(b), cptr(b), cptr(w1), cptr(w2), fwd, 1) - b ./= M - (c, b) -end - -# Lane by lane: gather, convolve, scatter. Element j of lane i at px + i + Sx·j, py + i + Sy·j. -# Uses the plan's own buffers, so one Bluestein kernel runs on one thread at a time. -function run_lanes!(py::Ptr{T}, px::Ptr{T}, ::Ptr{T}, ::Ptr{T}, P::Bluestein{T}, Sx::Int, Sy::Int, B::Int) where {T} - n = P.n; buf = P.buf; M = length(buf) - cx = Ptr{Complex{T}}(px); cy = Ptr{Complex{T}}(py) - GC.@preserve P begin - pb = cptr(buf); pw1 = cptr(P.work1); pw2 = cptr(P.work2) - for i in 1:B - @inbounds for j in 1:n - buf[j] = unsafe_load(cx, i + Sx * (j - 1)) * P.chirp[j] - end - fill!(view(buf, n+1:M), zero(Complex{T})) - run!(pb, pb, pw1, pw2, P.fwd, 1) - buf .*= P.kernel - run!(pb, pb, pw1, pw2, P.bwd, 1) - @inbounds for k in 1:n - unsafe_store!(cy, buf[k] * P.chirp[k], i + Sy * (k - 1)) - end - end - end - nothing -end - -run!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::Bluestein{T}, S::Int) where {T} = - run_lanes!(py, px, pw1, pw2, P, S, S, S) +const GRID_FOURSTEP_BYTES = Sys.ARCH === :aarch64 ? typemax(Int) : 1 << 23 """ Composite{T} @@ -447,45 +367,45 @@ run!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::Bluestein{T}, S::Int) DFT of a length n = s·q with s > 1 and q coprime to it (see `composite_split`), as one Cooley-Tukey step: with j = j2 + q·j1 and k = k1 + s·k2, a kernel of length s over j1 (a chain if s is smooth, else again a composite or -Bluestein kernel), the twiddles w_n^(j2·k1), and Bluestein's algorithm of length q over j2. The -convolutions have length about 2q instead of 2n and stay in cache. +resampling kernel), the twiddles w_n^(j2·k1), and a resampled DFT of length q over j2. The +transforms have length about 2q instead of 2n and stay in cache. """ struct Composite{T,C,B} <: Kernel{T} n::Int s::Int q::Int chain::C # length s, applied to the q lanes of every slab - blue::B # Bluestein of length q, with its own buffers + rs::B # length q over the lanes (RsBatch, or Resample for few lanes), with its own buffers tw::Vector{Complex{T}} # w_n^m, m = 0…n-1 z::Vector{Complex{T}} # output of the first stage, S·n values end function Composite{T}(n::Int, s::Int, q::Int, σ::Int, S::Int) where {T} chain = first_kernel(T, s, σ, S * q, S) - blue = Bluestein{T}(q, σ) + rs = leftover_kernel(T, q, σ, s) tw = cached(() -> Complex{T}.(unit_roots(σ, n)), CacheKey(:composite, T, n, σ, Int[]))::Vector{Complex{T}} - Composite{T,typeof(chain),typeof(blue)}(n, s, q, chain, blue, tw, Vector{Complex{T}}(undef, S * n)) + Composite{T,typeof(chain),typeof(rs)}(n, s, q, chain, rs, tw, Vector{Complex{T}}(undef, S * n)) end # Kernel of the first stage of a composite kernel: a chain over L lanes, or a kernel with buffers of -# its own (composite, Bluestein) that sees at most one lane block, or S lanes, at a time. +# its own (composite, resampling) that sees at most one lane block, or S lanes, at a time. first_kernel(::Type{T}, s::Int, σ::Int, L::Int, S::Int) where {T} = smooth(s) ? RChain{T}(s, σ, radices(s, L; esize = sizeof(Complex{T}))) : kernel(T, s, σ, max(composite_block(L, s), S)) -# Leftovers (prime factors above MAX_RADIX_PRIME) whose convolution takes at most this many bytes -# stay one Bluestein transform. Zen 2, one thread: splitting a leftover with an 8 MiB or larger -# convolution into its prime powers is 1.02–1.4× faster (591159 = 3·193·1021: 1.4×, 1194467 = 179·6673: 1.2×), -# splitting one of 128 KiB to 2.7 MiB is 1.15–1.5× slower (262145 = 65·37·109, 1588381 = 19·41·2039): -# each extra level is one more pass over memory. +# Leftovers (prime factors above MAX_RADIX_PRIME) whose grid takes at most this many bytes stay one +# resampled transform. Measured with Bluestein's algorithm (Zen 2, one thread): splitting a leftover +# with an 8 MiB or larger convolution into its prime powers is 1.02–1.4× faster (591159 = 3·193·1021: +# 1.4×, 1194467 = 179·6673: 1.2×), splitting one of 128 KiB to 2.7 MiB is 1.15–1.5× slower +# (262145 = 65·37·109, 1588381 = 19·41·2039): each extra level is one more pass over memory. const COMPOSITE_SPLIT_BYTES = 1 << 22 # The split n = s·q of the composite kernel: q the whole leftover r (the product of the prime -# factors above MAX_RADIX_PRIME) while its convolution is small, else the full power of the -# largest prime factor of r. (n, 1) for smooth n; s = 1 means a single Bluestein transform. +# factors above MAX_RADIX_PRIME) while its grid is small, else the full power of the largest prime +# factor of r. (n, 1) for smooth n; s = 1 means a single resampled transform. function composite_split(::Type{T}, n::Int) where {T} _, r = factorize(n) r == 1 && return (n, 1) - bluestein_length(T, 2r - 1) * sizeof(Complex{T}) <= COMPOSITE_SPLIT_BYTES && return (n ÷ r, r) + grid_length(T, 2r) * sizeof(Complex{T}) <= COMPOSITE_SPLIT_BYTES && return (n ÷ r, r) p = r; d = 37 while d * d <= r # largest prime factor of r (all its factors are ≥ 37) while r % d == 0 @@ -517,9 +437,7 @@ function run_lanes!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::Composi run_lanes!(pz + B * j2 * CT, px + Sx * j2 * CT, pw1, pw2, P.chain, Sx * q, B * q, B) end end - for i in 0:B-1, k1 in 0:P.s-1 - twiddled_bluestein!(py, pz, P.tw, P.blue, n, P.s, B, i, k1, Sy) - end + leftover!(py, pz, P.tw, P.rs, n, P.s, B, Sy) end nothing end @@ -536,31 +454,6 @@ function chain_blocks!(pz::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, C::Kern end composite_block(L, s) = min(L, max(LANE_ALIGN, LANE_ALIGN * fld(TUNING[].block_budget ÷ s, LANE_ALIGN))) -# Bluestein of length q on lane (i, k1) of the chain output z (element j2 at i + B·(j2 + q·k1)), -# after the twiddles w_n^(j2·k1); output element k2 at y + i + Sy·(k1 + s·k2). -function twiddled_bluestein!(py::Ptr{T}, pz::Ptr{T}, tw::Vector{Complex{T}}, bl::Bluestein{T}, n::Int, s::Int, - B::Int, i::Int, k1::Int, Sy::Int) where {T} - q = bl.n; buf = bl.buf; M = length(buf); chirp = bl.chirp - cz = Ptr{Complex{T}}(pz); cy = Ptr{Complex{T}}(py) - GC.@preserve bl tw begin - pb = cptr(buf) - zb = i + B * q * k1 + 1 - m = 0 # (j2·k1) mod n - @inbounds for j in 1:q - buf[j] = unsafe_load(cz, zb + B * (j - 1)) * tw[m+1] * chirp[j] - m += k1; m >= n && (m -= n) - end - fill!(view(buf, q+1:M), zero(Complex{T})) - run!(pb, pb, cptr(bl.work1), cptr(bl.work2), bl.fwd, 1) - buf .*= bl.kernel - run!(pb, pb, cptr(bl.work1), cptr(bl.work2), bl.bwd, 1) - @inbounds for k in 1:q - unsafe_store!(cy, buf[k] * chirp[k], i + Sy * (k1 + s * (k - 1)) + 1) - end - end - nothing -end - run!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::Composite{T}, S::Int) where {T} = run_lanes!(py, px, pw1, pw2, P, S, S, S) @@ -568,7 +461,7 @@ function kernel(::Type{T}, n::Int, σ::Int, S::Int = 1) where {T} smooth(n) && n <= SMALL_N && !has_lazy_factor(n) && return Chain{T}(n, σ, radices(n, S; esize = sizeof(Complex{T}))) smooth(n) && return RChain{T}(n, σ, radices(n, S; esize = sizeof(Complex{T}))) s, q = composite_split(T, n) - s == 1 && return Bluestein{T}(n, σ) + s == 1 && return Resample{T}(n, σ) Composite{T}(n, s, q, σ, S) end diff --git a/src/precompile.jl b/src/precompile.jl index 89407660..3997dd31 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -7,7 +7,7 @@ # # "lean" (default) 1024/1000/1009 in 1D, 16×16 and 8×8×8, both element types, all three # transforms, RPlan, plan_rfft and plan_pass. It covers the common -# codelets (radices 4, 5, 8, 9, 10 and 16, plus Bluestein) and the JIT +# codelets (radices 4, 5, 8, 9, 10 and 16, plus resampling) and the JIT # paths for both the compile-time and the run-time plans. # "broad" the measured shapes 10^5, 2^20, 256×256, 64×64×64 plus the radix # 12/6/3 lengths (360, 243), 128, and a four-step (threads = 4) case. diff --git a/src/r2r.jl b/src/r2r.jl index 7a68e3c6..7a6a2dff 100644 --- a/src/r2r.jl +++ b/src/r2r.jl @@ -58,7 +58,7 @@ function check_r2r_region(reg, N::Int) nothing end -own_buffers(K::R2RLine) = (i = K.inner; i isa RealLine ? i.inner isa Union{Bluestein,Composite} : i isa Union{Bluestein,Composite}) +own_buffers(K::R2RLine) = (i = K.inner; i isa RealLine ? i.inner isa Union{Resample,Composite} : i isa Union{Resample,Composite}) # exp(iπ(a j + b)/d) for j = 0…m-1 function cis_table(::Type{T}, m::Int, a::Int, b::Int, d::Int) where {T} @@ -70,7 +70,7 @@ function cis_table(::Type{T}, m::Int, a::Int, b::Int, d::Int) where {T} end # Inner kernel of the line for run-time plans (C = RChain): one of these types -const RInner{T} = Union{RealLine{T,RKern{T}},RChain{T},RBluestein{T}} +const RInner{T} = Union{RealLine{T,RKern{T}},RChain{T},RResample{T}} function R2RLine{T}(kind::R2RKind, n::Int, ::Type{C} = Chain) where {T,C} n >= 1 || throw(ArgumentError("length must be positive")) @@ -505,7 +505,7 @@ function half_dft(::Type{T}, sz, d, kind, z::Vector{Complex{T}}, nt) where {T} HalfDFT{T,typeof(cplan)}(S, outer, n, kind == RODFT10, kind == DHT, cplan, zz, w, cis_table(T, M + 1, -1, 0, 2n)) end -# One line kernel per task. A kernel with buffers of its own (Bluestein, composite) cannot run on +# One line kernel per task. A kernel with buffers of its own (resampling, composite) cannot run on # two tasks at once, so each task gets a copy; the other kernels are shared. task_lines(K::R2RLine{T}, threads::Integer) where {T} = own_buffers(K) && threads > 1 ? [K; [R2RLine{T}(K.kind, K.n) for _ in 2:threads]] : [K] diff --git a/src/real.jl b/src/real.jl index 399f8402..1a6187c1 100644 --- a/src/real.jl +++ b/src/real.jl @@ -24,10 +24,10 @@ struct RealLine{T,K<:Kernel{T}} w::Vector{Complex{T}} # u_k, k = 0…n÷2 (even n), see above end -# Run-time kernels (RChain, or Bluestein over RChain) for the run-time plans (rreal.jl). -const RKern{T} = Union{RChain{T},RBluestein{T}} +# Run-time kernels (RChain, or resampling over RChain) for the run-time plans (rreal.jl). +const RKern{T} = Union{RChain{T},RResample{T}} rkernel(::Type{T}, n::Int, σ::Int) where {T} = - (smooth(n) ? RChain{T}(n, σ, radices(n, 1; esize = sizeof(Complex{T}))) : Bluestein{T}(n, σ, RChain))::RKern{T} + (smooth(n) ? RChain{T}(n, σ, radices(n, 1; esize = sizeof(Complex{T}))) : Resample{T}(n, σ, RChain))::RKern{T} ckernel(::Type{T}, n::Int, σ::Int, ::Type{Chain}) where {T} = kernel(T, n, σ) ckernel(::Type{T}, n::Int, σ::Int, ::Type{RChain}) where {T} = rkernel(T, n, σ) @@ -169,7 +169,7 @@ line_fourstep(p) = nothing # Half-length DFT of one line on several threads, by the kernel type. run_half!(py::Ptr{T}, px::Ptr{T}, K::FourStep{T}, nt::Int) where {T} = run_fourstep!(py, px, K, nt) run_half!(py::Ptr{T}, px::Ptr{T}, K::CompositeMT{T}, nt::Int) where {T} = run_composite_mt!(py, px, K, nt) -run_half!(py::Ptr{T}, px::Ptr{T}, K::BluesteinMT{T}, nt::Int) where {T} = run_bluestein_mt!(py, px, K, nt) +run_half!(py::Ptr{T}, px::Ptr{T}, K::ResampleMT{T}, nt::Int) where {T} = run_resample_mt!(py, px, K, nt) # Lines along dimension d of an array of size `sz` (real side), as (lane stride S, line offsets # in elements of the real side and of the complex side). @@ -238,14 +238,14 @@ end sortregion(region, N) = Tuple(normregion(region, N)) # One real line of length n on several tasks: the half-length DFT on several threads (four-step for -# a smooth half length, composite or Bluestein for the others). The extra parallel pass for the +# a smooth half length, composite or resampling for the others). The extra parallel pass for the # pairs needs about 256 KB per line for smooth lengths and 128 KB for non-smooth ones (M1, 8 threads, # gain of the threaded over the serial line: Float32 smooth 0.7 at 128 KB, 1.5 at 256 KB; non-smooth # 1.05 at 120-130 KB, 1.8 at 147 KB; Float64 non-smooth 1.5 at 240 KB). use_real_threads(::Type{T}, n, nt) where {T} = (use_fourstep(inner_length(n), nt) && n * sizeof(T) >= TUNING[].real_fourstep_bytes) || - ((use_composite_mt(T, inner_length(n), nt) || use_bluestein_mt(inner_length(n), nt)) && - n * sizeof(T) >= TUNING[].real_bluestein_bytes) + ((use_composite_mt(T, inner_length(n), nt) || use_resample_mt(inner_length(n), nt)) && + n * sizeof(T) >= TUNING[].real_resample_bytes) # Odd n: the full-length DFT on several threads, in the line's buffers. Forward: the real line is # the real part of a complex line; backward: the half spectrum is completed by conjugate symmetry. @@ -285,7 +285,7 @@ inner_length(n) = iseven(n) ? n ÷ 2 : n function half_threads(::Type{T}, h, σ, nt) where {T} use_fourstep(h, nt) && return FourStep{T}(h, σ, nt) use_composite_mt(T, h, nt) && return CompositeMT{T}(h, σ, nt) - use_bluestein_mt(h, nt) && return BluesteinMT{T}(h, σ, nt) + use_resample_mt(h, nt) && return ResampleMT{T}(h, σ, nt) nothing end @@ -342,13 +342,13 @@ function Base.show(io::IO, p::BRFFTLPlan{T}) where {T} print(io, "FFTL c2r plan for $(join(p.sz, "×")) Complex{$T} → length $(p.d), dims $(p.region)") end -# Whether a line kernel uses work buffers of its own (a Bluestein kernel holds them in the kernel, +# Whether a line kernel uses work buffers of its own (a resampling kernel holds them in the kernel, # not in the task's buffer set): such lines run on one task. -own_buffers(k) = k isa Union{Bluestein,Composite} -own_buffers(K::RealLine) = K.inner isa Union{Bluestein,Composite} +own_buffers(k) = k isa Union{Resample,Composite} +own_buffers(K::RealLine) = K.inner isa Union{Resample,Composite} line_tasks(K, nt) = own_buffers(K) ? 1 : nt -# One line kernel per task. A kernel with buffers of its own (Bluestein, composite) cannot run on two +# One line kernel per task. A kernel with buffers of its own (resampling, composite) cannot run on two # tasks at once, so each task gets a copy; the other kernels are shared. # The copies have the kernel type of the line (C = RChain for the run-time plans). task_reallines(line::RealLine{T}, n::Int, σ::Int, nt::Int, ::Type{C} = Chain) where {T,C} = diff --git a/src/realpass.jl b/src/realpass.jl index 4f209a6d..0b5748e1 100644 --- a/src/realpass.jl +++ b/src/realpass.jl @@ -80,7 +80,7 @@ end """ plan_r2c_pass_runtime(v::StridedArray{<:Real}, d; threads = 1) -The same transform as [`plan_r2c_pass`](@ref), with run-time kernels (RChain, and Bluestein over RChain +The same transform as [`plan_r2c_pass`](@ref), with run-time kernels (RChain, and resampling over RChain for lengths with large prime factors). The plan type is `RealPassPlan{T,N,RealLine{T,RKern{T}}}`, the same for every length, so a library can compile it ahead of time (`juliac --trim`). """ diff --git a/src/resample.jl b/src/resample.jl new file mode 100644 index 00000000..42a21498 --- /dev/null +++ b/src/resample.jl @@ -0,0 +1,520 @@ +""" + Resample{T} + +DFT of any length n by oversampled interpolation (a type-2 NUFFT at the points k·m/n): the input, +divided by the Fourier transform of the kernel, is zero-padded to a smooth length m ≥ 2n and +transformed; output k is interpolated from the w grid values around k·m/n with the +exponential-of-semicircle kernel ψ(z) = exp(β(√(1-z²)-1)). One transform of length 2n where +Bluestein's algorithm takes two; the error is about exp(-πw√(1-n/m)) relative to ‖x‖. +""" +struct Resample{T,F,WS} <: Kernel{T} + n::Int + corr::Vector{T} # 1/Φ̂(j'/m) for input j (j' = j, or j - n for the upper half) + ws::WS # interpolation weights (RTable or RPoly) + fwd::F # length m, sign σ + buf::Vector{Complex{T}} # the grid between two margins (resample_grid) + work1::Vector{Complex{T}} + work2::Vector{Complex{T}} +end + +# Output k reads grid values l, …, l+w-1 (mod m) with l = q - h + 1, k·m = q·n + r; its weights +# depend only on δ = r/n. Either a table of the weights of every output, each stored twice for +# the re and im lanes of the interleaved grid (no shuffles to duplicate them), or P polynomial +# pieces of degree d in δ (Horner across the taps, then one duplication), coefficient c of tap i +# on piece p at C[(p(d+1) + c)w + i + 1]. Both are concrete per T (the run-time plans hold them). +struct RTable{T} + wt::Vector{T} + l0::Vector{Int32} +end +struct RPoly{T} + C::Vector{T} + P::Int + d::Int +end + +resample_tol(::Type{T}) where {T} = 10 * Float64(eps(T)) + +# Narrowest multiple of 4 meeting the tolerance at m ≥ 2n; a constant per T (16 for Float64, 8 for +# Float32), so the kernels take it from the element type. +resample_width(::Type{T}) where {T} = 4 * cld(ceil(Int, log(1 / resample_tol(T)) / (π * sqrt(0.5))), 4) + +# Largest weight table; above it the polynomial weights (no table streamed through the cache). At +# parity on Tiger Lake: F64 near 1 MiB, F32 (8 taps, less Horner work per byte saved) near 4 MiB. +resample_table_bytes(::Type{T}) where {T} = sizeof(T) >= 8 ? 1 << 20 : 1 << 22 + +# Margin on each side of a grid, in complex values. After the transform the margins hold copies of +# the other end of the grid (wrap_grid!), so the window of every output, w values from an index in +# -w/2+1 … m-w/2, is one contiguous load. Only w/2 values per side are used; the margin is a whole +# page so that the grid keeps the page offset of the two work buffers. With a 128-byte margin the +# prime length 100003 took 1.03–1.05 of the time without margins, with a page 0.98 (Float64, EPYC +# 9354, one thread). +resample_pad(::Type{T}) where {T} = 4096 ÷ sizeof(Complex{T}) +resample_grid(::Type{T}, m::Int) where {T} = Vector{Complex{T}}(undef, m + 2 * resample_pad(T)) +grid_ptr(buf::Vector{Complex{T}}) where {T} = cptr(buf) + resample_pad(T) * sizeof(Complex{T}) + +# Fill the margins of a transformed grid of length m ≥ w/2. +@inline function wrap_grid!(buf::Vector{Complex{T}}, m::Int) where {T} + H = resample_pad(T); h = resample_width(T) ÷ 2 + @inbounds for t in 1:h + buf[H-h+t] = buf[H+m-h+t] + buf[H+m+t] = buf[H+t] + end + nothing +end + +# C: the transform's chain type (Chain, or RChain for the run-time plans) +function Resample{T}(n::Int, σ::Int, ::Type{C} = Chain) where {T,C} + m, _, corr, ws = resample_parts(T, n) + fwd = resample_transform(T, m, σ, C) + Resample{T,typeof(fwd),typeof(ws)}(n, corr, ws, fwd, resample_grid(T, m), + Vector{Complex{T}}(undef, m), Vector{Complex{T}}(undef, m)) +end + +resample_transform(::Type{T}, m::Int, σ::Int, ::Type{C}) where {T,C} = C{T}(m, σ) +resample_transform(::Type{T}, m::Int, σ::Int, ::Type{Chain}) where {T} = + m * sizeof(Complex{T}) >= GRID_FOURSTEP_BYTES && fourstep_split(m) > 0 ? FourStep{T}(m, σ, 1) : RChain{T}(m, σ) + +# Grid length, width, deconvolution factors and weights of a length-n resampling (shared by plans). +function resample_parts(::Type{T}, n::Int, m::Int = grid_length(T, 2n)) where {T} + w = resample_width(T) + corr = cached(() -> resample_corr(T, n, m, w), CacheKey(:resample_corr, T, n, 0, [m, w]))::Vector{T} + ws = n * 2w * sizeof(T) <= resample_table_bytes(T) ? + cached(() -> resample_table(T, n, m, w), CacheKey(:resample_table, T, n, 0, [m, w]))::RTable{T} : + cached(() -> resample_poly(T, n, m, w), CacheKey(:resample_poly, T, n, 0, [m, w]))::RPoly{T} + (m, w, corr, ws) +end + +""" + RsBatch{T} + +The length-q DFTs of the s lanes of a composite kernel in one pass (FFTJ's batched resampling): +the grid holds rows of C lanes, lane fastest, so the transform is one chain batched over the +lanes and every output row is w FMAs of whole grid rows with broadcast weights (no reductions). +Lanes in chunks of C, so that a chunk's grid stays in cache. +""" +struct RsBatch{T,F,WS} + n::Int + m::Int + C::Int + corr::Vector{T} + ws::WS + fwd::F # length m over C lanes, sign σ + g::Vector{Complex{T}} # grid, m·C + work1::Vector{Complex{T}} + work2::Vector{Complex{T}} +end + +# grid bytes of one lane chunk (as FFTJ's chunk budget) +const RS_CHUNK_BYTES = 1 << 18 + +function RsBatch{T}(q::Int, σ::Int, s::Int) where {T} + m, _, corr, ws = resample_parts(T, q) + CU = ncomplex(first(vec_ladder(T))) + C = min(s, max(CU, CU * fld(RS_CHUNK_BYTES ÷ (m * sizeof(Complex{T})), CU))) + fwd = RChain{T}(m, σ, radices(m, C; esize = sizeof(Complex{T}))) + RsBatch{T,typeof(fwd),typeof(ws)}(q, m, C, corr, ws, fwd, Vector{Complex{T}}(undef, m * C), + Vector{Complex{T}}(undef, m * C), Vector{Complex{T}}(undef, m * C)) +end + +es_kernel(z, β) = abs(z) >= 1 ? zero(z) : exp(β * (sqrt(1 - z^2) - 1)) +resample_beta(n, m, w) = 0.97π * w * (1 - n / 2m) + +# Gauss-Legendre nodes and weights on [-1, 1] (Newton on the three-term recurrence) +function gauss_legendre(nq::Int) + z = zeros(nq); wq = zeros(nq) + for i in 1:nq + x = cospi((i - 0.25) / (nq + 0.5)) + dp = 1.0 + for _ in 1:20 + p0, p1 = 1.0, x + for k in 2:nq + p0, p1 = p1, ((2k - 1) * x * p1 - (k - 1) * p0) / k + end + dp = nq * (x * p1 - p0) / (x^2 - 1) + dx = p1 / dp + x -= dx + abs(dx) < 4eps() && break + end + z[i] = x + wq[i] = 2 / ((1 - x^2) * dp^2) + end + z, wq +end + +# Φ̂(ν) = w·∫₀¹ψ(z)cos(πνwz)dz by the positive half of a symmetric Gauss-Legendre rule +function resample_corr(::Type{T}, n::Int, m::Int, w::Int) where {T} + β = resample_beta(n, m, w) + zq, wq = gauss_legendre(2w + 20) + hp = cld(n, 2) + corr = Vector{T}(undef, n) + for j in 0:n-1 + ν = (j < hp ? j : j - n) / m + Φ = 0.0 + for q in eachindex(zq) + zq[q] > 0 && (Φ += wq[q] * es_kernel(zq[q], β) * cospi(ν * w * zq[q])) + end + corr[j+1] = T(1 / (w * Φ)) + end + corr +end + +function resample_table(::Type{T}, n::Int, m::Int, w::Int) where {T} + β = resample_beta(n, m, w) + h = w ÷ 2 + wt = Vector{T}(undef, n * 2w) + l0 = Vector{Int32}(undef, n) + for k in 0:n-1 + q, r = divrem(widemul(k, m), n) + l0[k+1] = Int32(q - h + 1) + δ = r / n + for i in 0:w-1 + wt[k*2w+2i+1] = wt[k*2w+2i+2] = T(es_kernel((δ + h - 1 - i) / h, β)) + end + end + RTable{T}(wt, l0) +end + +# Chebyshev interpolants of each tap's weight on each piece, in monomials of the local variable +# t ∈ [-1, 1]; the degree grows until the fit error is ≤ tol/16 or stops halving (coefficient +# rounding floors it). +function resample_poly(::Type{T}, n::Int, m::Int, w::Int) where {T} + β = resample_beta(n, m, w) + P = sizeof(T) >= 8 ? 32 : 64 + tolw = resample_tol(T) / 16 + prev, preverr = resample_fit(T, β, w, P, 1) + d = 1 + for dd in 2:30 + C, err = resample_fit(T, β, w, P, dd) + if err <= tolw + prev, d = C, dd + break + end + err > preverr / 2 && break + prev, preverr, d = C, err, dd + end + RPoly{T}(prev, P, d) +end + +function resample_fit(::Type{T}, β, w::Int, P::Int, d::Int) where {T} + h = w ÷ 2 + tj = [cospi((j + 0.5) / (d + 1)) for j in 0:d] + C = zeros(T, P * (d + 1) * w) + mono = zeros(d + 1) + err = 0.0 + for p in 0:P-1, i in 0:w-1 + f(t) = es_kernel(((p + (t + 1) / 2) / P + h - 1 - i) / h, β) + fj = f.(tj) + fill!(mono, 0) + Tkm, Tk = zeros(d + 1), zeros(d + 1) # T₀, T₁ in monomials + Tkm[1] = 1 + Tk[2] = 1 + for k in 0:d + a = sum(fj[j+1] * cospi(k * (j + 0.5) / (d + 1)) for j in 0:d) * 2 / (d + 1) + k == 0 && (a /= 2) + mono .+= a .* (k == 0 ? Tkm : Tk) + if k >= 1 # T_{k+1} = 2t·T_k - T_{k-1} + Tn = -Tkm + Tn[2:end] .+= 2 .* Tk[1:end-1] + Tkm, Tk = Tk, Tn + end + end + for c in 0:d + mono[c+1] = C[(p*(d+1)+c)*w+i+1] = T(mono[c+1]) + end + for t in range(-1.0, 1.0; length = 4d + 4) + err = max(err, abs(evalpoly(t, mono) - f(t))) + end + end + C, err +end + +# Grid start (0-based, not reduced mod m) and weights of output k. The table needs no state; the +# polynomials step st = (q, r), k·m = q·n + r, by (mq, mr) = divrem(m, n) per output, from output k. +rstart(::RTable, k, m, n) = nothing +rstart(::RPoly, k, m, n) = divrem(k * m, n) +@inline radvance(::RTable, st, mq, mr, n) = st +@inline function radvance(::RPoly, (q, r), mq, mr, n) + r += mr + c = r >= n + (q + mq + c, r - ifelse(c, n, 0)) +end +@inline rweights(ws::RTable{T}, k, st, n, ::Val{w}) where {T,w} = + @inbounds (Int(ws.l0[k+1]), vload(Vec{2w,T}, pointer(ws.wt) + k * 2w * sizeof(T))) +@inline function rweights(ws::RPoly{T}, k, (q, r), n, ::Val{w}) where {T,w} + d = ws.d + s = r * (ws.P / n) # δ·P ∈ [0, P) + pc = min(unsafe_trunc(Int, s), ws.P - 1) + t = T(2 * (s - pc) - 1) + pC = pointer(ws.C) + pc * (d + 1) * w * sizeof(T) + wv = vload(Vec{w,T}, pC + d * w * sizeof(T)) + for c in d-1:-1:0 + wv = muladd(wv, t, vload(Vec{w,T}, pC + c * w * sizeof(T))) + end + (q - w ÷ 2 + 1, shufflevector(wv, Val(ntuple(l -> (l - 1) ÷ 2, Val(2w))))) +end + +# Lane by lane: deconvolve and zero-pad, transform, interpolate. Element j of lane i at +# px + i + Sx·j, py + i + Sy·j. Uses the kernel's own buffers (one thread at a time). +function run_lanes!(py::Ptr{T}, px::Ptr{T}, ::Ptr{T}, ::Ptr{T}, P::Resample{T}, Sx::Int, Sy::Int, + B::Int) where {T} + n = P.n; buf = P.buf; H = resample_pad(T); m = length(buf) - 2H; hp = cld(n, 2); corr = P.corr + cx = Ptr{Complex{T}}(px); cy = Ptr{Complex{T}}(py) + GC.@preserve P begin + pb = grid_ptr(buf) + for i in 1:B + @inbounds for j in 1:hp + buf[H+j] = unsafe_load(cx, i + Sx * (j - 1)) * corr[j] + end + @inbounds for j in hp+1:n + buf[H+m-n+j] = unsafe_load(cx, i + Sx * (j - 1)) * corr[j] + end + fill!(view(buf, H+hp+1:H+m-n+hp), zero(Complex{T})) + run!(pb, pb, cptr(P.work1), cptr(P.work2), P.fwd, 1) + wrap_grid!(buf, m) + resample_interp!(cy, i, Sy, pb, P.ws, n, m, Val(resample_width(T)), 0, n) + end + end + nothing +end + +run!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::Resample{T}, S::Int) where {T} = + run_lanes!(py, px, pw1, pw2, P, S, S, S) + +# Fewest lanes s of a composite kernel that take the batched pass (RsBatch); below it every lane is +# one resampled transform (Resample) at the full vector width. The choice depends on the length only, +# so serial and threaded plans take the same path. EPYC 9354 (AVX-512), one thread, batched / lane by +# lane for q = 37, 1009, 10007: s = 2: 1.7–3.8 (F64), 2.1–7.3 (F32); 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. Lane +# counts that are not a multiple of the vector's complex count lose most (their remainder runs narrow). +const RS_BATCH_LANES = 32 + +# The length-q DFTs of the s lanes of a composite kernel: batched, or lane by lane for few lanes. +# C: the chain type of a lane's grid transform (Chain, or RChain for the run-time plans). +leftover_kernel(::Type{T}, q::Int, σ::Int, s::Int, ::Type{C} = Chain) where {T,C} = + s >= RS_BATCH_LANES ? RsBatch{T}(q, σ, s) : Resample{T}(q, σ, C) + +# Lane by lane (the batched pass, RsBatch, has its own method below). +function leftover!(py::Ptr{T}, pz::Ptr{T}, tw::Vector{Complex{T}}, rs::Resample{T}, n::Int, s::Int, B::Int, Sy::Int) where {T} + for i in 0:B-1, k1 in 0:s-1 + twiddled_leftover!(py, pz, tw, rs, n, s, B, i, k1, Sy) + end + nothing +end + +# The length-q DFT on lane (i, k1) of the chain output z (element j2 at i + B·(j2 + q·k1)), after the +# twiddles w_n^(j2·k1), folded into the deconvolution; output element k2 at y + i + Sy·(k1 + s·k2). +function twiddled_leftover!(py::Ptr{T}, pz::Ptr{T}, tw::Vector{Complex{T}}, rs::Resample{T}, n::Int, s::Int, + B::Int, i::Int, k1::Int, Sy::Int) where {T} + q = rs.n; buf = rs.buf; H = resample_pad(T); M = length(buf) - 2H; hp = cld(q, 2); corr = rs.corr + cz = Ptr{Complex{T}}(pz); cy = Ptr{Complex{T}}(py) + GC.@preserve rs tw begin + pb = grid_ptr(buf) + zb = i + B * q * k1 + 1 + e = 0 # (j2·k1) mod n + @inbounds for j in 1:hp + buf[H+j] = unsafe_load(cz, zb + B * (j - 1)) * tw[e+1] * corr[j] + e += k1; e >= n && (e -= n) + end + @inbounds for j in hp+1:q + buf[H+M-q+j] = unsafe_load(cz, zb + B * (j - 1)) * tw[e+1] * corr[j] + e += k1; e >= n && (e -= n) + end + fill!(view(buf, H+hp+1:H+M-q+hp), zero(Complex{T})) + run!(pb, pb, cptr(rs.work1), cptr(rs.work2), rs.fwd, 1) + wrap_grid!(buf, M) + resample_interp!(cy, i + Sy * k1 + 1, Sy * s, pb, rs.ws, q, M, Val(resample_width(T)), 0, q) + end + nothing +end + +# Items of the leftover pass of a threaded composite kernel (n = s·q, one slab lane) for ntasks tasks, +# and item b on the task's kernel: lanes one by one, or chunks of the batched pass (one per task where +# the lanes allow, multiples of the vector's complex count, see rs_chunk!). +leftover_chunk(R::Resample, s::Int, ntasks::Int) = 1 +leftover_chunk(R::RsBatch{T}, s::Int, ntasks::Int) where {T} = + (CU = ncomplex(first(vec_ladder(T))); min(R.C, CU * cld(cld(s, ntasks), CU))) +leftover_item!(py::Ptr{T}, pz::Ptr{T}, tw::Vector{Complex{T}}, R::Resample{T}, n::Int, s::Int, Cm::Int, b::Int) where {T} = + twiddled_leftover!(py, pz, tw, R, n, s, 1, 0, b, 1) +leftover_item!(py::Ptr{T}, pz::Ptr{T}, tw::Vector{Complex{T}}, R::RsBatch{T}, n::Int, s::Int, Cm::Int, b::Int) where {T} = + rs_chunk!(py, pz, tw, R, n, s, 1, 1, 0, b * Cm, min(Cm, s - b * Cm)) + +# Grid length of the leftover kernel of a composite kernel. +grid_size(R::RsBatch) = R.m +grid_size(R::Resample{T}) where {T} = length(R.buf) - 2 * resample_pad(T) + +# Tree reduction of U vectors of U complex values, each vector one output's partial sums, into one +# vector of the U outputs: per level, pairs of vectors halve their segments (FFTJ's interleaved +# reduction with complex values as units). Returns the statements and the result in output order. +function tree_reduce_units(vecs, U) + labels = [fill(j, U) for j in 0:U-1] + lanes(us) = Tuple(l for u in us for l in (2u, 2u + 1)) + stmts = Any[] + L = U + while L > 1 + half = L ÷ 2 + nvecs = Any[]; nlabels = Vector{Int}[] + for p in 1:2:length(vecs) + ip = Int[]; iq = Int[]; lab = Int[] + for u in 0:U-1 + seg, pos = divrem(u, L) + first = pos < half + sw = seg * L + (pos + half) % L + push!(ip, first ? u : U + u) + push!(iq, first ? sw : U + sw) + push!(lab, first ? labels[p][u+1] : labels[p+1][u+1]) + end + t = gensym(:t) + A, B = vecs[p], vecs[p+1] + push!(stmts, :($t = shufflevector($A, $B, Val($(lanes(ip)))) + shufflevector($A, $B, Val($(lanes(iq)))))) + push!(nvecs, t); push!(nlabels, lab) + end + vecs, labels = nvecs, nlabels + L = half + end + perm = [findfirst(==(u), labels[1]) - 1 for u in 0:U-1] + res = perm == 0:U-1 ? vecs[1] : :(shufflevector($(vecs[1]), Val($(lanes(perm))))) + stmts, res +end + +# Output k = Σ_t weight_t · grid[(l + t) mod m], read as grid[l + t] from a grid with filled margins +# (pb points at grid value 0, see wrap_grid!). Blocks of U = V/2 outputs (V lanes in the widest +# vector): per output one load of the w interleaved grid values, folded with the weights into V lanes +# by FMAs, then one tree reduction for the block and one store. The outputs left over at the end one +# at a time. Outputs k0 … k1-1 (a range of them per task on threads). +@generated function resample_interp!(cy::Ptr{Complex{T}}, i::Int, Sy::Int, pb::Ptr{T}, ws::WS, n::Int, + m::Int, ::Val{w}, k0::Int, k1::Int) where {T,WS,w} + V = 2 * ncomplex(first(vec_ladder(T))) + U = V ÷ 2 + ch(v, c) = :(shufflevector($v, Val($(Tuple(c*V:c*V+V-1))))) + fold(g, wv) = foldl((acc, c) -> :(muladd($(ch(g, c)), $(ch(wv, c)), $acc)), 1:2w÷V-1; + init = :($(ch(g, 0)) * $(ch(wv, 0)))) + blk = Any[] + fs = Any[] + for j in 0:U-1 + f = gensym(:f) + push!(blk, quote + l, wv = rweights(ws, k + $j, st, n, Val($w)) + g = vload(Vec{$(2w),T}, pb + l * $(2sizeof(T))) + $f = $(fold(:g, :wv)) + st = radvance(ws, st, mq, mr, n) + end) + push!(fs, f) + end + stmts, res = tree_reduce_units(fs, U) + blocked = U >= 2 && (2w) % V == 0 + quote + mq, mr = divrem(m, n) + st = rstart(ws, k0, m, n) + k = k0 + GC.@preserve ws @inbounds while k < k1 + if $blocked && k + $U <= k1 + $(blk...) + $(stmts...) + if Sy == 1 + vstore($res, Ptr{T}(cy) + (i - 1 + k) * $(2sizeof(T))) + else + v = $res + for u in 0:$(U - 1) + unsafe_store!(cy, Complex{T}(v[2u+1], v[2u+2]), i + Sy * (k + u)) + end + end + k += $U + continue + end + l, wv = rweights(ws, k, st, n, Val($w)) + p = vload(Vec{$(2w),T}, pb + l * $(2sizeof(T))) * wv + acc = Complex{T}(sum(shufflevector(p, Val($(Tuple(0:2:2w-1))))), sum(shufflevector(p, Val($(Tuple(1:2:2w-1)))))) + unsafe_store!(cy, acc, i + Sy * k) + st = radvance(ws, st, mq, mr, n) + k += 1 + end + nothing + end +end + +# The composite kernel's length-q DFTs by the batched pass: for each slab lane i, the s lanes k1 +# (z element j2 at i + B·(j2 + q·k1)) chunk by chunk into the grid with the twiddles w_n^(j2·k1) +# and the deconvolution, one batched transform, then rows of output (y + i + Sy·(k1 + s·k2)). +function leftover!(py::Ptr{T}, pz::Ptr{T}, tw::Vector{Complex{T}}, R::RsBatch{T}, n::Int, s::Int, + B::Int, Sy::Int) where {T} + for i in 0:B-1, c0 in 0:R.C:s-1 + rs_chunk!(py, pz, tw, R, n, s, B, Sy, i, c0, min(R.C, s - c0)) + end + nothing +end + +# Lanes c0 … c0+L-1 (L ≤ C, c0 a multiple of the vector's complex count CU) of slab lane i: one +# chunk of the batched pass, its multiple of CU lanes and its remainder as separate groups. The +# chain's vector width follows the lane count, and vectors of one complex value round differently; +# this way only the last s % CU lanes ever take the narrow path, so the threaded composite kernel +# (other chunk sizes, multiples of CU) agrees bitwise with the serial one. +function rs_chunk!(py::Ptr{T}, pz::Ptr{T}, tw::Vector{Complex{T}}, R::RsBatch{T}, n::Int, s::Int, B::Int, + Sy::Int, i::Int, c0::Int, L::Int) where {T} + r = L % ncomplex(first(vec_ladder(T))) + L > r && rs_group!(py, pz, tw, R, n, s, B, Sy, i, c0, L - r) + r > 0 && rs_group!(py, pz, tw, R, n, s, B, Sy, i, c0 + L - r, r) + nothing +end + +function rs_group!(py::Ptr{T}, pz::Ptr{T}, tw::Vector{Complex{T}}, R::RsBatch{T}, n::Int, s::Int, B::Int, + Sy::Int, i::Int, c0::Int, L::Int) where {T} + q = R.n; m = R.m; hp = cld(q, 2); corr = R.corr; g = R.g + cz = Ptr{Complex{T}}(pz); cy = Ptr{Complex{T}}(py) + GC.@preserve R tw begin + pg = cptr(g) + @inbounds for j2 in 0:q-1 + jj = j2 < hp ? j2 : m - q + j2 + e = (j2 * c0) % n # (j2·k1) mod n, stepped over the lanes + cr = corr[j2+1] + zb = i + B * (j2 + q * c0) + 1 + for c in 0:L-1 + g[jj*L+c+1] = unsafe_load(cz, zb + B * q * c) * (tw[e+1] * cr) + e += j2; e >= n && (e -= n) + end + end + fill!(view(g, hp*L+1:(m-q+hp)*L), zero(Complex{T})) + run!(pg, pg, cptr(R.work1), cptr(R.work2), R.fwd, L) + rs_rows!(cy, i + Sy * c0 + 1, Sy, Sy * s, pg, R.ws, q, m, L, Val(resample_width(T))) + end + nothing +end + +# Output rows of the batched pass: y[base + Ss·k2 + Sy·c] = Σ_t weight_t · grid row (l + t) mod m, +# lane c; whole vectors of lanes when Sy = 1, the rest one lane at a time. +@generated function rs_rows!(cy::Ptr{Complex{T}}, base::Int, Sy::Int, Ss::Int, pg::Ptr{T}, ws::WS, q::Int, + m::Int, L::Int, ::Val{w}) where {T,WS,w} + CU = ncomplex(first(vec_ladder(T))) + VT = Vec{2CU,T} + CT = 2sizeof(T) + o = [Symbol(:o, t) for t in 0:w-1] + b = [Symbol(:b, t) for t in 0:w-1] + offs = [:($(o[t+1]) = (g = l + $t; g += ifelse(g < 0, m, 0); g -= ifelse(g >= m, m, 0); g * L)) for t in 0:w-1] + bcast = [:($(b[t+1]) = $VT(wv[$(2t + 1)])) for t in 0:w-1] + vsum = foldl((acc, t) -> :(muladd($(b[t+1]), vload($VT, pg + ($(o[t+1]) + c) * $CT), $acc)), 1:w-1; + init = :($(b[1]) * vload($VT, pg + ($(o[1]) + c) * $CT))) + ssum = foldl((acc, t) -> :($acc + wv[$(2t + 1)] * unsafe_load(Ptr{Complex{T}}(pg), $(o[t+1]) + c + 1)), 1:w-1; + init = :(wv[1] * unsafe_load(Ptr{Complex{T}}(pg), $(o[1]) + c + 1))) + quote + mq, mr = divrem(m, q) + st = rstart(ws, 0, m, q) + GC.@preserve ws @inbounds for k2 in 0:q-1 + l, wv = rweights(ws, k2, st, q, Val($w)) + $(offs...) + yo = base + Ss * k2 + c = 0 + if Sy == 1 + $(bcast...) + while c + $CU <= L + vstore($vsum, Ptr{T}(cy) + (yo - 1 + c) * $CT) + c += $CU + end + end + while c < L + unsafe_store!(cy, $ssum, yo + Sy * c) + c += 1 + end + st = radvance(ws, st, mq, mr, q) + end + nothing + end +end diff --git a/src/rplan.jl b/src/rplan.jl index ecb6df92..12e3a740 100644 --- a/src/rplan.jl +++ b/src/rplan.jl @@ -1,11 +1,11 @@ -# Plans whose kernels are all run-time data: RChain for smooth lengths, Bluestein over RChain +# Plans whose kernels are all run-time data: RChain for smooth lengths, resampling over RChain # otherwise. Sizes, region and kernels live in vectors, so constructing a plan never builds a type # from run-time values and executing it compiles nothing per size; this is the path for code # compiled ahead of time (`juliac --trim`) and keeps the first plan of a new size cheap. # Same strategies as FFTLPlan: four-step for one long line on several threads, dimensions 1 and # 2 fused per plane when a plane fits in cache. -const RBluestein{T} = Bluestein{T,RChain{T},RChain{T}} +const RResample{T} = Union{Resample{T,RChain{T},RTable{T}},Resample{T,RChain{T},RPoly{T}}} const RFourStep{T} = FourStep{T,RChain{T},RChain{T}} """ @@ -20,15 +20,15 @@ mutable struct RPlan{T<:FFTLFloat} <: AbstractFFTs.Plan{Complex{T}} region::Vector{Int} σ::Int chains::Vector{RChain{T}} # per dimension in region (smooth lengths) - blue::Vector{Union{Nothing,RBluestein{T}}} # per dimension: Bluestein, or nothing + rs::Vector{Union{Nothing,RResample{T}}} # per dimension: resampling, or nothing four::Union{Nothing,RFourStep{T}} # one long contiguous line on several threads work1::Vector{Complex{T}} work2::Vector{Complex{T}} nthreads::Int tbuf::Vector{Vector{Complex{T}}} # two work buffers per thread pinv::AbstractFFTs.Plan{Complex{T}} - RPlan{T}(sz, region, σ, chains, blue, four, w1, w2, nt, tbuf) where {T} = - new{T}(sz, region, σ, chains, blue, four, w1, w2, nt, tbuf) + RPlan{T}(sz, region, σ, chains, rs, four, w1, w2, nt, tbuf) where {T} = + new{T}(sz, region, σ, chains, rs, four, w1, w2, nt, tbuf) end function RPlan{T}(sz::AbstractVector{<:Integer}, region::AbstractVector{<:Integer}, σ::Integer; @@ -38,7 +38,7 @@ function RPlan{T}(sz::AbstractVector{<:Integer}, region::AbstractVector{<:Intege reg = normregion(region, N) nt = check_plan(dims, σ, threads) chains = RChain{T}[] - blue = Union{Nothing,RBluestein{T}}[] + rs = Union{Nothing,RResample{T}}[] # one long line (the other transformed dimensions have length 1): the four-step algorithm long = filter(d -> dims[d] > 1, reg) four = length(long) == 1 && fourstep_dim(dims, long[1], nt) ? FourStep{T}(dims[long[1]], Int(σ), nt, RChain) : nothing @@ -47,14 +47,14 @@ function RPlan{T}(sz::AbstractVector{<:Integer}, region::AbstractVector{<:Intege S, _ = dim_extent(dims, d) if smooth(n) push!(chains, RChain{T}(n, Int(σ), radices(n, S; esize = sizeof(Complex{T})))) - push!(blue, nothing) + push!(rs, nothing) else push!(chains, RChain{T}(1, Int(σ), Int[])) - push!(blue, Bluestein{T}(n, Int(σ), RChain)) + push!(rs, Resample{T}(n, Int(σ), RChain)) end end nwork, ntwork = work_sizes(dims, reg, nt) - RPlan{T}(dims, reg, Int(σ), chains, blue, four, Vector{Complex{T}}(undef, nwork), Vector{Complex{T}}(undef, nwork), + RPlan{T}(dims, reg, Int(σ), chains, rs, four, Vector{Complex{T}}(undef, nwork), Vector{Complex{T}}(undef, nwork), nt, thread_buffers(T, ntwork, 2nt)) end @@ -78,7 +78,7 @@ function execute!(p::RPlan{T}, py::Ptr{Complex{T}}, px::Ptr{Complex{T}}) where { return nothing end i0 = 1 - if length(p.region) >= 2 && p.region[1] == 1 && p.region[2] == 2 && p.blue[1] === nothing && p.blue[2] === nothing + if length(p.region) >= 2 && p.region[1] == 1 && p.region[2] == 2 && p.rs[1] === nothing && p.rs[2] === nothing n1 = p.sz[1]; n2 = p.sz[2] outer = dim_extent(p.sz, 2)[2] if fuse_planes(n1, n2, outer, p.nthreads) @@ -97,7 +97,7 @@ function execute!(p::RPlan{T}, py::Ptr{Complex{T}}, px::Ptr{Complex{T}}) where { outer *= p.sz[j] end src = i == 1 ? rx : ry - b = p.blue[i] + b = p.rs[i] if b === nothing run_dim!(ry, src, p, p.chains[i], S, p.sz[d], outer) else diff --git a/src/rreal.jl b/src/rreal.jl index 87d5bddb..36c257f4 100644 --- a/src/rreal.jl +++ b/src/rreal.jl @@ -1,5 +1,5 @@ # Run-time plans of the real transforms (r2c, c2r, r2r) for arrays of any number of dimensions: -# sizes, region and kernels are data (RChain, Bluestein over RChain), only the element type is +# sizes, region and kernels are data (RChain, resampling over RChain), only the element type is # a type parameter. They run the line kernels of real.jl and r2r.jl and are the real-data # counterpart of RPlan for code compiled ahead of time (libfftl). diff --git a/src/threads.jl b/src/threads.jl index d51772de..78186e00 100644 --- a/src/threads.jl +++ b/src/threads.jl @@ -244,7 +244,7 @@ end const ROW_PAD = 4 -# One line on the calling thread (the convolution transforms of a large Bluestein kernel). +# One line on the calling thread (the grid transform of a large resampling kernel). function run!(py::Ptr{T}, px::Ptr{T}, ::Ptr{T}, ::Ptr{T}, P::FourStep{T}, S::Int) where {T} S == 1 || throw(ArgumentError("a four-step kernel transforms one line")) run_fourstep!(py, px, P, 1) @@ -287,46 +287,57 @@ function run_fourstep!(py::Ptr{T}, px::Ptr{T}, P::FourStep{T}, nt::Int) where {T end """ - BluesteinMT{T} + ResampleMT{T} -Bluestein's algorithm for one long contiguous line of a length with a prime factor above -`MAX_RADIX_PRIME`, on several threads: the two convolution transforms of length M run as -four-step transforms and the three element-wise passes (chirp product with zero padding, -filter product, chirp product) are split over the threads. +A resampled DFT (see `Resample`) for one long contiguous line of a length with a prime factor +above `MAX_RADIX_PRIME`, on several threads: the transform of length m runs as a four-step +transform, the deconvolution and the interpolation (ranges of outputs) are split over the threads. """ -struct BluesteinMT{T,F1,F2} <: Kernel{T} +struct ResampleMT{T,F,WS} <: Kernel{T} n::Int - chirp::Vector{Complex{T}} # c_j, j = 0…n-1 - kernel::Vector{Complex{T}} # DFT of the convolution kernel, scaled by 1/M - fwd::F1 - bwd::F2 - buf::Vector{Complex{T}} # length M + corr::Vector{T} + ws::WS + fwd::F + buf::Vector{Complex{T}} # the grid between two margins (resample_grid) end -# Smallest 5-smooth M ≥ 2n - 1 that the four-step algorithm can split. -function bluestein_mt_length(::Type{T}, n::Int) where {T} - m = bluestein_length(T, 2n - 1) +function ResampleMT{T}(n::Int, σ::Int, nt::Int) where {T} + m = grid_length(T, 2n) while fourstep_split(m) == 0 - m = bluestein_length(T, m + 1) + m = grid_length(T, m + 1) end - m + _, _, corr, ws = resample_parts(T, n, m) + fwd = FourStep{T}(m, σ, nt) + ResampleMT{T,typeof(fwd),typeof(ws)}(n, corr, ws, fwd, resample_grid(T, m)) end -function BluesteinMT{T}(n::Int, σ::Int, nt::Int) where {T} - M = bluestein_mt_length(T, n) - ser = RChain{T}(M, -1) - w1 = Vector{Complex{T}}(undef, M); w2 = Vector{Complex{T}}(undef, M) - chirp, b = cached(() -> bluestein_tables(T, n, σ, ser, w1, w2), CacheKey(:bluestein, T, n, σ, Int[]) - )::Tuple{Vector{Complex{T}},Vector{Complex{T}}} - fwd = FourStep{T}(M, -1, nt); bwd = FourStep{T}(M, 1, nt) - BluesteinMT{T,typeof(fwd),typeof(bwd)}(n, chirp, b, fwd, bwd, Vector{Complex{T}}(undef, M)) +function run_resample_mt!(py::Ptr{T}, px::Ptr{T}, P::ResampleMT{T}, nt::Int) where {T} + n = P.n; m = length(P.buf) - 2 * resample_pad(T); hp = cld(n, 2) + GC.@preserve P begin + cx = Ptr{Complex{T}}(px); cy = Ptr{Complex{T}}(py) + pg = grid_ptr(P.buf) + pb = Ptr{Complex{T}}(pg); pc = pointer(P.corr) + chunked(m, nt) do j + jj = j < hp ? j : j - (m - n) # input index of grid point j (or a padding zero) + v = j < hp || jj >= hp ? unsafe_load(cx, jj + 1) * unsafe_load(pc, jj + 1) : zero(Complex{T}) + unsafe_store!(pb, v, j + 1) + end + run_fourstep!(pg, pg, P.fwd, nt) + wrap_grid!(P.buf, m) + nch = min(BLOCKS_PER_THREAD * nt, cld(n, 1 << 12)) + parallel_items(nch, nt) do _, it + r = chunk_range(it + 1, nch, n) + isempty(r) || resample_interp!(cy, 1, 1, pg, P.ws, n, m, Val(resample_width(T)), first(r), last(r) + 1) + end + end + nothing end """ CompositeMT{T} -The composite kernel (first stage of length s, Bluestein of length q) for one long contiguous -line on several threads: the first stage in lane blocks over the tasks, then the s Bluestein +The composite kernel (first stage of length s, resampled DFT of length q) for one long contiguous +line on several threads: the first stage in lane blocks over the tasks, then the s resampled transforms spread over the tasks. Kernels with buffers of their own exist once per task. The blocks and the arithmetic are those of the serial kernel, so the results agree bitwise. """ @@ -335,7 +346,7 @@ struct CompositeMT{T,C,B} <: Kernel{T} s::Int q::Int chains::Vector{C} # first stage, one per task - blues::Vector{B} # one per task + rss::Vector{B} # leftover kernel (RsBatch, or Resample for few lanes), one per task tw::Vector{Complex{T}} # w_n^m, m = 0…n-1 z::Vector{Complex{T}} # output of the first stage tbuf::Vector{Vector{Complex{T}}} # two per task, one lane block of the first stage each @@ -345,16 +356,16 @@ function CompositeMT{T}(n::Int, σ::Int, nt::Int) where {T} s, q = composite_split(T, n) c1 = first_kernel(T, s, σ, q, 1) chains = [c1; [c1 isa AnyChain ? c1 : first_kernel(T, s, σ, q, 1) for _ in 2:nt]] - b1 = Bluestein{T}(q, σ) - blues = [b1; [Bluestein{T}(q, σ) for _ in 2:nt]] + r1 = leftover_kernel(T, q, σ, s) + rss = [r1; [leftover_kernel(T, q, σ, s) for _ in 2:nt]] tw = cached(() -> Complex{T}.(unit_roots(σ, n)), CacheKey(:composite, T, n, σ, Int[]))::Vector{Complex{T}} - CompositeMT{T,typeof(c1),typeof(b1)}(n, s, q, chains, blues, tw, Vector{Complex{T}}(undef, n), + CompositeMT{T,typeof(c1),typeof(r1)}(n, s, q, chains, rss, tw, Vector{Complex{T}}(undef, n), thread_buffers(T, composite_block(q, s) * s, 2nt)) end -# A long line n = s·q (see composite_split) with at least one Bluestein transform per thread. +# A long line n = s·q (see composite_split) with at least one leftover lane per thread. function use_composite_mt(::Type{T}, n, nt) where {T} - (nt > 1 && n >= TUNING[].bluestein_mt_min) || return false + (nt > 1 && n >= TUNING[].resample_mt_min) || return false s, q = composite_split(T, n) q > 1 && s >= nt end @@ -366,7 +377,8 @@ function run_composite_mt!(py::Ptr{T}, px::Ptr{T}, P::CompositeMT{T}, nt::Int) w CT = 2 * sizeof(T) Bb = composite_block(q, s) nb = cld(q, Bb) - ntasks = min(nt, length(P.blues)) + ntasks = min(nt, length(P.rss)) + Cm = leftover_chunk(P.rss[1], s, ntasks) GC.@preserve P begin pz = cptr(P.z) parallel_phases(ntasks, @@ -374,7 +386,7 @@ function run_composite_mt!(py::Ptr{T}, px::Ptr{T}, P::CompositeMT{T}, nt::Int) w local a = b * Bb run_lanes!(pz + a * CT, px + a * CT, cptr(P.tbuf[2c-1]), cptr(P.tbuf[2c]), P.chains[c], q, q, min(Bb, q - a)) end, - s, (c, k1) -> twiddled_bluestein!(py, pz, P.tw, P.blues[c], n, s, 1, 0, k1, 1)) + cld(s, Cm), (c, b) -> leftover_item!(py, pz, P.tw, P.rss[c], n, s, Cm, b)) end nothing end @@ -383,7 +395,7 @@ end LinesMT{T} Many lines of a length with a prime factor above `MAX_RADIX_PRIME` on several threads: the kernels -(Bluestein or composite) hold buffers of their own, so each task has one, and the tasks take blocks +(resampling or composite) hold buffers of their own, so each task has one, and the tasks take blocks of lanes of every slice from a shared counter. """ struct LinesMT{T,K} <: Kernel{T} @@ -403,11 +415,11 @@ end lines_mt_dim(sz, d::Int, nt::Int) = !smooth(sz[d]) && (so = dim_extent(sz, d); so[1] * so[2] > 1 && threaded(so[1], sz[d], so[2], nt)) -# From this length on, a Bluestein line runs on several threads. +# From this length on, a resampled line runs on several threads. -use_bluestein_mt(n, nt) = nt > 1 && n >= TUNING[].bluestein_mt_min && !smooth(n) +use_resample_mt(n, nt) = nt > 1 && n >= TUNING[].resample_mt_min && !smooth(n) -bluestein_mt_dim(sz, d::Int, nt::Int) = (so = dim_extent(sz, d); so[1] == 1 && so[2] == 1 && use_bluestein_mt(sz[d], nt)) +resample_mt_dim(sz, d::Int, nt::Int) = (so = dim_extent(sz, d); so[1] == 1 && so[2] == 1 && use_resample_mt(sz[d], nt)) # f(i) for i in chunk `it` of `nchunks` chunks of 0…N-1, the chunks over nt tasks. @inline function chunked(f::F, N, nt) where {F} @@ -419,27 +431,6 @@ bluestein_mt_dim(sz, d::Int, nt::Int) = (so = dim_extent(sz, d); so[1] == 1 && s end end -function run_bluestein_mt!(py::Ptr{T}, px::Ptr{T}, P::BluesteinMT{T}, nt::Int) where {T} - n = P.n; M = length(P.buf) - GC.@preserve P begin - cx = Ptr{Complex{T}}(px); cy = Ptr{Complex{T}}(py) - pb = Ptr{Complex{T}}(cptr(P.buf)) - pc = Ptr{Complex{T}}(cptr(P.chirp)); pk = Ptr{Complex{T}}(cptr(P.kernel)) - chunked(M, nt) do j - unsafe_store!(pb, j < n ? unsafe_load(cx, j + 1) * unsafe_load(pc, j + 1) : zero(Complex{T}), j + 1) - end - run_fourstep!(cptr(P.buf), cptr(P.buf), P.fwd, nt) - chunked(M, nt) do j - unsafe_store!(pb, unsafe_load(pb, j + 1) * unsafe_load(pk, j + 1), j + 1) - end - run_fourstep!(cptr(P.buf), cptr(P.buf), P.bwd, nt) - chunked(n, nt) do k - unsafe_store!(cy, unsafe_load(pb, k + 1) * unsafe_load(pc, k + 1), k + 1) - end - end - nothing -end - # A[k2 + (n2 + ROW_PAD)·(i0 + i)] = w[i + Bc·k2] · tw[k2 + n2·(i0 + i)] for i < Bc, k2 < n2: # contiguous writes along k2, strided reads from the cache buffer, one complex value per vector. function twiddle_transpose!(pa::Ptr{Complex{T}}, pw::Ptr{Complex{T}}, ptw::Ptr{Complex{T}}, i0, Bc, n2) where {T} @@ -513,8 +504,8 @@ end run_kernel_dim!(py::Ptr{T}, px::Ptr{T}, p, k::FourStep{T}, S::Int, n::Int, outer::Int) where {T} = run_fourstep!(py, px, k, p.nthreads) -run_kernel_dim!(py::Ptr{T}, px::Ptr{T}, p, k::BluesteinMT{T}, S::Int, n::Int, outer::Int) where {T} = - run_bluestein_mt!(py, px, k, p.nthreads) +run_kernel_dim!(py::Ptr{T}, px::Ptr{T}, p, k::ResampleMT{T}, S::Int, n::Int, outer::Int) where {T} = + run_resample_mt!(py, px, k, p.nthreads) run_kernel_dim!(py::Ptr{T}, px::Ptr{T}, p, k::CompositeMT{T}, S::Int, n::Int, outer::Int) where {T} = run_composite_mt!(py, px, k, p.nthreads) @@ -562,7 +553,7 @@ function run_kernel_dim!(py::Ptr{T}, px::Ptr{T}, p, k::Union{Chain{T},RChain{T}} nothing end -# Other kernels (Bluestein) and chains that are not threaded: one slice after the other. +# Other kernels (resampling, composite) and chains that are not threaded: one slice after the other. run_kernel_dim!(py::Ptr{T}, px::Ptr{T}, p, k::Kernel{T}, S::Int, n::Int, outer::Int) where {T} = run_serial_dim!(py, px, p, k, S, n, outer) diff --git a/src/tuning.jl b/src/tuning.jl index 1eb1d992..5171f9da 100644 --- a/src/tuning.jl +++ b/src/tuning.jl @@ -12,9 +12,10 @@ Size thresholds of FFTL (see `TUNING`): `clamp(prevpow2(L2 ÷ 32), 2^12, block_budget)`) - `fourstep_min`: length from which a single line uses the four-step algorithm on several threads - `factor_twiddles`: bytes above which the four-step uses factor twiddle tables -- `bluestein_mt_min`: length from which Bluestein's algorithm runs on several threads +- `resample_mt_min`: length from which a length with a large prime factor (resampling) runs on + several threads - `real_fourstep_bytes`: bytes of a real line from which it runs as a four-step on several threads -- `real_bluestein_bytes`: bytes of a real line from which a non-smooth one runs on several threads +- `real_resample_bytes`: bytes of a real line from which a non-smooth one runs on several threads - `plane_budget`: complex values of a plane (dimensions 1 and 2) that may be fused """ struct Tuning @@ -22,9 +23,9 @@ struct Tuning block_budget::Int fourstep_min::Int factor_twiddles::Int - bluestein_mt_min::Int + resample_mt_min::Int real_fourstep_bytes::Int - real_bluestein_bytes::Int + real_resample_bytes::Int plane_budget::Int end diff --git a/test/fftw_compat.jl b/test/fftw_compat.jl index df050ea1..513d7efb 100644 --- a/test/fftw_compat.jl +++ b/test/fftw_compat.jl @@ -1,26 +1,7 @@ # The interface under FFTW.jl's names, the scope rules and AbstractFFTs's generic functions. -# Needs FFTW.jl (test dependency); the first testset runs before FFTW is loaded. +# Needs FFTW.jl (test dependency); the test without FFTW runs first, in runtests.jl. import AbstractFFTs -fftw_loaded() = any(m -> nameof(m) === :FFTW, values(Base.loaded_modules)) - -@testset "AbstractFFTs without FFTW" begin - if fftw_loaded() - @test_skip "FFTW already loaded" - else - x = randn(ComplexF64, 12, 10) - @test AbstractFFTs.plan_fft(x) isa FFTL.FFTLPlan - @test AbstractFFTs.fft(x) ≈ FFTL.fft(x) - @test AbstractFFTs.fft(real(x)) ≈ FFTL.fft(real(x)) - @test AbstractFFTs.ifft(AbstractFFTs.fft(x)) ≈ x - r = randn(Float32, 9, 4) - @test AbstractFFTs.irfft(AbstractFFTs.rfft(r), 9) ≈ r - @test AbstractFFTs.plan_rfft(r) isa FFTL.RFFTLPlan - y = copy(x); AbstractFFTs.fft!(y) - @test y ≈ FFTL.fft(x) - end -end - @eval using FFTW @testset "Same results as FFTW.jl" begin diff --git a/test/inference.jl b/test/inference.jl index 14f97e64..ff39da8a 100644 --- a/test/inference.jl +++ b/test/inference.jl @@ -18,7 +18,7 @@ @test_broken @inferred FFTL.plan_fft(x) @test_broken @inferred FFTL.plan_fft(x, (1, 2)) end - # Bluestein length, the plan holds the work buffers + # a length with a large prime factor (resampling), the plan holds the work buffers let x = randn(Complex{T}, 37), y = similar(x) p = FFTL.plan_fft(x) mul!(y, p, x) @@ -159,7 +159,7 @@ exe_r2r!(p, y, x) = GC.@preserve x y FFTL.execute_r2r!(p, pointer(y), pointer(x) @test (@inferred exe!(p, y, x)) === nothing @test (@allocated exe!(p, y, x)) == 0 end - let p = (@inferred FFTL.RPlan{T}([37], [1], 1; threads = 1)), x = randn(Complex{T}, 37), y = similar(x) # Bluestein + let p = (@inferred FFTL.RPlan{T}([37], [1], 1; threads = 1)), x = randn(Complex{T}, 37), y = similar(x) # resampling exe!(p, y, x) @test (@allocated exe!(p, y, x)) == 0 end diff --git a/test/runtests.jl b/test/runtests.jl index cd5f6849..33378c34 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,13 +1,32 @@ using FFTL, Test, LinearAlgebra - -# Reference DFT in Float64 (BigFloat would be exact but slow; the tolerances below allow for it). -function dft(x, σ = -1) - n = length(x) - [sum(x[j+1] * cispi(2σ * mod(j * k, n) / n) for j in 0:n-1) for k in 0:n-1] +import AbstractFFTs + +# The AbstractFFTs generics with FFTL as the only backend: before FFTW.jl is loaded (below). +fftw_loaded() = any(m -> nameof(m) === :FFTW, values(Base.loaded_modules)) + +@testset "AbstractFFTs without FFTW" begin + if fftw_loaded() + @test_skip "FFTW already loaded" + else + x = randn(ComplexF64, 12, 10) + @test AbstractFFTs.plan_fft(x) isa FFTL.FFTLPlan + @test AbstractFFTs.fft(x) ≈ FFTL.fft(x) + @test AbstractFFTs.fft(real(x)) ≈ FFTL.fft(real(x)) + @test AbstractFFTs.ifft(AbstractFFTs.fft(x)) ≈ x + r = randn(Float32, 9, 4) + @test AbstractFFTs.irfft(AbstractFFTs.rfft(r), 9) ≈ r + @test AbstractFFTs.plan_rfft(r) isa FFTL.RFFTLPlan + y = copy(x); AbstractFFTs.fft!(y) + @test y ≈ FFTL.fft(x) + end end -# The DFT at a few frequencies only, each in O(n): the full reference is O(n²), too slow for lengths -# of 10⁵ and more. The frequencies cover both ends, the middle and a spread in between. +# Reference DFTs in Float64 from FFTW.jl: an independent implementation, O(n log n) at every size. +import FFTW +dft(x, σ = -1) = σ == -1 ? FFTW.fft(ComplexF64.(x)) : FFTW.bfft(ComplexF64.(x)) + +# The DFT at a few frequencies only, each in O(n), by direct summation. The frequencies cover both +# ends, the middle and a spread in between. function dft_at(x, σ = -1) n = length(x) ks = unique([0, 1, n ÷ 3, n ÷ 2, n - 2, n - 1, (7n) ÷ 11]) @@ -15,11 +34,11 @@ function dft_at(x, σ = -1) ks, ref end -# Direct DFT along the dimensions in `region`, one after the other. +# Along the dimensions in `region`, one after the other. function dftn(x, region, σ = -1) y = ComplexF64.(x) for d in region - y = mapslices(v -> dft(v, σ), y; dims = d) + y = σ == -1 ? FFTW.fft(y, d) : FFTW.bfft(y, d) end y end @@ -78,7 +97,7 @@ thin(T, v) = FULL || T == Float64 ? v : v[1:2:end] end end - @testset "Bluestein, $T" for T in (Float64, Float32) + @testset "prime factors above 31 (resampling), $T" for T in (Float64, Float32) for n in pick([37, 41, 74, 101, 37^2, 1009], [37, 74, 1009]) x = randn(Complex{T}, n) @test relerr(FFTL.fft(x), dft(ComplexF64.(x))) < 10tol(T) @@ -87,7 +106,7 @@ thin(T, v) = FULL || T == Float64 ? v : v[1:2:end] end @testset "smooth part times leftover, $T" for T in (Float64, Float32) - # n = s·q: a chain of length s, Bluestein on the leftover q (prime or composite) + # n = s·q: a chain of length s, a resampled DFT of the leftover q (prime or composite) for n in pick([2 * 37, 3 * 41, 16 * 43, 29 * 37, 64 * 37, 2 * 37^2, 5 * 13 * 37, 2 * 37 * 41, 6 * 1009, 2 * 37 * 41 * 43, 6 * 37^2 * 41], [2 * 37, 29 * 37, 2 * 37 * 41, 6 * 1009]) @@ -97,7 +116,7 @@ thin(T, v) = FULL || T == Float64 ? v : v[1:2:end] @test relerr(FFTL.bfft(x), dft(ComplexF64.(x), 1)) < 10tol(T) end # a leftover split into prime powers (the planner does this for large leftovers only): first - # stages that are composite and Bluestein kernels themselves + # stages that are composite and resampling kernels themselves for (n, s, q) in pick([(2 * 37 * 41 * 43, 2 * 37 * 41, 43), (37 * 41, 37, 41), (6 * 37^2 * 41, 6 * 37^2, 41)], [(2 * 37 * 41 * 43, 2 * 37 * 41, 43), (37 * 41, 37, 41)]) x = randn(Complex{T}, n) @@ -217,7 +236,7 @@ thin(T, v) = FULL || T == Float64 ? v : v[1:2:end] @test all(==(7), parent(w)[mask]) @test relerr(mul!(v, p, v), ref) < 10tol(T) end - # LinesMT owns separate Bluestein/composite kernels per task. Exercise first-axis + # LinesMT owns separate resampling/composite kernels per task. Exercise first-axis # slices, interior-axis lane blocks (including a partial block), and multiple passes. # The three-lane case has gather workspace with two planned tasks, but none with three. for (sz, region) in (((1009, 64), (1,)), ((65, 1009), (2,)), ((2018, 33), (1,)), @@ -268,24 +287,24 @@ thin(T, v) = FULL || T == Float64 ? v : v[1:2:end] end end - @testset "Bluestein on several threads, $T" for T in (Float64, Float32) - # one long line with a prime factor above MAX_RADIX_PRIME: threaded convolution transforms + @testset "resampling on several threads, $T" for T in (Float64, Float32) + # one long line with a prime factor above MAX_RADIX_PRIME: threaded grid transform for n in pick((16411, 20011, 100003), (20011,)) x = randn(Complex{T}, n) for (mk, σ) in ((FFTL.plan_fft, -1), (FFTL.plan_bfft, 1)) ref = mk(x; threads = 1) * x for nt in pick((2, 3, 4), (3,)) p = mk(x; threads = nt) - @test p.kernels[1] isa FFTL.BluesteinMT + @test p.kernels[1] isa FFTL.ResampleMT @test relerr(p * x, ref) < tol(T) z = copy(x); mul!(z, p, z) @test relerr(z, ref) < tol(T) end end end - # below the threshold the serial Bluestein kernel stays + # below the threshold the serial resampling kernel stays x = randn(Complex{T}, 4099) - @test FFTL.plan_fft(x; threads = 3).kernels[1] isa FFTL.Bluestein + @test FFTL.plan_fft(x; threads = 3).kernels[1] isa FFTL.Resample # smooth part times leftover: the composite kernel on threads, bitwise equal to the serial one # 3·193·1021: the leftover is split, the first stage is a composite kernel on each task for n in pick((65 * 4033, 7 * 11777, 1500 * 1013, 3 * 193 * 1021), (3 * 193 * 1021,)) @@ -302,7 +321,7 @@ thin(T, v) = FULL || T == Float64 ? v : v[1:2:end] end end # many lines with a prime factor above MAX_RADIX_PRIME: slices and lane blocks over the tasks, - # Bluestein kernels (127, 1009, 37) and composite kernels (129 = 3·43, 2018 = 2·1009) + # resampling kernels (127, 1009, 37) and composite kernels (129 = 3·43, 2018 = 2·1009) for (sz, region) in pick([((256, 127, 16), 1:3), ((16, 256, 127), (3,)), ((1009, 64), (1,)), ((64, 1009), (2,)), ((37, 37, 37), 1:3), ((129, 300), 1:2), ((300, 129), 1:2), ((2018, 32), (1,)), ((32, 2018), (2,))], [((256, 127, 16), 1:3), ((64, 1009), (2,)), ((1009, 64), (1,)), ((300, 129), 1:2)]) @@ -526,7 +545,7 @@ thin(T, v) = FULL || T == Float64 ? v : v[1:2:end] mask = trues(size(G)); mask[b...] .= false @test G[mask] == G0[mask] end - # out of place, Bluestein length along d, 2D + # out of place, a length with a large prime factor along d, 2D H0 = randn(Complex{T}, 45, 41); H1 = similar(H0) for d in 1:2 v = view(H0, 2:44, 3:39); w = view(H1, 2:44, 3:39) @@ -554,7 +573,7 @@ thin(T, v) = FULL || T == Float64 ? v : v[1:2:end] p = FFTL.plan_pass_runtime(randn(Complex{T}, sz), d) @test p isa FFTL.PassPlan{T,2,FFTL.RPassKern{T}} end - # against the length-specialised plan, every kernel kind: smooth, Bluestein, composite, four-step + # against the length-specialised plan, every kernel kind: smooth, resampling, composite, four-step for (n, sz_other, nt) in ((64, 5, 1), (100, 7, 1), (1000, 3, 2), (1009, 4, 1), (2 * 37 * 41, 2, 1), (2^17, 1, 3)) x = randn(Complex{T}, n, sz_other) @@ -582,7 +601,7 @@ thin(T, v) = FULL || T == Float64 ? v : v[1:2:end] # real passes: run-time plans have one type for all lengths; the results agree with the specialised ones for (sz, d) in (((100, 3), 1), ((1009, 3), 1), ((30, 2), 1)) x = randn(T, sz) - @test FFTL.plan_r2c_pass_runtime(x, d) isa FFTL.RealPassPlan{T,2,FFTL.RealLine{T,Union{FFTL.RChain{T},FFTL.Bluestein{T,FFTL.RChain{T},FFTL.RChain{T}}}}} + @test FFTL.plan_r2c_pass_runtime(x, d) isa FFTL.RealPassPlan{T,2,FFTL.RealLine{T,FFTL.RKern{T}}} ref = FFTL.plan_r2c_pass(x, d); rt = FFTL.plan_r2c_pass_runtime(x, d) ya = zeros(Complex{T}, ref.sz_out); yb = similar(ya) mul!(ya, ref, x); mul!(yb, rt, x) @@ -871,7 +890,7 @@ thin(T, v) = FULL || T == Float64 ? v : v[1:2:end] for (T, n) in ((Float64, 37888), (Float64, 92935), (Float64, 2 * 1009 * 37), (Float32, 33490), (Float32, 92935), (Float32, 100003)) x = randn(T, n) p = FFTL.plan_rfft(x; threads = 3) - @test p.four isa Union{FFTL.CompositeMT,FFTL.BluesteinMT} + @test p.four isa Union{FFTL.CompositeMT,FFTL.ResampleMT} X = p * x @test relerr(X, FFTL.plan_rfft(x; threads = 1) * x) < 10tol(T) @test relerr(FFTL.plan_brfft(X, n; threads = 3) * X, n .* x) < 10tol(T) @@ -882,7 +901,7 @@ thin(T, v) = FULL || T == Float64 ? v : v[1:2:end] end @testset "argument checks and fixed bugs" begin - # Bluestein lines in threaded real and r2r plans (the kernel's own buffers, one task at a time) + # resampled lines in threaded real and r2r plans (the kernel's own buffers, one task at a time) x = rand(1031, 40) @test FFTL.plan_rfft(x; threads = 4) * x == FFTL.plan_rfft(x; threads = 1) * x @test FFTL.plan_r2r(x, FFTL.REDFT10, 1; threads = 4) * x == FFTL.plan_r2r(x, FFTL.REDFT10, 1; threads = 1) * x