Parallelism and distribution

monoprop scales the same way on one node and on many: it partitions the operator into disjoint pieces and applies each gate to every partition in lock-step, exchanging only the terms that cross a partition boundary. Two axes compose:

• Across cores (default): one single-threaded partition per physical core, within a single process. This is on automatically, and can be configured with environment variables (see below).
• Across nodes (opt-in): MPI ranks, each of which partitions further across its own cores. Distributing the operator and its graph across ranks lets you simulate and variationally optimise systems larger than one node's memory.

Operator partitioning (single node)

By default monoprop splits the operator into one partition per physical core and gives each partition a pinned worker thread that runs it serially. Each partition keeps its own small, cache-resident term index; a gate is applied to all partitions at once, synchronised by a lightweight barrier, and anticommuting terms whose partner lives in another partition are resolved through a per-gate exchange.

A partition addresses its own terms with a 32-bit index, so one partition holds at most 2^32 (about 4.3 billion) terms. That is a few hundred gigabytes of resident operator in a single partition, so the ceiling is reached by running too few partitions rather than by problem size, and it is not a build option. An append that would cross it is refused before it happens — a RuntimeError naming the ceiling and the term counts, with the operator left as it was, rather than a silent wrap. The fix is more partitions or more MPI ranks — both divide the terms per partition — or a larger lower_atol (see Truncation and cutoffs).

Runtime environment variables

Three further variables belong with the distribution axis rather than this one and are documented under Rank routing below: two control which MPI rank owns a term, and one reports the gate exchange's wire volume.

# Run 8 partitions instead of one-per-core:
export monoprop_NUM_THREADS=8

Where the operator's bytes are

operator_memory_breakdown() names the resident bytes of one rank's operator. total_bytes() sums the fields that are live at quiescence; the d_ keys are counts, subsets of a field, or quantities no resting field can hold, so folding one into the total would double-count.

The row store and the anticommutation index's dense columns are held as fixed-size chunks taken from a pool, not as vectors: growth appends a chunk, nothing is copied, and the only slack is the tail of the last chunk. Four keys describe that storage. d_terms_slack_bytes is that tail. d_pool_mapped_bytes is what the pools have mapped -- a pool maps whole arenas and keeps one as long as a single chunk in it is live, so it sits at or above what the chunk-counting fields price, and d_pool_free_chunk_bytes is the part of it not currently handed out. Those bytes are faulted only as they are written, but the mapping is what the kernel's high-water mark is charged for once they have been, which is why the ledger's total can sit below the resident set with nothing missing from any named field.

d_row_wide_rows, d_row_inline_width and d_row_restrides describe the row layout. A row is a popcount slot followed by its set positions, at one width paid for by every row; rows wider than that but still within the cutoff's structural bound go to a second tier at a fixed wider stride, and rows over the bound to a lossless side-map. The width is predicted from the model, so a wrong guess costs one restride -- the store re-lays every row at the bound between gates, keeping every row index -- and never a wrong answer. A nonzero d_row_restrides says the prediction was too narrow for this model; a wide-row count that is a large share of the term count says the same before it happens.

indexing_bytes is the key -> row table over the stored terms: four-byte slots at a load of at most 0.7, so 5.7-11.4 bytes per term depending on where the term count sits between two doublings, and 0 until the first gate materialises it. See Rank routing for the key it hashes.

d_op_coeffs_slack_bytes is the most capacity the coefficient array held beyond its live rows at any point in the last propagate or build_graph call. The array grows at the row store's own 1.5x policy and is shrunk to fit when the call ends, which is the only moment a caller can read the ledger -- so it is reported as a high-water mark taken at the growth sites rather than as a measurement of the array as it stands, which would always be 0. Across partitions the marks are summed rather than maxed: the partitions grow together within a call, so the sum is the figure a per-process footprint wants, and it errs high.

d_wire_staging_bytes is not the operator's memory at all but the transport's, which is why it is outside the total: the in-process transports keep a payload staging buffer per direction, grown to the widest message the run has needed and never shrunk, plus the count and displacement tables that are fixed by the rank and partition counts. On a run whose widest gate came early those buffers are resident for the whole of it while nothing resting names them. One transport serves all the partitions of a rank, so it is a per-rank figure: exactly one partition reports it and a sum over partitions counts a rank's staging once. A pure-MPI rank reports 0, because there its payload buffers belong to the in-flight exchange handle and die with the gate that opened the round.

gate_scratch_bytes is the per-gate layer-build scratch, which the propagator owns so that its capacity survives the gate: the anticommuting fold's words, the protocol's per-row marks, the join's hit slots, the miss stage and the decoded incoming records. It carries no state between gates. d_gate_buffers_hwm_bytes is outside the total for the opposite reason to the fields above: it is not a subset of any of them but a quantity no resting field can hold. The records a gate puts on the wire, the ones it receives, the responses it stages and the answers it applies are all freed when the gate returns, so by the time a caller can ask, nothing of them is left to measure -- yet on a wide gate they are the largest thing the call adds to resident memory. The engine therefore stamps their combined size at the two instants where the set is widest, and the field carries the maximum over the gates of the last propagate or build_graph. Across partitions the values are summed, which is an upper bound rather than a simultaneous figure: each partition's peak is over its own timeline.

Peak memory of large runs: pin glibc's mmap threshold

Every gate allocates and frees a few large transient buffers: the records it stages, sends, receives and answers. glibc's allocator reacts to the first free of such a buffer by raising its dynamic mmap threshold (up to 32 MiB) and its trim threshold, so from then on the widest gate's transients are served from the heap and kept on free lists the allocator cannot return to the kernel. On a 250 M-term Hubbard run that retention is about a fifth of the resident set. Pinning the threshold disables the escalation, so every transient above it is a private mapping that is unmapped on free:

export GLIBC_TUNABLES=glibc.malloc.mmap_threshold=1048576   # 1 MiB
mpirun -x GLIBC_TUNABLES ...                                # MPI ranks must inherit it

Measured on a 16-core workstation at 250 M terms (paired interleaved repetitions): peak resident high-water mark 0.95x in-process with 16 partitions and 0.91x (ranks summed) at 4 ranks x 4 partitions, for about 1 % more wall time; a 128 KiB threshold takes a further 2 % of memory for 2-4 % more time. The cost is page faults on the re-mapped transients, so the setting pays only when a gate's compute amortises them: use it at or above ~10^8 terms with several partitions, and not on small single-core runs, where it costs ~20 % time. Results are bit-identical either way -- this changes where bytes live, not what is computed.

Placement report

Building a propagator writes one COMMPLACE line per rank to stderr, naming the CPUs the launcher gave that rank and whether co-located ranks got disjoint masks. It is report-only — no placement decision reads it — and there is no knob: redirect stderr to drop it.

COMMPLACE rank=0 node_rank=0 node_size=2 masks=private cpus=64 node_cpus=128 cpu_list=0-63

masks is private when the co-located ranks' affinity masks are pairwise disjoint, shared when two ranks can land on the same CPU, alone when this rank is the only one on its host (which is not evidence a multi-rank launcher bound correctly), and unknown when a mask did not fit the exchanged window.

Pinning each partition to a core is not configurable: leaving placement to the launcher measured propagate[hubbard] 2.90x slower, so the disabled arm is gone.

MPI distribution (multi-node)

MPI partitions the operator and graph across ranks, composing with per-rank partitioning into one flat world of R × S partitions (R ranks, S partitions each). MPI communication is serialised through each rank's first partition, bracketed by the intra-rank barriers.

Rank routing

For $N$ modes, monomials under a gate form the group $(\mathbb{F}_2^{2N}, \oplus)$: a generator $G$ sends $M \mapsto M \oplus G$. The rank index is chosen to be a homomorphism of that group — $h(M) = \bigoplus_{i \in \mathrm{supp}(M)} v_i$, over one fixed vector $v_i$ per Majorana slot, with $h_d$ its low $d = \log_2 R$ bits — so that

Three consequences. A rank owning the fibre $h_d^{-1}(r)$ sends every query for $G$ to $r \oplus h_d(G)$: one peer, independent of $R$ and of $|\mathrm{supp}\,G|$, where a full-avalanche hash sprays the same queries across all $R$ ranks. XOR is an involution, so the peer relation is symmetric and both sides derive the pairing without communicating. And the fibres are cosets of $\ker h_d$, all of size $2^{2N-d}$, so a uniformly drawn monomial is balanced by construction; imbalance can only come from the operator's support being non-uniform, which is why balance is measured rather than proved — rank occupancy max/mean 1.001 at $R = 128$, with every rank used.

That one peer is what the exchange is then allowed to spend. A gate's query round carries a PeerPlan -- the rank shift $h_d(G)$, derived once per generator where the generator is held -- and the transports read it rather than a rank count: one Isend/Irecv pair with the peer instead of an MPI_Alltoall on the counts and an MPI_Alltoallv on the payload, and a plain memcpy when the shift is zero and the peer is this rank itself. The reply round retraces the queries, and XOR is an involution, so it takes the same plan. Under the hybrid the larger win is not the messages but the serial sweeps: the count transpose, the two staging sizers, the recv column and the scatter are all $O(R S^2)$ and all run on partition 0 while the other $S - 1$ partitions park at a barrier, and restricting them to the reachable ranks makes them $O(S^2)$. Graph replay is narrowed the same way, without a plan: its count matrix is symmetric, so what a rank sends a peer is that peer's receive count and both ends drop the same legs on the same value.

Which transport a round takes must be the same on every rank, or a rank inside a collective waits forever on ranks that chose point-to-point. So it is never a predicate on a rank's own traffic -- rows vary, and any threshold on one straddles -- but a function of the resolved routing mode and the rank count alone, the same two things the propagator's construction-time agreement check reduces across the world. A wrong shift that every rank agrees on does not hang; it silently drops the blocks outside the peer set, and debug builds assert against exactly that.

The one condition is that the per-generator shifts $\{h_d(G)\}$ span $\mathbb{F}_2^{d}$. If they span only $\rho < d$ dimensions, the ranks a query can reach form a coset of a $\rho$-dimensional subspace: only $2^{\rho}$ of the $2^{d}$ are ever used and the rest stay empty. That is a load-balance failure and not a wrong answer, so it is reported as one COMMROUTE line per rank rather than enforced at runtime — measured $\rho = 32$ over the 60-site Hubbard's 416 distinct shifts, against the $d = 7$ that $R = 128$ needs.

The same image $h(M)$ — the term's fingerprint, computed by XOR-ing one label per set position, so a packed row need never be expanded into a bitset to be routed — also places the term within a rank. With $S$ a power of two the partition takes the next $\log_2 S$ bits of $h$, so the whole flat slot is linear in $M$ and every partition has exactly one peer partition per generator, inside the rank as well as across ranks:

With $S$ not a power of two only the rank level has that structure: the partition is the mixed fingerprint modulo $S$, balance only, and a generator's queries from one partition land on the $S$ partitions of one peer rank.

The fingerprint serves one more purpose. Its top 32 bits are the join key, $\mathrm{key}(M) = h(M) \gg 32$, which a persistent open-addressing table over every stored term is hashed by (indexing_bytes: four-byte slots holding a row index below a hash prefilter, at a load of at most 0.7). No key is resident: the table folds a row's key off the row's own positions when it indexes it -- streaming the store on a rebuild, and reading the rows a gate has just minted on an append -- and the emit path folds the partner's, off the positions the partner product has just produced. A record names its partner by that key, and the gate's join is one batched, prefetch-pipelined probe per record -- proportional to the records, not to the anticommuting set -- with every key match confirmed against the row's positions.

Taking the key as a plain projection rather than a mixed one keeps it linear, $\mathrm{key}(M \oplus G) = \mathrm{key}(M) \oplus \mathrm{key}(G)$, so a rank folds the same number off positions it decoded from the wire that the sender would fold off the term itself. It projects the high half because the low bits are the routing bits, so a partition's rows do not share a key prefix. Collisions become structured -- two terms share a key exactly when their symmetric difference has key zero -- but the labels are mix64 outputs, so the projection is uniform, and a key match was never more than a prefilter: the confirm against the positions is what decides a partner, so a collision costs a compare and never a wrong answer.

Every gate is then one exchange round for the terms that rotate. A term $M$ that anticommutes with $G$ and passes the rotation gate sends one record — the partner's positions, the rotation phase $\varphi(M, G)$, one bit saying whether $M$'s own side rotates, and in the fused path $M$'s pre-gate coefficient — to the partner's owner, all of them in one alltoallv. The owner of $M \oplus G$ joins the key against its own stored terms: a hit is the partner, and the pair rotates if either side asked for it, each owner applying the record it received ($\varphi(M \oplus G, G) = -\varphi(M, G)$, so the two adds are exactly the two-term rotation); a miss mints the partner. Graph mode records the same joins as the layer's in/out endpoint lists, so the replay pairs them positionally without exchanging indices.

A pair where both sides rotate is answered by that one round alone: each side receives the other's record. Where only one side rotates it is not, because a term below lower_atol still owns a coefficient its rotating partner needs, and only that term can supply it. Making every such term send a record of its own would put the whole anticommuting set on the wire and into the join — in the lower_atol regime, six times the records for the same rotations. So the hit carries the answer back instead: a record that lands on a silent row applies its own half there and returns that row's coefficient in a second, small message, naming the record by its position in the sender's stream. The sender applies its half from the answer. That makes three distinguishable outcomes for a record, and they are what tells a tracked partner from an absent one: the partner's own record arrived, or an answer arrived, or neither — and neither means absent, because only the partner's owner could have replied at all. Every add is the same add the two-pass protocol made, with the same value on the same side of the pair, so the round-and-a-half protocol reproduces it bit for bit at the volume of the rotations rather than of the anticommuting set.

Graph mode has no coefficients and therefore no silent terms to answer: it keeps the symmetric predicate, sending for every anticommuting term whose partner is structurally admissible, which is what its positional in/out pairing is proved on.

Linear routing is a switch and not a dial: the rank index takes every bit of $h$ or none of them, and none of them is splitmix, which reproduces the dense hash % (R × S) bit for bit. It therefore requires $R$ to be a power of two; any other rank count has no XOR structure to route by and is raised at propagator construction rather than silently routed onto a subspace of the ranks. $R = 1$ is a power of two, takes no rank bit, and so is the dense case already.

A related distributed scheme maps an index by summing its $k$-bit blocks modulo the rank count (Broers et al., 2026). That sum is additive modulo that count while the gate acts by XOR, so the carries make $f(I \oplus J)$ differ from $f(I)$ by $\pm 2^{J_j \bmod k}$ terms whose signs depend on the bits of $I$, bounding the destinations at $2^{2|J|+1}$ rather than collapsing them to one. Being linear over the same group the gate acts by is what turns that bound into an identity.

Single-node (MPI.COMM_SELF)

from mpi4py import MPI
sim = MajoranaPropagator(..., comm=MPI.COMM_SELF)

Multi-node (MPI.COMM_WORLD)

Replace the communicator and launch with mpiexec:

from mpi4py import MPI
sim = MajoranaPropagator(..., comm=MPI.COMM_WORLD)

mpiexec -n 8 uv run python your_script.py

A pure-MPI rank keeps its whole share in one partition unless monoprop_NUM_THREADS is set, so an MPI user who has not asked for threads gets one partition per rank.

Enabling MPI

MPI is off by default, so the prebuilt PyPI wheels run single-rank and the communicators above only distribute work after a from-source build with MPI enabled. See Building from source for the full build instructions.