perf(format): back large read buffers with transparent huge pages

A full read of a contiguous dataset is one memcpy from the mapped file,
yet ran at a quarter of h5py's speed on one thread: the fresh output Vec
took a page fault and a kernel page clear for every 4 KiB page written,
16384 per 64 MiB, costing several times the copy (the benchmark spent
6.2 s of 8 s in the kernel, 4.3M minor faults). numpy, so h5py, madvises
MADV_HUGEPAGE on allocations of 4 MiB or more; the typed readers' output,
the raw contiguous read and the chunk assembly buffer now do the same
(Linux only, libc as a Linux-only dependency; no-op otherwise).

New h5py comparison tests cover full and selection reads of contiguous
data for every 1-8-byte integer and float type, both byte orders, ranks
1-4, empty selections, and datasets past the 4 MiB threshold.

Co-Authored-By: Claude Opus 5.5 (1M context) <[email protected]>

CHANGELOG.md

@@ -2,6 +2,23 @@
 
 ## Unreleased
 
+### Contiguous read speed (2026-09-26)
+- **Large read buffers are backed by transparent huge pages.** A full read
+  of a contiguous dataset was one `memcpy` from the mapped file, yet ran at
+  a quarter of h5py's speed on one thread: the fresh output `Vec` took a
+  page fault (and a kernel page clear) for every 4 KiB page it was written
+  to, 16384 of them for 64 MiB, and those cost several times the copy.
+  numpy, and so h5py, asks for transparent huge pages on every allocation of
+  4 MiB or more; clawhdf5-format's read buffers now do too
+  (`madvise(MADV_HUGEPAGE)` on Linux, `libc` added as a Linux-only
+  dependency; a no-op elsewhere or when THP is disabled). It applies to the
+  typed readers' output (`read_f32`, `read_f64`, `read_i32`, `read_i64`,
+  `read_u64`, both byte orders), the raw contiguous read and the chunk
+  assembly buffer. Values are unchanged; new h5py comparison
+  `crates/clawhdf5/tests/contiguous_read_interop.rs` covers every 1-8-byte
+  integer and float type in both byte orders, ranks 1-4, and datasets past
+  the 4 MiB threshold.
+
 ### Plugin filters (2026-09-26)
 - **LZF, bitshuffle, bzip2 and Blosc read and write, in pure Rust.** Files
   written by h5py with `compression="lzf"`, or with hdf5plugin's

crates/clawhdf5-format/Cargo.toml

@@ -30,6 +30,10 @@ ruzstd = { version = "0.9", optional = true }
 bzip2 = { version = "0.6", optional = true }
 snap = { version = "1", optional = true }
 
+[target.'cfg(target_os = "linux")'.dependencies]
+# madvise(MADV_HUGEPAGE) for large read buffers (see src/bulk_alloc.rs).
+libc = { version = "0.2", default-features = false }
+
 [dev-dependencies]
 half = { workspace = true }
 serde_json = "1"

crates/clawhdf5-format/src/bulk_alloc.rs

@@ -0,0 +1,78 @@
+//! Large output buffers backed by transparent huge pages where the OS offers
+//! them.
+//!
+//! A fresh multi-megabyte `Vec` is mapped lazily by the kernel: the first
+//! write to each 4 KiB page takes a page fault, and the kernel zeroes the page
+//! before handing it over. For a 64 MiB read that is 16384 faults, and they
+//! cost far more than the copy that fills the buffer — single-threaded
+//! contiguous reads ran at about a quarter of h5py's speed because of them.
+//! numpy (so h5py) avoids this by asking for transparent huge pages
+//! (`madvise(MADV_HUGEPAGE)`) on every allocation of 4 MiB or more, which
+//! turns 512 faults into one; this module does the same.
+//!
+//! The advice only changes how the pages are backed, never their contents, so
+//! it is harmless when it cannot be honoured (THP disabled, not Linux, a
+//! region that is part of the heap): the buffer is then exactly what it would
+//! have been without it.
+
+#[cfg(not(feature = "std"))]
+use alloc::vec::Vec;
+
+/// Buffers smaller than this are left alone (numpy uses the same threshold).
+pub(crate) const HUGE_PAGE_THRESHOLD: usize = 4 << 20;
+
+/// Advise the kernel to back `[ptr, ptr + len)` with transparent huge pages,
+/// when `len` is large enough to benefit. Call it before the first write so
+/// the faults happen at huge-page granularity.
+#[inline]
+pub(crate) fn advise_huge_pages(ptr: *const u8, len: usize) {
+    #[cfg(target_os = "linux")]
+    if len >= HUGE_PAGE_THRESHOLD {
+        const PAGE: usize = 4096;
+        let start = (ptr as usize).next_multiple_of(PAGE);
+        let end = (ptr as usize + len) & !(PAGE - 1);
+        if end > start {
+            // SAFETY: `[start, end)` lies inside an allocation of `len` bytes
+            // at `ptr` that the caller owns, and is page aligned as madvise
+            // requires. MADV_HUGEPAGE does not change the memory's contents or
+            // validity; on failure (EINVAL when THP is compiled out, etc.) the
+            // region is simply left as it was, so the result is ignored.
+            unsafe {
+                libc::madvise(start as *mut libc::c_void, end - start, libc::MADV_HUGEPAGE);
+            }
+        }
+    }
+    #[cfg(not(target_os = "linux"))]
+    let _ = (ptr, len);
+}
+
+/// `Vec::with_capacity(count)` for a buffer about to be filled in bulk, with
+/// huge-page advice when it is large (see the module docs).
+#[inline]
+pub(crate) fn vec_for_bulk<T>(count: usize) -> Vec<T> {
+    let v: Vec<T> = Vec::with_capacity(count);
+    advise_huge_pages(
+        v.as_ptr().cast::<u8>(),
+        v.capacity().saturating_mul(core::mem::size_of::<T>()),
+    );
+    v
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn bulk_vec_is_an_ordinary_vec() {
+        for count in [0usize, 1, 1000, HUGE_PAGE_THRESHOLD / 4 + 3] {
+            let mut v: Vec<u32> = vec_for_bulk(count);
+            assert!(v.capacity() >= count);
+            v.extend((0..count as u32).map(|i| i.wrapping_mul(2654435761)));
+            assert!(
+                v.iter()
+                    .enumerate()
+                    .all(|(i, &x)| x == (i as u32).wrapping_mul(2654435761))
+            );
+        }
+    }
+}

crates/clawhdf5-format/src/chunked_read.rs

@@ -277,6 +277,8 @@ pub(crate) fn alloc_output(len: usize) -> Result<Vec<u8>, FormatError> {
     if ptr.is_null() {
         return Err(failed());
     }
+    // Before anything writes to it, so a large buffer faults in huge pages.
+    crate::bulk_alloc::advise_huge_pages(ptr, len);
     // SAFETY: `ptr` came from the global allocator with the layout of
     // `[u8; len]`, which is exactly what `Vec<u8>` with capacity `len` frees;
     // all `len` bytes are initialised (zero).

crates/clawhdf5-format/src/data_read.rs

@@ -180,7 +180,9 @@ fn read_raw_data_full_impl(
                 });
             }
             ensure_len(file_data, addr, sz)?;
-            Ok(file_data[addr..addr + sz].to_vec())
+            let mut out = crate::bulk_alloc::vec_for_bulk(sz);
+            out.extend_from_slice(&file_data[addr..addr + sz]);
+            Ok(out)
         }
         DataLayout::Chunked { .. } => read_chunked_data(
             file_data,
@@ -765,7 +767,7 @@ fn get_size(dt: &Datatype) -> usize {
 fn native_le_to_vec<T: Copy>(raw: &[u8], count: usize) -> Vec<T> {
     let bytes = count * core::mem::size_of::<T>();
     debug_assert!(bytes <= raw.len());
-    let mut result: Vec<T> = Vec::with_capacity(count);
+    let mut result: Vec<T> = crate::bulk_alloc::vec_for_bulk(count);
     // SAFETY: `result` has capacity for `count` values of `T`, i.e. `bytes`
     // bytes; `raw` holds at least `bytes` bytes (callers derive `count` from
     // `raw.len() / size_of::<T>()`); the regions cannot overlap because
@@ -803,7 +805,7 @@ pub fn read_as_f64(raw: &[u8], datatype: &Datatype) -> Result<Vec<f64>, FormatEr
     }
 
     let order = get_byte_order(datatype);
-    let mut result = Vec::with_capacity(count);
+    let mut result = crate::bulk_alloc::vec_for_bulk(count);
     if let Datatype::FloatingPoint { .. } = datatype {
         let format = FloatFormat::of(datatype)?;
         for chunk in raw.chunks_exact(elem_size) {
@@ -955,7 +957,7 @@ pub fn read_as_i64(raw: &[u8], datatype: &Datatype) -> Result<Vec<i64>, FormatEr
     }
 
     let order = get_byte_order(datatype);
-    let mut result = Vec::with_capacity(count);
+    let mut result = crate::bulk_alloc::vec_for_bulk(count);
     for i in 0..count {
         let chunk = &raw[i * elem_size..(i + 1) * elem_size];
         result.push(decode_scalar(chunk, datatype, &order)?.to_i64());
@@ -985,7 +987,7 @@ pub fn read_as_u64(raw: &[u8], datatype: &Datatype) -> Result<Vec<u64>, FormatEr
     }
     let count = raw.len() / elem_size;
     let order = get_byte_order(datatype);
-    let mut result = Vec::with_capacity(count);
+    let mut result = crate::bulk_alloc::vec_for_bulk(count);
     for i in 0..count {
         let chunk = &raw[i * elem_size..(i + 1) * elem_size];
         result.push(decode_scalar(chunk, datatype, &order)?.to_u64());
@@ -1018,14 +1020,17 @@ pub fn read_as_f32(raw: &[u8], datatype: &Datatype) -> Result<Vec<f32>, FormatEr
     // Little-endian IEEE half precision (numpy float16): widen directly.
     if is_native_le_float(datatype, FloatFormat::Half) {
         let (halves, _) = raw[..count * 2].as_chunks::<2>();
-        return Ok(halves
+        let mut result = crate::bulk_alloc::vec_for_bulk(count);
-            .iter()
+        result.extend(
-            .map(|&b| f16_bits_to_f32(u16::from_le_bytes(b)))
+            halves
-            .collect());
+                .iter()
+                .map(|&b| f16_bits_to_f32(u16::from_le_bytes(b))),
+        );
+        return Ok(result);
     }
 
     let order = get_byte_order(datatype);
-    let mut result = Vec::with_capacity(count);
+    let mut result = crate::bulk_alloc::vec_for_bulk(count);
     if let Datatype::FloatingPoint { .. } = datatype {
         let format = FloatFormat::of(datatype)?;
         for chunk in raw.chunks_exact(elem_size) {
@@ -1114,7 +1119,7 @@ pub fn read_as_i32(raw: &[u8], datatype: &Datatype) -> Result<Vec<i32>, FormatEr
     }
 
     let order = get_byte_order(datatype);
-    let mut result = Vec::with_capacity(count);
+    let mut result = crate::bulk_alloc::vec_for_bulk(count);
     for i in 0..count {
         let chunk = &raw[i * elem_size..(i + 1) * elem_size];
         result.push(decode_scalar(chunk, datatype, &order)?.to_i32());

crates/clawhdf5-format/src/lib.rs

@@ -61,6 +61,7 @@ pub mod attribute;
 pub mod attribute_info;
 pub mod btree_v1;
 pub mod btree_v2;
+mod bulk_alloc;
 pub mod checksum;
 pub mod chunk_cache;
 mod chunk_grid;

crates/clawhdf5/tests/contiguous_read_interop.rs

@@ -0,0 +1,403 @@
+//! Reads of contiguous datasets — full reads and hyperslab/point selections,
+//! through every typed reader — checked against h5py/libhdf5 for every
+//! integer and float width, both byte orders, ranks 1 to 4, and datasets
+//! larger than the huge-page threshold (4 MiB) the read buffers use.
+//!
+//! h5py writes the file and, for each selection, reads it with libhdf5's own
+//! hyperslab/point selection (`select_hyperslab` with stride and block,
+//! `select_elements`) and saves the raw bytes it gets back; the byte-level
+//! [`Dataset::read_selection`] must return exactly those bytes, and the typed
+//! readers the same values. Skipped when python3 with h5py is unavailable,
+//! unless `CLAWHDF5_REQUIRE_INTEROP=1`.
+
+use std::path::Path;
+use std::process::Command;
+
+use clawhdf5::File;
+use clawhdf5_format::selection::Selection;
+
+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)
+}
+
+macro_rules! skip_if_no_python {
+    () => {
+        if !python_available() {
+            assert!(
+                !interop_required(),
+                "CLAWHDF5_REQUIRE_INTEROP=1 but python3 with h5py is not available"
+            );
+            eprintln!("SKIP: python3 with h5py not available");
+            return;
+        }
+    };
+}
+
+fn run_python(script: &str) {
+    let output = Command::new(python())
+        .args(["-c", script])
+        .output()
+        .expect("failed to run python");
+    assert!(
+        output.status.success(),
+        "python failed:\n{}",
+        String::from_utf8_lossy(&output.stderr)
+    );
+}
+
+/// numpy type codes, with the modulus of the value pattern: every value is an
+/// integer exactly representable in the type and in every typed reader's
+/// output (f16 is exact below 2048, f32 below 2^24).
+const DTYPES: [(&str, i64, bool); 11] = [
+    ("i1", 201, true),
+    ("u1", 251, false),
+    ("i2", 2039, true),
+    ("u2", 2039, false),
+    ("i4", 1_000_003, true),
+    ("u4", 1_000_003, false),
+    ("i8", 1_000_003, true),
+    ("u8", 1_000_003, false),
+    ("f2", 2039, true),
+    ("f4", 1_000_003, true),
+    ("f8", 1_000_003, true),
+];
+
+/// Value of element `i` of a dataset of type `code` (the same formula as the
+/// Python side): a permutation-revealing pattern, centred on 0 when signed.
+fn value(code: &str, i: u64) -> i64 {
+    let (_, m, signed) = DTYPES.iter().find(|d| d.0 == code).unwrap();
+    let v = ((i as i128 * 7919) % *m as i128) as i64;
+    if *signed { v - m / 2 } else { v }
+}
+
+const SHAPES: [&[u64]; 4] = [&[1000], &[37, 53], &[7, 11, 13], &[3, 5, 7, 9]];
+
+/// Datasets past the 4 MiB huge-page threshold, as (type, shape).
+const BIG: [(&str, [u64; 2]); 4] = [
+    ("f4", [1100, 1024]),
+    ("i4", [1100, 1024]),
+    ("f8", [600, 1024]),
+    ("i8", [600, 1024]),
+];
+
+fn datasets() -> Vec<(String, String, Vec<u64>)> {
+    let mut out = Vec::new();
+    for (code, _, _) in DTYPES {
+        for (tag, _) in [("le", '<'), ("be", '>')] {
+            for shape in SHAPES {
+                out.push((
+                    format!("{code}{tag}_r{}", shape.len()),
+                    code.to_string(),
+                    shape.to_vec(),
+                ));
+            }
+        }
+    }
+    for (code, shape) in BIG {
+        for tag in ["le", "be"] {
+            out.push((format!("{code}{tag}_big"), code.to_string(), shape.to_vec()));
+        }
+    }
+    out
+}
+
+fn write_file(path: &Path) {
+    let script = format!(
+        r#"
+import h5py, numpy as np
+M = {{'i1': 201, 'u1': 251, 'i2': 2039, 'u2': 2039, 'i4': 1000003, 'u4': 1000003,
+     'i8': 1000003, 'u8': 1000003, 'f2': 2039, 'f4': 1000003, 'f8': 1000003}}
+def values(code, n):
+    v = (np.arange(n, dtype=np.int64) * 7919) % M[code]
+    if code[0] != 'u':
+        v -= M[code] // 2
+    return v
+shapes = [(1000,), (37, 53), (7, 11, 13), (3, 5, 7, 9)]
+big = [('f4', (1100, 1024)), ('i4', (1100, 1024)), ('f8', (600, 1024)), ('i8', (600, 1024))]
+with h5py.File("{path}", "w") as f:
+    for code in M:
+        for tag, e in (('le', '<'), ('be', '>')):
+            for shape in shapes:
+                n = int(np.prod(shape))
+                f.create_dataset(f"{{code}}{{tag}}_r{{len(shape)}}",
+                                 data=values(code, n).astype(e + code).reshape(shape))
+    for code, shape in big:
+        for tag, e in (('le', '<'), ('be', '>')):
+            n = int(np.prod(shape))
+            f.create_dataset(f"{{code}}{{tag}}_big",
+                             data=values(code, n).astype(e + code).reshape(shape))
+"#,
+        path = path.display()
+    );
+    run_python(&script);
+}
+
+#[test]
+fn full_reads_match_h5py_for_every_type_order_and_size() {
+    skip_if_no_python!();
+    let dir = tempfile::tempdir().unwrap();
+    let path = dir.path().join("contig.h5");
+    write_file(&path);
+    let file = File::open(&path).unwrap();
+    for (name, code, shape) in datasets() {
+        let ds = file.dataset(&name).unwrap();
+        assert!(ds.read_raw_ref().unwrap().is_some(), "{name} is contiguous");
+        let n: u64 = shape.iter().product();
+        let want: Vec<i64> = (0..n).map(|i| value(&code, i)).collect();
+        assert_eq!(
+            ds.read_f64().unwrap(),
+            want.iter().map(|&v| v as f64).collect::<Vec<_>>(),
+            "{name} read_f64"
+        );
+        assert_eq!(
+            ds.read_f32().unwrap(),
+            want.iter().map(|&v| v as f32).collect::<Vec<_>>(),
+            "{name} read_f32"
+        );
+        assert_eq!(ds.read_i64().unwrap(), want, "{name} read_i64");
+        assert_eq!(
+            ds.read_i32().unwrap(),
+            want.iter().map(|&v| v as i32).collect::<Vec<_>>(),
+            "{name} read_i32"
+        );
+        // libhdf5 saturates negative values to 0 when reading as unsigned.
+        assert_eq!(
+            ds.read_u64().unwrap(),
+            want.iter().map(|&v| v.max(0) as u64).collect::<Vec<_>>(),
+            "{name} read_u64"
+        );
+    }
+}
+
+struct Rng(u64);
+impl Rng {
+    fn next(&mut self) -> u64 {
+        self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15);
+        let mut z = self.0;
+        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
+        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
+        z ^ (z >> 31)
+    }
+    fn below(&mut self, n: u64) -> u64 {
+        self.next() % n.max(1)
+    }
+}
+
+/// A hyperslab from per-dimension `(start, stride, count, block)`.
+fn slab(dims: &[(u64, u64, u64, u64)]) -> Selection {
+    Selection::Hyperslab {
+        start: dims.iter().map(|d| d.0).collect(),
+        stride: dims.iter().map(|d| d.1).collect(),
+        count: dims.iter().map(|d| d.2).collect(),
+        block: dims.iter().map(|d| d.3).collect(),
+    }
+}
+
+/// Selections of every shape the read paths distinguish, all valid for `dims`.
+fn selections(rng: &mut Rng, dims: &[u64]) -> Vec<Selection> {
+    let mut out = Vec::new();
+    // Unit-stride box.
+    out.push(slab(
+        &dims
+            .iter()
+            .map(|&n| {
+                let c = 1 + rng.below(n);
+                (rng.below(n - c + 1), 1, c, 1)
+            })
+            .collect::<Vec<_>>(),
+    ));
+    // Strided (block 1), blocked (stride > block), and adjacent blocks
+    // (stride == block, which reads like a box): (block, stride - block).
+    for (block, gap) in [(1, 1), (2, 1), (2, 0)] {
+        out.push(slab(
+            &dims
+                .iter()
+                .map(|&n| {
+                    let b = (block + rng.below(2)).min(n);
+                    let st = b + gap + rng.below(2) * gap;
+                    let s = rng.below(n - b + 1);
+                    let c = 1 + rng.below((n - s - b) / st + 1);
+                    (s, st, c, b)
+                })
+                .collect::<Vec<_>>(),
+        ));
+    }
+    // Whole inner rows (one run across rows), and the whole dataset.
+    let r0 = rng.below(dims[0]);
+    let mut rows = vec![(r0, 1, 1 + rng.below(dims[0] - r0), 1)];
+    rows.extend(dims[1..].iter().map(|&n| (0, 1, n, 1)));
+    out.push(slab(&rows));
+    out.push(slab(
+        &dims.iter().map(|&n| (0, 1, n, 1)).collect::<Vec<_>>(),
+    ));
+    // One element.
+    out.push(slab(
+        &dims
+            .iter()
+            .map(|&n| (rng.below(n), 1, 1, 1))
+            .collect::<Vec<_>>(),
+    ));
+    // Distinct points in no particular order.
+    let mut points: Vec<Vec<u64>> = Vec::new();
+    for _ in 0..1 + rng.below(15) {
+        let p: Vec<u64> = dims.iter().map(|&n| rng.below(n)).collect();
+        if !points.contains(&p) {
+            points.push(p);
+        }
+    }
+    out.push(Selection::Points(points));
+    out
+}
+
+fn join(v: &[u64]) -> String {
+    v.iter().map(u64::to_string).collect::<Vec<_>>().join(",")
+}
+
+/// One element of type `code`, given as its raw file-order bytes, as an
+/// integer (every value in these files is one).
+fn decode(code: &str, big_endian: bool, bytes: &[u8]) -> i64 {
+    let mut b = bytes.to_vec();
+    if big_endian {
+        b.reverse();
+    }
+    let mut w = [0u8; 8];
+    w[..b.len()].copy_from_slice(&b);
+    let u = u64::from_le_bytes(w);
+    match code {
+        "i1" => u as u8 as i8 as i64,
+        "i2" => u as u16 as i16 as i64,
+        "i4" => u as u32 as i32 as i64,
+        "i8" => u as i64,
+        "u1" | "u2" | "u4" | "u8" => u as i64,
+        "f2" => clawhdf5_format::float16::f16_bits_to_f32(u as u16) as i64,
+        "f4" => f32::from_bits(u as u32) as i64,
+        "f8" => f64::from_bits(u) as i64,
+        _ => unreachable!(),
+    }
+}
+
+#[test]
+fn selection_reads_match_h5py_for_every_type_order_and_rank() {
+    skip_if_no_python!();
+    let dir = tempfile::tempdir().unwrap();
+    let path = dir.path().join("contig.h5");
+    write_file(&path);
+
+    let mut rng = Rng(2026);
+    let mut cases: Vec<(String, String, Selection)> = Vec::new();
+    for (name, code, shape) in datasets() {
+        for sel in selections(&mut rng, &shape) {
+            cases.push((name.clone(), code.clone(), sel));
+        }
+        // Empty: a zero count, and Selection::None.
+        let mut empty = shape.iter().map(|&n| (0, 1, n, 1)).collect::<Vec<_>>();
+        empty[shape.len() - 1].2 = 0;
+        cases.push((name.clone(), code.clone(), slab(&empty)));
+        cases.push((name, code, Selection::None));
+    }
+    let mut spec = String::new();
+    for (k, (name, _, sel)) in cases.iter().enumerate() {
+        let line = match sel {
+            Selection::Hyperslab { count, .. } if count.contains(&0) => "N".to_string(),
+            Selection::Hyperslab {
+                start,
+                stride,
+                count,
+                block,
+            } => format!(
+                "H {};{};{};{}",
+                join(start),
+                join(stride),
+                join(count),
+                join(block)
+            ),
+            Selection::Points(points) => format!(
+                "P {}",
+                points.iter().map(|p| join(p)).collect::<Vec<_>>().join(";")
+            ),
+            Selection::None => "N".to_string(),
+            Selection::All => unreachable!(),
+        };
+        spec.push_str(&format!("{k} {name} {line}\n"));
+    }
+    let spec_path = dir.path().join("cases.txt");
+    std::fs::write(&spec_path, spec).unwrap();
+    run_python(&format!(
+        r#"
+import h5py, numpy as np
+with h5py.File("{path}", "r") as f:
+    for line in open("{spec}"):
+        k, name, kind, *rest = line.split()
+        d = f[name]
+        space = d.id.get_space()
+        if kind == 'H':
+            start, stride, count, block = (tuple(int(x) for x in part.split(','))
+                                           for part in rest[0].split(';'))
+            space.select_hyperslab(start, count, stride, block)
+        elif kind == 'P':
+            pts = np.array([[int(x) for x in p.split(',')] for p in rest[0].split(';')],
+                           dtype=np.uint64)
+            space.select_elements(pts)
+        else:
+            space.select_none()
+        n = space.get_select_npoints()
+        out = np.empty(n, dtype=d.dtype)
+        if n:
+            d.id.read(h5py.h5s.create_simple((n,)), space, out)
+        open("{dir}/sel_" + k + ".bin", "wb").write(out.tobytes())
+"#,
+        path = path.display(),
+        spec = spec_path.display(),
+        dir = dir.path().display(),
+    ));
+
+    let file = File::open(&path).unwrap();
+    for (k, (name, code, sel)) in cases.iter().enumerate() {
+        let ds = file.dataset(name).unwrap();
+        let want_bytes = std::fs::read(dir.path().join(format!("sel_{k}.bin"))).unwrap();
+        let got_bytes = ds.read_selection(sel).unwrap();
+        assert!(
+            got_bytes == want_bytes,
+            "{name} {sel:?}: raw bytes differ from libhdf5's ({} vs {} bytes)",
+            got_bytes.len(),
+            want_bytes.len()
+        );
+        let size = ds.raw_datatype().unwrap().type_size() as usize;
+        let want: Vec<i64> = want_bytes
+            .chunks_exact(size)
+            .map(|e| decode(code, name.contains("be_"), e))
+            .collect();
+        assert_eq!(
+            ds.read_f64_selection(sel).unwrap(),
+            want.iter().map(|&v| v as f64).collect::<Vec<_>>(),
+            "{name} {sel:?} as f64"
+        );
+        assert_eq!(
+            ds.read_f32_selection(sel).unwrap(),
+            want.iter().map(|&v| v as f32).collect::<Vec<_>>(),
+            "{name} {sel:?} as f32"
+        );
+        assert_eq!(
+            ds.read_i64_selection(sel).unwrap(),
+            want,
+            "{name} {sel:?} as i64"
+        );
+        assert_eq!(
+            ds.read_i32_selection(sel).unwrap(),
+            want.iter().map(|&v| v as i32).collect::<Vec<_>>(),
+            "{name} {sel:?} as i32"
+        );
+    }
+}