BENCHMARKS.md

@@ -482,6 +482,85 @@ The rows and columns of the uncompressed layouts are within 20% (chunked
 column 0.45 -> 0.49 ms, contiguous column 2.55 -> 2.61 ms). This run does not
 explain the slower windows.
 
+## Concurrent reads
+
+**Not yet measured.** The harness exists; no numbers are published until it
+has been run on an idle machine. The smoke runs used while building it (tiny
+files, other jobs compiling on the box) are not results.
+
+The question: libhdf5's threadsafe build serialises every API call under one
+global mutex, and h5py holds a global lock around every call too, so threads
+reading through h5py cannot decode in parallel; h5py users scale with
+processes. A clawhdf5 `File` is `Send + Sync`, and nothing on the read paths
+this harness uses (`read_f32`, `read_f32_selection`) takes a library-wide
+lock: the one mutex is the `File`'s chunk cache (keyed per dataset), taken by
+full reads of chunked datasets for each chunk's O(1) lookup and insert, never
+across a decode; hyperslab reads do not use the cache. How does
+decoded throughput scale with threads on one open file, against h5py threads
+and h5py processes on the same files?
+
+Workload (`crates/clawhdf5-bench/src/bin/concurrent_read.rs`; the h5py script
+mirrors it): `<dir>/deflate.h5` and `<dir>/contiguous.h5`, each with 64 `f32`
+datasets of 64 MiB decoded (`[16384, 1024]`; the deflate file chunked
+`256 x 256`, level 4), written by clawhdf5 on first use and reused while
+`manifest.json` matches. The data is a slowly varying ramp plus 8 bits of
+noise per element, every value exact in `f32`, so both harnesses check what
+they read; it deflates about 3.1x (128 MiB -> 40.7 MiB for two 64 MiB
+datasets). For each layout and thread count
+(1, 2, 4, 8, 16; fixed total work per repetition, split among the threads):
+
+- `distinct`: every dataset read in full once, thread `t` taking datasets
+  `t, t + T, ...`;
+- `same`: 1024 random `256 x 256` hyperslabs of `d00` in total, from a seeded
+  splitmix64 stream that both harnesses generate identically.
+
+Reported per row: MB/s of decoded (selected) data from the median of the
+repetitions, and scaling efficiency `MB/s(T) / (T x MB/s(1))`. Each worker
+times itself from a start barrier; a repetition spans the earliest start to
+the latest finish. Page cache: warm by default (each file is read once before
+timing); `--cold` evicts the files with `posix_fadvise(POSIX_FADV_DONTNEED)`
+before every repetition (no root needed; best effort). clawhdf5 opens one
+`File` per repetition, shared by all threads; h5py threads share one
+`h5py.File`; h5py processes (spawned before timing) each open the file inside
+the timed region.
+
+Decode inside a single clawhdf5 read is itself parallel in this binary
+(clawhdf5-format's `parallel` feature, enabled here through clawhdf5-agent;
+it is off in the facade's default features), so a 1-thread clawhdf5 full read
+of the deflate file already uses the whole rayon pool. Run both
+`--decode-threads 1` (each read decodes on its calling thread, like h5py —
+this isolates the API's own scaling) and the default pool.
+
+> **Run** (from the repository root). The default files take about 5.4 GiB
+> of disk (4 GiB contiguous + about 1.3 GiB deflate). Generating them is
+> memory-hungry because `FileBuilder` holds a whole file in memory: peak RSS
+> was 676 MB for `--datasets 2 --mib 64` (2026-09-25, tank,
+> `/usr/bin/time -f %M`), about 5x one file's decoded size, so expect about
+> 21 GB at the defaults (once; later runs reuse the files). Put `--dir` on a
+> real disk, not tmpfs, if `--cold` is to mean anything.
+>
+> ```bash
+> DIR=/path/on/disk/concurrent-read
+> BENCH=crates/clawhdf5-bench/scripts
+> PY=.venv/bin/python   # h5py 3.16 / HDF5 2.0 in this repo
+> cargo build --release -p clawhdf5-bench --bin concurrent_read
+> B=target/release/concurrent_read
+> $B --dir $DIR --json claw-pool.json                    # generates on first run
+> $B --dir $DIR --decode-threads 1 --json claw-1.json
+> $PY $BENCH/concurrent_read_h5py.py --dir $DIR --executor threads   --json h5py-threads.json
+> $PY $BENCH/concurrent_read_h5py.py --dir $DIR --executor processes --json h5py-procs.json
+> $PY $BENCH/compare_concurrent_read.py claw-1.json    h5py-threads.json h5py-procs.json
+> $PY $BENCH/compare_concurrent_read.py claw-pool.json h5py-threads.json h5py-procs.json
+> ```
+>
+> Cold page cache: add `--cold` to every harness command. Smoke test (seconds):
+> `$B --dir /tmp/cr --datasets 4 --mib 1 --threads 1,2,4 --slabs 16 --reps 1`
+> and the same `--threads/--slabs/--reps` to the h5py script.
+
+Other flags (both harnesses): `--threads`, `--reps`, `--slab`, `--slabs`,
+`--seed`, `--modes distinct,same`, `--layouts deflate,contiguous`; sizes
+(`--datasets`, `--mib`) only on the Rust harness, which writes the files.
+
 ## Search harness baseline (v2.3.0)
 
 Produced by `cargo run --release -p clawhdf5-bench --bin search_harness -- --full`

CHANGELOG.md

@@ -195,6 +195,13 @@
   README claimed but nothing measured.
 - `footprint_bench` reports whether it built `float16` or `f32` stores and
   takes `--f32`; it had kept printing "f32" after the default changed.
+- New `concurrent_read` harness, with an h5py counterpart
+  (`crates/clawhdf5-bench/scripts/concurrent_read_h5py.py`, threads or
+  processes) and `compare_concurrent_read.py`: decoded read throughput and
+  scaling efficiency at 1-16 threads on one open file, full reads of distinct
+  datasets and random hyperslabs of one dataset, deflate and contiguous, warm
+  or `--cold` page cache, JSON output. Not yet measured — `BENCHMARKS.md`
+  ("Concurrent reads") has the commands and no numbers.
 
 ### Interop
 - **Conformance sweep in the repo** (`conformance/`, report in

crates/clawhdf5-bench/Cargo.toml

@@ -34,6 +34,10 @@ path = "src/bin/consolidation_efficiency.rs"
 name = "ephemeral_perf"
 path = "src/bin/ephemeral_perf.rs"
 
+[[bin]]
+name = "concurrent_read"
+path = "src/bin/concurrent_read.rs"
+
 [[bin]]
 name = "mpi_io_bench"
 path = "src/bin/mpi_io_bench.rs"
@@ -64,6 +68,10 @@ clawhdf5-io = { path = "../clawhdf5-io" }
 mpi = { version = "0.8", optional = true }
 serde = { workspace = true }
 serde_json = "1"
+# concurrent_read: size the decode pool (--decode-threads) and evict files
+# from the page cache (--cold, posix_fadvise). Both pure Rust / bindings only.
+rayon = "1"
+libc = "0.2"
 tempfile = { workspace = true }
 # Optional: libhdf5 C wrapper for side-by-side comparison (requires system libhdf5).
 # Enable with: cargo bench -p clawhdf5-bench --features libhdf5-compare

crates/clawhdf5-bench/scripts/compare_concurrent_read.py

@@ -0,0 +1,70 @@
+#!/usr/bin/env python3
+"""Tabulate concurrent_read JSON results (clawhdf5, h5py threads/processes).
+
+    python compare_concurrent_read.py clawhdf5.json h5py-threads.json h5py-procs.json
+
+Prints one Markdown table: for each layout, mode and thread count, every
+tool's MB/s and scaling efficiency, and the first file's MB/s relative to each
+of the others. Refuses to compare runs whose workload parameters differ.
+"""
+
+import json
+import sys
+
+COMPARED = ("datasets", "rows", "cols", "chunk", "deflate_level", "slab", "slabs", "seed")
+
+
+def main(paths):
+    if len(paths) < 2:
+        sys.exit(__doc__)
+    docs = []
+    for p in paths:
+        with open(p) as fh:
+            docs.append(json.load(fh))
+    ref = docs[0]
+    for d, p in zip(docs[1:], paths[1:]):
+        diff = [k for k in COMPARED if d["params"].get(k) != ref["params"].get(k)]
+        if diff:
+            sys.exit(f"{p}: workload differs from {paths[0]} in {', '.join(diff)}")
+        if d["cache"] != ref["cache"]:
+            print(f"warning: {p} ran {d['cache']!r}, {paths[0]} ran {ref['cache']!r}",
+                  file=sys.stderr)
+        if d.get("host") != ref.get("host"):
+            print(f"warning: {p} ran on {d.get('host')}, {paths[0]} on {ref.get('host')}",
+                  file=sys.stderr)
+
+    names = [d["tool"] for d in docs]
+    for d in docs:
+        extra = f", HDF5 {d['hdf5_version']}" if "hdf5_version" in d else ""
+        print(f"- {d['tool']} {d['version']}{extra}: host {d.get('host')}, "
+              f"{d.get('cpus')} CPUs, cache {d['cache']}, decode threads per read "
+              f"{d.get('decode_threads')}")
+    p = ref["params"]
+    print(f"\n{p['datasets']} datasets of {p['rows']} x {p['cols']} f32, chunks "
+          f"{p['chunk'][0]} x {p['chunk'][1]} (deflate {p['deflate_level']}); "
+          f"`same`: {p['slabs']} slabs of {p['slab']} x {p['slab']}\n")
+
+    index = [{(r["layout"], r["mode"], r["threads"]): r for r in d["results"]} for d in docs]
+    keys = [(r["layout"], r["mode"], r["threads"]) for r in ref["results"]]
+
+    head = ["layout", "mode", "threads"]
+    head += [f"{n} MB/s (eff)" for n in names]
+    head += [f"{names[0]} / {n}" for n in names[1:]]
+    print("| " + " | ".join(head) + " |")
+    print("|---|---|" + "---:|" * (len(head) - 2))
+    for key in keys:
+        cells = [key[0], key[1], str(key[2])]
+        rs = [ix.get(key) for ix in index]
+        for r in rs:
+            if r is None:
+                cells.append("-")
+            else:
+                eff = "-" if r["efficiency"] is None else f"{r['efficiency']:.2f}"
+                cells.append(f"{r['mb_s']:.0f} ({eff})")
+        for r in rs[1:]:
+            cells.append("-" if r is None else f"{rs[0]['mb_s'] / r['mb_s']:.2f}x")
+        print("| " + " | ".join(cells) + " |")
+
+
+if __name__ == "__main__":
+    main(sys.argv[1:])

crates/clawhdf5-bench/scripts/concurrent_read_h5py.py

@@ -0,0 +1,265 @@
+#!/usr/bin/env python3
+"""The concurrent_read workload with h5py, on the files concurrent_read wrote.
+
+libhdf5 serialises every API call under one global lock, and h5py holds its
+own global lock around every call as well, so h5py *threads* cannot decode in
+parallel. h5py users scale with *processes* instead; ``--executor processes``
+measures that (each worker opens the file itself).
+
+The workload mirrors ``crates/clawhdf5-bench/src/bin/concurrent_read.rs``:
+
+* ``distinct``: every dataset read in full once per repetition; worker ``t``
+  of ``T`` reads datasets ``t, t + T, ...``.
+* ``same``: ``--slabs`` random ``--slab`` x ``--slab`` hyperslabs of ``d00``
+  (slab ``j`` to worker ``j % T``), offsets from the same splitmix64 stream.
+
+Each worker times itself from a start barrier; a repetition spans the earliest
+start to the latest finish (CLOCK_MONOTONIC, comparable across processes).
+Threads share one ``h5py.File`` per repetition; process workers open the file
+inside the timed region (a few ms against reads of many MiB).
+
+Generate the files first with the Rust harness (it writes ``manifest.json``),
+then, for example::
+
+    python concurrent_read_h5py.py --dir DIR --executor threads --json h5py-threads.json
+    python concurrent_read_h5py.py --dir DIR --executor processes --json h5py-procs.json
+"""
+
+import argparse
+import json
+import multiprocessing as mp
+import os
+import platform
+import socket
+import sys
+import threading
+import time
+
+import h5py
+import numpy as np
+
+M64 = (1 << 64) - 1
+
+
+def splitmix64(state):
+    """Return (new_state, value); the same stream as the Rust harness."""
+    state = (state + 0x9E3779B97F4A7C15) & M64
+    z = state
+    z = ((z ^ (z >> 30)) * 0xBF58476D1CE4E5B9) & M64
+    z = ((z ^ (z >> 27)) * 0x94D049BB133111EB) & M64
+    return state, z ^ (z >> 31)
+
+
+def value(k, i):
+    """Element i (row-major) of dataset k, exactly as concurrent_read writes it."""
+    _, noise = splitmix64(i ^ (k << 40))
+    return np.float32((((i >> 6) % 16384) + k) + (noise & 0xFF) / 256.0)
+
+
+def slab_offsets(seed, count, rows, cols, slab):
+    s = seed
+    out = []
+    for _ in range(count):
+        s, r = splitmix64(s)
+        s, c = splitmix64(s)
+        out.append((r % (rows - slab + 1), c % (cols - slab + 1)))
+    return out
+
+
+def now():
+    return time.clock_gettime(time.CLOCK_MONOTONIC)
+
+
+def work(f, mode, t, threads, m, slabs, slab, verify):
+    """Worker t's share of one repetition on an open h5py.File."""
+    n = m["rows"] * m["cols"]
+    if mode == "distinct":
+        for k in range(t, m["datasets"], threads):
+            got = f[f"d{k:02d}"][...]
+            assert got.size == n
+            if verify:
+                flat = got.reshape(-1)
+                for i in (0, n // 3, n - 1):
+                    assert flat[i] == value(k, i), f"d{k:02d}[{i}]"
+    else:
+        ds = f["d00"]
+        cols = m["cols"]
+        for r, c in slabs[t::threads]:
+            got = ds[r : r + slab, c : c + slab]
+            assert got.shape == (slab, slab)
+            if verify:
+                assert got[0, 0] == value(0, r * cols + c)
+                last = (r + slab - 1) * cols + c + slab - 1
+                assert got[-1, -1] == value(0, last)
+
+
+# ----- process workers ------------------------------------------------------
+
+_barrier = None
+
+
+def _init(barrier):
+    global _barrier
+    _barrier = barrier
+
+
+def _proc_task(task):
+    path, mode, t, threads, m, slabs, slab = task
+    _barrier.wait()
+    start = now()
+    with h5py.File(path, "r") as f:
+        work(f, mode, t, threads, m, slabs, slab, False)
+    return start, now()
+
+
+def _noop(_):
+    return os.getpid()
+
+
+def run_threads(path, mode, threads, m, slabs, slab):
+    spans = [None] * threads
+    barrier = threading.Barrier(threads)
+    with h5py.File(path, "r") as f:
+
+        def body(t):
+            barrier.wait()
+            start = now()
+            work(f, mode, t, threads, m, slabs, slab, False)
+            spans[t] = (start, now())
+
+        ts = [threading.Thread(target=body, args=(t,)) for t in range(threads)]
+        for th in ts:
+            th.start()
+        for th in ts:
+            th.join()
+    return max(e for _, e in spans) - min(s for s, _ in spans)
+
+
+def run_processes(pool, path, mode, threads, m, slabs, slab):
+    tasks = [(path, mode, t, threads, m, slabs, slab) for t in range(threads)]
+    # One task per worker: each blocks in the barrier until all T have
+    # started, so no worker can take a second task.
+    spans = pool.map(_proc_task, tasks, chunksize=1)
+    return max(e for _, e in spans) - min(s for s, _ in spans)
+
+
+def warm(path):
+    with open(path, "rb") as fh:
+        while fh.read(1 << 24):
+            pass
+
+
+def evict(path):
+    fd = os.open(path, os.O_RDONLY)
+    try:
+        os.posix_fadvise(fd, 0, 0, os.POSIX_FADV_DONTNEED)
+    finally:
+        os.close(fd)
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__.split("\n\n")[0])
+    ap.add_argument("--dir", default="concurrent-read-data")
+    ap.add_argument("--executor", choices=["threads", "processes"], default="threads")
+    ap.add_argument("--threads", default="1,2,4,8,16")
+    ap.add_argument("--reps", type=int, default=3)
+    ap.add_argument("--slab", type=int, default=256)
+    ap.add_argument("--slabs", type=int, default=1024)
+    ap.add_argument("--seed", type=int, default=42)
+    ap.add_argument("--cold", action="store_true")
+    ap.add_argument("--modes", default="distinct,same")
+    ap.add_argument("--layouts", default="deflate,contiguous")
+    ap.add_argument("--json")
+    a = ap.parse_args()
+
+    # The Rust harness pins this value (splitmix64_reference).
+    assert splitmix64(42)[1] == 0xBDD732262FEB6E95, "splitmix64 port is wrong"
+
+    try:
+        with open(os.path.join(a.dir, "manifest.json")) as fh:
+            m = json.load(fh)
+    except FileNotFoundError:
+        sys.exit(f"{a.dir}/manifest.json not found: generate the files with "
+                 "`cargo run --release -p clawhdf5-bench --bin concurrent_read -- --dir ...` first")
+    threads_list = [int(x) for x in a.threads.split(",")]
+    modes = a.modes.split(",")
+    layouts = a.layouts.split(",")
+    if a.slab < 1 or a.slab > min(m["rows"], m["cols"]):
+        sys.exit(f"--slab must be 1..={min(m['rows'], m['cols'])}")
+    files = dict(m["files"])
+    slabs = slab_offsets(a.seed, a.slabs, m["rows"], m["cols"], a.slab)
+    dataset_bytes = m["rows"] * m["cols"] * 4
+    tool = f"h5py-{a.executor}"
+
+    ctx = mp.get_context("spawn")  # never fork a process holding HDF5 state
+    pools = {}
+    if a.executor == "processes":
+        for t in threads_list:
+            pool = ctx.Pool(t, initializer=_init, initargs=(ctx.Barrier(t),))
+            pool.map(_noop, range(t))  # start the workers outside the timing
+            pools[t] = pool
+
+    rows = []
+    print("| layout | mode | threads | MB/s | efficiency | median s |")
+    print("|---|---|---:|---:|---:|---:|")
+    try:
+        for layout in layouts:
+            path = os.path.join(a.dir, files[layout])
+            if not a.cold:
+                warm(path)
+            for mode in modes:
+                with h5py.File(path, "r") as f:  # untimed, checked pass
+                    work(f, mode, 0, 1, m, slabs, a.slab, True)
+                nbytes = (dataset_bytes * m["datasets"] if mode == "distinct"
+                          else a.slab * a.slab * 4 * a.slabs)
+                base = None
+                for t in threads_list:
+                    times = []
+                    for _ in range(a.reps):
+                        if a.cold:
+                            evict(path)
+                        if a.executor == "threads":
+                            times.append(run_threads(path, mode, t, m, slabs, a.slab))
+                        else:
+                            times.append(run_processes(pools[t], path, mode, t, m, slabs, a.slab))
+                    med = sorted(times)[len(times) // 2]
+                    mb_s = nbytes / (1 << 20) / med
+                    if t == 1:
+                        base = mb_s
+                    eff = mb_s / (t * base) if base else None
+                    print(f"| {layout} | {mode} | {t} | {mb_s:.0f} | "
+                          f"{'-' if eff is None else f'{eff:.2f}'} | {med:.4f} |")
+                    rows.append({
+                        "layout": layout, "mode": mode, "threads": t, "bytes": nbytes,
+                        "times_s": times, "median_s": med, "mb_s": mb_s, "efficiency": eff,
+                    })
+    finally:
+        for pool in pools.values():
+            pool.terminate()
+
+    if a.json:
+        doc = {
+            "tool": tool,
+            "version": h5py.__version__,
+            "hdf5_version": h5py.version.hdf5_version,
+            "python": platform.python_version(),
+            "host": socket.gethostname(),
+            "cpus": os.cpu_count(),
+            "unix_time": int(time.time()),
+            "cache": ("cold (posix_fadvise DONTNEED before each repetition)"
+                      if a.cold else "warm"),
+            "decode_threads": 1,
+            "params": {
+                "datasets": m["datasets"], "rows": m["rows"], "cols": m["cols"],
+                "chunk": m["chunk"], "deflate_level": m["deflate_level"],
+                "mib": dataset_bytes // (1 << 20), "slab": a.slab, "slabs": a.slabs,
+                "seed": a.seed, "reps": a.reps, "dir": a.dir,
+            },
+            "results": rows,
+        }
+        with open(a.json, "w") as fh:
+            json.dump(doc, fh, indent=2)
+
+
+if __name__ == "__main__":
+    main()

crates/clawhdf5-bench/src/bin/concurrent_read.rs

@@ -0,0 +1,523 @@
+//! Concurrent-read harness: how does decoded read throughput scale with the
+//! number of threads reading one open file?
+//!
+//! libhdf5 (threadsafe build) serialises every API call under one global
+//! mutex, and h5py holds it too, so threads cannot decode in parallel there.
+//! A clawhdf5 [`File`] is `Send + Sync`; this harness measures what that buys.
+//! `crates/clawhdf5-bench/scripts/concurrent_read_h5py.py` runs the same
+//! workload on the same files with h5py (threads, and processes), and
+//! `compare_concurrent_read.py` tabulates the JSON both write.
+//!
+//! Files (generated on first use, reused while `manifest.json` matches):
+//!
+//! * `<dir>/deflate.h5`: `--datasets` datasets `d00`, `d01`, ... of `f32`,
+//!   `--mib` MiB decoded each, shape `[mib * 256, 1024]`, chunks `256 x 256`,
+//!   deflate level 4.
+//! * `<dir>/contiguous.h5`: the same datasets, contiguous.
+//!
+//! Modes, for each layout and each thread count `T` (strong scaling: the total
+//! work per repetition is fixed, split among the threads):
+//!
+//! * `distinct`: every dataset is read in full once; thread `t` reads datasets
+//!   `t, t + T, t + 2T, ...`.
+//! * `same`: all threads read `d00`, `--slabs` random `--slab` x `--slab`
+//!   hyperslabs in total (slab `j` goes to thread `j % T`). The offsets come
+//!   from a splitmix64 stream seeded with `--seed`, identical in the h5py
+//!   script.
+//!
+//! One `File` per layout per repetition is shared by all threads (opened
+//! fresh each repetition, so no chunk cache carries over). Page cache:
+//! `warm` (default) reads every file once before timing; `--cold` evicts the
+//! files from the page cache with `posix_fadvise(POSIX_FADV_DONTNEED)` before
+//! every repetition (no root needed; it only evicts clean, unmapped pages, so
+//! it is best effort — the JSON says which was used).
+//!
+//! Decode inside one read is itself parallel when clawhdf5-format's `parallel`
+//! feature is on (it is in this binary, via clawhdf5-agent). `--decode-threads
+//! N` sizes that rayon pool; `--decode-threads 1` measures the API's own
+//! thread scaling, comparable with h5py where each call decodes on the
+//! calling thread.
+//!
+//! ```text
+//! cargo run --release -p clawhdf5-bench --bin concurrent_read -- \
+//!     --dir /data/concurrent-read --json clawhdf5.json
+//! cargo run --release -p clawhdf5-bench --bin concurrent_read -- \
+//!     --dir /tmp/cr --datasets 4 --mib 1 --threads 1,2 --slabs 16 --reps 1   # smoke
+//! ```
+
+use std::path::{Path, PathBuf};
+use std::sync::Barrier;
+use std::time::Instant;
+
+use clawhdf5::{File, FileBuilder, Selection};
+use serde::{Deserialize, Serialize};
+
+const COLS: u64 = 1024;
+const ROWS_PER_MIB: u64 = 256; // 256 rows x 1024 cols x 4 bytes = 1 MiB
+const CHUNK: u64 = 256;
+const DEFLATE_LEVEL: u32 = 4;
+const LAYOUTS: [&str; 2] = ["deflate", "contiguous"];
+const MANIFEST_VERSION: u32 = 1;
+
+/// splitmix64 — shared with the h5py script, which must produce the same
+/// stream (both the data and the hyperslab offsets depend on it).
+fn splitmix64(state: &mut u64) -> u64 {
+    *state = state.wrapping_add(0x9E37_79B9_7F4A_7C15);
+    let mut z = *state;
+    z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
+    z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
+    z ^ (z >> 31)
+}
+
+/// Element `i` (row-major) of dataset `k`: a slowly varying integer part plus
+/// 8 bits of noise, so deflate has real work to do (about 3.1x) and every value
+/// is exact in `f32` (< 2^15 with 8 fraction bits), which lets both harnesses
+/// check what they read against this formula.
+fn value(k: u64, i: u64) -> f32 {
+    let mut s = i ^ (k << 40);
+    let noise = splitmix64(&mut s) & 0xff;
+    (((i >> 6) % 16384) + k) as f32 + noise as f32 / 256.0
+}
+
+#[derive(Serialize, Deserialize, PartialEq, Debug, Clone)]
+struct Manifest {
+    version: u32,
+    datasets: u64,
+    rows: u64,
+    cols: u64,
+    chunk: [u64; 2],
+    deflate_level: u32,
+    files: Vec<(String, String)>, // (layout, file name)
+    writer: String,
+}
+
+fn manifest_for(datasets: u64, mib: u64) -> Manifest {
+    Manifest {
+        version: MANIFEST_VERSION,
+        datasets,
+        rows: mib * ROWS_PER_MIB,
+        cols: COLS,
+        chunk: [CHUNK, CHUNK],
+        deflate_level: DEFLATE_LEVEL,
+        files: LAYOUTS
+            .iter()
+            .map(|l| (l.to_string(), format!("{l}.h5")))
+            .collect(),
+        writer: format!("clawhdf5 {}", env!("CARGO_PKG_VERSION")),
+    }
+}
+
+fn dataset_values(k: u64, n: u64) -> Vec<f32> {
+    (0..n).map(|i| value(k, i)).collect()
+}
+
+/// Write the files unless `dir` already holds ones matching `want`.
+fn ensure_files(dir: &Path, want: &Manifest) -> std::io::Result<bool> {
+    let manifest_path = dir.join("manifest.json");
+    if let Ok(text) = std::fs::read_to_string(&manifest_path)
+        && let Ok(have) = serde_json::from_str::<Manifest>(&text)
+        && have.version == want.version
+        && have.datasets == want.datasets
+        && have.rows == want.rows
+        && have.cols == want.cols
+        && have.chunk == want.chunk
+        && have.deflate_level == want.deflate_level
+        && have.files == want.files
+        && want.files.iter().all(|(_, f)| dir.join(f).exists())
+    {
+        return Ok(false);
+    }
+    std::fs::create_dir_all(dir)?;
+    // A stale manifest must not survive a half-written regeneration.
+    let _ = std::fs::remove_file(&manifest_path);
+    let n = want.rows * want.cols;
+    for (layout, file) in &want.files {
+        // One layout at a time keeps the peak memory to about twice one
+        // file's decoded size.
+        let mut b = FileBuilder::new();
+        for k in 0..want.datasets {
+            let ds = b.create_dataset(&format!("d{k:02}"));
+            ds.with_f32_data(&dataset_values(k, n))
+                .with_shape(&[want.rows, want.cols]);
+            if layout == "deflate" {
+                ds.with_chunks(&[CHUNK.min(want.rows), CHUNK])
+                    .with_deflate(DEFLATE_LEVEL);
+            }
+        }
+        b.write(dir.join(file)).map_err(std::io::Error::other)?;
+    }
+    std::fs::write(
+        &manifest_path,
+        serde_json::to_string_pretty(want).map_err(std::io::Error::other)?,
+    )?;
+    Ok(true)
+}
+
+fn slab_offsets(seed: u64, count: usize, rows: u64, cols: u64, slab: u64) -> Vec<(u64, u64)> {
+    let mut s = seed;
+    (0..count)
+        .map(|_| {
+            let r = splitmix64(&mut s) % (rows - slab + 1);
+            let c = splitmix64(&mut s) % (cols - slab + 1);
+            (r, c)
+        })
+        .collect()
+}
+
+/// Warm the page cache by reading every byte of `path`.
+fn warm(path: &Path) -> std::io::Result<()> {
+    let mut f = std::fs::File::open(path)?;
+    std::io::copy(&mut f, &mut std::io::sink())?;
+    Ok(())
+}
+
+/// Ask the kernel to drop `path`'s pages from the page cache.
+fn evict(path: &Path) -> std::io::Result<()> {
+    use std::os::fd::AsRawFd;
+    let f = std::fs::File::open(path)?;
+    // SAFETY: plain syscall on a valid, open file descriptor.
+    let rc = unsafe { libc::posix_fadvise(f.as_raw_fd(), 0, 0, libc::POSIX_FADV_DONTNEED) };
+    if rc != 0 {
+        return Err(std::io::Error::from_raw_os_error(rc));
+    }
+    Ok(())
+}
+
+#[derive(Serialize)]
+struct Row {
+    layout: String,
+    mode: String,
+    threads: usize,
+    /// Decoded (selected) bytes read per repetition.
+    bytes: u64,
+    times_s: Vec<f64>,
+    median_s: f64,
+    mb_s: f64,
+    /// `mb_s / (threads * mb_s at threads = 1)`; null without a 1-thread row.
+    efficiency: Option<f64>,
+}
+
+struct Args {
+    dir: PathBuf,
+    datasets: u64,
+    mib: u64,
+    threads: Vec<usize>,
+    reps: usize,
+    slab: u64,
+    slabs: usize,
+    seed: u64,
+    cold: bool,
+    decode_threads: usize,
+    modes: Vec<String>,
+    layouts: Vec<String>,
+    json: Option<PathBuf>,
+}
+
+const USAGE: &str = "\
+usage: concurrent_read [--dir DIR] [--datasets N] [--mib N] [--threads 1,2,4,8,16]
+                       [--reps N] [--slab N] [--slabs N] [--seed N] [--cold]
+                       [--decode-threads N] [--modes distinct,same]
+                       [--layouts deflate,contiguous] [--json FILE]";
+
+fn parse_list<T: std::str::FromStr>(s: &str) -> Result<Vec<T>, String> {
+    s.split(',')
+        .map(|x| x.trim().parse().map_err(|_| format!("bad list item {x:?}")))
+        .collect()
+}
+
+fn parse_args() -> Result<Args, String> {
+    let mut a = Args {
+        dir: PathBuf::from("concurrent-read-data"),
+        datasets: 64,
+        mib: 64,
+        threads: vec![1, 2, 4, 8, 16],
+        reps: 3,
+        slab: 256,
+        slabs: 1024,
+        seed: 42,
+        cold: false,
+        decode_threads: 0,
+        modes: vec!["distinct".into(), "same".into()],
+        layouts: LAYOUTS.iter().map(|s| s.to_string()).collect(),
+        json: None,
+    };
+    let mut it = std::env::args().skip(1);
+    while let Some(flag) = it.next() {
+        if flag == "--cold" {
+            a.cold = true;
+            continue;
+        }
+        if flag == "-h" || flag == "--help" {
+            return Err(USAGE.into());
+        }
+        let v = it.next().ok_or(format!("{flag} needs a value\n{USAGE}"))?;
+        let num = |v: &str| {
+            v.parse::<u64>()
+                .map_err(|_| format!("{flag}: bad number {v:?}"))
+        };
+        match flag.as_str() {
+            "--dir" => a.dir = v.into(),
+            "--datasets" => a.datasets = num(&v)?,
+            "--mib" => a.mib = num(&v)?,
+            "--threads" => a.threads = parse_list(&v)?,
+            "--reps" => a.reps = num(&v)? as usize,
+            "--slab" => a.slab = num(&v)?,
+            "--slabs" => a.slabs = num(&v)? as usize,
+            "--seed" => a.seed = num(&v)?,
+            "--decode-threads" => a.decode_threads = num(&v)? as usize,
+            "--modes" => a.modes = parse_list(&v)?,
+            "--layouts" => a.layouts = parse_list(&v)?,
+            "--json" => a.json = Some(v.into()),
+            _ => return Err(format!("unknown flag {flag}\n{USAGE}")),
+        }
+    }
+    if a.datasets == 0 || a.datasets > 100 {
+        return Err("--datasets must be 1..=100".into());
+    }
+    if a.mib == 0 || a.reps == 0 || a.slabs == 0 || a.threads.contains(&0) {
+        return Err("--mib, --reps, --slabs and every --threads value must be > 0".into());
+    }
+    if a.slab == 0 || a.slab > COLS || a.slab > a.mib * ROWS_PER_MIB {
+        return Err(format!(
+            "--slab must be 1..={}",
+            COLS.min(a.mib * ROWS_PER_MIB)
+        ));
+    }
+    for m in &a.modes {
+        if m != "distinct" && m != "same" {
+            return Err(format!("unknown mode {m:?}"));
+        }
+    }
+    for l in &a.layouts {
+        if !LAYOUTS.contains(&l.as_str()) {
+            return Err(format!("unknown layout {l:?}"));
+        }
+    }
+    Ok(a)
+}
+
+/// One timed repetition: `T` threads on one shared `File`. Returns seconds.
+fn run_once(
+    path: &Path,
+    mode: &str,
+    threads: usize,
+    m: &Manifest,
+    slabs: &[(u64, u64)],
+    slab: u64,
+    verify: bool,
+) -> f64 {
+    let file = File::open(path).expect("open");
+    let barrier = Barrier::new(threads + 1); // + the spawning thread
+    let n = m.rows * m.cols;
+    // Each thread times itself from the barrier; the repetition spans the
+    // earliest start to the latest finish (timing on the spawning thread
+    // instead undercounts whenever it is scheduled after the workers ran).
+    let spans: Vec<(Instant, Instant)> = std::thread::scope(|s| {
+        let handles: Vec<_> = (0..threads)
+            .map(|t| {
+                let (file, barrier) = (&file, &barrier);
+                s.spawn(move || {
+                    barrier.wait();
+                    let start = Instant::now();
+                    match mode {
+                        "distinct" => {
+                            for k in (t as u64..m.datasets).step_by(threads) {
+                                let got = file.dataset(&format!("d{k:02}")).unwrap().read_f32();
+                                let got = got.unwrap();
+                                assert_eq!(got.len() as u64, n);
+                                if verify {
+                                    for i in [0, n / 3, n - 1] {
+                                        assert_eq!(got[i as usize], value(k, i), "d{k:02}[{i}]");
+                                    }
+                                }
+                                std::hint::black_box(got);
+                            }
+                        }
+                        _ => {
+                            let ds = file.dataset("d00").unwrap();
+                            for &(r, c) in slabs.iter().skip(t).step_by(threads) {
+                                let sel = Selection::Hyperslab {
+                                    start: vec![r, c],
+                                    stride: vec![1, 1],
+                                    count: vec![slab, slab],
+                                    block: vec![1, 1],
+                                };
+                                let got = ds.read_f32_selection(&sel).unwrap();
+                                assert_eq!(got.len() as u64, slab * slab);
+                                if verify {
+                                    let last = (r + slab - 1) * m.cols + c + slab - 1;
+                                    assert_eq!(got[0], value(0, r * m.cols + c));
+                                    assert_eq!(*got.last().unwrap(), value(0, last));
+                                }
+                                std::hint::black_box(got);
+                            }
+                        }
+                    }
+                    (start, Instant::now())
+                })
+            })
+            .collect();
+        barrier.wait();
+        handles.into_iter().map(|h| h.join().unwrap()).collect()
+    });
+    let start = spans.iter().map(|s| s.0).min().unwrap();
+    let end = spans.iter().map(|s| s.1).max().unwrap();
+    (end - start).as_secs_f64()
+}
+
+fn median(v: &[f64]) -> f64 {
+    let mut s = v.to_vec();
+    s.sort_by(f64::total_cmp);
+    s[s.len() / 2]
+}
+
+fn hostname() -> String {
+    std::fs::read_to_string("/proc/sys/kernel/hostname")
+        .map(|s| s.trim().to_string())
+        .unwrap_or_else(|_| "unknown".into())
+}
+
+fn main() {
+    let args = match parse_args() {
+        Ok(a) => a,
+        Err(e) => {
+            eprintln!("{e}");
+            std::process::exit(2);
+        }
+    };
+    if cfg!(debug_assertions) {
+        eprintln!("warning: debug build — numbers are meaningless. Use --release.");
+    }
+    if args.decode_threads > 0 {
+        rayon::ThreadPoolBuilder::new()
+            .num_threads(args.decode_threads)
+            .build_global()
+            .expect("configure rayon pool");
+    }
+
+    let manifest = manifest_for(args.datasets, args.mib);
+    let t = Instant::now();
+    match ensure_files(&args.dir, &manifest) {
+        Ok(true) => eprintln!(
+            "generated {} in {:.1} s",
+            args.dir.display(),
+            t.elapsed().as_secs_f64()
+        ),
+        Ok(false) => eprintln!("reusing {}", args.dir.display()),
+        Err(e) => {
+            eprintln!("cannot write test files in {}: {e}", args.dir.display());
+            std::process::exit(1);
+        }
+    }
+    let path_of = |layout: &str| args.dir.join(format!("{layout}.h5"));
+    let slabs = slab_offsets(
+        args.seed,
+        args.slabs,
+        manifest.rows,
+        manifest.cols,
+        args.slab,
+    );
+    let dataset_bytes = manifest.rows * manifest.cols * 4;
+
+    let mut rows: Vec<Row> = Vec::new();
+    println!("| layout | mode | threads | MB/s | efficiency | median s |");
+    println!("|---|---|---:|---:|---:|---:|");
+    for layout in &args.layouts {
+        let path = path_of(layout);
+        // Untimed pass: page cache warm (unless --cold), results checked.
+        if !args.cold {
+            warm(&path).expect("warm page cache");
+        }
+        for mode in &args.modes {
+            run_once(&path, mode, 1, &manifest, &slabs, args.slab, true);
+            let bytes = match mode.as_str() {
+                "distinct" => dataset_bytes * manifest.datasets,
+                _ => args.slab * args.slab * 4 * args.slabs as u64,
+            };
+            let mut base: Option<f64> = None;
+            for &threads in &args.threads {
+                let times: Vec<f64> = (0..args.reps)
+                    .map(|_| {
+                        if args.cold {
+                            evict(&path).expect("posix_fadvise");
+                        }
+                        run_once(&path, mode, threads, &manifest, &slabs, args.slab, false)
+                    })
+                    .collect();
+                let med = median(&times);
+                let mb_s = bytes as f64 / (1 << 20) as f64 / med;
+                if threads == 1 {
+                    base = Some(mb_s);
+                }
+                let efficiency = base.map(|b| mb_s / (threads as f64 * b));
+                println!(
+                    "| {layout} | {mode} | {threads} | {mb_s:.0} | {} | {med:.4} |",
+                    efficiency.map_or("-".into(), |e| format!("{e:.2}"))
+                );
+                rows.push(Row {
+                    layout: layout.clone(),
+                    mode: mode.clone(),
+                    threads,
+                    bytes,
+                    times_s: times,
+                    median_s: med,
+                    mb_s,
+                    efficiency,
+                });
+            }
+        }
+    }
+
+    if let Some(out) = &args.json {
+        let doc = serde_json::json!({
+            "tool": "clawhdf5",
+            "version": env!("CARGO_PKG_VERSION"),
+            "host": hostname(),
+            "cpus": std::thread::available_parallelism().map_or(0, |n| n.get()),
+            "unix_time": std::time::SystemTime::now()
+                .duration_since(std::time::UNIX_EPOCH)
+                .map_or(0, |d| d.as_secs()),
+            "cache": if args.cold { "cold (posix_fadvise DONTNEED before each repetition)" } else { "warm" },
+            "decode_threads": rayon::current_num_threads(),
+            "params": {
+                "datasets": manifest.datasets,
+                "mib": args.mib,
+                "rows": manifest.rows,
+                "cols": manifest.cols,
+                "chunk": manifest.chunk,
+                "deflate_level": manifest.deflate_level,
+                "slab": args.slab,
+                "slabs": args.slabs,
+                "seed": args.seed,
+                "reps": args.reps,
+                "dir": args.dir,
+            },
+            "results": rows,
+        });
+        std::fs::write(out, serde_json::to_string_pretty(&doc).unwrap()).expect("write json");
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn values_are_exact_in_f32() {
+        for k in [0, 7, 63] {
+            for i in [0u64, 1, 4095, 1 << 20, (1 << 24) - 1] {
+                let v = value(k, i);
+                assert_eq!(v, (v as f64) as f32);
+                assert!(v < 32768.0);
+                assert_eq!((v * 256.0).fract(), 0.0);
+            }
+        }
+    }
+
+    /// The h5py script hard-codes this vector to check its splitmix64 port.
+    #[test]
+    fn splitmix64_reference() {
+        let mut s = 42;
+        assert_eq!(splitmix64(&mut s), 0xBDD7_3226_2FEB_6E95);
+    }
+}

crates/clawhdf5-bench/tests/concurrent_read_smoke.rs

@@ -0,0 +1,148 @@
+//! Keeps the concurrent-read harnesses working: runs `concurrent_read`, the
+//! h5py script (threads and processes) and the comparison script end to end
+//! on tiny files. h5py reading the files also checks, element by element at
+//! spot positions, that both harnesses generate the same data and slabs.
+//!
+//! The h5py half is skipped when python3 with h5py is unavailable, unless
+//! `CLAWHDF5_REQUIRE_INTEROP=1`; `CLAWHDF5_PYTHON` picks the interpreter.
+
+use std::path::{Path, PathBuf};
+use std::process::Command;
+
+fn python() -> String {
+    std::env::var("CLAWHDF5_PYTHON").unwrap_or_else(|_| "python3".to_string())
+}
+
+fn interop_required() -> bool {
+    std::env::var("CLAWHDF5_REQUIRE_INTEROP").is_ok_and(|v| v == "1")
+}
+
+fn python_available() -> bool {
+    Command::new(python())
+        .args(["-c", "import h5py, numpy"])
+        .output()
+        .map(|o| o.status.success())
+        .unwrap_or(false)
+}
+
+fn scripts() -> PathBuf {
+    Path::new(env!("CARGO_MANIFEST_DIR")).join("scripts")
+}
+
+fn run(cmd: &mut Command) -> String {
+    let out = cmd.output().expect("spawn");
+    assert!(
+        out.status.success(),
+        "{cmd:?} failed\nSTDOUT:\n{}\nSTDERR:\n{}",
+        String::from_utf8_lossy(&out.stdout),
+        String::from_utf8_lossy(&out.stderr)
+    );
+    String::from_utf8_lossy(&out.stdout).into_owned()
+}
+
+const SMALL: [&str; 8] = [
+    "--threads",
+    "1,2",
+    "--slabs",
+    "8",
+    "--reps",
+    "1",
+    "--slab",
+    "64",
+];
+
+fn results(path: &Path) -> serde_json::Value {
+    serde_json::from_str(&std::fs::read_to_string(path).unwrap()).unwrap()
+}
+
+#[test]
+fn harnesses_run_end_to_end_on_tiny_files() {
+    let dir = tempfile::TempDir::new().unwrap();
+    let data = dir.path().join("data");
+    let claw = dir.path().join("claw.json");
+
+    let bin = env!("CARGO_BIN_EXE_concurrent_read");
+    run(Command::new(bin)
+        .arg("--dir")
+        .arg(&data)
+        .args(["--datasets", "3", "--mib", "1"])
+        .args(SMALL)
+        .arg("--json")
+        .arg(&claw));
+    // Second run reuses the files (and exercises --cold).
+    let out = Command::new(bin)
+        .arg("--dir")
+        .arg(&data)
+        .args(["--datasets", "3", "--mib", "1", "--cold"])
+        .args(SMALL)
+        .output()
+        .unwrap();
+    assert!(out.status.success());
+    assert!(String::from_utf8_lossy(&out.stderr).contains("reusing"));
+
+    let doc = results(&claw);
+    assert_eq!(doc["tool"], "clawhdf5");
+    // 2 layouts x 2 modes x 2 thread counts.
+    assert_eq!(doc["results"].as_array().unwrap().len(), 8);
+    for r in doc["results"].as_array().unwrap() {
+        assert!(r["mb_s"].as_f64().unwrap() > 0.0, "{r}");
+    }
+
+    if !python_available() {
+        assert!(
+            !interop_required(),
+            "CLAWHDF5_REQUIRE_INTEROP=1 but {} has no h5py",
+            python()
+        );
+        eprintln!("skipping the h5py half: no h5py in {}", python());
+        return;
+    }
+    let mut jsons = vec![claw];
+    for executor in ["threads", "processes"] {
+        let out = dir.path().join(format!("h5py-{executor}.json"));
+        run(Command::new(python())
+            .arg(scripts().join("concurrent_read_h5py.py"))
+            .arg("--dir")
+            .arg(&data)
+            .args(["--executor", executor])
+            .args(SMALL)
+            .arg("--json")
+            .arg(&out));
+        let doc = results(&out);
+        assert_eq!(doc["tool"], format!("h5py-{executor}"));
+        assert_eq!(doc["results"].as_array().unwrap().len(), 8);
+        jsons.push(out);
+    }
+    let table = run(Command::new(python())
+        .arg(scripts().join("compare_concurrent_read.py"))
+        .args(&jsons));
+    assert!(table.contains("| deflate | same | 2 |"), "{table}");
+    assert!(table.contains("clawhdf5 / h5py-processes"), "{table}");
+
+    // A different workload must not be compared.
+    let other = dir.path().join("other.json");
+    run(Command::new(python())
+        .arg(scripts().join("concurrent_read_h5py.py"))
+        .arg("--dir")
+        .arg(&data)
+        .args([
+            "--threads",
+            "1",
+            "--slabs",
+            "4",
+            "--reps",
+            "1",
+            "--slab",
+            "64",
+        ])
+        .arg("--json")
+        .arg(&other));
+    let out = Command::new(python())
+        .arg(scripts().join("compare_concurrent_read.py"))
+        .arg(&jsons[0])
+        .arg(&other)
+        .output()
+        .unwrap();
+    assert!(!out.status.success());
+    assert!(String::from_utf8_lossy(&out.stderr).contains("slabs"));
+}