CHANGELOG.md

@@ -16,6 +16,44 @@
   its lookups cost a few percent at 16 threads. Throughput with the default
   pool is unchanged, and still short of an h5py process pool.
 
+### 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.
+- **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
   written by h5py with `compression="lzf"`, or with hdf5plugin's

conformance/probe/Cargo.lock

@@ -64,6 +64,7 @@ dependencies = [
  "bzip2",
  "flate2",
  "libaec-sys",
+ "libc",
  "lz4_flex",
  "pco",
  "portable-atomic",

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,79 @@
+//! 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).
+#[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,
+/// 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

@@ -278,6 +278,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,
@@ -529,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())
@@ -538,96 +545,40 @@ pub fn extract_selection_from_buffer(
                 crate::chunked_read::checked_byte_len(output_elements, elem_size)?,
             )?;
 
-            // Compute dataset strides (row-major)
+            // One copy per run of elements contiguous in `full_data`
-            let mut ds_strides = vec![1usize; rank];
+            // (`gather`'s runs). Coordinates past the extent are skipped and
-            for i in (0..rank.saturating_sub(1)).rev() {
+            // runs past the end of `full_data` left as zeros, element by
-                ds_strides[i] = ds_strides[i + 1] * dims[i + 1] as usize;
+            // element, as this extractor always did; validated selections
-            }
+            // never hit either.
-
+            let mut out_at = 0usize;
-            // Compute output shape and strides
+            crate::gather::hyperslab_runs(dims, start, stride, count, block, |first, n| {
-            let output_dims: Vec<usize> = count
+                let big = |v: u64| usize::try_from(v).unwrap_or(usize::MAX);
-                .iter()
+                let (first, n) = (big(first), big(n));
-                .zip(block.iter())
+                let len = n.saturating_mul(elem_size);
-                .map(|(&c, &b)| (c * b) as usize)
+                let src = first.saturating_mul(elem_size);
-                .collect();
+                let out_end = out_at.saturating_add(len);
-            let mut out_strides = vec![1usize; rank];
+                if let (Some(from), Some(to)) = (
-            for i in (0..rank.saturating_sub(1)).rev() {
+                    full_data.get(src..src.saturating_add(len)),
-                out_strides[i] = out_strides[i + 1] * output_dims[i + 1];
+                    output.get_mut(out_at..out_end),
-            }
-
-            // 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 {
+                    to.copy_from_slice(from);
-                    // Copy one element
+                } else {
-                    let src = current_ds_offset * elem_size;
+                    for k in 0..n {
-                    let dst = *out_linear * elem_size;
+                        let s = first.saturating_add(k).saturating_mul(elem_size);
-                    if src + elem_size <= full_data.len() && dst + elem_size <= output.len() {
+                        let o = out_at.saturating_add(k.saturating_mul(elem_size));
-                        output[dst..dst + elem_size]
+                        if o >= output.len() {
-                            .copy_from_slice(&full_data[src..src + elem_size]);
+                            break;
                         }
-                    *out_linear += 1;
+                        if let (Some(from), Some(to)) = (
-                    return;
+                            full_data.get(s..s.saturating_add(elem_size)),
-                }
+                            output.get_mut(o..o.saturating_add(elem_size)),
-
+                        ) {
-                for bi in 0..count[d] {
+                            to.copy_from_slice(from);
-                    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],
-                            );
                         }
                     }
                 }
-            }
+                out_at = out_end;
-
+            });
-            iterate_hyperslab(
-                0,
-                rank,
-                start,
-                stride,
-                count,
-                block,
-                dims,
-                &ds_strides,
-                elem_size,
-                full_data,
-                &mut output,
-                &mut out_linear,
-                0,
-            );
 
             Ok(output)
         }
@@ -755,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
+mod sealed {
-/// little-endian target, in one copy.
+    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_to_vec<T: NativeElement>(raw: &[u8], count: usize) -> Vec<T> {
-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());
+    assert!(bytes <= raw.len(), "native_to_vec: source too short");
-    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
+    // bytes; `raw` holds at least `bytes` bytes (asserted); the regions
-    // `raw.len() / size_of::<T>()`); the regions cannot overlap because
+    // cannot overlap because `result` was just allocated. `T: NativeElement`
-    // `result` was just allocated. Every `T` used here (f32/f64/i32/i64) is
+    // is valid for any bit pattern, so after the copy all `count` values are
-    // 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::<u8>(), bytes);
@@ -779,6 +784,44 @@ fn native_le_to_vec<T: Copy>(raw: &[u8], count: usize) -> Vec<T> {
     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<T>`, 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<T: NativeElement>(
+    raw: &[u8],
+    dims: &[u64],
+    datatype: &Datatype,
+    selection: &crate::selection::Selection,
+) -> Result<Option<Vec<T>>, FormatError> {
+    if !T::is_native(datatype) {
+        return Ok(None);
+    }
+    let elem_size = core::mem::size_of::<T>();
+    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::<T>(raw, dims, elem_size, selection).map(Some)
+}
+
 /// Convert raw bytes to `f64` values.
 pub fn read_as_f64(raw: &[u8], datatype: &Datatype) -> Result<Vec<f64>, FormatError> {
     // Array datatypes read as a flat sequence of their base elements, and
@@ -797,13 +840,12 @@ pub fn read_as_f64(raw: &[u8], datatype: &Datatype) -> Result<Vec<f64>, FormatEr
     let count = raw.len() / elem_size;
 
     // Fast path: native-endian f64 — single bulk memcpy
-    #[cfg(target_endian = "little")]
+    if f64::is_native(datatype) {
-    if is_native_le_float(datatype, FloatFormat::Double) {
+        return Ok(native_to_vec::<f64>(raw, count));
-        return Ok(native_le_to_vec::<f64>(raw, count));
     }
 
     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) {
@@ -939,23 +981,12 @@ pub fn read_as_i64(raw: &[u8], datatype: &Datatype) -> Result<Vec<i64>, FormatEr
     let count = raw.len() / elem_size;
 
     // Fast path: native LE i64 — single bulk memcpy
-    #[cfg(target_endian = "little")]
+    if i64::is_native(datatype) {
-    if elem_size == 8
+        return Ok(native_to_vec::<i64>(raw, count));
-        && is_full_width(datatype)
-        && matches!(
-            datatype,
-            Datatype::FixedPoint {
-                byte_order: DatatypeByteOrder::LittleEndian,
-                signed: true,
-                ..
-            }
-        )
-    {
-        return Ok(native_le_to_vec::<i64>(raw, count));
     }
 
     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());
@@ -984,8 +1015,14 @@ pub fn read_as_u64(raw: &[u8], datatype: &Datatype) -> Result<Vec<u64>, FormatEr
         });
     }
     let count = raw.len() / elem_size;
+
+    // Fast path: native u64 — single bulk memcpy
+    if u64::is_native(datatype) {
+        return Ok(native_to_vec::<u64>(raw, count));
+    }
+
     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());
@@ -1011,21 +1048,23 @@ pub fn read_as_f32(raw: &[u8], datatype: &Datatype) -> Result<Vec<f32>, FormatEr
     let count = raw.len() / elem_size;
 
     // Fast path: native-endian f32 — single bulk memcpy
-    #[cfg(target_endian = "little")]
+    if f32::is_native(datatype) {
-    if is_native_le_float(datatype, FloatFormat::Single) {
+        return Ok(native_to_vec::<f32>(raw, count));
-        return Ok(native_le_to_vec::<f32>(raw, count));
     }
     // 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);
+        result.extend(
+            halves
                 .iter()
-            .map(|&b| f16_bits_to_f32(u16::from_le_bytes(b)))
+                .map(|&b| f16_bits_to_f32(u16::from_le_bytes(b))),
-            .collect());
+        );
+        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) {
@@ -1098,23 +1137,12 @@ pub fn read_as_i32(raw: &[u8], datatype: &Datatype) -> Result<Vec<i32>, FormatEr
     let count = raw.len() / elem_size;
 
     // Fast path: native LE i32 — single bulk memcpy
-    #[cfg(target_endian = "little")]
+    if i32::is_native(datatype) {
-    if elem_size == 4
+        return Ok(native_to_vec::<i32>(raw, count));
-        && is_full_width(datatype)
-        && matches!(
-            datatype,
-            Datatype::FixedPoint {
-                byte_order: DatatypeByteOrder::LittleEndian,
-                signed: true,
-                ..
-            }
-        )
-    {
-        return Ok(native_le_to_vec::<i32>(raw, count));
     }
 
     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/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<u64> {
+    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<F: FnMut(u64, u64)> {
+    start: u64,
+    len: u64,
+    emit: F,
+}
+
+impl<F: FnMut(u64, u64)> Coalesce<F> {
+    #[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<T>`, one `memcpy` per
+/// contiguous run, with no zero-filling of the output first.
+///
+/// For `T` other than `u8`, `elem_size` must equal `size_of::<T>()`. 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<T: NativeElement>(
+    src: &[u8],
+    dims: &[u64],
+    elem_size: usize,
+    selection: &Selection,
+) -> Result<Vec<T>, FormatError> {
+    let t_size = core::mem::size_of::<T>();
+    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<T> = crate::bulk_alloc::vec_for_bulk(out_len);
+    let dst = out.as_mut_ptr().cast::<u8>();
+    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<u8> = (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<u8> = gather(&src, &dims, 2, &sel).unwrap();
+        let want: Vec<u8> = [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<u8> = gather(&src, &dims, 2, &pts).unwrap();
+        let want: Vec<u8> = [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::<u8>(&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::<u8>(&src, &dims, 2, &past).is_err());
+        assert!(gather::<u8>(&src[..20], &dims, 2, &pts).is_err());
+    }
+}

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;
@@ -93,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;

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
+//! large dataset took about as long as reading all of it. A contiguous
-//! bounding box is materialised instead — only the rows of a contiguous
+//! dataset's selection is now copied straight out of the file, one `memcpy`
-//! dataset, or only the chunks, that overlap it — and the existing extractor
+//! per contiguous run of selected elements (`crate::gather`). For chunked
-//! runs over that small buffer with the selection translated to the box's
+//! data the selection's bounding box is materialised — only the chunks that
-//! origin. Extraction semantics are therefore exactly the full-read ones.
+//! 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::<u8>(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(_),
             ..

crates/clawhdf5/src/reader.rs

@@ -569,9 +569,7 @@ impl<'f> Dataset<'f> {
         &self,
         selection: &clawhdf5_format::selection::Selection,
     ) -> Result<Vec<f64>, Error> {
-        let raw = self.read_selection(selection)?;
+        self.read_typed_selection(selection, data_read::read_as_f64, || self.read_f64())
-        let dt = self.datatype()?;
-        Ok(data_read::read_as_f64(&raw, &dt)?)
     }
 
     /// Read selected elements as `f32` values.
@@ -579,9 +577,7 @@ impl<'f> Dataset<'f> {
         &self,
         selection: &clawhdf5_format::selection::Selection,
     ) -> Result<Vec<f32>, Error> {
-        let raw = self.read_selection(selection)?;
+        self.read_typed_selection(selection, data_read::read_as_f32, || self.read_f32())
-        let dt = self.datatype()?;
-        Ok(data_read::read_as_f32(&raw, &dt)?)
     }
 
     /// Read selected elements as `i32` values.
@@ -589,9 +585,7 @@ impl<'f> Dataset<'f> {
         &self,
         selection: &clawhdf5_format::selection::Selection,
     ) -> Result<Vec<i32>, Error> {
-        let raw = self.read_selection(selection)?;
+        self.read_typed_selection(selection, data_read::read_as_i32, || self.read_i32())
-        let dt = self.datatype()?;
-        Ok(data_read::read_as_i32(&raw, &dt)?)
     }
 
     /// Read selected elements as `i64` values.
@@ -599,9 +593,34 @@ impl<'f> Dataset<'f> {
         &self,
         selection: &clawhdf5_format::selection::Selection,
     ) -> Result<Vec<i64>, 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<T>`, one copy per contiguous run of selected elements; anything
+    /// else reads the selection's bytes and converts them with `convert`.
+    fn read_typed_selection<T: data_read::NativeElement>(
+        &self,
+        selection: &clawhdf5_format::selection::Selection,
+        convert: fn(&[u8], &Datatype) -> Result<Vec<T>, FormatError>,
+        full: impl FnOnce() -> Result<Vec<T>, Error>,
+    ) -> Result<Vec<T>, 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::<T>(raw, &dims, &dt, selection)?
+            {
+                return Ok(values);
+            }
+        }
+        let raw = self.read_selection(selection)?;
+        Ok(convert(&raw, &dt)?)
     }
 
     /// Zero-copy read of contiguous raw data.

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"
+        );
+    }
+}

docs/known-issues.md

@@ -9,7 +9,8 @@ deleting it.
 
 ## Concurrent and contiguous read performance (measured 2026-09-26)
 
-**Status:** open; one cause of the first bullet fixed (2026-09-26). Measured on
+**Status:** open for chunked full reads (one cause fixed 2026-09-26); 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"):
 - **Partly fixed 2026-09-26.** Full reads of chunked datasets from several threads
@@ -34,6 +35,12 @@ tank with `concurrent_read` against h5py 3.16 / HDF5 2.0 (`BENCHMARKS.md`,
   the `f32` copy of it, and a new buffer per decoded chunk).
 - 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