From 78c769f179dd26b85a675e8de436d9e83138b5ef Mon Sep 17 00:00:00 2001 From: osobh Date: Sat, 26 Sep 2026 08:12:55 -0500 Subject: [PATCH 1/3] 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) --- CHANGELOG.md | 17 + crates/clawhdf5-format/Cargo.toml | 4 + crates/clawhdf5-format/src/bulk_alloc.rs | 78 ++++ crates/clawhdf5-format/src/chunked_read.rs | 2 + crates/clawhdf5-format/src/data_read.rs | 27 +- crates/clawhdf5-format/src/lib.rs | 1 + .../clawhdf5/tests/contiguous_read_interop.rs | 403 ++++++++++++++++++ 7 files changed, 521 insertions(+), 11 deletions(-) create mode 100644 crates/clawhdf5-format/src/bulk_alloc.rs create mode 100644 crates/clawhdf5/tests/contiguous_read_interop.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index ef4bd2f..99f86f2 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/crates/clawhdf5-format/Cargo.toml b/crates/clawhdf5-format/Cargo.toml index 75fd233..1d1998c 100644 --- a/crates/clawhdf5-format/Cargo.toml +++ b/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" diff --git a/crates/clawhdf5-format/src/bulk_alloc.rs b/crates/clawhdf5-format/src/bulk_alloc.rs new file mode 100644 index 0000000..22a8bf2 --- /dev/null +++ b/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(count: usize) -> Vec { + let v: Vec = Vec::with_capacity(count); + advise_huge_pages( + v.as_ptr().cast::(), + v.capacity().saturating_mul(core::mem::size_of::()), + ); + 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 = 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)) + ); + } + } +} diff --git a/crates/clawhdf5-format/src/chunked_read.rs b/crates/clawhdf5-format/src/chunked_read.rs index cbbaa69..0baec25 100644 --- a/crates/clawhdf5-format/src/chunked_read.rs +++ b/crates/clawhdf5-format/src/chunked_read.rs @@ -277,6 +277,8 @@ pub(crate) fn alloc_output(len: usize) -> Result, 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` with capacity `len` frees; // all `len` bytes are initialised (zero). diff --git a/crates/clawhdf5-format/src/data_read.rs b/crates/clawhdf5-format/src/data_read.rs index b8f2ffe..1c53418 100644 --- a/crates/clawhdf5-format/src/data_read.rs +++ b/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(raw: &[u8], count: usize) -> Vec { let bytes = count * core::mem::size_of::(); debug_assert!(bytes <= raw.len()); - let mut result: Vec = Vec::with_capacity(count); + let mut result: Vec = 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::()`); the regions cannot overlap because @@ -803,7 +805,7 @@ pub fn read_as_f64(raw: &[u8], datatype: &Datatype) -> Result, 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, 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, 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, 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 - .iter() - .map(|&b| f16_bits_to_f32(u16::from_le_bytes(b))) - .collect()); + let mut result = crate::bulk_alloc::vec_for_bulk(count); + result.extend( + halves + .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, 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()); diff --git a/crates/clawhdf5-format/src/lib.rs b/crates/clawhdf5-format/src/lib.rs index 905ce84..05ebd46 100644 --- a/crates/clawhdf5-format/src/lib.rs +++ b/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; diff --git a/crates/clawhdf5/tests/contiguous_read_interop.rs b/crates/clawhdf5/tests/contiguous_read_interop.rs new file mode 100644 index 0000000..188856e --- /dev/null +++ b/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)> { + 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 = (0..n).map(|i| value(&code, i)).collect(); + assert_eq!( + ds.read_f64().unwrap(), + want.iter().map(|&v| v as f64).collect::>(), + "{name} read_f64" + ); + assert_eq!( + ds.read_f32().unwrap(), + want.iter().map(|&v| v as f32).collect::>(), + "{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::>(), + "{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::>(), + "{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 { + 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::>(), + )); + // 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::>(), + )); + } + // 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::>(), + )); + // One element. + out.push(slab( + &dims + .iter() + .map(|&n| (rng.below(n), 1, 1, 1)) + .collect::>(), + )); + // Distinct points in no particular order. + let mut points: Vec> = Vec::new(); + for _ in 0..1 + rng.below(15) { + let p: Vec = 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::>().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::>(); + 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::>().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 = 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::>(), + "{name} {sel:?} as f64" + ); + assert_eq!( + ds.read_f32_selection(sel).unwrap(), + want.iter().map(|&v| v as f32).collect::>(), + "{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::>(), + "{name} {sel:?} as i32" + ); + } +} From 2bc4cb46a67f0abfb56b2e6e14982799efc09118 Mon Sep 17 00:00:00 2001 From: osobh Date: Sat, 26 Sep 2026 08:21:47 -0500 Subject: [PATCH 2/3] perf: copy contiguous hyperslab and point reads run by run A 256 x 256 hyperslab of a contiguous f32 dataset read at an eighth of h5py's speed: partial_read copied the bounding box out of the file, the extractor then walked it element by element (a recursive call and two bounds checks per element) into a second buffer, and read_f32_selection converted that into a third. Selections of contiguous data are now copied straight from the file, one memcpy per run of elements contiguous in the file (gather.rs: a block along the last dimension, touching blocks as one range, whole rows merged), with no zero-filled intermediate and no full copy for large selections. The typed selection readers copy into their Vec directly when the dataset stores T natively (new data_read::read_selection_native and sealed NativeElement trait, which the read_as_* fast paths now share; read_as_u64 gains one) and convert as before otherwise. The general extractor used by the chunked paths runs on the same run walker, keeping its old handling of unvalidated selections. Checked against h5py (contiguous_read_interop.rs) for strided, blocked, adjacent-block and whole-row hyperslabs, points and empty selections of every 1-8-byte type in both byte orders, ranks 1-4. Also keeps the huge-page threshold constant out of no_std builds, where it was unused. Co-Authored-By: Claude Opus 5.5 (1M context) --- CHANGELOG.md | 21 ++ crates/clawhdf5-format/src/bulk_alloc.rs | 1 + crates/clawhdf5-format/src/data_read.rs | 279 +++++++++-------- crates/clawhdf5-format/src/gather.rs | 343 +++++++++++++++++++++ crates/clawhdf5-format/src/lib.rs | 1 + crates/clawhdf5-format/src/partial_read.rs | 61 ++-- crates/clawhdf5/src/reader.rs | 41 ++- docs/known-issues.md | 9 +- 8 files changed, 586 insertions(+), 170 deletions(-) create mode 100644 crates/clawhdf5-format/src/gather.rs diff --git a/CHANGELOG.md b/CHANGELOG.md index 99f86f2..bb70d01 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -18,6 +18,27 @@ `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. +- **Hyperslab and point reads of contiguous data copy runs, not elements.** + A 256 x 256 hyperslab of a contiguous `f32` dataset read at an eighth of + h5py's speed: the selection's bounding box was copied out of the file, + then walked element by element (a recursive call and two bounds checks per + element) into a second buffer, which `read_f32_selection` converted into + a third. Selections of contiguous data are now copied straight from the + file, one `memcpy` per run of elements that is contiguous in the file + (a block along the last dimension, blocks that touch, and whole rows when + the inner dimensions are selected in full, merged), with no zero-filled + intermediate; a selection covering most of the dataset no longer makes a + full copy first. The typed selection readers (`read_f32_selection`, + `read_f64_selection`, `read_i32_selection`, `read_i64_selection`) copy + directly into their output when the dataset stores that type natively, + and convert as before otherwise (big-endian, other widths). The chunked + paths use the same run-based extraction. New public + `clawhdf5_format::data_read::read_selection_native` and the sealed + `NativeElement` trait (also used by the `read_as_*` fast paths, which + gained one for native `u64`). Values are unchanged: checked against h5py + by `contiguous_read_interop.rs` (strided, blocked, adjacent-block and + whole-row hyperslabs, points, empty selections; every type, both byte + orders, ranks 1-4). ### Plugin filters (2026-09-26) - **LZF, bitshuffle, bzip2 and Blosc read and write, in pure Rust.** Files diff --git a/crates/clawhdf5-format/src/bulk_alloc.rs b/crates/clawhdf5-format/src/bulk_alloc.rs index 22a8bf2..085f125 100644 --- a/crates/clawhdf5-format/src/bulk_alloc.rs +++ b/crates/clawhdf5-format/src/bulk_alloc.rs @@ -19,6 +19,7 @@ use alloc::vec::Vec; /// Buffers smaller than this are left alone (numpy uses the same threshold). +#[cfg(any(target_os = "linux", test))] pub(crate) const HUGE_PAGE_THRESHOLD: usize = 4 << 20; /// Advise the kernel to back `[ptr, ptr + len)` with transparent huge pages, diff --git a/crates/clawhdf5-format/src/data_read.rs b/crates/clawhdf5-format/src/data_read.rs index 1c53418..633ecd7 100644 --- a/crates/clawhdf5-format/src/data_read.rs +++ b/crates/clawhdf5-format/src/data_read.rs @@ -531,6 +531,11 @@ pub fn extract_selection_from_buffer( block, } => { let rank = dims.len(); + if [start.len(), stride.len(), count.len(), block.len()] != [rank; 4] { + return Err(FormatError::SelectionOutOfBounds(format!( + "hyperslab rank does not match dataset rank {rank}" + ))); + } let output_elements = count .iter() .zip(block.iter()) @@ -540,96 +545,40 @@ pub fn extract_selection_from_buffer( crate::chunked_read::checked_byte_len(output_elements, elem_size)?, )?; - // Compute dataset strides (row-major) - let mut ds_strides = vec![1usize; rank]; - for i in (0..rank.saturating_sub(1)).rev() { - ds_strides[i] = ds_strides[i + 1] * dims[i + 1] as usize; - } - - // Compute output shape and strides - let output_dims: Vec = count - .iter() - .zip(block.iter()) - .map(|(&c, &b)| (c * b) as usize) - .collect(); - let mut out_strides = vec![1usize; rank]; - for i in (0..rank.saturating_sub(1)).rev() { - out_strides[i] = out_strides[i + 1] * output_dims[i + 1]; - } - - // Iterate over all selected elements - // For each block in the hyperslab, copy the elements - let mut out_linear = 0usize; - let _block_coords = vec![0u64; rank]; - - #[allow(clippy::too_many_arguments)] - fn iterate_hyperslab( - d: usize, - rank: usize, - start: &[u64], - stride: &[u64], - count: &[u64], - block: &[u64], - dims: &[u64], - ds_strides: &[usize], - elem_size: usize, - full_data: &[u8], - output: &mut [u8], - out_linear: &mut usize, - current_ds_offset: usize, - ) { - if d == rank { - // Copy one element - let src = current_ds_offset * elem_size; - let dst = *out_linear * elem_size; - if src + elem_size <= full_data.len() && dst + elem_size <= output.len() { - output[dst..dst + elem_size] - .copy_from_slice(&full_data[src..src + elem_size]); - } - *out_linear += 1; - return; - } - - for bi in 0..count[d] { - let block_start = start[d] + bi * stride[d]; - for bj in 0..block[d] { - let coord = block_start + bj; - if coord < dims[d] { - iterate_hyperslab( - d + 1, - rank, - start, - stride, - count, - block, - dims, - ds_strides, - elem_size, - full_data, - output, - out_linear, - current_ds_offset + coord as usize * ds_strides[d], - ); + // One copy per run of elements contiguous in `full_data` + // (`gather`'s runs). Coordinates past the extent are skipped and + // runs past the end of `full_data` left as zeros, element by + // element, as this extractor always did; validated selections + // never hit either. + let mut out_at = 0usize; + crate::gather::hyperslab_runs(dims, start, stride, count, block, |first, n| { + let big = |v: u64| usize::try_from(v).unwrap_or(usize::MAX); + let (first, n) = (big(first), big(n)); + let len = n.saturating_mul(elem_size); + let src = first.saturating_mul(elem_size); + let out_end = out_at.saturating_add(len); + if let (Some(from), Some(to)) = ( + full_data.get(src..src.saturating_add(len)), + output.get_mut(out_at..out_end), + ) { + to.copy_from_slice(from); + } else { + for k in 0..n { + let s = first.saturating_add(k).saturating_mul(elem_size); + let o = out_at.saturating_add(k.saturating_mul(elem_size)); + if o >= output.len() { + break; + } + if let (Some(from), Some(to)) = ( + full_data.get(s..s.saturating_add(elem_size)), + output.get_mut(o..o.saturating_add(elem_size)), + ) { + to.copy_from_slice(from); } } } - } - - iterate_hyperslab( - 0, - rank, - start, - stride, - count, - block, - dims, - &ds_strides, - elem_size, - full_data, - &mut output, - &mut out_linear, - 0, - ); + out_at = out_end; + }); Ok(output) } @@ -757,22 +706,76 @@ fn get_size(dt: &Datatype) -> usize { dt.type_size() as usize } -/// Reinterpret little-endian bytes as `count` native values of `T` on a -/// little-endian target, in one copy. +mod sealed { + pub trait Sealed {} +} + +/// A numeric type whose values can be copied straight out of a dataset's +/// bytes when the dataset stores exactly that type in the target's byte +/// order: `u8`, `i32`, `i64`, `u64`, `f32` and `f64`. +/// +/// # Safety +/// +/// Implementors have no padding and no invalid bit patterns, so a buffer of +/// them may be filled by copying bytes. The trait is sealed. +pub unsafe trait NativeElement: sealed::Sealed + Copy + 'static { + /// Whether `datatype`'s stored bytes are this type's native in-memory + /// representation (same size, byte order, signedness, full precision, + /// IEEE layout), so reading needs a copy and no conversion. + fn is_native(datatype: &Datatype) -> bool; +} + +/// A full-width fixed-point type of `size` bytes and the given signedness in +/// the target's byte order. +fn is_native_int(datatype: &Datatype, size: u32, want_signed: bool) -> bool { + let order = if cfg!(target_endian = "little") { + DatatypeByteOrder::LittleEndian + } else { + DatatypeByteOrder::BigEndian + }; + matches!( + datatype, + Datatype::FixedPoint { size: s, signed, byte_order, .. } + if *s == size && *signed == want_signed && (size == 1 || *byte_order == order) + ) && is_full_width(datatype) +} + +macro_rules! native_element { + ($($t:ty => |$dt:ident| $check:expr;)*) => {$( + impl sealed::Sealed for $t {} + // SAFETY: a primitive integer or float: no padding, and every bit + // pattern is a valid value. + unsafe impl NativeElement for $t { + fn is_native($dt: &Datatype) -> bool { + $check + } + } + )*}; +} + +native_element! { + u8 => |dt| is_native_int(dt, 1, false); + i32 => |dt| is_native_int(dt, 4, true); + i64 => |dt| is_native_int(dt, 8, true); + u64 => |dt| is_native_int(dt, 8, false); + f32 => |dt| cfg!(target_endian = "little") && is_native_le_float(dt, FloatFormat::Single); + f64 => |dt| cfg!(target_endian = "little") && is_native_le_float(dt, FloatFormat::Double); +} + +/// Copy `count` values of `T` out of `raw`, which holds them in `T`'s native +/// representation (see [`NativeElement::is_native`]), in one copy. /// /// The buffer is allocated uninitialised and filled by the copy. It used to be /// `vec![0; count]` first, which for a large dataset meant writing every page /// twice (zero it, then overwrite it) — about as expensive as the copy itself. -#[cfg(target_endian = "little")] -fn native_le_to_vec(raw: &[u8], count: usize) -> Vec { +fn native_to_vec(raw: &[u8], count: usize) -> Vec { let bytes = count * core::mem::size_of::(); - debug_assert!(bytes <= raw.len()); + assert!(bytes <= raw.len(), "native_to_vec: source too short"); let mut result: Vec = 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::()`); the regions cannot overlap because - // `result` was just allocated. Every `T` used here (f32/f64/i32/i64) is - // valid for any bit pattern, so after the copy all `count` values are + // bytes; `raw` holds at least `bytes` bytes (asserted); the regions + // cannot overlap because `result` was just allocated. `T: NativeElement` + // is valid for any bit pattern, so after the copy all `count` values are // initialised and `set_len` is sound. unsafe { core::ptr::copy_nonoverlapping(raw.as_ptr(), result.as_mut_ptr().cast::(), bytes); @@ -781,6 +784,44 @@ fn native_le_to_vec(raw: &[u8], count: usize) -> Vec { result } +/// Read `selection` of a dataset whose raw bytes (all of them, row-major, of +/// shape `dims`) are `raw` — typically a contiguous dataset's bytes borrowed +/// from the file — straight into a `Vec`, copying each contiguous run of +/// selected elements once. +/// +/// Returns `Ok(None)` when `datatype` is not `T`'s native representation +/// ([`NativeElement::is_native`]); the caller then converts through +/// [`read_raw_data_selection`] and the `read_as_*` functions. The selection is +/// validated like every selection read: out-of-range coordinates are +/// [`FormatError::SelectionOutOfBounds`]. +pub fn read_selection_native( + raw: &[u8], + dims: &[u64], + datatype: &Datatype, + selection: &crate::selection::Selection, +) -> Result>, FormatError> { + if !T::is_native(datatype) { + return Ok(None); + } + let elem_size = core::mem::size_of::(); + let total = dims + .iter() + .try_fold(1u64, |acc, &d| acc.checked_mul(d)) + .ok_or_else(|| FormatError::Overflow("dataset shape overflows".into()))?; + let expected = crate::chunked_read::checked_byte_len(total, elem_size)?; + if raw.len() != expected { + return Err(FormatError::DataSizeMismatch { + expected, + actual: raw.len(), + }); + } + if let crate::selection::Selection::All = selection { + return Ok(Some(native_to_vec(raw, expected / elem_size))); + } + crate::partial_read::validate(selection, dims)?; + crate::gather::gather::(raw, dims, elem_size, selection).map(Some) +} + /// Convert raw bytes to `f64` values. pub fn read_as_f64(raw: &[u8], datatype: &Datatype) -> Result, FormatError> { // Array datatypes read as a flat sequence of their base elements, and @@ -799,9 +840,8 @@ pub fn read_as_f64(raw: &[u8], datatype: &Datatype) -> Result, FormatEr let count = raw.len() / elem_size; // Fast path: native-endian f64 — single bulk memcpy - #[cfg(target_endian = "little")] - if is_native_le_float(datatype, FloatFormat::Double) { - return Ok(native_le_to_vec::(raw, count)); + if f64::is_native(datatype) { + return Ok(native_to_vec::(raw, count)); } let order = get_byte_order(datatype); @@ -941,19 +981,8 @@ pub fn read_as_i64(raw: &[u8], datatype: &Datatype) -> Result, FormatEr let count = raw.len() / elem_size; // Fast path: native LE i64 — single bulk memcpy - #[cfg(target_endian = "little")] - if elem_size == 8 - && is_full_width(datatype) - && matches!( - datatype, - Datatype::FixedPoint { - byte_order: DatatypeByteOrder::LittleEndian, - signed: true, - .. - } - ) - { - return Ok(native_le_to_vec::(raw, count)); + if i64::is_native(datatype) { + return Ok(native_to_vec::(raw, count)); } let order = get_byte_order(datatype); @@ -986,6 +1015,12 @@ pub fn read_as_u64(raw: &[u8], datatype: &Datatype) -> Result, FormatEr }); } let count = raw.len() / elem_size; + + // Fast path: native u64 — single bulk memcpy + if u64::is_native(datatype) { + return Ok(native_to_vec::(raw, count)); + } + let order = get_byte_order(datatype); let mut result = crate::bulk_alloc::vec_for_bulk(count); for i in 0..count { @@ -1013,9 +1048,8 @@ pub fn read_as_f32(raw: &[u8], datatype: &Datatype) -> Result, FormatEr let count = raw.len() / elem_size; // Fast path: native-endian f32 — single bulk memcpy - #[cfg(target_endian = "little")] - if is_native_le_float(datatype, FloatFormat::Single) { - return Ok(native_le_to_vec::(raw, count)); + if f32::is_native(datatype) { + return Ok(native_to_vec::(raw, count)); } // Little-endian IEEE half precision (numpy float16): widen directly. if is_native_le_float(datatype, FloatFormat::Half) { @@ -1103,19 +1137,8 @@ pub fn read_as_i32(raw: &[u8], datatype: &Datatype) -> Result, FormatEr let count = raw.len() / elem_size; // Fast path: native LE i32 — single bulk memcpy - #[cfg(target_endian = "little")] - if elem_size == 4 - && is_full_width(datatype) - && matches!( - datatype, - Datatype::FixedPoint { - byte_order: DatatypeByteOrder::LittleEndian, - signed: true, - .. - } - ) - { - return Ok(native_le_to_vec::(raw, count)); + if i32::is_native(datatype) { + return Ok(native_to_vec::(raw, count)); } let order = get_byte_order(datatype); diff --git a/crates/clawhdf5-format/src/gather.rs b/crates/clawhdf5-format/src/gather.rs new file mode 100644 index 0000000..22d2ae7 --- /dev/null +++ b/crates/clawhdf5-format/src/gather.rs @@ -0,0 +1,343 @@ +//! Copying a selection out of a row-major buffer one contiguous run at a time. +//! +//! A selection's elements, in output order, fall into runs that are adjacent +//! in the source: a whole block along the last dimension, blocks that touch +//! (`stride == block`), and whole rows when the inner dimensions are selected +//! in full. Copying run by run turns a 256 x 256 hyperslab of a 1024-wide +//! dataset into 256 `memcpy`s of 1 KiB, where the old extractor recursed and +//! bounds-checked once per element. + +#[cfg(not(feature = "std"))] +use alloc::{vec, vec::Vec}; + +use crate::data_read::NativeElement; +use crate::error::FormatError; +use crate::selection::Selection; + +/// Row-major element strides of `dims` (the last dimension has stride 1). +fn strides(dims: &[u64]) -> Vec { + let mut s = vec![1u64; dims.len()]; + for d in (0..dims.len().saturating_sub(1)).rev() { + s[d] = s[d + 1].wrapping_mul(dims[d + 1]); + } + s +} + +/// Merges adjacent runs before handing them on. +struct Coalesce { + start: u64, + len: u64, + emit: F, +} + +impl Coalesce { + #[inline] + fn push(&mut self, start: u64, len: u64) { + if len == 0 { + return; + } + if self.len > 0 && self.start.wrapping_add(self.len) == start { + self.len += len; + return; + } + self.flush(); + self.start = start; + self.len = len; + } + + fn flush(&mut self) { + if self.len > 0 { + (self.emit)(self.start, self.len); + self.len = 0; + } + } +} + +/// Call `emit(first_element, element_count)` for each run of a hyperslab's +/// elements that is contiguous in a row-major dataset of shape `dims`, in +/// the order the selection returns them. Adjacent runs are merged. +/// +/// Coordinates at or past a dimension's extent are skipped, as the +/// element-wise extractor always did; callers that want them to be an error +/// validate the selection first. The four vectors must have `dims.len()` +/// entries. +pub(crate) fn hyperslab_runs( + dims: &[u64], + start: &[u64], + stride: &[u64], + count: &[u64], + block: &[u64], + emit: impl FnMut(u64, u64), +) { + let rank = dims.len(); + let mut out = Coalesce { + start: 0, + len: 0, + emit, + }; + if rank == 0 { + out.push(0, 1); + out.flush(); + return; + } + if (0..rank).any(|d| count[d] == 0 || block[d] == 0) { + return; + } + let strides = strides(dims); + let last = rank - 1; + // Odometer over the outer dimensions: (block index, offset in block). + let mut ci = vec![0u64; last]; + let mut bi = vec![0u64; last]; + 'outer: loop { + // Base offset of this row, or skip it if a coordinate is out of range. + let mut base = 0u64; + let mut in_range = true; + for d in 0..last { + let coord = start[d] + .saturating_add(ci[d].saturating_mul(stride[d])) + .saturating_add(bi[d]); + if coord >= dims[d] { + in_range = false; + break; + } + base = base.wrapping_add(coord.wrapping_mul(strides[d])); + } + if in_range && (stride[last] == block[last] || count[last] == 1) { + // Blocks that touch (the common unit-stride case: block 1, + // stride 1) are one range; don't split it into per-element runs. + let s = start[last]; + let e = s + .saturating_add(count[last].saturating_mul(block[last])) + .min(dims[last]); + if s < e { + out.push(base.wrapping_add(s), e - s); + } + } else if in_range { + for c in 0..count[last] { + let s = start[last].saturating_add(c.saturating_mul(stride[last])); + if s >= dims[last] { + continue; + } + let e = s.saturating_add(block[last]).min(dims[last]); + out.push(base.wrapping_add(s), e - s); + } + } + // Advance the odometer, last outer dimension fastest. + let mut d = last; + loop { + if d == 0 { + break 'outer; + } + d -= 1; + bi[d] += 1; + if bi[d] < block[d] { + break; + } + bi[d] = 0; + ci[d] += 1; + if ci[d] < count[d] { + break; + } + ci[d] = 0; + } + } + out.flush(); +} + +/// The selected elements of `src` — a row-major dataset of shape `dims` and +/// `elem_size`-byte elements — copied into a fresh `Vec`, one `memcpy` per +/// contiguous run, with no zero-filling of the output first. +/// +/// For `T` other than `u8`, `elem_size` must equal `size_of::()`. The +/// selection must be a validated hyperslab, point list or `None` (`All` is the +/// caller's to handle); `src` must hold exactly the dataset. Anything that +/// would read outside `src` is an error, never a partial result. +pub(crate) fn gather( + src: &[u8], + dims: &[u64], + elem_size: usize, + selection: &Selection, +) -> Result, FormatError> { + let t_size = core::mem::size_of::(); + if elem_size == 0 || (t_size != 1 && t_size != elem_size) { + return Err(FormatError::DataSizeMismatch { + expected: t_size, + actual: elem_size, + }); + } + let n_elements = match selection { + Selection::None => 0, + Selection::Hyperslab { count, block, .. } => count + .iter() + .zip(block) + .try_fold(1u64, |acc, (&c, &b)| acc.checked_mul(c.checked_mul(b)?)) + .ok_or_else(|| FormatError::Overflow("hyperslab count x block overflows".into()))?, + Selection::Points(points) => points.len() as u64, + Selection::All => { + return Err(FormatError::SelectionOutOfBounds( + "gather does not take Selection::All".into(), + )); + } + }; + let out_bytes = crate::chunked_read::checked_byte_len(n_elements, elem_size)?; + let out_len = out_bytes / t_size; + let mut out: Vec = crate::bulk_alloc::vec_for_bulk(out_len); + let dst = out.as_mut_ptr().cast::(); + let mut written = 0usize; + let mut failed = false; + let mut copy_run = |first: u64, n: u64| { + if failed { + return; + } + let range = usize::try_from(first) + .ok() + .and_then(|f| f.checked_mul(elem_size)) + .zip( + usize::try_from(n) + .ok() + .and_then(|n| n.checked_mul(elem_size)), + ) + .and_then(|(at, len)| Some((at, len, at.checked_add(len)?))); + match range { + Some((at, len, end)) if end <= src.len() && written + len <= out_bytes => { + // SAFETY: `src[at..end]` is in bounds (checked above), and + // `dst + written .. + len` lies within `out`'s capacity of + // `out_bytes` bytes (checked above); `out` is a fresh + // allocation, so the regions do not overlap. + unsafe { + core::ptr::copy_nonoverlapping(src.as_ptr().add(at), dst.add(written), len) + }; + written += len; + } + _ => failed = true, + } + }; + let mut bad_point = false; + match selection { + Selection::Hyperslab { + start, + stride, + count, + block, + } => { + let rank = dims.len(); + if [start.len(), stride.len(), count.len(), block.len()] != [rank; 4] { + return Err(FormatError::SelectionOutOfBounds( + "hyperslab rank does not match dataset rank".into(), + )); + } + hyperslab_runs(dims, start, stride, count, block, &mut copy_run); + } + Selection::Points(points) => { + let strides = strides(dims); + let mut runs = Coalesce { + start: 0, + len: 0, + emit: &mut copy_run, + }; + for p in points { + if p.len() != dims.len() || p.iter().zip(dims).any(|(c, n)| c >= n) { + bad_point = true; + break; + } + let at = p + .iter() + .zip(&strides) + .fold(0u64, |acc, (c, s)| acc.wrapping_add(c.wrapping_mul(*s))); + runs.push(at, 1); + } + runs.flush(); + } + Selection::None | Selection::All => {} + } + if failed || bad_point || written != out_bytes { + return Err(FormatError::SelectionOutOfBounds( + "selection addresses elements outside the dataset".into(), + )); + } + // SAFETY: all `out_bytes` bytes, i.e. `out_len` values of `T`, were + // written above, and every bit pattern is a valid `T` (`NativeElement`). + unsafe { out.set_len(out_len) }; + Ok(out) +} + +#[cfg(test)] +mod tests { + use super::*; + + fn runs(dims: &[u64], sel: [&[u64]; 4]) -> Vec<(u64, u64)> { + let mut v = Vec::new(); + hyperslab_runs(dims, sel[0], sel[1], sel[2], sel[3], |s, n| v.push((s, n))); + v + } + + #[test] + fn runs_merge_blocks_and_whole_rows() { + // A box: one run per row. + assert_eq!( + runs(&[4, 10], [&[1, 2], &[1, 1], &[2, 3], &[1, 1]]), + vec![(12, 3), (22, 3)] + ); + // Whole rows: one run. + assert_eq!( + runs(&[4, 10], [&[1, 0], &[1, 1], &[3, 10], &[1, 1]]), + vec![(10, 30)] + ); + // stride == block: blocks merge. + assert_eq!( + runs(&[1, 10], [&[0, 1], &[1, 2], &[1, 4], &[1, 2]]), + vec![(1, 8)] + ); + // Strided with blocks along both dimensions. + assert_eq!( + runs(&[6, 10], [&[0, 1], &[3, 4], &[2, 2], &[2, 2]]), + vec![ + (1, 2), + (5, 2), + (11, 2), + (15, 2), + (31, 2), + (35, 2), + (41, 2), + (45, 2) + ] + ); + // Empty. + assert!(runs(&[4, 10], [&[0, 0], &[1, 1], &[0, 3], &[1, 1]]).is_empty()); + // Scalar. + assert_eq!(runs(&[], [&[], &[], &[], &[]]), vec![(0, 1)]); + } + + #[test] + fn gather_matches_element_order_and_rejects_out_of_range() { + let dims = [3u64, 4]; + let src: Vec = (0..12u16).flat_map(|v| v.to_le_bytes()).collect(); + let sel = Selection::Hyperslab { + start: vec![0, 1], + stride: vec![2, 2], + count: vec![2, 2], + block: vec![1, 1], + }; + let got: Vec = gather(&src, &dims, 2, &sel).unwrap(); + let want: Vec = [1u16, 3, 9, 11] + .iter() + .flat_map(|v| v.to_le_bytes()) + .collect(); + assert_eq!(got, want); + let pts = Selection::Points(vec![vec![2, 3], vec![0, 0], vec![0, 1]]); + let got: Vec = gather(&src, &dims, 2, &pts).unwrap(); + let want: Vec = [11u16, 0, 1].iter().flat_map(|v| v.to_le_bytes()).collect(); + assert_eq!(got, want); + // Past the extent, or a source shorter than the dataset: an error. + let bad = Selection::Points(vec![vec![3, 0]]); + assert!(gather::(&src, &dims, 2, &bad).is_err()); + let past = Selection::Hyperslab { + start: vec![2, 0], + stride: vec![1, 1], + count: vec![2, 4], + block: vec![1, 1], + }; + assert!(gather::(&src, &dims, 2, &past).is_err()); + assert!(gather::(&src[..20], &dims, 2, &pts).is_err()); + } +} diff --git a/crates/clawhdf5-format/src/lib.rs b/crates/clawhdf5-format/src/lib.rs index 05ebd46..23840b8 100644 --- a/crates/clawhdf5-format/src/lib.rs +++ b/crates/clawhdf5-format/src/lib.rs @@ -94,6 +94,7 @@ mod filters_szip; pub mod fixed_array; pub mod float16; pub mod fractal_heap; +mod gather; pub mod global_heap; pub mod group_info; pub mod group_v1; diff --git a/crates/clawhdf5-format/src/partial_read.rs b/crates/clawhdf5-format/src/partial_read.rs index 7d73599..9865c46 100644 --- a/crates/clawhdf5-format/src/partial_read.rs +++ b/crates/clawhdf5-format/src/partial_read.rs @@ -3,11 +3,13 @@ //! //! [`crate::data_read::read_raw_data_selection`] used to decode the *entire* //! dataset and then pick elements out of it, so reading a 64x64 window of a -//! large dataset took about as long as reading all of it. Here the selection's -//! bounding box is materialised instead — only the rows of a contiguous -//! dataset, or only the chunks, that overlap it — and the existing extractor -//! runs over that small buffer with the selection translated to the box's -//! origin. Extraction semantics are therefore exactly the full-read ones. +//! large dataset took about as long as reading all of it. A contiguous +//! dataset's selection is now copied straight out of the file, one `memcpy` +//! per contiguous run of selected elements (`crate::gather`). For chunked +//! data the selection's bounding box is materialised — only the chunks that +//! overlap it — and the extractor runs over that small buffer with the +//! selection translated to the box's origin. Extraction semantics are +//! therefore exactly the full-read ones. #[cfg(not(feature = "std"))] use alloc::string as alloc_or_std; @@ -250,10 +252,33 @@ pub fn read_selection( if dims.is_empty() || elem_size == 0 { return Ok(None); } + let total = dataspace.checked_num_elements()?; + // Contiguous data is addressable in place: copy the selection's runs + // straight out of it, whatever fraction of the dataset it covers, with no + // intermediate box (and no full copy for a large selection). + if let ( + DataLayout::Contiguous { + address: Some(address), + .. + }, + Selection::Hyperslab { .. } | Selection::Points(_), + ) = (layout, selection) + { + validate(selection, dims)?; + let base = usize::try_from(*address) + .map_err(|_| FormatError::Overflow("data address exceeds usize".into()))?; + let data = file_data + .get(base..) + .and_then(|d| d.get(..checked_byte_len(total, elem_size).ok()?)) + .ok_or(FormatError::UnexpectedEof { + expected: base, + available: file_data.len(), + })?; + return crate::gather::gather::(data, dims, elem_size, selection).map(Some); + } let Some((box_start, box_extent)) = bounding_box(selection, dims) else { return Ok(None); }; - let total = dataspace.checked_num_elements()?; let box_elements = box_extent .iter() .try_fold(1u64, |acc, &e| acc.checked_mul(e)) @@ -265,30 +290,6 @@ pub fn read_selection( let mut boxed = alloc_output(checked_byte_len(box_elements, elem_size)?)?; match layout { - DataLayout::Contiguous { - address: Some(address), - .. - } => { - let base = usize::try_from(*address) - .map_err(|_| FormatError::Overflow("data address exceeds usize".into()))?; - let data = file_data - .get(base..) - .and_then(|d| d.get(..checked_byte_len(total, elem_size).ok()?)) - .ok_or(FormatError::UnexpectedEof { - expected: base, - available: file_data.len(), - })?; - let origin = vec![0u64; dims.len()]; - copy_overlap( - data, - &origin, - dims, - &mut boxed, - &box_start, - &box_extent, - elem_size, - ); - } DataLayout::Chunked { btree_address: Some(_), .. diff --git a/crates/clawhdf5/src/reader.rs b/crates/clawhdf5/src/reader.rs index 83cbfa8..74c268d 100644 --- a/crates/clawhdf5/src/reader.rs +++ b/crates/clawhdf5/src/reader.rs @@ -569,9 +569,7 @@ impl<'f> Dataset<'f> { &self, selection: &clawhdf5_format::selection::Selection, ) -> Result, Error> { - let raw = self.read_selection(selection)?; - let dt = self.datatype()?; - Ok(data_read::read_as_f64(&raw, &dt)?) + self.read_typed_selection(selection, data_read::read_as_f64, || self.read_f64()) } /// Read selected elements as `f32` values. @@ -579,9 +577,7 @@ impl<'f> Dataset<'f> { &self, selection: &clawhdf5_format::selection::Selection, ) -> Result, Error> { - let raw = self.read_selection(selection)?; - let dt = self.datatype()?; - Ok(data_read::read_as_f32(&raw, &dt)?) + self.read_typed_selection(selection, data_read::read_as_f32, || self.read_f32()) } /// Read selected elements as `i32` values. @@ -589,9 +585,7 @@ impl<'f> Dataset<'f> { &self, selection: &clawhdf5_format::selection::Selection, ) -> Result, Error> { - let raw = self.read_selection(selection)?; - let dt = self.datatype()?; - Ok(data_read::read_as_i32(&raw, &dt)?) + self.read_typed_selection(selection, data_read::read_as_i32, || self.read_i32()) } /// Read selected elements as `i64` values. @@ -599,9 +593,34 @@ impl<'f> Dataset<'f> { &self, selection: &clawhdf5_format::selection::Selection, ) -> Result, Error> { - let raw = self.read_selection(selection)?; + self.read_typed_selection(selection, data_read::read_as_i64, || self.read_i64()) + } + + /// The typed selection readers. `All` is a full read. A contiguous dataset + /// that stores `T` natively is copied from the file straight into the + /// `Vec`, one copy per contiguous run of selected elements; anything + /// else reads the selection's bytes and converts them with `convert`. + fn read_typed_selection( + &self, + selection: &clawhdf5_format::selection::Selection, + convert: fn(&[u8], &Datatype) -> Result, FormatError>, + full: impl FnOnce() -> Result, Error>, + ) -> Result, Error> { + if matches!(selection, clawhdf5_format::selection::Selection::All) { + return full(); + } let dt = self.datatype()?; - Ok(data_read::read_as_i64(&raw, &dt)?) + if T::is_native(&dt) + && let Ok(Some(raw)) = self.read_raw_ref() + { + let dims = self.dataspace()?.dimensions; + if let Some(values) = data_read::read_selection_native::(raw, &dims, &dt, selection)? + { + return Ok(values); + } + } + let raw = self.read_selection(selection)?; + Ok(convert(&raw, &dt)?) } /// Zero-copy read of contiguous raw data. diff --git a/docs/known-issues.md b/docs/known-issues.md index bf8f77a..60fc9d8 100644 --- a/docs/known-issues.md +++ b/docs/known-issues.md @@ -9,7 +9,8 @@ deleting it. ## Concurrent and contiguous read performance (measured 2026-09-26) -**Status:** open. Measured on tank with `concurrent_read` against h5py +**Status:** open for chunked full reads; the contiguous item is fixed +(2026-09-26). Measured on tank with `concurrent_read` against h5py 3.16 / HDF5 2.0 (`BENCHMARKS.md`, "Concurrent reads"): - Full reads of chunked datasets from several threads through one `File` stop scaling at about 4 threads (880 MB/s on deflate data vs 4424 MB/s @@ -17,6 +18,12 @@ deleting it. scale to 1244 MB/s, so the `File`'s shared chunk cache is the suspect. - Contiguous datasets read 4x slower than h5py on one thread (2.5 vs 9.8 GB/s full, 0.12x for 256 x 256 hyperslabs). + **Fixed 2026-09-26** (not yet re-measured for `BENCHMARKS.md`): full + reads were dominated by 4 KiB page faults on the fresh output buffer, + which is now backed by transparent huge pages as numpy's is; hyperslab + reads copied the selection three times, element by element, and now copy + each contiguous run once, straight from the file into the output (see + `CHANGELOG.md`). The chunked-read scaling item above is still open. Values are correct; this is speed only. ## Silent wrong data found by the 2026-09-25 HDF5 audit From 6e8421a81e71bacddbe0a36e75857508f56f654f Mon Sep 17 00:00:00 2001 From: osobh Date: Sat, 26 Sep 2026 08:27:09 -0500 Subject: [PATCH 3/3] build: record libc in the conformance probe's lockfile clawhdf5-format now depends on libc on Linux (huge-page advice for read buffers); the probe's committed lockfile picks that up. Co-Authored-By: Claude Opus 5.5 (1M context) --- conformance/probe/Cargo.lock | 1 + 1 file changed, 1 insertion(+) diff --git a/conformance/probe/Cargo.lock b/conformance/probe/Cargo.lock index 964e0ab..36cfe85 100644 --- a/conformance/probe/Cargo.lock +++ b/conformance/probe/Cargo.lock @@ -64,6 +64,7 @@ dependencies = [ "bzip2", "flate2", "libaec-sys", + "libc", "lz4_flex", "pco", "portable-atomic",