Benchmark harness, figures and panel data in bench/

bench/.gitignore

@@ -0,0 +1,5 @@
+machines.local/
+vkfft/vkfft_bench
+cufft/cufft_bench
+ducc-ffi/libducc_ffi.*
+ducc-ffi/build.log

bench/Project.toml

@@ -0,0 +1,8 @@
+[deps]
+AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
+FFTL = "b53ef638-a7f0-46a4-97a9-16ee08a896af"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
+
+[sources]
+FFTL = {path = ".."}

bench/ab/bench_mt.jl

@@ -0,0 +1,18 @@
+# Multi-thread A/B timings of one package version: minimum over blocks of calls, threads = ARGS[1].
+# Pin to a fixed core set with taskset and start Julia with the same -t, one process per variant.
+# Usage: taskset -c 64-79 julia --project=<env> -t 16 ab/bench_mt.jl 16
+
+using FFTL, Printf
+nt = parse(Int, ARGS[1])
+function tm(f; reps = 10)
+  f(); minimum((t = time_ns(); for _ in 1:reps; f(); end; (time_ns() - t) / 1e9 / reps) for _ in 1:5)
+end
+cases = [("c2c", (37888,)), ("c2c", (64576,)), ("c2c", (92935,)), ("c2c", (1009, 16)), ("c2c", (101, 101)),
+         ("c2c", (1009, 1009)), ("c2c", (101, 101, 101)), ("c2c", (1024, 1024)), ("c2c", (2^20,)), ("c2c", (1000003,)),
+         ("r2c", (37888,)), ("r2c", (92935,)), ("r2c", (101, 101, 101)), ("r2c", (1024, 1024))]
+for (k, sh) in cases
+  x = k == "r2c" ? randn(Float64, sh...) : randn(ComplexF64, sh...)
+  p = k == "r2c" ? FFTL.plan_rfft(x; threads = nt) : FFTL.plan_fft(x; threads = nt)
+  t = tm(() -> p * x)
+  @printf("%s %-14s %10.1f us\n", k, string(sh), 1e6 * t)
+end

bench/ab/bench_pr.jl

@@ -0,0 +1,16 @@
+# Single-case A/B timings of one package version: one thread, minimum over blocks of calls.
+# Run one process per variant and thread count, pinned, and interleave the variants over rounds.
+# Usage: taskset -c N julia --project=<env with FFTL developed from the checkout> -t 1 ab/bench_pr.jl
+
+using FFTL, Printf
+function tm(f; reps = 30)
+  f(); minimum((t = time_ns(); for _ in 1:reps; f(); end; (time_ns() - t) / 1e9 / reps) for _ in 1:5)
+end
+cases = [("c2c", (1024,)), ("c2c", (4096,)), ("c2c", (1000,)), ("c2c", (101,)), ("c2c", (1009,)), ("c2c", (37888,)),
+         ("c2c", (64576,)), ("c2c", (1009, 1009)), ("c2c", (101, 101, 101)), ("r2c", (37888,)), ("r2c", (92935,)), ("c2c", (12000,))]
+for (k, sh) in cases
+  x = k == "r2c" ? randn(Float64, sh...) : randn(ComplexF64, sh...)
+  p = k == "r2c" ? FFTL.plan_rfft(x) : FFTL.plan_fft(x)
+  t = tm(() -> p * x)
+  @printf("%s %-14s %9.1f us\n", k, string(sh), 1e6 * t)
+end

bench/cufft/build.sh

@@ -0,0 +1,8 @@
+#!/bin/bash
+# Builds cufft_bench against the CUDA toolkit's cuFFT (CUDA_HOME, default /usr/local/cuda).
+#   bench/cufft/build.sh
+set -euo pipefail
+cd "$(dirname "$0")"
+CU=${CUDA_HOME:-/usr/local/cuda}
+g++ -O2 -std=c++17 -I"$CU/include" cufft_bench.cpp -o cufft_bench -L"$CU/lib64" -lcufft -lcudart -Wl,-rpath,"$CU/lib64"
+echo "built $(pwd)/cufft_bench"

bench/cufft/cufft.jl

@@ -0,0 +1,25 @@
+# Julia side of the native cuFFT harness (cufft_bench.cpp): runs one transform in the harness
+# process and returns its output and timings.
+#   include("bench/cufft/cufft.jl"); y, t = cufft_native(:r2c, x)   # t = (plan_ms, first_ms, exec_us)
+# Kinds: :c2c, :bc2c, :r2c, :c2r (needs n1 = size of the real array). ndims: transform along the
+# first ndims dimensions, batch over the others.
+const CUFFT_BIN = joinpath(@__DIR__, "cufft_bench")
+
+function cufft_native(kind::Symbol, x::AbstractArray; n1 = nothing, ndims = Base.ndims(x), reps = 50, samples = 7)
+  T = real(eltype(x))
+  sz = kind === :c2r ? (n1, size(x)[2:end]...) : size(x)
+  osz = kind === :r2c ? (size(x, 1) ÷ 2 + 1, size(x)[2:end]...) : kind === :c2r ? sz : size(x)
+  OT = kind in (:c2c, :bc2c, :r2c) ? Complex{T} : T
+  dir = mktempdir()
+  fin = joinpath(dir, "in.bin"); fout = joinpath(dir, "out.bin")
+  write(fin, collect(x))
+  prec = T === Float64 ? "f64" : "f32"
+  dims = sz[1:ndims]; b = prod(sz[ndims+1:end]; init = 1)
+  args = b > 1 ? [string.(dims)..., "b$b"] : string.(dims)
+  line = readchomp(`$CUFFT_BIN $kind $prec $fin $fout $reps $samples $args`)
+  y = Array{OT}(undef, osz)
+  read!(fout, y)
+  rm(dir; recursive = true)
+  v = parse.(Float64, split(line))
+  y, (plan_ms = v[1], first_ms = v[2], exec_us = v[3])
+end

bench/cufft/cufft_bench.cpp

@@ -0,0 +1,103 @@
+// cuFFT benchmark harness, the counterpart of bench/vkfft/vkfft_bench.cpp: one transform per
+// process, data through files, timings on stdout; the plan is made once and reused.
+//
+//   cufft_bench <kind> <f32|f64> <in.bin> <out.bin> <reps> <samples> n1 [n2 [n3]] [bB]
+//
+// kind: c2c (forward), bc2c (backward, unnormalised), r2c, c2r (unnormalised). Arrays are column
+// major (n1 contiguous, as in Julia), so cuFFT gets the dimensions in reverse order; r2c/c2r use
+// the unpadded layout (n1 real, n1/2+1 complex values per line). A last argument bB transforms B
+// arrays one after another (batch). All transforms are out of place; cuFFT's c2r may overwrite
+// its input for more than one dimension, so the input is reloaded before each sample (outside the
+// timing), as in the VkFFT harness.
+//
+// Output, one line: plan_ms first_ms exec_us; plan_ms is cufftPlanMany (including cuFFT's
+// kernel setup), first_ms the first execution including its synchronisation, exec_us the median
+// over <samples> of the mean time per transform of <reps> back-to-back executions with one
+// synchronisation.
+//
+// Build: bench/cufft/build.sh
+#include <algorithm>
+#include <chrono>
+#include <cstdio>
+#include <cstdlib>
+#include <string>
+#include <vector>
+#include <cuda_runtime.h>
+#include <cufft.h>
+
+using clk = std::chrono::steady_clock;
+static double ms_since(clk::time_point t) { return std::chrono::duration<double, std::milli>(clk::now() - t).count(); }
+static void die(const char* what, long code) {
+  fprintf(stderr, "cufft_bench: %s failed (%ld)\n", what, code);
+  exit(1);
+}
+#define CHECK(call, what) do { long r_ = (long)(call); if (r_ != 0) die(what, r_); } while (0)
+
+int main(int argc, char** argv) {
+  if (argc < 8) { fprintf(stderr, "usage: cufft_bench kind f32|f64 in out reps samples n1 [n2 [n3]] [bB]\n"); return 2; }
+  std::string kind = argv[1];
+  bool dbl = std::string(argv[2]) == "f64";
+  const char* fin = argv[3]; const char* fout = argv[4];
+  int reps = atoi(argv[5]), samples = atoi(argv[6]);
+  int dim = argc - 7;
+  long long batch = 1;
+  if (argv[argc - 1][0] == 'b') { batch = atoll(argv[argc - 1] + 1); dim--; }
+  if (dim < 1 || dim > 3) die("dimension", dim);
+  long long n[3] = {1, 1, 1};
+  for (int d = 0; d < dim; d++) n[d] = atoll(argv[7 + d]);
+  long long per = n[0] * n[1] * n[2], h = n[0] / 2 + 1, perh = h * (per / n[0]);
+  size_t es = dbl ? 8 : 4;
+  bool r2c = kind == "r2c", c2r = kind == "c2r";
+  if (!(r2c || c2r || kind == "c2c" || kind == "bc2c")) die("kind", 0);
+  size_t bytes_c = 2 * es * per * batch, bytes_r = es * per * batch, bytes_h = 2 * es * perh * batch;
+  size_t in_bytes = r2c ? bytes_r : c2r ? bytes_h : bytes_c;
+  size_t out_bytes = r2c ? bytes_h : c2r ? bytes_r : bytes_c;
+
+  std::vector<char> hin(in_bytes), hout(out_bytes);
+  FILE* f = fopen(fin, "rb"); if (!f || fread(hin.data(), 1, in_bytes, f) != in_bytes) die("reading input", 0); fclose(f);
+  CHECK(cudaSetDevice(0), "cudaSetDevice");
+  void *din, *dout;
+  CHECK(cudaMalloc(&din, in_bytes), "cudaMalloc"); CHECK(cudaMalloc(&dout, out_bytes), "cudaMalloc");
+  CHECK(cudaMemcpy(din, hin.data(), in_bytes, cudaMemcpyHostToDevice), "HtoD");
+  CHECK(cudaDeviceSynchronize(), "sync");
+
+  // cuFFT dimensions: slowest first (row-major), i.e. the column-major sizes reversed
+  long long nn[3];
+  for (int d = 0; d < dim; d++) nn[d] = n[dim - 1 - d];
+  cufftType type = r2c ? (dbl ? CUFFT_D2Z : CUFFT_R2C) : c2r ? (dbl ? CUFFT_Z2D : CUFFT_C2R) : (dbl ? CUFFT_Z2Z : CUFFT_C2C);
+  long long idist = r2c ? per : c2r ? perh : per, odist = r2c ? perh : c2r ? per : per;
+  cufftHandle plan;
+  size_t work = 0;
+  auto t0 = clk::now();
+  CHECK(cufftCreate(&plan), "cufftCreate");
+  CHECK(cufftMakePlanMany64(plan, dim, nn, nullptr, 1, idist, nullptr, 1, odist, type, batch, &work), "cufftMakePlanMany64");
+  double plan_ms = ms_since(t0);
+  int dir = kind == "bc2c" ? CUFFT_INVERSE : CUFFT_FORWARD;
+  auto run = [&]() {
+    if (type == CUFFT_C2C) CHECK(cufftExecC2C(plan, (cufftComplex*)din, (cufftComplex*)dout, dir), "exec");
+    else if (type == CUFFT_Z2Z) CHECK(cufftExecZ2Z(plan, (cufftDoubleComplex*)din, (cufftDoubleComplex*)dout, dir), "exec");
+    else if (type == CUFFT_R2C) CHECK(cufftExecR2C(plan, (cufftReal*)din, (cufftComplex*)dout), "exec");
+    else if (type == CUFFT_D2Z) CHECK(cufftExecD2Z(plan, (cufftDoubleReal*)din, (cufftDoubleComplex*)dout), "exec");
+    else if (type == CUFFT_C2R) CHECK(cufftExecC2R(plan, (cufftComplex*)din, (cufftReal*)dout), "exec");
+    else CHECK(cufftExecZ2D(plan, (cufftDoubleComplex*)din, (cufftDoubleReal*)dout), "exec");
+  };
+  t0 = clk::now();
+  run(); CHECK(cudaDeviceSynchronize(), "sync");
+  double first_ms = ms_since(t0);
+  CHECK(cudaMemcpy(hout.data(), dout, out_bytes, cudaMemcpyDeviceToHost), "DtoH");
+  f = fopen(fout, "wb"); fwrite(hout.data(), 1, out_bytes, f); fclose(f);
+
+  std::vector<double> ts;
+  for (int s = 0; s < samples; s++) {
+    if (c2r) CHECK(cudaMemcpy(din, hin.data(), in_bytes, cudaMemcpyHostToDevice), "HtoD");
+    run(); CHECK(cudaDeviceSynchronize(), "sync");
+    t0 = clk::now();
+    for (int r = 0; r < reps; r++) run();
+    CHECK(cudaDeviceSynchronize(), "sync");
+    ts.push_back(ms_since(t0) * 1e3 / reps);
+  }
+  std::sort(ts.begin(), ts.end());
+  printf("%.3f %.3f %.3f\n", plan_ms, first_ms, ts[ts.size() / 2]);
+  cufftDestroy(plan); cudaFree(din); cudaFree(dout);
+  return 0;
+}

bench/ducc-ffi/build.sh

@@ -0,0 +1,14 @@
+#!/bin/bash
+# Builds libducc_ffi against the ducc0 source tree, with native flags. Only ducc0/fft/fft.h is
+# included by ducc_ffi.cc; the instantiation units and infrastructure files are needed to link.
+#   bench/ducc-ffi/build.sh <ducc-src-dir (contains ducc0/)>
+set -euo pipefail
+src=$1
+here=$(cd "$(dirname "$0")" && pwd)
+arch=-march=native; [ "$(uname -m)" = arm64 ] && arch=-mcpu=native
+ext=so; [ "$(uname)" = Darwin ] && ext=dylib
+${CXX:-g++} -O3 $arch -std=c++17 -ffp-contract=fast -fPIC -shared -pthread -I"$src" \
+  "$here/ducc_ffi.cc" "$src/ducc0/fft/fft_inst1.cc" "$src/ducc0/fft/fft_inst2.cc" \
+  "$src/ducc0/infra/mav.cc" "$src/ducc0/infra/threading.cc" "$src/ducc0/infra/string_utils.cc" \
+  -o "$here/libducc_ffi.$ext"
+echo "built $here/libducc_ffi.$ext"