From 36043bf5c885f1aba2120b6cc886ba683775d9d5 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Sat, 7 Nov 2026 15:48:17 -0700 Subject: [PATCH 1/8] Full-width small-stride stages for single lines MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Stages of a single line with stride s below the complex count C of the widest vector (s = 1, 2, 4 for Float32 on AVX-512) ran on narrower vectors. stage_small! vectorises over f = q + s·p at full width: contiguous loads, per-lane twiddles (tw1 for s = 1, a new table twg otherwise) and an in-register interleave of the R outputs into contiguous stores, its permutation checked by the generator. Both chain forms use it. The kernels (radices 2, 4, 8, 16 × sign × stride, the widest vector type) take the stride as a constant and the length as run-time data, so they are a fixed set: plain methods in kernel_methods.jl, reached from RChain through cstage_small! and the new stage call OP_SMALL of the single-line schedule. Intel Core i7-1185G7, one line, Float32, RChain, radices of 16s and 8s (µs, without -> with; the same code against itself differs by up to 3 %): 128 0.093 -> 0.069 256 0.176 -> 0.131 1024 0.81 -> 0.61 2048 1.60 -> 1.35 4096 3.75 -> 3.29 8192 7.85 -> 6.69 16384 18.0 -> 15.4 Float64 has no such stage on this machine. Not done: m and s as constants of a single line. It needs one kernel per length (compiled Chain only), and in Chain it measured 0.90-1.04 of the time with run-time m and s (lengths 16-16384, both element types), inside the spread. Recorded in docs/src/explanation/design-decisions.md. Co-Authored-By: Claude Opus 5.5 Co-Authored-By: Claude Fable 5.2 --- docs/src/explanation/design-decisions.md | 21 +++++++ src/kernel_methods.jl | 12 ++++ src/kernels.jl | 76 ++++++++++++++++++++++++ src/plan.jl | 39 +++++++++++- src/rchain.jl | 36 +++++++++-- test/runtests.jl | 20 +++++++ 6 files changed, 196 insertions(+), 8 deletions(-) diff --git a/docs/src/explanation/design-decisions.md b/docs/src/explanation/design-decisions.md index ca1600e2..1f581f55 100644 --- a/docs/src/explanation/design-decisions.md +++ b/docs/src/explanation/design-decisions.md @@ -95,6 +95,27 @@ bodies the `@generated` methods build. The `@generated` methods remain as the fa **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"). +## Small strides of a single line run at full vector width + +**Decision (2026-11-11).** A stage of a single line whose stride s is below the complex count C of the +widest vector (s divides C, power-of-two radix, not the last stage) runs `stage_small!`: it is vectorised +over q + s·p, so every input is one contiguous load, the twiddles are stored per lane, and the R outputs +are interleaved in registers into contiguous stores. The generator checks the permutation of the +interleave. Both chain forms use it; the kernels (radices 2, 4, 8, 16 × sign × stride, widest vector +only) are part of the fixed set of plain methods, with the stride as a constant of the kernel and the +length as run-time data. + +**Why.** Such stages ran on narrower vectors (the grouped first stage on 4 complex values, later stages +on s). `Float32` on AVX-512 (C = 8) is the case that gains; `Float64` has no such stage there. Intel +Core i7-1185G7, one line, `RChain`, radix sequences of 16s and 8s, µs without → with: +128: 0.093 → 0.069, 256: 0.176 → 0.131, 1024: 0.81 → 0.61, 2048: 1.60 → 1.35, 4096: 3.75 → 3.29, +8192: 7.85 → 6.69, 16384: 18.0 → 15.4 (the same code against itself differs by up to 3 %). + +**Not done.** Passing m and s to the stage kernels as constants of a single line (immediates for the +stream offsets) needs one kernel per length, which only the compiled `Chain` can have. Measured in +`Chain` against the same kernels with run-time m and s, lengths 16–16384, both element types: ratios +0.90–1.04 without a trend, inside the spread of the measurement. Not verified on other machines. + ## Codelets of the radices 17–31 are compiled on first use **Decision (2026-11-10).** The codelets of the radices 17, 19, 23, 29 and 31 are not part of the package diff --git a/src/kernel_methods.jl b/src/kernel_methods.jl index e90b555d..85a2ff7d 100644 --- a/src/kernel_methods.jl +++ b/src/kernel_methods.jl @@ -40,6 +40,18 @@ for T in (Float32, Float64), (name, head, vts) in KERNEL_HEADS end end +# stage_small!: the widest vector type, power-of-two radices, every stride s < C that divides C. +for T in (Float32, Float64), σ in (-1, 1), R in RADICES + VT = first(vec_ladder(T)); C = ncomplex(VT) + (ispow2(R) && !(R in LAZY_RADICES)) || continue + for s in 1:C-1 + C % s == 0 || continue + h = :(stage_small!(py::Ptr{$T}, px::Ptr{$T}, pt::Ptr{$T}, ms::Int, ::Val{$R}, ::Val{$σ}, ::Type{$VT}, ::Val{$s})) + b = kernel_subst(stage_small_body(T, R, σ, VT, s), Dict{Symbol,Any}(:T => T, :VT => VT)) + @eval $(Expr(:function, h, b)) + end +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 diff --git a/src/kernels.jl b/src/kernels.jl index 8120237b..1529a44d 100644 --- a/src/kernels.jl +++ b/src/kernels.jl @@ -194,6 +194,82 @@ end stage_first_group_range_body(T, R, σ, VT) end +# Complex values of u ++ v (2C) in chunks of S, alternating u_0 v_0 u_1 v_1 …: the lane masks of +# the first C and the last C. +function chunk_masks(C, S) + idx = [(isodd(i ÷ S) ? C : 0) + (i ÷ S ÷ 2) * S + i % S for i in 0:2C-1] + lanes(ix) = Tuple(l for c in ix for l in (2c, 2c + 1)) + lanes(idx[1:C]), lanes(idx[C+1:2C]) +end + +# Interleave R = 2^r vectors in chunks of S complex values: pairs at S, then the lower and upper +# halves recursively at 2S; chunks of a whole vector are already in place. labs tracks +# (vector, complex value) of every complex value so the generator can check the permutation. +function ilv_net!(ex, ys, labs, S, C) + (length(ys) == 1 || S >= C) && return ys, labs + ml, mh = chunk_masks(C, S) + lo = Any[]; hi = Any[]; llo = Any[]; lhi = Any[] + for i in 1:2:length(ys) + cat = vcat(labs[i], labs[i+1]) + push!(llo, cat[collect(ml[1:2:end]) .÷ 2 .+ 1]); push!(lhi, cat[collect(mh[1:2:end]) .÷ 2 .+ 1]) + a = gensym(:lo); b = gensym(:hi) + push!(ex, :($a = shufflevector($(ys[i]), $(ys[i+1]), Val($ml))), + :($b = shufflevector($(ys[i]), $(ys[i+1]), Val($mh)))) + push!(lo, a); push!(hi, b) + end + o1, l1 = ilv_net!(ex, lo, llo, 2S, C) + o2, l2 = ilv_net!(ex, hi, lhi, 2S, C) + vcat(o1, o2), vcat(l1, l2) +end + +# Outputs X_k (lanes p·s + q) in the order of y: vector t, lane l holds position tC + l. +function chunk_interleave!(ex, ys, s, C) + R = length(ys) + out, labs = ilv_net!(ex, ys, [[(k, l) for l in 0:C-1] for k in 0:R-1], s, C) + for t in 0:R-1, l in 0:C-1 + pos = t * C + l + labs[t+1][l+1] == ((pos % (s * R)) ÷ s, (pos ÷ (s * R)) * s + pos % s) || + error("chunk interleave of $R vectors at stride $s, $C complex values is wrong") + end + out +end + +# Stage of a single line with stride s < C (C complex values per vector, s | C, R a power of +# two): vectorised over f = q + s·p, so each input is one contiguous load and the R outputs are +# interleaved in registers into contiguous stores. Twiddles per lane: w^(kp) at (k-1)ms + f, +# p = f ÷ s (the grouped first-stage table for s = 1). ms = m·s must be a multiple of C; s is a +# constant because the interleave depends on it. +function stage_small_body(::Type{T}, R::Int, σ::Int, ::Type{VT}, s::Int) where {T,VT} + stmts, a, X = codelet(R, σ) + C = ncomplex(VT); B = 2sizeof(T) + loads = [:($(a[j]) = vload($VT, px + (f + $(j - 1) * ms) * $B)) for j in 1:R] + outs = Any[X[1]] + twm = Any[] + for k in 1:R-1 + b = Symbol("b", k) + push!(twm, :($b = ccmul($(X[k+1]), vload($VT, pt + (f + $(k - 1) * ms) * $B)))) + push!(outs, b) + end + shuf = Any[] + ys = chunk_interleave!(shuf, outs, s, C) + stores = [:(vstore($(ys[t+1]), py + ($R * f + $(t * C)) * $B)) for t in 0:R-1] + quote + for f in 0:$C:ms-1 + $(loads...) + $(stmts...) + $(twm...) + $(shuf...) + $(stores...) + end + nothing + end +end + +@generated function stage_small!(py::Ptr{T}, px::Ptr{T}, pt::Ptr{T}, ms::Int, ::Val{R}, ::Val{σ}, + ::Type{VT}, ::Val{s}) where {T,R,σ,VT,s} + stage_small_body(T, R, σ, VT, s) +end + # Scalar tail of the grouped first stage: one complex value per vector, the twiddle table in # the grouped layout (k-1)m + p, one store per output. function stage_first_group_1_body(::Type{T}, R::Int, σ::Int, ::Type{VT}) where {T,VT} diff --git a/src/plan.jl b/src/plan.jl index 7c65e909..cb9401a2 100644 --- a/src/plan.jl +++ b/src/plan.jl @@ -14,6 +14,7 @@ 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 + twg::NTuple{K,Vector{Complex{T}}} # stage_small! with s > 1: w^(kp) at (k-1)ms + f, p = f ÷ s; else empty end # Read-only parts of kernels (chains, four-step twiddles, Bluestein chirp and filter) are shared @@ -77,9 +78,36 @@ Chain{T}(n::Int, σ::Int, rs::AbstractVector{Int} = radices(n)) where {T} = 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) + twg = small_tables(ncomplex(first(vec_ladder(T))), n, rs, tw) + Chain{T,Tuple(rs),σ,length(rs)}(n, Tuple(tw), tw1, Tuple(twg)) +end + +# Twiddles of the stage_small! stages with s > 1 (empty for every other stage), from the stage +# tables tw. C: complex values of the widest vector. +function small_tables(C::Int, n::Int, rs::AbstractVector{Int}, tw::Vector{Vector{Complex{T}}}) where {T} + twg = Vector{Complex{T}}[] + s = 1 + for (i, r) in enumerate(rs) + push!(twg, s > 1 && i < length(rs) && small_stride(C, r, n, s) ? small_twiddles(tw[i], r, s, n ÷ r) : + Complex{T}[]) + s *= r + end + twg end +# The stage table t (w^(kp) at (r-1)p + k) per lane: w^(kp) at (k-1)ms + f, p = f ÷ s. +function small_twiddles(t::Vector{Complex{T}}, r::Int, s::Int, ms::Int) where {T} + g = Vector{Complex{T}}(undef, (r - 1) * ms) + for k in 1:r-1, f in 0:ms-1 + g[(k-1)*ms+f+1] = t[(r-1)*(f÷s)+k] + end + g +end + +# Stages of a single line with stride s below the complex count C of the widest vector run +# stage_small! (full vectors, outputs interleaved in registers). +small_stride(C::Int, r::Int, n::Int, s::Int) = ispow2(r) && s < C && C % s == 0 && (n ÷ r) % C == 0 + # Stage twiddles and first-stage table of a radix sequence, shared by Chain and RChain. 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}}} @@ -220,7 +248,14 @@ function chain_body(T, Rs, σ, src, dsts, base, single) # 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 i == K + # small strides at full width, where the grouped first stage would not be + small = single && i < K && small_stride(ncomplex(ladder[1]), r, n, sst) && + (sst > 1 || group === nothing || ncomplex(ladder[group]) < ncomplex(ladder[1])) + if small + pt = sst == 1 ? :(cptr(P.tw1)) : :(cptr(P.twg[$i])) + push!(body, :(stage_small!($dst, $src, $pt, $(n ÷ r), Val($r), Val($σ), $(ladder[1]), Val($sst))), + :(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 diff --git a/src/rchain.jl b/src/rchain.jl index 5e8b5b8f..33e52ea0 100644 --- a/src/rchain.jl +++ b/src/rchain.jl @@ -46,6 +46,18 @@ all(in(RADICES), LAZY_RADICES) || error("LAZY_RADICES must be radices from RADIC # 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])) +# The stage_small! call of a codelet: a switch over the strides s < C that divide C (a constant +# of the kernel). A callable of a lazy radix (V a type parameter) takes s as a run-time Val. +function small_call(r, σ, V) + V isa Type || return :(stage_small!(py, px, pt, ms, Val($r), Val($σ), $V, Val(s))) + C = ncomplex(V) + ex = :(throw(ArgumentError("no small-stride stage for this stride"))) + for s in C-1:-1:1 + C % s == 0 && (ex = :(s == $s ? stage_small!(py, px, pt, ms, Val($r), Val($σ), $V, Val($s)) : $ex)) + end + ex +end + # 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 @@ -60,6 +72,7 @@ for T in (Float32, Float64) (:cstage_last!, (:py, :px, :s), (r, σ, V) -> :(stage_last!(py, px, s, Val($r), Val($σ), $V))), (:cstage_last_range!, (:py, :px, :s, :qlo, :qhi), (r, σ, V) -> :(stage_last!(py, px, s, qlo, qhi, Val($r), Val($σ), $V))), (:cstage_group!, (:py, :px, :pt, :m), (r, σ, V) -> :(stage_first_group!(py, px, pt, m, Val($r), Val($σ), $V))), + (:cstage_small!, (:py, :px, :pt, :ms, :s), (r, σ, V) -> small_call(r, σ, V)), (:cstage_lanes!, (:py, :px, :pt, :m, :s0, :sx, :sy, :B), (r, σ, V) -> :(stage_lanes!(py, px, pt, m, s0, sx, sy, B, Val($r), Val($σ), $V))), (:cstage_last_lanes!, (:py, :px, :s0, :sx, :sy, :B), (r, σ, V) -> :(stage_last_lanes!(py, px, s0, sx, sy, B, Val($r), Val($σ), $V)))) for (name, args, call) in kinds @@ -70,6 +83,8 @@ for T in (Float32, Float64) # never selects another vector type for it, and code compiled ahead of time must not # contain the call (its generator would fail) name === :cstage_group! && !(ncomplex(V) in (2, 4)) && continue + # stage_small! exists for the widest vector type and power-of-two radices + name === :cstage_small! && !(vi == 1 && ispow2(r)) && continue r in LAZY_RADICES && continue c = ((si - 1) * NL + vi - 1) * NR + ri push!(cases, :(c == $c && return $(call(r, σ, V)))) @@ -115,6 +130,7 @@ const OP_MID_TAIL = 2 # q in vectors up to s - rem, the scalar remainder aft const OP_LAST = 3 const OP_LAST_TAIL = 4 const OP_GROUP = 5 # first stage vectorised over the butterflies p +const OP_SMALL = 6 # stride below the widest vector: vectorised over q + s·p (stage_small!) struct StageCall op::Int @@ -136,6 +152,7 @@ struct RChain{T} <: Kernel{T} 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 + twg::Vector{Vector{Complex{T}}} # stage_small! with s > 1: w^(kp) at (k-1)ms + f, p = f ÷ s; else empty single::Vector{StageCall} # the stage calls of a single line end @@ -148,13 +165,15 @@ RChain{T}(n::Int, σ::Int, rs::AbstractVector{Int} = radices(n)) where {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))) + twg = small_tables(ladder_ncomplex(T)[1], n, rs, tw) + RChain{T}(n, σ, collect(Int, rs), tw, tw1, twg, single_schedule(T, n, σ, collect(Int, rs))) end # The stage calls of a single line, with the rules of `chain_body` for S = 1: 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. +# when its stride is 1; power-of-two stages with a stride below the widest vector run +# stage_small! at full width. function single_schedule(::Type{T}, n::Int, σ::Int, rs::Vector{Int}) where {T} ncs = ladder_ncomplex(T) nv = length(ncs) @@ -177,7 +196,10 @@ function single_schedule(::Type{T}, n::Int, σ::Int, rs::Vector{Int}) where {T} rem = s % ncs[vi] c = codelet_index(T, σ, vi, r); c1 = codelet_index(T, σ, nv, r) m = L ÷ r - if i == K + if i < K && small_stride(ncs[1], r, n, s) && (s > 1 || group == 0 || ncs[group] < ncs[1]) + # small strides at full width, where the grouped first stage would not be + push!(calls, StageCall(OP_SMALL, codelet_index(T, σ, 1, r), c1, m, s, 0)) + elseif i == K push!(calls, StageCall(rem == 0 ? OP_LAST : OP_LAST_TAIL, c, c1, m, s, rem)) elseif i == 1 && group != 0 && ncs[vi] == 1 push!(calls, StageCall(OP_GROUP, codelet_index(T, σ, group, r), c1, m, s, 0)) @@ -276,6 +298,8 @@ function run_single!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::RChain 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(st.s == 1 ? P.tw1 : P.twg[i]), st.m * st.s, st.s, st.c) elseif op == OP_LAST cstage_last!(dst, src, st.s, st.c) elseif op == OP_MID_TAIL @@ -349,13 +373,13 @@ 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) + c = C(P.n, ntuple(i -> P.tw[i], Val(K)), P.tw1, ntuple(i -> P.twg[i], Val(K))) 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) + c = C(P.n, ntuple(i -> P.tw[i], Val(K)), P.tw1, ntuple(i -> P.twg[i], Val(K))) GC.@preserve P run_lanes!(py, px, pw1, pw2, c, Sx, Sy, B) end @@ -404,7 +428,7 @@ end # 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) + c = C(P.n, ntuple(i -> P.tw[i], Val(K)), P.tw1, ntuple(i -> P.twg[i], Val(K))) 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) diff --git a/test/runtests.jl b/test/runtests.jl index a8d1a458..66e75ba8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -753,6 +753,26 @@ thin(T, v) = FULL || T == Float64 ? v : v[1:2:end] end end + @testset "small-stride stages, $T" for T in (Float64, Float32) + # stage_small!: every power-of-two radix at every stride below the widest vector, in both + # chain forms, against the direct DFT + p(a) = Ptr{T}(pointer(a)) + C = FFTL.ladder_ncomplex(T)[1] + for s in 1:C-1, r in (2, 4, 8, 16), σ in (-1, 1) + C % s == 0 || continue + rs = s == 1 ? [r, 2C] : [s, r, 2C] + n = prod(rs) + c = FFTL.build_chain(T, n, σ, rs); rc = FFTL.RChain{T}(n, σ, rs) + s == 1 || @test rc.single[2].op == FFTL.OP_SMALL + 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) + end + end + # run-time plans (RChain kernels): same results as the reference, four-step and plane fusion @testset "RPlan, $T" for T in pick((Float64, Float32), (Float64,)) for (sz, reg, nt) in pick([((1000,), [1], 1), ((2 * 3^11,), [1], 3), ((1009,), [1], 1), ((64, 48), [1, 2], 1), From 6dec53ad32424acf0c606f2a98062489d040a880 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Sat, 7 Nov 2026 15:48:17 -0700 Subject: [PATCH 2/8] Radix 16 for single lines on AVX-512 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Powers of two from 128 to 16384 take as many 16s as possible and 8s for the rest (512 stays 8 8 8). Same rule for Float32 and Float64. Intel Core i7-1185G7, one line, RChain (µs, 4s and 8s -> 16s and 8s; the same code against itself differs by up to 3 %): F64 128 0.145 -> 0.128 256 0.302 -> 0.256 1024 1.51 -> 1.34 2048 3.99 -> 3.53 4096 7.64 -> 7.13 8192 18.0 -> 16.3 16384 36.9 -> 36.2 F32 128 0.074 -> 0.069 256 0.131 -> 0.131 1024 0.632 -> 0.611 2048 1.27 -> 1.35 4096 3.25 -> 3.29 8192 7.26 -> 6.69 16384 16.0 -> 15.4 Float64 is 2-15 % faster, Float32 between 6 % slower (2048) and 8 % faster. 512 as 16 4 8 was not measured on the run-time chains; it stays 8 8 8. Co-Authored-By: Claude Opus 5.5 Co-Authored-By: Claude Fable 5.2 --- src/factor.jl | 8 ++++++++ test/data/radices-x86_64-32regs.csv | 14 +++++++------- 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/src/factor.jl b/src/factor.jl index 8d95ffe8..334cef30 100644 --- a/src/factor.jl +++ b/src/factor.jl @@ -7,6 +7,9 @@ # up to 4096 on the M1; on Tiger Lake (AVX-512) 7–30 % from 4096 on (S = 1) and up to 45 % in # batches, within noise below; on Zen 3 (AVX2, 16 # registers) it is 3–25 % slower from 256 on. +# Single lines on AVX-512 take radix 16 as often as possible and 8s for the rest (Tiger Lake, +# 128–16384: Float64 2–15 % faster, Float32 between 6 % slower and 8 % faster). Not measured +# on NEON. # Otherwise radix 4 throughout, one 8 (or 2) for an odd exponent; radix 16 and above lose # (register spills, cache-set conflicts at large strides). # - aarch64: a factor 3 is merged with a 4 into 12, else with a 2 into 6; a factor 5 with a spare @@ -105,6 +108,11 @@ function heuristic_radices(n::Int, S::Int) end if RADIX8_POW2 && ispow2(n) && 64 < n <= 16384 && S <= LARGE_STRIDE k = trailing_zeros(n) + if HAS_AVX512 && S == 1 + # as many 16s as possible, the rest in 8s: 4t ≤ k, t ≡ k (mod 3) + t = k % 3 + 3 * ((k ÷ 4 - k % 3) ÷ 3) + return vcat(fill(16, t), fill(8, (k - 4t) ÷ 3)) + end spare = k % 3 return vcat(spare == 2 ? [4] : spare == 1 ? [4, 4] : Int[], fill(8, (k - (spare == 2 ? 2 : spare == 1 ? 4 : 0)) ÷ 3)) end diff --git a/test/data/radices-x86_64-32regs.csv b/test/data/radices-x86_64-32regs.csv index 890ce426..291c0d04 100644 --- a/test/data/radices-x86_64-32regs.csv +++ b/test/data/radices-x86_64-32regs.csv @@ -192,7 +192,7 @@ 125,1024,5 5 5 126,1,2 7 9 126,1024,7 6 3 -128,1,4 4 8 +128,1,16 8 128,1024,8 4 4 130,1,2 5 13 130,1024,13 10 @@ -330,7 +330,7 @@ 253,1024,23 11 255,1,3 5 17 255,1024,17 5 3 -256,1,4 8 8 +256,1,16 16 256,1024,4 4 4 4 260,1,4 5 13 260,1024,4 5 13 @@ -878,15 +878,15 @@ 1020,1024,17 12 5 1023,1,3 11 31 1023,1024,31 11 3 -1024,1,4 4 8 8 +1024,1,16 8 8 1024,1024,4 4 4 4 4 -2048,1,4 8 8 8 +2048,1,16 16 8 2048,1024,8 4 4 4 4 -4096,1,8 8 8 8 +4096,1,16 16 16 4096,1024,4 4 4 4 4 4 -8192,1,4 4 8 8 8 +8192,1,16 8 8 8 8192,1024,8 4 4 4 4 4 -16384,1,4 8 8 8 8 +16384,1,16 16 8 8 16384,1024,4 4 4 4 4 4 4 32768,1,8 4 4 4 4 4 4 32768,1024,8 4 4 4 4 4 4 From fef0ff9387f8dd8b7780133f2a1f5e809c641f23 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Sat, 7 Nov 2026 16:18:09 -0700 Subject: [PATCH 3/8] Don't compile the run-time codelet table for every first plan kernel() compiles its Bluestein branch, whose filter transform calls run! with an uninferred Chain. Inference enumerated the few run! methods and compiled run!(::RChain), whose dispatch over every radix, sign and vector type instantiates all codelets. max_methods 1 on run! leaves such calls to run-time dispatch; concretely typed calls (the run-time plans) are unaffected. First fft of a fresh session, n = 1024: F64 19.6 -> 13.9 s (607 -> 19 kernel specializations), F32 afterwards 22.7 -> 5.0 s. The test suite takes as long (it compiles the run-time table in its own testsets). Co-Authored-By: Claude Opus 5.5 --- src/plan.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/plan.jl b/src/plan.jl index cb9401a2..24966bd5 100644 --- a/src/plan.jl +++ b/src/plan.jl @@ -286,6 +286,11 @@ function stride_dispatch(T, body) :(S == 1 ? $(body(1, true)) : $ex) end +# A call with an uninferred kernel (Bluestein's Chain{T}(M, σ) comes out of the cache as Any) must +# dispatch at run time: enumerating the few run! methods compiled run!(::RChain) and with it every +# codelet of the run-time dispatch (~20 s on the first plan of a session). +Base.Experimental.@max_methods 1 function run! end + """ run!(py, px, pw1, pw2, P::Kernel, S) From 761b8feb57aadfd6c6c49983ab4fe3ca4f2fa34f Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Wed, 11 Nov 2026 20:16:45 -0700 Subject: [PATCH 4/8] drop max_method --- src/plan.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/plan.jl b/src/plan.jl index 24966bd5..cb9401a2 100644 --- a/src/plan.jl +++ b/src/plan.jl @@ -286,11 +286,6 @@ function stride_dispatch(T, body) :(S == 1 ? $(body(1, true)) : $ex) end -# A call with an uninferred kernel (Bluestein's Chain{T}(M, σ) comes out of the cache as Any) must -# dispatch at run time: enumerating the few run! methods compiled run!(::RChain) and with it every -# codelet of the run-time dispatch (~20 s on the first plan of a session). -Base.Experimental.@max_methods 1 function run! end - """ run!(py, px, pw1, pw2, P::Kernel, S) From b60c8b9d237c6829df27ff62edabbad14f1bc70e Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 12 Nov 2026 05:06:00 +0000 Subject: [PATCH 5/8] Full-width stage_small! for the first stage only stage_small! runs the first stage of a single line (stride 1, power-of-two radix) on the widest vector where the grouped first stage is narrower (Float32 on AVX-512). Its twiddles are the grouped first-stage table tw1, so the chains carry no per-lane tables, the kernels have no stride parameter and the RChain dispatch has no stride switch. Stages with a stride of 2 or 4 below the widest vector run on the vector that fits their stride. Only Float32 lengths whose radices start with 4, 4 have such a stage (stride 4; 625 lengths up to 70000, 65536 the only power of two), none with stride 2, none for Float64. The lengths with radices of 16s and 8s (128-16384) never had one. AMD EPYC 9354, one line, Float32, RChain, time without the full-width stride-4 stage over time with it, medians of 10 interleaved rounds (65536 and 262144: 5 rounds); the same code against itself differs by 5-20 % between rounds: 192 1.08 320 1.06 768 1.03 960 1.03 3072 0.96 4800 0.99 20480 0.97 65536 0.96 262144 0.97 Only the 6-8 % at 192 and 320 repeated in a second series. Not measured on other machines. Recorded in docs/src/explanation/design-decisions.md. Co-Authored-By: Claude Fable 5.2 --- docs/src/explanation/design-decisions.md | 35 ++++++++++++-------- src/kernel_methods.jl | 13 +++----- src/kernels.jl | 37 ++++++++++----------- src/plan.jl | 42 +++++------------------- src/rchain.jl | 34 +++++-------------- test/runtests.jl | 11 +++---- 6 files changed, 66 insertions(+), 106 deletions(-) diff --git a/docs/src/explanation/design-decisions.md b/docs/src/explanation/design-decisions.md index 1f581f55..be313d78 100644 --- a/docs/src/explanation/design-decisions.md +++ b/docs/src/explanation/design-decisions.md @@ -95,22 +95,31 @@ bodies the `@generated` methods build. The `@generated` methods remain as the fa **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"). -## Small strides of a single line run at full vector width - -**Decision (2026-11-11).** A stage of a single line whose stride s is below the complex count C of the -widest vector (s divides C, power-of-two radix, not the last stage) runs `stage_small!`: it is vectorised -over q + s·p, so every input is one contiguous load, the twiddles are stored per lane, and the R outputs -are interleaved in registers into contiguous stores. The generator checks the permutation of the -interleave. Both chain forms use it; the kernels (radices 2, 4, 8, 16 × sign × stride, widest vector -only) are part of the fixed set of plain methods, with the stride as a constant of the kernel and the -length as run-time data. - -**Why.** Such stages ran on narrower vectors (the grouped first stage on 4 complex values, later stages -on s). `Float32` on AVX-512 (C = 8) is the case that gains; `Float64` has no such stage there. Intel -Core i7-1185G7, one line, `RChain`, radix sequences of 16s and 8s, µs without → with: +## The first stage of a single line runs at full vector width + +**Decision (2026-11-12).** The first stage of a single line (stride 1, power-of-two radix, not the last +stage) runs `stage_small!` where the widest vector holds more complex values than the grouped first +stage handles (C = 8, `Float32` on AVX-512): it is vectorised over the butterflies p, so every input is +one contiguous load, the twiddles come from the grouped first-stage table, and the R outputs are +interleaved in registers into contiguous stores. The generator checks the permutation of the +interleave. Both chain forms use it; the kernels (radices 2, 4, 8, 16 × sign, widest vector only) are +part of the fixed set of plain methods. + +**Why.** That stage ran the grouped first stage on 4 complex values. `Float64` has no such stage on +AVX-512. Intel Core i7-1185G7, one line, `Float32`, `RChain`, radix sequences of 16s and 8s, µs without +→ with: 128: 0.093 → 0.069, 256: 0.176 → 0.131, 1024: 0.81 → 0.61, 2048: 1.60 → 1.35, 4096: 3.75 → 3.29, 8192: 7.85 → 6.69, 16384: 18.0 → 15.4 (the same code against itself differs by up to 3 %). +**Not done.** The same scheme for later stages with a stride s = 2 or 4 below C (vectorised over +q + s·p). It needs a per-lane twiddle table for every such stage, one kernel per stride and a stride +switch in the `RChain` dispatch. Only `Float32` lengths whose radices start with 4, 4 have such a stage +(s = 4; 625 lengths up to 70000, 65536 the only power of two), none with s = 2. AMD EPYC 9354, one line, +`Float32`, `RChain`, time without the s = 4 stage over time with it, medians of 10 interleaved rounds: +192: 1.08, 320: 1.06, 768: 1.03, 960: 1.03, 3072: 0.96, 4800: 0.99, 20480: 0.97; 65536: 0.96 and +262144: 0.97 from 5 rounds. The same code against itself differs by 5–20 % between rounds on that +machine, so only the 6–8 % at 192 and 320 (the same in a second series) stands out. + **Not done.** Passing m and s to the stage kernels as constants of a single line (immediates for the stream offsets) needs one kernel per length, which only the compiled `Chain` can have. Measured in `Chain` against the same kernels with run-time m and s, lengths 16–16384, both element types: ratios diff --git a/src/kernel_methods.jl b/src/kernel_methods.jl index 85a2ff7d..5aa839c0 100644 --- a/src/kernel_methods.jl +++ b/src/kernel_methods.jl @@ -40,16 +40,13 @@ for T in (Float32, Float64), (name, head, vts) in KERNEL_HEADS end end -# stage_small!: the widest vector type, power-of-two radices, every stride s < C that divides C. +# stage_small!: the widest vector type, power-of-two radices. for T in (Float32, Float64), σ in (-1, 1), R in RADICES - VT = first(vec_ladder(T)); C = ncomplex(VT) + VT = first(vec_ladder(T)) (ispow2(R) && !(R in LAZY_RADICES)) || continue - for s in 1:C-1 - C % s == 0 || continue - h = :(stage_small!(py::Ptr{$T}, px::Ptr{$T}, pt::Ptr{$T}, ms::Int, ::Val{$R}, ::Val{$σ}, ::Type{$VT}, ::Val{$s})) - b = kernel_subst(stage_small_body(T, R, σ, VT, s), Dict{Symbol,Any}(:T => T, :VT => VT)) - @eval $(Expr(:function, h, b)) - end + h = :(stage_small!(py::Ptr{$T}, px::Ptr{$T}, pt::Ptr{$T}, m::Int, ::Val{$R}, ::Val{$σ}, ::Type{$VT})) + 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 diff --git a/src/kernels.jl b/src/kernels.jl index 1529a44d..ed295f13 100644 --- a/src/kernels.jl +++ b/src/kernels.jl @@ -222,39 +222,36 @@ function ilv_net!(ex, ys, labs, S, C) vcat(o1, o2), vcat(l1, l2) end -# Outputs X_k (lanes p·s + q) in the order of y: vector t, lane l holds position tC + l. -function chunk_interleave!(ex, ys, s, C) +# Outputs X_k (lane p) in the order of y: vector t, lane l holds position tC + l = Rp + k. +function chunk_interleave!(ex, ys, C) R = length(ys) - out, labs = ilv_net!(ex, ys, [[(k, l) for l in 0:C-1] for k in 0:R-1], s, C) + out, labs = ilv_net!(ex, ys, [[(k, l) for l in 0:C-1] for k in 0:R-1], 1, C) for t in 0:R-1, l in 0:C-1 - pos = t * C + l - labs[t+1][l+1] == ((pos % (s * R)) ÷ s, (pos ÷ (s * R)) * s + pos % s) || - error("chunk interleave of $R vectors at stride $s, $C complex values is wrong") + labs[t+1][l+1] == ((t * C + l) % R, (t * C + l) ÷ R) || + error("interleave of $R vectors of $C complex values is wrong") end out end -# Stage of a single line with stride s < C (C complex values per vector, s | C, R a power of -# two): vectorised over f = q + s·p, so each input is one contiguous load and the R outputs are -# interleaved in registers into contiguous stores. Twiddles per lane: w^(kp) at (k-1)ms + f, -# p = f ÷ s (the grouped first-stage table for s = 1). ms = m·s must be a multiple of C; s is a -# constant because the interleave depends on it. -function stage_small_body(::Type{T}, R::Int, σ::Int, ::Type{VT}, s::Int) where {T,VT} +# First stage of a single line (s = 1) at any vector width, for R a power of two and m a +# multiple of C: butterflies p, …, p+C-1 share a register (twiddles in the grouped first-stage +# table) and the R outputs are interleaved in registers into contiguous stores. +function stage_small_body(::Type{T}, R::Int, σ::Int, ::Type{VT}) where {T,VT} stmts, a, X = codelet(R, σ) C = ncomplex(VT); B = 2sizeof(T) - loads = [:($(a[j]) = vload($VT, px + (f + $(j - 1) * ms) * $B)) for j in 1:R] + loads = [:($(a[j]) = vload($VT, px + (p + $(j - 1) * m) * $B)) for j in 1:R] outs = Any[X[1]] twm = Any[] for k in 1:R-1 b = Symbol("b", k) - push!(twm, :($b = ccmul($(X[k+1]), vload($VT, pt + (f + $(k - 1) * ms) * $B)))) + push!(twm, :($b = ccmul($(X[k+1]), vload($VT, pt + (p + $(k - 1) * m) * $B)))) push!(outs, b) end shuf = Any[] - ys = chunk_interleave!(shuf, outs, s, C) - stores = [:(vstore($(ys[t+1]), py + ($R * f + $(t * C)) * $B)) for t in 0:R-1] + ys = chunk_interleave!(shuf, outs, C) + stores = [:(vstore($(ys[t+1]), py + ($R * p + $(t * C)) * $B)) for t in 0:R-1] quote - for f in 0:$C:ms-1 + for p in 0:$C:m-1 $(loads...) $(stmts...) $(twm...) @@ -265,9 +262,9 @@ function stage_small_body(::Type{T}, R::Int, σ::Int, ::Type{VT}, s::Int) where end end -@generated function stage_small!(py::Ptr{T}, px::Ptr{T}, pt::Ptr{T}, ms::Int, ::Val{R}, ::Val{σ}, - ::Type{VT}, ::Val{s}) where {T,R,σ,VT,s} - stage_small_body(T, R, σ, VT, s) +@generated function stage_small!(py::Ptr{T}, px::Ptr{T}, pt::Ptr{T}, m::Int, ::Val{R}, ::Val{σ}, + ::Type{VT}) where {T,R,σ,VT} + stage_small_body(T, R, σ, VT) end # Scalar tail of the grouped first stage: one complex value per vector, the twiddle table in diff --git a/src/plan.jl b/src/plan.jl index cb9401a2..f5d1bcd2 100644 --- a/src/plan.jl +++ b/src/plan.jl @@ -14,7 +14,6 @@ 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 - twg::NTuple{K,Vector{Complex{T}}} # stage_small! with s > 1: w^(kp) at (k-1)ms + f, p = f ÷ s; else empty end # Read-only parts of kernels (chains, four-step twiddles, Bluestein chirp and filter) are shared @@ -78,35 +77,14 @@ Chain{T}(n::Int, σ::Int, rs::AbstractVector{Int} = radices(n)) where {T} = function build_chain(::Type{T}, n::Int, σ::Int, rs::AbstractVector{Int}) where {T} tw, tw1 = chain_tables(T, n, σ, rs) - twg = small_tables(ncomplex(first(vec_ladder(T))), n, rs, tw) - Chain{T,Tuple(rs),σ,length(rs)}(n, Tuple(tw), tw1, Tuple(twg)) -end - -# Twiddles of the stage_small! stages with s > 1 (empty for every other stage), from the stage -# tables tw. C: complex values of the widest vector. -function small_tables(C::Int, n::Int, rs::AbstractVector{Int}, tw::Vector{Vector{Complex{T}}}) where {T} - twg = Vector{Complex{T}}[] - s = 1 - for (i, r) in enumerate(rs) - push!(twg, s > 1 && i < length(rs) && small_stride(C, r, n, s) ? small_twiddles(tw[i], r, s, n ÷ r) : - Complex{T}[]) - s *= r - end - twg -end - -# The stage table t (w^(kp) at (r-1)p + k) per lane: w^(kp) at (k-1)ms + f, p = f ÷ s. -function small_twiddles(t::Vector{Complex{T}}, r::Int, s::Int, ms::Int) where {T} - g = Vector{Complex{T}}(undef, (r - 1) * ms) - for k in 1:r-1, f in 0:ms-1 - g[(k-1)*ms+f+1] = t[(r-1)*(f÷s)+k] - end - g + Chain{T,Tuple(rs),σ,length(rs)}(n, Tuple(tw), tw1) end -# Stages of a single line with stride s below the complex count C of the widest vector run -# stage_small! (full vectors, outputs interleaved in registers). -small_stride(C::Int, r::Int, n::Int, s::Int) = ispow2(r) && s < C && C % s == 0 && (n ÷ r) % C == 0 +# 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. chain_tables(::Type{T}, n::Int, σ::Int, rs::AbstractVector{Int}) where {T} = @@ -248,12 +226,8 @@ function chain_body(T, Rs, σ, src, dsts, base, single) # 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) - # small strides at full width, where the grouped first stage would not be - small = single && i < K && small_stride(ncomplex(ladder[1]), r, n, sst) && - (sst > 1 || group === nothing || ncomplex(ladder[group]) < ncomplex(ladder[1])) - if small - pt = sst == 1 ? :(cptr(P.tw1)) : :(cptr(P.twg[$i])) - push!(body, :(stage_small!($dst, $src, $pt, $(n ÷ r), Val($r), Val($σ), $(ladder[1]), Val($sst))), + 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)) diff --git a/src/rchain.jl b/src/rchain.jl index 33e52ea0..583bee35 100644 --- a/src/rchain.jl +++ b/src/rchain.jl @@ -46,18 +46,6 @@ all(in(RADICES), LAZY_RADICES) || error("LAZY_RADICES must be radices from RADIC # 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])) -# The stage_small! call of a codelet: a switch over the strides s < C that divide C (a constant -# of the kernel). A callable of a lazy radix (V a type parameter) takes s as a run-time Val. -function small_call(r, σ, V) - V isa Type || return :(stage_small!(py, px, pt, ms, Val($r), Val($σ), $V, Val(s))) - C = ncomplex(V) - ex = :(throw(ArgumentError("no small-stride stage for this stride"))) - for s in C-1:-1:1 - C % s == 0 && (ex = :(s == $s ? stage_small!(py, px, pt, ms, Val($r), Val($σ), $V, Val($s)) : $ex)) - end - ex -end - # 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 @@ -72,7 +60,7 @@ for T in (Float32, Float64) (:cstage_last!, (:py, :px, :s), (r, σ, V) -> :(stage_last!(py, px, s, Val($r), Val($σ), $V))), (:cstage_last_range!, (:py, :px, :s, :qlo, :qhi), (r, σ, V) -> :(stage_last!(py, px, s, qlo, qhi, Val($r), Val($σ), $V))), (:cstage_group!, (:py, :px, :pt, :m), (r, σ, V) -> :(stage_first_group!(py, px, pt, m, Val($r), Val($σ), $V))), - (:cstage_small!, (:py, :px, :pt, :ms, :s), (r, σ, V) -> small_call(r, σ, V)), + (:cstage_small!, (:py, :px, :pt, :m), (r, σ, V) -> :(stage_small!(py, px, pt, m, Val($r), Val($σ), $V))), (:cstage_lanes!, (:py, :px, :pt, :m, :s0, :sx, :sy, :B), (r, σ, V) -> :(stage_lanes!(py, px, pt, m, s0, sx, sy, B, Val($r), Val($σ), $V))), (:cstage_last_lanes!, (:py, :px, :s0, :sx, :sy, :B), (r, σ, V) -> :(stage_last_lanes!(py, px, s0, sx, sy, B, Val($r), Val($σ), $V)))) for (name, args, call) in kinds @@ -130,7 +118,7 @@ const OP_MID_TAIL = 2 # q in vectors up to s - rem, the scalar remainder aft const OP_LAST = 3 const OP_LAST_TAIL = 4 const OP_GROUP = 5 # first stage vectorised over the butterflies p -const OP_SMALL = 6 # stride below the widest vector: vectorised over q + s·p (stage_small!) +const OP_SMALL = 6 # first stage over the butterflies p on the widest vector (stage_small!) struct StageCall op::Int @@ -152,7 +140,6 @@ struct RChain{T} <: Kernel{T} 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 - twg::Vector{Vector{Complex{T}}} # stage_small! with s > 1: w^(kp) at (k-1)ms + f, p = f ÷ s; else empty single::Vector{StageCall} # the stage calls of a single line end @@ -165,15 +152,13 @@ RChain{T}(n::Int, σ::Int, rs::AbstractVector{Int} = radices(n)) where {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) - twg = small_tables(ladder_ncomplex(T)[1], n, rs, tw) - RChain{T}(n, σ, collect(Int, rs), tw, tw1, twg, single_schedule(T, n, σ, collect(Int, rs))) + RChain{T}(n, σ, collect(Int, rs), tw, tw1, single_schedule(T, n, σ, collect(Int, rs))) end # The stage calls of a single line, with the rules of `chain_body` for S = 1: 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; power-of-two stages with a stride below the widest vector run -# stage_small! at full width. +# when its stride is 1, with stage_small! where that is wider. function single_schedule(::Type{T}, n::Int, σ::Int, rs::Vector{Int}) where {T} ncs = ladder_ncomplex(T) nv = length(ncs) @@ -196,8 +181,7 @@ function single_schedule(::Type{T}, n::Int, σ::Int, rs::Vector{Int}) where {T} rem = s % ncs[vi] c = codelet_index(T, σ, vi, r); c1 = codelet_index(T, σ, nv, r) m = L ÷ r - if i < K && small_stride(ncs[1], r, n, s) && (s > 1 || group == 0 || ncs[group] < ncs[1]) - # small strides at full width, where the grouped first stage would not be + if i == 1 && small_first(ncs, group, r, n) push!(calls, StageCall(OP_SMALL, codelet_index(T, σ, 1, r), c1, m, s, 0)) elseif i == K push!(calls, StageCall(rem == 0 ? OP_LAST : OP_LAST_TAIL, c, c1, m, s, rem)) @@ -299,7 +283,7 @@ function run_single!(py::Ptr{T}, px::Ptr{T}, pw1::Ptr{T}, pw2::Ptr{T}, P::RChain elseif op == OP_GROUP cstage_group!(dst, src, cptr(P.tw1), st.m, st.c) elseif op == OP_SMALL - cstage_small!(dst, src, cptr(st.s == 1 ? P.tw1 : P.twg[i]), st.m * st.s, st.s, st.c) + 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 @@ -373,13 +357,13 @@ 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, ntuple(i -> P.twg[i], Val(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, ntuple(i -> P.twg[i], Val(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 @@ -428,7 +412,7 @@ end # 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, ntuple(i -> P.twg[i], Val(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) diff --git a/test/runtests.jl b/test/runtests.jl index 66e75ba8..cd5f6849 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -754,16 +754,15 @@ thin(T, v) = FULL || T == Float64 ? v : v[1:2:end] end @testset "small-stride stages, $T" for T in (Float64, Float32) - # stage_small!: every power-of-two radix at every stride below the widest vector, in both - # chain forms, against the direct DFT + # stage_small!: every power-of-two radix as the first stage, in both chain forms, against + # the direct DFT p(a) = Ptr{T}(pointer(a)) C = FFTL.ladder_ncomplex(T)[1] - for s in 1:C-1, r in (2, 4, 8, 16), σ in (-1, 1) - C % s == 0 || continue - rs = s == 1 ? [r, 2C] : [s, r, 2C] + 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) - s == 1 || @test rc.single[2].op == FFTL.OP_SMALL + @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) From 3efcae67ead54887df0db5b7d1e45764eb7a2d2b Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 12 Nov 2026 06:16:59 +0000 Subject: [PATCH 6/8] Twiddles as memory operands and scaled stream addresses in the stage kernels MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The stages vectorised over q (stage!, stage_range!, stage_last!) need few registers besides the data in their inner loop: - ccmul_at multiplies by the twiddle at a pointer: re(w) and im(w) are scalar broadcasts, a·w = a·re(w) -/+ swap(a)·im(w) with the sign alternating over the lanes. On x86 both are memory operands. Every complex multiply outside the FCMA path has this form (ccmul_parts), so all kernels round alike and threaded and single-threaded transforms stay bit-identical. - The R input and R output streams are addressed as base + odd·stride·2^e with odd in (1, 3, 5, 7) and one base per 8 streams. The base pointer of the loop passes through llvm.launder.invariant.group and the scaling is a typed GEP; without both, loop strength reduction rewrites the addresses into one hoisted pointer per stream. A radix-16 stage has 15 twiddles and 32 streams; with twiddles in registers and one pointer per stream its inner loop reloaded vectors and pointers from the stack. Instructions per butterfly in the inner loop, AMD EPYC 9354 (AVX-512), Float64 (Float32 within 3), before -> after: radix 16 middle stage 284 -> 198 last stage 182 -> 151 radix 8 middle stage 100 -> 79 last stage 71 -> 58 Time, twiddle multiply alone (an earlier form with the x86 vfmaddsub intrinsics), one line, plan_fft, medians of 8 interleaved rounds, after/before: Float32 512-16384 0.88-0.93 (8192: 0.94 on the minima), Float64 1024, 4096, 8192, 16384 0.94-0.96, 32768 and above 0.98-1.01; the same code against itself differs by 1-5 % for most lengths. The time of both changes together is not verified. Not verified on aarch64. The test suite was not run on this commit; spot checks against a direct DFT and the threaded rfft equality pass. Recorded in docs/src/explanation/design-decisions.md. Co-Authored-By: Claude Fable 5.2 --- docs/src/explanation/design-decisions.md | 27 +++++++++++ src/complexvec.jl | 15 ++++-- src/kernels.jl | 59 +++++++++++++++++++----- 3 files changed, 85 insertions(+), 16 deletions(-) diff --git a/docs/src/explanation/design-decisions.md b/docs/src/explanation/design-decisions.md index be313d78..45324eb3 100644 --- a/docs/src/explanation/design-decisions.md +++ b/docs/src/explanation/design-decisions.md @@ -125,6 +125,33 @@ stream offsets) needs one kernel per length, which only the compiled `Chain` can `Chain` against the same kernels with run-time m and s, lengths 16–16384, both element types: ratios 0.90–1.04 without a trend, inside the spread of the measurement. Not verified on other machines. +## Stage kernels keep twiddles in memory and address streams from few registers + +**Decision (2026-11-12).** The stages vectorised over q (`stage!`, `stage_range!`, `stage_last!`) are +written so that the inner loop needs few registers besides the data: + +- The twiddle multiply reads re(w) and im(w) from the table as scalar broadcasts, + a·w = (a·re(w)) ∓ (swap(a)·im(w)) with the sign alternating over the lanes (`ccmul_at`). On x86 both + become memory operands of the multiply and the fused multiply. Every complex multiply outside the + FCMA path has this form (`ccmul_parts`), so all kernels round alike and threaded and single-threaded + transforms stay bit-identical. +- The R input and R output streams are addressed as base + odd·stride·2^e with odd in (1, 3, 5, 7) and + one base per 8 streams (`stream_addr`). The base pointer of the loop passes through + `llvm.launder.invariant.group` and the scaling is a typed GEP: without both, LLVM's loop strength + reduction rewrites the addresses into one hoisted pointer per stream. + +**Why.** A radix-16 stage has 15 twiddles, 16 input and 16 output streams. With the twiddles (and their +swapped copies) in registers and one pointer per stream, the inner loop reloaded vectors and pointers +from the stack. Instructions per butterfly in the inner loop, AMD EPYC 9354 (AVX-512), `Float64` +(`Float32` within 3), before → after: radix 16 middle stage 284 → 198, last stage 182 → 151; radix 8 +middle stage 100 → 79, last stage 71 → 58. + +Time: the twiddle multiply alone (in an earlier form written with the x86 `vfmaddsub` intrinsics), one +line, `plan_fft`, medians of 8 interleaved rounds, time after over time before: `Float32` 512–16384: +0.88–0.93 (8192: 0.94 on the minima), `Float64` 1024, 4096, 8192, 16384: 0.94–0.96, lengths of 32768 +and above 0.98–1.01; the same code against itself differs by 1–5 % for most lengths. The time of both +changes together is not verified, and neither change is verified on aarch64. + ## Codelets of the radices 17–31 are compiled on first use **Decision (2026-11-10).** The codelets of the radices 17, 19, 23, 29 and 31 are not part of the package diff --git a/src/complexvec.jl b/src/complexvec.jl index c02fdf86..3e7bd9ea 100644 --- a/src/complexvec.jl +++ b/src/complexvec.jl @@ -98,6 +98,7 @@ if HAS_FCMA @inline fcmla90(c::Vec{$N,$T}, a::Vec{$N,$T}, b::Vec{$N,$T}) = Vec(ccall($("llvm.aarch64.neon.vcmla.rot90." * suf), llvmcall, $VT, ($VT, $VT, $VT), c.data, a.data, b.data)) @inline ccmul(a::Vec{$N,$T}, w::Vec{$N,$T}) = fcmla90(fcmla0(zero(Vec{$N,$T}), a, w), a, w) + @inline ccmul_at(a::Vec{$N,$T}, p::Ptr{$T}) = ccmul(a, cbroadcast(Vec{$N,$T}, p)) end end else @@ -129,11 +130,17 @@ else end end end - # a·w = (a.re·w.re - a.im·w.im, a.re·w.im + a.im·w.re) - @inline function ccmul(a::Vec{N,T}, w::Vec{N,T}) where {N,T} - t = dup_im(a) * swap_reim(w) - blend_reim(muladd(dup_re(a), w, -t), muladd(dup_re(a), w, t)) + # a·w = (a.re·w.re - a.im·w.im, a.im·w.re + a.re·w.im) from re(w) and im(w), each repeated over + # the two lanes of its complex value + @inline function ccmul_parts(a::Vec{N,T}, wre::Vec{N,T}, wim::Vec{N,T}) where {N,T} + t = swap_reim(a) * wim + blend_reim(muladd(a, wre, -t), muladd(a, wre, t)) end + @inline ccmul(a::Vec{N,T}, w::Vec{N,T}) where {N,T} = ccmul_parts(a, dup_re(w), dup_im(w)) + # a·w for the complex value w at p, repeated over all complex lanes: re(w) and im(w) are scalar + # broadcasts, which x86 folds into the multiplies as memory operands + @inline ccmul_at(a::Vec{N,T}, p::Ptr{T}) where {N,T} = + ccmul_parts(a, Vec{N,T}(unsafe_load(p)), Vec{N,T}(unsafe_load(p, 2))) end # (re, im) repeated over all complex lanes of the register diff --git a/src/kernels.jl b/src/kernels.jl index ed295f13..46b26928 100644 --- a/src/kernels.jl +++ b/src/kernels.jl @@ -18,19 +18,22 @@ stage!(py::Ptr{T}, px::Ptr{T}, pt::Ptr{T}, m::Int, s::Int, ::Val{R}, ::Val{σ}, function stage_body(::Type{T}, R::Int, σ::Int, ::Type{VT}) where {T,VT} stmts, a, X = codelet(R, σ) C = ncomplex(VT); B = 2sizeof(T) - loads = [:($(a[j]) = vload($VT, px + (q + s * (p + $(j - 1) * m)) * $B)) for j in 1:R] - twl = [:($(Symbol("w", k)) = cbroadcast($VT, pt + ($(R - 1) * p + $(k - 1)) * $(2sizeof(T)))) for k in 1:R-1] - stores = Any[:(vstore($(X[1]), py + o * $B))] + loads = [:($(a[j]) = vload($VT, $(stream_addr(:xb, :xr, j - 1)))) for j in 1:R] + twl = [:($(Symbol("w", k)) = pt + ($(R - 1) * p + $(k - 1)) * $(2sizeof(T))) for k in 1:R-1] + stores = Any[:(vstore($(X[1]), yb0))] for k in 1:R-1 - push!(stores, :(vstore(ccmul($(X[k+1]), $(Symbol("w", k))), py + (o + $k * s) * $B))) + push!(stores, :(vstore(ccmul_at($(X[k+1]), $(Symbol("w", k))), $(stream_addr(:yb, :yr, k))))) end quote + $(stream_strides(:xr, R, :(s * m * $B))...) + $(stream_strides(:yr, R, :(s * $B))...) for p in plo:phi-1 $(twl...) for q in 0:$C:s-1 + $(stream_bases(:xb, :xr, R, :(px + (q + s * p) * $B))...) $(loads...) $(stmts...) - o = q + s * $R * p + $(stream_bases(:yb, :yr, R, :(py + (q + s * $R * p) * $B))...) $(stores...) end end @@ -60,22 +63,51 @@ 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")) + ir = "define ptr @entry(ptr %0, i64 %1) #0 {\ntop:\n %p = getelementptr $ty, ptr %0, i64 %1\n ret ptr %p\n}\nattributes #0 = { alwaysinline }" + @eval @inline scaled_ptr(b::Ptr{T}, i::Int, ::Val{$(e - 1)}) where {T} = Base.llvmcall(($ir, "entry"), Ptr{T}, Tuple{Ptr{T},Int}, b, i) +end + +# The pointer p, opaque to the loop optimisations: an induction pointer passed through here is +# not rewritten into one pointer per stream. +@inline opaque_ptr(p::Ptr{T}) where {T} = Base.llvmcall( + ("declare ptr @llvm.launder.invariant.group.p0(ptr)\ndefine ptr @entry(ptr %0) #0 {\ntop:\n %p = call ptr @llvm.launder.invariant.group.p0(ptr %0)\n ret ptr %p\n}\nattributes #0 = { alwaysinline }", "entry"), + Ptr{T}, Tuple{Ptr{T}}, p) + +# Address of stream j (zero-based) of streams r bytes apart: base + odd·r·2^e with odd in +# (1, 3, 5, 7) and one base per 8 streams, so that a few registers reach every stream. +function stream_addr(base, r, j) + b = Symbol(base, j ÷ 8); k = j % 8 + k == 0 && return b + e = trailing_zeros(k) + :(scaled_ptr($b, $(Symbol(r, k >> e)), Val($e))) +end +stream_bases(base, r, R, first) = + [i == 0 ? :($(Symbol(base, 0)) = opaque_ptr($first)) : :($(Symbol(base, i)) = scaled_ptr($(Symbol(base, i - 1)), $(Symbol(r, 1)), Val(3))) for i in 0:(R-1)÷8] +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]) + 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) - loads = [:($(a[j]) = vload($VT, px + (q + s * (p + $(j - 1) * m)) * $B)) for j in 1:R] - twl = [:($(Symbol("w", k)) = cbroadcast($VT, pt + ($(R - 1) * p + $(k - 1)) * $(2sizeof(T)))) for k in 1:R-1] - stores = Any[:(vstore($(X[1]), py + o * $B))] + loads = [:($(a[j]) = vload($VT, $(stream_addr(:xb, :xr, j - 1)))) for j in 1:R] + twl = [:($(Symbol("w", k)) = pt + ($(R - 1) * p + $(k - 1)) * $(2sizeof(T))) for k in 1:R-1] + stores = Any[:(vstore($(X[1]), yb0))] for k in 1:R-1 - push!(stores, :(vstore(ccmul($(X[k+1]), $(Symbol("w", k))), py + (o + $k * s) * $B))) + push!(stores, :(vstore(ccmul_at($(X[k+1]), $(Symbol("w", k))), $(stream_addr(:yb, :yr, k))))) end quote + $(stream_strides(:xr, R, :(s * m * $B))...) + $(stream_strides(:yr, R, :(s * $B))...) for p in plo:phi-1 $(twl...) for q in qlo:$C:qhi-1 + $(stream_bases(:xb, :xr, R, :(px + (q + s * p) * $B))...) $(loads...) $(stmts...) - o = q + s * $R * p + $(stream_bases(:yb, :yr, R, :(py + (q + s * $R * p) * $B))...) $(stores...) end end @@ -95,12 +127,15 @@ stage_last!(py::Ptr{T}, px::Ptr{T}, s::Int, ::Val{R}, ::Val{σ}, ::Type{VT}) whe function stage_last_body(::Type{T}, R::Int, σ::Int, ::Type{VT}) where {T,VT} stmts, a, X = codelet(R, σ) C = ncomplex(VT); B = 2sizeof(T) - loads = [:($(a[j]) = vload($VT, px + (q + s * $(j - 1)) * $B)) for j in 1:R] - stores = [:(vstore($(X[k+1]), py + (q + $k * s) * $B)) for k in 0:R-1] + loads = [:($(a[j]) = vload($VT, $(stream_addr(:xb, :r, j - 1)))) for j in 1:R] + stores = [:(vstore($(X[k+1]), $(stream_addr(:yb, :r, k)))) for k in 0:R-1] quote + $(stream_strides(:r, R, :(s * $B))...) for q in qlo:$C:qhi-1 + $(stream_bases(:xb, :r, R, :(px + q * $B))...) $(loads...) $(stmts...) + $(stream_bases(:yb, :r, R, :(py + q * $B))...) $(stores...) end nothing From 3598e95e42e058810b43e9319349a0f3ecd6ff44 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 12 Nov 2026 06:39:39 +0000 Subject: [PATCH 7/8] Shorten the design decisions on the single-line stage kernels The two entries state the decision, the measured evidence and what was not done, without the detail that the code and the commit log carry. The entry on twiddles and stream addresses now reports both changes together (AMD EPYC 9354, one line, plan_fft, medians of 8 interleaved rounds against the kernels before): lengths 512-16384 take 0.79-0.94 of the time for Float32 and 0.82-0.92 for Float64, lengths from 32768 on are unchanged; the same code against itself differs by 1-5 %. The quick test selection passes with these kernels. Co-Authored-By: Claude Fable 5.2 --- docs/src/explanation/design-decisions.md | 73 +++++++----------------- 1 file changed, 22 insertions(+), 51 deletions(-) diff --git a/docs/src/explanation/design-decisions.md b/docs/src/explanation/design-decisions.md index 45324eb3..e9a6ada4 100644 --- a/docs/src/explanation/design-decisions.md +++ b/docs/src/explanation/design-decisions.md @@ -97,60 +97,31 @@ specialisation, 130 MB of the 146 MB image. See [Package size and load time](@re ## The first stage of a single line runs at full vector width -**Decision (2026-11-12).** The first stage of a single line (stride 1, power-of-two radix, not the last -stage) runs `stage_small!` where the widest vector holds more complex values than the grouped first -stage handles (C = 8, `Float32` on AVX-512): it is vectorised over the butterflies p, so every input is -one contiguous load, the twiddles come from the grouped first-stage table, and the R outputs are -interleaved in registers into contiguous stores. The generator checks the permutation of the -interleave. Both chain forms use it; the kernels (radices 2, 4, 8, 16 × sign, widest vector only) are -part of the fixed set of plain methods. - -**Why.** That stage ran the grouped first stage on 4 complex values. `Float64` has no such stage on -AVX-512. Intel Core i7-1185G7, one line, `Float32`, `RChain`, radix sequences of 16s and 8s, µs without -→ with: -128: 0.093 → 0.069, 256: 0.176 → 0.131, 1024: 0.81 → 0.61, 2048: 1.60 → 1.35, 4096: 3.75 → 3.29, -8192: 7.85 → 6.69, 16384: 18.0 → 15.4 (the same code against itself differs by up to 3 %). - -**Not done.** The same scheme for later stages with a stride s = 2 or 4 below C (vectorised over -q + s·p). It needs a per-lane twiddle table for every such stage, one kernel per stride and a stride -switch in the `RChain` dispatch. Only `Float32` lengths whose radices start with 4, 4 have such a stage -(s = 4; 625 lengths up to 70000, 65536 the only power of two), none with s = 2. AMD EPYC 9354, one line, -`Float32`, `RChain`, time without the s = 4 stage over time with it, medians of 10 interleaved rounds: -192: 1.08, 320: 1.06, 768: 1.03, 960: 1.03, 3072: 0.96, 4800: 0.99, 20480: 0.97; 65536: 0.96 and -262144: 0.97 from 5 rounds. The same code against itself differs by 5–20 % between rounds on that -machine, so only the 6–8 % at 192 and 320 (the same in a second series) stands out. - -**Not done.** Passing m and s to the stage kernels as constants of a single line (immediates for the -stream offsets) needs one kernel per length, which only the compiled `Chain` can have. Measured in -`Chain` against the same kernels with run-time m and s, lengths 16–16384, both element types: ratios -0.90–1.04 without a trend, inside the spread of the measurement. Not verified on other machines. +**Decision (2026-11-12).** Where the widest vector holds 8 complex values (`Float32` on AVX-512), the +first stage of a single line runs `stage_small!`: vectorised over the butterflies, with the outputs +interleaved in registers into contiguous stores. It exists for the power-of-two radices and belongs to +the fixed set of plain methods. + +**Why.** The grouped first stage handles at most 4 complex values per vector. Intel Core i7-1185G7, +one line, `Float32`, lengths 128–16384 with radices of 16s and 8s: 12–26 % faster (the same code +against itself differs by up to 3 %). + +**Not done.** The same for later stages with a stride of 2 or 4. It needs a per-lane twiddle table and +one kernel per stride, and only `Float32` lengths whose radices start with 4, 4 have such a stage. On an +AMD EPYC 9354 it gained 6–8 % at 192 and 320 and nothing measurable at other lengths. ## Stage kernels keep twiddles in memory and address streams from few registers -**Decision (2026-11-12).** The stages vectorised over q (`stage!`, `stage_range!`, `stage_last!`) are -written so that the inner loop needs few registers besides the data: - -- The twiddle multiply reads re(w) and im(w) from the table as scalar broadcasts, - a·w = (a·re(w)) ∓ (swap(a)·im(w)) with the sign alternating over the lanes (`ccmul_at`). On x86 both - become memory operands of the multiply and the fused multiply. Every complex multiply outside the - FCMA path has this form (`ccmul_parts`), so all kernels round alike and threaded and single-threaded - transforms stay bit-identical. -- The R input and R output streams are addressed as base + odd·stride·2^e with odd in (1, 3, 5, 7) and - one base per 8 streams (`stream_addr`). The base pointer of the loop passes through - `llvm.launder.invariant.group` and the scaling is a typed GEP: without both, LLVM's loop strength - reduction rewrites the addresses into one hoisted pointer per stream. - -**Why.** A radix-16 stage has 15 twiddles, 16 input and 16 output streams. With the twiddles (and their -swapped copies) in registers and one pointer per stream, the inner loop reloaded vectors and pointers -from the stack. Instructions per butterfly in the inner loop, AMD EPYC 9354 (AVX-512), `Float64` -(`Float32` within 3), before → after: radix 16 middle stage 284 → 198, last stage 182 → 151; radix 8 -middle stage 100 → 79, last stage 71 → 58. - -Time: the twiddle multiply alone (in an earlier form written with the x86 `vfmaddsub` intrinsics), one -line, `plan_fft`, medians of 8 interleaved rounds, time after over time before: `Float32` 512–16384: -0.88–0.93 (8192: 0.94 on the minima), `Float64` 1024, 4096, 8192, 16384: 0.94–0.96, lengths of 32768 -and above 0.98–1.01; the same code against itself differs by 1–5 % for most lengths. The time of both -changes together is not verified, and neither change is verified on aarch64. +**Decision (2026-11-12).** In the stages vectorised over q, the twiddle multiply reads re(w) and im(w) +from the table (`ccmul_at`; memory operands on x86), and the R input and R output streams are addressed +as base + odd·stride·2^e from one base per 8 streams (`stream_addr`). The base pointer is hidden from +LLVM's loop optimiser (`opaque_ptr`), which otherwise goes back to one pointer per stream. + +**Why.** A radix-16 stage has 15 twiddles and 32 streams, more than fit in registers, so its inner loop +reloaded vectors and pointers from the stack. AMD EPYC 9354 (AVX-512), instructions per butterfly: +radix 16 middle stage 284 → 198, radix 8 100 → 79. One line, `plan_fft`, medians of 8 interleaved +rounds: lengths 512–16384 take 0.79–0.94 (`Float32`) and 0.82–0.92 (`Float64`) of the time, lengths +from 32768 on are unchanged (the same code against itself differs by 1–5 %). Not verified on aarch64. ## Codelets of the radices 17–31 are compiled on first use From 612b245585759812b086a69599f4849156327461 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Thu, 12 Nov 2026 07:08:19 +0000 Subject: [PATCH 8/8] Radix 16 for single lines only below an L1 set stride between streams MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit A radix-16 stage reads 16 streams n/16 elements apart. Once that spacing reaches the L1 set stride (4 KB) the streams share one cache set, and the sequences of 8s are faster. Single lines on AVX-512 therefore take radix 16 up to 4096 for Float32 and 2048 for Float64; 8192 and 16384 (and Float64 4096) use 8s with a 4 or two 4s again. AMD EPYC 9354, one line, all buffers at fixed page offsets (staggered by 1 KB), µs, sequence of 8s -> sequence of 16s: Float32 4096 3.56 -> 3.13 8192 7.48 -> 7.93 16384 15.95 -> 16.9 Float64 4096 6.59 -> 8.79 8192 14.81 -> 15.32 16384 31.76 -> 35.3 The ranking is the same with all buffers at one page offset. plan_fft against main (103b9b6), medians of 8 interleaved rounds, time of this branch over main: Float64 4096 0.86 (1.23 before this rule), 8192 0.94, 16384 0.95; Float32 8192 0.85, 16384 0.87. The same code against itself differs by 1-5 % for most lengths, up to 25 % for a few. The quick test selection passes. Not measured on Intel with this rule. Recorded in docs/src/explanation/design-decisions.md. Co-Authored-By: Claude Fable 5.2 --- docs/src/explanation/design-decisions.md | 6 ++++++ src/factor.jl | 13 +++++++------ test/data/radices-x86_64-32regs.csv | 6 +++--- 3 files changed, 16 insertions(+), 9 deletions(-) diff --git a/docs/src/explanation/design-decisions.md b/docs/src/explanation/design-decisions.md index e9a6ada4..437bef8a 100644 --- a/docs/src/explanation/design-decisions.md +++ b/docs/src/explanation/design-decisions.md @@ -35,6 +35,12 @@ and the large losses of an estimate-only planner are cache effects the algorithm A per-machine calibration or an optional measuring mode can follow if the rule leaves too much on the table. +One such cache effect is part of the rule: single lines on AVX-512 take radix 16 only while the 16 +streams of a stage are less than an L1 set stride (4 KB) apart, that is up to 4096 for `Float32` and +2048 for `Float64`. Beyond that the streams share one cache set and radix 8 wins. AMD EPYC 9354, one +line, all buffers at fixed page offsets, time of the radix-16 sequence over the radix-8 sequence: +`Float32` 4096: 0.88, 8192: 1.06, 16384: 1.06; `Float64` 4096: 1.33, 8192: 1.03, 16384: 1.11. + ## 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 diff --git a/src/factor.jl b/src/factor.jl index 334cef30..869fd86f 100644 --- a/src/factor.jl +++ b/src/factor.jl @@ -7,9 +7,10 @@ # up to 4096 on the M1; on Tiger Lake (AVX-512) 7–30 % from 4096 on (S = 1) and up to 45 % in # batches, within noise below; on Zen 3 (AVX2, 16 # registers) it is 3–25 % slower from 256 on. -# Single lines on AVX-512 take radix 16 as often as possible and 8s for the rest (Tiger Lake, -# 128–16384: Float64 2–15 % faster, Float32 between 6 % slower and 8 % faster). Not measured -# on NEON. +# Single lines on AVX-512 take radix 16 as often as possible and 8s for the rest while the 16 +# streams of a stage are less than an L1 set stride apart (Float32 up to 4096, Float64 up to +# 2048); beyond that they share one cache set and the 8s win (Zen 4: 3–33 % at 4096–16384). +# Not measured on NEON. # Otherwise radix 4 throughout, one 8 (or 2) for an odd exponent; radix 16 and above lose # (register spills, cache-set conflicts at large strides). # - aarch64: a factor 3 is merged with a 4 into 12, else with a 2 into 6; a factor 5 with a spare @@ -62,7 +63,7 @@ Radix sequence for a chain of length `n` applied to `S` interleaved sequences of `esize` bytes. """ function radices(n::Int, S::Int = 1; esize::Int = 16) - r = heuristic_radices(n, S) + r = heuristic_radices(n, S, esize) S > LARGE_STRIDE && set_conflict(r, n, S, esize) || return r alt = packed_radices(n) set_conflict(alt, n, S, esize) ? r : alt @@ -97,7 +98,7 @@ function packed_radices(n::Int) [reverse!(sort!(filter(ispow2, bins))); reverse!(sort!(filter(!ispow2, bins)))] end -function heuristic_radices(n::Int, S::Int) +function heuristic_radices(n::Int, S::Int, esize::Int = 16) n == 1 && return Int[] merge = merge_composite(S) f, rest = factorize(n) @@ -108,7 +109,7 @@ function heuristic_radices(n::Int, S::Int) end if RADIX8_POW2 && ispow2(n) && 64 < n <= 16384 && S <= LARGE_STRIDE k = trailing_zeros(n) - if HAS_AVX512 && S == 1 + if HAS_AVX512 && S == 1 && (n ÷ 16) * esize < 1 << L1_SET_LOG2 # as many 16s as possible, the rest in 8s: 4t ≤ k, t ≡ k (mod 3) t = k % 3 + 3 * ((k ÷ 4 - k % 3) ÷ 3) return vcat(fill(16, t), fill(8, (k - 4t) ÷ 3)) diff --git a/test/data/radices-x86_64-32regs.csv b/test/data/radices-x86_64-32regs.csv index 291c0d04..826e8734 100644 --- a/test/data/radices-x86_64-32regs.csv +++ b/test/data/radices-x86_64-32regs.csv @@ -882,11 +882,11 @@ 1024,1024,4 4 4 4 4 2048,1,16 16 8 2048,1024,8 4 4 4 4 -4096,1,16 16 16 +4096,1,8 8 8 8 4096,1024,4 4 4 4 4 4 -8192,1,16 8 8 8 +8192,1,4 4 8 8 8 8192,1024,8 4 4 4 4 4 -16384,1,16 16 8 8 +16384,1,4 8 8 8 8 16384,1024,4 4 4 4 4 4 4 32768,1,8 4 4 4 4 4 4 32768,1024,8 4 4 4 4 4 4