GDB: sharded determinant basis, numpy marshalling, drivers and first tests

[Sophia Wen (hfwen0502)]

Makes GDB usable for a sparse, distributed subspace. Experimental — GDB moves fast upstream and this will follow it.

Stacked on #44. Base is examples-reorg, so this diff is only the GDB work; GitHub retargets to main once #44 merges. Review #44 first.

Why

gdb_diag required every rank to pass the whole determinant list, so b_comm_size had to be 1. That cascaded: b is the only dimension that shards memory (h_comm is a row stride plus an allreduce; MakeHelpers ignores it), t <= b is structural, and Thrust throws unless h_comm_size == 1 — so multi-GPU GDB was unreachable and memory did not scale at all. Upstream's own run.sh uses --b_comm_size 2; the restriction was ours, not SBD's.

What changed

• gdb_diag takes a per-rank shard. b_comm_size == 1 keeps today's meaning (whole basis); > 1 means this rank's shard, verified globally rather than trusted — shard identity across h_comm, local sort, cross-shard disjointness. All validation votes through MPI_Allreduce before any rank throws, so a bad input cannot deadlock.
• numpy (N, words) in and out, via c_style | forcecast so existing list callers still work. Removes two of three copies; adds from_strings for bulk packing instead of one from_string call per determinant (issue Duplicate-determinant handling differs between TPB and GDB; cost of the #29 label conversion at scale #31 flagged that at 39k).
• All six placement schemes as one determinant_distribution string. The three legacy booleans stay unexposed: they encode a four-way choice with a priority rule, and do_redist_alpha_eq defaults true so the others are dead unless it is zeroed.
• Guardrails for two silent failures: t > b segfaults upstream (a starved rank dereferences an empty exidx) and t*b not dividing the rank count silently builds ragged communicators. Both now refuse up front.
• Two drivers — run_gdb_diag.py and run_gdb_heatbath.py (a cutoff ladder; carryover_type 2/3 returns parents with candidates, so one round's result is the next subspace). In-memory throughout, no file round-trip.
• Per-spin electron-count validation. Feeding qiskit-addon-sqd's concatenated [beta | alpha] strings where GDB wants them interleaved used to diagonalize to a plausible wrong energy and then abort inside the expansion with std::out_of_range. Now refused with the cause named. The occupation density cannot catch it — permuting bits preserves how many are set.

Verification

GDB had no test anywhere before this, in the wrapper or upstream. Now:

• TPB is equivalent to GDB on subspaces both can express (an interleaved product), at 576, 59,536 and 1,000,000 determinants — two independent solvers, same Hilbert space, no pinned value needed.
• Published references: h2o-1em3 -76.23594663, h2o-1em4 -76.24295848, n2-1em4, n2-1em5, up to 16.3M determinants.
• Sharding invariance: same energy at b 1/2/4/8 and t 1/2, across all six schemes — placement must not change the answer.
• Driver tests run the scripts as subprocesses, so argument parsing and the defaults are exercised rather than bypassed; the h2o case asserts the same anchor the library test does.
• Multi-GPU GDB confirmed on 8x H100 with the devices actually busy, not inferred from a matching energy.

28 passed serial, 12 with --run-slow, 13 under mpirun -n 2.

Known gaps

RDMs (do_rdm=1) and savename are untested at b > 1; no multi-node run; DetBasisCommunicator leaks four comms per call upstream, which a long loop would feel. gpu-omp has no GDB kernels at all, so GDB on AMD is CPU-only — documented rather than worked around.

🤖 Generated with Claude Code

[Jim Garrison (garrison)]

Related to Qiskit/qiskit-addon-sqd#389

[Jim Garrison (garrison)]

This is well-grounded work and I'd like to get it in rather than let it sit. I had Claude check the upstream claims in the description against vendor/sbd-upstream at 93ebabe and they all hold: the exidx[0].slide starvation (gdb/helper.h:736-761), the Thrust h_size != 1 throw (gdb/mult_thrust.h:310-314), energy assigned only inside the two method branches (gdb/sbdiag.h:418, :501), do_shuffle parsed then never read, and h_comm_size shadowed at gdb/sbdiag.h:206. It also modelled DetBasisCommunicator's splits for several (ranks, b, t) combinations, and the documented shard index rank % b_comm_size is correct in each — that is the load-bearing claim of the whole contract.

The collective discipline looks right: every per-rank check votes through sbd_all_ranks_ok before any rank throws, and the checks that throw directly are functions of the config alone, which every rank passes identically. Declining to expose h_comm_size and do_shuffle, with the upstream lines explaining why, is defensible. The TPB-equals-GDB equivalence test is a stronger anchor than a pinned number.

Three things I'd like addressed, then this looks good to me.

1. bit_length=64 is undefined behavior on the paths this PR newly enables

The drivers default to --bit_length 64 and test_gdb_equivalence.py sets BIT_LENGTH = 64. bitadvance() (framework/bit_manipulation.h:373) computes (((size_t) 1) << bit_length) - 1, so at 64 the shift is UB and the mask collapses to 0 in practice. It is reached from mpi_redistribution and mpi_sort_bitarray, i.e. from the new count and count-sorted schemes. This is the hazard established in #12.

Claude's reading of why the suite still passes is worth repeating: _shard() in the test reproduces get_mpi_range (framework/mpi_utility.h:22) exactly, so redistribution's fast path (caop/basic/basis.h:54-75) returns early on every rank and bitadvance never executes. So count is nominally covered but its redistribution code never runs. A caller who shards unevenly — the realistic case, and the one the count scheme exists for — would reach it.

Could the drivers default to 62, and could one test use a deliberately uneven split so the fast path is bypassed? 62 rather than 63 because of the interaction with check_spin_weights below, and rather than a smaller value because it keeps h2o's 48 interleaved bits in a single word — det_vector::init_elem_size fixes ceil(2*norb / bit_length) process-wide, so a value implying two words would collide with test_reference_energies.py and test_sqd_integration.py when they share a process. This came from tracing the call chain rather than from an observed failure, so treat it as a thing to confirm rather than a reproduction.

2. check_spin_weights skips silently on any bit_length it cannot handle

check_spin_weights (run_gdb_diag.py:117, run_gdb_heatbath.py:236) returns early unless bit_length == 64, printing nothing. --bit_length is user-settable and examples/gdb/README.md:323 demonstrates --bit_length 20, so that path drops the validation today — no default change needed to reach it.

That seems worth failing loudly over, because the function's own rationale is that this error class is invisible downstream: as the description notes, the occupation density cannot catch a permuted bit order since permuting preserves how many bits are set. A check premised on "nothing else will catch this" should not opt itself out quietly.

The masks are also more capable than the guard suggests. 0x5555…/0xAAAA… are correct whenever a word boundary preserves the alpha/beta alternation, which holds for any even bit_length — 20 and 32 included — and breaks only for odd ones. So the condition looks like it wants to be bit_length % 2 rather than != 64, and the odd case should raise or warn rather than return. Worth noting that 63 is an odd value, which is why I suggested 62 above: pairing a 63 default with this guard as written would turn a silent skip into a silent miscomputation.

3. The MPI tests do not follow the _standalone / _mpi split

test_gdb_equivalence.py:234, :350, :372, :400, :424 and :443 are single @pytest.mark.mpi functions. pytest-mpi filters on that marker in opposite directions, so these never run under tox -e py — including test_gdb_placement_does_not_change_the_energy, where four of the six schemes are meaningful at one rank. The convention elsewhere in the suite is one shared private body wrapped in two thin functions (test_sqd_integration.py:200-210).

Smaller things

• README.md:263 still shows gdb_diag without the three new keyword arguments, and the Returns list below it omits local_dim, global_dim and determinant_distribution, which the docstring and the release note both document.
• The savename change from one file to one per b_comm position is user-visible and currently documented only in the docstring. It seems to belong in the release note's upgrade section next to the method rejection.
• MPI_UNSIGNED_LONG at python/bindings.cpp:908 and :910 is correct on LP64 but is the only place size_t is sent as that type; everything else in the file uses MPI_UNSIGNED_LONG_LONG.

Neither Claude nor I built or ran this branch, so the reported 28 / 12 / 13 counts are unverified on our side.

---

This review was drafted by Claude Opus 5 under my guidance.

[Sophia Wen (hfwen0502)]

Thanks for the careful review. All three points are addressed in edfc1aa, along with the smaller ones. Point 1 is confirmed, and point 2 turned out to go deeper than the spin check.

1. bit_length=64

Confirmed. bitadvance() is reached on the path you traced. With deliberately uneven shards (384/192, rebalanced to 288/288), mpi_redistribution runs, and the last rank's exclusive upper bound ("last determinant + 1") collapses to all zeros at 64. On clang/arm64 that happens not to change the placement for h2o, so the energy was still right. It is undefined behavior on a real path all the same.

• Both drivers default to 62, and test_gdb_equivalence.py uses 62.
• New test_gdb_count_schemes_rebalance_uneven_shards (MPI) hands count and count-sorted a lopsided split, then asserts the counts come back balanced. That proves the slow path ran. It also asserts the energy matches TPB.

2. check_spin_weights, and odd bit_length in SBD itself

The masks are now per word: word j starts at global bit j * bit_length, so its alpha mask is 0x5555… or 0xAAAA… depending on the parity of that start. The check works at every bit_length and never skips. I verified it against direct string counts at 20, 31, 62, 63 and 64, up to 4 words, in both drivers.

While testing odd sizes I found that SBD's GDB heatbath expansion breaks whenever bit_length is odd and a determinant spans two or more words. Two upstream functions assume an even size:

• gdb::getHalfDets (chemistry/gdb/helper.h:160), compiled in every build: half_bl = bit_length / 2, plus even-position masks applied per word.
• DetFromAlphaBeta under SBD_TRADMODE (chemistry/basic/determinants.h:20), and the copy in omp_offload.h: half = bit_length / 2. The #else variant is correct.

Python ports of both mis-convert 100% of h2o determinants at 25, 31 and 33, partly at 47, and none at even sizes or at 63, where h2o fits in one word. Observed with the h2o heatbath:

bit_length	SBD_TRADMODE on (our CPU/gpu-omp build)	SBD_TRADMODE off
25, 31, 33	crash	crash
47	crash	wrong energy, exit 0 (−76.0714804798 vs −76.0723973374)
20, 62	ok	ok
24, 30, 63	ok	not run

On 8×H100 it also fails on gpu and gpu-omp (std::bad_alloc, heap corruption). Plain gdb_diag stays correct at odd sizes, apparently because it only uses the half-determinants as internally consistent keys. So the drivers now require an even bit_length in [2, 62] and refuse anything else up front. gdb_diag itself still passes the value through.

3. _standalone / _mpi

Split for the sharded-basis, placement (all six schemes) and rank-tiling tests, which now run under tox -e py too. Three stay MPI-only:

• test_gdb_task_dimension_works_once_the_ring_has_stations and test_gdb_rejects_shards_that_are_not_disjoint need 4 and 2 ranks, and already skip below that.
• test_gdb_under_mpi_matches_tpb_on_the_same_subspace would duplicate test_gdb_matches_tpb_on_the_same_subspace at one rank; its docstring says so.

Running on GPUs also showed that two of these tests could never pass on Thrust: they put ranks on the helper dimension, which the binding refuses there by design. On Thrust they now assert that refusal.

Smaller things

• README gdb_diag shows the three keyword arguments and the local_dim / global_dim / determinant_distribution result keys.
• The release note has an upgrade entry for savename writing one file per b_comm position. Before this PR, b_comm_size had to be 1, so it was always a single {savename}000000.bin.
• bindings.cpp:908/910 now send SBD_MPI_SIZE_T rather than MPI_UNSIGNED_LONG_LONG. The buffers are std::vector<size_t>, and upstream already defines that macro from SIZE_MAX. The other sites in the file send variables declared unsigned long long, which is why they use the _LONG_LONG type.

Verification

Built and run this time:

• Laptop (CPU): 47 single-process tests pass; the MPI suite passes at 2, 3 and 4 ranks.
• 8×H100, all three backends:
  • 47 single-process tests pass on each.
  • The MPI suite passes at 2, 4 and 8 ranks.
  • h2o_1em4 as a 2,380,849-determinant product basis: GDB on gpu under equal-bra-a, count, count-sorted and grid-cyclic, plus cpu and gpu-omp, all give −76.2429584823, equal to TPB on the same subspace at the same tolerance. All 8 GPUs were busy during the Thrust run.
  • Fe4S4, 72-bit determinants packed into two 62-bit words:
    • Upstream's bundled files (59,536 dets) give −326.6982518821 under all four placements on 8 GPUs.
    • A 1,000,000-determinant product basis gives −326.6135721292 from GDB on gpu, cpu and gpu-omp, the same as TPB.
    • The heatbath ladder on 4 GPUs grows 59,536 → 636,338 → 806,932 determinants, and the energy decreases every round.

---

This reply was drafted by Claude Opus 5.5 under my guidance.