crates/clawhdf5-format/src/file_writer.rs

@@ -307,7 +307,7 @@ pub(crate) fn build_single_block_fractal_heap(
     base_address: u64,
     max_heap_size: u16,
     heap_id_length: u16,
-) -> FractalHeapBlock {
+) -> Result<FractalHeapBlock, FormatError> {
     let os = OFFSET_SIZE as usize;
     let ls = LENGTH_SIZE as usize;
     let block_offset_bytes = (max_heap_size as usize).div_ceil(8);
@@ -435,183 +435,346 @@ pub(crate) fn build_single_block_fractal_heap(
     blob.extend_from_slice(&frhp);
     blob.extend_from_slice(&dblock);
 
-    FractalHeapBlock {
+    Ok(FractalHeapBlock {
         blob,
         frhp_addr,
         btree_addr,
         heap_ids,
         heap_id_length,
-    }
+    })
 }
 
-/// Build a multi-block fractal heap: a root indirect block (FHIB) over multiple
+/// Build a multi-block fractal heap: a root indirect block (FHIB) over direct
-/// direct blocks sized by the doubling table. Used when the objects don't fit
+/// blocks sized by the doubling table. Used when the objects don't fit in a
-/// in a single direct block. Objects do not span blocks (no huge-object path).
+/// single direct block.
+///
+/// Rows of the doubling table whose block size exceeds the maximum direct
+/// block size hold child indirect blocks, as the HDF5 spec (and libhdf5)
+/// reads them: a child in row `r` spans that row's block size of heap space
+/// and has `log2(size) - log2(start * width) + 1` rows of its own, which may
+/// in turn hold indirect blocks. Objects are packed into direct blocks in
+/// heap-offset order and never span blocks; a block too small for the next
+/// object is left unallocated (an undefined address), as libhdf5 skips rows
+/// when it needs a bigger block. There is no huge-object path.
 fn build_multiblock_fractal_heap(
     serialized: &[Vec<u8>],
     base_address: u64,
     max_heap_size: u16,
     heap_id_length: u16,
-) -> FractalHeapBlock {
+) -> Result<FractalHeapBlock, FormatError> {
     let os = OFFSET_SIZE as usize;
     let block_offset_bytes = (max_heap_size as usize).div_ceil(8);
-    let max_direct_block_size: u64 = 65536;
+    let geom = HeapGeometry {
-    let table_width: u16 = 4;
+        width: 4,
-    let starting_block_size: u64 = 512;
+        starting_block_size: 512,
-    let dblock_header_size = 4 + 1 + os + block_offset_bytes + 4;
+        max_direct_block_size: 65536,
-    let block_capacity =
+        dblock_header_size: 4 + 1 + os + block_offset_bytes + 4,
-        |row: usize| block_size_for_row(starting_block_size, row) - dblock_header_size as u64;
+        iblock_fixed_size: 5 + os + block_offset_bytes + 4,
+        max_heap_size,
+    };
 
-    // ---- Pack objects into direct blocks (row-major over the doubling table) ----
+    // ---- Pack objects into the doubling table ----
-    struct Blk {
+    let mut packer = HeapPacker {
-        row: usize,
+        geom: &geom,
-        size: u64,
+        objects: serialized,
-        heap_offset: u64,
+        next: 0,
-        data: Vec<u8>,
+        blocks: Vec::new(),
-    }
+        obj_loc: vec![(0, 0); serialized.len()],
-    let mut blocks: Vec<Blk> = Vec::new();
+    };
-    // Each object's (heap_offset, length) for the heap ID.
+    let root = packer.fill(0, None)?;
-    let mut obj_loc: Vec<(u64, u64)> = vec![(0, 0); serialized.len()];
+    let HeapPacker {
-
+        blocks, obj_loc, ..
-    let mut row = 0usize;
+    } = packer;
-    let mut col = 0u16;
-    let mut heap_off = 0u64;
-    let mut cur: Option<Blk> = None;
-
-    for (idx, s) in serialized.iter().enumerate() {
-        loop {
-            if cur.is_none() {
-                let size = block_size_for_row(starting_block_size, row);
-                cur = Some(Blk {
-                    row,
-                    size,
-                    heap_offset: heap_off,
-                    data: Vec::new(),
-                });
-            }
-            let blk = cur.as_mut().unwrap();
-            let cap = block_capacity(blk.row) as usize;
-            if !blk.data.is_empty() && blk.data.len() + s.len() > cap {
-                // Doesn't fit; finalize this block and advance to the next slot.
-                let finished = cur.take().unwrap();
-                heap_off += finished.size;
-                blocks.push(finished);
-                col += 1;
-                if col >= table_width {
-                    col = 0;
-                    row += 1;
-                }
-                continue;
-            }
-            // Place the object (a fresh block always accepts at least one object
-            // up to its capacity; objects larger than a max block are unsupported).
-            let pos_in_block = dblock_header_size + blk.data.len();
-            obj_loc[idx] = (blk.heap_offset + pos_in_block as u64, s.len() as u64);
-            blk.data.extend_from_slice(s);
-            break;
-        }
-    }
-    if let Some(b) = cur.take() {
-        blocks.push(b);
-    }
-
-    let cur_rows = (blocks.last().map(|b| b.row).unwrap_or(0) + 1) as u16;
 
     // ---- Addresses ----
     let frhp_size = frhp_header_size(os, LENGTH_SIZE as usize);
     let frhp_addr = base_address;
     let fhib_addr = frhp_addr + frhp_size as u64;
-    let fhib_entries = cur_rows as usize * table_width as usize;
+    let heap_len = root.subtree_size(&geom, &blocks);
-    let fhib_size = 5 + os + block_offset_bytes + fhib_entries * os + 4;
+    let btree_addr = fhib_addr + heap_len;
-    let first_dblock_addr = fhib_addr + fhib_size as u64;
 
-    // Assign each used block an address (laid out consecutively after the FHIB).
+    // Bookkeeping totals, as libhdf5 keeps them: the managed space is what
-    let mut blk_addrs: Vec<u64> = Vec::with_capacity(blocks.len());
+    // the root's rows span, the allocated space the direct blocks written,
-    let mut a = first_dblock_addr;
+    // and the allocation iterator the heap offset after the last of them.
-    for b in &blocks {
+    let cur_rows = root.nrows as u16;
-        blk_addrs.push(a);
+    let managed_space: u64 = (0..root.nrows).map(|r| geom.row_size(r) * geom.width).sum();
-        a += b.size;
-    }
-    let heap_end = a;
-    let btree_addr = heap_end;
-
-    // Bookkeeping totals.
-    let managed_space: u64 = (0..cur_rows as usize)
-        .map(|r| block_size_for_row(starting_block_size, r) * table_width as u64)
-        .sum();
     let alloc_space: u64 = blocks.iter().map(|b| b.size).sum();
     let used: u64 = blocks
         .iter()
-        .map(|b| dblock_header_size as u64 + b.data.len() as u64)
+        .map(|b| geom.dblock_header_size as u64 + b.data.len() as u64)
         .sum();
     let free_space = alloc_space.saturating_sub(used);
+    let alloc_iter = blocks.last().map_or(0, |b| b.heap_offset + b.size);
 
     // ---- FRHP header ----
-    let max_managed = max_direct_block_size as u32 - dblock_header_size as u32;
+    let max_managed = geom.max_managed();
     let frhp = write_frhp(WriteFrhp {
         heap_id_length,
         max_managed,
         free_space,
         managed_space,
         alloc_space,
+        alloc_iter,
         nobjects: serialized.len() as u64,
-        table_width,
+        table_width: geom.width as u16,
-        starting_block_size,
+        starting_block_size: geom.starting_block_size,
-        max_direct_block_size,
+        max_direct_block_size: geom.max_direct_block_size,
         max_heap_size,
         root_addr: fhib_addr,
         cur_rows,
     });
     debug_assert_eq!(frhp.len(), frhp_size);
 
-    // ---- Root indirect block (FHIB) ----
+    // ---- Indirect and direct blocks, depth first after the root ----
-    let mut fhib = Vec::with_capacity(fhib_size);
-    fhib.extend_from_slice(b"FHIB");
-    fhib.push(0); // version
-    write_offset(&mut fhib, frhp_addr, OFFSET_SIZE);
-    fhib.extend_from_slice(&vec![0u8; block_offset_bytes]); // block offset = 0 (root)
-    for &addr in &blk_addrs {
-        write_offset(&mut fhib, addr, OFFSET_SIZE);
-    }
-    // Remaining slots within the current rows are unallocated.
-    for _ in blk_addrs.len()..fhib_entries {
-        write_undef_offset(&mut fhib, OFFSET_SIZE);
-    }
-    let fhib_checksum = crate::checksum::jenkins_lookup3(&fhib);
-    fhib.extend_from_slice(&fhib_checksum.to_le_bytes());
-    debug_assert_eq!(fhib.len(), fhib_size);
-
-    // ---- Direct blocks ----
     let mut blob = frhp;
-    blob.extend_from_slice(&fhib);
+    root.emit(&geom, &blocks, frhp_addr, fhib_addr, &mut blob);
-    for b in &blocks {
+    debug_assert_eq!(blob.len() as u64, frhp_size as u64 + heap_len);
-        let mut dblock = Vec::with_capacity(b.size as usize);
-        dblock.extend_from_slice(b"FHDB");
-        dblock.push(0); // version
-        write_offset(&mut dblock, frhp_addr, OFFSET_SIZE);
-        let mut bo = b.heap_offset.to_le_bytes().to_vec();
-        bo.truncate(block_offset_bytes);
-        dblock.extend_from_slice(&bo);
-        let cksum_pos = dblock.len();
-        dblock.extend_from_slice(&[0u8; 4]); // checksum placeholder
-        dblock.extend_from_slice(&b.data);
-        dblock.resize(b.size as usize, 0);
-        let cksum = crate::checksum::jenkins_lookup3(&dblock);
-        dblock[cksum_pos..cksum_pos + 4].copy_from_slice(&cksum.to_le_bytes());
-        blob.extend_from_slice(&dblock);
-    }
 
     let heap_ids: Vec<Vec<u8>> = obj_loc
         .iter()
         .map(|(off, len)| encode_managed_id(*off, *len, max_heap_size, heap_id_length))
         .collect();
 
-    FractalHeapBlock {
+    Ok(FractalHeapBlock {
         blob,
         frhp_addr,
         btree_addr,
         heap_ids,
         heap_id_length,
+    })
+}
+
+/// The doubling table of a heap the writer builds.
+struct HeapGeometry {
+    width: u64,
+    starting_block_size: u64,
+    max_direct_block_size: u64,
+    dblock_header_size: usize,
+    /// An indirect block's size without its child entries.
+    iblock_fixed_size: usize,
+    max_heap_size: u16,
+}
+
+impl HeapGeometry {
+    fn row_size(&self, row: usize) -> u64 {
+        block_size_for_row(self.starting_block_size, row)
+    }
+
+    /// Rows holding direct blocks: `log2(max_direct / start) + 2`.
+    fn max_direct_rows(&self) -> usize {
+        (self.max_direct_block_size / self.starting_block_size).ilog2() as usize + 2
+    }
+
+    /// `log2(start * width)`, libhdf5's `first_row_bits`.
+    fn first_row_bits(&self) -> u32 {
+        (self.starting_block_size * self.width).ilog2()
+    }
+
+    /// Rows of an indirect block spanning `size` bytes of heap space
+    /// (libhdf5's `H5HF__dtable_size_to_rows`).
+    fn rows_for_size(&self, size: u64) -> usize {
+        (size.ilog2() - self.first_row_bits() + 1) as usize
+    }
+
+    /// Rows the root indirect block can have: enough to span the heap's
+    /// whole `2^max_heap_size` address space.
+    fn max_root_rows(&self) -> usize {
+        (u32::from(self.max_heap_size) - self.first_row_bits() + 1) as usize
+    }
+
+    /// The largest object a direct block holds.
+    fn max_managed(&self) -> u32 {
+        (self.max_direct_block_size - self.dblock_header_size as u64) as u32
+    }
+}
+
+/// A direct block the packer filled.
+struct HeapDirectBlock {
+    size: u64,
+    heap_offset: u64,
+    data: Vec<u8>,
+}
+
+/// One entry of an indirect block.
+enum HeapSlot {
+    /// Not allocated (undefined address).
+    Empty,
+    /// Index into the packer's direct blocks.
+    Direct(usize),
+    Indirect(HeapIndirectBlock),
+}
+
+struct HeapIndirectBlock {
+    heap_offset: u64,
+    nrows: usize,
+    /// `nrows * width` entries, row-major.
+    slots: Vec<HeapSlot>,
+}
+
+impl HeapIndirectBlock {
+    fn own_size(&self, geom: &HeapGeometry) -> u64 {
+        (geom.iblock_fixed_size + self.slots.len() * OFFSET_SIZE as usize) as u64
+    }
+
+    /// Bytes of this block and everything below it.
+    fn subtree_size(&self, geom: &HeapGeometry, blocks: &[HeapDirectBlock]) -> u64 {
+        self.own_size(geom)
+            + self
+                .slots
+                .iter()
+                .map(|s| match s {
+                    HeapSlot::Empty => 0,
+                    HeapSlot::Direct(i) => blocks[*i].size,
+                    HeapSlot::Indirect(ib) => ib.subtree_size(geom, blocks),
+                })
+                .sum::<u64>()
+    }
+
+    /// Append this block at `addr` (= `out`'s current end, relative to the
+    /// same base as `frhp_addr`), then its children in entry order.
+    fn emit(
+        &self,
+        geom: &HeapGeometry,
+        blocks: &[HeapDirectBlock],
+        frhp_addr: u64,
+        addr: u64,
+        out: &mut Vec<u8>,
+    ) {
+        let block_offset_bytes = (geom.max_heap_size as usize).div_ceil(8);
+        let start = out.len();
+        out.extend_from_slice(b"FHIB");
+        out.push(0); // version
+        write_offset(out, frhp_addr, OFFSET_SIZE);
+        out.extend_from_slice(&self.heap_offset.to_le_bytes()[..block_offset_bytes]);
+        let mut child = addr + self.own_size(geom);
+        for s in &self.slots {
+            match s {
+                HeapSlot::Empty => write_undef_offset(out, OFFSET_SIZE),
+                HeapSlot::Direct(i) => {
+                    write_offset(out, child, OFFSET_SIZE);
+                    child += blocks[*i].size;
+                }
+                HeapSlot::Indirect(ib) => {
+                    write_offset(out, child, OFFSET_SIZE);
+                    child += ib.subtree_size(geom, blocks);
+                }
+            }
+        }
+        let checksum = crate::checksum::jenkins_lookup3(&out[start..]);
+        out.extend_from_slice(&checksum.to_le_bytes());
+
+        let mut child = addr + self.own_size(geom);
+        for s in &self.slots {
+            match s {
+                HeapSlot::Empty => {}
+                HeapSlot::Direct(i) => {
+                    let b = &blocks[*i];
+                    let d = out.len();
+                    out.extend_from_slice(b"FHDB");
+                    out.push(0); // version
+                    write_offset(out, frhp_addr, OFFSET_SIZE);
+                    out.extend_from_slice(&b.heap_offset.to_le_bytes()[..block_offset_bytes]);
+                    let cksum_pos = out.len();
+                    out.extend_from_slice(&[0u8; 4]); // checksum placeholder
+                    out.extend_from_slice(&b.data);
+                    out.resize(d + b.size as usize, 0);
+                    let cksum = crate::checksum::jenkins_lookup3(&out[d..]);
+                    out[cksum_pos..cksum_pos + 4].copy_from_slice(&cksum.to_le_bytes());
+                    child += b.size;
+                }
+                HeapSlot::Indirect(ib) => {
+                    ib.emit(geom, blocks, frhp_addr, child, out);
+                    child += ib.subtree_size(geom, blocks);
+                }
+            }
+        }
+    }
+}
+
+/// Packs objects into a heap's doubling table in heap-offset order.
+struct HeapPacker<'a> {
+    geom: &'a HeapGeometry,
+    objects: &'a [Vec<u8>],
+    /// The next object to place.
+    next: usize,
+    blocks: Vec<HeapDirectBlock>,
+    /// Each object's (heap offset, length).
+    obj_loc: Vec<(u64, u64)>,
+}
+
+impl HeapPacker<'_> {
+    /// Fill an indirect block at `heap_offset` with `nrows` rows, or, for the
+    /// root (`None`), with as many rows as the objects need.
+    fn fill(
+        &mut self,
+        heap_offset: u64,
+        nrows: Option<usize>,
+    ) -> Result<HeapIndirectBlock, FormatError> {
+        let geom = self.geom;
+        let width = geom.width as usize;
+        let mut slots = Vec::new();
+        let mut off = heap_offset;
+        let mut row = 0usize;
+        while self.next < self.objects.len() && nrows.is_none_or(|n| row < n) {
+            if nrows.is_none() && row >= geom.max_root_rows() {
+                return Err(FormatError::SerializationError(format!(
+                    "fractal heap: {} objects do not fit its {}-bit address space",
+                    self.objects.len(),
+                    geom.max_heap_size
+                )));
+            }
+            let size = geom.row_size(row);
+            for _ in 0..width {
+                if self.next == self.objects.len() {
+                    slots.push(HeapSlot::Empty);
+                } else if row < geom.max_direct_rows() {
+                    slots.push(self.fill_direct(off, size));
+                } else {
+                    let child = self.fill(off, Some(geom.rows_for_size(size)))?;
+                    let used = child.slots.iter().any(|s| !matches!(s, HeapSlot::Empty));
+                    slots.push(if used {
+                        HeapSlot::Indirect(child)
+                    } else {
+                        HeapSlot::Empty
+                    });
+                }
+                off += size;
+            }
+            row += 1;
+        }
+        let nrows = nrows.unwrap_or(row);
+        slots.resize_with(nrows * width, || HeapSlot::Empty);
+        Ok(HeapIndirectBlock {
+            heap_offset,
+            nrows,
+            slots,
+        })
+    }
+
+    /// Fill the direct block at `heap_offset` with as many of the next
+    /// objects as fit; leave it unallocated if not even the next one does.
+    fn fill_direct(&mut self, heap_offset: u64, size: u64) -> HeapSlot {
+        let header = self.geom.dblock_header_size;
+        let capacity = size as usize - header;
+        let mut data = Vec::new();
+        while let Some(obj) = self.objects.get(self.next) {
+            if data.len() + obj.len() > capacity {
+                break;
+            }
+            self.obj_loc[self.next] =
+                (heap_offset + (header + data.len()) as u64, obj.len() as u64);
+            data.extend_from_slice(obj);
+            self.next += 1;
+        }
+        if data.is_empty() && self.objects.get(self.next).is_some_and(|o| !o.is_empty()) {
+            return HeapSlot::Empty;
+        }
+        self.blocks.push(HeapDirectBlock {
+            size,
+            heap_offset,
+            data,
+        });
+        HeapSlot::Direct(self.blocks.len() - 1)
     }
 }
 
@@ -661,6 +824,8 @@ struct WriteFrhp {
     free_space: u64,
     managed_space: u64,
     alloc_space: u64,
+    /// Heap offset of the next direct block to allocate.
+    alloc_iter: u64,
     nobjects: u64,
     table_width: u16,
     starting_block_size: u64,
@@ -685,7 +850,7 @@ fn write_frhp(p: WriteFrhp) -> Vec<u8> {
     write_undef_offset(&mut frhp, OFFSET_SIZE); // free_space_mgr_addr
     write_length(&mut frhp, p.managed_space, LENGTH_SIZE); // managed_space_in_heap
     write_length(&mut frhp, p.alloc_space, LENGTH_SIZE); // allocated_managed_space
-    write_length(&mut frhp, 0, LENGTH_SIZE); // dblock_alloc_iter
+    write_length(&mut frhp, p.alloc_iter, LENGTH_SIZE); // dblock_alloc_iter
     write_length(&mut frhp, p.nobjects, LENGTH_SIZE); // managed_objects_count
     write_length(&mut frhp, 0, LENGTH_SIZE); // huge_objects_size
     write_length(&mut frhp, 0, LENGTH_SIZE); // huge_objects_count
@@ -704,7 +869,10 @@ fn write_frhp(p: WriteFrhp) -> Vec<u8> {
 }
 
 /// Build dense attribute storage for a set of attributes.
-pub(crate) fn build_dense_attrs(attrs: &[AttributeMessage], base_address: u64) -> DenseAttrBlob {
+pub(crate) fn build_dense_attrs(
+    attrs: &[AttributeMessage],
+    base_address: u64,
+) -> Result<DenseAttrBlob, FormatError> {
     // Dense attrs use v3 attribute messages (adds character set encoding byte).
     let serialized: Vec<Vec<u8>> = attrs.iter().map(|a| a.serialize_v3(LENGTH_SIZE)).collect();
 
@@ -717,7 +885,7 @@ pub(crate) fn build_dense_attrs(attrs: &[AttributeMessage], base_address: u64) -
     let ls = LENGTH_SIZE as usize;
 
     // Attribute heaps use max_heap_size 40 / heap ID length 8 (matching libhdf5).
-    let heap = build_single_block_fractal_heap(&serialized, base_address, 40, 8);
+    let heap = build_single_block_fractal_heap(&serialized, base_address, 40, 8)?;
     let frhp_addr = heap.frhp_addr;
     let btree_addr = heap.btree_addr;
     let heap_id_length = heap.heap_id_length;
@@ -781,10 +949,10 @@ pub(crate) fn build_dense_attrs(attrs: &[AttributeMessage], base_address: u64) -
 
     let attr_info = serialize_attribute_info(frhp_addr, bthd_addr);
 
-    DenseAttrBlob {
+    Ok(DenseAttrBlob {
         attr_info_message: attr_info,
         blob,
-    }
+    })
 }
 
 // ---- Dense link blob ----
@@ -873,7 +1041,7 @@ pub(crate) fn build_dense_links(
 
     // libhdf5's link heap uses max_heap_size 32 / heap ID length 7 (vs 40/8 for
     // attributes), giving a 7-byte heap ID and an 11-byte type-5 record.
-    let heap = build_single_block_fractal_heap(&serialized, base_address, 32, 7);
+    let heap = build_single_block_fractal_heap(&serialized, base_address, 32, 7)?;
     let heap_id_length = heap.heap_id_length;
 
     // Type 5 records: hash(4) + heap_id. The B-tree's key is the name hash,
@@ -1405,7 +1573,9 @@ impl FileWriter {
             .enumerate()
             .map(|(gi, g)| {
                 let dummy_links = g.link_messages(&[], &[]);
-                let attr_blob = group_dense[gi].then(|| build_dense_attrs(&g.attrs, 0));
+                let attr_blob = group_dense[gi]
+                    .then(|| build_dense_attrs(&g.attrs, 0))
+                    .transpose()?;
                 let li = if group_links_dense[gi] {
                     serialize_link_info(
                         g.track_order.then_some(0),
@@ -1439,7 +1609,9 @@ impl FileWriter {
         let mut dummy_blobs: Vec<DataBlob> = Vec::new();
         let mut dummy_cursor = 0u64;
         for (i, d) in all_ds.iter().enumerate() {
-            let dense_blob = ds_dense[i].then(|| build_dense_attrs(&d.attrs, 0));
+            let dense_blob = ds_dense[i]
+                .then(|| build_dense_attrs(&d.attrs, 0))
+                .transpose()?;
             if is_vds[i] {
                 // VDS: dummy OH with address 0 to get the OH size. The global
                 // heap blob will be placed after the OHs in pass 2.
@@ -1563,7 +1735,7 @@ impl FileWriter {
                 group_link_blob_addrs.push(None);
             }
             if group_dense[gi] {
-                let blob = build_dense_attrs(&g.attrs, cursor2 as u64);
+                let blob = build_dense_attrs(&g.attrs, cursor2 as u64)?;
                 cursor2 += blob.blob.len();
                 group_dense_blobs.push(Some(blob));
             } else {
@@ -1580,15 +1752,15 @@ impl FileWriter {
                 let addr = cursor2 as u64;
                 cursor2 += sz;
                 if ds_dense[i] {
-                    let blob = build_dense_attrs(&all_ds[i].attrs, cursor2 as u64);
+                    let blob = build_dense_attrs(&all_ds[i].attrs, cursor2 as u64)?;
                     cursor2 += blob.blob.len();
                     ds_dense_blobs.push(Some(blob));
                 } else {
                     ds_dense_blobs.push(None);
                 }
-                addr
+                Ok(addr)
             })
-            .collect();
+            .collect::<Result<_, FormatError>>()?;
 
         let mut ds_blobs2: Vec<DataBlob> = Vec::new();
         let global_align_threshold = self.alignment_threshold;

crates/clawhdf5-tools/tests/h5rs_interop.rs

@@ -858,3 +858,34 @@ fn check_and_dump_files_with_nested_groups_and_links() {
         assert_eq!(stdout(&ours), r, "{name}");
     }
 }
+
+#[test]
+fn check_files_with_big_dense_storage() {
+    // Dense links and attributes past the 512 KiB the root indirect block's
+    // direct blocks hold: the heap then needs child indirect blocks, which
+    // the writer used to write as direct blocks ("fractal heap indirect
+    // block: bad signature").
+    use clawhdf5::{AttrValue, FileBuilder};
+    let dir = tempfile::tempdir().unwrap();
+    let mut b = FileBuilder::new();
+    let x = b.create_dataset("x");
+    x.with_i32_data(&[7]);
+    for i in 0..150usize {
+        let len = if i % 3 == 0 { 7_000 } else { 1 + i };
+        x.set_attr(
+            &format!("a{i:03}"),
+            AttrValue::F64Array(vec![i as f64; len]),
+        );
+    }
+    let mut g = b.create_group("g");
+    for i in 0..40_000 {
+        g.add_hard_link(&format!("link_{i:06}_{}", "x".repeat(88)), "/x");
+    }
+    b.add_group(g.finish());
+    let p = dir.path().join("big.h5").to_string_lossy().into_owned();
+    b.write(&p).unwrap();
+    // Structure only: `--data` looks every link up by a linear scan.
+    let o = h5rs(&["check", &p]);
+    assert_eq!(code(&o), 0, "{p}:\n{}", stdout(&o));
+    assert!(stdout(&o).contains("no problems found"), "{}", stdout(&o));
+}

crates/clawhdf5/tests/writer_groups_interop.rs

@@ -638,3 +638,130 @@ fn a_group_attribute_set_again_takes_the_new_value() {
     let f = File::open(&path).unwrap();
     assert!(matches!(f.root().attrs().unwrap()["v"], AttrValue::I64(2)));
 }
+
+// ---- big dense storage: child indirect blocks in the fractal heap ----
+
+/// A name `len` bytes long, unique per `i`.
+fn long_name(i: usize, len: usize) -> String {
+    let n = format!("link_{i:06}_");
+    format!("{n}{}", "x".repeat(len - n.len()))
+}
+
+#[test]
+fn dense_links_past_the_direct_blocks_of_the_root() {
+    skip_if_no_python!();
+    // A dense group's links live in a fractal heap whose root indirect
+    // block holds direct blocks up to 64 KiB: 512 KiB of link messages.
+    // Rows past that are child indirect blocks. The writer used to write
+    // them as direct blocks, which libhdf5 cannot read ("incorrect metadata
+    // checksum"), from about 17 000 links with 20-byte names.
+    // `g` crosses the first boundary (0.6 MB of links); `deep` has 65 535
+    // links of about 110 bytes (7 MB), so its heap reaches the child indirect
+    // blocks that hold indirect blocks themselves.
+    let dir = tempfile::tempdir().unwrap();
+    let mut b = FileBuilder::new();
+    b.create_dataset("x").with_i32_data(&[7]);
+    let mut g = b.create_group("g");
+    for i in 0..20_000 {
+        g.create_dataset(&format!("dataset_number_{i:06}"))
+            .with_i32_data(&[i]);
+    }
+    b.add_group(g.finish());
+    let mut g = b.create_group("deep");
+    g.track_order(true);
+    for i in 0..usize::from(u16::MAX) {
+        g.add_hard_link(&long_name(i, 100), "/x");
+    }
+    b.add_group(g.finish());
+    let path = write(&dir, "big_links.h5", b);
+
+    let out = h5py(
+        &path,
+        "with h5py.File(path, 'r') as f:\n\
+         \x20   g, d = f['g'], f['deep']\n\
+         \x20   names = list(g)\n\
+         \x20   dn = list(d)\n\
+         \x20   print(json.dumps([len(names), names[-1], int(g[names[-1]][0]),\n\
+         \x20       sum(int(g[n][0]) for n in names), len(dn), dn[0][:12], dn[-1][:12],\n\
+         \x20       int(d[dn[-1]][0]), h5py.h5o.get_info(f['x'].id).rc]))",
+    );
+    assert_eq!(
+        out,
+        r#"[20000, "dataset_number_019999", 19999, 199990000, 65535, "link_000000_", "link_065534_", 7, 65536]"#
+    );
+    h5dump_ok(&path);
+    let f = File::open(&path).unwrap();
+    let g = f.group("g").unwrap();
+    assert_eq!(g.datasets().unwrap().len(), 20_000);
+    assert_eq!(
+        g.dataset("dataset_number_019999")
+            .unwrap()
+            .read_i32()
+            .unwrap(),
+        [19999]
+    );
+    let d = f.group("deep").unwrap();
+    assert_eq!(d.datasets().unwrap().len(), usize::from(u16::MAX));
+    assert_eq!(
+        d.dataset(&long_name(65_534, 100))
+            .unwrap()
+            .read_i32()
+            .unwrap(),
+        [7]
+    );
+
+    // libhdf5 can add to and delete from the heap. It could not when the
+    // header's block allocation offset was 0: its next block overwrote the
+    // first ("bad version number for message").
+    let out = h5py(
+        &path,
+        "with h5py.File(path, 'r+') as f:\n\
+         \x20   f['g']['zz_new'] = np.arange(3)\n\
+         \x20   del f['g/dataset_number_000005']\n\
+         with h5py.File(path, 'r') as f:\n\
+         \x20   print(json.dumps([len(f['g']), int(f['g/zz_new'][2]),\n\
+         \x20       int(f['g/dataset_number_019998'][0])]))",
+    );
+    assert_eq!(out, r#"[20000, 2, 19998]"#);
+    h5dump_ok(&path);
+}
+
+#[test]
+fn dense_attributes_past_the_direct_blocks_of_the_root() {
+    skip_if_no_python!();
+    // Dense attributes share the heap writer. 150 attributes of up to 56 KB
+    // (8 MB) need child indirect blocks, and a big attribute after small
+    // ones must skip the small blocks rather than overrun one.
+    let dir = tempfile::tempdir().unwrap();
+    let mut b = FileBuilder::new();
+    let ds = b.create_dataset("x");
+    ds.with_i32_data(&[1]);
+    for i in 0..150usize {
+        let len = if i % 3 == 0 { 7_000 } else { 1 + i };
+        let v: Vec<f64> = (0..len).map(|k| (i * 100_000 + k) as f64).collect();
+        ds.set_attr(&format!("a{i:03}"), AttrValue::F64Array(v));
+    }
+    let path = write(&dir, "big_attrs.h5", b);
+
+    let out = h5py(
+        &path,
+        "with h5py.File(path, 'r') as f:\n\
+         \x20   a = f['x'].attrs\n\
+         \x20   ok = all(np.array_equal(a['a%03d' % i],\n\
+         \x20       np.arange(7000 if i % 3 == 0 else 1 + i) + i * 100000) for i in range(150))\n\
+         \x20   print(json.dumps([len(a), ok]))",
+    );
+    assert_eq!(out, "[150, true]");
+    h5dump_ok(&path);
+    let f = File::open(&path).unwrap();
+    let attrs = f.dataset("x").unwrap().attrs().unwrap();
+    assert_eq!(attrs.len(), 150);
+    for i in [0usize, 1, 147, 149] {
+        let len = if i % 3 == 0 { 7_000 } else { 1 + i };
+        let want: Vec<f64> = (0..len).map(|k| (i * 100_000 + k) as f64).collect();
+        match &attrs[&format!("a{i:03}")] {
+            AttrValue::F64Array(v) => assert_eq!(*v, want, "a{i:03}"),
+            other => panic!("a{i:03}: {other:?}"),
+        }
+    }
+}