From a204ed8f86dcc53a8eb6298c7f0a6aed7769edf0 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 13 Nov 2026 02:39:51 +0000 Subject: [PATCH] One chain form: every smooth length runs on RChain The radix sequence and the stage records of a Stockham chain are data, and the codelets come from the fixed set. The chain type with the radix sequence in its type (Chain, used for smooth lengths up to 128) is removed, with its generated run!/run_lanes!, the 583 compiled chain methods and the dispatch trees from an RChain of a short length into them. Nothing is compiled per length any more. The run-time chain used to cost 9-20 ns more per call. That was not the choice of the codelet at run time (about 1 ns per stage) but the code around it, which the stage loop now avoids: - a line of one stage is one codelet call; - the stage records are one flat vector that carries the twiddle addresses and is read through a raw pointer, so the loop needs no GC frame; - RChain is a mutable struct with const fields, passed as one pointer instead of being copied to the stack of each caller; - the codelet switch is inlined into the loop; - the rows of a plane run inside one call. The addresses in the stage records belong to one chain object. A chain remembers its own address when it is linked and links again when it finds itself elsewhere (a copy, or a plan loaded from a package image, which keeps Ptr values as they were). The tag of the default plans in Bluestein, RealLine and R2RLine is the abstract type TypedKernels. AMD EPYC 9354 (AVX-512), one thread, plan_fft and mul!, 92 cases, 8 rounds in alternating processes, median time over that of the previous code: geometric mean 0.973, 27 cases below 0.95, 8 above 1.05 (short lines of two or three stages, at most 1.13). 1796 methods instead of 2385, package image 61 MB instead of 72 MB, libfftl.so 49.5 MB instead of 55.6 MB. Not verified on aarch64. Co-Authored-By: Claude Fable 5.2 --- docs/AGENTS.md | 6 +- docs/src/explanation/design-decisions.md | 65 ++-- docs/src/explanation/image-size.md | 9 +- docs/src/explanation/implementation.md | 47 +-- src/interface.jl | 4 +- src/kernel_methods.jl | 14 +- src/kernels.jl | 34 +- src/pass.jl | 6 +- src/plan.jl | 234 ++----------- src/precompile.jl | 4 +- src/r2r.jl | 2 +- src/rchain.jl | 408 +++++++++++------------ src/real.jl | 8 +- src/threads.jl | 6 +- test/runtests.jl | 79 +++-- 15 files changed, 371 insertions(+), 555 deletions(-) diff --git a/docs/AGENTS.md b/docs/AGENTS.md index 8f2fc521..5a047b98 100644 --- a/docs/AGENTS.md +++ b/docs/AGENTS.md @@ -35,8 +35,8 @@ rules there apply to everything written here as well. Everything in this folder - **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 - 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 + reductions), implementation (type specialisation, SIMD, kernel cache, cache sizes, the stage loop of the + chain, why libfftl needs `RChain`), design decisions (what was decided and why, with the evidence), package size and load time, performance (measurement protocol). - **Reference** (`src/reference/`): Julia API (`@docs`), AbstractFFTs (what is implemented), libfftl C API (prose on FFTW3's interface as libfftl implements it; the header pages are generated from @@ -46,7 +46,7 @@ rules there apply to everything written here as well. Everything in this folder ## Decisions for the manual -The design decisions of the package (interface, chain forms, rules for lengths, package size) are on the +The design decisions of the package (interface, the chain, rules for lengths, package size) are on the 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 diff --git a/docs/src/explanation/design-decisions.md b/docs/src/explanation/design-decisions.md index 437bef8a..0babbe08 100644 --- a/docs/src/explanation/design-decisions.md +++ b/docs/src/explanation/design-decisions.md @@ -74,29 +74,49 @@ with one thread, and from ×0.77 to ×1.16 with 16 threads. measured on Zen 2 only, and the 4 MiB rule is checked on 12 lengths that also took part in its derivation. -## Two chain forms: compiled for short smooth lengths, run-time data for the rest - -**Decision (2026-11-09).** `kernel` builds the compiled `Chain` (radix sequence in the type) for smooth -lengths n ≤ 128 (`SMALL_N`) and the `RChain` (radix sequence as data, a fixed set of codelets) for -every other smooth length. The choice is made when the plan is built, so a call carries no dispatch. -Both forms share the twiddle tables, the radix rule and the vector-type rules (bitwise equal in 552 -checks). `RChain` is also the only form the trimmed C library can contain. - -**Why.** A plan of a new size with `Chain` is compilation, not planning: 0.31–0.48 s for plan and first -call against 0.013–0.044 s with the run-time plan (M1, one fresh process per case). But on x86 the -run-time choice of the codelet costs 9–20 ns per call (n = 16, `Float64`: `Chain` 22 ns, `RChain` -39 ns), which is 1.2–1.5× for short lines. With the split, short lengths are at parity or better. Warm execution of the run-time plan against the -compiled one was 1.00–1.08 on the cases measured, with one exception (Bluestein, 100003, `Float32`: 1.42, -cause open). - -**Replaced.** The idea of removing `Chain` altogether (2026-11-08) and the earlier decision to keep the -JIT and run-time plan families as separate types (2026-11-02). +## One chain form: the radix sequence is data + +**Decision (2026-11-13).** Every smooth length runs on `RChain`: the radix sequence and the stage +records are data, the codelets come from a fixed set. There is no chain type with the radix sequence +in its type, and nothing is compiled per length. `RChain` is also the only form the trimmed C library +can contain. + +**Why.** The compiled chain (`Chain{T,Rs,σ,K}`, used for smooth lengths up to 128 until this decision) +existed because the run-time chain cost 9–20 ns more per call on x86. That cost was not the choice of +the codelet at run time (about 1 ns per stage) but the code around it: GC frames for references loaded +from the chain, the copy of an immutable chain into each caller, bounds checks, and one more call level +per stage. With those removed (see [Implementation](@ref "Implementation")) the compiled chain has no +advantage left. AMD EPYC 9354 (AVX-512), one thread, `plan_fft` and `mul!`, 92 cases (`Float64` and +`Float32`; lines of 2 to 4096, planes of 8 × 8 to 128 × 128, cubes of 8³ to 64³, in-place, backward and +real transforms), 8 rounds in alternating processes, median time over the median time with the compiled +chain: + +- geometric mean 0.973; 27 cases below 0.95, 8 above 1.05; +- lines of one stage (2, 4, 16): 0.80–0.95; 8 × 8, 8 × 8 × 8, 7 × 12, 5 × 64: 0.85–0.91; +- the cases above 1.05 are short lines of two or three stages: `Float32` 20, 32, 96 and 100 (1.09–1.10), + 37 rows of 64 (1.13), `rfft` of 64 (1.07 and 1.09), and `Float64` 12 × 7 (1.05). The stage loop runs + about 50 instructions per stage around the codelet (no mispredicted branches), 2.4 ns per stage + against direct calls of the two kernels of `Float32` 32. +- The machine was shared during the measurement: over the 8 rounds the slowest time of the compiled + chain over its fastest is 1.04 in the median over the cases, 1.19 at the 90th percentile and 1.55 at + most. + +Not verified on aarch64. + +What the single form removes: 583 compiled chain methods and their generators, the dispatch trees from +an `RChain` of a short length into the compiled chains, and the compilation at the first plan of a +short length with a radix sequence outside the compiled set (a short dimension with a large stride). +The package has 1796 methods instead of 2385, the package image on this machine is 61 MB instead +of 72 MB, and `libfftl.so` (x86-64-v3, `--trim=safe`) is 49.5 MB instead of 55.6 MB. + +**Replaced.** The two chain forms (2026-11-09: compiled `Chain` for smooth lengths up to 128, `RChain` +for the rest), and before that the idea of removing `Chain` without changing the stage loop +(2026-11-08) and the separate JIT and run-time plan families (2026-11-02). ## Kernels of a fixed set are plain methods, not `@generated` -**Decision (2026-11-10).** The stage kernels of the codelet set that `RChain` dispatches to, and the -compiled chains of the short lengths, are plain methods defined in `src/kernel_methods.jl` from the -bodies the `@generated` methods build. The `@generated` methods remain as the fallback. +**Decision (2026-11-10).** The stage kernels of the codelet set that `RChain` dispatches to are plain +methods defined in `src/kernel_methods.jl` from the bodies the `@generated` methods build. The `@generated` methods remain as the fallback. **Why.** A package image keeps the complete optimised IR of every precompiled `@generated` specialisation, 130 MB of the 146 MB image. See [Package size and load time](@ref "Package size and load time"). @@ -157,9 +177,8 @@ and matches the direct DFT to 3·10⁻¹⁴ when called from C. Planning on 64 `d = 1` and 0.9 µs for `d = 2` on the M1 (the specialised plan: 4.7 µs and 2.8 µs). Applying the 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. +**Costs.** For real lengths with large prime factors, the run-time kernel is a Bluestein 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/image-size.md b/docs/src/explanation/image-size.md index 3b73cd81..0e4e78b8 100644 --- a/docs/src/explanation/image-size.md +++ b/docs/src/explanation/image-size.md @@ -43,10 +43,9 @@ later Julia versions, so re-run the count below after an upgrade. ## What keeps it small 1. **Plain methods for the fixed set.** The stage kernels that `RChain` dispatches to (radix × sign × - vector type, empty zero-input mask or full output mask) and the compiled chains of the smooth lengths - up to `SMALL_N` are defined in `src/kernel_methods.jl` with `@eval`, from the same body-building + vector type) are defined in `src/kernel_methods.jl` with `@eval`, from the same body-building functions the generators call. The `@generated` methods stay as the fallback for every other combination - (other radices, vector types and chains), so behaviour is unchanged. A plain method keeps its (compressed) source + (the lazy radices), so behaviour is unchanged. A plain method keeps its (compressed) source and its machine code, not its optimised IR. 2. **Lazy radices.** The codelets of the radices in `LAZY_RADICES` (17, 19, 23, 29, 31 by default; the preference `lazy_radices`, read in `src/rchain.jl`) are 54 % of the statements of the stage kernels and are rarely needed. They are not @@ -69,6 +68,10 @@ later Julia versions, so re-run the count below after an upgrade. To add a radix, put it in `RADICES`; `kernel_methods.jl` defines its plain methods. To make a radix lazy, add it to the default of the preference in `src/rchain.jl`. +The numbers on this page were measured while the smooth lengths up to 128 had compiled chains of their +own (583 more plain methods). Without them the image on an AMD EPYC 9354 (AVX-512, which has more vector +types than the M1) went from 72 MB to 61 MB; the M1 numbers have not been measured again. + ## What is left The 32 MB consist of about 3.8 MB for the package without the kernels, 8.6 MB for the definitions of the diff --git a/docs/src/explanation/implementation.md b/docs/src/explanation/implementation.md index 65e62587..65078e2b 100644 --- a/docs/src/explanation/implementation.md +++ b/docs/src/explanation/implementation.md @@ -3,18 +3,16 @@ The [algorithms](@ref "Algorithms") explain what is computed. This page explains how the code is organised and why, for someone who wants to change it. -## Specialisation by type, not by configuration +## Specialisation by type for the codelets, data for the plan -A plan is a Julia value, and its type says what it is. A Stockham chain is -`Chain{T,Rs,σ,K}`: the element type `T`, the radix sequence `Rs` as a tuple of integers, the sign -`σ` and the kernel type `K` are type parameters. The compiler therefore sees the complete stage -sequence when it compiles the transform, unrolls the loop over stages where it can, and emits -machine code for exactly that sequence. The price is one compilation per radix sequence, which -is what the [plan cache](@ref "Plans and reuse") and the [precompile workload](@ref "Precompilation") -are for. +A codelet is specialised by type: the stage kernels take the radix, the sign and the vector type as +type parameters, and the straight-line code of each combination is compiled once. The combinations +form a fixed set (radix × sign × vector type × element type), so all of them are compiled with the +package. -The codelets follow the same rule. `@generated` functions build the straight-line code for a -given radix, element type and vector width when they are first called with those types. +The stage sequence of a length is data. A Stockham chain (`RChain`, `src/rchain.jl`) holds its +radix sequence and one record per stage (codelet index, stride, twiddle table), and a loop calls the +codelet of each stage. Nothing is compiled when a plan of a new length is made. ## Codelets on SIMD vectors @@ -56,25 +54,28 @@ buffers per task are the minimum: the Stockham stages ping-pong between them. That buffer ownership is why a plan must not run on two tasks at once. The plan is a container for the buffers, not a stateless function. -## Two forms of the same chain +## The stage loop -There are two representations of a Stockham chain: +`RChain` calls codelets from the fixed set in a plain loop, selected by an integer index. A short line +is a few nanoseconds of arithmetic, so the loop is written to cost little next to it: -- `Chain` (`src/plan.jl`) has the radix sequence in its type and the stage sequence is unrolled at compile - time. It is the fastest form for short lines: a call carries no dispatch. -- `RChain` (`src/rchain.jl`) has the radix sequence as data and calls codelets from a fixed set (radix × - sign × vector type × element type) in a plain loop, selected by an integer index. +- A line that is one stage is one codelet call, without the loop. +- The stage records lie in one flat vector and carry the address of their twiddle table, and the loop + reads them through a raw pointer. A reference loaded from the chain inside the loop would need a GC + frame, which costs more than a short stage. +- The chain is a mutable struct, so it is passed as one pointer; an immutable chain is copied to the + stack of each caller. +- The switch over the codelet index is inlined into the loop: a stage is one jump and one call. +- The rows of a plane run inside one call of the loop function. -The plan chooses when it is built: `Chain` for smooth lengths up to `SMALL_N` = 128, `RChain` for every -other smooth length. A new radix sequence of `Chain` is a new type, so the first plan of such a size is -compilation; `RChain` has no per-size compilation, and its codelet set is finite, so it can be -precompiled completely. On x86 the run-time choice of the codelet costs 9–20 ns per call, which is why -the short lengths keep the compiled form. +The addresses in the stage records are valid for the chain object they were taken for. A chain +remembers its own address at that moment; a copy of a plan, or a plan stored in a package image and +loaded into another process, is at another address and links its records again before it uses them. `RChain` is also what the C library is built from. `juliac --trim` removes everything it cannot resolve -at build time, and a chain whose type is only known at run time cannot be compiled that way. The +at build time and compiles nothing at run time; the fixed codelet set is all a chain of any length needs. The [libfftl how-to](@ref "C library (libfftl)") describes what the library is built from, and -[Design decisions](@ref "Design decisions") gives the measurements behind the split. The stage kernels of +[Design decisions](@ref "Design decisions") gives the measurements behind the single chain form. The stage kernels of the fixed codelet set are plain methods, not `@generated` ones, to keep the package image small; see [Package size and load time](@ref "Package size and load time"). diff --git a/src/interface.jl b/src/interface.jl index 7f68cb26..afe5ba75 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -116,7 +116,7 @@ check_size(p, x) = size(x) == p.sz || # Dimensions 1 and 2 are fused per plane when a plane fits in cache (`run_planes!`). function apply_dims!(y::StridedArray{Complex{T}}, x::StridedArray{Complex{T}}, p::FFTLPlan{T}, ks::Tuple, ds::Tuple) where {T} isempty(ks) && return nothing - if length(ks) >= 2 && ds[1] == 1 && ds[2] == 2 && ks[1] isa AnyChain && ks[2] isa AnyChain + if length(ks) >= 2 && ds[1] == 1 && ds[2] == 2 && ks[1] isa RChain && ks[2] isa RChain n1, n2 = p.sz[1], p.sz[2] outer = dim_extent(p.sz, 2)[2] if fuse_planes(n1, n2, outer, p.nthreads) @@ -193,7 +193,7 @@ function run_strided_dim!(py::Ptr{T}, px::Ptr{T}, p, k::Kernel{T}, d::Int, groups = CartesianIndices(ntuple(j -> j <= m || j == d ? 1 : p.sz[j], N)) # LinesMT owns a kernel and buffers per task; a composite kernel holds at most k.B lanes. nt = k isa LinesMT ? min(p.nthreads, length(k.kernels)) : - k isa AnyChain && threaded(lanes, n, length(groups), p.nthreads) && + k isa RChain && threaded(lanes, n, length(groups), p.nthreads) && length(p.tbuf[1]) >= n ? p.nthreads : 1 nw = nt == 1 ? length(p.work1) : length(p.tbuf[1]) B = k isa LinesMT ? min(lanes, k.B) : min(lanes, lane_block(lanes, n, nt), nw ÷ n) diff --git a/src/kernel_methods.jl b/src/kernel_methods.jl index 5aa839c0..5a1fe9d0 100644 --- a/src/kernel_methods.jl +++ b/src/kernel_methods.jl @@ -1,5 +1,5 @@ # Compiled-ahead methods of the stage kernels for the fixed codelet set (radix × sign × vector -# type) that `RChain` dispatches to, and of the chains of the small lengths (rchain.jl). The +# type) that `RChain` dispatches to. The # `@generated` kernels in kernels.jl define every other combination. A package image keeps the # complete optimised IR of each precompiled `@generated` specialization but none of a plain # method, so the fixed set is defined as plain methods with the same bodies. The codelets of @@ -48,15 +48,3 @@ for T in (Float32, Float64), σ in (-1, 1), R in RADICES b = kernel_subst(stage_small_body(T, R, σ, VT), Dict{Symbol,Any}(:T => T, :VT => VT)) @eval $(Expr(:function, h, b)) end - -# The compiled chains of the small lengths (rchain.jl) run through `run!`/`run_lanes!` methods of -# plain definition for the same reason; any other chain type uses the generic methods. -for T in (Float32, Float64), (n, σ) in SMALL_ENTRIES - Rs = Tuple(radices(n, 1; esize = sizeof(Complex{T}))) - C = Chain{T,Rs,σ,length(Rs)} - @eval $(Expr(:function, :(run!(py::Ptr{$T}, px::Ptr{$T}, pw1::Ptr{$T}, pw2::Ptr{$T}, P::$C, S::Int)), - run_body(T, Rs, σ, length(Rs)))) - @eval $(Expr(:function, :(run_lanes!(py::Ptr{$T}, px::Ptr{$T}, pw1::Ptr{$T}, pw2::Ptr{$T}, P::$C, - Sx::Int, Sy::Int, B::Int)), - run_lanes_body(T, Rs, σ, length(Rs)))) -end diff --git a/src/kernels.jl b/src/kernels.jl index 46b26928..ef34bede 100644 --- a/src/kernels.jl +++ b/src/kernels.jl @@ -46,23 +46,6 @@ end stage_body(T, R, σ, VT) end -# Stage for strides that are not a multiple of the complex lanes of VT (odd batch sizes): the -# first s - s % C lanes in vectors of C complex values, the scalar remainder one complex value -# at a time. -@inline function stage_split!(py::Ptr{T}, px::Ptr{T}, pt::Ptr{T}, m::Int, s::Int, - ::Val{R}, ::Val{σ}, ::Type{VT}) where {T,R,σ,VT} - C = ncomplex(VT) - rem = s % C - if rem == 0 - stage!(py, px, pt, m, s, Val(R), Val(σ), VT) - else - stage_range!(py, px, pt, m, s, 0, s - rem, 0, m, Val(R), Val(σ), VT) - stage_range!(py, px, pt, m, s, s - rem, s, 0, m, Val(R), Val(σ), Vec{2,T}) - end -end - -# Stage over a q range: the kernel behind the split of a single line, whose q extent is known -# at generation time so the vectorised range and the scalar tail are emitted directly. # base + i·2^e as one scaled-index address computation (a GEP over an element of 2^e bytes), so # that the scaling stays in the addressing mode instead of becoming a value of its own. for (e, ty) in enumerate(("i8", "i16", "i32", "i64")) @@ -89,6 +72,9 @@ stream_bases(base, r, R, first) = stream_strides(r, R, unit) = vcat(:($(Symbol(r, 1)) = $unit), [:($(Symbol(r, o)) = $o * $(Symbol(r, 1))) for o in (3, 5, 7) if o < R]) +# Stage over a q range: a stride that is not a multiple of the complex lanes of VT runs the +# first s - s % C lanes in vectors of C complex values and the remainder one complex value at a +# time. function stage_range_body(::Type{T}, R::Int, σ::Int, ::Type{VT}) where {T,VT} stmts, a, X = codelet(R, σ) C = ncomplex(VT); B = 2sizeof(T) @@ -147,20 +133,6 @@ end stage_last_body(T, R, σ, VT) end -# Last stage for strides that are not a multiple of the complex lanes of VT (odd batch sizes): -# the first s - s % C lanes in vectors, the scalar remainder one complex value at a time. -@inline function stage_last_split!(py::Ptr{T}, px::Ptr{T}, s::Int, - ::Val{R}, ::Val{σ}, ::Type{VT}) where {T,R,σ,VT} - C = ncomplex(VT) - rem = s % C - if rem == 0 - stage_last!(py, px, s, Val(R), Val(σ), VT) - else - stage_last!(py, px, s, 0, s - rem, Val(R), Val(σ), VT) - stage_last!(py, px, s, s - rem, s, Val(R), Val(σ), Vec{2,T}) - end -end - # First stage with s = 1 for vectors of C = 2 or 4 complex values (m a multiple of C): # butterflies p, …, p+C-1 share a register. Groups of C outputs X_k…X_(k+C-1) are transposed # (complex value l of each output into one vector) so that each store writes C consecutive diff --git a/src/pass.jl b/src/pass.jl index f7843c54..3945bafb 100644 --- a/src/pass.jl +++ b/src/pass.jl @@ -88,7 +88,7 @@ end function pass_build(::Type{K}, v::StridedArray{Complex{T},N}, d::Int, σ::Int, nt::Int, k, lanes::Int, Sd::Int, offs::Vector{Int}) where {K,T,N} sz = size(v); st = strides(v); n = sz[d] - B = d == 1 ? 1 : lane_block(lanes, n, k isa AnyChain ? nt : 1) + B = d == 1 ? 1 : lane_block(lanes, n, k isa RChain ? nt : 1) nw = d == 1 ? n : B * n PassPlan{T,N,K}(sz, st, d, σ, k, lanes, Sd, offs, B, nt, thread_buffers(T, nw, 2nt)) end @@ -165,14 +165,14 @@ function run_pass!(py::Ptr{T}, px::Ptr{T}, p::PassPlan{T}) where {T} run_fourstep!(py + off, px + off, k, p.nthreads) elseif p.d == 1 # lines in groups on the tasks (kernels with buffers of their own on one task) - each_line(length(p.offs), p.sz[1], k isa AnyChain ? p.nthreads : 1) do c, i + each_line(length(p.offs), p.sz[1], k isa RChain ? p.nthreads : 1) do c, i local off = p.offs[i + 1] * CT run!(py + off, px + off, cptr(p.bufs[2c-1]), cptr(p.bufs[2c]), k, 1) end else B = p.B; nb = cld(p.lanes, B) nitems = nb * length(p.offs) - ntasks = (k isa AnyChain && p.nthreads > 1 && length(p.offs) * p.lanes * p.sz[p.d] >= TUNING[].thread_min_work) ? p.nthreads : 1 + ntasks = (k isa RChain && p.nthreads > 1 && length(p.offs) * p.lanes * p.sz[p.d] >= TUNING[].thread_min_work) ? p.nthreads : 1 parallel_items(nitems, ntasks) do c, it local g, b; g, b = divrem(it, nb) local i0 = b * B diff --git a/src/plan.jl b/src/plan.jl index f5d1bcd2..eaa09b85 100644 --- a/src/plan.jl +++ b/src/plan.jl @@ -1,20 +1,32 @@ -# 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 +# Kernels for one dimension: a Stockham chain (`RChain`, rchain.jl) 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. abstract type Kernel{T} end """ - Chain{T,Rs,σ,K} + TypedKernels -Stockham chain with radix sequence `Rs` (a tuple, part of the type so that the stage sequence is -compiled into one function) and sign `σ`. +Tag of the plans that carry the types of their kernels in their own type (the default plans), as +opposed to `RChain`, the tag of the run-time plans, whose kernels are of one fixed type. """ -struct Chain{T,Rs,σ,K} <: Kernel{T} - n::Int - tw::NTuple{K,Vector{Complex{T}}} # stage i: w^(kp) at (R-1)p + k - tw1::Vector{Complex{T}} # first stage: w^(kp) at (k-1)m + p -end +abstract type TypedKernels end + +""" + run!(py, px, pw1, pw2, P::Kernel, S) + +Transform S interleaved sequences of length `P.n` from `px` into `py` (which may equal `px`), +using the work buffers `pw1`, `pw2` (S·n complex values each). +""" +function run! end + +""" + run_lanes!(py, px, pw1, pw2, P::Kernel, Sx, Sy, B) + +Transform a block of B lanes; element j of lane i at `px + i + Sx·j` and `py + i + Sy·j` +(py may equal px). The work buffers hold B·n complex values each. +""" +function run_lanes! end # Read-only parts of kernels (chains, four-step twiddles, Bluestein chirp and filter) are shared # between plans and threads: a second plan of a size costs a lookup. Bounded; when full, an @@ -72,21 +84,13 @@ Drop the twiddle and filter tables that FFTL keeps for reuse by later plans. """ clear_cache!() = (@lock KCACHE_LOCK empty!(KCACHE); nothing) -Chain{T}(n::Int, σ::Int, rs::AbstractVector{Int} = radices(n)) where {T} = - cached(() -> build_chain(T, n, σ, rs), CacheKey(:chain, T, n, σ, collect(Int, rs))) - -function build_chain(::Type{T}, n::Int, σ::Int, rs::AbstractVector{Int}) where {T} - tw, tw1 = chain_tables(T, n, σ, rs) - Chain{T,Tuple(rs),σ,length(rs)}(n, Tuple(tw), tw1) -end - # The first of several stages of a single line runs stage_small! (the widest vector, C complex # values, outputs interleaved in registers) where the grouped first stage (vector index group # of the complex counts ncs, 0 for none) would be narrower. small_first(ncs, group::Int, r::Int, n::Int) = r < n && ispow2(r) && (n ÷ r) % ncs[1] == 0 && ncs[1] > (group == 0 ? 1 : ncs[group]) -# Stage twiddles and first-stage table of a radix sequence, shared by Chain and RChain. +# Stage twiddles and first-stage table of a radix sequence. chain_tables(::Type{T}, n::Int, σ::Int, rs::AbstractVector{Int}) where {T} = cached(() -> build_tables(T, n, σ, rs), CacheKey(:tables, T, n, σ, collect(Int, rs)))::Tuple{Vector{Vector{Complex{T}}},Vector{Complex{T}}} @@ -167,136 +171,6 @@ function first_twiddles(::Type{T}, W::Vector{ComplexF64}, r::Int) where {T} t end -# Stage calls of a chain from `src` through `dsts` (one destination per stage). `base` is a -# number of complex values known to divide S (1 if nothing is known); stage i has stride -# S·(R_1⋯R_(i-1)) and takes the widest vector type whose complex count divides base·R_1⋯R_(i-1). -# With `single` (S = 1) the first stage groups 2 or 4 butterflies in one vector when m allows. -# Stage calls with the q (or butterfly) remainder known at generation time: divisible stages -# take the plain kernels (the same calls as before the tail batching), single lines with a -# scalar remainder get the vectorised range plus the scalar tail emitted directly. -# rem === nothing (a runtime stride, e.g. an odd batch size) takes the -# splitting wrappers, which divide at run time. -function mid_call(dst, src, tw, r, σ, V, V1, rem) - rem === nothing && return :(stage_split!($dst, $src, $tw, m, s, Val($r), Val($σ), $V)) - rem == 0 && return :(stage!($dst, $src, $tw, m, s, Val($r), Val($σ), $V)) - quote - stage_range!($dst, $src, $tw, m, s, 0, s - $rem, 0, m, Val($r), Val($σ), $V) - stage_range!($dst, $src, $tw, m, s, s - $rem, s, 0, m, Val($r), Val($σ), $V1) - end -end - -function last_call(dst, src, r, σ, V, V1, rem) - rem === nothing && return :(stage_last_split!($dst, $src, s, Val($r), Val($σ), $V)) - rem == 0 && return :(stage_last!($dst, $src, s, Val($r), Val($σ), $V)) - quote - stage_last!($dst, $src, s, 0, s - $rem, Val($r), Val($σ), $V) - stage_last!($dst, $src, s, s - $rem, s, Val($r), Val($σ), $V1) - end -end - -function group_call(dst, src, tw1, r, σ, GV, V1, rem) - rem === nothing && return :(stage_first_group!($dst, $src, $tw1, m, Val($r), Val($σ), $GV)) - rem == 0 && return :(stage_first_group_range!($dst, $src, $tw1, m, 0, m, Val($r), Val($σ), $GV)) - quote - stage_first_group_range!($dst, $src, $tw1, m, 0, m - $rem, Val($r), Val($σ), $GV) - stage_first_group_1!($dst, $src, $tw1, m, m - $rem, m, Val($r), Val($σ), $V1) - end -end - -function chain_body(T, Rs, σ, src, dsts, base, single) - K = length(Rs) - n = prod(Rs; init = 1) - ladder = vec_ladder(T) - # vector type of a grouped first stage: the widest with 2 or 4 complex values (the scalar - # remainder of m is handled by the kernel) - group = findfirst(V -> ncomplex(V) in (2, 4), ladder) - # L and s are compile-time constants where provable (prod(Rs) == n is validated by - # build_tables; the single dispatch guarantees S == 1), so the remainders of the stages are - # known at generation time. - body = Any[:(L = $n), single ? :(s = 1) : :(s = S)] - sst = 1 - for i in 1:K - r = Rs[i] - dst = dsts[i] - # widest vector type whose complex count fits the stride (single lines take a scalar - # remainder; batched lines need divisibility) - V = base == 1 ? ladder[findfirst(V -> ncomplex(V) <= base * sst, ladder)] : - ladder[findfirst(V -> (base * sst) % ncomplex(V) == 0, ladder)] - # the q remainder of this stage: a compile-time constant for single lines, provably zero - # for batched lines (their stride is a multiple of the vector type), runtime for odd batch - # sizes - rem = single ? sst % ncomplex(V) : (base > 1 ? 0 : nothing) - if single && i == 1 && small_first(map(ncomplex, ladder), something(group, 0), r, n) - push!(body, :(stage_small!($dst, $src, cptr(P.tw1), $(n ÷ r), Val($r), Val($σ), $(ladder[1]))), - :(L = L ÷ $r), :(s *= $r)) - elseif i == K - push!(body, last_call(dst, src, r, σ, V, single_vec(T), rem)) - elseif single && i == 1 && group !== nothing && ncomplex(V) == 1 - # the first stage of a single line: vectorised over the butterflies (scalar remainder - # handled by the kernel) - push!(body, :(m = L ÷ $r)) - GV = ladder[group] - push!(body, group_call(dst, src, :(cptr(P.tw1)), r, σ, GV, single_vec(T), (n ÷ r) % ncomplex(GV))) - push!(body, :(L = m), :(s *= $r)) - else - push!(body, :(m = L ÷ $r)) - push!(body, mid_call(dst, src, :(cptr(P.tw[$i])), r, σ, V, single_vec(T), rem)) - push!(body, :(L = m), :(s *= $r)) - end - src = dst - sst *= r - end - body -end - -# Branches on S: S = 1, then the complex counts of the ladder (widest first) that divide S. -function stride_dispatch(T, body) - ex = body(1, false) - for V in reverse(vec_ladder(T)) - C = ncomplex(V) - C == 1 && continue - ex = :(S % $C == 0 ? $(body(C, false)) : $ex) - end - :(S == 1 ? $(body(1, true)) : $ex) -end - -""" - run!(py, px, pw1, pw2, P::Kernel, S) - -Transform S interleaved sequences of length `P.n` from `px` into `py` (which may equal `px`), -using the work buffers `pw1`, `pw2` (S·n complex values each). -""" -function run_body(::Type{T}, Rs, σ::Int, K::Int) where {T} - K == 0 && return :(py == px || unsafe_copyto!(py, px, 2S); nothing) - # Out of place: the stages alternate between py and pw1 and end in py; px is only read. - oop = [iseven(K - i) ? :py : :pw1 for i in 1:K] - # In place: the first stage reads px, the last writes px, the ones between alternate between - # the two work buffers. A single stage (m = 1) loads all R inputs of a butterfly before it - # stores, so it works on px directly. - ip = Any[:px for _ in 1:K] - for i in K-1:-1:1 - ip[i] = isodd(K - i) ? :pw1 : :pw2 - end - dispatch(dsts) = stride_dispatch(T, (base, single) -> quote - $(chain_body(T, Rs, σ, :px, dsts, base, single)...) - end) - inplace = dispatch(ip) - quote - GC.@preserve P begin - if py != px - $(dispatch(oop)) - else - $inplace - end - end - nothing - end -end - -@generated function run!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::Chain{T,Rs,σ,K}, S::Int) where {T,Rs,σ,K} - run_body(T, Rs, σ, K) -end - """ Bluestein{T} @@ -378,8 +252,8 @@ end # 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} +# C: TypedKernels, or RChain for run-time plans (their convolution transforms are chains) +function Bluestein{T}(n::Int, σ::Int, ::Type{C} = TypedKernels) 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) @@ -393,7 +267,7 @@ 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} +function convolution_transforms(::Type{T}, M::Int, ::Type{TypedKernels}) 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 @@ -565,68 +439,12 @@ run!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::Composite{T}, S::Int) run_lanes!(py, px, pw1, pw2, P, S, S, S) 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, σ) Composite{T}(n, s, q, σ, S) end -""" - run_lanes!(py, px, pw1, pw2, P::Chain, Sx, Sy, B) - -Transform a block of B lanes; element j of lane i at `px + i + Sx·j` and `py + i + Sy·j` -(py may equal px). The work buffers hold B·n complex values each. -""" -function run_lanes_body(::Type{T}, Rs, σ::Int, K::Int) where {T} - function body(V) - ex = Any[:(L = P.n), :(s0 = 1)] - src, srcS = :px, :Sx - for i in 1:K - r = Rs[i] - last = i == K - dst, dstS = last ? (:py, :Sy) : (isodd(i) ? :pw1 : :pw2, :B) - if last - push!(ex, :(stage_last_lanes!($dst, $src, s0, $srcS, $dstS, B, Val($r), Val($σ), $V))) - else - push!(ex, :(m = L ÷ $r)) - push!(ex, :(stage_lanes!($dst, $src, cptr(P.tw[$i]), m, s0, $srcS, $dstS, B, Val($r), Val($σ), $V))) - push!(ex, :(L = m), :(s0 *= $r)) - end - src, srcS = dst, dstS - end - ex - end - K == 0 && return :(lane_copy!(py, px, Sy, Sx, B, 1); nothing) - ladder = vec_ladder(T) - # widest vector type whose complex count divides B - function by_width(f) - ex = f(ladder[end]) - for V in reverse(ladder[1:end-1]) - ex = :(B % $(ncomplex(V)) == 0 ? $(f(V)) : $ex) - end - ex - end - if K == 1 - # a single stage loads all R inputs of a butterfly before it stores: py may equal px - r = Rs[1] - single(V) = :(stage_last_lanes!(py, px, 1, Sx, Sy, B, Val($r), Val($σ), $V)) - return quote - GC.@preserve P $(by_width(single)) - nothing - end - end - quote - GC.@preserve P $(by_width(V -> quote $(body(V)...) end)) - nothing - end -end - -@generated function run_lanes!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::Chain{T,Rs,σ,K}, - Sx::Int, Sy::Int, B::Int) where {T,Rs,σ,K} - run_lanes_body(T, Rs, σ, K) -end - # y[i + Sy·j] = x[i + Sx·j] for B lanes and n elements function lane_copy!(py::Ptr{T}, px::Ptr{T}, Sy::Int, Sx::Int, B::Int, n::Int) where {T} cy = Ptr{Complex{T}}(py); cx = Ptr{Complex{T}}(px) diff --git a/src/precompile.jl b/src/precompile.jl index 89407660..c563acc8 100644 --- a/src/precompile.jl +++ b/src/precompile.jl @@ -7,8 +7,8 @@ # # "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 -# paths for both the compile-time and the run-time plans. +# codelets (radices 4, 5, 8, 9, 10 and 16, plus Bluestein) and the +# plan types of both the default 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. # "none" no workload at all (for CI or when the image size is more important diff --git a/src/r2r.jl b/src/r2r.jl index 7a68e3c6..4fb5f8e2 100644 --- a/src/r2r.jl +++ b/src/r2r.jl @@ -72,7 +72,7 @@ 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}} -function R2RLine{T}(kind::R2RKind, n::Int, ::Type{C} = Chain) where {T,C} +function R2RLine{T}(kind::R2RKind, n::Int, ::Type{C} = TypedKernels) where {T,C} n >= 1 || throw(ArgumentError("length must be positive")) kind == REDFT00 && n < 2 && throw(ArgumentError("REDFT00 needs length ≥ 2")) base = kind in (RODFT10,) ? REDFT10 : kind == RODFT01 ? REDFT01 : kind == RODFT11 ? REDFT11 : kind diff --git a/src/rchain.jl b/src/rchain.jl index 583bee35..3260364d 100644 --- a/src/rchain.jl +++ b/src/rchain.jl @@ -1,5 +1,4 @@ -# Stockham chain with the radix sequence as run-time data. Same stages, twiddle tables, buffer -# schedule and vector-type rules as `Chain` (plan.jl), but the stage sequence is a loop over +# Stockham chain with the radix sequence as run-time data: the stage sequence is a loop over # plain code that calls the codelet of each stage from a fixed set (radix × sign × vector type), # selected by an integer index. Nothing is compiled per length, and the execution path contains # no call into a generated function whose body calls further generated functions, so it can be @@ -43,9 +42,6 @@ using Preferences: @load_preference const LAZY_RADICES = Tuple(Int.(@load_preference("lazy_radices", [17, 19, 23, 29, 31]))) all(in(RADICES), LAZY_RADICES) || error("LAZY_RADICES must be radices from RADICES") -# Whether the factorisation of n has a prime factor whose codelet is compiled on first use. -has_lazy_factor(n::Int) = any(in(LAZY_RADICES), keys(factorize(n)[1])) - # Codelet dispatch: one function per kind of stage call, a switch over the codelet index, # generated for both element types so that every call has constant type arguments. Codelets of # LAZY_RADICES are not in the switch; `_lazy` takes their index and calls the stage kernel @@ -120,13 +116,22 @@ const OP_LAST_TAIL = 4 const OP_GROUP = 5 # first stage vectorised over the butterflies p const OP_SMALL = 6 # first stage over the butterflies p on the widest vector (stage_small!) -struct StageCall +struct StageCall{T} op::Int c::Int # codelet of the vectorised part c1::Int # codelet of the scalar remainder (one complex value per vector) m::Int s::Int rem::Int + pt::Ptr{T} # twiddle table of the stage (C_NULL for a last stage), see link_calls! +end + +# One stage of a batch of lines or of a lane block: the vector type is chosen per call from the +# stride, so the record holds the codelet index of the widest vector type (the index of vector +# type vi is c + (vi - 1)·length(RADICES)) and the twiddle table. +struct BatchStage{T} + c::Int + pt::Ptr{T} # C_NULL for the last stage end """ @@ -134,28 +139,56 @@ end Stockham chain of length `n` and sign `σ` with the radix sequence `rs` stored as a vector. """ -struct RChain{T} <: Kernel{T} - n::Int - σ::Int - rs::Vector{Int} - tw::Vector{Vector{Complex{T}}} # stage i: w^(kp) at (R-1)p + k - tw1::Vector{Complex{T}} # first stage: w^(kp) at (k-1)m + p - single::Vector{StageCall} # the stage calls of a single line +mutable struct RChain{T} <: Kernel{T} + const n::Int + const σ::Int + const rs::Vector{Int} + const tw::Vector{Vector{Complex{T}}} # stage i: w^(kp) at (R-1)p + k + const tw1::Vector{Complex{T}} # first stage: w^(kp) at (k-1)m + p + const single::Vector{StageCall{T}} # the stage calls of a single line + const batch::Vector{BatchStage{T}} # the stages of a batch or a lane block + const one::Int # codelet of a single line that is one stage, else 0 + linked::UInt # address of the chain when its stage records were linked end -# Either kind of Stockham chain: the compiled Chain (short smooth lengths) or the run-time RChain. -const AnyChain = Union{Chain,RChain} - RChain{T}(n::Int, σ::Int, rs::AbstractVector{Int} = radices(n)) where {T} = cached(() -> build_rchain(T, n, σ, rs), CacheKey(:rchain, T, n, σ, collect(Int, rs)))::RChain{T} function build_rchain(::Type{T}, n::Int, σ::Int, rs::AbstractVector{Int}) where {T} all(in(RADICES), rs) || throw(ArgumentError("radix without a codelet in RADICES")) tw, tw1 = chain_tables(T, n, σ, rs) - RChain{T}(n, σ, collect(Int, rs), tw, tw1, single_schedule(T, n, σ, collect(Int, rs))) + calls = single_schedule(T, n, σ, collect(Int, rs)) + batch = [BatchStage{T}(codelet_index(T, σ, 1, r), Ptr{T}(C_NULL)) for r in rs] + P = RChain{T}(n, σ, collect(Int, rs), tw, tw1, calls, batch, length(calls) == 1 ? calls[1].c : 0, UInt(0)) + link_calls!(P) + P +end + +# The stage records carry the address of their twiddle table, so that a stage loop reads one +# flat vector and nothing else (a load through the chain inside the loop needs a GC frame, which +# costs more than a short stage). The tables belong to the chain and are never resized. The +# addresses are valid for the chain object they were taken for: a copy, or a chain that was stored +# in a package image and loaded into another process, is at another address than `linked` and is +# linked again before its records are used. +@inline linked(P::RChain) = P.linked == UInt(pointer_from_objref(P)) + +@noinline function link_calls!(P::RChain{T}) where {T} + calls = P.single + for i in eachindex(calls) + st = calls[i] + pt = st.op == OP_GROUP || st.op == OP_SMALL ? cptr(P.tw1) : + st.op == OP_MID || st.op == OP_MID_TAIL ? cptr(P.tw[i]) : Ptr{T}(C_NULL) + calls[i] = StageCall{T}(st.op, st.c, st.c1, st.m, st.s, st.rem, pt) + end + batch = P.batch + for i in 1:length(batch)-1 + batch[i] = BatchStage{T}(batch[i].c, cptr(P.tw[i])) + end + P.linked = UInt(pointer_from_objref(P)) + nothing end -# The stage calls of a single line, with the rules of `chain_body` for S = 1: stage i takes the +# The stage calls of a single line: stage i takes the # widest vector type whose complex count fits its stride s = R_1⋯R_(i-1), the remainder # s mod C runs one complex value at a time; the first stage is grouped over the butterflies # when its stride is 1, with stage_small! where that is wider. @@ -169,7 +202,7 @@ function single_schedule(::Type{T}, n::Int, σ::Int, rs::Vector{Int}) where {T} break end end - calls = StageCall[] + calls = StageCall{T}[] K = length(rs) L = n; s = 1 for i in 1:K @@ -182,13 +215,13 @@ function single_schedule(::Type{T}, n::Int, σ::Int, rs::Vector{Int}) where {T} c = codelet_index(T, σ, vi, r); c1 = codelet_index(T, σ, nv, r) m = L ÷ r if i == 1 && small_first(ncs, group, r, n) - push!(calls, StageCall(OP_SMALL, codelet_index(T, σ, 1, r), c1, m, s, 0)) + push!(calls, StageCall{T}(OP_SMALL, codelet_index(T, σ, 1, r), c1, m, s, 0, Ptr{T}(C_NULL))) elseif i == K - push!(calls, StageCall(rem == 0 ? OP_LAST : OP_LAST_TAIL, c, c1, m, s, rem)) + push!(calls, StageCall{T}(rem == 0 ? OP_LAST : OP_LAST_TAIL, c, c1, m, s, rem, Ptr{T}(C_NULL))) elseif i == 1 && group != 0 && ncs[vi] == 1 - push!(calls, StageCall(OP_GROUP, codelet_index(T, σ, group, r), c1, m, s, 0)) + push!(calls, StageCall{T}(OP_GROUP, codelet_index(T, σ, group, r), c1, m, s, 0, Ptr{T}(C_NULL))) else - push!(calls, StageCall(rem == 0 ? OP_MID : OP_MID_TAIL, c, c1, m, s, rem)) + push!(calls, StageCall{T}(rem == 0 ? OP_MID : OP_MID_TAIL, c, c1, m, s, rem, Ptr{T}(C_NULL))) end L = m; s *= r end @@ -201,21 +234,23 @@ end @inline stage_dst(py, px, pw1, pw2, i, K, inplace) = inplace ? (i == K ? px : (isodd(K - i) ? pw1 : pw2)) : (iseven(K - i) ? py : pw1) -function run!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::RChain{T}, S::Int) where {T} - K = length(P.rs) +# The single line is kept apart from the batch: the batch needs a GC frame, which a short single +# line must not pay for. +@inline run!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::RChain{T}, S::Int) where {T} = + S == 1 ? run_single!(py, px, pw1, pw2, P) : run_batch!(py, px, pw1, pw2, P, S) + +@noinline function run_batch!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::RChain{T}, S::Int) where {T} + K = length(P.batch) if K == 0 py == px || unsafe_copyto!(py, px, 2S) return nothing end - if P.n <= SMALL_N && SMALL_INDEX[P.n] > 0 - i = SMALL_INDEX[P.n] + (P.σ == 1 ? 1 : 0) - small_run!(py, px, pw1, pw2, P, i, S) && return nothing - end - S == 1 && return run_single!(py, px, pw1, pw2, P) + linked(P) || link_calls!(P) + pb = pointer(P.batch) ncs = ladder_ncomplex(T) nv = length(ncs) - σ = P.σ - # base: the widest complex count that divides S (as stride_dispatch in plan.jl) + NR = length(RADICES) + # base: the widest complex count that divides S base = 1 for C in ncs if C > 1 && S % C == 0 @@ -224,221 +259,172 @@ function run!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::RChain{T}, S: end end inplace = py == px - GC.@preserve P begin - L = P.n; s = S; sst = 1; src = px - for i in 1:K - r = P.rs[i] - dst = stage_dst(py, px, pw1, pw2, i, K, inplace) - # base > 1: the widest vector type dividing base·sst (no remainder); base = 1 (odd S): - # the widest one that fits sst, with the remainder of s split off at run time - vi = 1 - if base > 1 - while (base * sst) % ncs[vi] != 0 - vi += 1 - end + L = P.n; s = S; sst = 1; src = px + @inbounds for i in 1:K + st = unsafe_load(pb, i) + dst = stage_dst(py, px, pw1, pw2, i, K, inplace) + # base > 1: the widest vector type dividing base·sst (no remainder); base = 1 (odd S): + # the widest one that fits sst, with the remainder of s split off at run time + vi = 1 + if base > 1 + while (base * sst) % ncs[vi] != 0 + vi += 1 + end + else + while ncs[vi] > sst + vi += 1 + end + end + c = st.c + (vi - 1) * NR + rem = s % ncs[vi] + if i == K + r = L + if rem == 0 + cstage_last!(dst, src, s, c) else - while ncs[vi] > sst - vi += 1 - end + cstage_last_range!(dst, src, s, 0, s - rem, c) + cstage_last_range!(dst, src, s, s - rem, s, st.c + (nv - 1) * NR) end - c = codelet_index(T, σ, vi, r) - rem = s % ncs[vi] - if i == K - if rem == 0 - cstage_last!(dst, src, s, c) - else - cstage_last_range!(dst, src, s, 0, s - rem, c) - cstage_last_range!(dst, src, s, s - rem, s, codelet_index(T, σ, nv, r)) - end + else + r = RADICES[(st.c - 1) % NR + 1] + m = L ÷ r + if rem == 0 + cstage!(dst, src, st.pt, m, s, c) else - m = L ÷ r - pt = cptr(P.tw[i]) - if rem == 0 - cstage!(dst, src, pt, m, s, c) - else - cstage_range!(dst, src, pt, m, s, 0, s - rem, c) - cstage_range!(dst, src, pt, m, s, s - rem, s, codelet_index(T, σ, nv, r)) - end - L = m; s *= r + cstage_range!(dst, src, st.pt, m, s, 0, s - rem, c) + cstage_range!(dst, src, st.pt, m, s, s - rem, s, st.c + (nv - 1) * NR) end - src = dst - sst *= r + L = m; s *= r end + src = dst + sst *= r end nothing end -function run_single!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::RChain{T}) where {T} - calls = P.single - K = length(calls) +# A single line. One stage is one codelet call; several stages run the stage loop. +@inline function run_single!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::RChain{T}) where {T} + c = P.one + if c != 0 + cstage_last!(py, px, 1, c) + return nothing + end + run_calls!(py, px, pw1, pw2, P) +end + +# The stage loop of a single line. The calls are read through a raw pointer, so that no reference +# is live across a stage (the chain itself is kept alive by the caller); a reference loaded from +# the chain here or passed in by the caller costs a GC frame. The codelet switches are inlined: a +# stage is one jump and one call. +@inline function stage_loop!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, pc::Ptr{StageCall{T}}, K::Int) where {T} inplace = py == px - GC.@preserve P begin - src = px - for i in 1:K - st = calls[i] - dst = stage_dst(py, px, pw1, pw2, i, K, inplace) - op = st.op - if op == OP_MID - cstage!(dst, src, cptr(P.tw[i]), st.m, st.s, st.c) - elseif op == OP_GROUP - cstage_group!(dst, src, cptr(P.tw1), st.m, st.c) - elseif op == OP_SMALL - cstage_small!(dst, src, cptr(P.tw1), st.m, st.c) - elseif op == OP_LAST - cstage_last!(dst, src, st.s, st.c) - elseif op == OP_MID_TAIL - pt = cptr(P.tw[i]) - cstage_range!(dst, src, pt, st.m, st.s, 0, st.s - st.rem, st.c) - cstage_range!(dst, src, pt, st.m, st.s, st.s - st.rem, st.s, st.c1) - else - cstage_last_range!(dst, src, st.s, 0, st.s - st.rem, st.c) - cstage_last_range!(dst, src, st.s, st.s - st.rem, st.s, st.c1) - end - src = dst + src = px + for i in 1:K + st = unsafe_load(pc, i) + dst = stage_dst(py, px, pw1, pw2, i, K, inplace) + op = st.op + if op == OP_LAST + @inline cstage_last!(dst, src, st.s, st.c) + elseif op == OP_MID + @inline cstage!(dst, src, st.pt, st.m, st.s, st.c) + elseif op == OP_SMALL + @inline cstage_small!(dst, src, st.pt, st.m, st.c) + elseif op == OP_GROUP + @inline cstage_group!(dst, src, st.pt, st.m, st.c) + elseif op == OP_MID_TAIL + @inline cstage_range!(dst, src, st.pt, st.m, st.s, 0, st.s - st.rem, st.c) + @inline cstage_range!(dst, src, st.pt, st.m, st.s, st.s - st.rem, st.s, st.c1) + else + @inline cstage_last_range!(dst, src, st.s, 0, st.s - st.rem, st.c) + @inline cstage_last_range!(dst, src, st.s, st.s - st.rem, st.s, st.c1) + end + src = dst + end + nothing +end + +# The calls of a chain, linked. +@inline function linked_calls(P::RChain{T}) where {T} + linked(P) || link_calls!(P) + pointer(P.single) +end + +@noinline function run_calls!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::RChain{T}) where {T} + K = length(P.single) + if K == 0 + py == px || unsafe_copyto!(py, px, 2) + return nothing + end + stage_loop!(py, px, pw1, pw2, linked_calls(P), K) +end + +# `count` single lines, `bytes` apart. +@noinline function run_calls_rows!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::RChain{T}, count::Int, bytes::Int) where {T} + K = length(P.single) + if K == 0 + for j in 0:count-1 + py == px || unsafe_copyto!(py + j * bytes, px + j * bytes, 2) end + return nothing + end + pc = linked_calls(P) + for j in 0:count-1 + stage_loop!(py + j * bytes, px + j * bytes, pw1, pw2, pc, K) end nothing end function run_lanes!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::RChain{T}, Sx::Int, Sy::Int, B::Int) where {T} - rs = P.rs; K = length(rs); σ = P.σ + K = length(P.batch) if K == 0 lane_copy!(py, px, Sy, Sx, B, 1) return nothing end - if P.n <= SMALL_N && SMALL_INDEX[P.n] > 0 - i = SMALL_INDEX[P.n] + (P.σ == 1 ? 1 : 0) - small_lanes_run!(py, px, pw1, pw2, P, i, Sx, Sy, B) && return nothing - end + linked(P) || link_calls!(P) + pb = pointer(P.batch) ncs = ladder_ncomplex(T) # widest vector type whose complex count divides B, else the narrowest vi = length(ncs) - for j in 1:length(ncs)-1 + @inbounds for j in 1:length(ncs)-1 if B % ncs[j] == 0 vi = j break end end - GC.@preserve P begin - if K == 1 # one stage loads all inputs of a butterfly before it stores: py may equal px - cstage_last_lanes!(py, px, 1, Sx, Sy, B, codelet_index(T, σ, vi, rs[1])) - return nothing - end - L = P.n; s0 = 1; src = px; srcS = Sx - for i in 1:K - r = rs[i] - c = codelet_index(T, σ, vi, r) - if i == K - cstage_last_lanes!(py, src, s0, srcS, Sy, B, c) - else - dst = isodd(i) ? pw1 : pw2 - m = L ÷ r - cstage_lanes!(dst, src, cptr(P.tw[i]), m, s0, srcS, B, B, c) - L = m; s0 *= r - src = dst; srcS = B - end - end - end - nothing -end - -# Small lengths: the stage loop of a single line runs as compiled code for each smooth length up -# to SMALL_N and both signs. The entries of a length are consecutive (σ = -1, then σ = +1); one -# table lookup per call finds them. The compiled chain uses the tables of the RChain, and its -# radix sequence is checked against the compiled one, so any other plan takes the stage loop. -const SMALL_N = 128 -const SMALL_ENTRIES = [(n, σ) for n in 2:SMALL_N if smooth(n) && !has_lazy_factor(n) for σ in (-1, 1)] -const SMALL_INDEX = let t = zeros(Int, SMALL_N) - for (i, (n, σ)) in enumerate(SMALL_ENTRIES) - σ == -1 && (t[n] = i) - end - t -end - -@noinline function small_chain!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::RChain{T}, S::Int, - ::Type{C}, ::Val{K}) where {T,C<:Chain{T},K} - c = C(P.n, ntuple(i -> P.tw[i], Val(K)), P.tw1) - GC.@preserve P run!(py, px, pw1, pw2, c, S) -end - -@noinline function small_lanes_chain!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::RChain{T}, - Sx::Int, Sy::Int, B::Int, ::Type{C}, ::Val{K}) where {T,C<:Chain{T},K} - c = C(P.n, ntuple(i -> P.tw[i], Val(K)), P.tw1) - GC.@preserve P run_lanes!(py, px, pw1, pw2, c, Sx, Sy, B) -end - -# Binary tree of compares over the entry index (a flat chain of ~120 compares costs ~50 ns per -# call). Each leaf checks the radix sequence of the RChain against the compiled one. `leaf` gives -# the statement of a leaf for the chain type C with K stages. -function small_tree(T, lo, hi, leaf) - if lo == hi - n, σ = SMALL_ENTRIES[lo] - rs = Tuple(radices(n, 1; esize = sizeof(Complex{T}))) - K = length(rs) - same = Any[:(i == $lo), :(length(P.rs) == $K)] - for (j, r) in enumerate(rs) - push!(same, :(P.rs[$j] == $r)) - end - return :(if $(Expr(:&&, same...)) - $(leaf(:(Chain{$T,$rs,$σ,$K}), K)) - return true - end) + NR = length(RADICES) + off = (vi - 1) * NR + if K == 1 # one stage loads all inputs of a butterfly before it stores: py may equal px + cstage_last_lanes!(py, px, 1, Sx, Sy, B, unsafe_load(pb).c + off) + return nothing end - mid = (lo + hi) ÷ 2 - :(if i <= $mid - $(small_tree(T, lo, mid, leaf)) + L = P.n; s0 = 1; src = px; srcS = Sx + for i in 1:K + st = unsafe_load(pb, i) + c = st.c + off + if i == K + cstage_last_lanes!(py, src, s0, srcS, Sy, B, c) else - $(small_tree(T, mid + 1, hi, leaf)) - end) -end - -for T in (Float32, Float64) - @eval function small_run!(py::Ptr{$T}, px::Ptr{$T}, pw1::Ptr{$T}, pw2::Ptr{$T}, P::RChain{$T}, i::Int, S::Int) - $(small_tree(T, 1, length(SMALL_ENTRIES), - (C, K) -> :(small_chain!(py, px, pw1, pw2, P, S, $C, Val($K)))) - ) - false - end - @eval function small_lanes_run!(py::Ptr{$T}, px::Ptr{$T}, pw1::Ptr{$T}, pw2::Ptr{$T}, P::RChain{$T}, - i::Int, Sx::Int, Sy::Int, B::Int) - $(small_tree(T, 1, length(SMALL_ENTRIES), - (C, K) -> :(small_lanes_chain!(py, px, pw1, pw2, P, Sx, Sy, B, $C, Val($K)))) - ) - false - end -end - -# Rows of a plane: `count` single lines of the chain, row j at py + j·rowstride complex values. -# A compiled chain takes all rows in one call: one dispatch per plane instead of one per row. -@noinline function small_rows_chain!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::RChain{T}, - count::Int, rowstride::Int, ::Type{C}, ::Val{K}) where {T,C<:Chain{T},K} - c = C(P.n, ntuple(i -> P.tw[i], Val(K)), P.tw1) - bytes = rowstride * 2 * sizeof(T) - GC.@preserve P for j in 0:count-1 - run!(py + j * bytes, px + j * bytes, pw1, pw2, c, 1) + dst = isodd(i) ? pw1 : pw2 + r = @inbounds RADICES[(st.c - 1) % NR + 1] + m = L ÷ r + cstage_lanes!(dst, src, st.pt, m, s0, srcS, B, B, c) + L = m; s0 *= r + src = dst; srcS = B + end end nothing end -for T in (Float32, Float64) - @eval function small_rows_run!(py::Ptr{$T}, px::Ptr{$T}, pw1::Ptr{$T}, pw2::Ptr{$T}, P::RChain{$T}, - i::Int, count::Int, rowstride::Int) - $(small_tree(T, 1, length(SMALL_ENTRIES), - (C, K) -> :(small_rows_chain!(py, px, pw1, pw2, P, count, rowstride, $C, Val($K)))) - ) - false - end -end - -# Generic: one call of run! per row (any kernel), compiled chains through their dispatch. +# One single line per row. function run_rows!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::RChain{T}, count::Int, rowstride::Int) where {T} - if P.n <= SMALL_N && SMALL_INDEX[P.n] > 0 - i = SMALL_INDEX[P.n] + (P.σ == 1 ? 1 : 0) - small_rows_run!(py, px, pw1, pw2, P, i, count, rowstride) && return nothing - end bytes = rowstride * 2 * sizeof(T) - for j in 0:count-1 - run!(py + j * bytes, px + j * bytes, pw1, pw2, P, 1) + c = P.one + if c != 0 + for j in 0:count-1 + @inline cstage_last!(py + j * bytes, px + j * bytes, 1, c) + end + else + run_calls_rows!(py, px, pw1, pw2, P, count, bytes) end nothing end diff --git a/src/real.jl b/src/real.jl index 399f8402..2e859b18 100644 --- a/src/real.jl +++ b/src/real.jl @@ -28,11 +28,11 @@ end const RKern{T} = Union{RChain{T},RBluestein{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} -ckernel(::Type{T}, n::Int, σ::Int, ::Type{Chain}) where {T} = kernel(T, n, σ) +ckernel(::Type{T}, n::Int, σ::Int, ::Type{TypedKernels}) where {T} = kernel(T, n, σ) ckernel(::Type{T}, n::Int, σ::Int, ::Type{RChain}) where {T} = rkernel(T, n, σ) -# C: Chain (kernel type in the line's type) or RChain (RKern{T}, for run-time plans) -function RealLine{T}(n::Int, σ::Int, ::Type{C} = Chain) where {T,C} +# C: TypedKernels (kernel type in the line's type) or RChain (RKern{T}, for run-time plans) +function RealLine{T}(n::Int, σ::Int, ::Type{C} = TypedKernels) where {T,C} inner = ckernel(T, iseven(n) ? n ÷ 2 : n, σ, C) u = σ == -1 ? ComplexF64(0, -0.5) : ComplexF64(0, 1) w = Complex{T}[] @@ -351,7 +351,7 @@ 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 # 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} = +task_reallines(line::RealLine{T}, n::Int, σ::Int, nt::Int, ::Type{C} = TypedKernels) where {T,C} = own_buffers(line) && nt > 1 ? [line; [RealLine{T}(n, σ, C) for _ in 2:nt]] : [line] rline_tasks(L, nt) = length(L) > 1 ? length(L) : line_tasks(L[1], nt) rline_task(L, c) = L[length(L) > 1 ? c : 1] diff --git a/src/threads.jl b/src/threads.jl index d51772de..7fa8db8b 100644 --- a/src/threads.jl +++ b/src/threads.jl @@ -344,7 +344,7 @@ end 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]] + chains = [c1; [c1 isa RChain ? c1 : first_kernel(T, s, σ, q, 1) for _ in 2:nt]] b1 = Bluestein{T}(q, σ) blues = [b1; [Bluestein{T}(q, σ) for _ in 2:nt]] tw = cached(() -> Complex{T}.(unit_roots(σ, n)), CacheKey(:composite, T, n, σ, Int[]))::Vector{Complex{T}} @@ -487,7 +487,7 @@ end fuse_planes(n1, n2, outer, nt) = n1 * n2 <= TUNING[].plane_budget && (nt == 1 || outer >= nt) -function run_planes!(py::Ptr{T}, px::Ptr{T}, p, k1::Union{Chain{T},RChain{T}}, k2::Union{Chain{T},RChain{T}}, +function run_planes!(py::Ptr{T}, px::Ptr{T}, p, k1::RChain{T}, k2::RChain{T}, n1::Int, n2::Int, outer::Int) where {T} CT = 2 * sizeof(T) step = n1 * n2 * CT @@ -534,7 +534,7 @@ end # Stockham chains: lane blocks when there are enough lanes, slices in groups when there are enough # slices, otherwise one slice after the other. -function run_kernel_dim!(py::Ptr{T}, px::Ptr{T}, p, k::Union{Chain{T},RChain{T}}, S::Int, n::Int, outer::Int) where {T} +function run_kernel_dim!(py::Ptr{T}, px::Ptr{T}, p, k::RChain{T}, S::Int, n::Int, outer::Int) where {T} nt = p.nthreads step = 2 * S * n * sizeof(T) if use_lanes(S, n, outer, nt) diff --git a/test/runtests.jl b/test/runtests.jl index cd5f6849..f5d11e2b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -727,48 +727,77 @@ thin(T, v) = FULL || T == Float64 ? v : v[1:2:end] # quick run: Float64 only (the first RChain call compiles every codelet of its element type) @testset "runtime-radix chains, $T" for T in pick((Float64, Float32), (Float64,)) - # RChain runs the same codelets in the same order as Chain: results must be identical + # single lines out of place and in place, lane blocks and strided batches against the + # direct DFT at a few frequencies; the two forms of each pair run the same codelets p(a) = Ptr{T}(pointer(a)) - for n in pick([1, 2, 7, 8, 16, 60, 100, 128, 243, 1000, 31 * 27, 4096, 2^14, 3^8, 2 * 10007 - 14], - thin(T, [1, 8, 60, 128, 243, 1000, 31 * 27, 2^14])) - rs = FFTL.radices(n, 1; esize = sizeof(Complex{T})) - c = FFTL.build_chain(T, n, -1, rs); rc = FFTL.RChain{T}(n, -1, rs) + at(y, ks, B = 1, i = 0) = [y[i + B * k + 1] for k in ks] + for n in pick([2, 7, 8, 16, 60, 100, 128, 243, 1000, 31 * 27, 4096, 2^14, 3^8, 2 * 10007 - 14], + thin(T, [2, 8, 60, 128, 243, 1000, 31 * 27, 2^14])) + rc = FFTL.RChain{T}(n, -1, FFTL.radices(n, 1; esize = sizeof(Complex{T}))) x = randn(Complex{T}, n); w1 = similar(x); w2 = similar(x) - y1 = similar(x); y2 = similar(x); z1 = copy(x); z2 = copy(x) - GC.@preserve x y1 y2 z1 z2 w1 w2 begin - FFTL.run!(p(y1), p(x), p(w1), p(w2), c, 1); FFTL.run!(p(y2), p(x), p(w1), p(w2), rc, 1) - FFTL.run!(p(z1), p(z1), p(w1), p(w2), c, 1); FFTL.run!(p(z2), p(z2), p(w1), p(w2), rc, 1) + y = similar(x); z = copy(x) + GC.@preserve x y z w1 w2 begin + FFTL.run!(p(y), p(x), p(w1), p(w2), rc, 1) + FFTL.run!(p(z), p(z), p(w1), p(w2), rc, 1) end - @test y1 == y2 && z1 == z2 + ks, ref = dft_at(x, -1) + @test relerr(at(y, ks), ref) < 4tol(T) + @test y == z # three lanes, element j of lane i at i + 3j; strided batch S = 3 through run! B = 3; xs = randn(Complex{T}, B * n); wl1 = similar(xs); wl2 = similar(xs) - l1 = similar(xs); l2 = similar(xs); s1 = similar(xs); s2 = similar(xs) - cS = FFTL.build_chain(T, n, 1, FFTL.radices(n, B; esize = sizeof(Complex{T}))) + l = similar(xs); s = similar(xs) rS = FFTL.RChain{T}(n, 1, FFTL.radices(n, B; esize = sizeof(Complex{T}))) - GC.@preserve xs l1 l2 s1 s2 wl1 wl2 begin - FFTL.run_lanes!(p(l1), p(xs), p(wl1), p(wl2), cS, B, B, B); FFTL.run_lanes!(p(l2), p(xs), p(wl1), p(wl2), rS, B, B, B) - FFTL.run!(p(s1), p(xs), p(wl1), p(wl2), cS, B); FFTL.run!(p(s2), p(xs), p(wl1), p(wl2), rS, B) + GC.@preserve xs l s wl1 wl2 begin + FFTL.run_lanes!(p(l), p(xs), p(wl1), p(wl2), rS, B, B, B) + FFTL.run!(p(s), p(xs), p(wl1), p(wl2), rS, B) end - @test l1 == l2 && s1 == s2 + for i in 0:B-1 + ks, ref = dft_at(xs[i+1:B:end], 1) + @test relerr(at(l, ks, B, i), ref) < 4tol(T) + @test relerr(at(s, ks, B, i), ref) < 4tol(T) + end + end + # a chain whose stage records are not linked to it (a copy, or one loaded from a package + # image) links them again at the first run + rs = FFTL.radices(60, 1; esize = sizeof(Complex{T})) + ref = FFTL.RChain{T}(60, -1, rs); rc = FFTL.build_rchain(T, 60, -1, rs) + unlink!(c) = (for i in eachindex(c.single) + st = c.single[i] + c.single[i] = FFTL.StageCall{T}(st.op, st.c, st.c1, st.m, st.s, st.rem, C_NULL) + c.batch[i] = FFTL.BatchStage{T}(c.batch[i].c, C_NULL) + end; c.linked = 0) + x = randn(Complex{T}, 4 * 60); y1 = similar(x); y2 = similar(x); w1 = similar(x); w2 = similar(x) + GC.@preserve x y1 y2 w1 w2 begin + unlink!(rc); FFTL.run!(p(y1), p(x), p(w1), p(w2), ref, 1); FFTL.run!(p(y2), p(x), p(w1), p(w2), rc, 1) + @test y1[1:60] == y2[1:60] + unlink!(rc); FFTL.run!(p(y1), p(x), p(w1), p(w2), ref, 4); FFTL.run!(p(y2), p(x), p(w1), p(w2), rc, 4) + @test y1 == y2 + unlink!(rc); FFTL.run_lanes!(p(y1), p(x), p(w1), p(w2), ref, 4, 4, 4); FFTL.run_lanes!(p(y2), p(x), p(w1), p(w2), rc, 4, 4, 4) + @test y1 == y2 + unlink!(rc); FFTL.run_rows!(p(y1), p(x), p(w1), p(w2), ref, 4, 60); FFTL.run_rows!(p(y2), p(x), p(w1), p(w2), rc, 4, 60) + @test y1 == y2 + end + # a chain of length 1 copies + rc = FFTL.RChain{T}(1, -1, Int[]); x = randn(Complex{T}, 3); y = similar(x); w = similar(x) + GC.@preserve x y w begin + FFTL.run!(p(y), p(x), p(w), p(w), rc, 1); @test y[1] == x[1] + FFTL.run!(p(y), p(x), p(w), p(w), rc, 3); @test y == x + fill!(y, 0); FFTL.run_rows!(p(y), p(x), p(w), p(w), rc, 3, 1); @test y == x end end @testset "small-stride stages, $T" for T in (Float64, Float32) - # stage_small!: every power-of-two radix as the first stage, in both chain forms, against - # the direct DFT + # stage_small!: every power-of-two radix as the first stage, against the direct DFT p(a) = Ptr{T}(pointer(a)) C = FFTL.ladder_ncomplex(T)[1] for r in (2, 4, 8, 16), σ in (-1, 1) rs = [r, 2C] n = prod(rs) - c = FFTL.build_chain(T, n, σ, rs); rc = FFTL.RChain{T}(n, σ, rs) + rc = FFTL.RChain{T}(n, σ, rs) @test (rc.single[1].op == FFTL.OP_SMALL) == (C > 4) - x = randn(Complex{T}, n); w1 = similar(x); w2 = similar(x); y1 = similar(x); y2 = similar(x) - GC.@preserve x y1 y2 w1 w2 begin - FFTL.run!(p(y1), p(x), p(w1), p(w2), c, 1); FFTL.run_single!(p(y2), p(x), p(w1), p(w2), rc) - end - @test y1 == y2 - @test relerr(y2, dft(ComplexF64.(x), σ)) < tol(T) + x = randn(Complex{T}, n); w1 = similar(x); w2 = similar(x); y = similar(x) + GC.@preserve x y w1 w2 FFTL.run_single!(p(y), p(x), p(w1), p(w2), rc) + @test relerr(y, dft(ComplexF64.(x), σ)) < tol(T) end end