crates/clawhdf5-format/src/chunk_grid.rs

@@ -0,0 +1,200 @@
+//! Chunk-index linearisation shared by the Fixed Array and Extensible Array
+//! chunk indexes (reader and writer).
+//!
+//! Both indexes store one element per chunk at a *linear* index, and the
+//! library derives that index from the chunk's scaled coordinates
+//! (`offset / chunk_dim`) using the dataset's **maximum** dimensions, not its
+//! current ones (`H5D__farray_idx_get_addr` / `H5D__earray_idx_get_addr`,
+//! via `layout->max_down_chunks`). A dataset whose current shape is smaller
+//! than its maxshape therefore has gaps in the index, and laying it out by the
+//! current shape puts every chunk after the first row in the wrong place.
+//!
+//! The Extensible Array adds one more step: its one unlimited dimension has no
+//! finite chunk count, so the library *swizzles* the coordinates to make that
+//! dimension the slowest-varying one (`H5VM_swizzle_coords`, which moves
+//! `coords[unlim_dim]` to the front and shifts the dimensions before it right
+//! by one) before linearising with `swizzled_max_down_chunks`. When the
+//! unlimited dimension is already dimension 0 no swizzle happens.
+
+#[cfg(not(feature = "std"))]
+extern crate alloc;
+
+#[cfg(not(feature = "std"))]
+use alloc::{vec, vec::Vec};
+
+use crate::error::FormatError;
+
+/// How a chunk index maps linear element indexes to chunk coordinates.
+#[derive(Debug, Clone)]
+pub(crate) struct ChunkGrid {
+    /// Spatial chunk dimensions, in dataset order.
+    chunk_dims: Vec<u64>,
+    /// Chunks per dimension covering the *current* extent, in dataset order.
+    cur_chunks: Vec<u64>,
+    /// Dataset dimension stored at each linearisation position (slowest
+    /// first). The identity except for a swizzled Extensible Array.
+    order: Vec<usize>,
+    /// Linear stride of each linearisation position.
+    down: Vec<u64>,
+}
+
+impl ChunkGrid {
+    /// Grid for a Fixed Array index: row-major over the chunk counts of the
+    /// maximum dimensions (`max_dims`, falling back to the current dimensions
+    /// when the dataspace records none).
+    pub(crate) fn fixed_array(
+        cur_dims: &[u64],
+        max_dims: Option<&[u64]>,
+        chunk_dims: &[u64],
+    ) -> Result<Self, FormatError> {
+        Self::build(cur_dims, max_dims, chunk_dims, None)
+    }
+
+    /// Grid for an Extensible Array index: like the Fixed Array, but the
+    /// unlimited dimension (the one whose maximum is `H5S_UNLIMITED`) is moved
+    /// to the slowest-varying position first.
+    pub(crate) fn extensible_array(
+        cur_dims: &[u64],
+        max_dims: Option<&[u64]>,
+        chunk_dims: &[u64],
+    ) -> Result<Self, FormatError> {
+        let unlim = max_dims.and_then(|m| m.iter().position(|&d| d == u64::MAX));
+        Self::build(cur_dims, max_dims, chunk_dims, unlim)
+    }
+
+    fn build(
+        cur_dims: &[u64],
+        max_dims: Option<&[u64]>,
+        chunk_dims: &[u64],
+        unlim: Option<usize>,
+    ) -> Result<Self, FormatError> {
+        let rank = chunk_dims.len();
+        if cur_dims.len() != rank || max_dims.is_some_and(|m| m.len() != rank) {
+            return Err(FormatError::ChunkedReadError(
+                "chunk index rank does not match the dataspace".into(),
+            ));
+        }
+        if chunk_dims.contains(&0) {
+            return Err(FormatError::ChunkedReadError(
+                "chunk dimension is zero".into(),
+            ));
+        }
+        let cur_chunks: Vec<u64> = cur_dims
+            .iter()
+            .zip(chunk_dims)
+            .map(|(&d, &c)| d.div_ceil(c))
+            .collect();
+        // Chunk counts of the maximum extent. An unlimited dimension has no
+        // finite count; it only ever sits in the slowest position, where its
+        // count never enters a stride. A (corrupt) maximum smaller than the
+        // current extent is widened so no allocated chunk becomes unreachable.
+        let max_chunks: Vec<u64> = (0..rank)
+            .map(|d| {
+                let max = max_dims.map_or(cur_dims[d], |m| m[d]);
+                if max == u64::MAX {
+                    u64::MAX
+                } else {
+                    max.div_ceil(chunk_dims[d]).max(cur_chunks[d])
+                }
+            })
+            .collect();
+
+        let mut order: Vec<usize> = (0..rank).collect();
+        if let Some(u) = unlim {
+            order.remove(u);
+            order.insert(0, u);
+        }
+        let mut down = vec![1u64; rank];
+        for p in (0..rank.saturating_sub(1)).rev() {
+            let next = max_chunks[order[p + 1]];
+            if next == u64::MAX {
+                // Only reachable with more than one unlimited dimension, which
+                // neither index type can describe.
+                return Err(FormatError::ChunkedReadError(
+                    "array chunk index with more than one unlimited dimension".into(),
+                ));
+            }
+            down[p] = down[p + 1].checked_mul(next).ok_or_else(|| {
+                FormatError::Overflow("chunk index linear stride overflows u64".into())
+            })?;
+        }
+        Ok(Self {
+            chunk_dims: chunk_dims.to_vec(),
+            cur_chunks,
+            order,
+            down,
+        })
+    }
+
+    /// Dataset-space offsets of the chunk stored at linear `index`, or `None`
+    /// when that chunk lies outside the current extent (the index still has a
+    /// slot for it; the library ignores such chunks on read).
+    pub(crate) fn offsets(&self, index: u64) -> Option<Vec<u64>> {
+        let rank = self.chunk_dims.len();
+        let mut offsets = vec![0u64; rank];
+        let mut rem = index;
+        for p in 0..rank {
+            let d = self.order[p];
+            let scaled = rem / self.down[p];
+            rem %= self.down[p];
+            if scaled >= self.cur_chunks[d] {
+                return None;
+            }
+            offsets[d] = scaled * self.chunk_dims[d];
+        }
+        Some(offsets)
+    }
+
+    /// Linear index of the chunk with scaled coordinates `scaled`
+    /// (`offset / chunk_dim` per dimension, in dataset order).
+    pub(crate) fn linear_index(&self, scaled: &[u64]) -> u64 {
+        self.order
+            .iter()
+            .zip(&self.down)
+            .map(|(&d, &stride)| scaled[d] * stride)
+            .sum()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn fixed_array_uses_max_dims() {
+        // shape (4, 6), chunks (2, 3), maxshape (20, 10): 10 x 4 chunk grid.
+        let g = ChunkGrid::fixed_array(&[4, 6], Some(&[20, 10]), &[2, 3]).unwrap();
+        assert_eq!(g.offsets(0), Some(vec![0, 0]));
+        assert_eq!(g.offsets(1), Some(vec![0, 3]));
+        assert_eq!(g.offsets(2), None); // column chunk 2 is beyond the extent
+        assert_eq!(g.offsets(4), Some(vec![2, 0]));
+        assert_eq!(g.offsets(5), Some(vec![2, 3]));
+        assert_eq!(g.offsets(8), None); // row chunk 2 is beyond the extent
+        assert_eq!(g.linear_index(&[1, 1]), 5);
+    }
+
+    #[test]
+    fn extensible_array_swizzles_unlimited_dim() {
+        // maxshape (10, None): dim 1 is unlimited and becomes slowest.
+        let g = ChunkGrid::extensible_array(&[4, 6], Some(&[10, u64::MAX]), &[2, 3]).unwrap();
+        // max chunks of dim 0 = 5, so index = c1 * 5 + c0.
+        assert_eq!(g.linear_index(&[1, 0]), 1);
+        assert_eq!(g.linear_index(&[0, 1]), 5);
+        assert_eq!(g.offsets(5), Some(vec![0, 3]));
+        assert_eq!(g.offsets(6), Some(vec![2, 3]));
+        assert_eq!(g.offsets(2), None);
+    }
+
+    #[test]
+    fn extensible_array_unlimited_first_is_row_major() {
+        let g = ChunkGrid::extensible_array(&[4, 6], Some(&[u64::MAX, 30]), &[2, 3]).unwrap();
+        // max chunks of dim 1 = 10.
+        assert_eq!(g.linear_index(&[1, 1]), 11);
+        assert_eq!(g.offsets(11), Some(vec![2, 3]));
+    }
+
+    #[test]
+    fn rejects_two_unlimited_dims_after_the_first() {
+        assert!(ChunkGrid::fixed_array(&[4, 6], Some(&[u64::MAX, u64::MAX]), &[2, 3]).is_err());
+    }
+}

crates/clawhdf5-format/src/chunked_read.rs

@@ -593,6 +593,7 @@ pub fn list_chunks(
                 file_data,
                 &header,
                 &dataspace.dimensions,
+                dataspace.max_dimensions.as_deref(),
                 spatial_chunk_dims,
                 elem_size as u32,
                 offset_size,
@@ -608,6 +609,7 @@ pub fn list_chunks(
                 file_data,
                 &header,
                 &dataspace.dimensions,
+                dataspace.max_dimensions.as_deref(),
                 spatial_chunk_dims,
                 elem_size as u32,
                 offset_size,

crates/clawhdf5-format/src/chunked_write.rs

@@ -4,10 +4,11 @@
 extern crate alloc;
 
 #[cfg(not(feature = "std"))]
-use alloc::{vec, vec::Vec};
+use alloc::{format, vec, vec::Vec};
 
 use crate::checksum::jenkins_lookup3;
 use crate::chunk_cache::{CACHE_LINE_SIZE, align_to_cache_line};
+use crate::chunk_grid::ChunkGrid;
 use crate::ea_writer;
 use crate::error::FormatError;
 use crate::filter_pipeline::{
@@ -443,6 +444,27 @@ fn serialize_v4_fixed_array(
     element_size: u32,
     max_bits: u8,
 ) -> Vec<u8> {
+    let mut buf = layout_v4_chunked_prefix(chunk_dims, element_size);
+
+    // chunk index type = 3 (Fixed Array)
+    buf.push(3);
+
+    // max_dblk_page_nelmts_bits — must match FAHD max_nelmts_bits
+    buf.push(max_bits);
+
+    // Fixed Array header address
+    match offset_size {
+        4 => buf.extend_from_slice(&(fixed_array_address as u32).to_le_bytes()),
+        8 => buf.extend_from_slice(&fixed_array_address.to_le_bytes()),
+        _ => {}
+    }
+
+    buf
+}
+
+/// The part of a v4 chunked layout message before the chunk index type:
+/// version, class, flags and the chunk dimensions (plus the element size).
+fn layout_v4_chunked_prefix(chunk_dims: &[u32], element_size: u32) -> Vec<u8> {
     let mut buf = Vec::new();
     buf.push(4); // version
     buf.push(2); // class = chunked
@@ -482,124 +504,142 @@ fn serialize_v4_fixed_array(
         4 => buf.extend_from_slice(&element_size.to_le_bytes()),
         _ => {}
     }
-
-    // chunk index type = 3 (Fixed Array)
-    buf.push(3);
-
-    // max_dblk_page_nelmts_bits — must match FAHD max_nelmts_bits
-    buf.push(max_bits);
-
-    // Fixed Array header address
-    match offset_size {
-        4 => buf.extend_from_slice(&(fixed_array_address as u32).to_le_bytes()),
-        8 => buf.extend_from_slice(&fixed_array_address.to_le_bytes()),
-        _ => {}
-    }
-
     buf
 }
 
+/// log2 of the elements per Fixed Array data block page (the library's
+/// default, `H5D_FARRAY_MAX_DBLK_PAGE_NELMTS_BITS`).
+const FA_PAGE_BITS: u8 = 10;
+
+pub(crate) fn push_addr(buf: &mut Vec<u8>, addr: u64, offset_size: u8) {
+    match offset_size {
+        4 => buf.extend_from_slice(&(addr as u32).to_le_bytes()),
+        _ => buf.extend_from_slice(&addr.to_le_bytes()),
+    }
+}
+
+/// Width of the chunk-size field of a filtered chunk index element. Must
+/// match the library's `H5D_FARRAY_FILT_COMPUTE_CHUNK_SIZE_LEN` (the EA and
+/// B-tree v2 indexes use the same formula):
+/// `1 + ((log2(unfiltered chunk bytes) + 8) / 8)`, capped at 8.
+pub(crate) fn filtered_chunk_size_len(slots: &[Option<WrittenChunk>]) -> usize {
+    let max_raw = slots
+        .iter()
+        .flatten()
+        .map(|c| c.raw_size)
+        .max()
+        .unwrap_or(1);
+    let log2_val = if max_raw <= 1 {
+        0
+    } else {
+        63 - max_raw.leading_zeros()
+    };
+    (1 + ((log2_val + 8) / 8) as usize).min(8)
+}
+
+/// Append one chunk index element: the chunk's address, plus its stored size
+/// and filter mask when the dataset is filtered. `None` is an unallocated
+/// chunk (undefined address, zero size and mask).
+pub(crate) fn push_index_element(
+    buf: &mut Vec<u8>,
+    slot: Option<&WrittenChunk>,
+    offset_size: u8,
+    chunk_size_bytes: Option<usize>,
+) {
+    match slot {
+        Some(c) => {
+            push_addr(buf, c.address, offset_size);
+            if let Some(n) = chunk_size_bytes {
+                buf.extend_from_slice(&c.compressed_size.to_le_bytes()[..n]);
+                buf.extend_from_slice(&c.filter_mask.to_le_bytes());
+            }
+        }
+        None => {
+            buf.extend(core::iter::repeat_n(0xFF, offset_size as usize));
+            if let Some(n) = chunk_size_bytes {
+                buf.extend(core::iter::repeat_n(0x00, n + 4));
+            }
+        }
+    }
+}
+
 /// Build a complete Fixed Array at a known absolute address.
+///
+/// `slots` holds one entry per element of the array, i.e. per chunk of the
+/// dataset's *maximum* extent in the order [`crate::chunk_grid`] defines;
+/// `None` marks a chunk that is not allocated. An array with more elements
+/// than fit in one page (`2^FA_PAGE_BITS`) gets a paged data block: a
+/// page-init bitmap after the prefix, then one checksummed page per
+/// `2^FA_PAGE_BITS` elements, the last one short (`H5FA__dblock_create`).
 pub fn build_fixed_array_at(
-    chunks: &[WrittenChunk],
+    slots: &[Option<WrittenChunk>],
     offset_size: u8,
     length_size: u8,
     has_filters: bool,
     fa_base_address: u64,
 ) -> Vec<u8> {
     let os = offset_size as usize;
-    let num_elements = chunks.len();
+    let num_elements = slots.len();
-
-    // For filtered chunks, compute chunk_size encoding width.
-    // Must match the HDF5 C library's H5D_FARRAY_FILT_COMPUTE_CHUNK_SIZE_LEN macro:
-    //   chunk_size_len = 1 + ((H5VM_log2_gen(chunk.size) + 8) / 8)
-    // where chunk.size is the unfiltered chunk size in bytes (product of all chunk dims).
-    let chunk_size_bytes: usize = if has_filters {
-        let max_raw = chunks.iter().map(|c| c.raw_size).max().unwrap_or(1);
-        let log2_val = if max_raw <= 1 {
-            0
-        } else {
-            63 - max_raw.leading_zeros()
-        };
-        let len = 1 + ((log2_val + 8) / 8) as usize;
-        len.min(8)
-    } else {
-        0
-    };
-
-    let elem_size = if has_filters {
-        os + chunk_size_bytes + 4
-    } else {
-        os
-    };
 
+    let chunk_size_bytes = has_filters.then(|| filtered_chunk_size_len(slots));
+    let elem_size = os + chunk_size_bytes.map_or(0, |n| n + 4);
     let client_id: u8 = if has_filters { 1 } else { 0 };
 
     // FAHD total size
-    let nelmts_field_size = length_size as usize;
+    let fahd_total_size = 4 + 1 + 1 + 1 + 1 + length_size as usize + os + 4;
-    let fahd_total_size = 4 + 1 + 1 + 1 + 1 + nelmts_field_size + os + 4;
     let fadb_address = fa_base_address + fahd_total_size as u64;
 
-    // Build FAHD
     let mut fahd = Vec::with_capacity(fahd_total_size);
     fahd.extend_from_slice(b"FAHD");
     fahd.push(0); // version
     fahd.push(client_id);
     fahd.push(elem_size as u8);
-
+    fahd.push(FA_PAGE_BITS);
-    // max_nelmts_bits: use 10 as default (page_size = 1024), matching h5py convention
-    let max_bits: u8 = 10;
-    fahd.push(max_bits);
-
     match length_size {
         4 => fahd.extend_from_slice(&(num_elements as u32).to_le_bytes()),
-        8 => fahd.extend_from_slice(&(num_elements as u64).to_le_bytes()),
         _ => fahd.extend_from_slice(&(num_elements as u64).to_le_bytes()),
     }
-
+    push_addr(&mut fahd, fadb_address, offset_size);
-    match offset_size {
-        4 => fahd.extend_from_slice(&(fadb_address as u32).to_le_bytes()),
-        8 => fahd.extend_from_slice(&fadb_address.to_le_bytes()),
-        _ => fahd.extend_from_slice(&fadb_address.to_le_bytes()),
-    }
-
-    // Checksum
     let checksum = jenkins_lookup3(&fahd);
     fahd.extend_from_slice(&checksum.to_le_bytes());
-
     assert_eq!(fahd.len(), fahd_total_size);
 
-    // Build FADB
+    // FADB prefix
     let mut fadb = Vec::new();
     fadb.extend_from_slice(b"FADB");
     fadb.push(0); // version
     fadb.push(client_id);
+    push_addr(&mut fadb, fa_base_address, offset_size);
 
-    // header address
+    let page_nelmts = 1usize << FA_PAGE_BITS;
-    match offset_size {
+    if num_elements <= page_nelmts {
-        4 => fadb.extend_from_slice(&(fa_base_address as u32).to_le_bytes()),
+        // Unpaged: the elements follow the prefix, one checksum over both.
-        8 => fadb.extend_from_slice(&fa_base_address.to_le_bytes()),
+        for slot in slots {
-        _ => fadb.extend_from_slice(&fa_base_address.to_le_bytes()),
+            push_index_element(&mut fadb, slot.as_ref(), offset_size, chunk_size_bytes);
         }
-
-    // Element data
-    for chunk in chunks {
-        match offset_size {
-            4 => fadb.extend_from_slice(&(chunk.address as u32).to_le_bytes()),
-            8 => fadb.extend_from_slice(&chunk.address.to_le_bytes()),
-            _ => fadb.extend_from_slice(&chunk.address.to_le_bytes()),
-        }
-        if has_filters {
-            // Write compressed size using chunk_size_bytes (variable width)
-            let cs_bytes = chunk.compressed_size.to_le_bytes();
-            fadb.extend_from_slice(&cs_bytes[..chunk_size_bytes]);
-            fadb.extend_from_slice(&chunk.filter_mask.to_le_bytes());
-        }
-    }
-
-    // FADB checksum
         let fadb_checksum = jenkins_lookup3(&fadb);
         fadb.extend_from_slice(&fadb_checksum.to_le_bytes());
+    } else {
+        // Paged: every page is written, so every page-init bit is set
+        // (MSB-first, as `H5VM_bit_set` packs them). The prefix and bitmap
+        // share a checksum; each page carries its own.
+        let npages = num_elements.div_ceil(page_nelmts);
+        let mut bitmap = vec![0u8; npages.div_ceil(8)];
+        for p in 0..npages {
+            bitmap[p / 8] |= 0x80 >> (p % 8);
+        }
+        fadb.extend_from_slice(&bitmap);
+        let prefix_checksum = jenkins_lookup3(&fadb);
+        fadb.extend_from_slice(&prefix_checksum.to_le_bytes());
+        for page in slots.chunks(page_nelmts) {
+            let start = fadb.len();
+            for slot in page {
+                push_index_element(&mut fadb, slot.as_ref(), offset_size, chunk_size_bytes);
+            }
+            let page_checksum = jenkins_lookup3(&fadb[start..]);
+            fadb.extend_from_slice(&page_checksum.to_le_bytes());
+        }
+    }
 
     let mut combined = fahd;
     combined.extend_from_slice(&fadb);
@@ -667,7 +707,8 @@ pub fn build_chunked_data_from_precompressed(
     pre: &PrecompressedChunks,
     base_address: u64,
     maxshape: Option<&[u64]>,
-) -> ChunkedDataResult {
+) -> Result<ChunkedDataResult, FormatError> {
+    let index = ChunkIndexPlan::new(&pre.shape, maxshape, &pre.chunk_dims)?;
     let offset_size: u8 = 8;
     let length_size: u8 = 8;
     let num_chunks = pre.chunks.len();
@@ -693,17 +734,18 @@ pub fn build_chunked_data_from_precompressed(
     }
 
     let chunk_dims_u32: Vec<u32> = pre.chunk_dims.iter().map(|&d| d as u32).collect();
-    let use_extensible = maxshape.is_some_and(|ms| ms.contains(&u64::MAX));
 
     let aligned_idx = align_to_cache_line(data_buf.len());
     if aligned_idx > data_buf.len() {
         data_buf.resize(aligned_idx, 0u8);
     }
 
-    let layout_message = if use_extensible {
+    let layout_message = match &index {
+        ChunkIndexPlan::ExtensibleArray(grid) => {
             let ea_address = base_address + data_buf.len() as u64;
+            let slots = index_slots(grid, &pre.shape, &pre.chunk_dims, &written_chunks, None)?;
             let ea_bytes = ea_writer::build_extensible_array_at(
-            &written_chunks,
+                &slots,
                 offset_size,
                 length_size,
                 pre.has_filters,
@@ -716,7 +758,8 @@ pub fn build_chunked_data_from_precompressed(
                 offset_size,
                 element_size as u32,
             )
-    } else if num_chunks == 1 {
+        }
+        ChunkIndexPlan::SingleChunk => {
             let chunk_addr = written_chunks[0].address;
             let filtered_size = if pre.has_filters {
                 Some(written_chunks[0].compressed_size)
@@ -732,10 +775,18 @@ pub fn build_chunked_data_from_precompressed(
                 offset_size,
                 element_size as u32,
             )
-    } else {
+        }
+        ChunkIndexPlan::FixedArray(grid, nslots) => {
             let fa_address = base_address + data_buf.len() as u64;
-        let fa_bytes = build_fixed_array_at(
+            let slots = index_slots(
+                grid,
+                &pre.shape,
+                &pre.chunk_dims,
                 &written_chunks,
+                Some(*nslots),
+            )?;
+            let fa_bytes = build_fixed_array_at(
+                &slots,
                 offset_size,
                 length_size,
                 pre.has_filters,
@@ -747,15 +798,263 @@ pub fn build_chunked_data_from_precompressed(
                 fa_address,
                 offset_size,
                 element_size as u32,
-            10, // max_nelmts_bits — matches h5py convention
+                FA_PAGE_BITS,
             )
+        }
+        ChunkIndexPlan::BTreeV2 => {
+            let bt_address = base_address + data_buf.len() as u64;
+            let records: Vec<(Vec<u64>, &WrittenChunk)> = written_chunks
+                .iter()
+                .enumerate()
+                .map(|(i, c)| (scaled_coords(&pre.shape, &pre.chunk_dims, i), c))
+                .collect();
+            let (bt_bytes, node_size) = build_btree_v2_chunk_index_at(
+                pre.shape.len(),
+                &records,
+                offset_size,
+                length_size,
+                pre.has_filters,
+                bt_address,
+            )?;
+            data_buf.extend_from_slice(&bt_bytes);
+            serialize_v4_btree_v2(
+                &chunk_dims_u32,
+                bt_address,
+                offset_size,
+                element_size as u32,
+                node_size,
+            )
+        }
     };
 
-    ChunkedDataResult {
+    Ok(ChunkedDataResult {
         data_bytes: data_buf,
         layout_message,
         pipeline_message: pre.pipeline_message.clone(),
+    })
+}
+
+/// Most slots a Fixed Array index may have before we refuse to build it: its
+/// data block holds one element per chunk of the *maximum* extent, so a huge
+/// finite maxshape with small chunks would otherwise exhaust memory.
+const MAX_FIXED_ARRAY_SLOTS: u64 = 1 << 26;
+
+/// Which chunk index a dataset gets, following the library's choice in
+/// `H5D__layout_set_latest_indexing`: version-2 B-tree for more than one
+/// unlimited dimension, Extensible Array for exactly one, Fixed Array for a
+/// finite maxshape, Single Chunk when the whole maximum extent is one chunk.
+enum ChunkIndexPlan {
+    SingleChunk,
+    /// The grid and the number of array elements (chunks of the max extent).
+    FixedArray(ChunkGrid, usize),
+    ExtensibleArray(ChunkGrid),
+    BTreeV2,
+}
+
+impl ChunkIndexPlan {
+    fn new(
+        shape: &[u64],
+        maxshape: Option<&[u64]>,
+        chunk_dims: &[u64],
+    ) -> Result<Self, FormatError> {
+        let bad = |what: &str| FormatError::ChunkedReadError(format!("maxshape: {what}"));
+        if let Some(ms) = maxshape {
+            if ms.len() != shape.len() {
+                return Err(bad("rank differs from the shape"));
             }
+            if ms.iter().zip(shape).any(|(&m, &s)| m < s) {
+                return Err(bad("smaller than the shape"));
+            }
+        }
+        let max = maxshape.unwrap_or(shape);
+        let nunlim = max.iter().filter(|&&d| d == u64::MAX).count();
+        match nunlim {
+            0 => {
+                let nslots = max
+                    .iter()
+                    .zip(chunk_dims)
+                    .try_fold(1u64, |acc, (&m, &c)| acc.checked_mul(m.div_ceil(c.max(1))))
+                    .filter(|&n| n <= MAX_FIXED_ARRAY_SLOTS)
+                    .ok_or_else(|| {
+                        bad("too many chunks for a Fixed Array index; \
+                             use larger chunks or an unlimited dimension")
+                    })?;
+                // A Single Chunk index needs that one chunk to exist; an
+                // empty dataset gets an all-unallocated Fixed Array instead.
+                let empty = shape.contains(&0);
+                if nslots == 1 && !empty {
+                    Ok(Self::SingleChunk)
+                } else {
+                    let grid = ChunkGrid::fixed_array(shape, Some(max), chunk_dims)?;
+                    Ok(Self::FixedArray(grid, nslots as usize))
+                }
+            }
+            1 => Ok(Self::ExtensibleArray(ChunkGrid::extensible_array(
+                shape,
+                Some(max),
+                chunk_dims,
+            )?)),
+            _ => Ok(Self::BTreeV2),
+        }
+    }
+}
+
+/// Place each written chunk at its linear index in `grid`. `chunks` are in
+/// row-major order over the chunks of the current extent (`split_into_chunks`).
+/// `len` fixes the slot count (Fixed Array); otherwise it is one past the
+/// highest index used.
+fn index_slots(
+    grid: &ChunkGrid,
+    shape: &[u64],
+    chunk_dims: &[u64],
+    chunks: &[WrittenChunk],
+    len: Option<usize>,
+) -> Result<Vec<Option<WrittenChunk>>, FormatError> {
+    let mut placed: Vec<(usize, &WrittenChunk)> = Vec::with_capacity(chunks.len());
+    for (i, chunk) in chunks.iter().enumerate() {
+        let scaled = scaled_coords(shape, chunk_dims, i);
+        let idx = usize::try_from(grid.linear_index(&scaled))
+            .map_err(|_| FormatError::Overflow("chunk index slot".into()))?;
+        placed.push((idx, chunk));
+    }
+    let n = len.unwrap_or_else(|| placed.iter().map(|&(i, _)| i + 1).max().unwrap_or(0));
+    let mut slots = vec![None; n];
+    for (idx, chunk) in placed {
+        *slots
+            .get_mut(idx)
+            .ok_or_else(|| FormatError::Overflow("chunk index slot".into()))? = Some(chunk.clone());
+    }
+    Ok(slots)
+}
+
+/// Scaled coordinates (`offset / chunk_dim`) of the `i`-th chunk in the
+/// row-major order `split_into_chunks` produces over the current extent.
+fn scaled_coords(shape: &[u64], chunk_dims: &[u64], i: usize) -> Vec<u64> {
+    let rank = shape.len();
+    let mut scaled = vec![0u64; rank];
+    let mut rem = i as u64;
+    for d in (0..rank).rev() {
+        let n = shape[d].div_ceil(chunk_dims[d]);
+        scaled[d] = rem % n;
+        rem /= n;
+    }
+    scaled
+}
+
+/// Node size the library gives a chunk index B-tree (`H5D_BT2_NODE_SIZE`),
+/// with its split and merge percentages.
+const BT2_NODE_SIZE: u32 = 2048;
+const BT2_SPLIT_PERCENT: u8 = 100;
+const BT2_MERGE_PERCENT: u8 = 40;
+/// B-tree v2 record types for chunk indexes (`H5B2_CDSET_ID`,
+/// `H5B2_CDSET_FILT_ID`).
+const BT2_CHUNK_UNFILTERED: u8 = 10;
+const BT2_CHUNK_FILTERED: u8 = 11;
+
+/// Build a version-2 B-tree chunk index (the library's index for datasets
+/// with more than one unlimited dimension) at a known absolute address.
+///
+/// `records` are `(scaled coordinates, chunk)` in lexicographic order of the
+/// coordinates, which is the order the library's comparator
+/// (`H5VM_vector_cmp_u`) keeps them in. The tree is a single leaf: the
+/// library's 2048-byte node when the records fit, otherwise a leaf node
+/// sized to hold them all (the root's record count is 16-bit, so at most
+/// 65535 chunks). Returns the bytes and the node size the layout message
+/// must record.
+fn build_btree_v2_chunk_index_at(
+    rank: usize,
+    records: &[(Vec<u64>, &WrittenChunk)],
+    offset_size: u8,
+    length_size: u8,
+    has_filters: bool,
+    base_address: u64,
+) -> Result<(Vec<u8>, u32), FormatError> {
+    let os = offset_size as usize;
+    let nrec = u16::try_from(records.len()).map_err(|_| {
+        FormatError::ChunkedReadError(
+            "more than 65535 chunks with more than one unlimited dimension: \
+             use larger chunks"
+                .into(),
+        )
+    })?;
+    let chunk_size_bytes = has_filters.then(|| {
+        let slots: Vec<Option<WrittenChunk>> =
+            records.iter().map(|(_, c)| Some((*c).clone())).collect();
+        filtered_chunk_size_len(&slots)
+    });
+    let record_size = os + chunk_size_bytes.map_or(0, |n| n + 4) + 8 * rank;
+    // Leaf: signature, version, type, records, checksum.
+    let leaf_len = 4 + 1 + 1 + records.len() * record_size + 4;
+    let node_size = u32::try_from(leaf_len)
+        .map_err(|_| FormatError::Overflow("B-tree v2 leaf size".into()))?
+        .max(BT2_NODE_SIZE);
+    let tree_type = if has_filters {
+        BT2_CHUNK_FILTERED
+    } else {
+        BT2_CHUNK_UNFILTERED
+    };
+
+    let hdr_len = 4 + 1 + 1 + 4 + 2 + 2 + 1 + 1 + os + 2 + length_size as usize + 4;
+    let leaf_address = base_address + hdr_len as u64;
+
+    let mut out = Vec::with_capacity(hdr_len + node_size as usize);
+    out.extend_from_slice(b"BTHD");
+    out.push(0); // version
+    out.push(tree_type);
+    out.extend_from_slice(&node_size.to_le_bytes());
+    out.extend_from_slice(&(record_size as u16).to_le_bytes());
+    out.extend_from_slice(&0u16.to_le_bytes()); // depth
+    out.push(BT2_SPLIT_PERCENT);
+    out.push(BT2_MERGE_PERCENT);
+    if records.is_empty() {
+        out.extend(core::iter::repeat_n(0xFF, os));
+    } else {
+        push_addr(&mut out, leaf_address, offset_size);
+    }
+    out.extend_from_slice(&nrec.to_le_bytes());
+    match length_size {
+        4 => out.extend_from_slice(&(records.len() as u32).to_le_bytes()),
+        _ => out.extend_from_slice(&(records.len() as u64).to_le_bytes()),
+    }
+    let sum = jenkins_lookup3(&out);
+    out.extend_from_slice(&sum.to_le_bytes());
+    debug_assert_eq!(out.len(), hdr_len);
+    if records.is_empty() {
+        return Ok((out, node_size));
+    }
+
+    let leaf_start = out.len();
+    out.extend_from_slice(b"BTLF");
+    out.push(0); // version
+    out.push(tree_type);
+    for (scaled, chunk) in records {
+        push_index_element(&mut out, Some(chunk), offset_size, chunk_size_bytes);
+        for &c in scaled {
+            out.extend_from_slice(&c.to_le_bytes());
+        }
+    }
+    let sum = jenkins_lookup3(&out[leaf_start..]);
+    out.extend_from_slice(&sum.to_le_bytes());
+    // The library reads whole nodes; pad the leaf out to the node size.
+    out.resize(leaf_start + node_size as usize, 0);
+    Ok((out, node_size))
+}
+
+/// Serialize a v4 layout message for a version-2 B-tree chunk index.
+fn serialize_v4_btree_v2(
+    chunk_dims: &[u32],
+    btree_address: u64,
+    offset_size: u8,
+    element_size: u32,
+    node_size: u32,
+) -> Vec<u8> {
+    let mut buf = layout_v4_chunked_prefix(chunk_dims, element_size);
+    buf.push(5); // chunk index type = 5 (version-2 B-tree)
+    buf.extend_from_slice(&node_size.to_le_bytes());
+    buf.push(BT2_SPLIT_PERCENT);
+    buf.push(BT2_MERGE_PERCENT);
+    push_addr(&mut buf, btree_address, offset_size);
+    buf
 }
 
 /// Build chunked data with absolute addresses.
@@ -790,11 +1089,7 @@ pub fn build_chunked_data_at_ext(
     maxshape: Option<&[u64]>,
 ) -> Result<ChunkedDataResult, FormatError> {
     let pre = precompress_chunks(raw_data, shape, chunk_dims, element_size, options)?;
-    Ok(build_chunked_data_from_precompressed(
+    build_chunked_data_from_precompressed(&pre, base_address, maxshape)
-        &pre,
-        base_address,
-        maxshape,
-    ))
 }
 
 /// Write selected elements into an existing in-memory dataset buffer.
@@ -1382,7 +1677,8 @@ mod tests {
                 filter_mask: 0,
             },
         ];
-        let fa = build_fixed_array_at(&chunks, 8, 8, false, 0x2000);
+        let slots: Vec<_> = chunks.into_iter().map(Some).collect();
+        let fa = build_fixed_array_at(&slots, 8, 8, false, 0x2000);
         // Should start with FAHD
         assert_eq!(&fa[0..4], b"FAHD");
         // FAHD size = 4+1+1+1+1+8+8+4 = 28
@@ -1429,7 +1725,8 @@ mod tests {
                 filter_mask: 0,
             },
         ];
-        let ea = ea_writer::build_extensible_array_at(&chunks, 8, 8, false, 0x2000);
+        let slots: Vec<_> = chunks.into_iter().map(Some).collect();
+        let ea = ea_writer::build_extensible_array_at(&slots, 8, 8, false, 0x2000);
         assert_eq!(&ea[0..4], b"EAHD");
         // Find EAIB after EAHD: 12 fixed + 6*8 stats + 8 addr + 4 checksum = 72
         let aehd_size = 4 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 6 * 8 + 8 + 4;

crates/clawhdf5-format/src/ea_writer.rs

@@ -7,7 +7,7 @@ extern crate alloc;
 use alloc::{vec, vec::Vec};
 
 use crate::checksum::jenkins_lookup3;
-use crate::chunked_write::WrittenChunk;
+use crate::chunked_write::{WrittenChunk, filtered_chunk_size_len, push_addr, push_index_element};
 
 /// Serialize a v4 Extensible Array layout message.
 pub(crate) fn serialize_v4_extensible_array(
@@ -58,11 +58,11 @@ pub(crate) fn serialize_v4_extensible_array(
     buf.push(4);
 
     // EA creation parameters (must match AEHD and HDF5 C library defaults)
-    buf.push(32); // max_nelmts_bits
+    buf.push(MAX_NELMTS_BITS);
-    buf.push(4); // idx_blk_elmts
+    buf.push(IDX_BLK_ELMTS);
-    buf.push(4); // super_blk_min_data_ptrs
+    buf.push(SUP_BLK_MIN_DATA_PTRS);
-    buf.push(16); // data_blk_min_elmts
+    buf.push(DATA_BLK_MIN_ELMTS);
-    buf.push(10); // max_dblk_page_nelmts_bits
+    buf.push(MAX_DBLK_PAGE_NELMTS_BITS);
 
     // EA header address
     match offset_size {
@@ -74,304 +74,281 @@ pub(crate) fn serialize_v4_extensible_array(
     buf
 }
 
+// EA creation parameters — the HDF5 library's defaults for chunk indexes
+// (`H5D_EARRAY_*`); the layout message above and the header must agree.
+const MAX_NELMTS_BITS: u8 = 32;
+const IDX_BLK_ELMTS: u8 = 4;
+const SUP_BLK_MIN_DATA_PTRS: u8 = 4;
+const DATA_BLK_MIN_ELMTS: u8 = 16;
+const MAX_DBLK_PAGE_NELMTS_BITS: u8 = 10;
+
+/// One data block of the array: its first element (relative to the end of
+/// the index block's own elements), element count, and address when it is
+/// allocated.
+struct DataBlock {
+    start: usize,
+    nelmts: usize,
+    addr: Option<u64>,
+}
+
 /// Build a complete Extensible Array at a known absolute address.
 ///
-/// For simplicity, we put all elements inline in the index block when the
+/// `slots[i]` is the element at linear index `i` (see `chunk_grid`); `None`
-/// number of chunks is small (up to idx_blk_elmts), otherwise use inline +
+/// marks an unallocated chunk. The first `IDX_BLK_ELMTS` elements live in
-/// direct data blocks.
+/// the index block, the rest in data blocks grouped by super block level
+/// exactly as `H5EA__hdr_init` sizes them: level `u` has `2^(u/2)` data
+/// blocks of `DATA_BLK_MIN_ELMTS * 2^ceil(u/2)` elements. The data blocks of
+/// the first levels are addressed straight from the index block; later
+/// levels go through a super block (EASB). Data blocks larger than a page
+/// (`2^MAX_DBLK_PAGE_NELMTS_BITS` elements) are paged, with their page-init
+/// bits kept in the owning super block. Only blocks holding a defined element
+/// are allocated; the rest keep the undefined address, as in a file the
+/// library wrote.
 pub fn build_extensible_array_at(
-    chunks: &[WrittenChunk],
+    slots: &[Option<WrittenChunk>],
     offset_size: u8,
     length_size: u8,
     has_filters: bool,
     ea_base_address: u64,
 ) -> Vec<u8> {
     let os = offset_size as usize;
-    let num_elements = chunks.len();
+    let chunk_size_bytes = has_filters.then(|| filtered_chunk_size_len(slots));
-
+    let elem_size = os + chunk_size_bytes.map_or(0, |n| n + 4);
-    // Compute element encoding size (same logic as Fixed Array)
-    let chunk_size_bytes: usize = if has_filters {
-        let max_raw = chunks.iter().map(|c| c.raw_size).max().unwrap_or(1);
-        let log2_val = if max_raw <= 1 {
-            0
-        } else {
-            63 - max_raw.leading_zeros()
-        };
-        let len = 1 + ((log2_val + 8) / 8) as usize;
-        len.min(8)
-    } else {
-        0
-    };
-
-    let elem_size = if has_filters {
-        os + chunk_size_bytes + 4
-    } else {
-        os
-    };
-
     let client_id: u8 = if has_filters { 1 } else { 0 };
+    let arr_off_size = (MAX_NELMTS_BITS as usize).div_ceil(8);
+    let page_nelmts = 1usize << MAX_DBLK_PAGE_NELMTS_BITS;
+    let idx_blk = IDX_BLK_ELMTS as usize;
 
-    // EA creation parameters — must match HDF5 C library defaults exactly
+    // Elements past the last defined one are never realised
-    let max_nelmts_bits: u8 = 32;
+    // (`max_idx_set` is one past the highest index ever set).
-    let idx_blk_elmts: u8 = 4;
+    let max_idx_set = slots.iter().rposition(Option::is_some).map_or(0, |i| i + 1);
-    let min_dblk_nelmts: u8 = 16;
+    let slots = &slots[..max_idx_set];
-    let super_blk_min_nelmts: u8 = 4;
+    let defined_in = |start: usize, n: usize| -> bool {
-    let max_dblk_nelmts_bits: u8 = 10;
+        let lo = idx_blk.saturating_add(start).min(slots.len());
+        let hi = idx_blk
+            .saturating_add(start)
+            .saturating_add(n)
+            .min(slots.len());
+        slots[lo..hi].iter().any(Option::is_some)
+    };
 
-    // EAHD size: fixed(12) + 6 stats(6*length_size) + addr(offset_size) + checksum(4)
+    // Super block levels: (ndblks, dblk_nelmts, first element).
+    let log2_dmin = (DATA_BLK_MIN_ELMTS as u32).trailing_zeros() as usize;
+    let nsblks = 1 + MAX_NELMTS_BITS as usize - log2_dmin;
+    let ndblk_addrs = 2 * (SUP_BLK_MIN_DATA_PTRS as usize - 1);
+    let mut levels: Vec<(usize, usize, usize)> = Vec::with_capacity(nsblks);
+    let mut start = 0usize;
+    for u in 0..nsblks {
+        let ndblks = 1usize << (u / 2);
+        let nelmts = (DATA_BLK_MIN_ELMTS as usize) << u.div_ceil(2);
+        levels.push((ndblks, nelmts, start));
+        // Saturate: on 32-bit targets the last levels only need to compare
+        // as "beyond the end".
+        start = start.saturating_add(ndblks.saturating_mul(nelmts));
+    }
+    // Levels whose data blocks the index block addresses directly.
+    let mut direct_levels = 0;
+    let mut n = 0;
+    while n < ndblk_addrs {
+        n += levels[direct_levels].0;
+        direct_levels += 1;
+    }
+    let nsblk_addrs = nsblks - direct_levels;
+
+    let dblk_size = |nelmts: usize| -> usize {
+        let prefix = 4 + 1 + 1 + os + arr_off_size + 4;
+        if nelmts > page_nelmts {
+            prefix + (nelmts / page_nelmts) * (page_nelmts * elem_size + 4)
+        } else {
+            prefix + nelmts * elem_size
+        }
+    };
+    let sblk_bitmap_len = |ndblks: usize, nelmts: usize| -> usize {
+        if nelmts > page_nelmts {
+            ndblks * (nelmts / page_nelmts).div_ceil(8)
+        } else {
+            0
+        }
+    };
+
+    // Plan addresses: header, index block, the direct data blocks, then each
+    // allocated super block followed by its allocated data blocks.
     let aehd_size = 4 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 1 + 6 * length_size as usize + os + 4;
     let aeib_address = ea_base_address + aehd_size as u64;
+    let aeib_size = 4 + 1 + 1 + os + idx_blk * elem_size + ndblk_addrs * os + nsblk_addrs * os + 4;
+    let mut cursor = aeib_address + aeib_size as u64;
 
-    // Determine how many elements go inline vs data blocks
+    let mut ndata_blks = 0u64;
-    let n_inline = (idx_blk_elmts as usize).min(num_elements);
+    let mut data_blk_size = 0u64;
-    let remaining_after_inline = num_elements.saturating_sub(n_inline);
+    let mut nsuper_blks = 0u64;
+    let mut super_blk_size = 0u64;
+    let mut realized = idx_blk as u64;
 
-    // Compute super block layout per HDF5 spec
+    let mut plan_dblk = |cursor: &mut u64, start: usize, nelmts: usize| -> DataBlock {
-    let sblk_min = super_blk_min_nelmts as usize;
+        let addr = defined_in(start, nelmts).then(|| {
-    let log2_dblk_min = if min_dblk_nelmts <= 1 {
+            let a = *cursor;
-        0
+            let size = dblk_size(nelmts) as u64;
-    } else {
+            *cursor += size;
-        (min_dblk_nelmts as u32).trailing_zeros() as usize
+            ndata_blks += 1;
-    };
+            data_blk_size += size;
-    let nsblks = (max_nelmts_bits as usize).saturating_sub(log2_dblk_min) + 1;
+            realized += nelmts as u64;
-
+            a
-    // Direct data block addresses (from super blocks 0..sblk_min-1)
+        });
-    let mut dblk_sizes: Vec<usize> = Vec::new();
+        DataBlock {
-    for sblk_idx in 0..sblk_min.min(nsblks) {
+            start,
-        let ndblks = 1usize << (sblk_idx / 2);
+            nelmts,
-        let dblk_nelmts = (min_dblk_nelmts as usize) * (1 << sblk_idx.div_ceil(2));
+            addr,
-        for _ in 0..ndblks {
-            dblk_sizes.push(dblk_nelmts);
         }
-    }
-    let n_direct_dblks = dblk_sizes.len();
-
-    // Super block addresses (for super blocks sblk_min..nsblks-1)
-    let n_sblk_addrs = nsblks.saturating_sub(sblk_min);
-
-    // EAIB size
-    let aeib_size = 4
-        + 1
-        + 1
-        + os
-        + idx_blk_elmts as usize * elem_size
-        + n_direct_dblks * os
-        + n_sblk_addrs * os
-        + 4;
-
-    // Build AEHD
-    let mut aehd = Vec::with_capacity(aehd_size);
-    aehd.extend_from_slice(b"EAHD");
-    aehd.push(0); // version
-    aehd.push(client_id);
-    aehd.push(elem_size as u8);
-    aehd.push(max_nelmts_bits);
-    aehd.push(idx_blk_elmts);
-    aehd.push(min_dblk_nelmts);
-    aehd.push(super_blk_min_nelmts);
-    aehd.push(max_dblk_nelmts_bits);
-
-    // Count data blocks that will have chunks
-    let n_active_dblks: u64 = if remaining_after_inline > 0 {
-        let mut count = 0u64;
-        let mut ci = n_inline;
-        for &sz in &dblk_sizes {
-            if ci < num_elements {
-                count += 1;
-                ci += sz;
-            }
-        }
-        count
-    } else {
-        0
-    };
-    let blk_off_size = (max_nelmts_bits as usize).div_ceil(8);
-    let aedb_header_overhead = 4 + 1 + 1 + os + blk_off_size + 4;
-    let data_blk_total_size: u64 = if remaining_after_inline > 0 {
-        let mut total = 0u64;
-        let mut ci = n_inline;
-        for &sz in &dblk_sizes {
-            if ci < num_elements {
-                total += (aedb_header_overhead + sz * elem_size) as u64;
-                ci += sz;
-            }
-        }
-        total
-    } else {
-        0
-    };
-    let max_idx_set: u64 = if remaining_after_inline > 0 {
-        let mut max_set = idx_blk_elmts as u64;
-        let mut ci = n_inline;
-        for &sz in &dblk_sizes {
-            if ci < num_elements {
-                max_set += sz as u64;
-                ci += sz;
-            }
-        }
-        max_set
-    } else {
-        idx_blk_elmts as u64
     };
 
+    let mut direct: Vec<DataBlock> = Vec::with_capacity(ndblk_addrs);
+    for &(ndblks, nelmts, first) in &levels[..direct_levels] {
+        for k in 0..ndblks {
+            direct.push(plan_dblk(&mut cursor, first + k * nelmts, nelmts));
+        }
+    }
+    // (super block address, level, its data blocks)
+    let mut supers: Vec<(Option<u64>, usize, Vec<DataBlock>)> = Vec::with_capacity(nsblk_addrs);
+    for (u, &(ndblks, nelmts, first)) in levels.iter().enumerate().skip(direct_levels) {
+        if !defined_in(first, ndblks.saturating_mul(nelmts)) {
+            supers.push((None, u, Vec::new()));
+            continue;
+        }
+        let sb_size =
+            4 + 1 + 1 + os + arr_off_size + sblk_bitmap_len(ndblks, nelmts) + ndblks * os + 4;
+        let sb_addr = cursor;
+        cursor += sb_size as u64;
+        nsuper_blks += 1;
+        super_blk_size += sb_size as u64;
+        let dblks = (0..ndblks)
+            .map(|k| plan_dblk(&mut cursor, first + k * nelmts, nelmts))
+            .collect();
+        supers.push((Some(sb_addr), u, dblks));
+    }
+
+    let slot = |i: usize| slots.get(i).and_then(Option::as_ref);
     let write_length = |buf: &mut Vec<u8>, val: u64| match length_size {
         4 => buf.extend_from_slice(&(val as u32).to_le_bytes()),
         _ => buf.extend_from_slice(&val.to_le_bytes()),
     };
-    let write_addr = |buf: &mut Vec<u8>, val: u64| match offset_size {
+    let write_addr_opt = |buf: &mut Vec<u8>, addr: Option<u64>| match addr {
-        4 => buf.extend_from_slice(&(val as u32).to_le_bytes()),
+        Some(a) => push_addr(buf, a, offset_size),
-        _ => buf.extend_from_slice(&val.to_le_bytes()),
+        None => buf.extend(core::iter::repeat_n(0xFF, os)),
+    };
+    let block_prefix = |buf: &mut Vec<u8>, sig: &[u8; 4], block_off: usize| {
+        buf.extend_from_slice(sig);
+        buf.push(0); // version
+        buf.push(client_id);
+        push_addr(buf, ea_base_address, offset_size);
+        buf.extend_from_slice(&(block_off as u64).to_le_bytes()[..arr_off_size]);
+    };
+    // Serialise one data block (paged or not) onto `out`.
+    let write_dblk = |out: &mut Vec<u8>, db: &DataBlock| {
+        let at = out.len();
+        block_prefix(out, b"EADB", db.start);
+        let first = idx_blk + db.start;
+        if db.nelmts > page_nelmts {
+            // Paged: the prefix carries only its own checksum; each page
+            // follows with one of its own.
+            let sum = jenkins_lookup3(&out[at..]);
+            out.extend_from_slice(&sum.to_le_bytes());
+            for p in 0..db.nelmts / page_nelmts {
+                let page_at = out.len();
+                for e in 0..page_nelmts {
+                    let i = first + p * page_nelmts + e;
+                    push_index_element(out, slot(i), offset_size, chunk_size_bytes);
+                }
+                let sum = jenkins_lookup3(&out[page_at..]);
+                out.extend_from_slice(&sum.to_le_bytes());
+            }
+        } else {
+            for i in first..first + db.nelmts {
+                push_index_element(out, slot(i), offset_size, chunk_size_bytes);
+            }
+            let sum = jenkins_lookup3(&out[at..]);
+            out.extend_from_slice(&sum.to_le_bytes());
+        }
+        debug_assert_eq!(out.len() - at, dblk_size(db.nelmts));
     };
 
-    write_length(&mut aehd, 0);
+    // Header (EAHD). The six statistics are, in order: super blocks, their
-    write_length(&mut aehd, 0);
+    // bytes, data blocks, their bytes, max index set, elements realised.
-    write_length(&mut aehd, n_active_dblks);
+    let mut out = Vec::with_capacity((cursor - ea_base_address) as usize);
-    write_length(&mut aehd, data_blk_total_size);
+    out.extend_from_slice(b"EAHD");
-    write_length(&mut aehd, num_elements as u64);
+    out.push(0); // version
-    write_length(&mut aehd, max_idx_set);
+    out.push(client_id);
+    out.push(elem_size as u8);
+    out.push(MAX_NELMTS_BITS);
+    out.push(IDX_BLK_ELMTS);
+    out.push(DATA_BLK_MIN_ELMTS);
+    out.push(SUP_BLK_MIN_DATA_PTRS);
+    out.push(MAX_DBLK_PAGE_NELMTS_BITS);
+    write_length(&mut out, nsuper_blks);
+    write_length(&mut out, super_blk_size);
+    write_length(&mut out, ndata_blks);
+    write_length(&mut out, data_blk_size);
+    write_length(&mut out, max_idx_set as u64);
+    write_length(&mut out, realized);
+    push_addr(&mut out, aeib_address, offset_size);
+    let sum = jenkins_lookup3(&out);
+    out.extend_from_slice(&sum.to_le_bytes());
+    debug_assert_eq!(out.len(), aehd_size);
 
-    write_addr(&mut aehd, aeib_address);
+    // Index block (EAIB): inline elements, data block and super block
-
+    // addresses.
-    let aehd_checksum = jenkins_lookup3(&aehd);
+    let ib_start = out.len();
-    aehd.extend_from_slice(&aehd_checksum.to_le_bytes());
+    out.extend_from_slice(b"EAIB");
-    debug_assert_eq!(aehd.len(), aehd_size);
+    out.push(0);
-
+    out.push(client_id);
-    // Build AEIB
+    push_addr(&mut out, ea_base_address, offset_size);
-    let mut aeib = Vec::with_capacity(aeib_size);
+    for i in 0..idx_blk {
-    aeib.extend_from_slice(b"EAIB");
+        push_index_element(&mut out, slot(i), offset_size, chunk_size_bytes);
-    aeib.push(0);
-    aeib.push(client_id);
-
-    match offset_size {
-        4 => aeib.extend_from_slice(&(ea_base_address as u32).to_le_bytes()),
-        8 => aeib.extend_from_slice(&ea_base_address.to_le_bytes()),
-        _ => aeib.extend_from_slice(&ea_base_address.to_le_bytes()),
     }
-
+    for db in &direct {
-    // Inline elements
+        write_addr_opt(&mut out, db.addr);
-    #[allow(clippy::needless_range_loop)]
-    for i in 0..idx_blk_elmts as usize {
-        if i < n_inline {
-            write_chunk_element(
-                &mut aeib,
-                &chunks[i],
-                offset_size,
-                has_filters,
-                chunk_size_bytes,
-            );
-        } else {
-            write_undefined_element(&mut aeib, offset_size, has_filters, chunk_size_bytes);
     }
+    for (sb_addr, _, _) in &supers {
+        write_addr_opt(&mut out, *sb_addr);
     }
+    let sum = jenkins_lookup3(&out[ib_start..]);
+    out.extend_from_slice(&sum.to_le_bytes());
+    debug_assert_eq!(out.len() - ib_start, aeib_size);
 
-    // Data block addresses + build data blocks
+    for db in direct.iter().filter(|d| d.addr.is_some()) {
-    let mut data_blocks_buf = Vec::new();
+        write_dblk(&mut out, db);
-    let dblks_base = aeib_address + aeib_size as u64;
-    let mut dblk_cursor = dblks_base;
-    let mut chunk_idx = n_inline;
-
-    for &nelmts in &dblk_sizes {
-        if chunk_idx >= num_elements {
-            match offset_size {
-                4 => aeib.extend_from_slice(&u32::MAX.to_le_bytes()),
-                8 => aeib.extend_from_slice(&u64::MAX.to_le_bytes()),
-                _ => aeib.extend_from_slice(&u64::MAX.to_le_bytes()),
     }
+    for (sb_addr, u, dblks) in &supers {
+        if sb_addr.is_none() {
             continue;
         }
-
+        let (ndblks, nelmts, first) = levels[*u];
-        match offset_size {
+        let sb_start = out.len();
-            4 => aeib.extend_from_slice(&(dblk_cursor as u32).to_le_bytes()),
+        block_prefix(&mut out, b"EASB", first);
-            8 => aeib.extend_from_slice(&dblk_cursor.to_le_bytes()),
+        if nelmts > page_nelmts {
-            _ => aeib.extend_from_slice(&dblk_cursor.to_le_bytes()),
+            // Page-init bits, `npages` per data block, packed MSB-first
-        }
+            // (`H5VM_bit_set`): every page of an allocated data block is
-
+            // written.
-        // Build EADB
+            let npages = nelmts / page_nelmts;
-        let mut aedb = Vec::new();
+            let mut bitmap = vec![0u8; sblk_bitmap_len(ndblks, nelmts)];
-        aedb.extend_from_slice(b"EADB");
+            for (k, db) in dblks.iter().enumerate() {
-        aedb.push(0);
+                if db.addr.is_some() {
-        aedb.push(client_id);
+                    for p in 0..npages {
-        match offset_size {
+                        let bit = k * npages + p;
-            4 => aedb.extend_from_slice(&(ea_base_address as u32).to_le_bytes()),
+                        bitmap[bit / 8] |= 0x80 >> (bit % 8);
-            8 => aedb.extend_from_slice(&ea_base_address.to_le_bytes()),
-            _ => aedb.extend_from_slice(&ea_base_address.to_le_bytes()),
-        }
-
-        let blk_off_size = (max_nelmts_bits as usize).div_ceil(8);
-        let blk_off_val = (chunk_idx - n_inline) as u64;
-        aedb.extend_from_slice(&blk_off_val.to_le_bytes()[..blk_off_size]);
-
-        for slot in 0..nelmts {
-            if chunk_idx + slot < num_elements {
-                write_chunk_element(
-                    &mut aedb,
-                    &chunks[chunk_idx + slot],
-                    offset_size,
-                    has_filters,
-                    chunk_size_bytes,
-                );
-            } else {
-                write_undefined_element(&mut aedb, offset_size, has_filters, chunk_size_bytes);
                     }
                 }
-
-        let aedb_checksum = jenkins_lookup3(&aedb);
-        aedb.extend_from_slice(&aedb_checksum.to_le_bytes());
-
-        dblk_cursor += aedb.len() as u64;
-        data_blocks_buf.extend_from_slice(&aedb);
-        chunk_idx += nelmts;
             }
-
+            out.extend_from_slice(&bitmap);
-    // Super block addresses (all undefined)
+        }
-    for _ in 0..n_sblk_addrs {
+        for db in dblks {
-        match offset_size {
+            write_addr_opt(&mut out, db.addr);
-            4 => aeib.extend_from_slice(&u32::MAX.to_le_bytes()),
+        }
-            8 => aeib.extend_from_slice(&u64::MAX.to_le_bytes()),
+        let sum = jenkins_lookup3(&out[sb_start..]);
-            _ => aeib.extend_from_slice(&u64::MAX.to_le_bytes()),
+        out.extend_from_slice(&sum.to_le_bytes());
+        for db in dblks.iter().filter(|d| d.addr.is_some()) {
+            write_dblk(&mut out, db);
         }
     }
-
+    debug_assert_eq!(out.len() as u64, cursor - ea_base_address);
-    let aeib_checksum = jenkins_lookup3(&aeib);
+    out
-    aeib.extend_from_slice(&aeib_checksum.to_le_bytes());
-    debug_assert_eq!(aeib.len(), aeib_size);
-
-    let mut combined = aehd;
-    combined.extend_from_slice(&aeib);
-    combined.extend_from_slice(&data_blocks_buf);
-    combined
-}
-
-fn write_chunk_element(
-    buf: &mut Vec<u8>,
-    chunk: &WrittenChunk,
-    offset_size: u8,
-    has_filters: bool,
-    chunk_size_bytes: usize,
-) {
-    match offset_size {
-        4 => buf.extend_from_slice(&(chunk.address as u32).to_le_bytes()),
-        8 => buf.extend_from_slice(&chunk.address.to_le_bytes()),
-        _ => buf.extend_from_slice(&chunk.address.to_le_bytes()),
-    }
-    if has_filters {
-        let cs_bytes = chunk.compressed_size.to_le_bytes();
-        buf.extend_from_slice(&cs_bytes[..chunk_size_bytes]);
-        buf.extend_from_slice(&chunk.filter_mask.to_le_bytes());
-    }
-}
-
-fn write_undefined_element(
-    buf: &mut Vec<u8>,
-    offset_size: u8,
-    has_filters: bool,
-    chunk_size_bytes: usize,
-) {
-    let os = offset_size as usize;
-    // Use extend with repeat to avoid heap-allocating a temporary Vec on each call.
-    buf.extend(core::iter::repeat_n(0xFF, os));
-    if has_filters {
-        buf.extend(core::iter::repeat_n(0x00, chunk_size_bytes));
-        buf.extend_from_slice(&0u32.to_le_bytes());
-    }
 }

crates/clawhdf5-format/src/extensible_array.rs

@@ -9,6 +9,7 @@ extern crate alloc;
 #[cfg(not(feature = "std"))]
 use alloc::{format, vec, vec::Vec};
 
+use crate::chunk_grid::ChunkGrid;
 use crate::chunked_read::ChunkInfo;
 use crate::error::FormatError;
 
@@ -203,8 +204,7 @@ fn read_element(
     offset_size: u8,
     chunk_byte_size: u64,
     linear_index: usize,
-    num_chunks_per_dim: &[u64],
+    grid: &ChunkGrid,
-    chunk_dimensions: &[u32],
 ) -> Result<(Option<ChunkInfo>, usize), FormatError> {
     let os = offset_size as usize;
 
@@ -220,7 +220,10 @@ fn read_element(
             return Ok((None, os));
         }
         let address = read_offset(data, pos, offset_size)?;
-        let offsets = index_to_chunk_offsets(linear_index, num_chunks_per_dim, chunk_dimensions);
+        // A slot beyond the current extent is ignored, as the library does.
+        let Some(offsets) = grid.offsets(linear_index as u64) else {
+            return Ok((None, os));
+        };
         Ok((
             Some(ChunkInfo {
                 chunk_size: chunk_byte_size as u32,
@@ -261,7 +264,9 @@ fn read_element(
             data[fm_off + 2],
             data[fm_off + 3],
         ]);
-        let offsets = index_to_chunk_offsets(linear_index, num_chunks_per_dim, chunk_dimensions);
+        let Some(offsets) = grid.offsets(linear_index as u64) else {
+            return Ok((None, elem_total));
+        };
         Ok((
             Some(ChunkInfo {
                 chunk_size: chunk_size as u32,
@@ -274,27 +279,6 @@ fn read_element(
     }
 }
 
-/// Convert a linear chunk index to N-dimensional chunk offsets in dataset space.
-fn index_to_chunk_offsets(
-    index: usize,
-    num_chunks_per_dim: &[u64],
-    chunk_dimensions: &[u32],
-) -> Vec<u64> {
-    let rank = num_chunks_per_dim.len();
-    let mut offsets = vec![0u64; rank];
-    let mut remaining = index as u64;
-    for d in (0..rank).rev() {
-        let nchunks = num_chunks_per_dim[d];
-        if nchunks == 0 {
-            continue;
-        }
-        let chunk_idx = remaining % nchunks;
-        remaining /= nchunks;
-        offsets[d] = chunk_idx * chunk_dimensions[d] as u64;
-    }
-    offsets
-}
-
 /// Collect elements from a data block at the given offset.
 #[allow(clippy::too_many_arguments)]
 /// Layout of super block `u`, per the HDF5 spec: the number of data blocks it
@@ -339,8 +323,7 @@ fn read_data_block_elements(
     offset_size: u8,
     chunk_byte_size: u64,
     start_index: usize,
-    num_chunks_per_dim: &[u64],
+    grid: &ChunkGrid,
-    chunk_dimensions: &[u32],
     page_init: &[u8],
     first_page: usize,
 ) -> Result<Vec<ChunkInfo>, FormatError> {
@@ -376,8 +359,7 @@ fn read_data_block_elements(
                 offset_size,
                 chunk_byte_size,
                 first_index + i,
-                num_chunks_per_dim,
+                grid,
-                chunk_dimensions,
             )?;
             if let Some(ci) = info {
                 chunks.push(ci);
@@ -449,25 +431,19 @@ pub fn read_extensible_array_chunks(
     file_data: &[u8],
     header: &ExtensibleArrayHeader,
     dataset_dims: &[u64],
+    max_dims: Option<&[u64]>,
     chunk_dimensions: &[u32],
     element_size: u32,
     offset_size: u8,
     _length_size: u8,
 ) -> Result<Vec<ChunkInfo>, FormatError> {
-    let rank = chunk_dimensions.len();
     let os = offset_size as usize;
 
-    let mut num_chunks_per_dim = Vec::with_capacity(rank);
+    // Linear indexes follow the maximum dimensions, with the unlimited
-    for d in 0..rank {
+    // dimension swizzled to the slowest position (see `chunk_grid`).
-        let ch_dim = chunk_dimensions[d] as u64;
+    let dims_u64: Vec<u64> = chunk_dimensions.iter().map(|&d| d as u64).collect();
-        if ch_dim == 0 {
+    let grid = ChunkGrid::extensible_array(dataset_dims, max_dims, &dims_u64)?;
-            return Err(FormatError::ChunkedReadError(
+    let grid = &grid;
-                "chunk dimension is zero".into(),
-            ));
-        }
-        let ds_dim = dataset_dims[d];
-        num_chunks_per_dim.push(ds_dim.div_ceil(ch_dim));
-    }
 
     let chunk_byte_size: u64 =
         chunk_dimensions.iter().map(|&d| d as u64).product::<u64>() * element_size as u64;
@@ -557,8 +533,7 @@ pub fn read_extensible_array_chunks(
             offset_size,
             chunk_byte_size,
             i,
-            &num_chunks_per_dim,
+            grid,
-            chunk_dimensions,
         )?;
         if let Some(ci) = info {
             chunks.push(ci);
@@ -594,8 +569,7 @@ pub fn read_extensible_array_chunks(
                 offset_size,
                 chunk_byte_size,
                 global_index,
-                &num_chunks_per_dim,
+                grid,
-                chunk_dimensions,
                 &[],
                 0,
             )?);
@@ -625,8 +599,7 @@ pub fn read_extensible_array_chunks(
                 offset_size,
                 chunk_byte_size,
                 global_index,
-                &num_chunks_per_dim,
+                grid,
-                chunk_dimensions,
             )?);
         }
         global_index =
@@ -653,8 +626,7 @@ fn read_super_block(
     offset_size: u8,
     chunk_byte_size: u64,
     start_index: usize,
-    num_chunks_per_dim: &[u64],
+    grid: &ChunkGrid,
-    chunk_dimensions: &[u32],
 ) -> Result<Vec<ChunkInfo>, FormatError> {
     let os = offset_size as usize;
     let sb_header_size = 4 + 1 + 1 + os + arr_off_size(header);
@@ -710,8 +682,7 @@ fn read_super_block(
                 offset_size,
                 chunk_byte_size,
                 global_idx,
-                num_chunks_per_dim,
+                grid,
-                chunk_dimensions,
                 bitmap,
                 i * npages,
             )?);
@@ -735,35 +706,18 @@ mod tests {
     }
     #[test]
     fn index_to_offsets_1d() {
-        let num_chunks = vec![5u64];
+        let g = ChunkGrid::fixed_array(&[100], None, &[20]).unwrap();
-        let chunk_dims = vec![20u32];
+        assert_eq!(g.offsets(0).unwrap(), vec![0]);
-        assert_eq!(index_to_chunk_offsets(0, &num_chunks, &chunk_dims), vec![0]);
+        assert_eq!(g.offsets(1).unwrap(), vec![20]);
-        assert_eq!(
+        assert_eq!(g.offsets(4).unwrap(), vec![80]);
-            index_to_chunk_offsets(1, &num_chunks, &chunk_dims),
-            vec![20]
-        );
-        assert_eq!(
-            index_to_chunk_offsets(4, &num_chunks, &chunk_dims),
-            vec![80]
-        );
     }
 
     #[test]
     fn index_to_offsets_2d() {
-        let num_chunks = vec![3u64, 2];
+        let g = ChunkGrid::fixed_array(&[10, 6], None, &[4, 3]).unwrap();
-        let chunk_dims = vec![4u32, 3];
+        assert_eq!(g.offsets(0).unwrap(), vec![0, 0]);
-        assert_eq!(
+        assert_eq!(g.offsets(1).unwrap(), vec![0, 3]);
-            index_to_chunk_offsets(0, &num_chunks, &chunk_dims),
+        assert_eq!(g.offsets(2).unwrap(), vec![4, 0]);
-            vec![0, 0]
-        );
-        assert_eq!(
-            index_to_chunk_offsets(1, &num_chunks, &chunk_dims),
-            vec![0, 3]
-        );
-        assert_eq!(
-            index_to_chunk_offsets(2, &num_chunks, &chunk_dims),
-            vec![4, 0]
-        );
     }
 
     #[test]
@@ -830,7 +784,7 @@ mod tests {
             index_block_address: (usize::MAX - 4) as u64,
         };
         let buf = vec![0u8; 64];
-        let r = read_extensible_array_chunks(&buf, &header, &[100], &[20], 8, 8, 8);
+        let r = read_extensible_array_chunks(&buf, &header, &[100], None, &[20], 8, 8, 8);
         assert!(r.is_err());
     }
 
@@ -913,8 +867,16 @@ mod tests {
         let header = ExtensibleArrayHeader::parse(&file_data, aehd_offset, os, ls).unwrap();
         let ds_dims = vec![40u64]; // 2 chunks × 20 elements
         let chunk_dims = vec![20u32];
-        let chunks =
+        let chunks = read_extensible_array_chunks(
-            read_extensible_array_chunks(&file_data, &header, &ds_dims, &chunk_dims, 8, os, ls)
+            &file_data,
+            &header,
+            &ds_dims,
+            None,
+            &chunk_dims,
+            8,
+            os,
+            ls,
+        )
         .unwrap();
 
         assert_eq!(chunks.len(), 2);
@@ -1023,8 +985,16 @@ mod tests {
         let header = ExtensibleArrayHeader::parse(&file_data, aehd_offset, os, ls).unwrap();
         let ds_dims = vec![40u64];
         let chunk_dims = vec![10u32];
-        let chunks =
+        let chunks = read_extensible_array_chunks(
-            read_extensible_array_chunks(&file_data, &header, &ds_dims, &chunk_dims, 8, os, ls)
+            &file_data,
+            &header,
+            &ds_dims,
+            None,
+            &chunk_dims,
+            8,
+            os,
+            ls,
+        )
         .unwrap();
 
         assert_eq!(chunks.len(), 4);
@@ -1047,10 +1017,8 @@ mod tests {
     #[test]
     fn read_element_unallocated() {
         let data = vec![0xFFu8; 16];
-        let num_chunks = vec![5u64];
+        let grid = ChunkGrid::fixed_array(&[50], None, &[10]).unwrap();
-        let chunk_dims = vec![10u32];
+        let (info, consumed) = read_element(&data, 0, 0, 8, 8, 80, 0, &grid).unwrap();
-        let (info, consumed) =
-            read_element(&data, 0, 0, 8, 8, 80, 0, &num_chunks, &chunk_dims).unwrap();
         assert!(info.is_none());
         assert_eq!(consumed, 8);
     }
@@ -1069,20 +1037,9 @@ mod tests {
         // Filter mask
         data[12..16].copy_from_slice(&0u32.to_le_bytes());
 
-        let num_chunks = vec![5u64];
+        let grid = ChunkGrid::fixed_array(&[50], None, &[10]).unwrap();
-        let chunk_dims = vec![10u32];
+        let (info, consumed) =
-        let (info, consumed) = read_element(
+            read_element(&data, 0, 1, elem_size as u8, os, 80, 2, &grid).unwrap();
-            &data,
-            0,
-            1,
-            elem_size as u8,
-            os,
-            80,
-            2,
-            &num_chunks,
-            &chunk_dims,
-        )
-        .unwrap();
         let ci = info.unwrap();
         assert_eq!(ci.address, 0x2000);
         assert_eq!(ci.chunk_size, 120);

crates/clawhdf5-format/src/file_writer.rs

@@ -1124,7 +1124,12 @@ impl FileWriter {
         let is_chunked: Vec<bool> = all_ds
             .iter()
             .enumerate()
-            .map(|(i, d)| !is_vds[i] && (d.chunk_options.is_chunked() || d.maxshape.is_some()))
+            .map(|(i, d)| {
+                // Only a dataset that can grow needs chunks; a maxshape equal
+                // to the shape is as fixed as no maxshape at all.
+                let resizable = d.maxshape.as_ref().is_some_and(|m| *m != d.ds.dimensions);
+                !is_vds[i] && (d.chunk_options.is_chunked() || resizable)
+            })
             .collect();
         // Determine which datasets use compact storage
         let is_compact: Vec<bool> = all_ds
@@ -1244,7 +1249,7 @@ impl FileWriter {
                     &pre,
                     dummy_cursor,
                     d.maxshape.as_deref(),
-                );
+                )?;
                 dummy_cursor += result.data_bytes.len() as u64;
                 let dense_blob = if ds_dense[i] {
                     Some(build_dense_attrs(&d.attrs, 0))
@@ -1424,7 +1429,7 @@ impl FileWriter {
                         .expect("chunked dataset missing precompressed cache"),
                     base_address,
                     d.maxshape.as_deref(),
-                );
+                )?;
                 cursor2 += result.data_bytes.len();
                 let oh = build_chunked_dataset_oh(
                     &d.dt,

crates/clawhdf5-format/src/fixed_array.rs

@@ -6,6 +6,7 @@ extern crate alloc;
 #[cfg(not(feature = "std"))]
 use alloc::{format, vec, vec::Vec};
 
+use crate::chunk_grid::ChunkGrid;
 use crate::chunked_read::ChunkInfo;
 use crate::error::FormatError;
 
@@ -151,13 +152,13 @@ pub fn read_fixed_array_chunks(
     file_data: &[u8],
     header: &FixedArrayHeader,
     dataset_dims: &[u64],
+    max_dims: Option<&[u64]>,
     chunk_dimensions: &[u32],
     element_size: u32,
     offset_size: u8,
     _length_size: u8,
 ) -> Result<Vec<ChunkInfo>, FormatError> {
     let db_offset = header.data_block_address as usize;
-    let rank = chunk_dimensions.len();
 
     // Parse data block header: FADB(4) + version(1) + client_id(1) + header_address(offset_size)
     let db_header_size = 4 + 1 + 1 + offset_size as usize;
@@ -198,19 +199,10 @@ pub fn read_fixed_array_chunks(
             ))
     };
 
-    // Compute chunk offsets based on index.
+    // The index is laid out over the chunk grid of the *maximum* dimensions
-    // Chunks are stored in row-major order within the dataset space.
+    // (row-major), so a dataset smaller than its maxshape has gaps.
-    let mut num_chunks_per_dim = Vec::with_capacity(rank);
+    let dims_u64: Vec<u64> = chunk_dimensions.iter().map(|&d| d as u64).collect();
-    for d_idx in 0..rank {
+    let grid = ChunkGrid::fixed_array(dataset_dims, max_dims, &dims_u64)?;
-        let ch_dim = chunk_dimensions[d_idx] as u64;
-        if ch_dim == 0 {
-            return Err(FormatError::ChunkedReadError(
-                "chunk dimension is zero".into(),
-            ));
-        }
-        let ds_dim = dataset_dims[d_idx];
-        num_chunks_per_dim.push(ds_dim.div_ceil(ch_dim));
-    }
 
     let chunk_byte_size: u64 =
         chunk_dimensions.iter().map(|&d| d as u64).product::<u64>() * element_size as u64;
@@ -226,7 +218,11 @@ pub fn read_fixed_array_chunks(
                 header.element_size,
                 chunk_byte_size,
             )? {
-                let offsets = index_to_chunk_offsets(i, &num_chunks_per_dim, chunk_dimensions);
+                // A slot beyond the current extent is ignored, as the
+                // library does.
+                let Some(offsets) = grid.offsets(i as u64) else {
+                    return Ok(());
+                };
                 chunks.push(ChunkInfo {
                     chunk_size,
                     filter_mask,
@@ -367,27 +363,6 @@ fn parse_fa_element(
     }
 }
 
-/// Convert a linear chunk index to N-dimensional chunk offsets in dataset space.
-fn index_to_chunk_offsets(
-    index: usize,
-    num_chunks_per_dim: &[u64],
-    chunk_dimensions: &[u32],
-) -> Vec<u64> {
-    let rank = num_chunks_per_dim.len();
-    let mut offsets = vec![0u64; rank];
-    let mut remaining = index as u64;
-    for d in (0..rank).rev() {
-        let nchunks = num_chunks_per_dim[d];
-        if nchunks == 0 {
-            continue;
-        }
-        let chunk_idx = remaining % nchunks;
-        remaining /= nchunks;
-        offsets[d] = chunk_idx * chunk_dimensions[d] as u64;
-    }
-    offsets
-}
-
 /// Read a variable-length little-endian unsigned integer.
 fn read_variable_length(data: &[u8], size: usize) -> Result<u64, FormatError> {
     if size > 8 || data.len() < size {
@@ -416,44 +391,21 @@ mod tests {
 
     #[test]
     fn index_to_offsets_1d() {
-        let num_chunks = vec![5u64];
+        let g = ChunkGrid::fixed_array(&[100], None, &[20]).unwrap();
-        let chunk_dims = vec![20u32];
+        assert_eq!(g.offsets(0).unwrap(), vec![0]);
-        assert_eq!(index_to_chunk_offsets(0, &num_chunks, &chunk_dims), vec![0]);
+        assert_eq!(g.offsets(1).unwrap(), vec![20]);
-        assert_eq!(
+        assert_eq!(g.offsets(4).unwrap(), vec![80]);
-            index_to_chunk_offsets(1, &num_chunks, &chunk_dims),
-            vec![20]
-        );
-        assert_eq!(
-            index_to_chunk_offsets(4, &num_chunks, &chunk_dims),
-            vec![80]
-        );
     }
 
     #[test]
     fn index_to_offsets_2d() {
         // 10x6 dataset with 4x3 chunks => ceil(10/4)=3, ceil(6/3)=2 => 6 chunks
-        let num_chunks = vec![3u64, 2];
+        let g = ChunkGrid::fixed_array(&[10, 6], None, &[4, 3]).unwrap();
-        let chunk_dims = vec![4u32, 3];
+        assert_eq!(g.offsets(0).unwrap(), vec![0, 0]);
-        assert_eq!(
+        assert_eq!(g.offsets(1).unwrap(), vec![0, 3]);
-            index_to_chunk_offsets(0, &num_chunks, &chunk_dims),
+        assert_eq!(g.offsets(2).unwrap(), vec![4, 0]);
-            vec![0, 0]
+        assert_eq!(g.offsets(3).unwrap(), vec![4, 3]);
-        );
+        assert_eq!(g.offsets(5).unwrap(), vec![8, 3]);
-        assert_eq!(
-            index_to_chunk_offsets(1, &num_chunks, &chunk_dims),
-            vec![0, 3]
-        );
-        assert_eq!(
-            index_to_chunk_offsets(2, &num_chunks, &chunk_dims),
-            vec![4, 0]
-        );
-        assert_eq!(
-            index_to_chunk_offsets(3, &num_chunks, &chunk_dims),
-            vec![4, 3]
-        );
-        assert_eq!(
-            index_to_chunk_offsets(5, &num_chunks, &chunk_dims),
-            vec![8, 3]
-        );
     }
 
     #[test]
@@ -517,7 +469,7 @@ mod tests {
 
         let read = |f: &[u8], fahd: usize| -> Result<Vec<ChunkInfo>, FormatError> {
             let h = FixedArrayHeader::parse(f, fahd, 8, 8)?;
-            read_fixed_array_chunks(f, &h, &[60], &[20], 8, 8, 8)
+            read_fixed_array_chunks(f, &h, &[60], None, &[20], 8, 8, 8)
         };
 
         let (clean, fahd) = build();
@@ -562,7 +514,7 @@ mod tests {
         let db = 0x100usize;
         buf[db..db + 4].copy_from_slice(b"FADB");
         let header = FixedArrayHeader::parse(&buf, fahd, 8, 8).unwrap();
-        let r = read_fixed_array_chunks(&buf, &header, &[100], &[20], 8, 8, 8);
+        let r = read_fixed_array_chunks(&buf, &header, &[100], None, &[20], 8, 8, 8);
         assert!(r.is_err());
     }
 
@@ -579,7 +531,7 @@ mod tests {
         stamp_checksum(&mut buf, fahd, fahd + 24);
         buf[0x80..0x84].copy_from_slice(b"FADB");
         let header = FixedArrayHeader::parse(&buf, fahd, 8, 8).unwrap();
-        let r = read_fixed_array_chunks(&buf, &header, &[100], &[20], 8, 8, 8);
+        let r = read_fixed_array_chunks(&buf, &header, &[100], None, &[20], 8, 8, 8);
         assert!(r.is_err());
     }
 
@@ -602,7 +554,7 @@ mod tests {
             data_block_address: (usize::MAX - 4) as u64,
         };
         let buf = vec![0u8; 64];
-        let r = read_fixed_array_chunks(&buf, &header, &[100], &[20], 8, 8, 8);
+        let r = read_fixed_array_chunks(&buf, &header, &[100], None, &[20], 8, 8, 8);
         assert!(r.is_err());
     }
 
@@ -664,6 +616,7 @@ mod tests {
             &file_data,
             &header,
             &ds_dims,
+            None,
             &chunk_dims,
             8,
             offset_size,
@@ -740,6 +693,7 @@ mod tests {
             &file_data,
             &header,
             &ds_dims,
+            None,
             &chunk_dims,
             8,
             offset_size,
@@ -840,6 +794,7 @@ mod tests {
             &file_data,
             &header,
             &ds_dims,
+            None,
             &chunk_dims,
             8,
             offset_size,

crates/clawhdf5-format/src/lib.rs

@@ -54,6 +54,7 @@ pub mod btree_v1;
 pub mod btree_v2;
 pub mod checksum;
 pub mod chunk_cache;
+mod chunk_grid;
 pub mod chunk_index;
 pub mod chunked_read;
 pub mod chunked_write;

crates/clawhdf5/tests/chunk_index_interop.rs

@@ -0,0 +1,593 @@
+//! Fixed Array / Extensible Array chunk-index interop with libhdf5 (via h5py).
+//!
+//! Both indexes place each chunk at a linear index computed from the
+//! dataset's *maximum* dimensions, and the Extensible Array additionally
+//! moves its unlimited dimension to the slowest-varying position. Getting
+//! either wrong reads (or writes) every chunk after the first row in the
+//! wrong place, silently, so these tests compare every value.
+//!
+//! Skipped when python3 with h5py is unavailable, unless
+//! `CLAWHDF5_REQUIRE_INTEROP=1`.
+
+use std::process::Command;
+
+use clawhdf5::{File, FileBuilder};
+
+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"])
+        .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) -> String {
+    let output = Command::new(python())
+        .args(["-c", script])
+        .output()
+        .expect("failed to run python");
+    if !output.status.success() {
+        panic!(
+            "Python script failed:\nSTDOUT: {}\nSTDERR: {}",
+            String::from_utf8_lossy(&output.stdout),
+            String::from_utf8_lossy(&output.stderr)
+        );
+    }
+    String::from_utf8_lossy(&output.stdout).trim().to_string()
+}
+
+/// Row-major `arange` of `shape`, cropped to `crop` (the current extent).
+fn arange_cropped(full: &[usize], crop: &[usize]) -> Vec<i32> {
+    let n: usize = crop.iter().product();
+    let mut out = Vec::with_capacity(n);
+    for flat in 0..n {
+        let mut rem = flat;
+        let mut src = 0usize;
+        let mut stride = 1usize;
+        let mut coords = vec![0usize; crop.len()];
+        for d in (0..crop.len()).rev() {
+            coords[d] = rem % crop[d];
+            rem /= crop[d];
+        }
+        for d in (0..full.len()).rev() {
+            src += coords[d] * stride;
+            stride *= full[d];
+        }
+        out.push(src as i32);
+    }
+    out
+}
+
+/// One `i4` dataset, filled with `arange` over `full` and then resized to
+/// `shape` (equal to `full` unless the case shrinks it).
+struct Case {
+    name: &'static str,
+    full: Vec<usize>,
+    shape: Vec<usize>,
+    chunks: Vec<usize>,
+    maxshape: &'static str,
+    extra: &'static str,
+    index: &'static str,
+}
+
+fn py_tuple(v: &[usize]) -> String {
+    let parts: Vec<String> = v.iter().map(|x| x.to_string()).collect();
+    format!("({},)", parts.join(","))
+}
+
+/// Have h5py (`libver="latest"`, so Fixed/Extensible Array indexes) write
+/// every case to one file, then read each back and compare every value.
+fn check_h5py_written(cases: &[Case]) {
+    let dir = tempfile::tempdir().unwrap();
+    let path = dir.path().join("h5py_chunk_index.h5");
+    let path_str = path.display().to_string();
+
+    let mut script =
+        format!("import h5py, numpy as np\nf = h5py.File(r'{path_str}', 'w', libver='latest')\n");
+    for c in cases {
+        script += &format!(
+            "d = f.create_dataset('{name}', data=np.arange({n}, dtype='i4').reshape({full}), \
+             chunks={chunks}, maxshape={maxshape}{extra})\n\
+             d.resize({shape})\n",
+            name = c.name,
+            n = c.full.iter().product::<usize>(),
+            full = py_tuple(&c.full),
+            chunks = py_tuple(&c.chunks),
+            maxshape = c.maxshape,
+            extra = c.extra,
+            shape = py_tuple(&c.shape),
+        );
+    }
+    script += "f.close()\n";
+    run_python(&script);
+
+    let file = File::open(&path).unwrap();
+    for c in cases {
+        let ds = file.dataset(c.name).unwrap();
+        let shape: Vec<usize> = ds.shape().unwrap().iter().map(|&d| d as usize).collect();
+        assert_eq!(shape, c.shape, "{}: shape", c.name);
+        let got = ds.read_i32().unwrap();
+        let want = arange_cropped(&c.full, &c.shape);
+        let bad = got.iter().zip(&want).filter(|(a, b)| a != b).count();
+        assert_eq!(
+            got,
+            want,
+            "{}: {bad} of {} values differ (index {})",
+            c.name,
+            want.len(),
+            c.index
+        );
+    }
+}
+
+/// h5py-written Extensible Array whose unlimited dimension is not the first,
+/// with the current shape smaller than the finite maximum: the library
+/// swizzles the unlimited dimension to the slowest position and strides the
+/// rest by their maximum chunk counts.
+#[test]
+fn h5py_extensible_array_partial_extent_reads_correctly() {
+    skip_if_no_python!();
+    check_h5py_written(&[
+        // The `ea_fa_partial.h5` repro from the conformance sweep.
+        Case {
+            name: "ea_10_none",
+            full: vec![4, 6],
+            shape: vec![4, 6],
+            chunks: vec![2, 3],
+            maxshape: "(10, None)",
+            extra: "",
+            index: "EA, unlimited dim 1",
+        },
+        Case {
+            name: "ea_none_10",
+            full: vec![4, 6],
+            shape: vec![4, 6],
+            chunks: vec![2, 3],
+            maxshape: "(None, 10)",
+            extra: "",
+            index: "EA, unlimited dim 0",
+        },
+        Case {
+            name: "ea_3d_mid",
+            full: vec![3, 4, 5],
+            shape: vec![3, 4, 5],
+            chunks: vec![2, 3, 2],
+            maxshape: "(5, None, 7)",
+            extra: "",
+            index: "EA, unlimited dim 1 of 3",
+        },
+        Case {
+            name: "ea_3d_last_gzip",
+            full: vec![3, 4, 5],
+            shape: vec![3, 4, 5],
+            chunks: vec![2, 3, 2],
+            maxshape: "(5, 9, None)",
+            extra: ", compression='gzip'",
+            index: "EA, unlimited dim 2 of 3, filtered",
+        },
+        // Many chunks: crosses data blocks, super blocks and paging.
+        Case {
+            name: "ea_many",
+            full: vec![3, 1500],
+            shape: vec![3, 1500],
+            chunks: vec![1, 1],
+            maxshape: "(4, None)",
+            extra: "",
+            index: "EA, 4500 slots",
+        },
+        // Shrunk after writing: chunks beyond the extent must be ignored.
+        Case {
+            name: "ea_shrunk",
+            full: vec![8, 9],
+            shape: vec![3, 4],
+            chunks: vec![2, 3],
+            maxshape: "(10, None)",
+            extra: "",
+            index: "EA, shrunk",
+        },
+    ]);
+}
+
+/// h5py-written Fixed Array with the current shape smaller than a finite
+/// maxshape: the index has one slot per chunk of the *maximum* extent.
+#[test]
+fn h5py_fixed_array_partial_extent_reads_correctly() {
+    skip_if_no_python!();
+    check_h5py_written(&[
+        Case {
+            name: "fa_20_10",
+            full: vec![4, 6],
+            shape: vec![4, 6],
+            chunks: vec![2, 3],
+            maxshape: "(20, 10)",
+            extra: "",
+            index: "FA",
+        },
+        Case {
+            name: "fa_3d_gzip",
+            full: vec![3, 4, 5],
+            shape: vec![3, 4, 5],
+            chunks: vec![2, 3, 2],
+            maxshape: "(6, 8, 10)",
+            extra: ", compression='gzip'",
+            index: "FA, filtered",
+        },
+        // Paged (> 1024 slots) with most of them beyond the extent.
+        Case {
+            name: "fa_paged",
+            full: vec![30, 50],
+            shape: vec![30, 50],
+            chunks: vec![1, 1],
+            maxshape: "(40, 60)",
+            extra: "",
+            index: "FA, 2400 slots, paged",
+        },
+        Case {
+            name: "fa_shrunk",
+            full: vec![8, 9],
+            shape: vec![5, 2],
+            chunks: vec![2, 3],
+            maxshape: "(20, 10)",
+            extra: "",
+            index: "FA, shrunk",
+        },
+    ]);
+}
+
+// ===========================================================================
+// Files we write, read back by libhdf5 (h5py and h5dump) and by us
+// ===========================================================================
+
+/// One `i4` dataset we write, filled with `arange` over `shape`.
+struct WriteCase {
+    name: String,
+    shape: Vec<u64>,
+    chunks: Vec<u64>,
+    maxshape: Option<Vec<u64>>,
+    deflate: bool,
+}
+
+fn wcase(name: &str, shape: &[u64], chunks: &[u64], maxshape: Option<&[u64]>) -> WriteCase {
+    WriteCase {
+        name: name.to_string(),
+        shape: shape.to_vec(),
+        chunks: chunks.to_vec(),
+        maxshape: maxshape.map(<[u64]>::to_vec),
+        deflate: false,
+    }
+}
+
+fn h5dump_available() -> bool {
+    Command::new("h5dump")
+        .arg("--version")
+        .output()
+        .map(|o| o.status.success())
+        .unwrap_or(false)
+}
+
+/// Write every case into one file with our writer, then check that our own
+/// reader, h5py and h5dump (when installed) all return every value. Only the
+/// libhdf5 half is skipped without h5py.
+fn check_we_write(cases: &[WriteCase]) {
+    let dir = tempfile::tempdir().unwrap();
+    let path = dir.path().join("ours_chunk_index.h5");
+    let path_str = path.display().to_string();
+
+    let mut b = FileBuilder::new();
+    for c in cases {
+        let n: u64 = c.shape.iter().product();
+        let data: Vec<i32> = (0..n as i32).collect();
+        let ds = b.create_dataset(&c.name);
+        ds.with_i32_data(&data)
+            .with_shape(&c.shape)
+            .with_chunks(&c.chunks);
+        if let Some(ms) = &c.maxshape {
+            ds.with_maxshape(ms);
+        }
+        if c.deflate {
+            ds.with_deflate(4);
+        }
+    }
+    b.write(&path).unwrap();
+
+    // Our reader.
+    let file = File::open(&path).unwrap();
+    for c in cases {
+        let got = file.dataset(&c.name).unwrap().read_i32().unwrap();
+        let n: u64 = c.shape.iter().product();
+        let bad = got
+            .iter()
+            .enumerate()
+            .filter(|&(i, &v)| v != i as i32)
+            .count();
+        assert!(
+            got.len() == n as usize && bad == 0,
+            "{}: our reader: {bad} of {n} values wrong",
+            c.name
+        );
+    }
+
+    // libhdf5 via h5py.
+    skip_if_no_python!();
+    let mut script =
+        format!("import h5py, numpy as np\nbad = []\nf = h5py.File(r'{path_str}', 'r')\n");
+    for c in cases {
+        let shape: Vec<String> = c.shape.iter().map(u64::to_string).collect();
+        let maxshape: Vec<String> = c
+            .maxshape
+            .as_ref()
+            .unwrap_or(&c.shape)
+            .iter()
+            .map(|&d| {
+                if d == u64::MAX {
+                    "None".to_string()
+                } else {
+                    d.to_string()
+                }
+            })
+            .collect();
+        script += &format!(
+            "d = f['{name}']\n\
+             want = np.arange({n}, dtype='i4').reshape(({shape},))\n\
+             got = d[()]\n\
+             if d.maxshape != ({maxshape},): bad.append(('{name}', 'maxshape', d.maxshape))\n\
+             elif not np.array_equal(got, want): \
+             bad.append(('{name}', int((got != want).sum()), 'of', got.size))\n",
+            name = c.name,
+            n = c.shape.iter().product::<u64>(),
+            shape = shape.join(","),
+            maxshape = maxshape.join(","),
+        );
+    }
+    script += "print(bad if bad else 'OK')\n";
+    let out = run_python(&script);
+    assert_eq!(out, "OK", "h5py disagrees");
+
+    // libhdf5's own tool, when installed.
+    if h5dump_available() {
+        let o = Command::new("h5dump").arg(&path).output().unwrap();
+        let stderr = String::from_utf8_lossy(&o.stderr);
+        assert!(
+            o.status.success() && !stderr.to_lowercase().contains("error"),
+            "h5dump failed: {stderr}"
+        );
+    }
+
+    // Let libhdf5 grow every resizable dataset by two chunks per dimension
+    // (capped at the maxshape) and rewrite it, which updates our index in
+    // place and inserts new chunks into it. Then both readers must agree.
+    let script = format!(
+        r#"
+import h5py, numpy as np
+grown = {{}}
+with h5py.File(r'{path_str}', 'r+') as f:
+    for name in f:
+        d = f[name]
+        if d.chunks is None:
+            continue
+        new = tuple(s + 2 * c if m is None else min(m, s + 2 * c)
+                    for s, m, c in zip(d.shape, d.maxshape, d.chunks))
+        if new == d.shape:
+            continue
+        old = d[()]
+        full = np.full(new, -7, 'i4')
+        full[tuple(slice(0, s) for s in old.shape)] = old
+        d.resize(new)
+        d[...] = full
+        grown[name] = (list(old.shape), list(new))
+with h5py.File(r'{path_str}', 'r') as f:
+    for name, (old, new) in grown.items():
+        want = np.full(new, -7, 'i4')
+        want[tuple(slice(0, s) for s in old)] = np.arange(int(np.prod(old)), dtype='i4').reshape(old)
+        assert np.array_equal(f[name][()], want), name
+for name, (old, new) in grown.items():
+    print(name, ','.join(map(str, old)), ','.join(map(str, new)))
+"#
+    );
+    let out = run_python(&script);
+    let growable = cases
+        .iter()
+        .filter(|c| c.maxshape.as_ref().is_some_and(|m| *m != c.shape))
+        .count();
+    assert_eq!(out.lines().count(), growable, "libhdf5 grew: {out}");
+    let dims = |s: &str| -> Vec<usize> { s.split(',').map(|x| x.parse().unwrap()).collect() };
+    let file = File::open(&path).unwrap();
+    for line in out.lines() {
+        let mut parts = line.split(' ');
+        let (name, old, new) = (
+            parts.next().unwrap(),
+            dims(parts.next().unwrap()),
+            dims(parts.next().unwrap()),
+        );
+        let got = file.dataset(name).unwrap().read_i32().unwrap();
+        let n: usize = new.iter().product();
+        let mut want = vec![-7i32; n];
+        for (flat, w) in want.iter_mut().enumerate() {
+            let mut rem = flat;
+            let mut coords = vec![0usize; new.len()];
+            for d in (0..new.len()).rev() {
+                coords[d] = rem % new[d];
+                rem /= new[d];
+            }
+            if coords.iter().zip(&old).all(|(c, o)| c < o) {
+                *w = coords.iter().zip(&old).fold(0, |acc, (c, o)| acc * o + c) as i32;
+            }
+        }
+        let bad = got.iter().zip(&want).filter(|(a, b)| a != b).count();
+        assert!(
+            got.len() == n && bad == 0,
+            "{name}: after libhdf5 grew it, our reader got {bad} of {n} values wrong"
+        );
+    }
+}
+
+/// A Fixed Array with more than 1024 elements must be paged, or libhdf5
+/// rejects the data block's checksum.
+#[test]
+fn we_write_paged_fixed_array() {
+    let mut cases: Vec<WriteCase> = [1023u64, 1024, 1025, 2048, 5000]
+        .iter()
+        .map(|&n| wcase(&format!("fa_{n}"), &[n * 4], &[4], None))
+        .collect();
+    // Filtered elements are wider; a 2-D grid pages the same way.
+    let mut filtered = wcase("fa_1500_deflate", &[1500 * 4], &[4], None);
+    filtered.deflate = true;
+    cases.push(filtered);
+    cases.push(wcase("fa_2d_1100", &[110, 40], &[1, 4], None));
+    check_we_write(&cases);
+}
+
+/// An Extensible Array holds 4 elements in its index block and 240 in the
+/// data blocks the index block addresses; everything after that lives under
+/// super blocks, and from ~131K elements on in paged data blocks. Chunks past
+/// index 243 used to be written but never indexed (read back as fill by us
+/// and by libhdf5).
+#[test]
+fn we_write_extensible_array_past_index_block() {
+    let unl: &[u64] = &[u64::MAX];
+    let mut cases: Vec<WriteCase> = [1u64, 4, 5, 243, 244, 245, 300, 1000, 5000]
+        .iter()
+        .map(|&n| wcase(&format!("ea_{n}"), &[n * 4], &[4], Some(unl)))
+        .collect();
+    let mut filtered = wcase("ea_300_deflate", &[300 * 4], &[4], Some(unl));
+    filtered.deflate = true;
+    cases.push(filtered);
+    // Several super blocks and paged data blocks (level 13, the first with
+    // data blocks over 1024 elements, starts at element 4 + 131056).
+    cases.push(wcase("ea_140000", &[140_000], &[1], Some(unl)));
+    check_we_write(&cases);
+}
+
+/// A maxshape larger than the shape: the index must be laid out over the
+/// chunks of the maximum extent (libhdf5 read our Fixed Array past its end:
+/// "addr overflow"), and an Extensible Array whose unlimited dimension is not
+/// the first must swizzle it to the slowest position (libhdf5 read our
+/// `(20, None)` dataset scrambled).
+#[test]
+fn we_write_maxshape_larger_than_shape() {
+    const U: u64 = u64::MAX;
+    let mut cases = vec![
+        // Fixed Array over the maximum extent.
+        wcase("fa2d_finite_max", &[20, 30], &[5, 5], Some(&[40, 60])),
+        wcase("fa1d_finite_max", &[40], &[4], Some(&[100])),
+        wcase("fa3d_edges", &[6, 7, 8], &[4, 3, 5], Some(&[10, 9, 20])),
+        wcase("fa_paged_max", &[30, 50], &[1, 1], Some(&[40, 60])),
+        wcase("fa_one_chunk_now", &[5], &[5], Some(&[50])),
+        // Extensible Array, unlimited dimension first (no swizzle) ...
+        wcase("ea2d_unl_fin", &[20, 30], &[5, 5], Some(&[U, 30])),
+        wcase("ea2d_unl_fin_max", &[20, 30], &[5, 5], Some(&[U, 60])),
+        // ... and not first (swizzled).
+        wcase("ea2d_fin_unl", &[20, 30], &[5, 5], Some(&[20, U])),
+        wcase("ea2d_fin_max_unl", &[20, 30], &[5, 5], Some(&[40, U])),
+        wcase("ea3d_mid", &[6, 7, 8], &[4, 3, 5], Some(&[10, U, 20])),
+        // Past the index block and into super blocks, swizzled.
+        wcase("ea2d_many", &[3, 2000], &[1, 1], Some(&[4, U])),
+    ];
+    let mut filtered = wcase(
+        "ea3d_last_deflate",
+        &[6, 7, 8],
+        &[4, 3, 5],
+        Some(&[6, 8, U]),
+    );
+    filtered.deflate = true;
+    cases.push(filtered);
+    check_we_write(&cases);
+}
+
+/// More than one unlimited dimension needs a version-2 B-tree chunk index,
+/// as the library uses; an Extensible Array for `(None, None)` made libhdf5
+/// refuse the whole file ("already found unlimited dimension").
+#[test]
+fn we_write_btree_v2_for_several_unlimited_dims() {
+    const U: u64 = u64::MAX;
+    let mut cases = vec![
+        wcase("unl_unl", &[20, 30], &[5, 5], Some(&[U, U])),
+        wcase("unl_fin_unl", &[6, 7, 8], &[4, 3, 5], Some(&[U, 9, U])),
+        // More records than the library's 2048-byte node holds (84 here).
+        wcase("unl_unl_2400", &[40, 60], &[1, 1], Some(&[U, U])),
+        wcase("unl_unl_empty", &[0, 0], &[4, 4], Some(&[U, U])),
+    ];
+    let mut filtered = wcase("unl_unl_deflate", &[6, 7, 8], &[4, 3, 5], Some(&[U, U, U]));
+    filtered.deflate = true;
+    cases.push(filtered);
+    check_we_write(&cases);
+}
+
+/// A single-leaf B-tree has a 16-bit record count; beyond it the writer
+/// refuses rather than writing a tree libhdf5 would misread.
+#[test]
+fn btree_v2_index_past_one_leaf_is_refused() {
+    let mut b = FileBuilder::new();
+    b.create_dataset("d")
+        .with_i32_data(&vec![0i32; 70_000])
+        .with_shape(&[70_000, 1])
+        .with_chunks(&[1, 1])
+        .with_maxshape(&[u64::MAX, u64::MAX]);
+    let dir = tempfile::tempdir().unwrap();
+    assert!(b.write(dir.path().join("too_many.h5")).is_err());
+}
+
+/// A maxshape equal to the shape cannot grow, so it needs no chunks: the
+/// dataset stays contiguous (as h5py makes it) unless chunks are requested.
+#[test]
+fn maxshape_equal_to_shape_stays_contiguous() {
+    let dir = tempfile::tempdir().unwrap();
+    let path = dir.path().join("ms_eq.h5");
+    let data: Vec<i32> = (0..40).collect();
+    let mut b = FileBuilder::new();
+    b.create_dataset("plain")
+        .with_i32_data(&data)
+        .with_shape(&[40])
+        .with_maxshape(&[40]);
+    b.create_dataset("chunked")
+        .with_i32_data(&data)
+        .with_shape(&[40])
+        .with_maxshape(&[40])
+        .with_chunks(&[8]);
+    b.write(&path).unwrap();
+
+    let file = File::open(&path).unwrap();
+    let plain = file.dataset("plain").unwrap();
+    assert_eq!(plain.read_i32().unwrap(), data);
+    assert_eq!(plain.max_dimensions().unwrap(), Some(vec![40]));
+    assert!(
+        plain.read_raw_ref().unwrap().is_some(),
+        "maxshape == shape should be contiguous"
+    );
+    let chunked = file.dataset("chunked").unwrap();
+    assert_eq!(chunked.read_i32().unwrap(), data);
+    assert!(chunked.read_raw_ref().unwrap().is_none());
+
+    skip_if_no_python!();
+    let out = run_python(&format!(
+        "import h5py, numpy as np\n\
+         f = h5py.File(r'{}', 'r')\n\
+         for n in ('plain', 'chunked'):\n\
+         \x20   d = f[n]\n\
+         \x20   assert np.array_equal(d[()], np.arange(40, dtype='i4')), n\n\
+         \x20   print(n, d.chunks, d.maxshape)\n",
+        path.display()
+    ));
+    assert_eq!(out, "plain None (40,)\nchunked (8,) (40,)");
+}