Remove the pruning of pass plans (nonzero, needed)

README.md

@@ -103,10 +103,9 @@ only dependency is the Julia runtime shipped next to it. On macOS arm64 the libr
 the runtime 10 MB. The run-time plans the library uses (`FFTL.RPlan`, `FFTL.RRFFTPlan`,
 `FFTL.RR2RPlan`) are public, so Julia programs that use them can be trimmed the same way.
 
-**Entry points for special variants.** `FFTL.plan_pass(v, d; nonzero, needed)` transforms all
-lines of a strided view along one dimension and can skip zero inputs and unneeded outputs
-(pruning); `FFTL.plan_r2r_pass` does the same for real-to-real kinds. Libraries that build their
-own transforms from one-dimensional passes use these; NFFT.jl is one of them.
+**Entry points for special variants.** `FFTL.plan_pass(v, d)` transforms all lines of a strided view
+along one dimension, and `FFTL.plan_r2r_pass` does the same for real-to-real kinds. Libraries that build
+their own transforms from one-dimensional passes use these; NFFT.jl is one of them.
 
 ## Performance
 
@@ -157,7 +156,7 @@ describes the setup and what a change must come with.
 
 ## Tests
 
-`] test FFTL` runs a quick selection that covers every feature (about 1.5 min). Set
+`] test FFTL` runs a quick selection that covers every feature (about 5 min with two threads). Set
 `FFTL_TEST_FULL=true` for the full lists and `FFTL_TEST_VERBOSE=true` for the time of every test
 set. CI runs the quick tests on Linux and macOS for pushes to `main` and pull requests, the full
 lists and the oldest supported Julia (1.10) once a week, and builds and tests the C library.

docs/AGENTS.md

@@ -29,12 +29,12 @@ rules there apply to everything written here as well. Everything in this folder
   with the DCT, Poisson with the DST (second-order convergence shown), fast convolution. Each has
   at least one figure (CairoMakie) next to its printed results.
 - **How-to Guides** (`src/howto/`): coming from FFTW.jl, plans and reuse, transforms along dimensions
-  (including `plan_pass` with pruning, which NFFT.jl uses), real data, cosine/sine/Hartley transforms,
+  (including `plan_pass`, which NFFT.jl uses), real data, cosine/sine/Hartley transforms,
   threads, number types, precompilation, GPU, C library (libfftl), compiling ahead of time (the run-time
   plans `RPlan`, `RRFFTPlan`, `RR2RPlan`), contribute.
 - **Explanation** (`src/explanation/`): the DFT and its conventions, FFTW and FFTL (history, papers, what
   FFTL takes from FFTW; written from papers and documentation only), algorithms (Stockham chains,
-  codelets, pruning, Bluestein, multidimensional without transposition, four-step, real and real-to-real
+  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
   size and load time, performance (measurement protocol).

docs/src/explanation/algorithms.md

@@ -38,13 +38,6 @@ The generated code runs on complex SIMD vectors (`src/complexvec.jl`), with the
 the element type. On aarch64 the complex add-and-rotate and the complex multiply are single
 instructions (`FCADD`, `FCMLA`).
 
-## Pruned codelets
-
-When the input is known to be zero outside some range, or only some outputs are used, the
-generator drops the terms that are zero and the statements that feed no needed output. A pass
-that transforms a zero-padded array therefore does not read the padding, and an NFFT that needs
-only part of the spectrum does not compute the rest. See [transforms along dimensions](@ref "Transforms along dimensions").
-
 ## Bluestein's algorithm
 
 Lengths with a prime factor above 31 go through Bluestein's identity. With ``jk = (j^2 + k^2 - (k-j)^2)/2``,
@@ -156,6 +149,6 @@ listed on [the cosine, sine and Hartley page](@ref "Cosine, sine and Hartley tra
 ## Passes
 
 `plan_pass` applies one 1D transform to all lines of a strided array with the same kernel. It
-lets a library assemble its own multidimensional transform one dimension at a time and declare,
-per dimension, which inputs are zero and which outputs it needs; both declarations use the
-pruned codelets above. NFFT.jl transforms its oversampled grid this way.
+lets a library assemble its own multidimensional transform one dimension at a time. NFFT.jl
+transforms its oversampled grid this way; it skips the boxes of the grid that it does not need
+itself, so the passes themselves carry no pruning.

docs/src/explanation/image-size.md

@@ -45,8 +45,8 @@ later Julia versions, so re-run the count below after an upgrade.
 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
-   functions the generators call. The `@generated` methods stay as the fallback for everything else
-   (pruned stages, other chains), so behaviour is unchanged. A plain method keeps its (compressed) source
+   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
    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
@@ -90,8 +90,8 @@ Further steps that were considered:
 | no `Float32` in the workload | 16 MB (at 73 MB) | 6 s at the first `Float32` call |
 
 Replacing the preference by environment variables did not shorten the load, and Julia does not notice a
-changed variable, so a stale image would stay in use. Removing the pruned kernels from the workload was not
-needed: no pruned kernel is in the image.
+changed variable, so a stale image would stay in use. The pass plans do not take zero masks any more (their
+kernels were `@generated` with a mask per segment), so no pruned kernel exists to store.
 
 ## Checking it
 

docs/src/howto/contribute.md

@@ -33,9 +33,9 @@ describe what the code does now, not how it got there.
 using Pkg; Pkg.test("FFTL")
 ```
 
-This runs a quick selection that covers every feature. On a 256-thread AMD EPYC 7702 with two Julia
-threads it takes about 40 minutes, so run it in the background or on the CI runner. Two environment
-variables change it:
+This runs a quick selection that covers every feature. With two Julia threads it takes about 5 minutes
+on a 256-thread AMD EPYC 7702, and it is the check to run before every push. Two environment variables
+change it:
 
 | variable | effect |
 |:--|:--|

docs/src/howto/dimensions.md

@@ -38,27 +38,8 @@ a 3D array costs little more than one for a 1D line.
 
 `FFTL.plan_pass` transforms along one dimension of a strided view, applied to all lines at once.
 It is the building block of the transforms above, and the entry point for libraries that build
-special variants of the transform themselves, such as NFFT.jl on its oversampled grid. Two
-optional declarations let it skip work:
-
-- `nonzero` says that each line is zero outside the given ranges of indices along `d`. The first
-  stage then does not read the zeros.
-- `needed` says that only the given output indices are used. The last stage computes only those.
-
-The result on the declared outputs is identical to the full transform. Here, the input is zero
-outside its first 16 of 64 entries, and only the first 8 outputs are needed:
-
-```@example dims
-n = 64
-V = zeros(ComplexF64, n, 10)
-V[1:16, :] .= randn(ComplexF64, 16, 10)
-
-full = fft(V, 1)
-pass = FFTL.plan_pass(V, 1; nonzero = (1:16,), needed = (1:8,))
-W = similar(V)
-mul!(W, pass, V)
-maximum(abs, W[1:8, :] - full[1:8, :])
-```
+special variants of the transform themselves, such as NFFT.jl on its oversampled grid. The plan is
+made for the size and strides of the view, and applied with `mul!` or `apply!`.
 
 Passes are single-threaded by default. A caller that already threads over many arrays keeps
 `threads = 1`; see [Threads](@ref "Threads").

docs/src/howto/fftw.md

@@ -83,7 +83,7 @@ roundtrip(p, x) ≈ x
 - Number types beyond `Float32` and `Float64`, computed in their own precision
   ([Number types](@ref "Number types")).
 - GPU arrays with the same functions ([GPU](gpu.md)).
-- Passes over one dimension of a strided array, with pruning
+- Passes over one dimension of a strided array
   ([Transforms along dimensions](@ref "Passes: one dimension of a strided array")).
 - A C library compiled from the Julia code, with FFTW3's C and Fortran interfaces
   ([libfftl](@ref "C library (libfftl)")).

docs/src/index.md

@@ -45,8 +45,8 @@ library needs no Julia installation. The run-time plans it is built from are pub
 programs that use them can be trimmed as well.
 
 **Entry points for special variants.** [Passes](@ref "Passes: one dimension of a strided array")
-transform one dimension of a strided array and can skip zero inputs and unneeded outputs.
-Libraries that assemble their own transforms from such passes, such as NFFT.jl, use them.
+transform one dimension of a strided array. Libraries that assemble their own transforms from such
+passes, such as NFFT.jl, use them.
 
 **The interface of the ecosystem.** FFTL exports FFTW.jl's names with FFTW.jl's conventions and
 scalings, so `using FFTL` can take the place of `using FFTW`, and it is a backend of

docs/src/reference/platforms.md

@@ -15,7 +15,7 @@ available.
 | DST-II | ✓ | ✓ | ✓ | ✓ |
 | DST-I, DST-III, DST-IV | ✓ | ✓ | ◐ partly tuned | ✓ |
 | DHT, `R2HC`, `HC2R` | ✓ | ✓ | ◐ not tuned | ✓ |
-| batched 1D passes (`plan_pass`, `plan_r2r_pass`) with pruning | ✓ | – | – | – |
+| batched 1D passes (`plan_pass`, `plan_r2r_pass`) | ✓ | – | – | – |
 | any rank and any subset of dimensions | ✓ | ✓ | ✓ | ✓ |
 | any length (Bluestein above prime factor 31) | ✓ | ✓ | ✓ | ✓ |
 

src/codegen.jl

@@ -11,8 +11,6 @@
 # Odd primes use the symmetric form
 #   X_k, X_{r-k} = A ± σ·i·B,  A = a_0 + Σ_j (a_j + a_{r-j}) cos(2πjk/r),  B = Σ_j (a_j - a_{r-j}) sin(2πjk/r).
 #
-# Pruned codelets: inputs known to be zero are the term `nothing` (sums and products simplify),
-# and only the needed outputs are computed (statements they do not depend on are dropped).
 
 mutable struct CodeGen
   stmts::Vector{Any}
@@ -29,13 +27,9 @@ function emit!(g::CodeGen, ex)
 end
 
 rotate(t, k) = (t[1], mod(t[2] + k, 4))
-rotate(::Nothing, k) = nothing
 # rotation units of one factor σ·i: -i = (-i)^1, +i = (-i)^3
 unit(g::CodeGen) = g.σ == -1 ? 1 : 3
 
-plus!(g, ::Nothing, b) = b
-plus!(g, a, ::Nothing) = a
-plus!(g, ::Nothing, ::Nothing) = nothing
 function plus!(g, a, b)
   d = mod(b[2] - a[2], 4)
   ex = d == 0 ? :($(a[1]) + $(b[1])) :
@@ -46,7 +40,6 @@ function plus!(g, a, b)
 end
 minus!(g, a, b) = plus!(g, a, rotate(b, 2))
 
-force!(g, ::Nothing) = emit!(g, :(zero(VT)))
 function force!(g, t)
   t[2] == 0 && return t[1]
   t[2] == 2 && return emit!(g, :(-$(t[1])))
@@ -55,7 +48,6 @@ function force!(g, t)
 end
 
 # t · W_N^e with constant e, N
-twiddle!(g, ::Nothing, e, N) = nothing
 function twiddle!(g, t, e, N)
   e = mod(e, N)
   e == 0 && return t
@@ -70,37 +62,34 @@ function twiddle!(g, t, e, N)
   (emit!(g, :(ccmul($(t[1]), cconst(VT, $c, $s)))), t[2])
 end
 
-# sums of plain symbols where `nothing` is zero
-radd!(g, x, y) = x === nothing ? y : y === nothing ? x : emit!(g, :($x + $y))
-rsub!(g, x, y) = y === nothing ? x : x === nothing ? emit!(g, :(-$y)) : emit!(g, :($x - $y))
-rmuladd!(g, x, c, acc) = x === nothing ? acc : acc === nothing ? emit!(g, :($c * $x)) : emit!(g, :(muladd($x, $c, $acc)))
+radd!(g, x, y) = emit!(g, :($x + $y))
+rsub!(g, x, y) = emit!(g, :($x - $y))
+rmuladd!(g, x, c, acc) = emit!(g, :(muladd($x, $c, $acc)))
+scale!(g, c, x) = emit!(g, :($c * $x))
 
 function prime!(g, a)
   r = length(a); h = (r - 1) ÷ 2
-  a = [x === nothing ? nothing : force!(g, x) for x in a]
+  a = [force!(g, x) for x in a]
   sp = [radd!(g, a[j+1], a[r-j+1]) for j in 1:h]
   dm = [rsub!(g, a[j+1], a[r-j+1]) for j in 1:h]
   X = Vector{Any}(undef, r)
-  acc = nothing
-  for j in 1:h
+  acc = sp[1]
+  for j in 2:h
     acc = radd!(g, acc, sp[j])
   end
-  x0 = radd!(g, a[1], acc)
-  X[1] = x0 === nothing ? nothing : (x0, 0)
+  X[1] = (radd!(g, a[1], acc), 0)
   for k in 1:h
     A = a[1]
     for j in 1:h
       A = rmuladd!(g, sp[j], cospi(2j * k / r), A)
     end
-    B = nothing
-    for j in 1:h
+    B = scale!(g, sinpi(2 * k / r), dm[1])
+    for j in 2:h
       B = rmuladd!(g, dm[j], sinpi(2j * k / r), B)
     end
     # X_k = A + σ·i·B and X_(r-k) = A - σ·i·B; -i·B is the term (B, 1), +i·B the term (B, 3)
-    tA = A === nothing ? nothing : (A, 0)
-    tB = B === nothing ? nothing : (B, 0)
-    X[k+1] = plus!(g, tA, rotate(tB, g.σ == -1 ? 1 : 3))
-    X[r-k+1] = plus!(g, tA, rotate(tB, g.σ == -1 ? 3 : 1))
+    X[k+1] = plus!(g, (A, 0), rotate((B, 0), g.σ == -1 ? 1 : 3))
+    X[r-k+1] = plus!(g, (A, 0), rotate((B, 0), g.σ == -1 ? 3 : 1))
   end
   X
 end
@@ -133,33 +122,15 @@ function dft!(g::CodeGen, a::Vector)
 end
 
 """
-    codelet(r, σ, zin = 0, need = all) -> (stmts, inputs, outputs)
+    codelet(r, σ) -> (stmts, inputs, outputs)
 
 Statements computing a size-`r` DFT with sign `σ` from the input symbols `a1…ar` into the
-returned output symbols. Inputs whose bit is set in `zin` (bit j-1 for input j) are zero and not
-read; outputs whose bit is clear in `need` are not computed (their entry is `nothing`).
+returned output symbols.
 """
-function codelet(r::Int, σ::Int, zin::Integer = 0, need::Integer = (1 << r) - 1)
+function codelet(r::Int, σ::Int)
   g = CodeGen(σ)
   a = [Symbol("a", j) for j in 1:r]
-  X = dft!(g, [((zin >> (j - 1)) & 1 == 1 ? nothing : (s, 0)) for (j, s) in enumerate(a)])
-  out = [(need >> (k - 1)) & 1 == 1 ? force!(g, X[k]) : nothing for k in 1:r]
-  prune_statements!(g.stmts, out), a, out
+  X = dft!(g, [(s, 0) for s in a])
+  out = [force!(g, X[k]) for k in 1:r]
+  g.stmts, a, out
 end
-
-# Drop statements that no needed output depends on.
-function prune_statements!(stmts, out)
-  live = Set{Symbol}(o for o in out if o isa Symbol)
-  keep = falses(length(stmts))
-  for i in length(stmts):-1:1
-    lhs, rhs = stmts[i].args
-    lhs in live || continue
-    keep[i] = true
-    collect_symbols!(live, rhs)
-  end
-  stmts[keep]
-end
-
-collect_symbols!(set, s::Symbol) = push!(set, s)
-collect_symbols!(set, ex::Expr) = foreach(x -> collect_symbols!(set, x), ex.args)
-collect_symbols!(set, x) = set

src/kernel_methods.jl

@@ -1,56 +1,53 @@
 # Compiled-ahead methods of the stage kernels for the fixed codelet set (radix × sign × vector
-# type, empty zero-input mask or full output mask) that `RChain` dispatches to. The `@generated`
-# kernels in kernels.jl handle every other combination (pruned stages); a package image keeps
-# the complete optimised IR of every 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 LAZY_RADICES are not in the set.
+# type) that `RChain` dispatches to, and of the chains of the small lengths (rchain.jl). 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
+# LAZY_RADICES are not in the set.
 
 kernel_subst(x, d) = x isa Symbol ? get(d, x, x) :
                      x isa Expr ? Expr(x.head, map(a -> kernel_subst(a, d), x.args)...) : x
 
 const KERNEL_HEADS = (
   (:stage!, :(stage!(py::Ptr{T}, px::Ptr{T}, pt::Ptr{T}, m::Int, s::Int, plo::Int, phi::Int,
-                     ::Val{R}, ::Val{σ}, ::Type{VT}, ::Val{mask})), :zin, :all),
+                     ::Val{R}, ::Val{σ}, ::Type{VT})), :all),
   (:stage_range!, :(stage_range!(py::Ptr{T}, px::Ptr{T}, pt::Ptr{T}, m::Int, s::Int, qlo::Int, qhi::Int,
-                                 plo::Int, phi::Int, ::Val{R}, ::Val{σ}, ::Type{VT}, ::Val{mask})), :zin, :all),
+                                 plo::Int, phi::Int, ::Val{R}, ::Val{σ}, ::Type{VT})), :all),
   (:stage_last!, :(stage_last!(py::Ptr{T}, px::Ptr{T}, s::Int, qlo::Int, qhi::Int,
-                               ::Val{R}, ::Val{σ}, ::Type{VT}, ::Val{mask})), :need, :all),
+                               ::Val{R}, ::Val{σ}, ::Type{VT})), :all),
   (:stage_lanes!, :(stage_lanes!(py::Ptr{T}, px::Ptr{T}, pt::Ptr{T}, m::Int, s0::Int, Sx::Int, Sy::Int,
-                                 B::Int, plo::Int, phi::Int, ::Val{R}, ::Val{σ}, ::Type{VT}, ::Val{mask})), :zin, :all),
+                                 B::Int, plo::Int, phi::Int, ::Val{R}, ::Val{σ}, ::Type{VT})), :all),
   (:stage_last_lanes!, :(stage_last_lanes!(py::Ptr{T}, px::Ptr{T}, s0::Int, Sx::Int, Sy::Int, B::Int,
-                                           qlo::Int, qhi::Int, ::Val{R}, ::Val{σ}, ::Type{VT}, ::Val{mask})), :need, :all),
+                                           qlo::Int, qhi::Int, ::Val{R}, ::Val{σ}, ::Type{VT})), :all),
   # the grouped first stage exists for 2 or 4 complex values per vector only; its scalar tail
   # runs on one complex value
   (:stage_first_group_range!, :(stage_first_group_range!(py::Ptr{T}, px::Ptr{T}, pt::Ptr{T}, m::Int,
-                                plo::Int, phi::Int, ::Val{R}, ::Val{σ}, ::Type{VT}, ::Val{mask})), :zin, :group),
+                                plo::Int, phi::Int, ::Val{R}, ::Val{σ}, ::Type{VT})), :group),
   (:stage_first_group_1!, :(stage_first_group_1!(py::Ptr{T}, px::Ptr{T}, pt::Ptr{T}, m::Int,
-                            plo::Int, phi::Int, ::Val{R}, ::Val{σ}, ::Type{VT}, ::Val{mask})), :zin, :tail))
+                            plo::Int, phi::Int, ::Val{R}, ::Val{σ}, ::Type{VT})), :tail))
 
-for T in (Float32, Float64), (name, head, kind, vts) in KERNEL_HEADS
+for T in (Float32, Float64), (name, head, vts) in KERNEL_HEADS
   body = getfield(@__MODULE__, Symbol(string(name)[1:end-1], "_body"))
   vts === :all ? (types = vec_ladder(T)) :
   vts === :group ? (types = filter(V -> ncomplex(V) in (2, 4), vec_ladder(T))) :
   (types = any(V -> ncomplex(V) in (2, 4), vec_ladder(T)) ? (Vec{2,T},) : ())
   for σ in (-1, 1), VT in types, R in RADICES
     R in LAZY_RADICES && continue
-    mask = kind === :need ? (1 << R) - 1 : 0
-    h = kernel_subst(head, Dict{Symbol,Any}(:T => T, :R => R, :σ => σ, :VT => VT, :mask => mask))
+    h = kernel_subst(head, Dict{Symbol,Any}(:T => T, :R => R, :σ => σ, :VT => VT))
     # the codelet statements name the type parameters T and VT, which a plain method does not bind
-    b = kernel_subst(body(T, R, σ, VT, mask), Dict{Symbol,Any}(:T => T, :VT => VT))
+    b = kernel_subst(body(T, R, σ, VT), 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, and every pruned call, uses the
-# @generated methods.
+# 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)}
-  Pr = typeof(NOPRUNE)
-  @eval $(Expr(:function, :(run!(py::Ptr{$T}, px::Ptr{$T}, pw1::Ptr{$T}, pw2::Ptr{$T}, P::$C, S::Int, Pr::$Pr)),
-               run_body(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, Pr::$Pr)),
-               run_lanes_body(T, Rs, σ, length(Rs), (), ())))
+                                       Sx::Int, Sy::Int, B::Int)),
+               run_lanes_body(T, Rs, σ, length(Rs))))
 end