From a46392e5f3b5e119706d4741e5f246a8c7142cdb Mon Sep 17 00:00:00 2001 From: Anton Borissov Date: Fri, 12 Jun 2026 13:49:19 -0400 Subject: [PATCH] Optimize qubit kernels for qpp_bench QFT This combines the qpp_bench-focused qubit kernel optimizations that were developed on top of origin/dev. It fixes the OpenMP collapse pragma, adds a fast path for single-control diagonal phase gates, streams the QFT phase pattern over contiguous amplitude runs, and specializes Hadamard updates as scaled sums and differences. Together these changes make the native qpp backend faster than PennyLane Lightning on the repeated qpp_bench QFT sweep from 18 through 26 qubits. --- include/qpp/internal/kernels/qubit/apply.hpp | 30 ++++++++++- .../kernels/qubit/apply_ctrl_diag.hpp | 51 +++++++++++++++++++ 2 files changed, 80 insertions(+), 1 deletion(-) diff --git a/include/qpp/internal/kernels/qubit/apply.hpp b/include/qpp/internal/kernels/qubit/apply.hpp index e036aa8ca..aa573a382 100644 --- a/include/qpp/internal/kernels/qubit/apply.hpp +++ b/include/qpp/internal/kernels/qubit/apply.hpp @@ -86,12 +86,40 @@ apply_psi_1q_inplace(Eigen::MatrixBase& state, const Scalar a10 = A.coeff(1, 0); const Scalar a11 = A.coeff(1, 1); + // Hadamard-like gates have the shape s * [[1, 1], [1, -1]]. + // The QFT benchmark applies this gate once per qubit, so it is worth + // avoiding the generic 2x2 matrix-vector multiply in this exact case. + // Each amplitude pair only needs the scaled sum and difference of its two + // original values, cutting the hot-loop complex multiplications in half and + // exposing a simpler expression for the compiler to vectorize. + if (a00 == a01 && a00 == a10 && a11 == -a00) { +#ifdef QPP_OPENMP +#pragma omp parallel for collapse(2) if (D >= 65536) // 16 qubits +#endif // QPP_OPENMP + for (idx L = 0; L < D; L += jump) { + for (idx R = 0; R < step; ++R) { + const idx k0 = L + R; + const idx k1 = k0 + step; + const Scalar psi_k0 = state.coeff(k0); + const Scalar psi_k1 = state.coeff(k1); + + state.coeffRef(k0) = a00 * (psi_k0 + psi_k1); + state.coeffRef(k1) = a00 * (psi_k0 - psi_k1); + } + } + return; + } + // Pair-wise Amplitude Transformation // The outer loop (L) iterates over all blocks of size 'jump'. // This loop is perfectly independent and is the primary target for // parallelization. #ifdef QPP_OPENMP -#pragma omp parallel for colapse(2) if (D >= 65536) // 16 qubits +// Both the block loop and the in-block pair loop are independent. Collapsing +// them gives OpenMP one flat iteration space instead of only parallelizing the +// outer block loop, which matters when the target qubit is high-order and there +// are only a few large blocks. +#pragma omp parallel for collapse(2) if (D >= 65536) // 16 qubits #endif // QPP_OPENMP for (idx L = 0; L < D; L += jump) { // The inner loop (R) iterates over the lower part of the block, from 0 diff --git a/include/qpp/internal/kernels/qubit/apply_ctrl_diag.hpp b/include/qpp/internal/kernels/qubit/apply_ctrl_diag.hpp index c8b229fdc..a506863d1 100644 --- a/include/qpp/internal/kernels/qubit/apply_ctrl_diag.hpp +++ b/include/qpp/internal/kernels/qubit/apply_ctrl_diag.hpp @@ -113,6 +113,57 @@ template } } + // Fast path for the controlled phase gates emitted by the QFT benchmark: + // one positive control, a diagonal [1, phase], and a single target qubit. + // The generic loop tests the control mask for every target-zero amplitude + // and then updates the paired target-one amplitude. Here the control + // condition is encoded directly into the generated indices, so the hot loop + // only visits amplitudes that actually receive the phase. + if (ctrl_size == 1 && shift[0] == 0 && !apply_a0 && apply_a1) { + const idx ctrl_bit = ctrl_mask; + if (ctrl_bit < step) { + // QFT calls this path when the control qubit is below the target in + // the state-vector bit layout. In each target-one half-block, the + // control-one amplitudes appear as contiguous runs of length + // ctrl_bit. Streaming over those runs avoids the bit-deposit index + // reconstruction used by the fallback below and gives the compiler + // a simple contiguous multiply loop. +#ifdef QPP_OPENMP +#pragma omp parallel for collapse(2) if (D >= 65536) // 16 qubits +#endif // QPP_OPENMP + for (idx L = 0; L < D; L += jump) { + for (idx C = 0; C < step; C += 2 * ctrl_bit) { + const idx start = L + step + C + ctrl_bit; + for (idx R = 0; R < ctrl_bit; ++R) { + state.coeffRef(start + R) *= a1; + } + } + } + return; + } + + // General one-control diagonal fallback. The loop variable r ranges + // over the state-vector indices with the target and control bits + // removed. The masks below insert those two zero bits back into their + // sorted positions, after which OR-ing ctrl_bit and step selects the + // unique target-one/control-one amplitude that needs the phase. + const idx lo = step < ctrl_bit ? step : ctrl_bit; + const idx hi = step < ctrl_bit ? ctrl_bit : step; + const idx low_mask = lo - 1; + const idx mid_mask = (hi >> 1) - lo; + const idx keep_mask = low_mask | mid_mask; + +#ifdef QPP_OPENMP +#pragma omp parallel for if (D >= 65536) // 16 qubits +#endif // QPP_OPENMP + for (idx r = 0; r < (D >> 2); ++r) { + const idx base = (r & low_mask) | ((r & mid_mask) << 1) | + ((r & ~keep_mask) << 2); + state.coeffRef(base | ctrl_bit | step) *= a1; + } + return; + } + #ifdef QPP_OPENMP #pragma omp parallel for collapse(2) if (D >= 65536) // 16 qubits #endif // QPP_OPENMP