diff --git a/apps/BUILD b/apps/BUILD index 76c45e83a..c85f631dc 100644 --- a/apps/BUILD +++ b/apps/BUILD @@ -32,6 +32,16 @@ cc_binary( ], ) +cc_binary( + name = "qsim_gate_batch", + srcs = ["qsim_gate_batch.cc"], + copts = gsframe_copts, + deps = [ + "//lib:run_qsim_gate_batch", + "//lib:run_qsim_lib", + ], +) + cc_binary( name = "qsim_von_neumann", srcs = ["qsim_von_neumann.cc"], diff --git a/apps/Makefile b/apps/Makefile index 19ccbc422..bd731cdf8 100644 --- a/apps/Makefile +++ b/apps/Makefile @@ -1,6 +1,8 @@ CXX_TARGETS = $(shell find . -maxdepth 1 -name '*.cc') CXX_TARGETS := $(CXX_TARGETS:%.cc=%.x) +OUTPUT_DIR ?= . + CUDA_TARGETS = $(shell find . -maxdepth 1 -name '*cuda.cu') CUDA_TARGETS := $(CUDA_TARGETS:%cuda.cu=%cuda.x) @@ -16,6 +18,12 @@ HIP_TARGETS := $(HIP_TARGETS:%cuda.cu=%hip.x) .PHONY: qsim qsim: $(CXX_TARGETS) +.PHONY: gate-batch +gate-batch: + mkdir -p "$(OUTPUT_DIR)" + $(CXX) -o "$(OUTPUT_DIR)/qsim_gate_batch.x" qsim_gate_batch.cc $(CXXFLAGS) $(ARCHFLAGS) + $(CXX) -o "$(OUTPUT_DIR)/qsim_base.x" qsim_base.cc $(CXXFLAGS) $(ARCHFLAGS) + .PHONY: qsim-cuda qsim-cuda: $(CUDA_TARGETS) diff --git a/apps/qsim_gate_batch.cc b/apps/qsim_gate_batch.cc new file mode 100644 index 000000000..3ef01d9e9 --- /dev/null +++ b/apps/qsim_gate_batch.cc @@ -0,0 +1,124 @@ +// Benchmark app for the proposal-faithful gate-batch runner. +// -f is the per-batch max_fused_size (0 = apply +// raw gates without fusing, 2-3 = proposal's suggestion). + +#include + +#ifdef _OPENMP +# include +#endif + +#include +#include + +#include "../lib/circuit_qsim_parser.h" +#include "../lib/formux.h" +#include "../lib/fuser_mqubit.h" +#include "../lib/gates_qsim.h" +#include "../lib/io_file.h" +#include "../lib/operation.h" +#include "../lib/run_qsim_gate_batch.h" +#include "../lib/seqfor.h" +#include "../lib/simmux.h" +#include "../lib/util_cpu.h" + +struct Options { + std::string circuit_file; + unsigned maxtime = std::numeric_limits::max(); + unsigned seed = 1; + unsigned num_threads = 1; + unsigned max_fused_size = 3; + unsigned block_qubits = 19; + unsigned verbosity = 0; +}; + +Options GetOptions(int argc, char* argv[]) { + constexpr char usage[] = "usage:\n ./qsim_gate_batch -c circuit " + "-d maxtime -s seed -t threads " + "-f max_fused_size -l block_qubits -v verbosity\n"; + Options opt; + int k; + while ((k = getopt(argc, argv, "c:d:s:t:f:l:v:")) != -1) { + switch (k) { + case 'c': opt.circuit_file = optarg; break; + case 'd': opt.maxtime = std::atoi(optarg); break; + case 's': opt.seed = std::atoi(optarg); break; + case 't': opt.num_threads = std::atoi(optarg); break; + case 'f': opt.max_fused_size = std::atoi(optarg); break; + case 'l': opt.block_qubits = std::atoi(optarg); break; + case 'v': opt.verbosity = std::atoi(optarg); break; + default: qsim::IO::errorf(usage); exit(1); + } + } + if (opt.circuit_file.empty()) { + qsim::IO::errorf(usage); + exit(1); + } + return opt; +} + +template +void PrintAmplitudes(unsigned num_qubits, const StateSpace& state_space, + const QubitMappedState& state) { + static constexpr char const* bits[8] = { + "000", "001", "010", "011", "100", "101", "110", "111", + }; + uint64_t size = std::min(uint64_t{8}, uint64_t{1} << num_qubits); + unsigned s = 3 - std::min(unsigned{3}, num_qubits); + for (uint64_t i = 0; i < size; ++i) { + auto a = state.GetAmpl(state_space, i); + qsim::IO::messagef("%s:%16.8g%16.8g%16.8g\n", + bits[i] + s, std::real(a), std::imag(a), std::norm(a)); + } +} + +int main(int argc, char* argv[]) { + using namespace qsim; + + auto opt = GetOptions(argc, argv); + +#ifdef _OPENMP + omp_set_num_threads(opt.num_threads); +#endif + + Circuit> circuit; + if (!CircuitQsimParser::FromFile(opt.maxtime, opt.circuit_file, + circuit)) { + return 1; + } + + struct Factory { + Factory(unsigned num_threads) : num_threads(num_threads) {} + using Simulator = qsim::Simulator; + using StateSpace = Simulator::StateSpace; + StateSpace CreateStateSpace() const { return StateSpace(num_threads); } + Simulator CreateSimulator() const { return Simulator(num_threads); } + unsigned num_threads; + }; + + using Simulator = Factory::Simulator; + using StateSpace = Simulator::StateSpace; + using State = StateSpace::State; + using Fuser = MultiQubitGateFuser; + using SeqSimulator = qsim::Simulator; + using Runner = QSimGateBatchRunner; + + StateSpace state_space = Factory(opt.num_threads).CreateStateSpace(); + QubitMappedState state(state_space.Create(circuit.num_qubits)); + if (state_space.IsNull(state.state)) { + IO::errorf("not enough memory: is the number of qubits too large?\n"); + return 1; + } + state_space.SetStateZero(state.state); + + Runner::Parameter param; + param.max_fused_size = opt.max_fused_size; + param.block_qubits = opt.block_qubits; + param.verbosity = opt.verbosity; + + if (Runner::Run(param, Factory(opt.num_threads), circuit, state)) { + PrintAmplitudes(circuit.num_qubits, state_space, state); + } + + return 0; +} diff --git a/lib/BUILD b/lib/BUILD index ddc2ed2bc..785cef2bb 100644 --- a/lib/BUILD +++ b/lib/BUILD @@ -564,6 +564,21 @@ cc_library( ], ) +cc_library( + name = "run_qsim_gate_batch", + hdrs = [ + "qubit_mapped_state.h", + "qubit_remap.h", + "run_qsim_gate_batch.h", + ], + deps = [ + ":gate", + ":matrix", + ":operation_base", + ":util", + ], +) + cc_library( name = "run_qsimh", hdrs = ["run_qsimh.h"], diff --git a/lib/qubit_mapped_state.h b/lib/qubit_mapped_state.h new file mode 100644 index 000000000..1ee839a6b --- /dev/null +++ b/lib/qubit_mapped_state.h @@ -0,0 +1,122 @@ +// Copyright 2026 Google LLC. All Rights Reserved. +// +// 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 +// +// https://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. + +#ifndef QUBIT_MAPPED_STATE_H_ +#define QUBIT_MAPPED_STATE_H_ + +#include +#include +#include +#include + +#include "qubit_remap.h" + +namespace qsim { + +// Bidirectional logical<->physical qubit map; the two arrays are inverse +// permutations of each other at all times. "Physical position p" means +// amplitude-index bit p of the state storage. +class QubitLayout { + public: + explicit QubitLayout(unsigned num_qubits) + : logical_to_physical_(num_qubits), + physical_to_logical_(num_qubits) { + for (unsigned i = 0; i < num_qubits; ++i) { + logical_to_physical_[i] = i; + physical_to_logical_[i] = i; + } + } + + std::size_t NumQubits() const { return logical_to_physical_.size(); } + + unsigned PhysicalPositionOf(unsigned logical_qubit) const { + return logical_to_physical_[logical_qubit]; + } + + unsigned LogicalQubitAt(unsigned physical_position) const { + return physical_to_logical_[physical_position]; + } + + const std::vector& LogicalToPhysical() const { + return logical_to_physical_; + } + + uint64_t PhysicalIndexOf(uint64_t logical_index) const { + uint64_t physical_index = 0; + for (unsigned q = 0; q < NumQubits(); ++q) { + if ((logical_index >> q) & 1) { + physical_index |= uint64_t{1} << PhysicalPositionOf(q); + } + } + return physical_index; + } + + bool IsIdentityAt(unsigned position) const { + return physical_to_logical_[position] == position; + } + + void SwapPositions(unsigned position_a, unsigned position_b) { + const auto qubit_a = physical_to_logical_[position_a]; + const auto qubit_b = physical_to_logical_[position_b]; + std::swap(physical_to_logical_[position_a], + physical_to_logical_[position_b]); + std::swap(logical_to_physical_[qubit_a], + logical_to_physical_[qubit_b]); + } + + // Builds one disjoint restoration pass and updates this layout to match. + bool BuildIdentityRestorationSwaps(std::vector& swaps) { + swaps.clear(); + std::vector position_touched(NumQubits(), 0); + + for (unsigned position = 0; position < NumQubits(); ++position) { + if (IsIdentityAt(position)) continue; + const auto paired_position = PhysicalPositionOf(position); + if (position_touched[position] || + position_touched[paired_position]) { + continue; + } + + swaps.emplace_back(position, paired_position); + position_touched[position] = 1; + position_touched[paired_position] = 1; + SwapPositions(position, paired_position); + } + return !swaps.empty(); + } + + private: + std::vector logical_to_physical_; + std::vector physical_to_logical_; +}; + +// Owns physical state storage together with its logical-qubit interpretation. +// A newly created mapped state starts in canonical one-to-one qubit order. +template +struct QubitMappedState { + explicit QubitMappedState(State&& state_storage) + : state(std::move(state_storage)), layout(state.num_qubits()) {} + + template + auto GetAmpl(const StateSpace& state_space, uint64_t logical_index) const { + return state_space.GetAmpl(state, layout.PhysicalIndexOf(logical_index)); + } + + State state; + QubitLayout layout; +}; + +} // namespace qsim + +#endif // QUBIT_MAPPED_STATE_H_ diff --git a/lib/qubit_remap.h b/lib/qubit_remap.h new file mode 100644 index 000000000..c929686c9 --- /dev/null +++ b/lib/qubit_remap.h @@ -0,0 +1,201 @@ +// Qubit remapping for cache-blocked simulation (refactoring plan, Step 1). +// +// In-place, single-pass application of a set of DISJOINT qubit-position +// transpositions (an involution) to a state stored in a chunked SIMD +// layout. Disjoint transpositions commute, so one pass over the state +// (1 read + 1 write of the moved amplitudes) applies them all. +// +// Layout: a chunked SIMD state space stores amplitudes as chunks of +// +// chunk = [re0 .. re(N-1), im0 .. im(N-1)] N = 2^chunk_qubits +// +// i.e. the low `chunk_qubits` amplitude-index bits select the lane inside a +// chunk (NEON/SSE: chunk_qubits=2, N=4; AVX2: 3, N=8; AVX512: 4, N=16; a +// non-chunked/fully-interleaved layout is the degenerate chunk_qubits=0 +// case, N=1). All swapped positions must therefore be >= chunk_qubits: the +// permutation then moves whole chunks and never has to touch lanes. +// `chunk_qubits` is an explicit property of the caller's state space. +// +// Execution shape: for the lowest swapped amplitude bit b_min, amplitudes +// move as contiguous spans of 2^(b_min - chunk_qubits) chunks. The pass +// enumerates chunk-span starts, computes each span's partner by flipping +// the bit pairs whose two bits differ, and swaps the two spans (each span +// is swapped exactly once via the partner > self check; spans whose pair +// bits all match are fixed points and are skipped). +// +// Reading guide: ApplyBitPairSwaps (the only public entry point) is +// BitSwapPlan construction followed by a parallel loop of +// FirstChunkOfSpan -> PartnerChunk -> SwapChunkSpans. + +#ifndef QUBIT_REMAP_H_ +#define QUBIT_REMAP_H_ + +#include +#include +#include +#include +#include + +namespace qsim { + +// A transposition of two amplitude-index bit positions. +using QubitSwap = std::pair; + +namespace remap_internal { + +// One amplitude-bit swap translated to chunk-index bit positions. +class ChunkBitSwap { + public: + ChunkBitSwap(unsigned first_amplitude_bit, unsigned second_amplitude_bit, + unsigned chunk_qubits) { + assert(first_amplitude_bit >= chunk_qubits); + assert(second_amplitude_bit >= chunk_qubits); + lower_chunk_bit_ = + std::min(first_amplitude_bit, second_amplitude_bit) - chunk_qubits; + upper_chunk_bit_ = + std::max(first_amplitude_bit, second_amplitude_bit) - chunk_qubits; + flip_mask_ = (uint64_t{1} << lower_chunk_bit_) | + (uint64_t{1} << upper_chunk_bit_); + } + + unsigned LowerChunkBit() const { return lower_chunk_bit_; } + + // A bit-position swap moves the chunk only when the two bits differ. + bool ChangesChunkIndex(uint64_t chunk_index) const { + return ((chunk_index >> lower_chunk_bit_) & 1) != + ((chunk_index >> upper_chunk_bit_) & 1); + } + + uint64_t FlipMask() const { return flip_mask_; } + + private: + unsigned lower_chunk_bit_; + unsigned upper_chunk_bit_; + uint64_t flip_mask_; +}; + +// Precomputed chunk geometry for applying a set of disjoint qubit swaps. +class BitSwapPlan { + public: + BitSwapPlan(unsigned num_qubits, unsigned chunk_qubits, + const std::vector& qubit_swaps) + : floats_per_chunk_(uint64_t{2} << chunk_qubits), + chunk_index_bits_(num_qubits - chunk_qubits), + chunk_span_bits_(chunk_index_bits_) { + assert(num_qubits >= chunk_qubits); +#ifndef NDEBUG + std::vector is_qubit_swapped(num_qubits, 0); +#endif + chunk_bit_swaps_.reserve(qubit_swaps.size()); + for (const QubitSwap& qubit_swap : qubit_swaps) { + assert(qubit_swap.first < num_qubits); + assert(qubit_swap.second < num_qubits); + assert(qubit_swap.first != qubit_swap.second); +#ifndef NDEBUG + assert(!is_qubit_swapped[qubit_swap.first]); + assert(!is_qubit_swapped[qubit_swap.second]); + is_qubit_swapped[qubit_swap.first] = 1; + is_qubit_swapped[qubit_swap.second] = 1; +#endif + AddQubitSwap(qubit_swap, chunk_qubits); + } + } + + // Floats in one chunk (2^chunk_qubits amplitudes, 2 floats each). + uint64_t FloatsPerChunk() const { return floats_per_chunk_; } + + // Amplitudes move in contiguous spans of 2^chunk_span_bits chunks. + uint64_t FloatsPerChunkSpan() const { + return floats_per_chunk_ << chunk_span_bits_; + } + + int64_t NumChunkSpans() const { + return int64_t{1} << (chunk_index_bits_ - chunk_span_bits_); + } + + uint64_t FirstChunkOfSpan(int64_t span_index) const { + return uint64_t(span_index) << chunk_span_bits_; + } + + // Applies every chunk-bit swap to find the involution partner. A chunk + // whose paired bits all match maps to itself. + uint64_t PartnerChunk(uint64_t chunk_index) const { + auto partner_chunk = chunk_index; + for (const ChunkBitSwap& bit_swap : chunk_bit_swaps_) { + if (bit_swap.ChangesChunkIndex(chunk_index)) { + partner_chunk ^= bit_swap.FlipMask(); + } + } + return partner_chunk; + } + + private: + void AddQubitSwap(const QubitSwap& qubit_swap, unsigned chunk_qubits) { + chunk_bit_swaps_.emplace_back(qubit_swap.first, qubit_swap.second, + chunk_qubits); + chunk_span_bits_ = std::min( + chunk_span_bits_, chunk_bit_swaps_.back().LowerChunkBit()); + } + + uint64_t floats_per_chunk_; + std::vector chunk_bit_swaps_; + unsigned chunk_index_bits_; + unsigned chunk_span_bits_; +}; + +// First float of the chunk with this index. +inline float* FirstFloatOfChunk(float* state, uint64_t floats_per_chunk, + uint64_t chunk_index) { + return state + floats_per_chunk * chunk_index; +} + +// Swap two contiguous chunk spans. Compilers vectorize this to SIMD-width +// loads and stores; the loop is generally memory-bandwidth-bound. +inline void SwapChunkSpans(float* __restrict first_span, + float* __restrict second_span, + uint64_t num_floats) { + for (uint64_t i = 0; i < num_floats; ++i) { + float temporary = first_span[i]; + first_span[i] = second_span[i]; + second_span[i] = temporary; + } +} + +} // namespace remap_internal + +// Applies all `qubit_swaps` transpositions of amplitude-bit positions to the +// state, in place, in a single pass. The pairs must be disjoint and every +// position must be >= chunk_qubits (see layout comment above). +inline void ApplyBitPairSwaps( + float* state, unsigned num_qubits, unsigned chunk_qubits, + const std::vector& qubit_swaps) { + namespace ri = remap_internal; + + if (qubit_swaps.empty()) { + return; + } + + const auto plan = ri::BitSwapPlan(num_qubits, chunk_qubits, qubit_swaps); + const auto floats_per_chunk = plan.FloatsPerChunk(); + const auto floats_per_chunk_span = plan.FloatsPerChunkSpan(); + const auto num_chunk_spans = plan.NumChunkSpans(); + +#pragma omp parallel for schedule(dynamic, 32) + for (int64_t span_index = 0; span_index < num_chunk_spans; ++span_index) { + const auto first_chunk = plan.FirstChunkOfSpan(span_index); + const auto partner_chunk = plan.PartnerChunk(first_chunk); + + // Each moved pair is visited from both sides; act on one of them + // (fixed points have partner == chunk and fall through). + if (partner_chunk > first_chunk) { + ri::SwapChunkSpans( + ri::FirstFloatOfChunk(state, floats_per_chunk, first_chunk), + ri::FirstFloatOfChunk(state, floats_per_chunk, partner_chunk), + floats_per_chunk_span); + } + } +} + +} // namespace qsim + +#endif // QUBIT_REMAP_H_ diff --git a/lib/run_qsim_gate_batch.h b/lib/run_qsim_gate_batch.h new file mode 100644 index 000000000..e4899b1e5 --- /dev/null +++ b/lib/run_qsim_gate_batch.h @@ -0,0 +1,1008 @@ +// Gate-batched cache-local simulation runner. Follows the original +// proposal as closely as possible (see original_proposal.md). +// +// Data hierarchy: +// Circuit: the full program, as an ordered list of raw gates. +// State block: a cache-sized slice of 2^block_qubits amplitudes; one +// gate batch runs against every state block. +// SIMD chunk: the vectorized unit inside a state block - 2^chunk_qubits +// complex amplitudes, architecture-dependent width (NEON/SSE: +// 4, AVX2: 8, AVX512: 16). Remapping moves whole chunks. +// +// A gate batch ties this together: the planner selects a compatible set +// of circuit gates, remaps their logical qubits into physical block +// positions, fuses them, and runs the result on every state block. +// +// Algorithm: +// Move the most-used logical qubits into the low block up front +// (PlaceHotQubitsInFixedZone) - free, since the all-zero initial state +// is permutation-invariant. +// +// while gates remain: +// Select the next gate batch (GateBatchPlanner::PlanNextGateBatch): +// score a bounded set of candidate qubit sets - the greedy +// first-fit scan, the qubits already resident, and a few sets +// seeded from upcoming gates - and keep the best. +// Remap the batch's qubits into physical positions [0, block_qubits) +// with one pass of disjoint transpositions (BuildBatchSwapPairs, +// ApplySwapsToState). The low chunk_qubits positions are pinned and +// never swapped; every other qubit draws from a bounded remap budget +// (GateBatchPlanner::RemapSlotCapacity). +// Fuse the batch's gates on physical qubits with the standard fuser +// (FuseBatchGates); max_fused_size == 0 disables fusion. +// Execute every fused gate on every state block, blocks in parallel, +// each with its own sequential simulator (ExecuteGateBatchOnBlocks). +// +// Restore identity qubit order with a final sequence of disjoint- +// transposition passes (RestoreIdentityQubitOrder). +// +// Gates are scheduled raw: there is no global pre-fusion, 1q-chain +// peephole, or diagonal-gate phase kernel/bundle. All fusion happens +// inside the per-batch fuser above. + +#ifndef RUN_QSIM_GATE_BATCH_H_ +#define RUN_QSIM_GATE_BATCH_H_ + +#include +#include +#include +#include +#include +#include + +#include "gate.h" +#include "qubit_mapped_state.h" +#include "matrix.h" +#include "operation_base.h" +#include "qubit_remap.h" +#include "util.h" + +namespace qsim { + +namespace gate_batch_internal { + +// Internal scheduling representation of one raw circuit gate. +template +struct PendingGate { + PendingGate(std::vector qubits, Matrix gate_matrix) + : logical_qubits(std::move(qubits)), + matrix(std::move(gate_matrix)) {} + + bool IsPending() const { return !applied_; } + void MarkApplied() { applied_ = true; } + std::size_t Arity() const { return logical_qubits.size(); } + + std::vector logical_qubits; // ascending + Matrix matrix; + + private: + bool applied_ = false; +}; + +// A fused (or passthrough) gate ready to execute inside a state block. +template +struct ExecutableGate { + std::vector physical_qubits; // ascending + Matrix matrix; +}; + +// Geometry of the cache-blocked state: an n-qubit state is viewed as +// 2^(n - L) contiguous blocks of 2^L amplitudes. +template +struct BlockPartition { + BlockPartition(unsigned state_qubits, unsigned requested_block_qubits) + : num_state_qubits(state_qubits), + block_qubits(std::min(requested_block_qubits, state_qubits)), + floats_per_block(StateSpace::MinSize(block_qubits)), + num_blocks(int64_t{1} << (state_qubits - block_qubits)) {} + + unsigned num_state_qubits; // n + unsigned block_qubits; // L + uint64_t floats_per_block; // state floats per block (SIMD layout) + int64_t num_blocks; // 2^(n - L) +}; + +// Counters and phase timings accumulated across gate batches. Gate batches +// and identity-restoration passes are counted separately because only gate +// batches make circuit progress; restore passes are pure overhead. Timers +// are filled only at verbosity > 1. +struct SimulationStats { + unsigned num_gate_batches = 0; + unsigned num_restore_passes = 0; + unsigned num_swaps = 0; + unsigned num_executable_gates = 0; + double plan_seconds = 0.0; + double swap_seconds = 0.0; + double fuse_seconds = 0.0; + double gate_seconds = 0.0; +}; + +// Buffers reused by every gate batch. +template +struct GateBatchWorkspace { + std::vector swap_pairs; + std::vector> executable_gates; +}; + +// Logical qubits selected for one gate batch, with eviction-aware remap-slot +// accounting. Logical qubits resident below the eviction floor join for free. +struct GateBatchQubitSet { + bool Contains(unsigned logical_qubit) const { + return contains_qubit[logical_qubit]; + } + + std::vector contains_qubit; + unsigned remap_slots_used = 0; + unsigned remap_slot_capacity = 0; +}; + +// The outcome of planning one gate batch, built without mutating gates, +// layout, or state. +struct GateBatchPlan { + bool HasGates() const { return !gate_indices.empty(); } + std::size_t NumGates() const { return gate_indices.size(); } + bool UsesQubit(unsigned logical_qubit) const { + return uses_qubit[logical_qubit]; + } + bool operator<(const GateBatchPlan& other) const { + return score < other.score; + } + + std::vector gate_indices; + std::vector uses_qubit; + unsigned required_swaps = 0; + double score = 0.0; +}; + +} // namespace gate_batch_internal + +template +class QSimGateBatchRunner final { + public: + using StateSpace = typename Factory::StateSpace; + using State = typename StateSpace::State; + using QubitMappedState = qsim::QubitMappedState; + using fp_type = typename StateSpace::fp_type; + using SeqStateSpace = typename SeqSimulator::StateSpace; + + static_assert(std::is_same_v, + "QSimGateBatchRunner requires a float state space."); + static_assert(std::is_same_v, + "State spaces must use the same floating-point type."); + static_assert(StateSpace::kChunkQubits == SeqStateSpace::kChunkQubits, + "State spaces must use the same SIMD chunk layout."); + + struct Parameter : public Fuser::Parameter { + Parameter() { this->max_fused_size = 3; } + + unsigned block_qubits = 19; + }; + + template + static bool Run(const Parameter& param, const Factory& factory, + const Circuit& circuit, State& state) { + // Keep the conventional qsim static entry point while putting all + // mutable execution state in a private, one-shot runner instance. + (void) factory; // Present for API compatibility with other runners. + QubitLayout layout(circuit.num_qubits); + QSimGateBatchRunner runner(param, circuit.num_qubits, state, layout); + return runner.SimulateCircuit(circuit, true); + } + + // Runs directly on a mapped state and leaves its physical storage in the + // returned logical-to-physical layout, avoiding final restoration. + template + static bool Run(const Parameter& param, const Factory& factory, + const Circuit& circuit, QubitMappedState& state) { + (void) factory; + QSimGateBatchRunner runner(param, circuit.num_qubits, state.state, + state.layout); + return runner.SimulateCircuit(circuit, false); + } + + private: + using PendingGate = gate_batch_internal::PendingGate; + using ExecutableGate = gate_batch_internal::ExecutableGate; + using BlockPartition = gate_batch_internal::BlockPartition; + using SimulationStats = gate_batch_internal::SimulationStats; + using GateBatchWorkspace = + gate_batch_internal::GateBatchWorkspace; + using GateBatchQubitSet = gate_batch_internal::GateBatchQubitSet; + using GateBatchPlan = gate_batch_internal::GateBatchPlan; + + QSimGateBatchRunner(const Parameter& param, unsigned num_qubits, + State& state, QubitLayout& layout) + : param_(param), + partition_(num_qubits, param.block_qubits), + state_data_(state.get()), + chunk_qubits_(StateSpace::kChunkQubits), + seq_sim_(1), + layout_(layout), + gate_batch_planner_(partition_.num_state_qubits, + partition_.block_qubits, chunk_qubits_) { + assert(partition_.block_qubits >= chunk_qubits_); + assert(layout_.NumQubits() == partition_.num_state_qubits); + } + + template + bool SimulateCircuit(const Circuit& circuit, bool restore_qubit_order) { + using Op = typename std::decay_t::value_type; + + const auto prepare_start = StartPhaseTimer(); + if (!PreparePendingGates(circuit)) return false; + LogPreparationTime(prepare_start); + + const auto simulation_start = param_.verbosity > 0 ? GetTime() : 0.0; + PlaceHotQubitsInFixedZone(); + + // Op depends on Circuit, so this buffer remains local rather than + // forcing the whole runner type to acquire another template parameter. + std::vector batch_operations; + + std::size_t applied_count = 0; + while (applied_count < pending_gates_.size()) { + const auto num_applied = ExecuteNextGateBatch(batch_operations); + if (num_applied == 0) return false; // fuser error or stuck schedule + applied_count += num_applied; + } + + if (restore_qubit_order) RestoreIdentityQubitOrder(); + LogSimulationSummary(simulation_start); + return true; + } + + // ======== Phase timing ======== + + double StartPhaseTimer() const { + return param_.verbosity > 1 ? GetTime() : 0.0; + } + + void AccumulatePhaseSeconds(double start, double& seconds) const { + if (param_.verbosity > 1) seconds += GetTime() - start; + } + + void LogPreparationTime(double prepare_start) const { + if (param_.verbosity <= 1) return; + IO::messagef("prepare time is %g seconds.\n", + GetTime() - prepare_start); + } + + // ======== Shared gate utilities ======== + + // Sorts a gate's qubits ascending, permuting `matrix` to match. + static void NormalizeGateQubitOrder(std::vector& qubits, + Matrix& matrix) { + if (qubits.size() < 2) return; + auto permutation = NormalToGateOrderPermutation(qubits); + if (!permutation.empty()) { + MatrixShuffle(permutation, unsigned(qubits.size()), matrix); + std::sort(qubits.begin(), qubits.end()); + } + } + + // ======== Pending-gate preparation ======== + + // The proposal starts from the raw circuit: every gate becomes one + // pending gate, verbatim (qubit order normalized). All fusion is deferred + // to the per-batch fuser. Controlled gates, measurements, and channels + // are not supported. + template + bool PreparePendingGates(const Circuit& circuit) { + pending_gates_.reserve(circuit.ops.size()); + + for (const auto& operation : circuit.ops) { + const auto* raw_gate = + OpGetAlternative>(operation); + if (raw_gate == nullptr) { + IO::errorf("qsim_gate_batch: unsupported operation " + "(controlled gate or measurement).\n"); + return false; + } + + auto logical_qubits = raw_gate->qubits; + auto matrix = raw_gate->matrix; + NormalizeGateQubitOrder(logical_qubits, matrix); + pending_gates_.emplace_back(std::move(logical_qubits), + std::move(matrix)); + } + return true; + } + + // ======== Gate-batch planning ======== + + // Chooses the gates for the next gate batch. First-fit selection is greedy + // maximal, not maximum: an early gate can fill the qubit set and shut out + // a larger family of later gates. Rather than solve that exactly (it is + // combinatorial), PlanNextGateBatch tries a bounded set of candidate + // qubit sets and keeps the best-scoring plan: + // - the first-fit greedy set (baseline: never worse than a single + // greedy scan, because re-evaluating a finished set can only + // admit more gates); + // - the currently resident set (a zero-swap batch); + // - a set seeded by each of the next few pending gates, for when + // the first pending gate is the one poisoning the set. + // PlanNextGateBatch reads gates and layout but never modifies them. + class GateBatchPlanner { + public: + GateBatchPlanner(unsigned num_state_qubits, unsigned block_qubits, + unsigned chunk_qubits) + : num_logical_qubits_(num_state_qubits), + block_qubits_(block_qubits), + chunk_qubits_(chunk_qubits), + is_qubit_blocked_(num_state_qubits) {} + + GateBatchPlan PlanNextGateBatch(const std::vector& gates, + const QubitLayout& layout) { + std::vector candidate_plans; + candidate_plans.reserve(2 + kMaxGateSeeds); + + // First-fit greedy candidate. + candidate_plans.push_back(PlanGateBatchFromQubitSet( + gates, layout, MakeEmptyQubitSet())); + + // Zero-swap resident candidate. + candidate_plans.push_back(PlanGateBatchFromQubitSet( + gates, layout, MakeResidentQubitSet(layout))); + + for (std::size_t seed_idx : CollectSeedGateIndices(gates)) { + auto qubit_set = MakeEmptyQubitSet(); + if (!TryAdmitQubits(gates[seed_idx].logical_qubits, + layout, qubit_set)) { + continue; + } + // Candidate seeded by a later pending gate. + candidate_plans.push_back(PlanGateBatchFromQubitSet( + gates, layout, std::move(qubit_set))); + } + + // max_element returns the first maximum, so ties preserve candidate + // priority: greedy, resident, then seeded plans in circuit order. + return std::move(*std::max_element( + candidate_plans.begin(), candidate_plans.end())); + } + + unsigned EvictionFloor() const { return ComputeEvictionFloor(); } + + private: + // How many pending gates beyond the first (which already leads the + // greedy baseline) get to seed their own candidate qubit set. + static constexpr unsigned kMaxGateSeeds = 4; + + // Keep enough remap positions for an ordinary two-qubit gate even when + // the block is too small to reach kMinEvictionFloor. + static constexpr unsigned kMinRemapSlots = 2; + + // Lowest eviction position allowed in a multi-block gate batch: swap-pass + // spans are 2^(floor - chunk_qubits) chunks, so floor 11 keeps every + // span at least 16 KiB for float states and at streaming bandwidth. + // See RemapSlotCapacity. + static constexpr unsigned kMinEvictionFloor = 11; + + // Starting any remapping incurs a full-state pass. Additional pairs share + // that pass and therefore carry a smaller marginal cost. With these + // weights, eight pairs offset one additional gate. + static constexpr double kSwapPassCost = 0.5; + static constexpr double kSwapPairCost = 1.0 / 16.0; + + GateBatchQubitSet MakeEmptyQubitSet() const { + GateBatchQubitSet qubit_set; + qubit_set.contains_qubit.assign(num_logical_qubits_, 0); + qubit_set.remap_slot_capacity = RemapSlotCapacity(); + return qubit_set; + } + + // Positions below the eviction floor are fixed residents and join a set + // for free. Positions at or above it consume one remap slot: a low + // resident protects a potential victim, while a high qubit requires a + // victim. Keeping the floor at kMinEvictionFloor makes every remap span + // at least 16 KiB for the default L=19. Small blocks retain at least two + // remap slots so an ordinary two-qubit gate can make progress. + unsigned RemapSlotCapacity() const { + return block_qubits_ - ComputeEvictionFloor(); + } + + unsigned ComputeEvictionFloor() const { + if (block_qubits_ == num_logical_qubits_) return block_qubits_; + + const auto full = block_qubits_ > chunk_qubits_ + ? block_qubits_ - chunk_qubits_ + : 0u; + const auto depth_capped = block_qubits_ > kMinEvictionFloor + ? block_qubits_ - kMinEvictionFloor + : 0u; + const auto capacity = + std::min(full, std::max(depth_capped, kMinRemapSlots)); + return block_qubits_ - capacity; + } + + // The zero-swap candidate: exactly the qubits already resident in + // the low block. Its occupancy may exceed the depth-capped capacity; + // that is safe because none of these qubits needs a swap (no + // evictions happen), and growth beyond them is still capacity-bound. + GateBatchQubitSet MakeResidentQubitSet( + const QubitLayout& layout) const { + auto qubit_set = MakeEmptyQubitSet(); + for (unsigned q = 0; q < num_logical_qubits_; ++q) { + if (layout.PhysicalPositionOf(q) < block_qubits_) { + qubit_set.remap_slots_used += + AdditionalRemapSlotsFor(q, qubit_set, layout); + qubit_set.contains_qubit[q] = 1; + } + } + return qubit_set; + } + + // Gate indices of the next kMaxGateSeeds pending gates after the + // first pending one. + std::vector CollectSeedGateIndices( + const std::vector& gates) const { + std::vector seeds; + bool skipped_first_pending = false; + for (std::size_t idx = 0; + idx < gates.size() && seeds.size() < kMaxGateSeeds; ++idx) { + if (!gates[idx].IsPending()) continue; + if (!skipped_first_pending) { + skipped_first_pending = true; + continue; + } + seeds.push_back(idx); + } + return seeds; + } + + GateBatchPlan PlanGateBatchFromQubitSet( + const std::vector& gates, const QubitLayout& layout, + GateBatchQubitSet qubit_set) { + GrowQubitSetGreedily(gates, layout, qubit_set); + return EvaluateQubitSet(gates, layout, qubit_set); + } + + // First-fit growth: scan pending gates in circuit order, admitting + // each gate's qubits when they fit and blocking them otherwise (the + // same causal rule EvaluateQubitSet applies to the finished set). + void GrowQubitSetGreedily(const std::vector& gates, + const QubitLayout& layout, + GateBatchQubitSet& qubit_set) { + ClearBlockedQubits(); + for (const PendingGate& gate : gates) { + if (!gate.IsPending()) continue; + if (AnyQubitBlocked(gate.logical_qubits) || + !TryAdmitQubits(gate.logical_qubits, layout, qubit_set)) { + BlockQubits(gate.logical_qubits); + } + } + } + + // Exact evaluation of a fixed qubit set: one causal-blocking scan over + // the pending gates collects every gate the set can execute, then the + // plan is scored. + GateBatchPlan EvaluateQubitSet( + const std::vector& gates, const QubitLayout& layout, + const GateBatchQubitSet& qubit_set) { + GateBatchPlan plan; + plan.uses_qubit.assign(num_logical_qubits_, 0); + ClearBlockedQubits(); + + for (std::size_t idx = 0; idx < gates.size(); ++idx) { + const PendingGate& gate = gates[idx]; + if (!gate.IsPending()) continue; + if (!AnyQubitBlocked(gate.logical_qubits) && + QubitSetContainsAll(gate.logical_qubits, qubit_set)) { + plan.gate_indices.push_back(idx); + for (unsigned q : gate.logical_qubits) { + plan.uses_qubit[q] = 1; + } + } else { + BlockQubits(gate.logical_qubits); + } + } + + plan.required_swaps = CountSwapsNeeded(plan, layout); + const auto swap_cost = + plan.required_swaps == 0 + ? 0.0 + : kSwapPassCost + kSwapPairCost * plan.required_swaps; + plan.score = double(plan.NumGates()) - swap_cost; + return plan; + } + + // Remap-slot cost of admitting q. Fixed residents below the eviction + // floor cost nothing. A resident in the eviction zone consumes a victim + // position by protecting it, and a high qubit consumes one by entering. + unsigned AdditionalRemapSlotsFor( + unsigned q, const GateBatchQubitSet& qubit_set, + const QubitLayout& layout) const { + if (qubit_set.Contains(q)) return 0; + return layout.PhysicalPositionOf(q) < ComputeEvictionFloor() ? 0 : 1; + } + + // Admits all of a gate's qubits into the set, or none of them. + bool TryAdmitQubits(const std::vector& qubits, + const QubitLayout& layout, + GateBatchQubitSet& qubit_set) const { + unsigned additional_remap_slots = 0; + for (unsigned q : qubits) { + additional_remap_slots += + AdditionalRemapSlotsFor(q, qubit_set, layout); + } + if (qubit_set.remap_slots_used + additional_remap_slots > + qubit_set.remap_slot_capacity) { + return false; + } + for (unsigned q : qubits) qubit_set.contains_qubit[q] = 1; + qubit_set.remap_slots_used += additional_remap_slots; + return true; + } + + static bool QubitSetContainsAll(const std::vector& qubits, + const GateBatchQubitSet& qubit_set) { + return std::all_of( + qubits.begin(), qubits.end(), + [&qubit_set](unsigned q) { return qubit_set.Contains(q); }); + } + + unsigned CountSwapsNeeded(const GateBatchPlan& plan, + const QubitLayout& layout) const { + unsigned count = 0; + for (unsigned q = 0; q < num_logical_qubits_; ++q) { + if (plan.UsesQubit(q) && + layout.PhysicalPositionOf(q) >= block_qubits_) { + ++count; + } + } + return count; + } + + void ClearBlockedQubits() { + std::fill(is_qubit_blocked_.begin(), is_qubit_blocked_.end(), 0); + } + + bool AnyQubitBlocked(const std::vector& qubits) const { + return std::any_of( + qubits.begin(), qubits.end(), + [this](unsigned q) { return is_qubit_blocked_[q]; }); + } + + void BlockQubits(const std::vector& qubits) { + for (unsigned q : qubits) is_qubit_blocked_[q] = 1; + } + + unsigned num_logical_qubits_; + unsigned block_qubits_; + unsigned chunk_qubits_; + std::vector is_qubit_blocked_; + }; + + // ======== Gate-batch pipeline ======== + + // Executes one complete gate batch following the proposal's loop body. + // Returns the number of gates applied; 0 signals failure (fuser error, or + // a schedule that cannot make progress). Planning is read-only; gates are + // marked applied only after fusion succeeds, so a failed batch never + // corrupts the schedule. + template + std::size_t ExecuteNextGateBatch(std::vector& batch_operations) { + const auto block_qubits = partition_.block_qubits; + + // pick_maximum_number_of_gates_acting_on_low_qubits + const auto plan_start = StartPhaseTimer(); + const auto plan = gate_batch_planner_.PlanNextGateBatch( + pending_gates_, layout_); + AccumulatePhaseSeconds(plan_start, simulation_stats_.plan_seconds); + if (!plan.HasGates()) { + IO::errorf("qsim_gate_batch: a gate does not fit the qubit set " + "of %u block qubits; use a larger block_qubits.\n", + block_qubits); + return 0; + } + LogPlannedGateBatch(plan); + + // select_qubits_for_swap + swap_low_and_high_qubits + BuildBatchSwapPairs(plan); + LogGateBatchSwapSummary(plan); + LogAppliedRemap(); + ApplySwapsToState(gate_batch_workspace_.swap_pairs); + + // fused_low_gates = fuse(low_gates) + const auto fuse_start = StartPhaseTimer(); + BuildBatchOperations(plan, batch_operations); + if (!FuseBatchGates(batch_operations)) { + return 0; + } + AccumulatePhaseSeconds(fuse_start, simulation_stats_.fuse_seconds); + simulation_stats_.num_executable_gates += + unsigned(gate_batch_workspace_.executable_gates.size()); + MarkGatesApplied(plan); + + // for i in 0..(2^num_high_qubits): apply all fused gates to block i + const auto gates_start = StartPhaseTimer(); + ExecuteGateBatchOnBlocks(gate_batch_workspace_.executable_gates, + state_data_, partition_, seq_sim_); + AccumulatePhaseSeconds(gates_start, simulation_stats_.gate_seconds); + + ++simulation_stats_.num_gate_batches; + return plan.NumGates(); + } + + // ======== Gate-batch diagnostics ======== + + // The last pair holds the lowest eviction position; a low floor means + // short scattered spans in the swap pass (see qubit_remap.h). + void LogGateBatchSwapSummary(const GateBatchPlan& plan) const { + if (param_.verbosity <= 2 || + gate_batch_workspace_.swap_pairs.empty()) { + return; + } + + IO::messagef("gate batch %u: %u gates, %u swaps, evict floor %u\n", + simulation_stats_.num_gate_batches, + unsigned(plan.NumGates()), + unsigned(gate_batch_workspace_.swap_pairs.size()), + gate_batch_workspace_.swap_pairs.back().first); + } + + // The planned gate batch before any remapping: every planned gate with its + // logical qubits, then every distinct qubit the batch uses with its + // current physical position; '*' marks qubits outside the block that + // the remap is about to swap in. + void LogPlannedGateBatch(const GateBatchPlan& plan) const { + if (param_.verbosity <= 3) return; + + IO::messagef("gate batch %u plan: %u gates, %u swaps needed\n gates:", + simulation_stats_.num_gate_batches, + unsigned(plan.NumGates()), + plan.required_swaps); + for (std::size_t idx : plan.gate_indices) { + const PendingGate& gate = pending_gates_[idx]; + IO::messagef(" #%u[", unsigned(idx)); + for (std::size_t i = 0; i < gate.Arity(); ++i) { + IO::messagef(i == 0 ? "q%u" : ",q%u", + gate.logical_qubits[i]); + } + IO::messagef("]"); + } + IO::messagef("\n qubits:"); + for (unsigned q = 0; q < unsigned(plan.uses_qubit.size()); ++q) { + if (!plan.UsesQubit(q)) continue; + const auto position = layout_.PhysicalPositionOf(q); + IO::messagef(" q%u@p%u%s", q, position, + position >= partition_.block_qubits ? "*" : ""); + } + IO::messagef("\n"); + } + + // The remap just applied (the layout is already updated): each + // transposition as "incoming qubit, its new<-old position, outgoing + // qubit", then the low-block layout the batch's gates will use. + void LogAppliedRemap() const { + if (param_.verbosity <= 3) return; + + IO::messagef(" swaps:"); + if (gate_batch_workspace_.swap_pairs.empty()) IO::messagef(" none"); + for (const QubitSwap& pair : gate_batch_workspace_.swap_pairs) { + IO::messagef(" [q%u in p%u<-p%u, q%u out]", + layout_.LogicalQubitAt(pair.first), pair.first, + pair.second, layout_.LogicalQubitAt(pair.second)); + } + IO::messagef("\n block:"); + for (unsigned p = 0; p < partition_.block_qubits; ++p) { + IO::messagef(" q%u", layout_.LogicalQubitAt(p)); + } + IO::messagef("\n"); + } + + // Extends the layout so every qubit the plan uses lands in physical + // positions [0, L), recording the transpositions to apply. Eviction + // walks down from L-1, skipping positions that hold used qubits; the + // planner's slot accounting guarantees it never reaches the pinned + // lane positions. + void BuildBatchSwapPairs(const GateBatchPlan& plan) { + gate_batch_workspace_.swap_pairs.clear(); + unsigned evict_position = partition_.block_qubits - 1; + + for (unsigned q = 0; q < unsigned(plan.uses_qubit.size()); ++q) { + if (!plan.UsesQubit(q)) continue; + const auto current_position = layout_.PhysicalPositionOf(q); + if (current_position < partition_.block_qubits) continue; + + while (plan.UsesQubit(layout_.LogicalQubitAt(evict_position))) { + --evict_position; + } + gate_batch_workspace_.swap_pairs.emplace_back(evict_position, + current_position); + layout_.SwapPositions(evict_position, current_position); + --evict_position; + } + } + + // Applies the transpositions to the state in one involution pass. + void ApplySwapsToState(const std::vector& swap_pairs) { + const auto swap_start = StartPhaseTimer(); + ApplyBitPairSwaps(state_data_, partition_.num_state_qubits, chunk_qubits_, + swap_pairs); + AccumulatePhaseSeconds(swap_start, simulation_stats_.swap_seconds); + simulation_stats_.num_swaps += unsigned(swap_pairs.size()); + } + + // Scores qubits by gate participation, weighted by gate arity. + std::vector ComputeQubitUsageScores() const { + std::vector scores(partition_.num_state_qubits, 0); + for (const PendingGate& gate : pending_gates_) { + const auto weight = gate.Arity(); + for (unsigned q : gate.logical_qubits) scores[q] += weight; + } + return scores; + } + + // Returns the count highest-scoring qubits in deterministic rank order. + static std::vector SelectHighestScoringQubits( + const std::vector& scores, unsigned first_qubit, + unsigned count) { + std::vector candidates; + candidates.reserve(scores.size() - first_qubit); + for (unsigned q = first_qubit; q < scores.size(); ++q) { + candidates.push_back(q); + } + assert(count <= candidates.size()); + std::partial_sort(candidates.begin(), candidates.begin() + count, + candidates.end(), + [&scores](unsigned a, unsigned b) { + return scores[a] != scores[b] + ? scores[a] > scores[b] : a < b; + }); + candidates.resize(count); + return candidates; + } + + // Builds disjoint swaps that put the selected logical qubits in the fixed + // physical zone and updates layout to match the state permutation. + std::vector + BuildFixedZonePlacementSwaps(const std::vector& fixed_qubits, + unsigned eviction_floor) { + std::vector is_fixed_qubit(layout_.NumQubits(), 0); + for (unsigned q : fixed_qubits) is_fixed_qubit[q] = 1; + + std::vector victim_positions; + std::vector incoming_qubits; + victim_positions.reserve(fixed_qubits.size()); + incoming_qubits.reserve(fixed_qubits.size()); + + for (unsigned p = chunk_qubits_; p < eviction_floor; ++p) { + if (!is_fixed_qubit[layout_.LogicalQubitAt(p)]) { + victim_positions.push_back(p); + } + } + for (unsigned q : fixed_qubits) { + if (layout_.PhysicalPositionOf(q) >= eviction_floor) { + incoming_qubits.push_back(q); + } + } + + assert(victim_positions.size() == incoming_qubits.size()); + std::vector swap_pairs; + swap_pairs.reserve(victim_positions.size()); + for (std::size_t i = 0; i < victim_positions.size(); ++i) { + const auto victim_position = victim_positions[i]; + const auto incoming_position = + layout_.PhysicalPositionOf(incoming_qubits[i]); + swap_pairs.emplace_back(victim_position, incoming_position); + layout_.SwapPositions(victim_position, incoming_position); + } + return swap_pairs; + } + + void LogFixedZonePlacement(const std::vector& usage_scores, + unsigned eviction_floor, + std::size_t num_initial_swaps) const { + if (param_.verbosity <= 1) return; + + IO::messagef("fixed hot zone [%u,%u):", chunk_qubits_, + eviction_floor); + for (unsigned p = chunk_qubits_; p < eviction_floor; ++p) { + const auto q = layout_.LogicalQubitAt(p); + IO::messagef(" q%u(%llu)", q, + static_cast(usage_scores[q])); + } + IO::messagef("; %u initial swaps.\n", unsigned(num_initial_swaps)); + } + + // Places the most frequently used logical qubits outside the in-chunk + // physical zone [chunk_qubits, eviction_floor). In-chunk positions + // remain untouched. Canonical runs restore qubit order at the end. + void PlaceHotQubitsInFixedZone() { + const auto eviction_floor = gate_batch_planner_.EvictionFloor(); + if (partition_.num_blocks == 1 || + eviction_floor <= chunk_qubits_) { + return; + } + + const auto num_fixed_qubits = + eviction_floor - chunk_qubits_; + const auto usage_scores = ComputeQubitUsageScores(); + const auto fixed_qubits = + SelectHighestScoringQubits(usage_scores, chunk_qubits_, + num_fixed_qubits); + auto swap_pairs = + BuildFixedZonePlacementSwaps(fixed_qubits, eviction_floor); + LogFixedZonePlacement(usage_scores, eviction_floor, swap_pairs.size()); + ApplySwapsToState(swap_pairs); + } + + // Rebuilds the plan's pending gates as ordinary gates on PHYSICAL + // qubits, with fresh sequential times, so the standard fuser can + // consume them as if they were a small L-qubit circuit. + template + void BuildBatchOperations(const GateBatchPlan& plan, + std::vector& batch_operations) const { + batch_operations.clear(); + batch_operations.reserve(plan.NumGates()); + + unsigned time = 0; + for (std::size_t idx : plan.gate_indices) { + const PendingGate& gate = pending_gates_[idx]; + + std::vector physical_qubits(gate.Arity()); + for (std::size_t i = 0; i < gate.Arity(); ++i) { + physical_qubits[i] = + layout_.PhysicalPositionOf(gate.logical_qubits[i]); + } + auto physical_matrix = gate.matrix; + NormalizeGateQubitOrder(physical_qubits, physical_matrix); + + // Synthetic gate: kind is irrelevant outside qsimh; sequential + // times satisfy the fuser's ordering contract. + auto physical_gate = Gate{ + 0, time++, std::move(physical_qubits), {}, + std::move(physical_matrix), false}; + batch_operations.push_back(Op{std::move(physical_gate)}); + } + } + + // Deferred until fusion has succeeded, so a failed batch leaves the + // schedule untouched. + void MarkGatesApplied(const GateBatchPlan& plan) { + for (std::size_t idx : plan.gate_indices) { + pending_gates_[idx].MarkApplied(); + } + } + + // fuse(low_gates): runs the standard fuser over the batch's physical + // gates (max_fused_size from param; the proposal suggests 2 or 3) and + // copies the results into ExecutableGate records. max_fused_size == 0 + // is the no-fusion control mode. + template + bool FuseBatchGates(const std::vector& batch_operations) { + auto& executable_gates = + gate_batch_workspace_.executable_gates; + executable_gates.clear(); + + if (param_.max_fused_size == 0) { + executable_gates.reserve(batch_operations.size()); + for (const auto& operation : batch_operations) { + if (!AppendExecutableGate(operation)) return false; + } + return true; + } + + auto fused_ops = + Fuser::FuseGates(param_, partition_.block_qubits, batch_operations); + if (fused_ops.size() == 0 && batch_operations.size() > 0) { + IO::errorf("qsim_gate_batch: fuser failed on a gate batch.\n"); + return false; + } + + executable_gates.reserve(fused_ops.size()); + for (const auto& operation : fused_ops) { + if (!AppendExecutableGate(operation)) return false; + } + return true; + } + + // Copies one (possibly fused) operation into an ExecutableGate, + // normalizing qubit order. The copy must happen while the fuser's input + // container is still alive (fused ops may point into it). + template + bool AppendExecutableGate(const Op& operation) { + const std::vector* physical_qubits = nullptr; + const Matrix* matrix = nullptr; + + if (const auto* fused_gate = + OpGetAlternative>(operation)) { + physical_qubits = &fused_gate->qubits; + matrix = &fused_gate->matrix; + } else if (const auto* raw_gate = + OpGetAlternative>(operation)) { + physical_qubits = &raw_gate->qubits; + matrix = &raw_gate->matrix; + } else { + IO::errorf("qsim_gate_batch: unsupported operation in gate batch.\n"); + return false; + } + + ExecutableGate executable_gate{*physical_qubits, *matrix}; + NormalizeGateQubitOrder(executable_gate.physical_qubits, + executable_gate.matrix); + gate_batch_workspace_.executable_gates.push_back( + std::move(executable_gate)); + return true; + } + + // The proposal's inner loops: for every block i, apply every fused + // gate to the block while it is cache-resident. Blocks in parallel. + // Unit-sized dynamic scheduling balances heterogeneous cores and gate + // workloads consistently, including when the binary is run directly. + static void ExecuteGateBatchOnBlocks( + const std::vector& executable_gates, + fp_type* state_data, + const BlockPartition& partition, SeqSimulator& seq_sim) { + const int64_t num_blocks = partition.num_blocks; + const auto floats_per_block = partition.floats_per_block; + const auto block_qubits = partition.block_qubits; + +#pragma omp parallel for schedule(dynamic, 1) + for (int64_t b = 0; b < num_blocks; ++b) { + fp_type* block_data = state_data + uint64_t(b) * floats_per_block; + auto block_view = SeqStateSpace::Create(block_data, block_qubits); + + for (const ExecutableGate& gate : executable_gates) { + seq_sim.ApplyGate(gate.physical_qubits, gate.matrix.data(), + block_view); + } + } + } + + // ======== Final cleanup ======== + + void RestoreIdentityQubitOrder() { + std::vector swap_pairs; + + while (layout_.BuildIdentityRestorationSwaps(swap_pairs)) { + ApplySwapsToState(swap_pairs); + ++simulation_stats_.num_restore_passes; + } + } + + // ======== Reporting ======== + + void LogSimulationSummary(double simulation_start) const { + if (param_.verbosity == 0) return; + + IO::messagef("simu time is %g seconds.\n", + GetTime() - simulation_start); + IO::messagef("gate-batch runner: %u gate batches, %u restore passes, " + "%u raw gates, %u executable gates, " + "%u qubit swaps.\n", + simulation_stats_.num_gate_batches, + simulation_stats_.num_restore_passes, + unsigned(pending_gates_.size()), + simulation_stats_.num_executable_gates, + simulation_stats_.num_swaps); + if (param_.verbosity > 1) { + IO::messagef("breakdown: plan %g s, swap passes %g s, fuse %g s, " + "gate passes %g s.\n", + simulation_stats_.plan_seconds, + simulation_stats_.swap_seconds, + simulation_stats_.fuse_seconds, + simulation_stats_.gate_seconds); + } + } + + // State owned or borrowed by one invocation of the public static Run. + // The instance never escapes Run, so the borrowed parameter and state + // storage remain valid for its complete lifetime. + const Parameter& param_; + const BlockPartition partition_; + fp_type* const state_data_; + const unsigned chunk_qubits_; // Low amplitude bits inside a state chunk. + SeqSimulator seq_sim_; + std::vector pending_gates_; + QubitLayout& layout_; + GateBatchPlanner gate_batch_planner_; + GateBatchWorkspace gate_batch_workspace_; + SimulationStats simulation_stats_; +}; + +} // namespace qsim + +#endif // RUN_QSIM_GATE_BATCH_H_ diff --git a/lib/statespace_avx.h b/lib/statespace_avx.h index f68733f0f..ae221ceff 100644 --- a/lib/statespace_avx.h +++ b/lib/statespace_avx.h @@ -74,6 +74,9 @@ class StateSpaceAVX : using State = typename Base::State; using fp_type = typename Base::fp_type; + // Number of low amplitude-index bits stored inside one SIMD chunk. + static constexpr unsigned kChunkQubits = 3; + template explicit StateSpaceAVX(ForArgs&&... args) : Base(args...) {} diff --git a/lib/statespace_avx512.h b/lib/statespace_avx512.h index e7f431183..35638ddb2 100644 --- a/lib/statespace_avx512.h +++ b/lib/statespace_avx512.h @@ -86,6 +86,9 @@ class StateSpaceAVX512 : using State = typename Base::State; using fp_type = typename Base::fp_type; + // Number of low amplitude-index bits stored inside one SIMD chunk. + static constexpr unsigned kChunkQubits = 4; + template explicit StateSpaceAVX512(ForArgs&&... args) : Base(args...) {} diff --git a/lib/statespace_basic.h b/lib/statespace_basic.h index 41a7f8bf1..d06334aab 100644 --- a/lib/statespace_basic.h +++ b/lib/statespace_basic.h @@ -41,6 +41,9 @@ class StateSpaceBasic : using State = typename Base::State; using fp_type = typename Base::fp_type; + // Number of low amplitude-index bits stored inside one state chunk. + static constexpr unsigned kChunkQubits = 0; + template explicit StateSpaceBasic(ForArgs&&... args) : Base(args...) {} diff --git a/lib/statespace_neon.h b/lib/statespace_neon.h index 14480cd0b..f6501ffd7 100644 --- a/lib/statespace_neon.h +++ b/lib/statespace_neon.h @@ -66,6 +66,9 @@ class StateSpaceNEON : using State = typename Base::State; using fp_type = typename Base::fp_type; + // Number of low amplitude-index bits stored inside one SIMD chunk. + static constexpr unsigned kChunkQubits = 2; + template explicit StateSpaceNEON(ForArgs&&... args) : Base(args...) {} diff --git a/lib/statespace_sse.h b/lib/statespace_sse.h index a3ed52c95..6ad0978e1 100644 --- a/lib/statespace_sse.h +++ b/lib/statespace_sse.h @@ -71,6 +71,9 @@ class StateSpaceSSE : using State = typename Base::State; using fp_type = typename Base::fp_type; + // Number of low amplitude-index bits stored inside one SIMD chunk. + static constexpr unsigned kChunkQubits = 2; + template explicit StateSpaceSSE(ForArgs&&... args) : Base(args...) {}