Add permutedims! benchmark suite

benchmark/Project.toml

@@ -0,0 +1,10 @@
+[deps]
+ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+Strided = "5e0ebb24-38b0-5f93-81fe-25c709ecae67"
+TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[compat]
+ArgParse = "1"
+BenchmarkTools = "1"
+TOML = "1"

benchmark/README.md

@@ -0,0 +1,45 @@
+# Strided.jl benchmarks
+
+This directory contains a [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl) suite (`benchmarks.jl`, defining `SUITE`) for the `permutedims!` machinery, in the format expected by [AirSpeedVelocity.jl](https://github.com/MilesCranmer/AirSpeedVelocity.jl) and PkgBenchmark.jl.
+
+Each case fixes a *shape* (dimension ratios) and a *permutation* and is swept over a list of total lengths, so its time can be divided by that of a plain memory copy of the same length — the auto-generated `copy` group — to get a machine-independent efficiency (`≥ 1`, approaching `1` when bandwidth-bound).
+The default groups (`balanced`, `skewed`, `strided`) live in `cases.toml`; the `strided` group makes the input/output non-contiguous strided views to exercise the paths a dense array never hits.
+See the comments in `cases.toml` for the format; a different cases file can be passed with `--cases`.
+
+The copy baseline is single-threaded, so multi-threaded permutation ratios are relative to a single-threaded memory copy.
+
+## Running manually
+
+`benchmarks.jl` only defines `SUITE`, but its contents are configurable through command-line arguments, which can be passed after a `--` separator:
+
+```bash
+julia --project=benchmark benchmark/benchmarks.jl --help
+julia --project=benchmark -t 8 \
+    -e 'include("benchmark/benchmarks.jl"); display(run(SUITE; verbose=true))' \
+    -- -g balanced -T Float64 -f L=1048576
+```
+
+The first use requires instantiating the environment:
+
+```bash
+julia --project=benchmark -e 'using Pkg; Pkg.develop(path="."); Pkg.instantiate()'
+```
+
+Multithreaded variants are generated only when Julia is started with more than one thread (`-t` / `JULIA_NUM_THREADS`).
+The full unfiltered grid takes on the order of an hour; restrict it with `-g`/`-T`/`-n`/`-f` (or `benchpkg --filter`).
+
+## Comparing revisions with AirSpeedVelocity
+
+```bash
+julia -e 'using Pkg; Pkg.add("AirSpeedVelocity")'  # once, in the global env
+
+export JULIA_NUM_THREADS=8
+benchpkg Strided --rev=v2.6.1,main,dirty --path=. --bench-on=main --filter=T=Float64
+```
+
+from the repository root.
+`benchpkgtable` / `benchpkgplot` format the results.
+AirSpeedVelocity includes the script with empty `ARGS`, so the default (full) configuration applies; restrict it with `--filter`, which matches substrings of the benchmark names (element type, thread count, permutation, and size are all part of the name).
+
+Note that `--bench-on` requires a revision that already contains this suite; to benchmark with an uncommitted version of the script, pass it explicitly along with its non-standard dependencies, e.g. `--script=benchmark/benchmarks.jl --add=ArgParse`.
+The `--output-dir` must exist beforehand, and the current AirSpeedVelocity version (0.6.5) crashes on `--bench-on=dirty`, so prefer a committed revision there.

benchmark/benchmarks.jl

@@ -0,0 +1,234 @@
+# Benchmark suite for Strided.jl, in the standard BenchmarkTools.jl format
+# expected by AirSpeedVelocity.jl / PkgBenchmark.jl: this file defines a
+# `SUITE::BenchmarkGroup` which the harness loads and runs for each revision.
+# This file never runs the suite itself.
+#
+# The suite is configured through command-line arguments (see --help). When
+# included by a harness (empty ARGS), the default configuration is used: the
+# full suite, which takes on the order of an hour per revision. Use the
+# harness' own filtering (e.g. `benchpkg --filter`) to restrict it. For a
+# manual run, arguments can be passed after `--`:
+#
+#   julia --project=benchmark -t 8 \
+#       -e 'include("benchmark/benchmarks.jl"); display(run(SUITE; verbose=true))' \
+#       -- --groups=balanced --eltypes=Float64
+#
+# The suite focuses on `permutedims!` of strided arrays. Each case fixes a
+# *shape* (dimension ratios) and a *permutation*, and is swept over a list of
+# total lengths; the cases are defined in a TOML file, by default the cases.toml
+# next to this script (see its comments for the format and the groups). A
+# different file can be passed with --cases.
+#
+# Alongside `permutedims!`, the suite auto-generates a "copy" group: a plain
+# `copyto!` on `Vector`s of each swept length. Since a permutation and a copy
+# both read and write the same number of elements, dividing a permutation's
+# time by the copy time of the same length gives a machine-independent
+# efficiency (>= 1, approaching 1 when bandwidth-bound). The copy baseline is
+# single-threaded, so multi-threaded permutation ratios are relative to a
+# single-threaded memory copy.
+#
+# Structure of the suite:
+#   SUITE["permutedims!"][group]["T=$T"]["nthreads=$nt"][case]
+#   SUITE["copy"]["T=$T"]["L=$L"]
+
+using ArgParse
+using BenchmarkTools
+using Strided
+using Strided: StridedView
+using TOML
+
+function parse_config(args)
+    s = ArgParseSettings(;
+        prog = "benchmark/benchmarks.jl",
+        description = "Benchmark suite for permutedims! on strided arrays; " *
+            "defines `SUITE` without running it.",
+    )
+    #! format: off
+    @add_arg_table! s begin
+        "--eltypes", "-T"
+            help = "comma-separated element types to benchmark"
+            default = "Float64,ComplexF64"
+        "--cases", "-c"
+            help = "TOML file defining the benchmark cases"
+            default = joinpath(@__DIR__, "cases.toml")
+        "--groups", "-g"
+            help = "comma-separated case groups to benchmark; " *
+                   "defaults to all groups in the cases file"
+            default = ""
+        "--nthreads", "-n"
+            help = "comma-separated Strided thread counts; values above the " *
+                   "number of Julia threads (set with julia -t) are clamped"
+            default = "1,$(Threads.nthreads())"
+        "--filter", "-f"
+            help = "only include benchmarks whose full name contains this substring"
+            default = ""
+    end
+    #! format: on
+    return parse_args(args, s)
+end
+
+const CONFIG = parse_config(ARGS)
+
+const ELTYPES = map(split(CONFIG["eltypes"], ',')) do s
+    T = getfield(Base, Symbol(strip(s)))
+    T isa Type || error("--eltypes: `$s` is not a type")
+    return T
+end
+const NTHREADS = sort!(
+    unique(
+        clamp.(
+            parse.(Int, strip.(split(CONFIG["nthreads"], ','))), 1, Threads.nthreads()
+        ),
+    ),
+)
+
+# Arrays shorter than this are never threaded by Strided's kernel, so threaded
+# variants of such cases would just duplicate the single-threaded numbers.
+const MINTHREADLENGTH = isdefined(Strided, :MINTHREADLENGTH) ? Strided.MINTHREADLENGTH : 1 << 15
+
+const CASEGROUPS = TOML.parsefile(CONFIG["cases"])
+const GROUPS = if isempty(strip(CONFIG["groups"]))
+    sort!(collect(keys(CASEGROUPS)))
+else
+    strip.(split(CONFIG["groups"], ','))
+end
+
+# Column-major strides of a dense array of size `sz`.
+colstrides(sz::NTuple{N, Int}) where {N} = ntuple(i -> prod(ntuple(j -> sz[j], i - 1)), N)
+
+# Dimensions realizing (approximately) `total` elements at the given shape ratio.
+function shapedims(shape::NTuple{N, <:Real}, total::Int) where {N}
+    c = (total / prod(shape))^(1 / N)
+    return ntuple(i -> max(1, round(Int, shape[i] * c)), N)
+end
+
+# A StridedView of logical size `sz` whose axes are spaced by `mult` inside a
+# freshly allocated (page-faulted) backing buffer. `mult` all ones gives a dense
+# view (the unit-stride fast path); a leading value > 1 gives a non-unit
+# innermost stride, a later value > 1 gives gaps between higher dimensions.
+function make_view(::Type{T}, sz::NTuple{N, Int}, mult::NTuple{N, Int}) where {T, N}
+    all(isone, mult) && return StridedView(rand(T, sz))
+    bufsz = ntuple(i -> sz[i] * mult[i], N)
+    strides = ntuple(i -> mult[i] * colstrides(bufsz)[i], N)
+    return StridedView(rand(T, prod(bufsz)), sz, strides)
+end
+
+# Expand a group spec into concrete case instances (one per swept total).
+function expand(spec, group)
+    totals = haskey(spec, "totals") ? Int.(spec["totals"]) : Int[]
+    insts = NamedTuple[]
+    for c in spec["cases"]
+        p = (Int.(c["p"])...,)
+        isperm(p) || error("$group: `p = $(c["p"])` is not a permutation")
+        N = length(p)
+        min = haskey(c, "stride_in") ? (Int.(c["stride_in"])...,) : ntuple(one, N)
+        mout = haskey(c, "stride_out") ? (Int.(c["stride_out"])...,) : ntuple(one, N)
+        (length(min) == N == length(mout)) ||
+            error("$group: `stride_in`/`stride_out` must have length $N")
+        if haskey(c, "size")
+            dims = (Int.(c["size"])...,)
+            push!(insts, (; p, dims, min, mout, L = prod(dims), shape = nothing))
+        elseif haskey(c, "shape")
+            shape = (Float64.(c["shape"])...,)
+            length(shape) == N || error("$group: `shape` must have length $N")
+            isempty(totals) &&
+                error("$group: case uses `shape` but the group has no `totals`")
+            for total in totals
+                push!(insts, (; p, dims = shapedims(shape, total), min, mout, L = total, shape))
+            end
+        else
+            error("$group: each case needs a `size` or a `shape`")
+        end
+    end
+    return insts
+end
+
+function casename(inst)
+    s = "p=" * join(inst.p) * "_L=" * string(inst.L)
+    inst.shape === nothing || (s *= "_shape=" * join(round.(Int, inst.shape), 'x'))
+    any(!isone, inst.min) && (s *= "_sin=" * join(inst.min, 'x'))
+    any(!isone, inst.mout) && (s *= "_sout=" * join(inst.mout, 'x'))
+    return s
+end
+
+# Arrays are allocated (and touched) in the per-sample setup rather than at
+# suite-construction time: the full grid would otherwise keep every array alive
+# at once. `evals` is fixed so that tuning never re-runs the setup, and is > 1
+# only for cases too short to time reliably in a single evaluation.
+function addcase!(group::BenchmarkGroup, inst, ::Type{T}, nt; samples, seconds) where {T}
+    L = prod(inst.dims)
+    nt > 1 && L < MINTHREADLENGTH && return group
+    dsz = ntuple(i -> inst.dims[inst.p[i]], length(inst.p))
+    evals = max(1, (1 << 15) ÷ L)
+    p, sz, min, mout = inst.p, inst.dims, inst.min, inst.mout
+    group[casename(inst)] = @benchmarkable(
+        permutedims!(dst, src, $p),
+        setup = (
+            Strided.set_num_threads($nt);
+            src = make_view($T, $sz, $min);
+            dst = make_view($T, $dsz, $mout)
+        ),
+        evals = evals,
+        samples = samples,
+        seconds = seconds,
+    )
+    return group
+end
+
+function addcopy!(group::BenchmarkGroup, L, ::Type{T}) where {T}
+    evals = max(1, (1 << 15) ÷ L)
+    group["L=$L"] = @benchmarkable(
+        copyto!(dst, src),
+        setup = (src = rand($T, $L); dst = Vector{$T}(undef, $L)),
+        evals = evals,
+        samples = 100,
+        seconds = 5,
+    )
+    return group
+end
+
+function filtered(group::BenchmarkGroup, pattern::AbstractString, prefix::String = "")
+    out = BenchmarkGroup(group.tags)
+    for (key, value) in group
+        name = prefix * "/" * key
+        if value isa BenchmarkGroup
+            sub = filtered(value, pattern, name)
+            isempty(sub.data) || (out[key] = sub)
+        elseif occursin(pattern, name)
+            out[key] = value
+        end
+    end
+    return out
+end
+
+const SUITE = let
+    suite = BenchmarkGroup()
+    permute = addgroup!(suite, "permutedims!")
+    lengths = Set{Int}()
+    for name in GROUPS
+        haskey(CASEGROUPS, name) ||
+            error("--groups: group `$name` not found in $(CONFIG["cases"])")
+        spec = CASEGROUPS[name]
+        insts = expand(spec, name)
+        samples, seconds = spec["samples"], spec["seconds"]
+        g = addgroup!(permute, name)
+        for T in ELTYPES
+            tg = addgroup!(g, "T=$T")
+            for nt in NTHREADS
+                ng = addgroup!(tg, "nthreads=$nt")
+                for inst in insts
+                    addcase!(ng, inst, T, nt; samples, seconds)
+                    push!(lengths, inst.L)
+                end
+            end
+        end
+    end
+    copies = addgroup!(suite, "copy")
+    for T in ELTYPES
+        tg = addgroup!(copies, "T=$T")
+        for L in sort!(collect(lengths))
+            addcopy!(tg, L, T)
+        end
+    end
+    isempty(CONFIG["filter"]) ? suite : filtered(suite, CONFIG["filter"])
+end

benchmark/cases.toml

@@ -0,0 +1,77 @@
+# Benchmark cases for benchmarks.jl.
+#
+# Each case is a permutation applied to an array whose *shape* (the ratio of its
+# dimensions) is fixed while its *total length* is swept, so a case can be
+# compared against a plain memory copy of the same length (the auto-generated
+# "copy" group) to get a machine-independent efficiency. Dimensions are derived
+# from `shape` and each `total` as d_i = round(s_i * (total / prod(s))^(1/N)),
+# so the realized length is approximately (not exactly) `total`.
+#
+# [group] fields:
+#   samples, seconds   BenchmarkTools budget (per case; `seconds` caps the big
+#                      cases, evals stabilize the tiny ones)
+#   totals             list of target element counts to sweep
+#   cases              list of cases, each a table with:
+#       p              1-based permutation (output axis i comes from input p[i])
+#       shape          relative dimension ratios, length N == length(p)
+#                      (alternatively `size` for explicit dims, ignoring totals)
+#       stride_in      optional per-input-axis extent multiplier (default all 1)
+#       stride_out     optional per-output-axis extent multiplier (default all 1)
+#                      A multiplier > 1 makes the array a non-contiguous strided
+#                      view: a leading >1 gives a non-unit innermost stride, a
+#                      later >1 gives gaps between higher dimensions.
+
+# Balanced shapes (all dimensions equal) across ranks 2-6, plus a couple of
+# non-reversing permutations. The bandwidth-bound baseline for the transpose.
+[balanced]
+samples = 100
+seconds = 5
+totals = [1024, 32768, 1048576, 8388608, 67108864]
+cases = [
+    { p = [2, 1],             shape = [1, 1] },
+    { p = [3, 2, 1],          shape = [1, 1, 1] },
+    { p = [2, 3, 1],          shape = [1, 1, 1] },
+    { p = [4, 3, 2, 1],       shape = [1, 1, 1, 1] },
+    { p = [3, 4, 1, 2],       shape = [1, 1, 1, 1] },
+    { p = [5, 4, 3, 2, 1],    shape = [1, 1, 1, 1, 1] },
+    { p = [6, 5, 4, 3, 2, 1], shape = [1, 1, 1, 1, 1, 1] },
+]
+
+# Skewed shapes: one dimension much larger than the rest, either leading (fat
+# first axis) or trailing (fat last axis), as in the HPTT/TTC reference set.
+[skewed]
+samples = 100
+seconds = 5
+totals = [1024, 32768, 1048576, 8388608, 67108864]
+cases = [
+    { p = [2, 1],             shape = [32, 1] },
+    { p = [2, 1],             shape = [1, 32] },
+    { p = [3, 2, 1],          shape = [16, 1, 1] },
+    { p = [3, 2, 1],          shape = [1, 1, 16] },
+    { p = [4, 3, 2, 1],       shape = [8, 1, 1, 1] },
+    { p = [4, 3, 2, 1],       shape = [1, 1, 1, 8] },
+    { p = [6, 5, 4, 3, 2, 1], shape = [8, 1, 1, 1, 1, 1] },
+    { p = [6, 5, 4, 3, 2, 1], shape = [1, 1, 1, 1, 1, 8] },
+]
+
+# Truly strided views: the input or output is a non-contiguous slice of a larger
+# buffer. `stride_in/out = [2, 1, ...]` gives a non-unit innermost stride (the
+# strongly non-contiguous case that defeats the unit-stride fast path);
+# `[1, 2, ...]` keeps the innermost dimension contiguous but adds a gap between
+# successive second-axis slices.
+[strided]
+samples = 100
+seconds = 5
+totals = [1024, 32768, 1048576, 8388608, 67108864]
+cases = [
+    { p = [2, 1],       shape = [1, 1],       stride_in  = [2, 1] },
+    { p = [2, 1],       shape = [1, 1],       stride_in  = [1, 2] },
+    { p = [2, 1],       shape = [1, 1],       stride_out = [2, 1] },
+    { p = [2, 1],       shape = [1, 1],       stride_out = [1, 2] },
+    { p = [3, 2, 1],    shape = [1, 1, 1],    stride_in  = [2, 1, 1] },
+    { p = [3, 2, 1],    shape = [1, 1, 1],    stride_in  = [1, 2, 1] },
+    { p = [3, 2, 1],    shape = [1, 1, 1],    stride_out = [2, 1, 1] },
+    { p = [4, 3, 2, 1], shape = [1, 1, 1, 1], stride_in  = [2, 1, 1, 1] },
+    { p = [4, 3, 2, 1], shape = [1, 1, 1, 1], stride_in  = [1, 2, 1, 1] },
+    { p = [4, 3, 2, 1], shape = [1, 1, 1, 1], stride_out = [2, 1, 1, 1] },
+]