benchmark: compare3.jl — FFTW vs registry FFTA vs checkouts, side by side

.gitignore

@@ -2,3 +2,5 @@
 docs/build
 benchmark/Manifest.toml
 benchmark/*/Manifest.toml
+benchmark/envs/
+benchmark/compare3_results/

benchmark/KERNEL_STAGES.md

@@ -0,0 +1,108 @@
+# Kernel stage breakdown (ComplexF64, single thread, neoverse-n1, Julia 1.12.6)
+
+Stage times are minimum times of the kernel re-run with the other stages removed; `copy` is one `copyto!` of the array (the memory-traffic floor of a single pass).
+
+Branch `integration/all` @ 901692b (before the SIMD-pass / leaves-first / codelet branches). Each block re-runs the kernel with stages removed, so the stage rows add up to the total; the last column relates each stage to one `copyto!` of the array (the memory-traffic floor of a single pass). Hardware counters were not available on this host.
+### n = 1024 = 2^10 (radix-4 recursion, 2 butterfly levels, 16 leaves of 64 points)
+
+| stage | time | share of FFTA total | time of one `copyto!` of the array: 346 ns |
+|:--|--:|--:|--:|
+| FFTA total | 9.52 µs | 100% | 27.5× copy |
+| leaves: 16 × 64-point codelet, input stride 16 | 5.32 µs | 56% | 15.4× copy |
+| (same codelets on contiguous input) | 5.28 µs | 55% | strided-load penalty 1.01× |
+| butterfly pass, blocks of 1024 (stride 256 between the 4 legs) | 2.11 µs | 22% | 6.1× copy; memory-only version of the pass 620 ns |
+| butterfly pass, blocks of 256 (stride 64 between the 4 legs) | 2.14 µs | 22% | 6.2× copy; memory-only version of the pass 631 ns |
+| leaves + passes (sum of the rows) | 9.56 µs | 100% | recursion/dispatch overhead = total − sum = -43 ns |
+| FFTW ESTIMATE / MEASURE | 6.32 µs / 5.28 µs | 66% / 55% | 18.2× / 15.2× copy |
+
+### n = 16384 = 2^14 (radix-4 recursion, 4 butterfly levels, 256 leaves of 64 points)
+
+| stage | time | share of FFTA total | time of one `copyto!` of the array: 6.04 µs |
+|:--|--:|--:|--:|
+| FFTA total | 244.68 µs | 100% | 40.5× copy |
+| leaves: 256 × 64-point codelet, input stride 256 | 97.64 µs | 40% | 16.2× copy |
+| (same codelets on contiguous input) | 86.64 µs | 35% | strided-load penalty 1.13× |
+| butterfly pass, blocks of 16384 (stride 4096 between the 4 legs) | 34.16 µs | 14% | 5.7× copy; memory-only version of the pass 14.08 µs |
+| butterfly pass, blocks of 4096 (stride 1024 between the 4 legs) | 34.40 µs | 14% | 5.7× copy; memory-only version of the pass 14.14 µs |
+| butterfly pass, blocks of 1024 (stride 256 between the 4 legs) | 33.84 µs | 14% | 5.6× copy; memory-only version of the pass 12.52 µs |
+| butterfly pass, blocks of 256 (stride 64 between the 4 legs) | 34.72 µs | 14% | 5.7× copy; memory-only version of the pass 12.86 µs |
+| leaves + passes (sum of the rows) | 234.76 µs | 96% | recursion/dispatch overhead = total − sum = 9.92 µs |
+| FFTW ESTIMATE / MEASURE | 176.64 µs / 157.64 µs | 72% / 64% | 29.2× / 26.1× copy |
+
+### n = 262144 = 2^18 (radix-4 recursion, 6 butterfly levels, 4096 leaves of 64 points)
+
+| stage | time | share of FFTA total | time of one `copyto!` of the array: 334.44 µs |
+|:--|--:|--:|--:|
+| FFTA total | 6.733 ms | 100% | 20.1× copy |
+| leaves: 4096 × 64-point codelet, input stride 4096 | 3.198 ms | 47% | 9.6× copy |
+| (same codelets on contiguous input) | 1.465 ms | 22% | strided-load penalty 2.18× |
+| butterfly pass, blocks of 262144 (stride 65536 between the 4 legs) | 560.61 µs | 8% | 1.7× copy; memory-only version of the pass 241.20 µs |
+| butterfly pass, blocks of 65536 (stride 16384 between the 4 legs) | 557.29 µs | 8% | 1.7× copy; memory-only version of the pass 237.44 µs |
+| butterfly pass, blocks of 16384 (stride 4096 between the 4 legs) | 553.53 µs | 8% | 1.7× copy; memory-only version of the pass 232.48 µs |
+| butterfly pass, blocks of 4096 (stride 1024 between the 4 legs) | 556.89 µs | 8% | 1.7× copy; memory-only version of the pass 233.80 µs |
+| butterfly pass, blocks of 1024 (stride 256 between the 4 legs) | 561.65 µs | 8% | 1.7× copy; memory-only version of the pass 219.52 µs |
+| butterfly pass, blocks of 256 (stride 64 between the 4 legs) | 600.41 µs | 9% | 1.8× copy; memory-only version of the pass 245.24 µs |
+| leaves + passes (sum of the rows) | 6.588 ms | 98% | recursion/dispatch overhead = total − sum = 145.52 µs |
+| FFTW ESTIMATE / MEASURE | 8.897 ms / 4.631 ms | 132% / 69% | 26.6× / 13.8× copy |
+
+### n = 4194304 = 2^22 (radix-4 recursion, 8 butterfly levels, 65536 leaves of 64 points)
+
+| stage | time | share of FFTA total | time of one `copyto!` of the array: 6.120 ms |
+|:--|--:|--:|--:|
+| FFTA total | 210.500 ms | 100% | 34.4× copy |
+| leaves: 65536 × 64-point codelet, input stride 65536 | 124.993 ms | 59% | 20.4× copy |
+| (same codelets on contiguous input) | 33.230 ms | 16% | strided-load penalty 3.76× |
+| butterfly pass, blocks of 4194304 (stride 1048576 between the 4 legs) | 9.916 ms | 5% | 1.6× copy; memory-only version of the pass 7.003 ms |
+| butterfly pass, blocks of 1048576 (stride 262144 between the 4 legs) | 9.250 ms | 4% | 1.5× copy; memory-only version of the pass 4.768 ms |
+| butterfly pass, blocks of 262144 (stride 65536 between the 4 legs) | 9.076 ms | 4% | 1.5× copy; memory-only version of the pass 4.336 ms |
+| butterfly pass, blocks of 65536 (stride 16384 between the 4 legs) | 9.051 ms | 4% | 1.5× copy; memory-only version of the pass 4.255 ms |
+| butterfly pass, blocks of 16384 (stride 4096 between the 4 legs) | 9.064 ms | 4% | 1.5× copy; memory-only version of the pass 4.171 ms |
+| butterfly pass, blocks of 4096 (stride 1024 between the 4 legs) | 9.362 ms | 4% | 1.5× copy; memory-only version of the pass 4.570 ms |
+| butterfly pass, blocks of 1024 (stride 256 between the 4 legs) | 10.913 ms | 5% | 1.8× copy; memory-only version of the pass 5.828 ms |
+| butterfly pass, blocks of 256 (stride 64 between the 4 legs) | 11.518 ms | 5% | 1.9× copy; memory-only version of the pass 6.324 ms |
+| leaves + passes (sum of the rows) | 203.144 ms | 97% | recursion/dispatch overhead = total − sum = 7.357 ms |
+| FFTW ESTIMATE / MEASURE | 242.198 ms / 104.928 ms | 115% / 50% | 39.6× / 17.1× copy |
+### n = 1000 = 8 × 125 (composite step: right = 125 COMPOSITE_FFT, left = 8 POW2RADIX4_FFT)
+
+| stage | time | share of FFTA total | time of one `copyto!` of the array: 338 ns |
+|:--|--:|--:|--:|
+| FFTA total | 35.32 µs | 100% | 104.4× copy |
+| 8 right sub-transforms of 125 (input stride 8) | 30.80 µs | 87% | 91.0× copy |
+| twiddle multiply pass | 1.27 µs | 4% | 3.7× copy |
+| 125 left sub-transforms of 8 (output stride 125) | 3.09 µs | 9% | 9.1× copy |
+| sum of the rows | 35.15 µs | 100% | overhead = 169 ns |
+| FFTW ESTIMATE / MEASURE | 10.28 µs / 7.36 µs | 29% / 21% | 30.4× / 21.7× copy |
+
+### n = 46305 = 27 × 1715 (composite step: right = 1715 COMPOSITE_FFT, left = 27 POW3_FFT)
+
+| stage | time | share of FFTA total | time of one `copyto!` of the array: 27.56 µs |
+|:--|--:|--:|--:|
+| FFTA total | 3.326 ms | 100% | 120.7× copy |
+| 27 right sub-transforms of 1715 (input stride 27) | 2.521 ms | 76% | 91.5× copy |
+| twiddle multiply pass | 68.76 µs | 2% | 2.5× copy |
+| 1715 left sub-transforms of 27 (output stride 1715) | 693.69 µs | 21% | 25.2× copy |
+| sum of the rows | 3.284 ms | 99% | overhead = 42.32 µs |
+| FFTW ESTIMATE / MEASURE | 1.052 ms / 865.49 µs | 32% / 26% | 38.2× / 31.4× copy |
+
+### n = 1000000 = 64 × 15625 (composite step: right = 15625 COMPOSITE_FFT, left = 64 POW2RADIX4_FFT)
+
+| stage | time | share of FFTA total | time of one `copyto!` of the array: 1.274 ms |
+|:--|--:|--:|--:|
+| FFTA total | 100.257 ms | 100% | 78.7× copy |
+| 64 right sub-transforms of 15625 (input stride 64) | 83.168 ms | 83% | 65.3× copy |
+| twiddle multiply pass | 2.005 ms | 2% | 1.6× copy |
+| 15625 left sub-transforms of 64 (output stride 15625) | 11.154 ms | 11% | 8.8× copy |
+| sum of the rows | 96.327 ms | 96% | overhead = 3.930 ms |
+| FFTW ESTIMATE / MEASURE | 55.945 ms / 24.771 ms | 56% / 25% | 43.9× / 19.4× copy |
+
+### n = 65537 (Bluestein, padded to 262144)
+
+| stage | time | share of FFTA total | time of one `copyto!` of the array: 60.56 µs |
+|:--|--:|--:|--:|
+| FFTA total | 15.552 ms | 100% | 256.8× copy |
+| chirp multiply + zero pad | 191.36 µs | 1% | |
+| forward FFT of 262144 | 6.390 ms | 41% | |
+| pointwise product | 540.45 µs | 3% | |
+| second FFT of 262144 | 6.401 ms | 41% | |
+| output chirp multiply | 131.68 µs | 1% | |
+| FFTW ESTIMATE / MEASURE | 3.362 ms / 2.660 ms | 22% / 17% | |

benchmark/README.md

@@ -63,6 +63,26 @@ The committed [`REPORT.md`](REPORT.md) was produced on the machine described
 in its *Environment* section; re-run the suite to obtain numbers for your own
 hardware.
 
+## Several implementations side by side (`compare3.jl`)
+
+```bash
+cd benchmark
+julia --project=. compare3.jl                        # FFTW vs FFTA 0.3.1 (registry) vs this checkout
+julia --project=. compare3.jl --impl fftw=fftw --impl int=/path/to/base --impl new=/path/to/branch --ref int \
+                              --threads 1,4,16 --only nd,batched
+```
+
+Each `--impl NAME=SPEC` column runs in its own process and environment
+(`envs/NAME/`, git-ignored): `fftw`, `@0.3.1` (an FFTA version from the
+registry) or a path to an FFTA checkout / worktree. The same case list as
+`suite.jl` (`cases.jl`; the `--only/--kinds/--maxlog2/--sizes/--seconds`
+options apply) is run per thread count in `--threads`, and one markdown
+table (`compare3_results/COMPARE3.md`) shows the times, each column's ratio
+to the first column and its speedup over the `--ref` column, with per-class
+geometric means and, for several thread counts, the thread scaling.
+Results are checksummed against the first column. `--render-only`
+re-renders the table from the JSON files of a previous run.
+
 ## Original 1D suite (`run_benchmarks.jl`)
 
 ## Structure

benchmark/cases.jl

@@ -0,0 +1,134 @@
+# Size classes and the case list shared by `suite.jl` and `compare3.jl`.
+#
+# A case is a NamedTuple (kind, T, sz, dims, shape, nthreads, measure); `T` is
+# the *real* element type (complex cases use `Complex{T}`).
+
+using Primes
+
+function size_classes(maxlog2)
+    nmax = 1 << maxlog2
+    pow2   = [1 << k for k in 3:maxlog2]
+    # 2^a 3^b 5^c 7^d, not powers of two, spread log-uniformly
+    smooth_all = Int[]
+    for a in 0:maxlog2, b in 0:14, c in 0:9, d in 0:8
+        n = 2^a * 3^b * 5^c * 7^d
+        (n <= nmax && n >= 8 && !ispow2(n)) && push!(smooth_all, n)
+    end
+    sort!(unique!(smooth_all))
+    smooth = Int[]
+    for t in exp.(range(log(12), log(nmax), length = 18))
+        push!(smooth, smooth_all[argmin(abs.(log.(smooth_all) .- log(t)))])
+    end
+    # explicitly include some "classic" smooth sizes
+    for n in (12, 60, 120, 360, 720, 1000, 1000_000)
+        n <= nmax && push!(smooth, n)
+    end
+    sort!(unique!(smooth))
+    # primes: below (DFT path) and above (Bluestein path) the cutoff
+    prime = Int[]
+    for t in exp.(range(log(7), log(nmax), length = 16))
+        push!(prime, nextprime(round(Int, t)))
+    end
+    # straddle the DFT/Bluestein crossover: the O(n²) leaf is used below the
+    # cutoff (73 originally, 29 since the twiddle-table work), so sample the
+    # band 23–46 as well as 61–79
+    for n in (23, 29, 31, 37, 43, 61, 71, 73, 79)
+        push!(prime, n)
+    end
+    filter!(<=(nmax), prime); sort!(unique!(prime))
+    # awkward: (prime > cutoff) × small factor
+    awk = Int[]
+    for p in (101, 1009, 4099, 16411, 65537, 262147), f in (2, 3, 4, 6, 16)
+        n = p * f
+        n <= nmax && push!(awk, n)
+    end
+    sort!(unique!(awk))
+    return (pow2 = pow2, smooth = smooth, prime = prime, awkward = awk)
+end
+
+function class_of(classes, n)
+    for (k, v) in pairs(classes)
+        n in v && return String(k)
+    end
+    return "other"
+end
+
+"""
+    case_list(; only, kinds, maxlog2, sizes, nthreads, measure) -> Vector{NamedTuple}
+
+The sweep of `suite.jl`: 1D classes (restricted with `classes`) × {Float64, Float32} × {fft, rfft}, a
+FFTW-`MEASURE` column for 1D pow2 ComplexF64, 2D/3D, batched `dims`, and the
+threaded subset. `nthreads > 1` selects the threaded subset only when `only`
+contains `"threads"`; the other sections are always emitted with `nthreads`.
+"""
+function case_list(; only = ("1d", "nd", "batched", "threads"), kinds = (:fft, :rfft),
+                     maxlog2 = 22, sizes = Int[], nthreads = 1, measure = true,
+                     classes = ("pow2", "smooth", "prime", "awkward"))
+    sizes_by_class = size_classes(maxlog2)
+    nmax = 1 << maxlog2
+    cases = NamedTuple[]
+    add!(kind, T, sz, dims, shape; measure = false, thr = nthreads) =
+        kind in kinds && push!(cases, (kind = kind, T = T, sz = sz, dims = dims, shape = shape,
+                                        nthreads = thr, measure = measure,
+                                        class = length(sz) == 1 ? class_of(sizes_by_class, sz[1]) : "nd"))
+    if "1d" in only
+        for cls in (:pow2, :smooth, :prime, :awkward), n in sizes_by_class[cls]
+            String(cls) in classes || continue
+            isempty(sizes) || n in sizes || continue
+            for T in (Float64, Float32)
+                add!(:fft,  T, (n,), 1, "1d")
+                add!(:rfft, T, (n,), 1, "1d")
+            end
+        end
+        if measure
+            for n in sizes_by_class.pow2
+                isempty(sizes) || n in sizes || continue
+                add!(:fft, Float64, (n,), 1, "1d"; measure = true)
+            end
+        end
+    end
+    if "nd" in only
+        nd2 = filter(n -> n * n <= nmax, [8, 16, 32, 64, 128, 256, 512, 1024, 2048])
+        for n in nd2, T in (Float64, Float32)
+            add!(:fft,  T, (n, n), (1, 2), "2d")
+            add!(:rfft, T, (n, n), (1, 2), "2d")
+        end
+        for sz in ((1000, 1000), (720, 480), (1009, 64), (64, 1009), (127, 257))
+            prod(sz) <= nmax || continue
+            add!(:fft,  Float64, sz, (1, 2), "2d")
+            add!(:rfft, Float64, sz, (1, 2), "2d")
+        end
+        nd3 = filter(n -> n^3 <= nmax, [8, 16, 32, 64, 128])
+        for n in nd3, T in (Float64, Float32)
+            add!(:fft,  T, (n, n, n), (1, 2, 3), "3d")
+            add!(:rfft, T, (n, n, n), (1, 2, 3), "3d")
+        end
+    end
+    if "batched" in only
+        for n in filter(n -> n * 64 <= nmax, [64, 256, 1024, 4096, 16384, 65536]), T in (Float64, Float32)
+            add!(:fft,  T, (n, 64), 1, "batched_dim1")
+            add!(:fft,  T, (64, n), 2, "batched_dim2")
+            add!(:rfft, T, (n, 64), 1, "batched_dim1")
+            add!(:rfft, T, (64, n), 2, "batched_dim2")
+        end
+        for (n, m) in ((256, 4096), (1024, 1024), (4096, 256))
+            n * m <= nmax || continue
+            add!(:rfft, Float64, (n, m), 1, "batched_dim1")
+        end
+    end
+    if "threads" in only && nthreads > 1
+        for n in filter(n -> n >= 1 << 16, sizes_by_class.pow2)
+            add!(:fft, Float64, (n,), 1, "1d")
+        end
+        for n in filter(n -> n * n <= nmax && n >= 256, [256, 512, 1024, 2048])
+            add!(:fft, Float64, (n, n), (1, 2), "2d")
+        end
+        for n in filter(n -> n * 64 <= nmax && n >= 1024, [1024, 4096, 16384, 65536])
+            add!(:fft, Float64, (n, 64), 1, "batched_dim1")
+        end
+    end
+    return unique(cases)
+end
+
+casekey(c) = string(c.kind, " ", c.T, " ", join(c.sz, "×"), " dims=", join(c.dims, ","),
+                    c.nthreads > 1 ? " ($(c.nthreads) thr)" : "", c.measure ? " MEASURE" : "")