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

Two further variables control which MPI rank owns a term; they belong with the distribution axis and are documented under Rank routing below.

# 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.

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.

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.

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.