diff --git a/cpp/include/monoprop/MonomialPropagator.h b/cpp/include/monoprop/MonomialPropagator.h index b684ff58..72ead85f 100644 --- a/cpp/include/monoprop/MonomialPropagator.h +++ b/cpp/include/monoprop/MonomialPropagator.h @@ -365,6 +365,8 @@ class MonomialPropagator { // it captures this and rejects a later call once it moves, as it does for a rebuilt graph. size_t initial_operator_epoch_{0}; + bool routing_coverage_reported_{false}; // report_routing_coverage_ speaks once per propagator + size_t logical_num_modes_{NumModes}; CutoffType cutoff_type_; @@ -449,6 +451,10 @@ class MonomialPropagator { const VecD ¶meters, std::optional only_rotate_len_k) -> void; + // Do this call's generator shifts span log2(R)? If not, ranks receive nothing (routing::gf2_rank). + // Not beside check_routing_agreement: at construction the gate list does not exist yet. + auto report_routing_coverage_(const std::vector &majoranas) -> void; + template auto run_gate_loop_(const std::vector &majoranas, std::optional only_rotate_len_k, diff --git a/cpp/monoprop/Evolution.cpp b/cpp/monoprop/Evolution.cpp index 296082cf..bc0bb531 100644 --- a/cpp/monoprop/Evolution.cpp +++ b/cpp/monoprop/Evolution.cpp @@ -26,6 +26,7 @@ #include "monoprop/detail/evolution/CosineRecomputeCallbacks.h" #include "monoprop/detail/mpi/Exchange.h" #include "monoprop/detail/mpi/MPICompat.h" +#include "monoprop/detail/mpi/Routing.h" namespace monoprop { namespace { @@ -73,10 +74,23 @@ auto &acquire_flat_exchange_buffers() { } // A property of the communicator, not the layer: all ranks participate even at local total_count 0. +// Still true with the pairwise arm: `wire_bits` decides the transport for the whole communicator, so a +// rank that skipped the round strands the others whichever transport they are on. auto layer_exchange_participates(const mpi::Comm &comm) -> bool { return mpi::size(comm) != 1; } +// The transport gate for post_flat_alltoallv, and the ONLY thing allowed to choose it: rank-uniform by +// construction, because linear_bits_for reads the environment alone and check_routing_agreement +// allreduces exactly this number at construction and throws rather than proceed on a mismatch. Fanout 1 +// -- every generator's queries land on one destination rank -- is both why a layer's other legs are +// empty and why no rank can be on the other side of the branch. Any other geometry keeps the collective. +auto layer_exchange_wire_bits(const mpi::Comm &comm) -> int { + const auto ranks = static_cast(mpi::geometry(comm).ranks); + const size_t bits = routing::linear_bits_for(ranks); + return (ranks >> bits) == 1 ? static_cast(bits) : 0; +} + // Derives both sides at once: the count matrix is symmetric, so the recv layout is the send layout. auto derive_layer_exchange(const LayerTraversal &layer, const mpi::Comm &comm, int scale, LayerExchangeLayout &layout) -> void { @@ -115,7 +129,8 @@ inline auto begin_flat_exchange(FlatExchangeBuffers &buffers, const mpi::Comm &c .recv_counts = layout.counts.data(), .recv_displs = layout.displs.data()}, mpi::size(comm), - comm); + comm, + layer_exchange_wire_bits(comm)); return handle; } diff --git a/cpp/monoprop/detail/EnvConfig.h b/cpp/monoprop/detail/EnvConfig.h index 30120417..81b0786a 100644 --- a/cpp/monoprop/detail/EnvConfig.h +++ b/cpp/monoprop/detail/EnvConfig.h @@ -14,18 +14,35 @@ #pragma once +#include #include +#include #include #include +#include +#include +#include // Single home for runtime environment configuration. Kept dependency-free by design, because it is // pulled into hot-path headers. // // monoprop_NUM_THREADS positive int (1..1e6), else ignored → num_threads // monoprop_PARTITIONS int N | "auto" | "off"; parsed where it is used (resolve_partition_count_) +// monoprop_ROUTING "splitmix" | "linear" → routing_mode +// monoprop_ROUTE_SEED decimal uint64 basis seed → route_seed +// +// Both routing knobs THROW on a malformed value instead of falling back: each silently changes the +// transport, so a typo that defaulted would stay invisible until a performance postmortem. namespace monoprop::config { +class EnvConfigError : public std::runtime_error { +public: + using std::runtime_error::runtime_error; +}; + +enum class RoutingMode : std::uint8_t { Splitmix, Linear }; + namespace detail { inline auto parse_positive_int(const char *text) -> std::optional { @@ -43,10 +60,49 @@ inline auto parse_positive_int(const char *text) -> std::optional { return static_cast(value); } +[[noreturn]] inline auto reject_env(std::string_view name, const char *text, std::string_view expected) -> void { + throw EnvConfigError(std::string{name} + "=\"" + text + "\" is not " + std::string{expected} + + "; correct it or leave the variable unset."); +} + +// Unset and empty are both nullopt; a value that is present but unparseable throws. +inline auto parse_uint64(std::string_view name, const char *text) -> std::optional { + if (text == nullptr || *text == '\0') { + return std::nullopt; + } + // strtoull WRAPS a negative literal to a huge unsigned rather than failing, so '-' is rejected here. + if (std::string_view{text}.find('-') != std::string_view::npos) { + reject_env(name, text, "a decimal uint64"); + } + errno = 0; + char *end = nullptr; + const unsigned long long value = std::strtoull(text, &end, 10); + if (end == text || *end != '\0' || errno == ERANGE) { + reject_env(name, text, "a decimal uint64"); + } + return static_cast(value); +} + +inline auto parse_routing_mode(std::string_view name, const char *text) -> std::optional { + if (text == nullptr || *text == '\0') { + return std::nullopt; + } + const std::string_view value{text}; + if (value == "splitmix") { + return RoutingMode::Splitmix; + } + if (value == "linear") { + return RoutingMode::Linear; + } + reject_env(name, text, "one of \"splitmix\" or \"linear\""); +} + } // namespace detail struct Settings { std::optional num_threads; + std::optional routing_mode; + std::optional route_seed; }; // Parse the environment once; the Settings are cached and shared across TUs. @@ -54,6 +110,8 @@ inline auto get() -> const Settings & { static const Settings settings = [] { Settings s; s.num_threads = detail::parse_positive_int(std::getenv("monoprop_NUM_THREADS")); + s.routing_mode = detail::parse_routing_mode("monoprop_ROUTING", std::getenv("monoprop_ROUTING")); + s.route_seed = detail::parse_uint64("monoprop_ROUTE_SEED", std::getenv("monoprop_ROUTE_SEED")); return s; }(); return settings; diff --git a/cpp/monoprop/detail/evolution/layer_build/Engine.h b/cpp/monoprop/detail/evolution/layer_build/Engine.h index 0f03f4af..3c93ecfa 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Engine.h +++ b/cpp/monoprop/detail/evolution/layer_build/Engine.h @@ -84,7 +84,8 @@ struct GraphSink { std::vector acc; size_t def_in_base_ = 0; // deferred self-miss bases into acc[my_rank] size_t def_out_base_ = 0; - std::vector in_base_; // cross-rank per-rank base into acc[s].in_entries (set in prepare) + // Cross-rank base into acc[slot].in_entries, over the query window (set in prepare). + mpi::WindowVec in_base_; GraphSink(size_t R_, size_t my_rank_) : R(R_), my_rank(my_rank_), acc(R_) {} @@ -105,33 +106,35 @@ struct GraphSink { // Cross-rank (R>1). Send buffer = the plain query stream (no value fusion). The exchange is positional: // responses[s][q] must answer incoming[s][q], one resolution per query. - auto send_buffer(std::vector &queries, - std::vector> & /*vals*/, - std::vector & /*scratch*/) -> std::vector & { + auto send_buffer(mpi::WindowVec &queries, + mpi::WindowVec> & /*vals*/, + mpi::WindowVec & /*scratch*/) -> mpi::WindowVec & { return queries; } + // `acc` stays flat [P] -- finalize hands it to build_layer_storage_unified, which is P-shaped -- so + // the window index is turned back into a slot here rather than re-basing it. auto prepare(const IncomingProbe & /*pr*/, - size_t rank_count, MPOperator & /*op*/, - const std::vector> &responses) -> void { - in_base_.assign(rank_count, 0); - for (size_t s = 0; s < rank_count; ++s) { - in_base_[s] = acc[s].in_entries.size(); - acc[s].in_entries.resize(in_base_[s] + responses[s].size()); + const mpi::WindowVec> &responses) -> void { + const mpi::SlotWindow w = responses.window(); + in_base_.reset(w); + for (size_t k = 0; k < w.count; ++k) { + const mpi::WindowIndex wi{k}; + PartnerAcc &a = acc[w.slot(wi)]; + in_base_[wi] = a.in_entries.size(); + a.in_entries.resize(in_base_[wi] + responses[wi].size()); } } auto on_resolved(size_t g, - size_t s, + mpi::WindowIndex s, size_t q, size_t ip, const IncomingProbe &pr, - const std::vector & /*incoming*/) -> Response { - acc[s].in_entries[in_base_[s] + q] = {ip, pr.phase_of[g]}; + const mpi::WindowVec & /*incoming*/) -> Response { + acc[pr.window.slot(s)].in_entries[in_base_[s] + q] = {ip, pr.phase_of[g]}; return static_cast(ip); } - auto process_reserve(const std::vector> & /*inc_r*/, - size_t /*rank_count*/, - size_t /*my_rank*/) -> void {} + auto process_reserve(const mpi::WindowVec> & /*inc_r*/, size_t /*my_rank*/) -> void {} auto on_response_block(size_t r, const std::vector &resp, const std::vector &srcs, @@ -236,28 +239,30 @@ struct ContractSink { // Cross-rank (R>1). Send buffer = queries interleaved with their v_src stream into `scratch` // (combined_qv_), so one alltoallv carries query + value. - auto send_buffer(std::vector &queries, std::vector> &vals, std::vector &scratch) - -> std::vector & { - scratch.resize(queries.size()); - for (size_t r = 0; r < queries.size(); ++r) { - QueryCodec::build_fused(queries[r], vals[r], scratch[r]); + auto send_buffer(mpi::WindowVec &queries, + mpi::WindowVec> &vals, + mpi::WindowVec &scratch) -> mpi::WindowVec & { + const mpi::SlotWindow w = queries.window(); + scratch.reset(w); + for (size_t k = 0; k < w.count; ++k) { + const mpi::WindowIndex wi{k}; + QueryCodec::build_fused(queries[wi], vals[wi], scratch[wi]); } return scratch; } auto prepare(const IncomingProbe &pr, - size_t /*rank_count*/, MPOperator &op, - const std::vector> & /*responses*/) -> void { + const mpi::WindowVec> & /*responses*/) -> void { state_mask_ = schrodinger ? initial_state_mask(op.initial_state) : Monomial{}; cross_base_ = fc.cross_half.size(); fc.cross_half.resize(cross_base_ + pr.nq_total); } auto on_resolved(size_t g, - size_t s, + mpi::WindowIndex s, size_t /*q*/, size_t ip, const IncomingProbe &pr, - const std::vector &incoming) -> Response { + const mpi::WindowVec &incoming) -> Response { double v_tgt; if (ip < pr.base) { v_tgt = fused_scale ? op_coeffs[ip] * inv_cos : op_coeffs[ip]; @@ -278,11 +283,13 @@ struct ContractSink { /*is_insert=*/ip >= pr.base}; return v_tgt; } - auto process_reserve(const std::vector> &inc_r, size_t rank_count, size_t my_rank_) -> void { + auto process_reserve(const mpi::WindowVec> &inc_r, size_t my_rank_) -> void { + const mpi::SlotWindow w = inc_r.window(); size_t incoming = 0; - for (size_t r = 0; r < rank_count; ++r) { - if (r != my_rank_) { - incoming += inc_r[r].size(); + for (size_t k = 0; k < w.count; ++k) { + const mpi::WindowIndex wi{k}; + if (w.slot(wi) != my_rank_) { + incoming += inc_r[wi].size(); } } fc.cross_half.reserve(fc.cross_half.size() + incoming); @@ -338,18 +345,24 @@ struct LayerBuildEngine { // → distinct found, so each slot is marked once. MatchedEpochSet &matched; size_t combined_size; - std::vector queries_r; - std::vector> src_idx_r; + // The destination slots this generator can reach: `plan`'s window for my_rank. Every per-slot array + // below is sized to it, so a flat slot only ever enters through WindowVec::at_slot. + mpi::SlotWindow window; + mpi::WindowVec queries_r; + mpi::WindowVec> src_idx_r; std::vector deferred_self_misses; // Deferred-miss positions, concatenated in miss order; parallel to deferred_self_misses. std::vector deferred_pos_flat_; // This pass's self-owned queries as positions, straight from the scan: never encoded, so the resolve - // below has nothing to decode. Parallel to src_idx_r[my_rank]. + // below has nothing to decode. Parallel to src_idx_r's self slot. SelfQueryStage self_stage_; // Scan-captured v_src per query (ContractSink only via Sink::wants_values; empty for GraphSink). - std::vector> src_val_r; + mpi::WindowVec> src_val_r; // Fused query+value send scratch (ContractSink, R>1): shared by a gate's two exchange passes. - std::vector combined_qv_; + mpi::WindowVec combined_qv_; + // Which destination ranks this gate's queries can reach. Dense unless the router is GF(2)-linear; + // see mpi::PeerPlan. Derived once per layer in build_layer, never per query. + mpi::PeerPlan plan; Sink sink; LayerBuildEngine(MPOperator &local_op_, @@ -358,34 +371,46 @@ struct LayerBuildEngine { size_t my_rank_, MatchedEpochSet &matched_scratch, size_t combined_size_, - Sink &&sink_) + Sink &&sink_, + mpi::PeerPlan plan_ = {}) // dense by default: the tests build the engine directly : local_op(local_op_), comm(comm_), R(R_), my_rank(my_rank_), matched(matched_scratch), combined_size(combined_size_), - queries_r(R_), - src_idx_r(R_), + plan(plan_), sink(std::move(sink_)) { + const auto geom = mpi::geometry(comm); + window = plan.window(my_rank, static_cast(geom.ranks), static_cast(geom.partitions)); + assert(window.stop() <= R && window.count != 0); + queries_r.reset(window); + src_idx_r.reset(window); matched.begin_gate(combined_size); } // Resolve this rank's own query stream inline, then clear it so the alltoallv never sends to self. + // Self is inside the window only when this generator's rank shift is zero; otherwise the window names + // another rank outright and the scan cannot have staged a self-owned partner. auto resolve_self_queries(bool is_leader_pass) -> void { - std::vector &ls = src_idx_r[my_rank]; + if (!window.contains(my_rank)) { + assert(self_stage_.size() == 0 && "a self-owned partner outside this generator's peer window"); + self_stage_.clear(); + return; + } + std::vector &ls = src_idx_r.at_slot(my_rank); std::vector *lv = nullptr; if constexpr (Sink::wants_values) { - lv = &src_val_r[my_rank]; + lv = &src_val_r.at_slot(my_rank); } // The scan routes a self-owned partner to the stage, never to the wire buffer. - assert(queries_r[my_rank].empty() && "a self-owned query was encoded instead of staged"); + assert(queries_r.at_slot(my_rank).empty() && "a self-owned query was encoded instead of staged"); assert(ls.size() == self_stage_.size() && "the self stage does not hold exactly one query per source"); resolve_range_(ls, lv, is_leader_pass); self_stage_.clear(); ls.clear(); if constexpr (Sink::wants_values) { - src_val_r[my_rank].clear(); + src_val_r.at_slot(my_rank).clear(); } } @@ -396,10 +421,14 @@ struct LayerBuildEngine { // pass must also drop the queries a leader already matched, and that only holds once the leader pass // has run. auto run_exchange(bool is_leader_pass, - std::vector &&queries, - std::vector> &&src_idx, - std::vector> &&src_val, + mpi::WindowVec &&queries, + mpi::WindowVec> &&src_idx, + mpi::WindowVec> &&src_val, SelfQueryStage &&self_stage) -> void { + // The scan sized its arrays to the same plan, so the two windows must agree exactly -- a mismatch + // would re-base every slot against the wrong base. + assert(queries.window().base == window.base && queries.window().count == window.count); + assert(src_idx.window().base == window.base && src_idx.window().count == window.count); queries_r = std::move(queries); src_idx_r = std::move(src_idx); self_stage_ = std::move(self_stage); @@ -412,30 +441,32 @@ struct LayerBuildEngine { if (R <= 1) { return; } - std::vector &send = sink.send_buffer(queries_r, src_val_r, combined_qv_); - std::vector> inc_q; - mpi::begin_alltoallv(send, comm).wait_into(inc_q); - auto resp = resolve_incoming(inc_q, local_op, R, is_leader_pass, matched, combined_size, sink); + mpi::WindowVec &send = sink.send_buffer(queries_r, src_val_r, combined_qv_); + mpi::WindowVec inc_q; + mpi::begin_alltoallv(send, comm, /*skip_self=*/false, /*known_recv_counts=*/nullptr, plan).wait_into(inc_q); + auto resp = resolve_incoming(inc_q, local_op, is_leader_pass, matched, combined_size, sink); std::vector resp_recv = response_recv_counts(); - std::vector> inc_r; - mpi::begin_alltoallv(resp, comm, /*skip_self=*/false, &resp_recv).wait_into(inc_r); - process_responses(inc_r, src_idx_r, queries_r, R, my_rank, sink); + mpi::WindowVec> inc_r; + // The answers retrace the queries, and the pairing is an XOR involution, so the same plan holds. + mpi::begin_alltoallv(resp, comm, /*skip_self=*/false, &resp_recv, plan).wait_into(inc_r); + process_responses(inc_r, src_idx_r, queries_r, my_rank, sink); } // Followers a leader already matched must not be re-resolved over the wire, so compact them out. auto drop_matched_cross_rank_followers() -> void { using QC = QueryCodec; const QueryLayout layout = sink.querier_layout(); - for (size_t r = 0; r < R; ++r) { - if (r == my_rank) { + for (size_t k = 0; k < window.count; ++k) { + const mpi::WindowIndex wi{k}; + if (window.slot(wi) == my_rank) { continue; } - VecZ &q = queries_r[r]; - std::vector &s = src_idx_r[r]; + VecZ &q = queries_r[wi]; + std::vector &s = src_idx_r[wi]; // Fused: the v_src stream is parallel to the query/source streams, so compact it in lockstep. std::vector *v = nullptr; if constexpr (Sink::wants_values) { - v = &src_val_r[r]; + v = &src_val_r[wi]; } const size_t nq = s.size(); size_t kept = 0; @@ -495,13 +526,16 @@ struct LayerBuildEngine { private: // Response counts are the transpose of the query counts (one answer per query), so passing them as // known_recv_counts skips the response count-Alltoall round. + // FLAT [P], which is what begin_alltoallv's known_recv_counts is indexed by; only the window's slots + // can be non-zero, and the window is what masks the rest. auto response_recv_counts() const -> std::vector { - std::vector counts(R); - for (size_t r = 0; r < R; ++r) { - // One response per QUERY, and src_idx_r[r] holds one source per query: no walk, no division. - assert(src_idx_r[r].size() == QueryCodec::count_queries(queries_r[r], sink.querier_layout()) + std::vector counts(R, 0); + for (size_t k = 0; k < window.count; ++k) { + const mpi::WindowIndex wi{k}; + // One response per QUERY, and src_idx_r's block holds one source per query: no walk, no division. + assert(src_idx_r[wi].size() == QueryCodec::count_queries(queries_r[wi], sink.querier_layout()) && "a querier buffer does not hold exactly one query per source"); - counts[r] = static_cast(src_idx_r[r].size()); + counts[window.slot(wi)] = static_cast(src_idx_r[wi].size()); } return counts; } @@ -600,6 +634,17 @@ auto build_layer(MPOperator &local_op, validate_only_rotate_len_k_(only_rotate_len_k, 2 * NumModes); const size_t my_rank = static_cast(mpi::rank(comm)); const size_t R = static_cast(mpi::size(comm)); + // R is the FLAT world (ranks x partitions); the router is what splits it back into the two levels. + const routing::Router router = router_for(comm); + assert(router.flat_world() == R); + // Under linear routing every query for THIS generator lands on the rank this rank's own index XOR + // rank_shift(gen), so the exchange knows its peer before it starts. Dense otherwise, which is + // today's collective. + const size_t gen_shift = router.rank_shift(gen); + const auto plan = mpi::PeerPlan{.sparse = router.is_linear(), .shift = static_cast(gen_shift)}; + // The reachable slots, once per generator: S of the P=R*S world under linear routing, all P otherwise. + // Every per-slot structure from the scan to the wire is sized to this run. + const mpi::SlotWindow scan_window = plan.window(my_rank, router.ranks(), router.partitions()); // Fused contraction runs at all rank counts (R>1 via the cross-rank half-rotation exchange). const bool use_fused = (fused_contract != nullptr); const auto cut_st = build_majorana_evolution_cutoff_state(atol, local_coeffs, upper_atol, param); @@ -619,36 +664,45 @@ auto build_layer(MPOperator &local_op, } assert(fused_scale_coeffs == nullptr || (local_coeffs && &local_coeffs->get() == fused_scale_coeffs)); - FusedScanResult fused = [&] { + // An identity generator anticommutes with nothing: the scan returns on its empty fold-column set + // with no query, no cosine block and no coefficient swept, and run_exchange's three collectives per + // pass would carry no payload. The generator list is replicated, so skipping needs no agreement. + // (A zero chemical potential alone contributes 60 of the 60-site Hubbard's 476 generators per + // Trotter layer.) No gate is merged: a no-op gate is simply not exchanged for. + const bool identity_gen = !gen.any(); + + FusedScanResult fused; + CosMask cos_all; + if (!identity_gen) { double *const sweep_ptr = fused_scale ? fused_scale_coeffs->data() : nullptr; - return with_algebra(basis, [&]() { + fused = with_algebra(basis, [&]() { return fused_find_and_collect(local_op, gen, cut_eval, cut_st, coeffs, only_rotate_len_k, - R, + scan_window, my_rank, + router, + gen_shift, /*capture_values=*/use_fused, sweep_ptr, cos_build); }); - }(); - - CosMask cos_all; - if (fused.cos_blocks.size() == 1) { - // The serial scan produces a single cosine block set — take it wholesale. - cos_all = std::move(fused.cos_blocks[0]); - } - else { - // Cosine block sets are disjoint and ascending; concatenate in order. - for (const auto &block : fused.cos_blocks) { - cos_all.total_count += block.total_count; - cos_all.blocks.insert(cos_all.blocks.end(), block.blocks.begin(), block.blocks.end()); + if (fused.cos_blocks.size() == 1) { + // The serial scan produces a single cosine block set — take it wholesale. + cos_all = std::move(fused.cos_blocks[0]); } + else { + // Cosine block sets are disjoint and ascending; concatenate in order. + for (const auto &block : fused.cos_blocks) { + cos_all.total_count += block.total_count; + cos_all.blocks.insert(cos_all.blocks.end(), block.blocks.begin(), block.blocks.end()); + } + } + fused.cos_blocks = std::vector{}; } - fused.cos_blocks = std::vector{}; auto run = [&](Sink sink) -> std::shared_ptr { LayerBuildEngine eng(local_op, @@ -657,17 +711,20 @@ auto build_layer(MPOperator &local_op, my_rank, matched_scratch, /*combined_size=*/local_op.store->size(), - std::move(sink)); - eng.run_exchange(/*is_leader_pass=*/true, - std::move(fused.leader_queries), - std::move(fused.leader_src), - std::move(fused.leader_val), - std::move(fused.leader_self)); - eng.run_exchange(/*is_leader_pass=*/false, - std::move(fused.follower_queries), - std::move(fused.follower_src), - std::move(fused.follower_val), - std::move(fused.follower_self)); + std::move(sink), + plan); + if (!identity_gen) { + eng.run_exchange(/*is_leader_pass=*/true, + std::move(fused.leader_queries), + std::move(fused.leader_src), + std::move(fused.leader_val), + std::move(fused.leader_self)); + eng.run_exchange(/*is_leader_pass=*/false, + std::move(fused.follower_queries), + std::move(fused.follower_src), + std::move(fused.follower_val), + std::move(fused.follower_self)); + } return eng.finish(std::move(cos_all), out_cos); }; diff --git a/cpp/monoprop/detail/evolution/layer_build/Resolve.h b/cpp/monoprop/detail/evolution/layer_build/Resolve.h index 0c49706f..7e6d5902 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Resolve.h +++ b/cpp/monoprop/detail/evolution/layer_build/Resolve.h @@ -24,6 +24,7 @@ #include "monoprop/detail/evolution/CutoffContext.h" #include "monoprop/detail/evolution/layer_build/Common.h" #include "monoprop/detail/evolution/layer_build/QueryCodec.h" +#include "monoprop/detail/mpi/Comm.h" #include "monoprop/detail/operator/MPOperator.h" #include "monoprop/detail/operator/RowAccess.h" @@ -38,10 +39,12 @@ struct IncomingProbe { // The STORE's position width, not the wire's: these positions exist to become rows. using PosT = typename OperatorIndex::PosT; - std::vector goff; // rank_count+1 flat offsets: g = goff[s] + q - DefaultInitVector sender_of; // g → sender rank + // The slots `incoming` covers; senders are named by their index into it, never by a flat slot. + mpi::SlotWindow window; + std::vector goff; // window.count+1 flat offsets: g = goff[k] + q + DefaultInitVector sender_wi; // g → sender's WINDOW index (see sender_index/sender_slot) DefaultInitVector phase_of; // g → query phase - // g → WORD offset of that query inside incoming[sender_of[g]]; a query ordinal names no position. + // g → WORD offset of that query inside its sender's buffer; a query ordinal names no position. DefaultInitVector off_of; DefaultInitVector idx_of; // g → resolved index (hit: < base; miss: base+j) std::vector miss_g; // j → the g that became miss j (Phase 4 reads the key of miss_g[j]) @@ -55,6 +58,10 @@ struct IncomingProbe { // g → fold_hash of the query key, folded by the probe and reused by the insert. DefaultInitVector hash_of; + // The two ways to name query g's sender. sender_wi is re-based, so nothing else reads it. + [[nodiscard]] auto sender_index(size_t g) const -> mpi::WindowIndex { return mpi::WindowIndex{sender_wi[g]}; } + [[nodiscard]] auto sender_slot(size_t g) const -> size_t { return window.slot(sender_index(g)); } + // BUILDS a bitset, so cold consumers only -- the fully paired minority, never anything per-term. [[nodiscard]] auto mono_at(size_t g) const -> Monomial { Monomial m; @@ -77,28 +84,29 @@ struct IncomingProbe { // fused for the ContractSink resolver, plain for GraphSink. The caller runs Phase 3, then // insert_incoming_misses. Counts and offsets come from the decode walk; there is no record stride. template -auto probe_incoming_queries(const std::vector &incoming, // serialized, one VecZ per sender +auto probe_incoming_queries(const mpi::WindowVec &incoming, // serialized, one VecZ per sender slot MPOperator &op, - size_t rank_count, QueryLayout layout) -> IncomingProbe { using QC = QueryCodec; IncomingProbe pr; + pr.window = incoming.window(); + const size_t senders = pr.window.count; - pr.goff.assign(rank_count + 1, 0); - for (size_t s = 0; s < rank_count; ++s) { - const size_t nq = QC::count_queries(incoming[s], layout); - pr.goff[s + 1] = pr.goff[s] + nq; + pr.goff.assign(senders + 1, 0); + for (size_t k = 0; k < senders; ++k) { + const size_t nq = QC::count_queries(incoming[mpi::WindowIndex{k}], layout); + pr.goff[k + 1] = pr.goff[k] + nq; } - pr.nq_total = pr.goff[rank_count]; + pr.nq_total = pr.goff[senders]; if (pr.nq_total == 0) { return pr; } - pr.sender_of.resize(pr.nq_total); - for (size_t s = 0; s < rank_count; ++s) { - std::fill(pr.sender_of.begin() + static_cast(pr.goff[s]), - pr.sender_of.begin() + static_cast(pr.goff[s + 1]), - static_cast(s)); + pr.sender_wi.resize(pr.nq_total); + for (size_t k = 0; k < senders; ++k) { + std::fill(pr.sender_wi.begin() + static_cast(pr.goff[k]), + pr.sender_wi.begin() + static_cast(pr.goff[k + 1]), + static_cast(k)); } // Phase 1 (read-only): deserialize, then probe with the group-prefetch batch find. One walk per sender. @@ -111,20 +119,21 @@ auto probe_incoming_queries(const std::vector &incoming, // serialized, on pr.pos_flat.clear(); // A hint only: the measured mean is 5.33 positions, so this is one allocation but for an outlier. pr.pos_flat.reserve(pr.nq_total * QueryCodec::kReservePositionsPerQuery); - for (size_t s = 0; s < rank_count; ++s) { + for (size_t si = 0; si < senders; ++si) { + const VecZ &buf = incoming[mpi::WindowIndex{si}]; size_t off = 0; - for (size_t g = pr.goff[s]; g < pr.goff[s + 1]; ++g) { + for (size_t g = pr.goff[si]; g < pr.goff[si + 1]; ++g) { int ph = 0; - const size_t k = QC::k_at(incoming[s], off); + const size_t k = QC::k_at(buf, off); const size_t at = pr.pos_flat.size(); pr.pos_flat.resize(at + k); // default-init grow: read_positions writes every element pr.pos_off[g] = at; pr.k_of[g] = static_cast(k); pr.off_of[g] = off; - off = QC::read_positions(incoming[s], layout, off, pr.pos_flat.data() + at, ph); + off = QC::read_positions(buf, layout, off, pr.pos_flat.data() + at, ph); pr.phase_of[g] = ph; } - assert(off == incoming[s].size() && "the query walk did not consume the sender's whole buffer"); + assert(off == buf.size() && "the query walk did not consume the sender's whole buffer"); } { const size_t op_size = op.store->size(); @@ -179,19 +188,19 @@ auto insert_incoming_misses(MPOperator &op, const IncomingProbe -auto resolve_incoming(const std::vector &incoming, // serialized, one VecZ per sender +auto resolve_incoming(const mpi::WindowVec &incoming, // serialized, one VecZ per sender slot MPOperator &op, - size_t rank_count, bool is_leader_pass, MatchedEpochSet &matched, size_t combined_size, // pre-layer op size: bounds the matched set - Sink &sink) -> std::vector> { + Sink &sink) -> mpi::WindowVec> { using Resp = typename Sink::Response; - const IncomingProbe pr = - probe_incoming_queries(incoming, op, rank_count, sink.incoming_layout()); - std::vector> responses(rank_count); - for (size_t s = 0; s < rank_count; ++s) { - responses[s].assign(pr.goff[s + 1] - pr.goff[s], Sink::init_response()); + const IncomingProbe pr = probe_incoming_queries(incoming, op, sink.incoming_layout()); + // The response window is the query window: the pairing is an XOR involution, so a rank answers + // exactly the slots it queried. + mpi::WindowVec> responses(pr.window); + for (size_t k = 0; k < pr.window.count; ++k) { + responses[mpi::WindowIndex{k}].assign(pr.goff[k + 1] - pr.goff[k], Sink::init_response()); } if (pr.nq_total == 0) { return responses; @@ -199,10 +208,10 @@ auto resolve_incoming(const std::vector &incoming, // serialized, one VecZ // Phase 3 (scatter): responses + sink records + matched-follower marks. Freshly inserted partners // (ip ≥ combined_size) skip the mark. - sink.prepare(pr, rank_count, op, responses); + sink.prepare(pr, op, responses); for (size_t g = 0; g < pr.nq_total; ++g) { - const size_t s = pr.sender_of[g]; - const size_t q = g - pr.goff[s]; + const mpi::WindowIndex s = pr.sender_index(g); + const size_t q = g - pr.goff[s.value]; const size_t ip = pr.idx_of[g]; responses[s][q] = sink.on_resolved(g, s, q, ip, pr, incoming); if (is_leader_pass && ip < combined_size) { @@ -215,20 +224,25 @@ auto resolve_incoming(const std::vector &incoming, // serialized, one VecZ } // Querier rank (any cross-rank sink): fold each resolver response into a querier-side record. The self/ -// local rank was already resolved inline, so it is skipped here. inc_r[r][q] answers query q from rank r. +// local slot was already resolved inline, so it is skipped here (and is in the window only when this +// generator's rank shift is zero). inc_r[k][q] answers query q sent to the window's k-th slot. template -auto process_responses(const std::vector> &inc_r, - const std::vector> &src_idx, - const std::vector &queries, // serialized query buffers (for phase recovery) - size_t rank_count, +auto process_responses(const mpi::WindowVec> &inc_r, + const mpi::WindowVec> &src_idx, + const mpi::WindowVec &queries, // serialized query buffers (for phase recovery) size_t my_rank, Sink &sink) -> void { - sink.process_reserve(inc_r, rank_count, my_rank); - for (size_t r = 0; r < rank_count; ++r) { + const mpi::SlotWindow w = inc_r.window(); + assert(src_idx.window().base == w.base && src_idx.window().count == w.count); + assert(queries.window().base == w.base && queries.window().count == w.count); + sink.process_reserve(inc_r, my_rank); + for (size_t k = 0; k < w.count; ++k) { + const mpi::WindowIndex wi{k}; + const size_t r = w.slot(wi); if (r == my_rank) { continue; } - sink.on_response_block(r, inc_r[r], src_idx[r], queries[r]); + sink.on_response_block(r, inc_r[wi], src_idx[wi], queries[wi]); } } diff --git a/cpp/monoprop/detail/evolution/layer_build/Scan.h b/cpp/monoprop/detail/evolution/layer_build/Scan.h index cce972b0..2cd41f95 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Scan.h +++ b/cpp/monoprop/detail/evolution/layer_build/Scan.h @@ -33,6 +33,7 @@ #include "monoprop/detail/evolution/layer_build/PartnerMerge.h" #include "monoprop/detail/evolution/layer_build/QueryCodec.h" #include "monoprop/detail/graph_encoding/MPGraphEncodingTypes.h" +#include "monoprop/detail/mpi/Comm.h" #include "monoprop/detail/mpi/MPIUtils.h" #include "monoprop/detail/operator/InvertedIndex.h" #include "monoprop/detail/operator/MPOperator.h" @@ -220,27 +221,37 @@ template template struct FusedScanResult { - std::vector cos_blocks; // ascending, disjoint, chunk order - std::vector leader_queries; // size R: serialized leader queries per owner rank - std::vector> leader_src; // size R: parallel to leader_queries (source op idx) - std::vector follower_queries; // size R: serialized follower queries per owner rank - std::vector> follower_src; // size R: parallel to follower_queries + std::vector cos_blocks; // ascending, disjoint, chunk order + // The six arrays below are indexed by DESTINATION SLOT through WindowVec::at_slot, and cover only + // the slots this generator can reach: S of the P=R*S world under linear routing, all P under + // splitmix. See mpi::PeerPlan::window. + mpi::SlotWindow window; + mpi::WindowVec leader_queries; // serialized leader queries per owner slot + mpi::WindowVec> leader_src; // parallel to leader_queries (source op idx) + mpi::WindowVec follower_queries; // serialized follower queries per owner slot + mpi::WindowVec> follower_src; // parallel to follower_queries // Fused-contraction only (capture_values): signed pre-cos source coeff (v_src) parallel to // leader_src / follower_src. Empty when capture_values is false. - std::vector> leader_val; - std::vector> follower_val; - // Self-owned queries, staged as positions instead of encoded into leader_queries[my_rank]: they are - // resolved inline and never reach a wire, so the codec is not on this leg. Order matches - // leader_src[my_rank] / follower_src[my_rank], which is the accumulation order. + mpi::WindowVec> leader_val; + mpi::WindowVec> follower_val; + // Self-owned queries, staged as positions instead of encoded into the window's self slot: they are + // resolved inline and never reach a wire, so the codec is not on this leg. Order matches that slot's + // leader_src / follower_src, which is the accumulation order. Empty unless the window contains + // my_rank, i.e. unless this generator's rank shift is zero. SelfQueryStage leader_self; SelfQueryStage follower_self; }; // Classify, cut off and emit in one pass over the anticommuting terms. Queries go to the owner of -// M'=M⊕G (hash%R; self at R==1) in ascending source-index order, so resolve and index assignment are +// M'=M⊕G (routing::Router; self at R==1) in ascending source-index order, so resolve and index assignment are // deterministic. `fused_scale_coeffs` (no length cap only; must alias coeffs.data()) scales every anticommuting // coeff in place by `fused_scale_cos`=cos(2·build_angle), so no cosine set is built and a hit's stored // value is post-cos (resolve recovers it via 1/cos). +// +// `gen_shift` is router.rank_shift(gen), and `op` must hold only terms `my_rank` owns -- then the owner of +// M⊕G is rank(M) ^ gen_shift and the linear planes never run per term. Both are what mpi::PeerPlan +// already assumes; a violation moves ownership silently, so the fast path asserts against dest(). +// `window` must be that plan's window for `my_rank`: it is what the six query arrays are sized to. template auto fused_find_and_collect(const MPOperator &op, const Monomial &gen, @@ -248,25 +259,30 @@ auto fused_find_and_collect(const MPOperator &op, const CutoffContext &cut_st, const VecD &coeffs, std::optional only_rotate_len_k, - size_t rank_count, + mpi::SlotWindow window, size_t my_rank, + const routing::Router &router, + size_t gen_shift, bool capture_values = false, double *fused_scale_coeffs = nullptr, double fused_scale_cos = 1.0) -> FusedScanResult { validate_only_rotate_len_k_(only_rotate_len_k, 2 * NumModes); const size_t gen_pop = gen.count(); + const size_t rank_count = router.flat_world(); const auto ectx = A::make_gen_context(gen); + assert(window.stop() <= rank_count && window.count != 0); FusedScanResult res; - res.leader_queries.assign(rank_count, VecZ{}); - res.leader_src.assign(rank_count, std::vector{}); - res.follower_queries.assign(rank_count, VecZ{}); - res.follower_src.assign(rank_count, std::vector{}); - // Sized to R even on the early-return paths below so the fused engine's per-rank src_val_r access - // is always in bounds (parallel to leader_src / follower_src). + res.window = window; + res.leader_queries.reset(window); + res.leader_src.reset(window); + res.follower_queries.reset(window); + res.follower_src.reset(window); + // Sized on the early-return paths below too, so the fused engine's per-slot src_val_r access is + // always in bounds (parallel to leader_src / follower_src). if (capture_values) { - res.leader_val.assign(rank_count, std::vector{}); - res.follower_val.assign(rank_count, std::vector{}); + res.leader_val.reset(window); + res.follower_val.reset(window); } { @@ -343,21 +359,28 @@ auto fused_find_and_collect(const MPOperator &op, double v_src, bool is_follower) { // Single rank: every partner is self-owned, skip the O(W) hash; multi-rank routes by owner. - // Must be the SAME function find_rank computes (MPIUtils.h) or a term is placed and queried - // on different ranks, which duplicates a row silently; mpi_utils_tests.cpp asserts it. + // routing::Router is the ONLY owner function -- find_rank (MPIUtils.h) must stay in step + // with it, or a term is placed and queried on different ranks, which duplicates a row + // silently; mpi_utils_tests.cpp asserts the agreement. size_t r_prime = my_rank; if (rank_count != 1) { - r_prime = monomial_hash(dense) % rank_count; + r_prime = router.dest_from_shift(dense, my_rank, gen_shift); + assert(r_prime == router.dest(dense)); // an identity, not an approximation } + // at_slot is the only re-basing door and asserts membership: a destination outside this + // generator's window means the shift is wrong, and would otherwise land on another peer. if (r_prime == my_rank) { (is_follower ? res.follower_self : res.leader_self).push(pos, k, phase); } else { - QueryCodec::push_positions(is_follower ? fq[r_prime] : lq[r_prime], pos, k, phase); + QueryCodec::push_positions(is_follower ? fq.at_slot(r_prime) : lq.at_slot(r_prime), + pos, + k, + phase); } - (is_follower ? fs[r_prime] : ls[r_prime]).push_back(i); + (is_follower ? fs : ls).at_slot(r_prime).push_back(i); if (capture_values) { - (is_follower ? fv[r_prime] : lv[r_prime]).push_back(v_src); + (is_follower ? fv : lv).at_slot(r_prime).push_back(v_src); } }; @@ -406,9 +429,9 @@ auto fused_find_and_collect(const MPOperator &op, // A hint only, off the measured mean of 5.33 positions; wider terms grow the buffer. const size_t pq = QueryCodec::kReservePositionsPerQuery; res.leader_self.reserve(n_anti - n_foll, pq); - ls[my_rank].reserve(n_anti - n_foll); + ls.at_slot(my_rank).reserve(n_anti - n_foll); res.follower_self.reserve(n_foll, pq); - fs[my_rank].reserve(n_foll); + fs.at_slot(my_rank).reserve(n_foll); } auto derive_coeff = [&](size_t i) -> std::pair { if (capture_values) { diff --git a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl index b2e0aae1..e415d6ff 100644 --- a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl +++ b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl @@ -143,6 +143,8 @@ MonomialPropagator::MonomialPropagator(const OperatorDict &initial_ope const size_t num_ranks = static_cast(mpi::size(comm_)); const size_t my_rank = static_cast(mpi::rank(comm_)); + check_routing_agreement(comm_); // a disagreement here would hang the first exchange, not corrupt it + const routing::Router router = router_for(comm_); // hoisted: geometry() can hit MPI, so never per term MonomialList local_heisenberg_terms; double core_term = 0.0; @@ -155,7 +157,7 @@ MonomialPropagator::MonomialPropagator(const OperatorDict &initial_ope core_term = encoded_coeff; continue; } - if (my_rank == find_rank(majorana_bitset, num_ranks)) { + if (my_rank == find_rank(majorana_bitset, router)) { mp_op_.init_op_map[majorana_bitset] = encoded_coeff; local_heisenberg_terms.push_back(majorana_bitset); } @@ -177,7 +179,7 @@ MonomialPropagator::MonomialPropagator(const OperatorDict &initial_ope // The initial monomials are distinct, so emplace (insert-if-absent) is an assigning insert here. for (size_t r = 0; r < op.size(); ++r) { const auto &mono = materialize_row(op, r); - if (my_rank == find_rank(mono, num_ranks)) { + if (my_rank == find_rank(mono, router)) { mp_op_.append_term(mono); mp_op_.store->emplace(mono, i++); } @@ -205,6 +207,7 @@ MonomialPropagator::MonomialPropagator(const MonomialPropagator &other upper_atol_(other.upper_atol_), core_term_(other.core_term_), initial_operator_epoch_(other.initial_operator_epoch_), + routing_coverage_reported_(other.routing_coverage_reported_), logical_num_modes_(other.logical_num_modes_), cutoff_type_(other.cutoff_type_), basis_change_(other.basis_change_), @@ -366,7 +369,7 @@ auto MonomialPropagator::apply_initial_operator_(const OperatorDict &o for_each_partition_([&](MonomialPropagator &s) { s.update_initial_operator(op_dict); }); return {}; } - const size_t num_ranks = static_cast(mpi::size(comm_)); + const routing::Router router = router_for(comm_); // hoisted: geometry() can hit MPI, so never per term const size_t my_rank = static_cast(mpi::rank(comm_)); OperatorDict new_op; @@ -376,7 +379,7 @@ auto MonomialPropagator::apply_initial_operator_(const OperatorDict &o core_term_ = algebra_encode_coeff(basis_, coeff, mono); continue; } - if (my_rank == find_rank(mono, num_ranks)) { + if (my_rank == find_rank(mono, router)) { const auto mono_indices = bitset_to_indices(mono); new_op[mono_indices] = coeff; } @@ -718,11 +721,51 @@ auto MonomialPropagator::propagate(const std::vector &majoranas, evolve_mode_contract_immediately_(majoranas, parameter_mapping, gen_coeffs, parameters, only_rotate_len_k); } +template +auto MonomialPropagator::report_routing_coverage_(const std::vector &majoranas) -> void { + // One report per rank, so only its partition 0 speaks, and only once whatever the outcome. + if (routing_coverage_reported_ || comm_.shm_rank != 0) { + return; + } + routing_coverage_reported_ = true; + const routing::Router router = router_for(comm_); + if (!router.is_linear()) { + return; // splitmix: no subspace to fall short of + } + std::vector shifts; + shifts.reserve(majoranas.size()); + for (const auto &gate : majoranas) { + // An out-of-range index is build_evolve_result_'s to reject, gate by gate: converting the whole + // list up front would pre-empt that throw. + if (std::ranges::any_of(gate, [this](size_t i) { return i >= 2 * logical_num_modes_; })) { + return; + } + shifts.push_back(static_cast(router.rank_shift(indices_to_bitset(gate)))); + } + std::ranges::sort(shifts); + shifts.erase(std::ranges::unique(shifts).begin(), shifts.end()); + const size_t span = routing::gf2_rank(shifts); + if (span >= router.linear_bits()) { + return; + } + // A warning, not a throw: every term still lands on one owner, they just do not cover the ranks. + // COMMPLACE's shape -- greppable prefix, rank-identified, one line. + const auto line = std::format("COMMROUTE rank={} linear_bits={} shift_rank={} shifts={} idle_ranks={}\n", + static_cast(mpi::rank(comm_)) / router.partitions(), + router.linear_bits(), + span, + shifts.size(), + router.ranks() - (size_t{1} << span)); + std::fputs(line.c_str(), stderr); + std::fflush(stderr); +} + template template auto MonomialPropagator::run_gate_loop_(const std::vector &majoranas, std::optional only_rotate_len_k, EvolutionFunc evolution_func) -> void { + report_routing_coverage_(majoranas); // Serial per partition; parallelism comes from partitioning the operator across cores. for (size_t i = 0; i < majoranas.size(); ++i) { const auto idx = !schrodinger_ ? majoranas.size() - 1 - i : i; diff --git a/cpp/monoprop/detail/mpi/CMakeLists.txt b/cpp/monoprop/detail/mpi/CMakeLists.txt index 5b9969b0..9128c8ef 100644 --- a/cpp/monoprop/detail/mpi/CMakeLists.txt +++ b/cpp/monoprop/detail/mpi/CMakeLists.txt @@ -11,7 +11,9 @@ target_sources( "HybridComm.h" "MPICompat.h" "MPIUtils.h" + "Pairwise.h" "PartitionBarrier.h" + "Routing.h" "ShmComm.h" ) diff --git a/cpp/monoprop/detail/mpi/Comm.h b/cpp/monoprop/detail/mpi/Comm.h index 9ff3a14c..4b45fff4 100644 --- a/cpp/monoprop/detail/mpi/Comm.h +++ b/cpp/monoprop/detail/mpi/Comm.h @@ -14,6 +14,7 @@ #pragma once +#include #include #include #include @@ -65,6 +66,111 @@ struct Comm { } }; +// A window-relative index. Distinct from a flat slot on purpose: the two are the same number only when +// the window starts at 0, so a swap addresses the wrong peer while staying in bounds. +struct WindowIndex { + size_t value = 0; + + constexpr WindowIndex() = default; + explicit constexpr WindowIndex(size_t v) noexcept : value(v) {} +}; + +// The contiguous run of flat destination slots a round can reach. Slots are rank-major +// (slot = rank * S + partition), so one rank's S partitions are contiguous and the single peer sparse +// routing leaves is exactly one such run; dense is the count == P value of the same run, not a second +// shape. See PeerPlan::window. +struct SlotWindow { + size_t base = 0; // first reachable flat slot + size_t count = 0; // slots in the run + + [[nodiscard]] constexpr auto stop() const -> size_t { return base + count; } + [[nodiscard]] constexpr auto contains(size_t slot) const -> bool { return slot >= base && slot < stop(); } + // The one flat-slot door: it asserts membership, so a slot from outside cannot become another's entry. + [[nodiscard]] constexpr auto index(size_t slot) const -> WindowIndex { + assert(contains(slot) && "flat slot outside the window it is being re-based into"); + return WindowIndex{slot - base}; + } + [[nodiscard]] constexpr auto slot(WindowIndex i) const -> size_t { + assert(i.value < count); + return base + i.value; + } +}; + +// A vector over a SlotWindow, addressed by flat slot through at_slot(); operator[] takes a WindowIndex, +// so a flat slot used as a raw index does not compile. Re-basing an array is only safe if every index +// site shifts together, and these two accessors are the only sites. +template +class WindowVec { +public: + using value_type = T; + + WindowVec() = default; + explicit WindowVec(SlotWindow w) : win_(w), v_(w.count) {} + + auto reset(SlotWindow w) -> void { + win_ = w; + v_.assign(w.count, T{}); + } + + [[nodiscard]] auto window() const -> SlotWindow { return win_; } + [[nodiscard]] auto size() const -> size_t { return v_.size(); } + + [[nodiscard]] auto operator[](WindowIndex i) -> T & { return v_[i.value]; } + [[nodiscard]] auto operator[](WindowIndex i) const -> const T & { return v_[i.value]; } + [[nodiscard]] auto at_slot(size_t slot) -> T & { return v_[win_.index(slot).value]; } + [[nodiscard]] auto at_slot(size_t slot) const -> const T & { return v_[win_.index(slot).value]; } + + [[nodiscard]] auto begin() { return v_.begin(); } + [[nodiscard]] auto end() { return v_.end(); } + [[nodiscard]] auto begin() const { return v_.begin(); } + [[nodiscard]] auto end() const { return v_.end(); } + +private: + SlotWindow win_{}; + std::vector v_; +}; + +// Which destination RANKS a round can touch, when the caller knows. Two states, matching +// routing::Router: dense, or sparse over the single peer GF(2)-linear routing implies. +// +// Sparse means the destination rank of every block is determined by the generator: it is this rank's +// own index XOR `shift`, so +// +// peer = me ^ shift, count == 1 +// +// -- one peer instead of all `ranks`, and the relation is symmetric (XOR is an involution), so every +// rank derives the same pairing with no communication. That is what lets a verb replace a dense +// collective with point-to-point. Linear routing takes ALL log2(ranks) rank bits, so there is no +// intermediate fanout to express here. +// +// Dense is the default: peer(k) == k and count == ranks, so the same loops walk every rank and the +// verbs take their collective path. Every single-rank run is dense (Router::is_linear is false at +// R == 1), so the collectives are not a fallback but the common case. +// +// Two distinct failure modes if `shift` is wrong, which is why the plan is derived in one place. Ranks +// that DISAGREE deadlock: the pairing stops being symmetric and someone waits on a send never posted. +// Ranks that all agree on the same wrong shift stay symmetric and never hang -- they silently DROP the +// blocks outside the peer set, because pack_count_matrix_ / size_staging_send_ / pack_send_ only ever +// touch peers. pack_count_matrix_ asserts the non-peer remainder is empty to catch that one. +struct PeerPlan { + bool sparse = false; + int shift = 0; + + [[nodiscard]] constexpr auto dense() const -> bool { return !sparse; } + [[nodiscard]] constexpr auto count(int ranks) const -> int { return sparse ? 1 : ranks; } + // `k` indexes the peer set, which is a singleton when sparse. + [[nodiscard]] constexpr auto peer(int me, int k) const -> int { return sparse ? (me ^ shift) : k; } + [[nodiscard]] constexpr auto contains(int me, int b) const -> bool { return !sparse || b == (me ^ shift); } + // The flat slots reachable from `me_flat` over a `ranks` x `parts` world. One expression per field: + // sparse names the peer rank's `parts` slots, dense is the same with peer rank 0 and count(ranks) + // == ranks, i.e. the whole world. + [[nodiscard]] constexpr auto window(size_t me_flat, size_t ranks, size_t parts) const -> SlotWindow { + const size_t peer_rank = sparse ? ((me_flat / parts) ^ static_cast(shift)) : 0; + return SlotWindow{.base = peer_rank * parts, + .count = static_cast(count(static_cast(ranks))) * parts}; + } +}; + // Argument bundles for the variable all-to-all verbs, deliberately here rather than in HybridComm.h: // ShmComm.h takes the resolve bundle and compiles in non-MPI builds, so neither bundle may name an // MPI-only type. MPI_Datatype therefore stays a separate parameter on the HybridComm verbs that need diff --git a/cpp/monoprop/detail/mpi/Exchange.h b/cpp/monoprop/detail/mpi/Exchange.h index 9ec32530..2ecf306c 100644 --- a/cpp/monoprop/detail/mpi/Exchange.h +++ b/cpp/monoprop/detail/mpi/Exchange.h @@ -14,11 +14,15 @@ #pragma once +#include #include #include #include #include "monoprop/detail/mpi/MPICompat.h" +#ifdef monoprop_ENABLE_MPI +#include "monoprop/detail/mpi/Pairwise.h" +#endif // Keeps #ifdef monoprop_ENABLE_MPI out of the consumers; non-MPI builds get self-copy stubs. @@ -28,11 +32,21 @@ namespace monoprop::mpi { // whatever the span holds, so a layout built for a differently sized communicator reads out of bounds. auto check_exchange_layout_width(std::span send_counts, const Comm &comm) -> void; +// Legs carrying a payload in either direction, i.e. an upper bound on what the pairwise path posts. +[[nodiscard]] inline auto active_leg_count(const int *send_counts, const int *recv_counts, int num_ranks) -> int { + int legs = 0; + for (int i = 0; i < num_ranks; ++i) { + legs += static_cast(send_counts[i] != 0 || recv_counts[i] != 0); + } + return legs; +} + // Idempotent completion handle for a posted payload transfer; move-only, so a request is waited on -// exactly once. wait() is a no-op on the blocking path and in non-MPI builds. Owns its request: the -// destructor completes anything still in flight, because a dropped in-flight MPI_Ialltoallv -- what an +// exactly once. wait() is a no-op on the blocking path and in non-MPI builds. Owns its requests: the +// destructor completes anything still in flight, because a dropped in-flight transfer -- what an // exception between post and wait does -- keeps writing into a thread_local buffer the next exchange -// reallocates. +// reallocates. The pairwise arm's request vector lives here for the same reason: MPI reads it until the +// wait, and a vector move keeps its heap block, so the handle can travel. class [[nodiscard("call wait() on the Ticket to complete the posted transfer")]] Ticket { public: Ticket() = default; @@ -45,6 +59,8 @@ class [[nodiscard("call wait() on the Ticket to complete the posted transfer")]] wait(); // never drop a request this handle already owns request_ = other.request_; other.request_ = MPI_REQUEST_NULL; + requests_ = std::move(other.requests_); + posted_ = std::exchange(other.posted_, 0); } #endif (void)other; @@ -58,22 +74,46 @@ class [[nodiscard("call wait() on the Ticket to complete the posted transfer")]] MPI_Wait(&request_, MPI_STATUS_IGNORE); request_ = MPI_REQUEST_NULL; } + if (posted_ != 0) { + MPI_Waitall(posted_, requests_.data(), MPI_STATUSES_IGNORE); + posted_ = 0; + } +#endif + } + + // Requests wait() still has to drain: 1 for the collective, two per pairwise leg, 0 for nothing + // posted. The only handle on which transport a post took. + [[nodiscard]] auto in_flight() const -> int { +#ifdef monoprop_ENABLE_MPI + return static_cast(request_ != MPI_REQUEST_NULL) + posted_; +#else + return 0; #endif } #ifdef monoprop_ENABLE_MPI explicit Ticket(MPI_Request request) : request_(request) {} + Ticket(std::vector requests, int posted) : requests_(std::move(requests)), posted_(posted) {} private: - MPI_Request request_ = MPI_REQUEST_NULL; + MPI_Request request_ = MPI_REQUEST_NULL; // the dense collective + std::vector requests_; // the pairwise arm's pairs; `posted_` of them are live + int posted_ = 0; #endif }; -// Never skipped on zero total: all ranks must participate or the collective deadlocks. Non-blocking -// (MPI_Ialltoallv) in an MPI build (the Ticket completes it); non-MPI build does a per-rank self-copy -// (recv layout == send layout). +// Never skipped on zero total: the collective arm needs all ranks or it deadlocks. Non-blocking in an +// MPI build -- MPI_Ialltoallv, or Isend/Irecv over the legs that carry a payload (the Ticket completes +// either); non-MPI build does a per-rank self-copy (recv layout == send layout). +// +// `wire_bits` picks the transport and MUST be RANK-UNIFORM: a rank choosing MPI_Ialltoallv waits forever +// on ranks that chose point-to-point, and no predicate over a rank's OWN row can promise that (rows vary, +// so any threshold on one straddles). It is the resolved linear-routing bit count when that routing gives +// fanout 1 -- one destination rank per generator, which is what empties the other legs -- and 0, today's +// collective, for every other geometry. The caller owns that derivation; see Evolution.cpp. template -inline auto post_flat_alltoallv(const FlatAlltoallvArgs &args, int num_ranks, Comm comm) -> Ticket { +inline auto post_flat_alltoallv(const FlatAlltoallvArgs &args, int num_ranks, Comm comm, int wire_bits = 0) + -> Ticket { // The in-process transports address the buffers as raw bytes; MPI_Ialltoallv below still takes the // typed pointers plus a datatype. Offsets stay in elements on both paths. if (comm.kind == Comm::Kind::Shm) { @@ -91,10 +131,34 @@ inline auto post_flat_alltoallv(const FlatAlltoallvArgs &args, int num_ranks, } #ifdef monoprop_ENABLE_MPI if (comm.kind == Comm::Kind::Hybrid) { - comm.hyb->alltoallv(comm.shm_rank, args.bytes(), datatype::get()); + // The wire is narrowed inside the verb, not here: only partition 0 reaches MPI, and it cannot + // name the rank's peer from its own row alone (its row may be the empty one). See HybridComm. + comm.hyb->alltoallv(comm.shm_rank, args.bytes(), datatype::get(), PeerPlan{}, wire_bits); return Ticket{}; } - (void)num_ranks; + if (wire_bits > 0) { + // Which legs to drop needs no plan and no count round: the count matrix is symmetric, so what + // this rank sends a peer IS that peer's recv count and both ends drop the same legs on the same + // value. The plan stays DENSE -- it walks all N and posts only the non-zero legs, which is + // exactly that -- so a mis-derived shift cannot drop a block here; `wire_bits` only chooses the + // transport. `legs` sizes the request vector, which a dense plan would otherwise take to 2N. + const int legs = active_leg_count(args.send_counts, args.recv_counts, num_ranks); + std::vector requests; + const int posted = sparse_pairwise(PeerPlan{}, + rank(comm), + num_ranks, + comm.mpi, + kFlatReplayTag, + datatype::get(), + sizeof(T), + reinterpret_cast(args.send), + PeerLayout{.counts = args.send_counts, .displs = args.send_displs}, + reinterpret_cast(args.recv), + PeerLayout{.counts = args.recv_counts, .displs = args.recv_displs}, + requests, + legs); + return Ticket(std::move(requests), posted); + } MPI_Request request = MPI_REQUEST_NULL; MPI_Ialltoallv(args.send, args.send_counts, diff --git a/cpp/monoprop/detail/mpi/HybridComm.h b/cpp/monoprop/detail/mpi/HybridComm.h index 47a79880..a373735f 100644 --- a/cpp/monoprop/detail/mpi/HybridComm.h +++ b/cpp/monoprop/detail/mpi/HybridComm.h @@ -32,6 +32,7 @@ #include "monoprop/detail/mpi/CheckedCount.h" #include "monoprop/detail/mpi/Comm.h" +#include "monoprop/detail/mpi/Pairwise.h" #include "monoprop/detail/mpi/PartitionBarrier.h" // Composes R MPI ranks x S in-process partitions into one flat P=R*S SPMD world. Global id is rank-major @@ -97,29 +98,47 @@ class HybridComm { auto operator=(const HybridComm &) -> HybridComm & = delete; auto size() const -> int { return r_ * s_; } + auto ranks() const -> int { return r_; } + auto partitions() const -> int { return s_; } auto global_rank(int local_partition) const -> int { return mpi_rank_ * s_ + local_partition; } - auto alltoall_counts(int local_partition, const int *send_counts /*[P]*/, int *recv_counts /*[P]*/) -> void { - guard_partition0_(local_partition, "alltoall_counts", [this, local_partition, send_counts, recv_counts] { - alltoall_counts_impl_(local_partition, send_counts, recv_counts); + auto alltoall_counts(int local_partition, + const int *send_counts /*[P]*/, + int *recv_counts /*[P]*/, + PeerPlan plan = {}) -> void { + guard_partition0_(local_partition, "alltoall_counts", [this, local_partition, send_counts, recv_counts, plan] { + alltoall_counts_impl_(local_partition, send_counts, recv_counts, plan); }); } // See AlltoallvArgs for the send-buffer lifetime and the element-vs-byte convention; `dt` is the MPI // datatype whose extent is args.elem, and it stays a separate argument because the bundle is shared // with the non-MPI-capable transport. - auto alltoallv(int local_partition, const AlltoallvArgs &args, MPI_Datatype dt) -> void { - guard_partition0_(local_partition, "alltoallv", [this, local_partition, &args, dt] { - alltoallv_impl_(local_partition, args, dt); + // + // `derive_wire_bits` > 0 asks partition 0 to narrow the WIRE itself, to that many linear bits, from + // the destination ranks the whole rank actually uses. A caller cannot supply that plan: only + // partition 0 reaches MPI, and its own row may be the empty one while a sibling partition has the + // rank's only traffic. Legal only for a SYMMETRIC layout (recv counts are the send counts), which is + // what lets the peer set be read off the published recv rows; asserted against the send rows. + auto alltoallv(int local_partition, + const AlltoallvArgs &args, + MPI_Datatype dt, + PeerPlan plan = {}, + int derive_wire_bits = 0) -> void { + guard_partition0_(local_partition, "alltoallv", [this, local_partition, &args, dt, plan, derive_wire_bits] { + alltoallv_impl_(local_partition, args, dt, plan, derive_wire_bits); }); } // See AlltoallvResolveArgs: the recv side is an output, and args.recv is resized here. template - auto alltoallv_resolve(int local_partition, const AlltoallvResolveArgs &args, MPI_Datatype dt) -> void { + auto alltoallv_resolve(int local_partition, + const AlltoallvResolveArgs &args, + MPI_Datatype dt, + PeerPlan plan = {}) -> void { // `args` by reference, not by value: the impl resizes args.recv and then writes through it. - guard_partition0_(local_partition, "alltoallv_resolve", [this, local_partition, &args, dt] { - alltoallv_resolve_impl_(local_partition, args, dt); + guard_partition0_(local_partition, "alltoallv_resolve", [this, local_partition, &args, dt, plan] { + alltoallv_resolve_impl_(local_partition, args, dt, plan); }); } @@ -170,18 +189,25 @@ class HybridComm { } } - // recv_counts[g] = amount global partition g sends to this partition. 2 barriers + one S*S-int MPI_Alltoall. - auto alltoall_counts_impl_(int local_partition, const int *send_counts /*[P]*/, int *recv_counts /*[P]*/) -> void { + // recv_counts[g] = amount global partition g sends to this partition. 2 barriers + one count round. + auto alltoall_counts_impl_(int local_partition, + const int *send_counts /*[P]*/, + int *recv_counts /*[P]*/, + PeerPlan plan) -> void { publish_counts_row_(local_partition, send_counts); sync(); if (local_partition == 0) { - pack_count_matrix_(); - MPI_Alltoall(counts_send_.data(), s_ * s_, MPI_INT, counts_recv_.data(), s_ * s_, MPI_INT, parent_); + fill_peers_(plan); + pack_count_matrix_(plan); + exchange_count_blocks_(plan); } sync(); // Partition t extracts its row: recv from (rank a, partition su) is contiguous per source rank a. + // Under a plan only the f peer ranks were exchanged, so the rest of the row is zero by definition + // (a non-peer cannot own the partner of any term this rank owns). const int t = local_partition; - for (int a = 0; a < r_; ++a) { + std::fill(recv_counts, recv_counts + static_cast(r_) * static_cast(s_), 0); + for (const int a : peers_) { for (int su = 0; su < s_; ++su) { recv_counts[a * s_ + su] = counts_recv_[counts_idx_(a, t, su)]; } @@ -192,19 +218,31 @@ class HybridComm { } // Flat variable all-to-all over caller-owned buffers; see AlltoallvArgs for the conventions. - auto alltoallv_impl_(int local_partition, const AlltoallvArgs &args, MPI_Datatype dt) -> void { + auto alltoallv_impl_(int local_partition, + const AlltoallvArgs &args, + MPI_Datatype dt, + PeerPlan plan, + int derive_wire_bits = 0) -> void { const size_t u = static_cast(local_partition); Slot &me = slots_[u]; me.ptr = args.send; me.send_displs = args.send_displs; publish_counts_row_(local_partition, args.send_counts); - publish_recv_rows_(local_partition, args.recv_counts); + publish_recv_rows_(local_partition, args.recv_counts, plan); sync(); // B1 // B2: partition 0 sizes/reallocates staging; must finish before any partition packs into stage_send_. if (local_partition == 0) { + // Written here, read again in the B3->B4 window: partition 0 is this member's only toucher. + wire_plan_ = plan; + if (derive_wire_bits > 0) { + wire_plan_ = derived_wire_plan_(derive_wire_bits); + // The recv rows it was read off against the send rows: the symmetry the parameter needs. + assert(narrowing_is_lossless_(wire_plan_)); + } + fill_peers_(wire_plan_); size_staging_send_(args.elem); - fill_recv_col_from_rows_(); + fill_recv_col_([this](int a, int t) { return row_recv_(t)[a]; }); size_staging_recv_(args.elem); } sync(); // B2 @@ -213,46 +251,29 @@ class HybridComm { pack_send_(local_partition, args.elem); sync(); // B3 - // B4: partition 0 runs the single MPI_Alltoallv while peers park at the barrier. + // B4: partition 0 moves the payload while peers park at the barrier. if (local_partition == 0) { - MPI_Alltoallv(stage_send_.data(), - mpi_send_counts_.data(), - mpi_send_displs_.data(), - dt, - stage_recv_.data(), - mpi_recv_counts_.data(), - mpi_recv_displs_.data(), - dt, - parent_); + exchange_payload_(dt, args.elem, wire_plan_); } sync(); // B4 - // Scatter each global source's contiguous run from stage_recv_ to recv_displs[g] (all legs, incl. - // self-rank, go through staging). Walks (a, su) in accumulation order, so `cur` re-derives the - // block starts from base_recv_ and this partition's own counts. - std::byte *dst = args.recv; - const int t = local_partition; - for (int a = 0; a < r_; ++a) { - size_t cur = base_recv_[static_cast(a) * static_cast(s_) + static_cast(t)]; - for (int su = 0; su < s_; ++su) { - const int g = a * s_ + su; - const int cnt = args.recv_counts[g]; - if (cnt != 0) { - std::memcpy(dst + static_cast(args.recv_displs[g]) * args.elem, - stage_recv_.data() + cur * args.elem, - static_cast(cnt) * args.elem); - } - cur += static_cast(cnt); - } - } + scatter_recv_(local_partition, args.recv, args.recv_counts, args.recv_displs, args.elem); // No trailing barrier: base_recv_ is rewritten only in a later verb's B1→B2 window. } // Fused count-resolve + payload alltoallv: folds the standalone count exchange into this verb's - // B1→B2 window (4 syncs instead of 6). recv_counts / recv_displs and `recv` (resized) are outputs. - // Bit-identical to alltoall_counts + alltoallv. + // barriered windows (4 syncs instead of 6). recv_counts / recv_displs and `recv` (resized) are + // outputs. Bit-identical to alltoall_counts + alltoallv. + // + // The count round is POSTED in B1→B2 and only waited on in B3→B4, so the whole B2→B3 packing runs + // underneath it. Split, not fused, because nothing before fill_recv_col_ reads counts_recv_: the + // send side sizes from the locally published rows, and pack_send_ from that sizing. The dense arm + // is MPI_Alltoall, blocking, and completes inside post_count_blocks_ regardless. template - auto alltoallv_resolve_impl_(int local_partition, const AlltoallvResolveArgs &args, MPI_Datatype dt) -> void { + auto alltoallv_resolve_impl_(int local_partition, + const AlltoallvResolveArgs &args, + MPI_Datatype dt, + PeerPlan plan) -> void { // Typed verb: element bytes are sizeof(T) by construction, so they are derived rather than passed. constexpr size_t elem = sizeof(T); const size_t u = static_cast(local_partition); @@ -261,22 +282,40 @@ class HybridComm { // never reconstructs T, so the slot stays type-erased for the untyped alltoallv_impl_ above. me.ptr = reinterpret_cast(args.send); me.send_displs = args.send_displs; - // Count row only: the recv counts do not exist until the count Alltoall in B1→B2. + // Count row only: the recv counts do not exist until the count round is drained in B3→B4. publish_counts_row_(local_partition, args.send_counts); sync(); // B1 if (local_partition == 0) { - pack_count_matrix_(); - MPI_Alltoall(counts_send_.data(), s_ * s_, MPI_INT, counts_recv_.data(), s_ * s_, MPI_INT, parent_); + fill_peers_(plan); + pack_count_matrix_(plan); + post_count_blocks_(plan); size_staging_send_(elem); - fill_recv_col_from_counts_recv_(); - size_staging_recv_(elem); } sync(); // B2 + // B3: each partition packs its own cross-rank blocks into stage_send_, the count round in flight. + pack_send_(local_partition, elem); + sync(); // B3 + + // B4: the counts land here -- fill_recv_col_ is their first reader -- then the payload moves. + if (local_partition == 0) { + wait_count_blocks_(); + fill_recv_col_([this](int a, int t) { return block_sum_(a, t); }); + size_staging_recv_(elem); + exchange_payload_(dt, elem, plan); + } + sync(); // B4 + + // Past B4 now, not before B3: counts_recv_ does not exist until the wait above. Partition 0 + // cannot rewrite it before a later verb's B1→B2 window, unreachable until every reader here + // has arrived at that verb's B1. const int t = local_partition; long long total = 0; - for (int a = 0; a < r_; ++a) { + const size_t p = static_cast(r_) * static_cast(s_); + std::fill(args.recv_counts, args.recv_counts + p, 0); + std::fill(args.recv_displs, args.recv_displs + p, 0); + for (const int a : peers_) { for (int su = 0; su < s_; ++su) { const int g = a * s_ + su; const int c = counts_recv_[counts_idx_(a, t, su)]; @@ -287,36 +326,11 @@ class HybridComm { } args.recv.resize(static_cast(checked_mpi_count(total, "Total recv count"))); - pack_send_(local_partition, elem); - sync(); // B3 - - if (local_partition == 0) { - MPI_Alltoallv(stage_send_.data(), - mpi_send_counts_.data(), - mpi_send_displs_.data(), - dt, - stage_recv_.data(), - mpi_recv_counts_.data(), - mpi_recv_displs_.data(), - dt, - parent_); - } - sync(); // B4 - - std::byte *dst = reinterpret_cast(args.recv.data()); // after the resize: it may reallocate - for (int a = 0; a < r_; ++a) { - size_t cur = base_recv_[static_cast(a) * static_cast(s_) + static_cast(t)]; - for (int su = 0; su < s_; ++su) { - const int g = a * s_ + su; - const int cnt = args.recv_counts[g]; - if (cnt != 0) { - std::memcpy(dst + static_cast(args.recv_displs[g]) * elem, - stage_recv_.data() + cur * elem, - static_cast(cnt) * elem); - } - cur += static_cast(cnt); - } - } + scatter_recv_(local_partition, + reinterpret_cast(args.recv.data()), // after the resize: it may reallocate + args.recv_counts, + args.recv_displs, + elem); // No trailing barrier: same discipline as alltoallv_impl_. } @@ -408,6 +422,35 @@ class HybridComm { + static_cast(su); } + // The plan's peer ranks, materialised once per verb: the sizing sweeps index them S times each. + // Written by partition 0 in the B1→B2 window, like base_recv_, so every reader past B2 sees it. + auto fill_peers_(PeerPlan plan) -> void { + const int f = plan.count(r_); + peers_.resize(static_cast(f)); + for (int k = 0; k < f; ++k) { + peers_[static_cast(k)] = plan.peer(mpi_rank_, k); + } + } + + // Zero only the peer ranks' slots of a [P] table. Every sweep that writes one of these tables and + // every sweep that reads it walks the same peer set, so the rest of the row is never looked at -- + // and these run SERIALLY on partition 0 while S-1 partitions park, so the width is the cost. + auto zero_peer_slots_(std::vector &table) -> void { + for (const int b : peers_) { + std::fill_n(table.begin() + (static_cast(b) * s_), s_, 0LL); + } + } + + // What rank a's block of the count message holds for partition t, summed over source partitions. + auto block_sum_(int a, int t) const -> long long { + const int *blk = counts_recv_.data() + counts_idx_(a, t, 0); + long long sum = 0; + for (int su = 0; su < s_; ++su) { + sum += blk[su]; + } + return sum; + } + // [S x P] payload offset table: row u is source partition u's staging starts, one per destination g. auto pack_idx_(int u, int g) const -> size_t { return static_cast(u) * static_cast(r_) * static_cast(s_) + static_cast(g); @@ -425,10 +468,16 @@ class HybridComm { static_cast(r_) * static_cast(s_) * sizeof(int)); } - // long long: it sums S int counts. - auto publish_recv_rows_(int local_partition, const int *recv_counts) -> void { + // long long: it sums S int counts. Masked through the plan, symmetric with its one reader + // (fill_recv_col_): a non-peer's row is zero by definition, and a count left there would size + // staging for a block no receive is posted for. + auto publish_recv_rows_(int local_partition, const int *recv_counts, PeerPlan plan) -> void { long long *rr = row_recv_(local_partition); - for (int a = 0; a < r_; ++a) { + // Every entry is written, not just the peers': derived_wire_plan_ reads the whole row. + std::fill_n(rr, r_, 0LL); + const int f = plan.count(r_); + for (int k = 0; k < f; ++k) { + const int a = plan.peer(mpi_rank_, k); long long sum = 0; for (int su = 0; su < s_; ++su) { sum += recv_counts[a * s_ + su]; @@ -440,10 +489,11 @@ class HybridComm { // Transpose the published count rows into counts_send_, dest-major then source-minor, for the one // S*S-int MPI_Alltoall. Partition 0 only, inside a barriered window. Source partition outer, so the // peer-owned side streams. Every element is written here, so no pre-zeroing. - auto pack_count_matrix_() -> void { + auto pack_count_matrix_([[maybe_unused]] PeerPlan plan) -> void { + assert(narrowing_is_lossless_(plan)); // a wrong shift every rank agrees on drops blocks silently for (int su = 0; su < s_; ++su) { const int *row = counts_row_(su); - for (int b = 0; b < r_; ++b) { + for (const int b : peers_) { for (int t = 0; t < s_; ++t) { counts_send_[counts_idx_(b, t, su)] = row[b * s_ + t]; } @@ -451,6 +501,118 @@ class HybridComm { } } + // The rank-level peer set, read off the recv rows every partition published before B1 -- the first + // point with a view wider than one partition's row. Under fanout-1 routing a layer's traffic is all + // on ONE rank, so this resolves to a shift; with nothing occupied it resolves to the self peer, whose + // legs are then all zero, and that keeps the collective-vs-pairwise branch a function of + // `derive_wire_bits` alone rather than of a rank's data (a data-dependent branch straddles and hangs). + // A set wider than one rank contradicts the caller's fanout claim: dense, so nothing is dropped. + auto derived_wire_plan_(int bits) -> PeerPlan { + int found = -1; + for (int u = 0; u < s_; ++u) { + const long long *rr = row_recv_(u); + for (int a = 0; a < r_; ++a) { + if (rr[a] != 0 && a != found) { + if (found >= 0) { + assert(false && "fanout claimed 1, but this rank's layer spans several peer ranks"); + return PeerPlan{}; + } + found = a; + } + } + } + // The plan is a boolean now, so every rank bit is a linear bit and the mask is the rank index. + return PeerPlan{.sparse = bits > 0, .shift = found < 0 ? 0 : (mpi_rank_ ^ found)}; + } + + // Do the published rows put anything outside the plan's peers? If so the narrowing silently drops it. + auto narrowing_is_lossless_(PeerPlan plan) const -> bool { + for (int su = 0; su < s_; ++su) { + const int *row = counts_matrix_ + static_cast(su) * counts_stride_; + for (int g = 0; g < r_ * s_; ++g) { + if (row[g] != 0 && !plan.contains(mpi_rank_, g / s_)) { + return false; + } + } + } + return true; + } + + // The count blocks: one S*S-int MPI_Alltoall when dense, else a pair per peer (with the full linear + // bits a zero rank shift keeps the whole round on-rank). Partition 0 only, in a barriered window. + // + // Sparse arm POSTS ONLY, so counts_recv_ and counts_send_ must stay put until wait_count_blocks_; + // the dense MPI_Alltoall is blocking and has completed on return, which is why plan.dense() leaves + // nothing live. The requests go in count_reqs_, never reqs_: see the member declaration. + auto post_count_blocks_(PeerPlan plan) -> void { + assert(count_posted_ == 0); // an un-drained round would be waited on twice + const int block = s_ * s_; + if (plan.dense()) { + MPI_Alltoall(counts_send_.data(), block, MPI_INT, counts_recv_.data(), block, MPI_INT, parent_); + return; + } + const PeerLayout blocks{.block = block}; + count_posted_ = sparse_pairwise(plan, + mpi_rank_, + r_, + parent_, + kHybridCountTag, + MPI_INT, + sizeof(int), + reinterpret_cast(counts_send_.data()), + blocks, + reinterpret_cast(counts_recv_.data()), + blocks, + count_reqs_); + } + + // Drains post_count_blocks_. A no-op on the dense arm, and on a plan whose only peer is this rank + // itself -- a count block is a fixed S*S ints, so no other peer can be skipped for a zero count. + auto wait_count_blocks_() -> void { + if (count_posted_ != 0) { + MPI_Waitall(count_posted_, count_reqs_.data(), MPI_STATUSES_IGNORE); + count_posted_ = 0; + } + } + + // Post and drain in one step, for callers with no work to overlap. + auto exchange_count_blocks_(PeerPlan plan) -> void { + post_count_blocks_(plan); + wait_count_blocks_(); + } + + // The staged payload: one MPI_Alltoallv when dense, else a pair per peer over the same per-rank + // counts and displacements (a non-peer's count is zero, so nothing is dropped). `elem` is dt's + // extent: needed to reach a block, and known exactly to both callers. + auto exchange_payload_(MPI_Datatype dt, size_t elem, PeerPlan plan) -> void { + if (plan.dense()) { + MPI_Alltoallv(stage_send_.data(), + mpi_send_counts_.data(), + mpi_send_displs_.data(), + dt, + stage_recv_.data(), + mpi_recv_counts_.data(), + mpi_recv_displs_.data(), + dt, + parent_); + return; + } + const int posted = + sparse_pairwise(plan, + mpi_rank_, + r_, + parent_, + kHybridPayloadTag, + dt, + elem, + stage_send_.data(), + PeerLayout{.counts = mpi_send_counts_.data(), .displs = mpi_send_displs_.data()}, + stage_recv_.data(), + PeerLayout{.counts = mpi_recv_counts_.data(), .displs = mpi_recv_displs_.data()}, + reqs_); + MPI_Waitall(posted, reqs_.data(), MPI_STATUSES_IGNORE); + } + template static auto grow_(V &v, size_t need) -> void { if (v.size() < need) { @@ -460,17 +622,25 @@ class HybridComm { // Partition 0's send-side staging sizing, between B1 and B2, in two sweeps of counts_matrix_. Wire // block order is destination major, source minor: a per-destination base plus a per-source prefix. + // Both sweeps run SERIALLY here while S-1 partitions park, so narrowing them to the peers matters as + // much as the message count does. auto size_staging_send_(size_t elem) -> void { - const size_t p = static_cast(r_) * static_cast(s_); // Pass A: the column sums W over source partitions, u outer so both sides sweep in address order. - std::fill(col_sum_.begin(), col_sum_.end(), 0LL); + zero_peer_slots_(col_sum_); for (int u = 0; u < s_; ++u) { const int *row = counts_row_(u); - for (size_t g = 0; g < p; ++g) { - col_sum_[g] += row[g]; + for (const int b : peers_) { + const size_t base = static_cast(b) * static_cast(s_); + for (int t = 0; t < s_; ++t) { + col_sum_[base + static_cast(t)] += row[base + static_cast(t)]; + } } } - for (int b = 0; b < r_; ++b) { + // The [R] arrays stay fully zeroed -- MPI_Alltoallv reads every entry on the dense arm, and the + // zeros are what make the prefix below equal the full 0..R one at the peers' positions. + std::ranges::fill(mpi_send_counts_, 0); + std::ranges::fill(mpi_send_displs_, 0); + for (const int b : peers_) { const long long *col = col_sum_.data() + static_cast(b) * static_cast(s_); long long send_sum = 0; for (int t = 0; t < s_; ++t) { @@ -478,13 +648,15 @@ class HybridComm { } mpi_send_counts_[static_cast(b)] = checked_mpi_count(send_sum, "Per-rank send count"); } + // peers_ is ascending (dense is 0..R-1, sparse is a singleton), so a prefix over it takes the + // same value at every peer as a prefix over all R: a non-peer contributes zero. long long send_running = 0; - for (int b = 0; b < r_; ++b) { + for (const int b : peers_) { mpi_send_displs_[static_cast(b)] = checked_mpi_count(send_running, "Send displacement"); send_running += mpi_send_counts_[static_cast(b)]; } const size_t total_send = static_cast(checked_mpi_count(send_running, "Total send count")); - for (int b = 0; b < r_; ++b) { + for (const int b : peers_) { size_t cur = static_cast(mpi_send_displs_[static_cast(b)]); for (int t = 0; t < s_; ++t) { const size_t g = static_cast(b) * static_cast(s_) + static_cast(t); @@ -493,37 +665,31 @@ class HybridComm { } } // Pass B: the exclusive prefix over source partitions; col_sum_ is free to be reused for it here. - std::fill(col_sum_.begin(), col_sum_.end(), 0LL); + zero_peer_slots_(col_sum_); for (int u = 0; u < s_; ++u) { const int *row = counts_row_(u); size_t *off = pack_off_.data() + pack_idx_(u, 0); - for (size_t g = 0; g < p; ++g) { - off[g] = base_send_[g] + static_cast(col_sum_[g]); - col_sum_[g] += row[g]; + for (const int b : peers_) { + const size_t base = static_cast(b) * static_cast(s_); + for (int t = 0; t < s_; ++t) { + const size_t g = base + static_cast(t); + off[g] = base_send_[g] + static_cast(col_sum_[g]); + col_sum_[g] += row[g]; + } } } // Grow-only, no zero-fill: pack_send_'s blocks tile [0, total_send) exactly. grow_(stage_send_, total_send * elem); } - // recv_col_[a*S + t] = what partition t receives from rank a: the rows published in Phase P0. - auto fill_recv_col_from_rows_() -> void { - for (int a = 0; a < r_; ++a) { + // recv_col_[a*S + t] = what partition t receives from rank a. `value(a, t)` reads it from the rows + // published in Phase P0 (alltoallv) or from the count blocks just exchanged (the fused resolve). + template + auto fill_recv_col_(Value &&value) -> void { + zero_peer_slots_(recv_col_); + for (const int a : peers_) { for (int t = 0; t < s_; ++t) { - recv_col_[static_cast(a) * static_cast(s_) + static_cast(t)] = row_recv_(t)[a]; - } - } - } - - auto fill_recv_col_from_counts_recv_() -> void { - for (int a = 0; a < r_; ++a) { - for (int t = 0; t < s_; ++t) { - const int *blk = counts_recv_.data() + counts_idx_(a, t, 0); - long long sum = 0; - for (int su = 0; su < s_; ++su) { - sum += blk[su]; - } - recv_col_[static_cast(a) * static_cast(s_) + static_cast(t)] = sum; + recv_col_[static_cast(a) * static_cast(s_) + static_cast(t)] = value(a, t); } } } @@ -531,20 +697,23 @@ class HybridComm { // Partition 0's recv-side staging sizing, from recv_col_. Only the per-(rank, partition) base; the // post-B4 scatter re-derives the per-source offsets as it walks (a, su). auto size_staging_recv_(size_t elem) -> void { - for (int a = 0; a < r_; ++a) { + std::ranges::fill(mpi_recv_counts_, 0); + std::ranges::fill(mpi_recv_displs_, 0); + for (const int a : peers_) { long long recv_sum = 0; for (int t = 0; t < s_; ++t) { recv_sum += recv_col_[static_cast(a) * static_cast(s_) + static_cast(t)]; } mpi_recv_counts_[static_cast(a)] = checked_mpi_count(recv_sum, "Per-rank recv count"); } + // Same ascending-peers prefix as the send side. long long recv_running = 0; - for (int a = 0; a < r_; ++a) { + for (const int a : peers_) { mpi_recv_displs_[static_cast(a)] = checked_mpi_count(recv_running, "Recv displacement"); recv_running += mpi_recv_counts_[static_cast(a)]; } const size_t total_recv = static_cast(checked_mpi_count(recv_running, "Total recv count")); - for (int a = 0; a < r_; ++a) { + for (const int a : peers_) { size_t cur = static_cast(mpi_recv_displs_[static_cast(a)]); for (int t = 0; t < s_; ++t) { const size_t g = static_cast(a) * static_cast(s_) + static_cast(t); @@ -563,13 +732,37 @@ class HybridComm { const int *my_send_counts = counts_row_(u); const int *my_send_displs = slots_[static_cast(u)].send_displs; const size_t *off = pack_off_.data() + pack_idx_(u, 0); - const int p = r_ * s_; - for (int g = 0; g < p; ++g) { - const int cnt = my_send_counts[g]; - if (cnt != 0) { - std::memcpy(stage_send_.data() + off[g] * elem, - src + static_cast(my_send_displs[g]) * elem, - static_cast(cnt) * elem); + for (const int b : peers_) { + const int base = b * s_; + for (int t = 0; t < s_; ++t) { + const int g = base + t; + const int cnt = my_send_counts[g]; + if (cnt != 0) { + std::memcpy(stage_send_.data() + off[g] * elem, + src + static_cast(my_send_displs[g]) * elem, + static_cast(cnt) * elem); + } + } + } + } + + // Scatter each global source's contiguous run out of stage_recv_ to recv_displs[g] (all legs, incl. + // self-rank, go through staging). Walks (a, su) in accumulation order, so `cur` re-derives the block + // starts from base_recv_ and this partition's own counts. + auto scatter_recv_(int local_partition, std::byte *dst, const int *recv_counts, const int *recv_displs, size_t elem) + -> void { + const int t = local_partition; + for (const int a : peers_) { + size_t cur = base_recv_[static_cast(a) * static_cast(s_) + static_cast(t)]; + for (int su = 0; su < s_; ++su) { + const int g = a * s_ + su; + const int cnt = recv_counts[g]; + if (cnt != 0) { + std::memcpy(dst + static_cast(recv_displs[g]) * elem, + stage_recv_.data() + cur * elem, + static_cast(cnt) * elem); + } + cur += static_cast(cnt); } } } @@ -627,6 +820,17 @@ class HybridComm { double red_f64_ = 0.0; uint64_t red_u64_ = 0; std::vector red_vec_; + // Point-to-point request scratch for the sparse payload round; grown on demand, partition 0 only. + std::vector reqs_; + // The count round's own scratch, separate from reqs_ by construction and not merely by the current + // ordering: it stays live across B2 and B3, and Pairwise.h's resize would move the buffer MPI holds + // pointers into the moment a payload post ever preceded the count wait. + std::vector count_reqs_; + int count_posted_ = 0; // live requests in count_reqs_; always 0 on the dense (blocking) arm + // This verb's peer ranks; see fill_peers_. + std::vector peers_; + // alltoallv's wire plan, partition 0 only: written in B1->B2, read in B3->B4. See derived_wire_plan_. + PeerPlan wire_plan_; PartitionBarrier barrier_; }; diff --git a/cpp/monoprop/detail/mpi/MPICompat.cpp b/cpp/monoprop/detail/mpi/MPICompat.cpp index 133a9fb7..c9285c40 100644 --- a/cpp/monoprop/detail/mpi/MPICompat.cpp +++ b/cpp/monoprop/detail/mpi/MPICompat.cpp @@ -14,10 +14,15 @@ #include "monoprop/detail/mpi/Exchange.h" +#include #include #include #include +#ifdef monoprop_ENABLE_MPI +#include "monoprop/detail/mpi/Pairwise.h" +#endif + namespace monoprop::mpi { #ifdef monoprop_ENABLE_MPI @@ -84,6 +89,20 @@ auto size(const Comm &comm) -> int { #endif } +auto geometry(const Comm &comm) -> Geometry { + if (comm.kind == Comm::Kind::Shm) { + return {.ranks = 1, .partitions = comm.shm->size()}; + } +#ifdef monoprop_ENABLE_MPI + if (comm.kind == Comm::Kind::Hybrid) { + return {.ranks = comm.hyb->ranks(), .partitions = comm.hyb->partitions()}; + } + return {.ranks = size(comm), .partitions = 1}; +#else + return {.ranks = 1, .partitions = 1}; +#endif +} + auto allreduce_sum_inplace(VecD &values, Comm comm) -> void { if (comm.kind == Comm::Kind::Shm) { comm.shm->allreduce_sum_inplace(comm.shm_rank, values.data(), values.size()); @@ -100,19 +119,45 @@ auto allreduce_sum_inplace(VecD &values, Comm comm) -> void { #endif } -auto alltoall_counts(const int *send_counts, int *recv_counts, int n, Comm comm) -> void { +auto alltoall_counts(const int *send_counts, int *recv_counts, int n, Comm comm, PeerPlan plan) -> void { if (comm.kind == Comm::Kind::Shm) { comm.shm->alltoall_counts(comm.shm_rank, send_counts, recv_counts); return; } #ifdef monoprop_ENABLE_MPI if (comm.kind == Comm::Kind::Hybrid) { - comm.hyb->alltoall_counts(comm.shm_rank, send_counts, recv_counts); + comm.hyb->alltoall_counts(comm.shm_rank, send_counts, recv_counts, plan); + return; + } + if (!plan.dense()) { + // S == 1 world: exchange one int with each reachable peer; the rest of the row is zero by + // definition, so it must be cleared rather than left from a previous round. + int me = 0; + MPI_Comm_rank(comm.mpi, &me); + std::fill(recv_counts, recv_counts + n, 0); + const PeerLayout one{.block = 1}; + // Eager by contract: recv_counts is caller memory the caller reads on return, so unlike the + // payload round this one cannot be handed on in a handle. + std::vector reqs; + const int posted = sparse_pairwise(plan, + me, + n, + comm.mpi, + kFlatCountTag, + MPI_INT, + sizeof(int), + reinterpret_cast(send_counts), + one, + reinterpret_cast(recv_counts), + one, + reqs); + MPI_Waitall(posted, reqs.data(), MPI_STATUSES_IGNORE); return; } (void)n; MPI_Alltoall(send_counts, 1, MPI_INT, recv_counts, 1, MPI_INT, comm.mpi); #else + (void)plan; // single participant: nothing to narrow for (int i = 0; i < n; ++i) { recv_counts[i] = send_counts[i]; } diff --git a/cpp/monoprop/detail/mpi/MPICompat.h b/cpp/monoprop/detail/mpi/MPICompat.h index 288ee6c8..11aabb6e 100644 --- a/cpp/monoprop/detail/mpi/MPICompat.h +++ b/cpp/monoprop/detail/mpi/MPICompat.h @@ -15,6 +15,7 @@ #pragma once #include +#include #include #include #include @@ -32,6 +33,7 @@ #include "monoprop/detail/mpi/ShmComm.h" #ifdef monoprop_ENABLE_MPI #include "monoprop/detail/mpi/HybridComm.h" +#include "monoprop/detail/mpi/Pairwise.h" #endif // These includes are here on purpose and should not be moved to the top @@ -90,6 +92,15 @@ inline auto finalize() -> void {} auto rank(const Comm &comm) -> int; auto size(const Comm &comm) -> int; +// How the flat world of size() is actually built: ranks * partitions. Routing needs the split, because +// an inter-rank message costs a network hop while an inter-partition one is a shared-memory copy -- +// size() alone cannot tell them apart. ranks * partitions == size() for every Kind. +struct Geometry { + int ranks = 1; + int partitions = 1; +}; +auto geometry(const Comm &comm) -> Geometry; + template inline auto allreduce_sum(T local_val, Comm comm) -> T { if (comm.kind == Comm::Kind::Shm) { @@ -109,8 +120,52 @@ inline auto allreduce_sum(T local_val, Comm comm) -> T { auto allreduce_sum_inplace(VecD &values, Comm comm) -> void; -// `n` is the comm size. -auto alltoall_counts(const int *send_counts, int *recv_counts, int n, Comm comm) -> void; +// `n` is the comm size. `plan` narrows the exchange to the destination ranks it can reach (see PeerPlan); +// the default is dense, i.e. today's collective. +auto alltoall_counts(const int *send_counts, int *recv_counts, int n, Comm comm, PeerPlan plan = {}) -> void; + +// The per-slot block arrays begin_alltoallv and wait_into accept: a plain [P] vector-of-vectors, or a +// WindowVec over the slots a PeerPlan can reach. These four overload pairs are the whole difference -- +// the verbs below are one code path walking one SlotWindow. +template +using SlotBlockValue = typename Blocks::value_type::value_type; + +template +inline auto slot_window_of(const std::vector> &v) -> SlotWindow { + return SlotWindow{.base = 0, .count = v.size()}; +} +template +inline auto slot_window_of(const WindowVec> &v) -> SlotWindow { + return v.window(); +} + +template +inline auto slot_block(const std::vector> &v, size_t slot) -> const std::vector & { + return v[slot]; +} +template +inline auto slot_block(std::vector> &v, size_t slot) -> std::vector & { + return v[slot]; +} +template +inline auto slot_block(const WindowVec> &v, size_t slot) -> const std::vector & { + return v.at_slot(slot); +} +template +inline auto slot_block(WindowVec> &v, size_t slot) -> std::vector & { + return v.at_slot(slot); +} + +// A plain destination keeps the full-world shape (a non-peer's block is empty, not absent); a WindowVec +// takes the round's window. +template +inline auto reset_slots(std::vector> &v, SlotWindow /*w*/, size_t world) -> void { + v.assign(world, std::vector{}); +} +template +inline auto reset_slots(WindowVec> &v, SlotWindow w, size_t /*world*/) -> void { + v.reset(w); +} // In-flight variable-size all-to-all owning its buffers + layout, so several can be in flight. // recv_counts is valid on return from begin_alltoallv; wait_into completes the payload transfer (a @@ -118,6 +173,8 @@ auto alltoall_counts(const int *send_counts, int *recv_counts, int n, Comm comm) template struct PendingAlltoallv { int num_ranks = 0; + // The slots this round touches; counts/displs are zero outside it. Set by begin_alltoallv. + SlotWindow window; std::vector send_counts; std::vector send_displs; std::vector recv_counts; @@ -125,74 +182,111 @@ struct PendingAlltoallv { std::vector send_buffer; std::vector recv_buffer; #ifdef monoprop_ENABLE_MPI - MPI_Request request = MPI_REQUEST_NULL; // set only on the Kind::Mpi async path + MPI_Request request = MPI_REQUEST_NULL; // set only on the Kind::Mpi dense async path + std::vector requests; // the Kind::Mpi sparse path's pairs; `posted` of them live + int posted = 0; // MPI reads send_buffer/recv_buffer until these complete, + // and both move with the handle, so the pointers hold #endif - auto wait_into(std::vector> &recv_data) -> void { + template + auto wait_into(Dest &recv_data) -> void { #ifdef monoprop_ENABLE_MPI if (request != MPI_REQUEST_NULL) { MPI_Wait(&request, MPI_STATUS_IGNORE); request = MPI_REQUEST_NULL; } + if (posted != 0) { + MPI_Waitall(posted, requests.data(), MPI_STATUSES_IGNORE); + posted = 0; + } #endif - recv_data.resize(static_cast(num_ranks)); - for (int i = 0; i < num_ranks; ++i) { - const auto lo = recv_buffer.begin() + recv_displs[static_cast(i)]; - recv_data[static_cast(i)].assign(lo, lo + recv_counts[static_cast(i)]); + reset_slots(recv_data, window, static_cast(num_ranks)); + for (size_t k = 0; k < window.count; ++k) { + const size_t i = window.slot(WindowIndex{k}); + const auto lo = recv_buffer.begin() + recv_displs[i]; + slot_block(recv_data, i).assign(lo, lo + recv_counts[i]); } } }; +// Debug-only: a caller may supply more slots than the plan reaches, and anything it left outside the +// window is DROPPED rather than refused -- the silent failure mode a wrong-but-agreed shift produces. +template +inline auto assert_outside_window_is_empty_([[maybe_unused]] const Blocks &send_data, + [[maybe_unused]] SlotWindow supplied, + [[maybe_unused]] SlotWindow window) -> void { +#ifndef NDEBUG + for (size_t i = supplied.base; i < supplied.stop(); ++i) { + assert((window.contains(i) || slot_block(send_data, i).empty()) + && "a block outside the plan's peer window would be dropped in silence"); + } +#endif +} + // The count exchange runs eagerly (recv_counts known on return); the Kind::Mpi payload is non-blocking // (wait_into completes it), Shm / single-process transfer here. // skip_self: do not send the self slot (the caller handles self inline) — self send/recv = 0. // known_recv_counts: recv counts already known (e.g. the transpose of the query counts), so skip the // count exchange. The self slot is also zeroed when skip_self is set. -template -inline auto begin_alltoallv(const std::vector> &send_data, +template > +inline auto begin_alltoallv(const Blocks &send_data, Comm comm, bool skip_self = false, - const std::vector *known_recv_counts = nullptr) -> PendingAlltoallv { + const std::vector *known_recv_counts = nullptr, + PeerPlan plan = {}) -> PendingAlltoallv { const int num_ranks = size(comm); - if (static_cast(send_data.size()) != num_ranks) { + const int me = rank(comm); + const auto geom = geometry(comm); + PendingAlltoallv h; + h.num_ranks = num_ranks; + // The plan IS the mask, dense included -- it is the count == P value of the same window. A caller may + // hand a whole [P] array under a sparse plan (the tests do), so the supplied array only has to COVER + // the window; assert_outside_window_is_empty_ catches what it leaves outside, which is the silent + // drop a wrong-but-agreed shift produces. + h.window = + plan.window(static_cast(me), static_cast(geom.ranks), static_cast(geom.partitions)); + const SlotWindow supplied = slot_window_of(send_data); + if (h.window.stop() > static_cast(num_ranks) || supplied.base > h.window.base + || supplied.stop() < h.window.stop()) { throw CollectiveArgumentError( - std::format("begin_alltoallv: send_data size ({}) must equal number of ranks ({})", - send_data.size(), + std::format("begin_alltoallv: send_data covers slots [{}, {}), which does not cover the plan's " + "[{}, {}) in a {}-slot world", + supplied.base, + supplied.stop(), + h.window.base, + h.window.stop(), num_ranks)); } - PendingAlltoallv h; - h.num_ranks = num_ranks; - h.send_counts.resize(static_cast(num_ranks)); - h.send_displs.resize(static_cast(num_ranks)); - h.recv_displs.resize(static_cast(num_ranks)); + assert_outside_window_is_empty_(send_data, supplied, h.window); + h.send_counts.assign(static_cast(num_ranks), 0); + h.send_displs.assign(static_cast(num_ranks), 0); + h.recv_displs.assign(static_cast(num_ranks), 0); - const int self = skip_self ? rank(comm) : -1; + const int self = skip_self ? me : -1; + // Counts and their prefix in ONE sweep over the window; the rest stay zero from the assign above. // Wide accumulator + checked narrowing: a wrapped count would size send_buffer short and then feed // MPI a negative count/displacement. - long long total_send = 0; - for (int i = 0; i < num_ranks; ++i) { - const size_t n = (i == self) ? 0 : send_data[static_cast(i)].size(); - const int c = checked_mpi_count(n, "Send count"); - h.send_counts[static_cast(i)] = c; - total_send += c; - } long long running_send = 0; - for (int i = 0; i < num_ranks; ++i) { - h.send_displs[static_cast(i)] = checked_mpi_count(running_send, "Send displacement"); - running_send += h.send_counts[static_cast(i)]; + for (size_t k = 0; k < h.window.count; ++k) { + const size_t i = h.window.slot(WindowIndex{k}); + const size_t n = (static_cast(i) == self) ? 0 : slot_block(send_data, i).size(); + const int c = checked_mpi_count(n, "Send count"); + h.send_counts[i] = c; + h.send_displs[i] = checked_mpi_count(running_send, "Send displacement"); + running_send += c; } - h.send_buffer.resize(static_cast(checked_mpi_count(total_send, "Total send count"))); - for (int i = 0; i < num_ranks; ++i) { - const int c = h.send_counts[static_cast(i)]; + h.send_buffer.resize(static_cast(checked_mpi_count(running_send, "Total send count"))); + for (size_t k = 0; k < h.window.count; ++k) { + const size_t i = h.window.slot(WindowIndex{k}); + const int c = h.send_counts[i]; if (c == 0) { continue; } - std::copy(send_data[static_cast(i)].begin(), - send_data[static_cast(i)].begin() + c, - h.send_buffer.begin() + h.send_displs[static_cast(i)]); + const auto &block = slot_block(send_data, i); + std::copy(block.begin(), block.begin() + c, h.send_buffer.begin() + h.send_displs[i]); } - h.recv_counts.resize(static_cast(num_ranks)); + h.recv_counts.assign(static_cast(num_ranks), 0); const AlltoallvResolveArgs resolve_args{.send = h.send_buffer.data(), .send_counts = h.send_counts.data(), @@ -209,29 +303,36 @@ inline auto begin_alltoallv(const std::vector> &send_data, } #ifdef monoprop_ENABLE_MPI if (known_recv_counts == nullptr && comm.kind == Comm::Kind::Hybrid) { - comm.hyb->alltoallv_resolve(comm.shm_rank, resolve_args, datatype::get()); + comm.hyb->alltoallv_resolve(comm.shm_rank, resolve_args, datatype::get(), plan); return h; } #endif if (known_recv_counts != nullptr) { - std::copy( - known_recv_counts->begin(), - known_recv_counts->begin() + std::min(known_recv_counts->size(), static_cast(num_ranks)), - h.recv_counts.begin()); + // The caller's array is FLAT [P]; the window is the mask, as alltoall_counts already masks the + // counts it exchanges. No receive is ever posted for a non-peer, so a non-zero count there sizes + // recv_buffer for bytes nothing writes and wait_into hands the caller uninitialised memory. + const size_t avail = known_recv_counts->size(); + for (size_t k = 0; k < h.window.count; ++k) { + const size_t i = h.window.slot(WindowIndex{k}); + if (i < avail) { + h.recv_counts[i] = (*known_recv_counts)[i]; + } + } if (self >= 0) { h.recv_counts[static_cast(self)] = 0; } } else { - alltoall_counts(h.send_counts.data(), h.recv_counts.data(), num_ranks, comm); + alltoall_counts(h.send_counts.data(), h.recv_counts.data(), num_ranks, comm, plan); } // Wide accumulator + checked narrowing: see checked_mpi_count. long long running = 0; - for (int i = 0; i < num_ranks; ++i) { - h.recv_displs[static_cast(i)] = checked_mpi_count(running, "Recv displacement"); - running += h.recv_counts[static_cast(i)]; + for (size_t k = 0; k < h.window.count; ++k) { + const size_t i = h.window.slot(WindowIndex{k}); + h.recv_displs[i] = checked_mpi_count(running, "Recv displacement"); + running += h.recv_counts[i]; } h.recv_buffer.resize(static_cast(checked_mpi_count(running, "Total recv count"))); @@ -255,21 +356,41 @@ inline auto begin_alltoallv(const std::vector> &send_data, } #ifdef monoprop_ENABLE_MPI else if (comm.kind == Comm::Kind::Hybrid) { - comm.hyb->alltoallv(comm.shm_rank, flat, datatype::get()); + comm.hyb->alltoallv(comm.shm_rank, flat, datatype::get(), plan); } #endif else { #ifdef monoprop_ENABLE_MPI - MPI_Ialltoallv(h.send_buffer.data(), - h.send_counts.data(), - h.send_displs.data(), - datatype::get(), - h.recv_buffer.data(), - h.recv_counts.data(), - h.recv_displs.data(), - datatype::get(), - comm.mpi, - &h.request); + if (plan.dense()) { + MPI_Ialltoallv(h.send_buffer.data(), + h.send_counts.data(), + h.send_displs.data(), + datatype::get(), + h.recv_buffer.data(), + h.recv_counts.data(), + h.recv_displs.data(), + datatype::get(), + comm.mpi, + &h.request); + } + else { + // S == 1 world: the same pairing as the Hybrid path, one message per reachable peer, left + // in flight in the handle exactly as MPI_Ialltoallv is. The buffers MPI holds live in `h` + // and travel with it: a vector move keeps its heap block, so returning `h` moves nothing + // MPI is reading. + h.posted = sparse_pairwise(plan, + rank(comm), + num_ranks, + comm.mpi, + kFlatPayloadTag, + datatype::get(), + sizeof(T), + reinterpret_cast(h.send_buffer.data()), + PeerLayout{.counts = h.send_counts.data(), .displs = h.send_displs.data()}, + reinterpret_cast(h.recv_buffer.data()), + PeerLayout{.counts = h.recv_counts.data(), .displs = h.recv_displs.data()}, + h.requests); + } #else h.recv_buffer = h.send_buffer; // single participant: self round-trip (layouts identical) #endif diff --git a/cpp/monoprop/detail/mpi/MPIUtils.h b/cpp/monoprop/detail/mpi/MPIUtils.h index 3a2d20ab..6022d2a6 100644 --- a/cpp/monoprop/detail/mpi/MPIUtils.h +++ b/cpp/monoprop/detail/mpi/MPIUtils.h @@ -16,12 +16,15 @@ #include #include +#include +#include #include #include "monoprop/MPGraph.h" #include "monoprop/TypeAliases.h" #include "monoprop/core/Monomial.h" #include "monoprop/detail/mpi/MPICompat.h" +#include "monoprop/detail/mpi/Routing.h" namespace monoprop::mpi_detail { @@ -49,12 +52,67 @@ inline auto read_monomial_from_words(const VecZ &buffer, size_t start) -> Monomi namespace monoprop { // Stateless and identical on every rank, so all ranks agree on a term's owner without communication. +// Goes through routing::Router::dest and nothing else: this and Scan.h's query emission must return the +// same slot for the same monomial, and a divergence splits ownership silently. There is deliberately no +// rank-count overload -- it would answer splitmix during a linear run, which is exactly that split. template -auto find_rank(const Monomial &mono, const size_t n_ranks) -> size_t { - if (n_ranks == 0) { - return 0; +auto find_rank(const Monomial &mono, const routing::Router &router) -> size_t { + return router.dest(mono); +} + +// The router this communicator's geometry implies, honouring monoprop_ROUTING. +// Templated because the router binds the transposed basis for this monomial width. +template +inline auto router_for(const mpi::Comm &comm) -> routing::Router { + const auto geom = mpi::geometry(comm); + return routing::make_router(static_cast(geom.ranks), static_cast(geom.partitions)); +} + +class RoutingDisagreement : public std::runtime_error { +public: + using std::runtime_error::runtime_error; +}; + +// Every participant must resolve the SAME router, and the failure mode if they do not is a hang, not a +// wrong answer: linear routing makes each rank post receives from the peers its own bits imply, so a +// rank whose monoprop_ROUTING or _ROUTE_SEED did not reach it waits forever on a peer that never sends. +// Turning that into an exception at construction costs two allreduces, called once and never per gate. +// +// TWO independent digests, not one: allreduce_sum is the only collective in the tree, and a sum is not +// an equality test -- differing values can add up to mine*world. Both must agree, so a disagreement +// survives at ~2^-128 rather than ~2^-64. Partitions are in the digest because S enters Router::dest: +// two ranks differing only in S agree on the mode and the seed and still route apart. +inline auto check_routing_agreement(const mpi::Comm &comm) -> void { + const auto world = static_cast(mpi::size(comm)); + if (world <= 1) { + return; + } + const auto geom = mpi::geometry(comm); + const auto parts = static_cast(geom.partitions); + const auto linear = static_cast(routing::linear_requested()); + const uint64_t seed = routing::seed_from_env(); + const auto digest = [&](uint64_t salt) { + return routing::mix64(routing::mix64(routing::mix64(salt ^ linear) ^ parts) ^ seed); + }; + const uint64_t first = digest(0x9E37'79B9'7F4A'7C15ULL); + const uint64_t second = digest(0xC2B2'AE3D'27D4'EB4FULL); + const auto agrees = [&](uint64_t mine) { + return mpi::allreduce_sum(mine, comm) == mine * static_cast(world); + }; + // Both allreduces run on every participant: short-circuiting the second would itself deadlock. + const bool ok_first = agrees(first); + const bool ok_second = agrees(second); + if (!ok_first || !ok_second) { + throw RoutingDisagreement( + std::format("routing configuration differs across the {} participants (this one: linear={}, " + "partitions={}, seed={}). monoprop_ROUTING / monoprop_ROUTE_SEED must reach every rank " + "identically -- under linear routing a disagreement deadlocks the exchange rather than " + "corrupting it.", + world, + linear, + parts, + seed)); } - return monomial_hash(mono) % n_ranks; } } // namespace monoprop diff --git a/cpp/monoprop/detail/mpi/Pairwise.h b/cpp/monoprop/detail/mpi/Pairwise.h new file mode 100644 index 00000000..84b66293 --- /dev/null +++ b/cpp/monoprop/detail/mpi/Pairwise.h @@ -0,0 +1,112 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +#include +#include +#include +#include + +#include + +#include "monoprop/detail/mpi/Comm.h" + +namespace monoprop::mpi { + +// One tag per (transport, verb), all five here so no two can collide unseen: one thread per rank calls +// MPI, so the tag is all that keeps a count round in flight from being matched by a payload receive. +// +// Why Engine.h's run_exchange may post BOTH its begin_alltoallv rounds under kFlatPayloadTag on one +// communicator: MPI does not overtake within a (src, dst, tag, comm); the query round's MPI_Waitall +// completes before round 2 posts; and both ends skip a zero-count leg on the same value -- what one +// sends a peer IS that peer's recv count, by transpose -- so the two posted sequences match element +// for element and a round-2 receive cannot match a round-1 send. +inline constexpr int kHybridCountTag = 0x6D70; // 'mp' +inline constexpr int kHybridPayloadTag = 0x6D71; +inline constexpr int kFlatPayloadTag = 0x6D72; +inline constexpr int kFlatCountTag = 0x6D73; +// Graph REPLAY payload (Exchange.h). Its own value, so a replay leg can match neither the build path's +// counts nor its payload even though both run over the same communicator. +inline constexpr int kFlatReplayTag = 0x6D74; + +// Per-peer element counts and offsets. Null `counts` is the fixed-block case: `block` each, at b*block. +struct PeerLayout { + const int *counts = nullptr; + const int *displs = nullptr; + int block = 0; + + [[nodiscard]] auto count(int b) const -> int { return counts != nullptr ? counts[b] : block; } + [[nodiscard]] auto displ(int b) const -> size_t { + return static_cast(displs != nullptr ? displs[b] : b * block); + } +}; + +// A variable all-to-all as point-to-point over `plan`'s peers: one Irecv/Isend pair each, the self peer +// copied in place. Counts and displacements are in ELEMENTS of `dt`, whose extent must be `elem`. +// `reqs` is caller storage, grown then INDEXED: MPI holds these pointers until the wait, so a +// reallocating push_back would dangle them. +// +// POSTS ONLY, and returns how many of `reqs` are live. The caller waits, so `send`, `recv` and `reqs` +// must all outlive that wait -- which is what lets a caller hold the round open (PendingAlltoallv) the +// same way the dense branch holds an MPI_Ialltoallv. +// +// `active_legs` is an UPPER BOUND on the peers that will post, for a caller that already knows it (a +// dense plan over a mostly-empty layout sizes `reqs` at 2 * n_ranks otherwise); negative means "assume +// every peer posts". Too small an upper bound is caught by the assert below, not silently. +[[nodiscard]] inline auto sparse_pairwise(PeerPlan plan, + int me, + int n_ranks, + MPI_Comm comm, + int tag, + MPI_Datatype dt, + size_t elem, + const std::byte *send, + PeerLayout send_lay, + std::byte *recv, + PeerLayout recv_lay, + std::vector &reqs, + int active_legs = -1) -> int { + const int f = plan.count(n_ranks); + const auto cap = static_cast(2 * (active_legs < 0 || active_legs > f ? f : active_legs)); + if (reqs.size() < cap) { + reqs.resize(cap); + } + int n_req = 0; + for (int k = 0; k < f; ++k) { + const int b = plan.peer(me, k); + const int sc = send_lay.count(b); + const int rc = recv_lay.count(b); + std::byte *rbuf = recv + recv_lay.displ(b) * elem; + const std::byte *sbuf = send + send_lay.displ(b) * elem; + if (b == me) { + // The self slot is a copy, not a message: its two counts are each other's transpose. + assert(sc == rc); + if (rc != 0) { + std::memcpy(rbuf, sbuf, static_cast(rc) * elem); + } + continue; + } + assert(static_cast(n_req) + 2 <= reqs.size() || (rc == 0 && sc == 0)); // active_legs too small + if (rc != 0) { + MPI_Irecv(rbuf, rc, dt, b, tag, comm, &reqs[static_cast(n_req++)]); + } + if (sc != 0) { + MPI_Isend(sbuf, sc, dt, b, tag, comm, &reqs[static_cast(n_req++)]); + } + } + return n_req; +} + +} // namespace monoprop::mpi diff --git a/cpp/monoprop/detail/mpi/Routing.h b/cpp/monoprop/detail/mpi/Routing.h new file mode 100644 index 00000000..c9a2ae04 --- /dev/null +++ b/cpp/monoprop/detail/mpi/Routing.h @@ -0,0 +1,327 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "monoprop/core/Monomial.h" +#include "monoprop/detail/EnvConfig.h" + +// The single home for "which flat slot owns this monomial". Two call sites depend on agreeing exactly +// (Scan.h emits queries by it, MonomialPropagator seeds the operator by it), and a disagreement splits +// ownership silently rather than crashing -- so both go through this Router and nothing else. Scan.h +// enters at dest_from_shift, which is dest with the per-generator rank bits taken from the shift. +// +// Two-level, because the levels cost differently: across MPI ranks the message COUNT is what hurts, so +// the rank index is GF(2)-linear in the support and a generator maps every query to one peer; within a +// rank partitions talk through shared memory, where fanout is free and only balance matters. +// +// part = q % S q = monomial_hash(M) (splitmix, unchanged) +// rank = a & (R - 1) a = linear_hash(M); all log2(R) rank bits, so fanout is 1 +// flat = rank * S + part +// +// Linear or not is a switch and not a dial: the rank takes every bit from the linear hash or none of +// them. None of them is the splitmix router, which is `q % (R*S)` bit for bit, and R == 1 is that case +// by construction. R > 1 must then be a power of two, or there is no XOR structure to route by and the +// geometry is rejected (UnroutableGeometry) rather than silently falling back. +// +// The derivation, what linear routing buys and what it costs: docs/content/docs/features/parallelism.mdx, +// under "Rank routing". +// +// Knobs, parsed and validated in EnvConfig.h: +// monoprop_ROUTING linear (default) | splitmix +// monoprop_ROUTE_SEED uint64 seed for the linear basis (default kDefaultSeed) + +namespace monoprop::routing { + +inline constexpr uint64_t kDefaultSeed = 0x5DEE'CE66'D0C6'2517ULL; + +// A rank count linear routing cannot serve. Thrown at Router construction, not at the first term: the +// alternative is a silent fallback to splitmix on some ranks and a deadlocked exchange. +class UnroutableGeometry : public std::runtime_error { +public: + using std::runtime_error::runtime_error; +}; + +inline constexpr auto mix64(uint64_t x) noexcept -> uint64_t { + x += 0x9E37'79B9'7F4A'7C15ULL; + x = (x ^ (x >> 30)) * 0xBF58'476D'1CE4'E5B9ULL; + x = (x ^ (x >> 27)) * 0x94D0'49BB'1331'11EBULL; + return x ^ (x >> 31); +} + +inline auto seed_from_env() -> uint64_t { + return config::get().route_seed.value_or(kDefaultSeed); +} + +// One 64-bit vector per Majorana mode. Deterministic from the seed alone, so every rank builds the +// same table with no communication -- the property find_rank's contract rests on. +template +inline auto linear_basis() -> const std::array & { + static const auto table = [] { + std::array v{}; + const uint64_t seed = seed_from_env(); + for (size_t i = 0; i < NumBits; ++i) { + v[i] = mix64(mix64(seed) + (static_cast(i) * 0x9E37'79B9'7F4A'7C15ULL)); + } + return v; + }(); + return table; +} + +// The full 64-bit image, by walking the set bits. Router::dest does NOT use this -- it needs only the +// low log2(R) bits and takes the transposed form below -- but the map is defined here, and the +// GF(2)-linearity test and the plane build both read it as the definition. +template +[[nodiscard]] inline auto linear_hash(const monoprop::Bitset &bits) noexcept -> uint64_t { + const auto &v = linear_basis(); + uint64_t h = 0; + for (size_t i = bits.find_first(); i < NumBits; i = bits.find_next(i)) { + h ^= v[i]; + } + return h; +} + +// One per output bit; a Router reads log2(R) of them. Not trimmed to that: the count is geometry, and +// keying the table on it would build one table per Router instead of one per (seed, width). +inline constexpr size_t kLinearPlanes = 64; + +template +inline constexpr size_t kPlaneWords = monoprop::Bitset::num_words(); + +// The basis transposed: plane j, bit i, is bit j of basis vector i. Then bit j of linear_hash(M) is +// popcount(M & plane_j) & 1 -- log2(R) * words branch-free AND/XOR/popcount ops instead of a gather +// whose length is the term's popcount (~20-28 under the production cutoff). Same map, bit for bit. +// +// Keyed on the seed and the width alone, never on the geometry, so one table serves every Router; the +// Router binds a pointer to it at construction and dest() never reaches the static-init guard. +template +inline auto linear_planes() -> const std::array> & { + static const auto table = [] { + std::array> planes{}; + const auto &v = linear_basis(); + for (size_t i = 0; i < NumBits; ++i) { + const size_t word = i / 64; + const uint64_t bit = uint64_t{1} << (i % 64); + for (size_t j = 0; j < kLinearPlanes; ++j) { + if (((v[i] >> j) & 1U) != 0) { + planes[(j * kPlaneWords)+word] |= bit; + } + } + } + return planes; + }(); + return table; +} + +// GF(2) rank of a set of 64-bit vectors, by Gaussian elimination over the bit columns. The per-generator +// rank shifts must span at least log2(R) dimensions or the reachable destination ranks form a strict +// subspace and R - 2^rank ranks stay empty. A balance failure, not a correctness one, so its caller +// (MonomialPropagator::report_routing_coverage_) warns on stderr rather than gating. +// Measured: rank 32 for the 60-site Hubbard's 416 distinct shifts, against the 7 bits R = 128 needs. +[[nodiscard]] inline auto gf2_rank(std::vector vectors) noexcept -> size_t { + size_t rank = 0; + for (size_t bit = 0; bit < 64 && rank < vectors.size(); ++bit) { + const uint64_t probe = uint64_t{1} << bit; + size_t pivot = vectors.size(); + for (size_t i = rank; i < vectors.size(); ++i) { + if ((vectors[i] & probe) != 0) { + pivot = i; + break; + } + } + if (pivot == vectors.size()) { + continue; + } + std::swap(vectors[rank], vectors[pivot]); + for (size_t i = 0; i < vectors.size(); ++i) { + if (i != rank && (vectors[i] & probe) != 0) { + vectors[i] ^= vectors[rank]; + } + } + ++rank; + } + return rank; +} + +// Trivially copyable and cheap to build; hold one per build_layer call rather than per term. +class Router final { +public: + // Bound to a monomial width: dest() reads the transposed basis for THAT width, and binding the + // pointer here is what keeps the static-init guard out of the per-term path. The only way to a + // linear router, so an unbound one cannot reach dest(). Throws if `linear` and ranks is not 2^k. + template + [[nodiscard]] static auto for_modes(size_t ranks, size_t partitions, bool linear) -> Router { + Router r{ranks, partitions, linear}; + r.planes_ = linear_planes<2 * NumModes>().data(); + r.plane_words_ = kPlaneWords<2 * NumModes>; + return r; + } + + // Today's routing: one flat world, full avalanche, no linear bits. Width-free: no plane is read. + static constexpr auto splitmix(size_t flat_world) noexcept -> Router { return Router{flat_world, 1, false}; } + + [[nodiscard]] constexpr auto ranks() const noexcept -> size_t { return ranks_; } + [[nodiscard]] constexpr auto partitions() const noexcept -> size_t { return parts_; } + [[nodiscard]] constexpr auto flat_world() const noexcept -> size_t { return flat_; } + // False for a splitmix router AND for R == 1, which has no rank bit to take: both route densely. + [[nodiscard]] constexpr auto is_linear() const noexcept -> bool { return linear_; } + // Rank bits read off the linear hash: log2(R), or 0 when not linear. The span a generator set must + // cover for every rank to be reachable. + [[nodiscard]] constexpr auto linear_bits() const noexcept -> size_t { + return linear_ ? static_cast(std::countr_zero(ranks_)) : 0; + } + + // The same number for a geometry alone, with no monomial width bound: it IS the constructor, so the + // resolution and the non-power-of-two throw cannot drift from the router's. + [[nodiscard]] static auto bits_for(size_t ranks, bool linear) -> size_t { + return Router{ranks, 1, linear}.linear_bits(); + } + + // Flat destination slot in [0, flat_world). Branch is on a member, so it is perfectly predicted. + template + [[nodiscard]] [[gnu::always_inline]] inline auto dest(const Monomial &mono) const noexcept -> size_t { + if (!linear_) { + return static_cast(monomial_hash(mono) % flat_); // bit-for-bit today's `hash % P` + } + const size_t rank = static_cast(linear_low_(mono)); + if (parts_ == 1) { + return rank; // q % 1 == 0 for every q, so S == 1 owes the hash nothing + } + return (rank * parts_) + part_of_(monomial_hash(mono)); + } + + // The emit path's dest, for a query M^G raised on a term M THIS slot owns. rank(M^G) == rank(M) ^ + // shift(G) is exact under linear routing, so the planes are a per-generator constant already in hand + // and only the partition index is per term. `my_flat` must be this slot's own index and `shift` this + // generator's rank_shift, or ownership moves silently -- the same precondition mpi::PeerPlan carries. + // Splitmix has no such identity and falls back to dest(), bit for bit. + template + [[nodiscard]] [[gnu::always_inline]] inline auto dest_from_shift(const Monomial &mono, + size_t my_flat, + size_t shift) const noexcept -> size_t { + if (!linear_) { + return dest(mono); + } + const size_t rank = rank_of_slot_(my_flat) ^ shift; + if (parts_ == 1) { + return rank; + } + return (rank * parts_) + part_of_(monomial_hash(mono)); + } + + // The rank-level shift a generator induces: rank(M^G) == rank(M) ^ shift(G). Zero for every G when + // the router is not linear. This is what makes the destination predictable. + template + [[nodiscard]] auto rank_shift(const Monomial &gen) const noexcept -> size_t { + return static_cast(linear_low_(gen)); + } + +private: + // ranks x partitions == the flat world the destinations index. Private: a linear router must go + // through for_modes, which binds the basis linear_low_ reads. + constexpr Router(size_t ranks, size_t partitions, bool linear) + : ranks_(ranks == 0 ? 1 : ranks), + parts_(partitions == 0 ? 1 : partitions), + flat_(ranks_ * parts_), + linear_(linear && ranks_ > 1), // R == 1 takes no rank bit, so it IS the dense case + parts_pow2_(std::has_single_bit(parts_)), + parts_mask_(parts_ - 1), + parts_log2_(static_cast(std::countr_zero(parts_))) { + if (linear && !std::has_single_bit(ranks_)) { + throw UnroutableGeometry( + std::format("linear routing needs a power-of-two rank count, got {}. Launch 2^k ranks, or set " + "monoprop_ROUTING=splitmix to keep the dense all-to-all.", + ranks_)); + } + } + + // q % S. Every production layout has S in {1,2,4,8,16}, where the modulo is a mask -- and parts_ is a + // runtime member, so the compiler cannot strength-reduce it on our behalf. Identical bit for bit. + [[nodiscard]] [[gnu::always_inline]] inline auto part_of_(uint64_t q) const noexcept -> size_t { + if (parts_pow2_) { + assert(static_cast(q & parts_mask_) == static_cast(q % parts_)); + return static_cast(q & parts_mask_); + } + return static_cast(q % parts_); + } + + // Flat slot -> rank index. Flat slots are rank * S + partition (mpi::rank under the hybrid comm). + [[nodiscard]] [[gnu::always_inline]] inline auto rank_of_slot_(size_t flat_slot) const noexcept -> size_t { + if (parts_pow2_) { + assert((flat_slot >> parts_log2_) == flat_slot / parts_); + return flat_slot >> parts_log2_; + } + return flat_slot / parts_; + } + + // linear_hash(M) & (R - 1), one output bit per plane: parity(popcount(M & plane_j)). Folding the + // words with XOR before the popcount is the same parity (popcount(x)+popcount(y) == popcount(x^y) + // mod 2) for one popcount per bit instead of one per word. A non-linear router reads no plane. + template + [[nodiscard]] [[gnu::always_inline]] inline auto linear_low_(const Monomial &m) const noexcept + -> uint64_t { + constexpr size_t kW = kPlaneWords<2 * NumModes>; + const size_t bits = linear_bits(); + assert(bits == 0 || (planes_ != nullptr && plane_words_ == kW)); // bound at a different width + uint64_t acc = 0; + for (size_t j = 0; j < bits; ++j) { + const uint64_t *plane = planes_ + (j * kW); + uint64_t fold = 0; + for (size_t w = 0; w < kW; ++w) { + fold ^= m.word(w) & plane[w]; + } + acc |= static_cast(std::popcount(fold) & 1) << j; + } + return acc; + } + + size_t ranks_; + size_t parts_; + size_t flat_; + bool linear_; + bool parts_pow2_; // S is 2^k, so `% S` is a mask and `/ S` a shift + size_t parts_mask_; // S - 1, and parts_log2_ == log2(S); both meaningless unless parts_pow2_ + size_t parts_log2_; + const uint64_t *planes_ = nullptr; // [kLinearPlanes x plane_words_], owned by linear_planes() + size_t plane_words_ = 0; +}; + +// The mode, before any geometry. Linear unless asked otherwise: measured at the production point it +// costs nothing on balance (rank occupancy max/mean 1.001 at R=128, all ranks used) and takes messages +// per rank per layer from 362,712 to 1,397, i.e. from proportional-to-R to flat. +inline auto linear_requested() -> bool { + return config::get().routing_mode.value_or(config::RoutingMode::Linear) == config::RoutingMode::Linear; +} + +// Resolved rank bits for a geometry, without a router: the replay transport gates on the number. +inline auto linear_bits_for(size_t ranks) -> size_t { + return Router::bits_for(ranks, linear_requested()); +} + +template +inline auto make_router(size_t ranks, size_t partitions) -> Router { + return Router::for_modes(ranks, partitions, linear_requested()); +} + +} // namespace monoprop::routing diff --git a/cpp/tests/README.md b/cpp/tests/README.md index 6cd35461..dc4c1e2b 100644 --- a/cpp/tests/README.md +++ b/cpp/tests/README.md @@ -87,8 +87,10 @@ name and cannot address suite-nested cases, tests use flat `pauli_algebra_tests.cpp`, `majorana_cutoff_tests.cpp` (length/support cutoff, CutoffEvaluator, interleave phase, coeff encode/decode), `validation_tests.cpp` (parameter validators), `mpi_utils_tests.cpp` (find_rank, word serialization, - scan routing agreement), `evolution_detail_tests.cpp` (MatchedEpochSet + - CutoffContext), + scan routing agreement), `routing_tests.cpp` (Router term -> flat-slot map: + the `d = 0` splitmix equivalence, the linear shift identity, the gf2_rank + coverage diagnostic, and the shipped default), + `evolution_detail_tests.cpp` (MatchedEpochSet + CutoffContext), `row_accessor_tests.cpp` (dense vs OperatorIndex row accessors), `sparse_monomial_tests.cpp` (the `(k, d)` cutoff predicates vs their bitset forms). diff --git a/cpp/tests/boost-test.cmake b/cpp/tests/boost-test.cmake index 7d0d52c9..19bf4828 100644 --- a/cpp/tests/boost-test.cmake +++ b/cpp/tests/boost-test.cmake @@ -1,6 +1,8 @@ +# 4 as well as 2: at R = 2 the peer plan can only ever resolve one peer, so every `for k in [0, f)` in +# the sparse transport stays a single iteration and the f > 1 cases never run. set( monoprop_MPI_TEST_PROCS - "2" + "2;4" CACHE STRING "Semicolon-separated list of ranks for MPI test variants" ) diff --git a/cpp/tests/env_config_tests.cpp b/cpp/tests/env_config_tests.cpp index fe46afe7..9d7b80e7 100644 --- a/cpp/tests/env_config_tests.cpp +++ b/cpp/tests/env_config_tests.cpp @@ -14,11 +14,16 @@ #include +#include #include #include "monoprop/detail/EnvConfig.h" +using monoprop::config::EnvConfigError; +using monoprop::config::RoutingMode; using monoprop::config::detail::parse_positive_int; +using monoprop::config::detail::parse_routing_mode; +using monoprop::config::detail::parse_uint64; BOOST_AUTO_TEST_CASE(env_config_parse_positive_int_null_and_malformed) { BOOST_CHECK(parse_positive_int(nullptr) == std::nullopt); @@ -44,3 +49,26 @@ BOOST_AUTO_TEST_CASE(env_config_settings_cached_singleton) { // Touch a field so the Settings aggregate is actually read. BOOST_CHECK(a.num_threads == std::nullopt || *a.num_threads >= 1); } + +// Both routing parsers throw where parse_positive_int returns nullopt: a routing knob that +// defaulted silently would change the transport with no diagnostic. +BOOST_AUTO_TEST_CASE(env_config_parse_uint64_unset_valid_and_rejected) { + BOOST_CHECK(parse_uint64("k", nullptr) == std::nullopt); + BOOST_CHECK(parse_uint64("k", "") == std::nullopt); + BOOST_CHECK(parse_uint64("k", "0") == std::optional(0)); + BOOST_CHECK(parse_uint64("k", "18446744073709551615") == std::optional(~std::uint64_t{0})); + BOOST_CHECK_THROW(parse_uint64("k", "abc"), EnvConfigError); + BOOST_CHECK_THROW(parse_uint64("k", "12x"), EnvConfigError); + BOOST_CHECK_THROW(parse_uint64("k", "-1"), EnvConfigError); // strtoull would WRAP it + BOOST_CHECK_THROW(parse_uint64("k", "18446744073709551616"), EnvConfigError); // ERANGE +} + +BOOST_AUTO_TEST_CASE(env_config_parse_routing_mode_rejects_a_typo) { + BOOST_CHECK(parse_routing_mode("k", nullptr) == std::nullopt); + BOOST_CHECK(parse_routing_mode("k", "") == std::nullopt); + BOOST_CHECK(parse_routing_mode("k", "splitmix") == std::optional(RoutingMode::Splitmix)); + BOOST_CHECK(parse_routing_mode("k", "linear") == std::optional(RoutingMode::Linear)); + BOOST_CHECK_THROW(parse_routing_mode("k", "dense"), EnvConfigError); // all three used to mean linear + BOOST_CHECK_THROW(parse_routing_mode("k", "off"), EnvConfigError); + BOOST_CHECK_THROW(parse_routing_mode("k", "Linear"), EnvConfigError); +} diff --git a/cpp/tests/evolution_detail_tests.cpp b/cpp/tests/evolution_detail_tests.cpp index bade37f4..1cc0c90f 100644 --- a/cpp/tests/evolution_detail_tests.cpp +++ b/cpp/tests/evolution_detail_tests.cpp @@ -165,7 +165,7 @@ BOOST_AUTO_TEST_CASE(self_resolve_mark_bounded_by_combined_size) { }; stage_self(terms[1], 1); stage_self(terms[5], -1); - eng.src_idx_r[0] = {0, 2}; + eng.src_idx_r.at_slot(0) = {0, 2}; eng.resolve_self_queries(/*is_leader_pass=*/true); diff --git a/cpp/tests/flat_exchange_tests.cpp b/cpp/tests/flat_exchange_tests.cpp new file mode 100644 index 00000000..7621e315 --- /dev/null +++ b/cpp/tests/flat_exchange_tests.cpp @@ -0,0 +1,230 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// post_flat_alltoallv's transport choice (the graph REPLAY path). `wire_bits` alone picks it -- never the +// layout -- because a data-dependent branch straddles: a rank inside MPI_Ialltoallv waits forever on one +// that chose point-to-point. Ticket::in_flight() is what tells the two apart: 1 for the collective, two +// requests per posted leg, 0 for a round with nothing to move. + +#include + +#include +#include +#include +#include + +#include "monoprop/detail/mpi/Comm.h" +#include "monoprop/detail/mpi/Exchange.h" + +#ifdef monoprop_ENABLE_MPI +#include + +#include +#endif + +using monoprop::mpi::active_leg_count; +using monoprop::mpi::Comm; +using monoprop::mpi::post_flat_alltoallv; + +namespace { + +// The prefix sum post_flat_alltoallv takes on both sides of a symmetric layout. +auto displs_of(const std::vector &counts) -> std::vector { + std::vector displs(counts.size()); + int running = 0; + for (size_t i = 0; i < counts.size(); ++i) { + displs[i] = running; + running += counts[i]; + } + return displs; +} + +auto total_of(const std::vector &counts) -> int { + return std::accumulate(counts.begin(), counts.end(), 0); +} + +#ifdef monoprop_ENABLE_MPI +auto world_size() -> int { + int n = 0; + MPI_Comm_size(MPI_COMM_WORLD, &n); + return n; +} +auto world_rank() -> int { + int r = 0; + MPI_Comm_rank(MPI_COMM_WORLD, &r); + return r; +} +// What Evolution's gate resolves to when routing gives fanout 1, which is the default at a power-of-two +// rank count. Tests must not derive it from their own layout -- that is the straddle. +auto wire_bits_for(int n) -> int { + return std::countr_zero(static_cast(n)); +} +#endif + +} // namespace + +// The upper bound the request vector is sized from: a leg counts if EITHER side carries a payload, so a +// bound too small can never be handed to sparse_pairwise. +BOOST_AUTO_TEST_CASE(flat_exchange_active_legs_take_either_side) { + constexpr int n = 8; + std::vector send(static_cast(n), 0); + std::vector recv(static_cast(n), 0); + BOOST_CHECK_EQUAL(active_leg_count(send.data(), recv.data(), n), 0); + send[1] = 4; + recv[6] = 4; + BOOST_CHECK_EQUAL(active_leg_count(send.data(), recv.data(), n), 2); + recv[1] = 4; // same leg, both sides + BOOST_CHECK_EQUAL(active_leg_count(send.data(), recv.data(), n), 2); +} + +// The self leg is a copy, not a message, so a one-rank world posts nothing at all and wait() is a no-op. +// Same on the non-MPI build, where the fallback self-copy runs instead. +BOOST_AUTO_TEST_CASE(flat_exchange_self_leg_is_copied_with_nothing_posted) { + Comm c{MPI_COMM_SELF}; + const std::vector counts{3}; + const auto displs = displs_of(counts); + const std::vector send{7, 8, 9}; + std::vector out(static_cast(total_of(counts)), -1); + auto ticket = post_flat_alltoallv({.send = send.data(), + .send_counts = counts.data(), + .send_displs = displs.data(), + .recv = out.data(), + .recv_counts = counts.data(), + .recv_displs = displs.data()}, + 1, + c, + /*wire_bits=*/1); + BOOST_CHECK_EQUAL(ticket.in_flight(), 0); + ticket.wait(); + BOOST_CHECK(out == send); +} + +#ifdef monoprop_ENABLE_MPI + +namespace { + +// One round over a symmetric layout, returning what wait() had to drain. `counts` is both sides. +auto run_round(int n, + int wire_bits, + const std::vector &counts, + const std::vector &send, + std::vector &out) -> int { + Comm c{MPI_COMM_WORLD}; + const auto displs = displs_of(counts); + auto ticket = post_flat_alltoallv({.send = send.data(), + .send_counts = counts.data(), + .send_displs = displs.data(), + .recv = out.data(), + .recv_counts = counts.data(), + .recv_displs = displs.data()}, + n, + c, + wire_bits); + const int drained = ticket.in_flight(); + ticket.wait(); + return drained; +} + +} // namespace + +// wire_bits == 0 is today's collective whatever the layout holds -- including the sparse layout the +// pairwise arm exists for. The layout must not be able to move the branch. +BOOST_AUTO_TEST_CASE(flat_exchange_zero_wire_bits_always_takes_the_collective) { + const int n = world_size(); + if (n < 2 || (n % 2) != 0) { + return; + } + const int me = world_rank(); + const int peer = me ^ 1; + for (const bool dense : {true, false}) { + std::vector counts(static_cast(n), 0); + if (dense) { + std::fill(counts.begin(), counts.end(), 1); + } + else { + counts[static_cast(peer)] = 1; // one leg: the shape the pairwise arm is for + } + const auto displs = displs_of(counts); + std::vector send(static_cast(total_of(counts)), me); + std::vector out(static_cast(total_of(counts)), -1); + BOOST_CHECK_EQUAL(run_round(n, 0, counts, send, out), 1); + for (int i = 0; i < n; ++i) { + if (counts[static_cast(i)] != 0) { + BOOST_CHECK_EQUAL(out[static_cast(displs[static_cast(i)])], i); + } + } + } +} + +// One peer, me ^ 1 -- an involution, so the count matrix stays symmetric and both ends drop the same +// n - 1 legs. Two requests (one Irecv, one Isend) and the same bytes delivered. +BOOST_AUTO_TEST_CASE(flat_exchange_single_leg_takes_the_pairwise_path) { + const int n = world_size(); + if (n < 2 || (n % 2) != 0) { + return; + } + const int me = world_rank(); + const int peer = me ^ 1; + constexpr int len = 3; + std::vector counts(static_cast(n), 0); + counts[static_cast(peer)] = len; + std::vector send(static_cast(len)); + for (int j = 0; j < len; ++j) { + send[static_cast(j)] = (me * 1000) + j; + } + std::vector out(static_cast(len), -1); + BOOST_CHECK_EQUAL(run_round(n, wire_bits_for(n), counts, send, out), 2); + for (int j = 0; j < len; ++j) { + BOOST_CHECK_EQUAL(out[static_cast(j)], (peer * 1000) + j); + } +} + +// Every leg active on the pairwise arm: it posts a pair per non-zero leg rather than falling back, so +// the transport really is the gate's choice and not the layout's. The self leg is a copy, not a pair. +BOOST_AUTO_TEST_CASE(flat_exchange_dense_layout_on_the_pairwise_arm) { + const int n = world_size(); + if (n < 2 || (n & (n - 1)) != 0) { + return; // off a power of two the gate resolves to 0 bits, i.e. the collective + } + const int me = world_rank(); + const std::vector counts(static_cast(n), 1); + const auto displs = displs_of(counts); + std::vector send(static_cast(n)); + for (int d = 0; d < n; ++d) { + send[static_cast(d)] = (me * 1000) + d; + } + std::vector out(static_cast(n), -1); + BOOST_CHECK_EQUAL(run_round(n, wire_bits_for(n), counts, send, out), 2 * (n - 1)); + for (int src = 0; src < n; ++src) { + BOOST_CHECK_EQUAL(out[static_cast(src)], (src * 1000) + me); + } +} + +// No active leg anywhere: nothing is posted, wait() drains nothing, and the recv buffer is left as the +// caller sized it. A layer with no cross-rank partners takes this shape on EVERY rank at once -- the +// symmetric layout is what makes that true -- so the pairwise arm is still the agreed one. +BOOST_AUTO_TEST_CASE(flat_exchange_empty_layout_posts_nothing) { + const int n = world_size(); + if (n < 2 || (n & (n - 1)) != 0) { + return; + } + const std::vector counts(static_cast(n), 0); + BOOST_REQUIRE_EQUAL(total_of(counts), 0); + const std::vector send(1, 0); + std::vector out(1, -1); // Evolution sizes an empty round to 1, not 0 + BOOST_CHECK_EQUAL(run_round(n, wire_bits_for(n), counts, send, out), 0); + BOOST_CHECK_EQUAL(out[0], -1); +} + +#endif // monoprop_ENABLE_MPI diff --git a/cpp/tests/hybrid_comm_tests.cpp b/cpp/tests/hybrid_comm_tests.cpp index bb6ca426..f69ef203 100644 --- a/cpp/tests/hybrid_comm_tests.cpp +++ b/cpp/tests/hybrid_comm_tests.cpp @@ -23,6 +23,8 @@ #ifdef monoprop_ENABLE_MPI #include +#include +#include #include #include #include @@ -204,8 +206,8 @@ BOOST_AUTO_TEST_CASE(hybrid_comm_repeated_alltoallv_varying_sizes) { BOOST_CHECK_EQUAL(failures.load(), 0); } -// alltoallv_resolve driven directly: it folds the count MPI_Alltoall into the payload verb's B1→B2 -// window and sizes recv itself. +// alltoallv_resolve driven directly: it folds the count round into the payload verb's barriered +// windows and sizes recv itself. BOOST_AUTO_TEST_CASE(hybrid_comm_alltoallv_resolve_fused) { if (world_size() < 2) { return; @@ -484,4 +486,670 @@ BOOST_AUTO_TEST_CASE(hybrid_comm_poison_releases_waiters) { } } +// The pairing is derived independently on both ends, so it must be an involution: whoever I send to +// sends back to me. A plan on which it fails deadlocks rather than mis-delivers, and no transport case +// below can distinguish the two. +BOOST_AUTO_TEST_CASE(hybrid_comm_sparse_plan_peer_is_an_involution) { + for (const int shift : {0, 1, 2, 3, 5, 8, 13, 255}) { + const monoprop::mpi::PeerPlan plan{.sparse = true, .shift = shift}; + BOOST_REQUIRE(!plan.dense()); + BOOST_REQUIRE_EQUAL(plan.count(256), 1); + for (int me = 0; me < 256; ++me) { + const int peer = plan.peer(me, 0); + BOOST_REQUIRE_EQUAL(plan.peer(peer, 0), me); + BOOST_REQUIRE(plan.contains(me, peer)); + BOOST_REQUIRE(plan.contains(peer, me)); + BOOST_REQUIRE(!plan.contains(me, peer ^ 1)); // the peer set is a singleton + } + } + const monoprop::mpi::PeerPlan dense_plan{}; + BOOST_REQUIRE(dense_plan.dense()); + BOOST_REQUIRE_EQUAL(dense_plan.count(7), 7); + for (int k = 0; k < 7; ++k) { + BOOST_REQUIRE_EQUAL(dense_plan.peer(3, k), k); // dense ignores `me`: every rank is a peer + BOOST_REQUIRE(dense_plan.contains(3, k)); + } +} + +// A sparse PeerPlan replaces the collectives with point-to-point over the one peer the plan names, so +// the two failure modes it can have are DROPPED data and a HANG -- neither of which a dense-path test +// can see. Every rank derives the same pairing from the same shift, and a block whose destination is +// not the peer must be empty: send only to the plan's peer and check the delivery is exactly that. +BOOST_AUTO_TEST_CASE(hybrid_comm_sparse_plan_delivers_only_to_its_peers) { + const int R = world_size(); + if (R < 2 || (R & (R - 1)) != 0) { + return; // the XOR pairing needs a power-of-two rank count + } + for (const int S : {1, 2, 3}) { + const int P = R * S; + for (int shift = 0; shift < R; ++shift) { + const monoprop::mpi::PeerPlan plan{.sparse = true, .shift = shift}; + const int peer = plan.peer(world_rank(), 0); + BOOST_REQUIRE_EQUAL(plan.count(R), 1); // sparse is always pairwise + std::vector>> recv(static_cast(S)); + auto errs = run_hybrid(S, [&](HybridComm &hyb, int u) { + Comm c = Comm::make_hybrid(&hyb, u); + const int g = monoprop::mpi::rank(c); + std::vector> send(static_cast(P)); + for (int t = 0; t < S; ++t) { + auto &blk = send[static_cast(peer * S + t)]; + for (int j = 0; j <= t; ++j) { + blk.push_back(g * 1000 + t * 10 + j); + } + } + auto h = monoprop::mpi::begin_alltoallv(send, + c, + /*skip_self=*/false, + /*known_recv_counts=*/nullptr, + plan); + std::vector> out; + h.wait_into(out); + recv[static_cast(u)] = out; + }); + for (const auto &e : errs) { + BOOST_CHECK(e == nullptr); + } + // XOR is an involution, so whoever I send to sends to me: my sources are exactly `peer`'s + // partitions, and the block from (peer, su) to my partition t has t+1 entries. + for (int t = 0; t < S; ++t) { + const auto &out = recv[static_cast(t)]; + BOOST_REQUIRE_EQUAL(static_cast(out.size()), P); + for (int src = 0; src < P; ++src) { + const auto &blk = out[static_cast(src)]; + if (src / S != peer) { + BOOST_CHECK(blk.empty()); // a non-peer must not appear at all + continue; + } + BOOST_REQUIRE_EQUAL(static_cast(blk.size()), t + 1); + for (int j = 0; j <= t; ++j) { + BOOST_CHECK_EQUAL(blk[static_cast(j)], src * 1000 + t * 10 + j); + } + } + } + } + } +} + +// Same contract on the S == 1 world, which takes the pure-MPI branch (MPI_Ialltoallv vs Isend/Irecv) +// rather than HybridComm's staged one, and on the KNOWN-recv-counts path the response round uses. +BOOST_AUTO_TEST_CASE(hybrid_comm_sparse_plan_on_the_plain_mpi_path) { + const int R = world_size(); + if (R < 2 || (R & (R - 1)) != 0) { + return; + } + Comm c{MPI_COMM_WORLD}; + for (int shift = 0; shift < R; ++shift) { + const monoprop::mpi::PeerPlan plan{.sparse = true, .shift = shift}; + const int peer = plan.peer(world_rank(), 0); + std::vector> send(static_cast(R)); + for (int j = 0; j < 4; ++j) { + send[static_cast(peer)].push_back(world_rank() * 1000 + j); + } + // Unknown recv layout: the counts round is point-to-point too. + std::vector> out; + monoprop::mpi::begin_alltoallv(send, c, false, nullptr, plan).wait_into(out); + BOOST_REQUIRE_EQUAL(static_cast(out.size()), R); + for (int src = 0; src < R; ++src) { + if (src != peer) { + BOOST_CHECK(out[static_cast(src)].empty()); + continue; + } + BOOST_REQUIRE_EQUAL(static_cast(out[static_cast(src)].size()), 4); + for (int j = 0; j < 4; ++j) { + BOOST_CHECK_EQUAL(out[static_cast(src)][static_cast(j)], src * 1000 + j); + } + } + // Known recv layout (the response round): counts are the transpose, so peer-only again. + std::vector known(static_cast(R), 0); + known[static_cast(peer)] = 4; + std::vector> out2; + monoprop::mpi::begin_alltoallv(send, c, false, &known, plan).wait_into(out2); + BOOST_REQUIRE_EQUAL(static_cast(out2.size()), R); + BOOST_CHECK(out2[static_cast(peer)] == out[static_cast(peer)]); + } +} + +// known_recv_counts is CALLER-supplied, so it can carry a count for a rank the plan does not name -- +// the response round's transpose is only as masked as whatever produced it. No receive is ever posted +// for a non-peer, so an unmasked count sizes recv_buffer for bytes nothing writes and wait_into would +// hand that slot to the caller as data. Both the plain-MPI and the staged HybridComm path must drop it. +BOOST_AUTO_TEST_CASE(hybrid_comm_known_recv_counts_are_masked_through_the_plan) { + const int R = world_size(); + if (R < 2 || (R & (R - 1)) != 0) { + return; + } + constexpr int kReal = 4; + constexpr int kBogus = 7; // what a stale or unmasked transpose would claim a non-peer is sending + for (int shift = 0; shift < R; ++shift) { + const monoprop::mpi::PeerPlan plan{.sparse = true, .shift = shift}; + const int peer = plan.peer(world_rank(), 0); + const int bad = (peer + 1) % R; // the peer set is exactly {peer} + BOOST_REQUIRE(bad != peer); + + // S == 1: the plain-MPI Isend/Irecv branch of begin_alltoallv. + { + Comm c{MPI_COMM_WORLD}; + std::vector> send(static_cast(R)); + for (int j = 0; j < kReal; ++j) { + send[static_cast(peer)].push_back(world_rank() * 1000 + j); + } + std::vector known(static_cast(R), 0); + known[static_cast(peer)] = kReal; + known[static_cast(bad)] = kBogus; + std::vector> out; + monoprop::mpi::begin_alltoallv(send, c, false, &known, plan).wait_into(out); + BOOST_REQUIRE_EQUAL(static_cast(out.size()), R); + BOOST_CHECK(out[static_cast(bad)].empty()); // unmasked, this holds kBogus elements + BOOST_REQUIRE_EQUAL(static_cast(out[static_cast(peer)].size()), kReal); + for (int j = 0; j < kReal; ++j) { + BOOST_CHECK_EQUAL(out[static_cast(peer)][static_cast(j)], peer * 1000 + j); + } + } + + // S == 2: the same array through HybridComm's staged alltoallv. + constexpr int S = 2; + const int P = R * S; + std::vector>> recv(static_cast(S)); + auto errs = run_hybrid(S, [&](HybridComm &hyb, int u) { + Comm c = Comm::make_hybrid(&hyb, u); + const int g = monoprop::mpi::rank(c); + std::vector> send(static_cast(P)); + std::vector known(static_cast(P), 0); + for (int t = 0; t < S; ++t) { + for (int j = 0; j < kReal; ++j) { + send[static_cast((peer * S) + t)].push_back((g * 1000) + j); + } + known[static_cast((peer * S) + t)] = kReal; + known[static_cast((bad * S) + t)] = kBogus; + } + std::vector> out; + monoprop::mpi::begin_alltoallv(send, c, false, &known, plan).wait_into(out); + recv[static_cast(u)] = out; + }); + for (const auto &e : errs) { + BOOST_CHECK(e == nullptr); + } + for (int t = 0; t < S; ++t) { + const auto &out = recv[static_cast(t)]; + BOOST_REQUIRE_EQUAL(static_cast(out.size()), P); + for (int su = 0; su < S; ++su) { + BOOST_CHECK(out[static_cast((bad * S) + su)].empty()); + const auto &blk = out[static_cast((peer * S) + su)]; + BOOST_REQUIRE_EQUAL(static_cast(blk.size()), kReal); + for (int j = 0; j < kReal; ++j) { + BOOST_CHECK_EQUAL(blk[static_cast(j)], (((peer * S) + su) * 1000) + j); + } + } + } + } +} + +namespace { + +// Varies along BOTH ends and hits 0, so a block landing on the wrong peer or the wrong partition +// changes a length, not just a value. +auto sparse_count(int src, int dst) -> int { + return ((src * 3) + (dst * 5)) % 4; +} +auto sparse_tag(int src, int dst, int j) -> int { + return (((src * 128) + dst) * 1000) + j; +} + +} // namespace + +// A zero-count leg is where a send/recv posting asymmetry deadlocks rather than mis-delivers: both ends +// must skip on the SAME value. Nothing above ever sends an empty block over a real message, so force +// one -- the lower-numbered end of every pair sends nothing while its peer sends four. +BOOST_AUTO_TEST_CASE(hybrid_comm_sparse_plan_with_an_empty_leg) { + const int R = world_size(); + if (R < 2 || (R & (R - 1)) != 0) { + return; + } + const int me = world_rank(); + constexpr int kLen = 4; + for (int shift = 1; shift < R; ++shift) { // shift 0 is the self peer, covered separately + const monoprop::mpi::PeerPlan plan{.sparse = true, .shift = shift}; + const int peer = plan.peer(me, 0); + BOOST_REQUIRE(peer != me); + const int my_len = me < peer ? 0 : kLen; // exactly one end of the pair is silent + const int peer_len = peer < me ? 0 : kLen; + + { + Comm c{MPI_COMM_WORLD}; + std::vector> send(static_cast(R)); + for (int j = 0; j < my_len; ++j) { + send[static_cast(peer)].push_back(sparse_tag(me, peer, j)); + } + std::vector> out; + monoprop::mpi::begin_alltoallv(send, c, false, nullptr, plan).wait_into(out); + BOOST_REQUIRE_EQUAL(static_cast(out.size()), R); + BOOST_REQUIRE_EQUAL(static_cast(out[static_cast(peer)].size()), peer_len); + for (int j = 0; j < peer_len; ++j) { + BOOST_CHECK_EQUAL(out[static_cast(peer)][static_cast(j)], sparse_tag(peer, me, j)); + } + } + + for (const int S : {1, 2}) { + const int P = R * S; + std::vector>> recv(static_cast(S)); + auto errs = run_hybrid(S, [&](HybridComm &hyb, int u) { + Comm c = Comm::make_hybrid(&hyb, u); + const int g = monoprop::mpi::rank(c); + std::vector> send(static_cast(P)); + for (int t = 0; t < S; ++t) { + const int d = (peer * S) + t; + for (int j = 0; j < my_len; ++j) { + send[static_cast(d)].push_back(sparse_tag(g, d, j)); + } + } + std::vector> out; + monoprop::mpi::begin_alltoallv(send, c, false, nullptr, plan).wait_into(out); + recv[static_cast(u)] = out; + }); + for (const auto &e : errs) { + BOOST_CHECK(e == nullptr); + } + for (int t = 0; t < S; ++t) { + const int g = (me * S) + t; + const auto &out = recv[static_cast(t)]; + BOOST_REQUIRE_EQUAL(static_cast(out.size()), P); + for (int su = 0; su < S; ++su) { + const auto &blk = out[static_cast((peer * S) + su)]; + BOOST_REQUIRE_EQUAL(static_cast(blk.size()), peer_len); + for (int j = 0; j < peer_len; ++j) { + BOOST_CHECK_EQUAL(blk[static_cast(j)], sparse_tag((peer * S) + su, g, j)); + } + } + } + } + } +} + +// shift == 0 makes every rank its own and only peer, so the sparse path's self slot is the whole +// exchange -- and skip_self then removes the one leg that would have moved anything for a partition. +// The remaining S-1 in-rank legs must still arrive. +BOOST_AUTO_TEST_CASE(hybrid_comm_sparse_plan_skip_self_at_shift_zero) { + const int R = world_size(); + if (R < 2 || (R & (R - 1)) != 0) { + return; + } + const int me = world_rank(); + const monoprop::mpi::PeerPlan plan{.sparse = true, .shift = 0}; + BOOST_REQUIRE_EQUAL(plan.peer(me, 0), me); + + { + Comm c{MPI_COMM_WORLD}; // the self peer is the ONLY peer, and skip_self drops it + std::vector> send(static_cast(R)); + for (int j = 0; j < 5; ++j) { + send[static_cast(me)].push_back(sparse_tag(me, me, j)); + } + std::vector> out; + monoprop::mpi::begin_alltoallv(send, c, /*skip_self=*/true, nullptr, plan).wait_into(out); + BOOST_REQUIRE_EQUAL(static_cast(out.size()), R); + for (const auto &blk : out) { + BOOST_CHECK(blk.empty()); + } + } + + for (const int S : {2, 3}) { + const int P = R * S; + std::vector>> recv(static_cast(S)); + auto errs = run_hybrid(S, [&](HybridComm &hyb, int u) { + Comm c = Comm::make_hybrid(&hyb, u); + const int g = monoprop::mpi::rank(c); + std::vector> send(static_cast(P)); + for (int t = 0; t < S; ++t) { + const int d = (me * S) + t; + for (int j = 0; j <= t; ++j) { + send[static_cast(d)].push_back(sparse_tag(g, d, j)); + } + } + std::vector> out; + monoprop::mpi::begin_alltoallv(send, c, /*skip_self=*/true, nullptr, plan).wait_into(out); + recv[static_cast(u)] = out; + }); + for (const auto &e : errs) { + BOOST_CHECK(e == nullptr); + } + for (int t = 0; t < S; ++t) { + const int g = (me * S) + t; + const auto &out = recv[static_cast(t)]; + BOOST_REQUIRE_EQUAL(static_cast(out.size()), P); + for (int src = 0; src < P; ++src) { + const bool in_rank = (src / S) == me; + const int want = (in_rank && src != g) ? t + 1 : 0; // own slot dropped by skip_self + BOOST_REQUIRE_EQUAL(static_cast(out[static_cast(src)].size()), want); + for (int j = 0; j < want; ++j) { + BOOST_CHECK_EQUAL(out[static_cast(src)][static_cast(j)], sparse_tag(src, g, j)); + } + } + } + } +} + +// Engine.h's run_exchange posts both its rounds on ONE communicator under kFlatPayloadTag, and the +// non-overtaking argument in Pairwise.h is what says a round-2 receive cannot match a round-1 send. +// Two back-to-back sparse rounds, the second on the transpose of the first's counts, assert it. +BOOST_AUTO_TEST_CASE(hybrid_comm_sparse_plan_back_to_back_rounds) { + const int R = world_size(); + if (R < 2 || (R & (R - 1)) != 0) { + return; + } + const int me = world_rank(); + Comm c{MPI_COMM_WORLD}; + for (int shift = 0; shift < R; ++shift) { + const monoprop::mpi::PeerPlan plan{.sparse = true, .shift = shift}; + const int peer = plan.peer(me, 0); + const int len = 3 + (me % 2); // asymmetric, so a swapped round is a length mismatch + + std::vector> q(static_cast(R)); + for (int j = 0; j < len; ++j) { + q[static_cast(peer)].push_back(sparse_tag(me, peer, j)); + } + std::vector> q_out; + monoprop::mpi::begin_alltoallv(q, c, false, nullptr, plan).wait_into(q_out); + + // Round 2 immediately, same comm and tag, sized from round 1's transpose -- the response shape. + std::vector known(static_cast(R), 0); + known[static_cast(peer)] = len; + std::vector> r(static_cast(R)); + const int back = static_cast(q_out[static_cast(peer)].size()); + for (int j = 0; j < back; ++j) { + r[static_cast(peer)].push_back(q_out[static_cast(peer)][static_cast(j)] + 7); + } + std::vector> r_out; + monoprop::mpi::begin_alltoallv(r, c, false, &known, plan).wait_into(r_out); + + BOOST_REQUIRE_EQUAL(back, 3 + (peer % 2)); + BOOST_REQUIRE_EQUAL(static_cast(r_out[static_cast(peer)].size()), len); + for (int j = 0; j < len; ++j) { + BOOST_CHECK_EQUAL(r_out[static_cast(peer)][static_cast(j)], sparse_tag(me, peer, j) + 7); + } + for (int src = 0; src < R; ++src) { + if (src != peer) { + BOOST_CHECK(r_out[static_cast(src)].empty()); + } + } + } +} + +namespace { + +// One fused-resolve round under `plan`: flattens send_blocks, runs the verb, and re-splits the payload +// by global source using the resolved counts, so a case compares delivered BLOCKS rather than offsets. +auto resolve_round(HybridComm &hyb, + int u, + int P, + const std::vector> &send_blocks, + monoprop::mpi::PeerPlan plan) -> std::vector> { + std::vector send; + std::vector sc(static_cast(P)), sd(static_cast(P)); + for (int d = 0; d < P; ++d) { + const auto &blk = send_blocks[static_cast(d)]; + sc[static_cast(d)] = static_cast(blk.size()); + sd[static_cast(d)] = static_cast(send.size()); + send.insert(send.end(), blk.begin(), blk.end()); + } + std::vector recv; + std::vector rc(static_cast(P)), rd(static_cast(P)); + hyb.alltoallv_resolve(u, + {.send = send.data(), + .send_counts = sc.data(), + .send_displs = sd.data(), + .recv = recv, + .recv_counts = rc.data(), + .recv_displs = rd.data()}, + monoprop::mpi::datatype::get(), + plan); + std::vector> out(static_cast(P)); + for (int src = 0; src < P; ++src) { + const auto off = static_cast(rd[static_cast(src)]); + const auto n = static_cast(rc[static_cast(src)]); + out[static_cast(src)].assign(recv.begin() + static_cast(off), + recv.begin() + static_cast(off + n)); + } + return out; +} + +// Every leg INTO an odd rank is empty, so a multi-peer plan always has a peer whose payload legs are +// skipped entirely while its S*S-int count block still travels. Both ends read the same function of +// (src, dst), which is what keeps the skip symmetric. +auto muted_count(int src, int dst, int S) -> int { + return (dst / S) % 2 == 1 ? 0 : sparse_count(src, dst); +} + +} // namespace + +// The fused resolve POSTS its count round in B1→B2 and drains it in B3→B4, with every partition's +// pack_send_ running underneath. Only the sparse arm splits -- plan.dense() is a blocking MPI_Alltoall +// -- so drive the SAME peer-masked layout through both arms and require the delivered blocks to agree +// element for element. A dropped or torn count block changes a length here, not just a value. +BOOST_AUTO_TEST_CASE(hybrid_comm_resolve_split_count_round_matches_the_dense_arm) { + const int R = world_size(); + if (R < 2 || (R & (R - 1)) != 0) { + return; // the XOR pairing needs a power-of-two rank count + } + const int me = world_rank(); + int cases = 0; + { + constexpr int f = 1; // a sparse plan is fanout 1 by construction; f > 1 is no longer expressible + for (int shift = 0; shift < R; ++shift) { + const monoprop::mpi::PeerPlan plan{.sparse = true, .shift = shift}; + BOOST_REQUIRE_EQUAL(plan.count(R), f); + ++cases; + for (const int S : {1, 2, 3}) { + const int P = R * S; + // Peer-masked, so the dense plan carries the identical bytes: its non-peer blocks are + // empty rather than absent, and a dense count of zero and a sparse absence must agree. + const auto fill = [&](int g) { + std::vector> send(static_cast(P)); + for (int k = 0; k < f; ++k) { + const int b = plan.peer(me, k); + for (int t = 0; t < S; ++t) { + const int d = (b * S) + t; + for (int j = 0; j < muted_count(g, d, S); ++j) { + send[static_cast(d)].push_back(sparse_tag(g, d, j)); + } + } + } + return send; + }; + std::vector>> dense(static_cast(S)); + std::vector>> sparse(static_cast(S)); + auto errs = run_hybrid(S, [&](HybridComm &hyb, int u) { + const int g = (me * S) + u; + dense[static_cast(u)] = resolve_round(hyb, u, P, fill(g), {}); + sparse[static_cast(u)] = resolve_round(hyb, u, P, fill(g), plan); + }); + for (const auto &e : errs) { + BOOST_CHECK(e == nullptr); + } + for (int t = 0; t < S; ++t) { + const int g = (me * S) + t; + const auto &got = sparse[static_cast(t)]; + const auto &want_dense = dense[static_cast(t)]; + BOOST_REQUIRE_EQUAL(static_cast(got.size()), P); + BOOST_REQUIRE_EQUAL(static_cast(want_dense.size()), P); + for (int src = 0; src < P; ++src) { + const int want = plan.contains(me, src / S) ? muted_count(src, g, S) : 0; + const auto &blk = got[static_cast(src)]; + BOOST_REQUIRE_EQUAL(static_cast(blk.size()), want); + BOOST_REQUIRE_EQUAL(static_cast(want_dense[static_cast(src)].size()), want); + for (int j = 0; j < want; ++j) { + BOOST_CHECK_EQUAL(blk[static_cast(j)], sparse_tag(src, g, j)); + BOOST_CHECK_EQUAL(blk[static_cast(j)], + want_dense[static_cast(src)][static_cast(j)]); + } + } + } + } + } + } + BOOST_TEST(cases > 0); // at R < 2 the case is a no-op and must not read as coverage +} + +// shift == 0 at full bits makes this rank its own and only peer, so the count round posts NOTHING: +// sparse_pairwise memcpys the S*S-int block in place and the B3→B4 wait must be a no-op rather than a +// wait on stale requests. The S^2 in-rank payload legs still have to arrive. +BOOST_AUTO_TEST_CASE(hybrid_comm_resolve_split_count_round_self_peer_only) { + const int R = world_size(); + if (R < 2 || (R & (R - 1)) != 0) { + return; + } + const int me = world_rank(); + const monoprop::mpi::PeerPlan plan{.sparse = true, .shift = 0}; + BOOST_REQUIRE_EQUAL(plan.peer(me, 0), me); + for (const int S : {1, 2, 3}) { + const int P = R * S; + std::vector>> recv(static_cast(S)); + auto errs = run_hybrid(S, [&](HybridComm &hyb, int u) { + const int g = (me * S) + u; + std::vector> send(static_cast(P)); + for (int t = 0; t < S; ++t) { + const int d = (me * S) + t; + for (int j = 0; j <= t; ++j) { + send[static_cast(d)].push_back(sparse_tag(g, d, j)); + } + } + recv[static_cast(u)] = resolve_round(hyb, u, P, send, plan); + }); + for (const auto &e : errs) { + BOOST_CHECK(e == nullptr); + } + for (int t = 0; t < S; ++t) { + const int g = (me * S) + t; + const auto &out = recv[static_cast(t)]; + BOOST_REQUIRE_EQUAL(static_cast(out.size()), P); + for (int src = 0; src < P; ++src) { + const int want = (src / S) == me ? t + 1 : 0; // block (su -> t) has t+1 entries + const auto &blk = out[static_cast(src)]; + BOOST_REQUIRE_EQUAL(static_cast(blk.size()), want); + for (int j = 0; j < want; ++j) { + BOOST_CHECK_EQUAL(blk[static_cast(j)], sparse_tag(src, g, j)); + } + } + } + } +} + +// A peer with nothing to send is where the split can hang: pack_send_ copies zero bytes in B2→B3 while +// the count round is still live, and the payload round posts no leg for it at all. Back-to-back rounds, +// the second silent, so the second also runs over the first's staging high-water bytes. +BOOST_AUTO_TEST_CASE(hybrid_comm_resolve_split_count_round_zero_count_peer) { + const int R = world_size(); + if (R < 2 || (R & (R - 1)) != 0) { + return; + } + const int me = world_rank(); + for (int shift = 1; shift < R; ++shift) { // shift 0 is the self peer, covered above + const monoprop::mpi::PeerPlan plan{.sparse = true, .shift = shift}; + const int peer = plan.peer(me, 0); + BOOST_REQUIRE(peer != me); + const int my_len = me < peer ? 0 : 4; // exactly one end of the pair is silent + const int peer_len = peer < me ? 0 : 4; + for (const int S : {1, 2}) { + const int P = R * S; + std::vector>> loud(static_cast(S)); + std::vector>> quiet(static_cast(S)); + auto errs = run_hybrid(S, [&](HybridComm &hyb, int u) { + const int g = (me * S) + u; + std::vector> send(static_cast(P)); + for (int t = 0; t < S; ++t) { + const int d = (peer * S) + t; + for (int j = 0; j < my_len; ++j) { + send[static_cast(d)].push_back(sparse_tag(g, d, j)); + } + } + loud[static_cast(u)] = resolve_round(hyb, u, P, send, plan); + quiet[static_cast(u)] = + resolve_round(hyb, u, P, std::vector>(static_cast(P)), plan); + }); + for (const auto &e : errs) { + BOOST_CHECK(e == nullptr); + } + for (int t = 0; t < S; ++t) { + const int g = (me * S) + t; + const auto &out = loud[static_cast(t)]; + BOOST_REQUIRE_EQUAL(static_cast(out.size()), P); + for (int su = 0; su < S; ++su) { + const auto &blk = out[static_cast((peer * S) + su)]; + BOOST_REQUIRE_EQUAL(static_cast(blk.size()), peer_len); + for (int j = 0; j < peer_len; ++j) { + BOOST_CHECK_EQUAL(blk[static_cast(j)], sparse_tag((peer * S) + su, g, j)); + } + } + // The all-silent round: every resolved count is zero, so a stale staged byte cannot hide. + for (const auto &blk : quiet[static_cast(t)]) { + BOOST_CHECK(blk.empty()); + } + } + } + } +} + +// derived_wire_plan_ is what lets alltoallv narrow its own wire when the caller cannot: only partition 0 +// reaches MPI, and ITS row may be the empty one while a sibling partition holds the rank's only traffic. +// Reading the peer set off partition 0's row alone resolves to the SELF peer, whose legs are all zero -- +// every block is then dropped with no hang to show for it. So this case puts the traffic where partition +// 0 cannot see it, and checks the payload arrives carrying its source's global id. +// +// alltoallv's derive_wire_bits has no library caller today, so this is its only exercise. +BOOST_AUTO_TEST_CASE(hybrid_comm_derived_wire_plan_reads_every_partitions_row) { + const int R = world_size(); + if (R < 2 || (R & (R - 1)) != 0) { + return; // the XOR pairing needs a power-of-two rank count + } + const int bits = std::countr_zero(static_cast(R)); + const int me = world_rank(); + const int peer = me ^ 1; // shift 1, so every rank derives the same pairing + constexpr int kLen = 5; + + for (const int S : {2, 4}) { + const int P = R * S; + const int carrier = S - 1; // never partition 0 -- that is the whole point of the case + std::vector> got(static_cast(S)); + auto errs = run_hybrid(S, [&](HybridComm &hyb, int u) { + const int g = (me * S) + u; + std::vector counts(static_cast(P), 0); + const std::vector displs(static_cast(P), 0); + std::vector send; + if (u == carrier) { + counts[static_cast((peer * S) + carrier)] = kLen; + for (int j = 0; j < kLen; ++j) { + send.push_back((g * 1000) + j); + } + } + // Symmetric layout: my peer's carrier partition sends me exactly what I send it, which is + // what derive_wire_bits requires and what lets the recv rows name the peer set. One block, so + // both displacements are zero. + std::vector recv(static_cast(u == carrier ? kLen : 0), -1); + const monoprop::mpi::AlltoallvArgs args{.send = reinterpret_cast(send.data()), + .send_counts = counts.data(), + .send_displs = displs.data(), + .recv = reinterpret_cast(recv.data()), + .recv_counts = counts.data(), + .recv_displs = displs.data(), + .elem = sizeof(int)}; + hyb.alltoallv(u, args, MPI_INT, monoprop::mpi::PeerPlan{}, /*derive_wire_bits=*/bits); + got[static_cast(u)] = recv; + }); + for (const auto &e : errs) { + BOOST_CHECK(e == nullptr); + } + for (int u = 0; u < S; ++u) { + if (u != carrier) { + BOOST_CHECK(got[static_cast(u)].empty()); + continue; + } + // The values name their sender, so this fails on a plan that named the wrong peer as well as + // on one that named none. + const int src = (peer * S) + carrier; + BOOST_REQUIRE_EQUAL(static_cast(got[static_cast(u)].size()), kLen); + for (int j = 0; j < kLen; ++j) { + BOOST_CHECK_EQUAL(got[static_cast(u)][static_cast(j)], (src * 1000) + j); + } + } + } +} + #endif // monoprop_ENABLE_MPI diff --git a/cpp/tests/mpi_utils_tests.cpp b/cpp/tests/mpi_utils_tests.cpp index 45fdaf9e..3fabc300 100644 --- a/cpp/tests/mpi_utils_tests.cpp +++ b/cpp/tests/mpi_utils_tests.cpp @@ -27,14 +27,17 @@ #include "monoprop/algebra/MajoranaAlgebra.h" #include "monoprop/detail/evolution/CutoffContext.h" #include "monoprop/detail/evolution/layer_build/Scan.h" +#include "monoprop/detail/mpi/Comm.h" #include "monoprop/detail/mpi/MPIUtils.h" +#include "monoprop/detail/mpi/Routing.h" #include "monoprop/detail/operator/MPOperator.h" #include "monoprop/detail/operator/OperatorIndex.h" using namespace monoprop; -// find_rank is splitmix over the dense words modulo the rank count, and nothing else, so the oracle -// is asserted unconditionally rather than as one of several permitted hashes. +// Under the splitmix router find_rank is the dense words modulo the rank count and nothing else, so the +// oracle is asserted unconditionally rather than as one of several permitted hashes. The router is +// constructed explicitly: there is no rank-count overload to reach it by accident. BOOST_AUTO_TEST_CASE(mpi_utils_find_rank_range_and_hash_mod) { constexpr size_t N = 32; std::mt19937_64 rng(0x9E3779B9ULL); @@ -46,19 +49,21 @@ BOOST_AUTO_TEST_CASE(mpi_utils_find_rank_range_and_hash_mod) { } const auto mono = indices_to_bitset(inds); for (size_t n_ranks : {size_t{1}, size_t{2}, size_t{3}, size_t{7}}) { - const size_t r = find_rank(mono, n_ranks); + const auto router = routing::Router::splitmix(n_ranks); + const size_t r = find_rank(mono, router); BOOST_TEST(r == monomial_hash(mono) % n_ranks); BOOST_TEST(r < n_ranks); - BOOST_TEST(r == find_rank(mono, n_ranks)); // deterministic + BOOST_TEST(r == find_rank(mono, router)); // deterministic } } } -// n_ranks == 0 is degenerate: owner is rank 0, not a modulo by zero. +// A zero-rank world is degenerate: Router clamps it to one slot, so the owner is rank 0 rather than a +// modulo by zero. BOOST_AUTO_TEST_CASE(mpi_utils_find_rank_zero_ranks) { constexpr size_t N = 32; const auto mono = indices_to_bitset(VecZ{0, 3, 5}); - BOOST_TEST(find_rank(mono, 0) == 0U); + BOOST_TEST(find_rank(mono, routing::Router::splitmix(0)) == 0U); } BOOST_AUTO_TEST_CASE(mpi_utils_monomial_words_roundtrip) { @@ -111,35 +116,42 @@ auto build_op(const std::vector> &terms) -> detail::MPOperator<32> return op; } -auto check_bucket_ownership(const std::vector &buckets, size_t ranks, size_t &checked) -> void { +auto check_bucket_ownership(const mpi::WindowVec &buckets, const routing::Router &router, size_t &checked) + -> void { // Every offset comes from the codec's walk: the record is VARIABLE WIDTH, so a hardcoded stride - // would compare a monomial decoded at the wrong offset against the wrong rank. + // would compare a monomial decoded at the wrong offset against the wrong rank. The bucket index is + // re-based, so the rank compared against is window.slot(k) -- a wrong window base shows up here. using QC = detail::QueryCodec<32>; const detail::QueryLayout layout{/*fused=*/false}; - for (size_t r = 0; r < buckets.size(); ++r) { + const mpi::SlotWindow w = buckets.window(); + for (size_t k = 0; k < w.count; ++k) { + const mpi::WindowIndex wi{k}; + const VecZ &bucket = buckets[wi]; size_t off = 0; - while (off < buckets[r].size()) { + while (off < bucket.size()) { Monomial<32> mono; int phase = 0; - QC::read_mono(buckets[r], off, mono, phase); - BOOST_REQUIRE_EQUAL(find_rank<32>(mono, ranks), r); - off = QC::next_off(buckets[r], layout, off); + QC::read_mono(bucket, off, mono, phase); + BOOST_REQUIRE_EQUAL(find_rank<32>(mono, router), w.slot(wi)); + off = QC::next_off(bucket, layout, off); ++checked; } - BOOST_REQUIRE_EQUAL(off, buckets[r].size()); + BOOST_REQUIRE_EQUAL(off, bucket.size()); } } // The self-owned bucket is staged as positions, not encoded, so it is invisible to the walk above -- // without this the r == my_rank arm of the routing decision goes unchecked. -auto check_self_ownership(const detail::SelfQueryStage<32> &stage, size_t ranks, size_t my_rank, size_t &checked) - -> void { +auto check_self_ownership(const detail::SelfQueryStage<32> &stage, + const routing::Router &router, + size_t my_rank, + size_t &checked) -> void { for (size_t q = 0; q < stage.size(); ++q) { Monomial<32> mono; for (size_t j = 0; j < stage.k_of[q]; ++j) { mono.set(static_cast(stage.pos_flat[stage.pos_off[q] + j])); } - BOOST_REQUIRE_EQUAL(find_rank<32>(mono, ranks), my_rank); + BOOST_REQUIRE_EQUAL(find_rank<32>(mono, router), my_rank); ++checked; } } @@ -157,47 +169,105 @@ BOOST_AUTO_TEST_CASE(mpi_utils_scan_routing_agrees_with_find_rank) { for (size_t i = 0; i < 2000; ++i) { terms.push_back(draw_well_formed(rng, kLogical, 1 + (rng() % 6))); } - auto op = build_op(terms); const Monomial gen = draw_well_formed(rng, kLogical, 4); - VecD coeffs(op.store->size(), 1.0); const CutoffFn fn = detail::LengthCutoff{10, kLogical}; const detail::CutoffEvaluator eval(fn); - const auto cut = detail::build_majorana_evolution_cutoff_state(std::nullopt, - std::cref(coeffs), - std::nullopt, - std::optional{0.3}); - size_t checked = 0; - size_t self_checked = 0; + // Each rank is handed exactly the terms it owns, as MonomialPropagator seeds it: the emit path routes + // by rank(M) ^ rank_shift(G), which is only the owner of M^G when M really is local. + // + // The partner count is a property of the operator and the gate, not of where the partners live, so + // the sum over the ranks is the same for every router; routing only moves a partner between ranks and + // between the encoded (cross-rank) and staged (self-owned) sides. Pinning that invariance is stronger + // than a floor: a routing bug that drops partners moves the total, and one that misroutes them moves + // the split. + size_t routers = 0; + std::optional first_total; for (const size_t ranks : {2U, 4U, 8U}) { - const auto res = detail::fused_find_and_collect>(op, - gen, - eval, - cut, - coeffs, - std::nullopt, - ranks, - 0, - false, - nullptr, - 1.0); - BOOST_REQUIRE_EQUAL(res.leader_queries.size(), ranks); - // The scan routes a self-owned partner to the stage, so bucket 0 must be empty here. - BOOST_REQUIRE(res.leader_queries[0].empty()); - BOOST_REQUIRE(res.follower_queries[0].empty()); - check_bucket_ownership(res.leader_queries, ranks, checked); - check_bucket_ownership(res.follower_queries, ranks, checked); - check_self_ownership(res.leader_self, ranks, /*my_rank=*/0, self_checked); - check_self_ownership(res.follower_self, ranks, /*my_rank=*/0, self_checked); + // BOTH routers, because the agreement is a property of the pair and not of either hash: the scan + // calls Router::dest_from_shift and find_rank calls Router::dest, so a divergence introduced by + // one of them shows up here whichever routing the geometry resolves to. + for (const bool linear : {false, true}) { + const auto router = routing::Router::for_modes(ranks, /*partitions=*/1, linear); + const size_t shift = router.rank_shift(gen); + const mpi::PeerPlan plan{.sparse = router.is_linear(), .shift = static_cast(shift)}; + // Per router, not summed over them: the floors are what stops the loop passing on an empty + // scan, and a sum lets one router carry the other. + size_t checked = 0; + size_t self_checked = 0; + for (size_t my_rank = 0; my_rank < ranks; ++my_rank) { + std::vector> owned; + for (const auto &t : terms) { + if (find_rank(t, router) == my_rank) { + owned.push_back(t); + } + } + BOOST_REQUIRE(!owned.empty()); + auto op = build_op(owned); + const mpi::SlotWindow window = plan.window(my_rank, ranks, /*parts=*/1); + VecD coeffs(op.store->size(), 1.0); + const auto cut = detail::build_majorana_evolution_cutoff_state(std::nullopt, + std::cref(coeffs), + std::nullopt, + std::optional{0.3}); + const auto res = detail::fused_find_and_collect>(op, + gen, + eval, + cut, + coeffs, + std::nullopt, + window, + my_rank, + router, + shift, + false, + nullptr, + 1.0); + BOOST_REQUIRE_EQUAL(res.leader_queries.size(), window.count); + if (window.contains(my_rank)) { + // The scan routes a self-owned partner to the stage, so my own bucket must be empty. + BOOST_REQUIRE(res.leader_queries.at_slot(my_rank).empty()); + BOOST_REQUIRE(res.follower_queries.at_slot(my_rank).empty()); + } + else { + // A non-zero shift puts self outside the window entirely, so nothing may be staged. + BOOST_REQUIRE_EQUAL(res.leader_self.size(), 0U); + BOOST_REQUIRE_EQUAL(res.follower_self.size(), 0U); + } + check_bucket_ownership(res.leader_queries, router, checked); + check_bucket_ownership(res.follower_queries, router, checked); + check_self_ownership(res.leader_self, router, my_rank, self_checked); + check_self_ownership(res.follower_self, router, my_rank, self_checked); + } + BOOST_TEST_MESSAGE("ranks=" << ranks << " linear=" << router.is_linear() << " shift=" << shift + << " encoded=" << checked << " staged=" << self_checked); + const size_t total = checked + self_checked; + if (first_total.has_value()) { + BOOST_TEST(total == *first_total); // routing moves partners, it does not create or lose them + } + else { + first_total = total; + } + // Measured, all six routers: total 387 every time, with the splitmix split running from + // 212/175 at R=2 to 343/44 at R=8 as more partners fall cross-rank, and the linear arm all + // encoded (shift 1, 3, 7). The floors sit below the observed minimum of each arm and only + // catch a scan that emitted nothing. + BOOST_TEST(total > 300U); + if (!router.is_linear()) { + // Splitmix re-hashes the partner, so both sides are populated at every rank count. + BOOST_TEST(checked > 150U); + BOOST_TEST(self_checked > 40U); + } + else if (shift != 0) { + // Fanout 1: every partner leaves for my_rank ^ shift, so nothing stays self-owned. + BOOST_TEST(self_checked == 0U); + } + else { + BOOST_TEST(checked == 0U); // a shift of zero keeps every partner on its own rank + } + ++routers; + } } - // Without this the loop above passes trivially if the scan emitted nothing. The floor is on the SUM - // because that is what is invariant across the split: the encoded counter alone fell to 797 of 1161 - // when the self-owned partners moved into the stage, with nothing going unchecked. Each arm still - // carries its own floor -- a routing bug sending everything one way leaves the sum intact -- and the - // message prints the measured 797/364 so those can be re-grounded rather than guessed. - BOOST_TEST_MESSAGE("encoded=" << checked << " staged=" << self_checked); - BOOST_TEST(checked + self_checked > 1000U); - BOOST_TEST(checked > 500U); - BOOST_TEST(self_checked > 200U); + BOOST_TEST(routers == 6U); // the floors above are per router, so the router count is part of them } diff --git a/cpp/tests/routing_tests.cpp b/cpp/tests/routing_tests.cpp new file mode 100644 index 00000000..f7a413d0 --- /dev/null +++ b/cpp/tests/routing_tests.cpp @@ -0,0 +1,447 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +// routing::Router -- the term -> flat-slot map. Two states, two properties that carry the whole design: +// * splitmix is bit-for-bit today's `monomial_hash % P`, so the refactor cannot move a single term; +// * under linear routing the destination RANK of M^G is rank(M) ^ shift(G), which is what turns the +// dense all-to-all into a pairwise exchange. A break here is silent -- wrong owner, not a crash -- +// so the shift identity and the Scan/find_rank agreement are both asserted explicitly. +// +// Flat cases with a shared prefix, no suite nesting (suites break Boost's ctest discovery here). + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "monoprop/algebra/MajoranaAlgebra.h" +#include "monoprop/detail/mpi/Comm.h" +#include "monoprop/detail/mpi/MPIUtils.h" +#include "monoprop/detail/mpi/Routing.h" + +using namespace monoprop; +using monoprop::routing::Router; + +namespace { + +constexpr size_t kN = 64; // 2N = 128 bits -> 2 words, so the multi-word hash path is exercised + +// Exactly `weight` distinct bits, so a case can pin the popcount the old bit-walk paid for. +auto monomials_of_weight(size_t count, size_t weight, uint64_t seed) -> std::vector> { + std::mt19937_64 rng(seed); + std::uniform_int_distribution slot(0, 2 * kN - 1); + std::vector> out; + out.reserve(count); + for (size_t i = 0; i < count; ++i) { + Monomial m; + while (m.count() < weight) { + m.set(slot(rng)); + } + out.push_back(m); + } + return out; +} + +// The map Router::dest computed before the basis was transposed: load one 64-bit vector per SET bit and +// XOR. Rebuilt from mix64 and the seed rather than read from linear_basis(), so the plane build and this +// share nothing but the specification. +template +auto linear_hash_reference(const monoprop::Bitset &bits) -> uint64_t { + const uint64_t seed = routing::seed_from_env(); + uint64_t h = 0; + for (size_t i = bits.find_first(); i < NumBits; i = bits.find_next(i)) { + h ^= routing::mix64(routing::mix64(seed) + (static_cast(i) * 0x9E37'79B9'7F4A'7C15ULL)); + } + return h; +} + +auto random_monomials(size_t count, size_t weight, uint64_t seed) -> std::vector> { + std::mt19937_64 rng(seed); + std::uniform_int_distribution slot(0, 2 * kN - 1); + std::vector> out; + out.reserve(count); + for (size_t i = 0; i < count; ++i) { + VecZ inds; + for (size_t k = 0; k < weight; ++k) { + inds.push_back(slot(rng)); + } + out.push_back(indices_to_bitset(inds)); + } + return out; +} + +} // namespace + +// The splitmix router must not move a single term relative to `monomial_hash % P`: this is the +// regression gate that licenses everything else. +BOOST_AUTO_TEST_CASE(routing_splitmix_is_bit_identical_to_hash_mod_p) { + const auto monos = random_monomials(500, 5, 0xC0FFEEULL); + for (const size_t flat : {size_t{1}, size_t{2}, size_t{7}, size_t{112}, size_t{1792}}) { + const auto router = Router::splitmix(flat); + BOOST_TEST(!router.is_linear()); + for (const auto &m : monos) { + const size_t expected = monomial_hash(m) % flat; + BOOST_TEST(router.dest(m) == expected); + BOOST_TEST(find_rank(m, router) == expected); + } + } +} + +// A two-level splitmix router is still today's routing, even though it knows about partitions: +// ((q/S) % R)*S + q%S == q % (R*S). 12 ranks is here too: splitmix carries no power-of-two condition. +BOOST_AUTO_TEST_CASE(routing_splitmix_two_level_collapses_to_flat_modulo) { + const auto monos = random_monomials(300, 6, 0xBEEF01ULL); + for (const auto [r, s] : {std::pair{8, 14}, {4, 28}, {128, 14}, {64, 28}, {12, 5}}) { + const auto router = Router::for_modes(r, s, false); + for (const auto &m : monos) { + BOOST_TEST(router.dest(m) == monomial_hash(m) % (r * s)); + } + } +} + +BOOST_AUTO_TEST_CASE(routing_linear_hash_is_gf2_linear) { + const auto a = random_monomials(200, 5, 0x1234ULL); + const auto b = random_monomials(200, 3, 0x5678ULL); + for (size_t i = 0; i < a.size(); ++i) { + const uint64_t ha = routing::linear_hash<2 * kN>(a[i]); + const uint64_t hb = routing::linear_hash<2 * kN>(b[i]); + BOOST_TEST(routing::linear_hash<2 * kN>(a[i] ^ b[i]) == (ha ^ hb)); + } + // The empty monomial is the identity of the group, so its hash must be the identity of the codomain. + BOOST_TEST(routing::linear_hash<2 * kN>(Monomial{}) == 0ULL); +} + +// The load-bearing identity: under linear routing the destination RANK of M^G is rank(M) ^ shift(G), +// so one rank's queries for one generator all land on one peer. +BOOST_AUTO_TEST_CASE(routing_shift_identity_holds_under_linear_routing) { + constexpr size_t kRanks = 16; + constexpr size_t kParts = 14; + const auto router = Router::for_modes(kRanks, kParts, true); + BOOST_REQUIRE(router.is_linear()); + BOOST_REQUIRE(router.linear_bits() == 4U); // every rank bit, log2(16) + + const auto terms = random_monomials(400, 6, 0xAAAA01ULL); + const auto gens = random_monomials(40, 4, 0xBBBB02ULL); + for (const auto &g : gens) { + const size_t shift = router.rank_shift(g); + for (const auto &m : terms) { + const size_t rank_m = router.dest(m) / kParts; + const size_t rank_mg = router.dest(m ^ g) / kParts; + BOOST_TEST(rank_mg == (rank_m ^ shift)); + } + } +} + +// What the emit path actually calls. rank(M^G) == rank(M) ^ shift(G) is an identity, so the check is +// exact equality against dest() over every geometry and both modes -- including splitmix, where +// dest_from_shift is a pass-through and must not move a term either. S == 1 (no partition index at all), +// S a power of two (the masked modulo) and S = 3, 14 (the division) are all here. +BOOST_AUTO_TEST_CASE(routing_dest_from_shift_agrees_with_dest) { + const auto terms = random_monomials(400, 6, 0x5E1F7A11ULL); + const auto gens = random_monomials(25, 4, 0x9110F7E5ULL); + const std::vector> + geometries{{1, 1}, {8, 1}, {16, 1}, {2, 2}, {32, 16}, {64, 8}, {4, 3}, {8, 14}, {1, 14}}; + size_t checked = 0; + for (const auto &[r, s] : geometries) { + for (const bool linear : {false, true}) { + const auto router = Router::for_modes(r, s, linear); + for (const auto &g : gens) { + const size_t shift = router.rank_shift(g); + for (const auto &m : terms) { + // dest(m) stands in for the flat slot the owning rank would pass in. + const auto partner = m ^ g; + BOOST_REQUIRE_EQUAL(router.dest_from_shift(partner, router.dest(m), shift), + router.dest(partner)); + ++checked; + } + } + } + } + BOOST_TEST_MESSAGE("dest_from_shift checks: " << checked); + BOOST_TEST(checked >= 100000U); +} + +// The S == 1 skip: the partition index is 0 for every term, so under linear routing the destination is +// the rank index alone and monomial_hash is not on the path at all. +BOOST_AUTO_TEST_CASE(routing_single_partition_destination_is_the_rank_index) { + constexpr size_t kRanks = 64; + const auto router = Router::for_modes(kRanks, 1, true); + for (const auto &m : random_monomials(500, 5, 0x0FAE7101ULL)) { + BOOST_REQUIRE_EQUAL(router.dest(m), routing::linear_hash<2 * kN>(m) & (kRanks - 1)); + } +} + +// The consequence that the transport will rely on: every term a rank owns sends its query for one +// generator to exactly ONE peer rank -- and the partition index within that peer still spreads. +BOOST_AUTO_TEST_CASE(routing_fanout_is_one_under_linear_routing) { + constexpr size_t kRanks = 8; + constexpr size_t kParts = 14; + const auto router = Router::for_modes(kRanks, kParts, true); + const auto terms = random_monomials(4000, 6, 0xCCCC03ULL); + const auto gens = random_monomials(12, 4, 0xDDDD04ULL); + + std::vector>> owned(kRanks); + for (const auto &m : terms) { + owned[router.dest(m) / kParts].push_back(m); + } + for (size_t src = 0; src < kRanks; ++src) { + BOOST_REQUIRE(!owned[src].empty()); + for (const auto &g : gens) { + std::set dest_ranks; + std::set dest_parts; + for (const auto &m : owned[src]) { + const size_t flat = router.dest(m ^ g); + dest_ranks.insert(flat / kParts); + dest_parts.insert(flat % kParts); + } + BOOST_TEST(dest_ranks.size() == 1U); + BOOST_TEST(dest_parts.size() > 1U); // partitions keep full avalanche + } + } +} + +// mpi::PeerPlan::window -- the slots a plan can reach, which is what lets the per-generator structures +// be S long instead of P. Dense is the count == P value of the same two expressions, so it is checked +// against the same formula rather than against a second one. +BOOST_AUTO_TEST_CASE(routing_slot_window_is_the_peer_ranks_partition_run) { + for (const auto [r, s] : {std::pair{8, 16}, {128, 16}, {16, 1}, {1, 1}, {4, 3}}) { + const size_t p = r * s; + for (size_t me = 0; me < p; ++me) { + const mpi::SlotWindow dense = mpi::PeerPlan{}.window(me, r, s); + BOOST_REQUIRE_EQUAL(dense.base, 0U); + BOOST_REQUIRE_EQUAL(dense.count, p); + BOOST_REQUIRE(dense.contains(me)); + + for (size_t shift = 0; shift < r; ++shift) { + const mpi::PeerPlan plan{.sparse = true, .shift = static_cast(shift)}; + const mpi::SlotWindow w = plan.window(me, r, s); + BOOST_REQUIRE_EQUAL(w.count, s); + BOOST_REQUIRE_EQUAL(w.base, ((me / s) ^ shift) * s); + // Every slot in the run belongs to the one peer rank the plan names, and no other. + for (size_t k = 0; k < w.count; ++k) { + const size_t slot = w.slot(mpi::WindowIndex{k}); + BOOST_REQUIRE(w.contains(slot)); + BOOST_REQUIRE_EQUAL(w.index(slot).value, k); + BOOST_REQUIRE_EQUAL(slot / s, ((me / s) ^ shift)); + BOOST_REQUIRE(plan.contains(static_cast(me / s), static_cast(slot / s))); + } + BOOST_REQUIRE(!w.contains(w.base + w.count)); + // Self is reachable only at shift 0; the engine's self-resolve leg turns on that. + BOOST_REQUIRE_EQUAL(w.contains(me), shift == 0); + } + } + } +} + +// The window must be symmetric, or the two ends of one exchange size different arrays: XOR is an +// involution, so the peer's own window points back at this rank's run. +BOOST_AUTO_TEST_CASE(routing_slot_window_pairing_is_symmetric) { + constexpr size_t kRanks = 32; + constexpr size_t kParts = 8; + for (size_t me = 0; me < kRanks * kParts; me += 3) { + for (size_t shift = 0; shift < kRanks; ++shift) { + const mpi::PeerPlan plan{.sparse = true, .shift = static_cast(shift)}; + const mpi::SlotWindow mine = plan.window(me, kRanks, kParts); + const mpi::SlotWindow theirs = plan.window(mine.base, kRanks, kParts); + BOOST_REQUIRE_EQUAL(theirs.base, (me / kParts) * kParts); + BOOST_REQUIRE_EQUAL(theirs.count, kParts); + } + } +} + +// The property the re-basing rests on: every destination the emit path can produce for one generator +// lies inside that generator's window, so `slot - base` is always a legal index. +BOOST_AUTO_TEST_CASE(routing_every_dest_lands_inside_the_generators_window) { + const auto terms = random_monomials(600, 6, 0x5107500DULL); + const auto gens = random_monomials(20, 4, 0x1CE0FF1CEULL); + const std::vector> geometries{{8, 16}, {16, 1}, {32, 4}, {1, 14}, {4, 3}}; + size_t checked = 0; + for (const auto &[r, s] : geometries) { + for (const bool linear : {false, true}) { + const auto router = Router::for_modes(r, s, linear); + for (const auto &g : gens) { + const size_t shift = router.rank_shift(g); + const mpi::PeerPlan plan{.sparse = router.is_linear(), .shift = static_cast(shift)}; + for (const auto &m : terms) { + const size_t me = router.dest(m); + const mpi::SlotWindow w = plan.window(me, r, s); + BOOST_REQUIRE(w.contains(router.dest_from_shift(m ^ g, me, shift))); + ++checked; + } + } + } + } + BOOST_TEST_MESSAGE("window containment checks: " << checked); + BOOST_TEST(checked >= 100000U); +} + +// WindowVec re-bases in exactly one place, so the flat slot the writer used is the flat slot the reader +// gets back -- including when the run does not start at 0. +BOOST_AUTO_TEST_CASE(routing_window_vec_round_trips_flat_slots) { + const mpi::SlotWindow w{.base = 48, .count = 16}; + mpi::WindowVec v(w); + BOOST_REQUIRE_EQUAL(v.size(), w.count); + for (size_t slot = w.base; slot < w.stop(); ++slot) { + v.at_slot(slot) = slot * 7; + } + for (size_t k = 0; k < w.count; ++k) { + BOOST_REQUIRE_EQUAL(v[mpi::WindowIndex{k}], (w.base + k) * 7); + BOOST_REQUIRE_EQUAL(v.at_slot(w.slot(mpi::WindowIndex{k})), (w.base + k) * 7); + } + v.reset(mpi::SlotWindow{.base = 0, .count = 4}); + BOOST_REQUIRE_EQUAL(v.size(), 4U); + BOOST_REQUIRE_EQUAL(v.at_slot(3), 0U); +} + +// Without a power-of-two rank count there is no XOR structure to exploit, and there is no partial dial +// to fall back to, so the geometry is rejected at construction rather than routed on a subspace. +BOOST_AUTO_TEST_CASE(routing_non_power_of_two_ranks_throw_under_linear_routing) { + for (const size_t r : {size_t{3}, size_t{7}, size_t{12}, size_t{112}}) { + BOOST_CHECK_THROW(static_cast(Router::for_modes(r, 14, true)), routing::UnroutableGeometry); + BOOST_CHECK_NO_THROW(static_cast(Router::for_modes(r, 14, false))); // splitmix has no condition + } + BOOST_CHECK_NO_THROW(static_cast(Router::for_modes(64, 14, true))); +} + +// R = 1 is a power of two, so it must NOT throw: it has no rank bit to take, which makes it the dense +// router and keeps the collective transport that serves every single-rank run. +BOOST_AUTO_TEST_CASE(routing_single_rank_is_dense_and_not_an_error) { + const auto monos = random_monomials(200, 5, 0x51A61EULL); + for (const size_t s : {size_t{1}, size_t{14}, size_t{112}}) { + const auto router = Router::for_modes(1, s, true); + BOOST_TEST(!router.is_linear()); + BOOST_TEST(router.linear_bits() == 0U); + for (const auto &m : monos) { + BOOST_TEST(router.dest(m) == monomial_hash(m) % s); + BOOST_TEST(router.rank_shift(m) == 0U); // no bit to shift, so every generator is on-rank + } + } +} + +BOOST_AUTO_TEST_CASE(routing_dest_is_in_range_and_deterministic) { + const auto monos = random_monomials(1000, 7, 0x9999ULL); + for (const auto [r, s] : {std::pair{1, 1}, {1, 112}, {8, 14}, {128, 14}, {64, 28}}) { + for (const bool linear : {false, true}) { + const auto router = Router::for_modes(r, s, linear); + for (const auto &m : monos) { + const size_t slot = router.dest(m); + BOOST_TEST(slot < r * s); + BOOST_TEST(slot == router.dest(m)); // stateless + } + } + } +} + +// The transposed basis must be the SAME map, not merely a faster one: a divergence is a silently wrong +// owner. Pin dest() and rank_shift() against the old bit-walk over both routers, every geometry, and +// popcounts from empty to full support. +BOOST_AUTO_TEST_CASE(routing_transposed_basis_is_bit_identical_to_the_bit_walk) { + std::vector> monos; + for (const size_t w : {size_t{0}, + size_t{1}, + size_t{2}, + size_t{3}, + size_t{5}, + size_t{8}, + size_t{13}, + size_t{21}, + size_t{34}, + size_t{2 * kN}}) { + const auto batch = monomials_of_weight(500, w, 0xD15EA5E0ULL + w); + monos.insert(monos.end(), batch.begin(), batch.end()); + } + BOOST_REQUIRE_EQUAL(monos.size(), 5000U); + + // (R, S), each run under both routers: log2(R) runs 0, 1, 3, 4, 5, 6, 7, 8, 10, 11, 12, so the plane + // count varies from none to twelve and R = 1 pins the geometry that is dense under either mode. + const std::vector> geometries{{128, 14}, + {16, 1}, + {64, 28}, + {8, 14}, + {1024, 1}, + {2, 112}, + {4096, 16}, + {32, 3}, + {1, 112}, + {256, 2}, + {2048, 1}}; + + size_t checked = 0; + for (const auto &[r, s] : geometries) { + for (const bool linear : {false, true}) { + const auto router = Router::for_modes(r, s, linear); + const size_t d = router.linear_bits(); + BOOST_REQUIRE_EQUAL(d, linear ? static_cast(std::countr_zero(r)) : 0U); + const uint64_t lin_mask = d == 0 ? 0ULL : (uint64_t{1} << d) - 1; + for (const auto &m : monos) { + const uint64_t q = monomial_hash(m); + // d == 0 is the splitmix arm, and R = 1 under linear routing lands there too. + const size_t expected = + d == 0 ? static_cast(q % (r * s)) + : static_cast(((linear_hash_reference<2 * kN>(m) & lin_mask) * s) + (q % s)); + BOOST_REQUIRE_EQUAL(router.dest(m), expected); + BOOST_REQUIRE_EQUAL(router.rank_shift(m), + static_cast(linear_hash_reference<2 * kN>(m) & lin_mask)); + ++checked; + } + } + } + BOOST_TEST_MESSAGE("bit-identity checks: " << checked); + BOOST_TEST(checked >= 100000U); +} + +// gf2_rank is the coverage diagnostic: shifts that span fewer than log2(R) dimensions leave ranks empty. +BOOST_AUTO_TEST_CASE(routing_gf2_rank_detects_a_degenerate_shift_set) { + BOOST_TEST(routing::gf2_rank({}) == 0U); + BOOST_TEST(routing::gf2_rank({0ULL, 0ULL}) == 0U); + BOOST_TEST(routing::gf2_rank({0b001ULL, 0b010ULL, 0b011ULL}) == 2U); // third is the XOR of the first two + BOOST_TEST(routing::gf2_rank({0b001ULL, 0b010ULL, 0b100ULL}) == 3U); + + // The real generator shifts must span at least log2(R) dimensions or the reachable ranks are a + // strict subspace of the rank space. + constexpr size_t kRanks = 128; + const auto router = Router::for_modes(kRanks, 14, true); + std::vector shifts; + for (const auto &g : random_monomials(200, 4, 0x7777ULL)) { + shifts.push_back(static_cast(router.rank_shift(g))); + } + // Seed-independent, so this is NOT skipped when monoprop_ROUTE_SEED is overridden: 200 vectors fail + // to span F_2^7 with probability ~2^-194, whatever basis the seed picks. + BOOST_TEST(routing::gf2_rank(shifts) == 7U); // == log2(128): every rank is reachable +} + +// The SHIPPED default, pinned by a test rather than left to a comment: with no environment override a +// power-of-two rank count routes linearly at fanout 1, a single rank routes densely, and a geometry +// with no XOR structure is refused. Skipped when the environment does override it, because then the +// default is not what is under test. +BOOST_AUTO_TEST_CASE(routing_default_is_linear_where_the_geometry_allows_it) { + const char *mode = std::getenv("monoprop_ROUTING"); + if (mode != nullptr && *mode != '\0') { + BOOST_TEST_MESSAGE("routing overridden in the environment; default not under test"); + return; + } + BOOST_TEST(routing::make_router(8, 14).is_linear()); + BOOST_TEST(routing::make_router(128, 14).linear_bits() == 7U); + BOOST_TEST(!routing::make_router(1, 112).is_linear()); // single rank: nothing to route between + + // 6 and 12 are not powers of two: no XOR structure and no partial dial to retreat to, so the + // geometry is rejected instead of silently routing onto a subspace of the ranks. + BOOST_CHECK_THROW(routing::make_router(6, 14), routing::UnroutableGeometry); + BOOST_CHECK_THROW(routing::make_router(12, 28), routing::UnroutableGeometry); +} diff --git a/cpp/tests/sparse_resolve_tests.cpp b/cpp/tests/sparse_resolve_tests.cpp index 37b70873..31161367 100644 --- a/cpp/tests/sparse_resolve_tests.cpp +++ b/cpp/tests/sparse_resolve_tests.cpp @@ -27,6 +27,7 @@ #include "monoprop/core/Monomial.h" #include "monoprop/detail/evolution/layer_build/QueryCodec.h" #include "monoprop/detail/evolution/layer_build/Resolve.h" +#include "monoprop/detail/mpi/Comm.h" #include "monoprop/detail/operator/MPOperator.h" #include "monoprop/detail/operator/OperatorIndex.h" @@ -106,14 +107,16 @@ auto draw_distinct(std::mt19937_64 &rng, size_t n) -> std::vector -auto serialize(const std::vector>> &queries, bool fused) -> std::vector { - std::vector incoming(queries.size()); +auto serialize(const std::vector>> &queries, bool fused, mpi::SlotWindow window) + -> mpi::WindowVec { + mpi::WindowVec incoming(window); for (size_t s = 0; s < queries.size(); ++s) { + VecZ &buf = incoming[mpi::WindowIndex{s}]; for (size_t q = 0; q < queries[s].size(); ++q) { const int phase = ((q % 2) == 0) ? 1 : -1; - detail::QueryCodec::push(incoming[s], queries[s][q], phase); + detail::QueryCodec::push(buf, queries[s][q], phase); if (fused) { - detail::QueryCodec::push_value(incoming[s], 0.5 + static_cast(q)); + detail::QueryCodec::push_value(buf, 0.5 + static_cast(q)); } } } @@ -121,8 +124,14 @@ auto serialize(const std::vector>> &queries, bool } template -auto check_probe_matches_the_queries(std::mt19937_64 &rng, size_t n_seed, size_t n_query, size_t rank_count, bool fused) - -> void { +auto check_probe_matches_the_queries(std::mt19937_64 &rng, + size_t n_seed, + size_t n_query, + size_t rank_count, + bool fused, + size_t window_base = 0) -> void { + // A non-zero base is the case a re-basing bug survives: sender_slot must still name the flat slot. + const mpi::SlotWindow window{.base = window_base, .count = rank_count}; const auto seed_terms = draw_distinct(rng, n_seed); const auto fresh_terms = draw_distinct(rng, n_query); @@ -159,11 +168,13 @@ auto check_probe_matches_the_queries(std::mt19937_64 &rng, size_t n_seed, size_t } } - const auto incoming = serialize(queries, fused); + const auto incoming = serialize(queries, fused, window); const detail::QueryLayout layout{fused}; auto op = make_op(seed_terms); - const auto pr = detail::probe_incoming_queries(incoming, op, rank_count, layout); + const auto pr = detail::probe_incoming_queries(incoming, op, layout); + BOOST_REQUIRE_EQUAL(pr.window.base, window.base); + BOOST_REQUIRE_EQUAL(pr.window.count, window.count); BOOST_REQUIRE_EQUAL(pr.nq_total, expect_mono.size()); BOOST_REQUIRE(pr.nq_total > 0); @@ -186,7 +197,8 @@ auto check_probe_matches_the_queries(std::mt19937_64 &rng, size_t n_seed, size_t BOOST_TEST((pr.mono_at(g) == want)); BOOST_TEST(pr.k_of[g] == want.count()); BOOST_TEST(pr.phase_of[g] == expect_phase[g]); - BOOST_TEST(pr.sender_of[g] == expect_sender[g]); + BOOST_TEST(pr.sender_index(g).value == expect_sender[g]); + BOOST_TEST(pr.sender_slot(g) == window.base + expect_sender[g]); BOOST_TEST(pr.is_paired_at(g) == monoprop::is_paired(want)); std::vector key; @@ -263,12 +275,22 @@ BOOST_AUTO_TEST_CASE(sparse_resolve_probe_matches_narrow_positions) { BOOST_AUTO_TEST_CASE(sparse_resolve_probe_matches_wide_positions) { std::mt19937_64 rng(20260815); static_assert(sizeof(detail::OperatorIndex<250>::PosT) == 2, "this case exists to cover the wide store"); - check_probe_matches_the_queries<250>(rng, /*n_seed=*/60, /*n_query=*/140, /*rank_count=*/4, /*fused=*/false); + check_probe_matches_the_queries<250>(rng, + /*n_seed=*/60, + /*n_query=*/140, + /*rank_count=*/4, + /*fused=*/false, + /*window_base=*/16); } BOOST_AUTO_TEST_CASE(sparse_resolve_probe_matches_fused_layout) { std::mt19937_64 rng(20260816); - check_probe_matches_the_queries<250>(rng, /*n_seed=*/50, /*n_query=*/120, /*rank_count=*/2, /*fused=*/true); + check_probe_matches_the_queries<250>(rng, + /*n_seed=*/50, + /*n_query=*/120, + /*rank_count=*/2, + /*fused=*/true, + /*window_base=*/6); } BOOST_AUTO_TEST_CASE(sparse_resolve_probe_matches_single_sender) { diff --git a/docs/bibliography.bib b/docs/bibliography.bib index 290f83b3..79964d64 100644 --- a/docs/bibliography.bib +++ b/docs/bibliography.bib @@ -89,3 +89,16 @@ @ARTICLE{chakraborty2026scalablequantumcircuitgeneration primaryClass={quant-ph}, url={https://arxiv.org/abs/2603.23444}, } + +@ARTICLE{Broers2025-or, + title = "{Scalable Simulation of Quantum Many-Body Dynamics with Or-Represented Quantum Algebra}", + author = "Broers, Lukas and Sun, Rong-Yang and Yunoki, Seiji", + journal = "Phys. Rev. Applied", + volume = 26, + pages = "024046", + year = 2026, + url = "https://arxiv.org/abs/2506.13241", + archivePrefix = "arXiv", + primaryClass = "quant-ph", + eprint = "2506.13241" +} diff --git a/docs/content/docs/features/parallelism.mdx b/docs/content/docs/features/parallelism.mdx index e9e56ed5..068c9a06 100644 --- a/docs/content/docs/features/parallelism.mdx +++ b/docs/content/docs/features/parallelism.mdx @@ -30,6 +30,9 @@ whose partner lives in another partition are resolved through a per-gate exchang | `monoprop_NUM_THREADS` | one partition per physical core | Caps the number of partitions. Set it to run fewer partitions than cores. | | `monoprop_PARTITIONS` | `auto` | `auto` = one partition per core (capped by `monoprop_NUM_THREADS`); an integer `N` = exactly `N` partitions; `off` = one partition holding the whole operator. | +Three further variables control which MPI rank owns a term; they belong with the +distribution axis and are documented under [Rank routing](#rank-routing) below. + ```bash # Run 8 partitions instead of one-per-core: export monoprop_NUM_THREADS=8 @@ -61,6 +64,60 @@ partitioning into one flat world of `R × S` partitions (`R` ranks, `S` partitio 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$ bits — so that + +$$ +h(M \oplus G) = h(M) \oplus h(G). +$$ + +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. + +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 left as a diagnostic +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. + +Within a rank, partitions keep the full-avalanche `splitmix` hash: fanout across shared +memory is free, so only balance matters there. Routing is therefore two-level. Writing $q$ +for the `splitmix` hash of $M$, the partition is $q \bmod S$ and the rank is $h_d(M)$ with +$d = \log_2 R$ — every rank bit, so the fanout is 1: + +$$ +\mathrm{flat}(M) = h_d(M)\,S + (q \bmod S). +$$ + +Linear routing is a switch and not a dial: the rank takes all $\log_2 R$ bits from $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. $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 [@Broers2025-or]. 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. + +| Variable | Default | Meaning | +| --- | --- | --- | +| `monoprop_ROUTING` | `linear` | `splitmix` selects the dense all-to-all; `linear`, or unset, takes every rank bit from $h$ and requires a power-of-two $R$. Any other value is rejected at startup rather than silently defaulting. | +| `monoprop_ROUTE_SEED` | `6768574230969066775` | Decimal `uint64` from which every rank derives the same basis $\{v_i\}$ with no communication. The same value must reach every rank: a mismatch in either of these variables, or in the partition count, is caught by two allreduces at propagator construction and raised, because under linear routing it deadlocks the exchange instead of corrupting it. | + ### Single-node (`MPI.COMM_SELF`) ```python notest diff --git a/docs/content/docs/references.mdx b/docs/content/docs/references.mdx index 887489ee..410f3596 100644 --- a/docs/content/docs/references.mdx +++ b/docs/content/docs/references.mdx @@ -50,4 +50,11 @@ described in [@Miller2023-bonsai]. Other representative large-scale applications arXiv [quant-ph], 2026. arXiv:2605.04025 +
  • + L. Broers, R.-Y. Sun, and S. Yunoki, + "Scalable simulation of quantum many-body dynamics with or-represented + quantum algebra," + Phys. Rev. Applied, vol. 26, p. 024046, 2026. + arXiv:2506.13241 +