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

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,

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<usize> = 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,
+            // 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),
                 ) {
-                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]);
+                    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;
                         }
-                    *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],
-                            );
+                        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<T: Copy>(raw: &[u8], count: usize) -> Vec<T> {
+fn native_to_vec<T: NativeElement>(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> = crate::bulk_alloc::vec_for_bulk(count);
     // SAFETY: `result` has capacity for `count` values of `T`, i.e. `bytes`
-    // bytes; `raw` holds at least `bytes` bytes (callers derive `count` from
-    // `raw.len() / size_of::<T>()`); the regions cannot overlap because
-    // `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::<u8>(), bytes);
@@ -781,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
@@ -799,9 +840,8 @@ 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 is_native_le_float(datatype, FloatFormat::Double) {
-        return Ok(native_le_to_vec::<f64>(raw, count));
+    if f64::is_native(datatype) {
+        return Ok(native_to_vec::<f64>(raw, count));
     }
 
     let order = get_byte_order(datatype);
@@ -941,19 +981,8 @@ 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 elem_size == 8
-        && is_full_width(datatype)
-        && matches!(
-            datatype,
-            Datatype::FixedPoint {
-                byte_order: DatatypeByteOrder::LittleEndian,
-                signed: true,
-                ..
-            }
-        )
-    {
-        return Ok(native_le_to_vec::<i64>(raw, count));
+    if i64::is_native(datatype) {
+        return Ok(native_to_vec::<i64>(raw, count));
     }
 
     let order = get_byte_order(datatype);
@@ -986,6 +1015,12 @@ 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 = 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<Vec<f32>, 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::<f32>(raw, count));
+    if f32::is_native(datatype) {
+        return Ok(native_to_vec::<f32>(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<Vec<i32>, 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::<i32>(raw, count));
+    if i32::is_native(datatype) {
+        return Ok(native_to_vec::<i32>(raw, count));
     }
 
     let order = get_byte_order(datatype);

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

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

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::<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)?;
-        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<Vec<f32>, 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<Vec<i32>, 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<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.

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