diff --git a/benchmark/Project.toml b/benchmark/Project.toml new file mode 100644 index 0000000..c75c69c --- /dev/null +++ b/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" diff --git a/benchmark/README.md b/benchmark/README.md new file mode 100644 index 0000000..0545750 --- /dev/null +++ b/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. diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl new file mode 100644 index 0000000..7398ad6 --- /dev/null +++ b/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 diff --git a/benchmark/cases.toml b/benchmark/cases.toml new file mode 100644 index 0000000..c290890 --- /dev/null +++ b/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] }, +]