From 69f40d11aa553fd8b576463be250945eb17e6140 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Thu, 1 Oct 2026 14:18:25 -0700 Subject: [PATCH 01/27] Extend auto-multithreading --- quest/src/cpu/cpu_subroutines.cpp | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/quest/src/cpu/cpu_subroutines.cpp b/quest/src/cpu/cpu_subroutines.cpp index c1ba9b3aa..a9b7c61de 100644 --- a/quest/src/cpu/cpu_subroutines.cpp +++ b/quest/src/cpu/cpu_subroutines.cpp @@ -829,7 +829,7 @@ void cpu_statevec_allTargDiagMatr_sub(Qureg qureg, FullStateDiagMatr matr, qcomp // every iteration modifies one amp, using one element qindex numIts = qureg.numAmpsPerNode; - #pragma omp parallel for if(qureg.isMultithreaded||qureg.isMultithreaded) + #pragma omp parallel for if(qureg.isMultithreaded||matr.isMultithreaded) for (qindex n=0; n coeffs // use template param to compile-time unroll inner loop below int numQuregs = accel_tryUseCompileTimeValue(inQuregs.size()); + // below, we choose to multithread only if outQureg is multithreaded, + // regardless of whether any of the inQuregs are multithreaded + #pragma omp parallel for if(outQureg.isMultithreaded) for (qindex n=0; n Date: Thu, 1 Oct 2026 14:23:27 -0700 Subject: [PATCH 02/27] Patch distributed error propagation We require a validation failure on one node to reach all other nodes; this requires an OR across nodes. The code previously erroneously performed an AND, requiring the failure occur on all nodes in order to be propagated. --- quest/src/api/channels.cpp | 2 +- quest/src/api/matrices.cpp | 2 +- quest/src/api/paulis.cpp | 2 +- quest/src/api/qureg.cpp | 2 +- quest/src/comm/comm_routines.cpp | 16 ++++++++++++++++ quest/src/comm/comm_routines.hpp | 2 ++ 6 files changed, 22 insertions(+), 4 deletions(-) diff --git a/quest/src/api/channels.cpp b/quest/src/api/channels.cpp index c6702438a..72aa22772 100644 --- a/quest/src/api/channels.cpp +++ b/quest/src/api/channels.cpp @@ -108,7 +108,7 @@ void freeAllMemoryIfAnyAllocsFailed(T& obj) { // determine whether any node experienced a failure bool anyFail = didAnyLocalAllocsFail(obj); if (comm_isActive()) - anyFail = comm_isTrueOnAllNodes(anyFail); + anyFail = comm_isTrueOnAnyNode(anyFail); // if so, free all memory before subsequent validation if (anyFail) diff --git a/quest/src/api/matrices.cpp b/quest/src/api/matrices.cpp index 07e37025c..a741cea40 100644 --- a/quest/src/api/matrices.cpp +++ b/quest/src/api/matrices.cpp @@ -166,7 +166,7 @@ void freeAllMemoryIfAnyAllocsFailed(T matr) { // ascertain whether any allocs failed on any node bool anyFail = didAnyLocalAllocsFail(matr); if (comm_isActive()) - anyFail = comm_isTrueOnAllNodes(anyFail); + anyFail = comm_isTrueOnAnyNode(anyFail); // if so, free all heap fields if (anyFail) diff --git a/quest/src/api/paulis.cpp b/quest/src/api/paulis.cpp index d5ac4d8e7..4875bb542 100644 --- a/quest/src/api/paulis.cpp +++ b/quest/src/api/paulis.cpp @@ -39,7 +39,7 @@ bool didAnyAllocsFailOnAnyNode(PauliStrSum sum) { ! mem_isAllocated(sum.isApproxHermitian) ); if (comm_isActive()) - anyFail = comm_isTrueOnAllNodes(anyFail); + anyFail = comm_isTrueOnAnyNode(anyFail); return anyFail; } diff --git a/quest/src/api/qureg.cpp b/quest/src/api/qureg.cpp index 683ca183d..677827b46 100644 --- a/quest/src/api/qureg.cpp +++ b/quest/src/api/qureg.cpp @@ -117,7 +117,7 @@ bool didAnyAllocsFailOnAnyNode(Qureg qureg) { bool anyFail = didAnyLocalAllocsFail(qureg); if (comm_isActive()) - anyFail = comm_isTrueOnAllNodes(anyFail); + anyFail = comm_isTrueOnAnyNode(anyFail); return anyFail; } diff --git a/quest/src/comm/comm_routines.cpp b/quest/src/comm/comm_routines.cpp index cf6956454..852316625 100644 --- a/quest/src/comm/comm_routines.cpp +++ b/quest/src/comm/comm_routines.cpp @@ -777,6 +777,22 @@ bool comm_isTrueOnAllNodes(bool val) { } +bool comm_isTrueOnAnyNode(bool val) { +#if QUEST_COMPILE_MPI + + // perform global OR and broadcast result back to all nodes + int local = (int) val; + int global; + MPI_Allreduce(&local, &global, 1, MPI_INT, MPI_LOR, comm_getMpiComm()); + return (bool) global; + +#else + error_commButEnvNotDistributed(); + return false; +#endif +} + + bool comm_isTrueOnRootNode(bool val) { #if QUEST_COMPILE_MPI diff --git a/quest/src/comm/comm_routines.hpp b/quest/src/comm/comm_routines.hpp index e75e889f6..d2bceae5c 100644 --- a/quest/src/comm/comm_routines.hpp +++ b/quest/src/comm/comm_routines.hpp @@ -69,6 +69,8 @@ void comm_reduceReals(qreal* localReals, qindex numLocalReals); bool comm_isTrueOnAllNodes(bool val); +bool comm_isTrueOnAnyNode(bool val); + bool comm_isTrueOnRootNode(bool val); From 9850367d53c846b2b0496dc7d2ae7d4c865970fe Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Thu, 1 Oct 2026 14:25:48 -0700 Subject: [PATCH 03/27] Add missing validation --- quest/src/api/calculations.cpp | 5 +++++ quest/src/api/decoherence.cpp | 1 + quest/src/api/qureg.cpp | 12 +++++++++++- 3 files changed, 17 insertions(+), 1 deletion(-) diff --git a/quest/src/api/calculations.cpp b/quest/src/api/calculations.cpp index 47e5d8a63..e18227cd5 100644 --- a/quest/src/api/calculations.cpp +++ b/quest/src/api/calculations.cpp @@ -428,6 +428,11 @@ qreal calcProbOfMultiQubitOutcome(Qureg qureg, vector qubits, vector o vector calcProbsOfAllMultiQubitOutcomes(Qureg qureg, vector qubits) { + // validate before exponentiating the caller-controlled vector size + validate_quregFields(qureg, __func__); + validate_targets(qureg, qubits.data(), qubits.size(), __func__); + validate_measurementOutcomesFitInGpuMem(qureg, qubits.size(), __func__); + // allocate temp vector, and validate successful (since it's exponentially large!) vector out; qindex numOut = powerOf2(qubits.size()); diff --git a/quest/src/api/decoherence.cpp b/quest/src/api/decoherence.cpp index 4e1901f25..3bd393e20 100644 --- a/quest/src/api/decoherence.cpp +++ b/quest/src/api/decoherence.cpp @@ -143,6 +143,7 @@ void mixQureg(Qureg outQureg, Qureg inQureg, qreal inProb) { void mixSuperOp(Qureg qureg, int* targets, int numTargets, SuperOp superop) { validate_quregFields(qureg, __func__); + validate_quregIsDensityMatrix(qureg, __func__); validate_targets(qureg, targets, numTargets, __func__); validate_superOpFields(superop, __func__); validate_superOpIsSynced(superop, __func__); diff --git a/quest/src/api/qureg.cpp b/quest/src/api/qureg.cpp index 677827b46..de726d7fa 100644 --- a/quest/src/api/qureg.cpp +++ b/quest/src/api/qureg.cpp @@ -534,6 +534,11 @@ extern "C" void _wrap_getDensityQuregAmp(qcomp* out, Qureg qureg, qindex row, qi vector getQuregAmps(Qureg qureg, qindex startInd, qindex numAmps) { + // validate before interpreting numAmps as an allocation size + validate_quregFields(qureg, __func__); + validate_quregIsStateVector(qureg, __func__); + validate_basisStateIndices(qureg, startInd, numAmps, __func__); + // allocate the output vector, and validate successful vector out; auto callback = [&]() { validate_tempListAllocSucceeded(false, numAmps, sizeof(qcomp), __func__); }; @@ -547,9 +552,14 @@ vector getQuregAmps(Qureg qureg, qindex startInd, qindex numAmps) { vector> getDensityQuregAmps(Qureg qureg, qindex startRow, qindex startCol, qindex numRows, qindex numCols) { + // validate before multiplying dimensions or interpreting them as allocation sizes + validate_quregFields(qureg, __func__); + validate_quregIsDensityMatrix(qureg, __func__); + validate_basisStateRowCols(qureg, startRow, startCol, numRows, numCols, __func__); + // allocate the output matrix, and validate successful vector> out; - qindex numElems = numRows * numCols; // never overflows (else Qureg alloc would fail) + qindex numElems = numRows * numCols; auto callback1 = [&]() { validate_tempListAllocSucceeded(false, numElems, sizeof(qcomp), __func__); }; util_tryAllocMatrix(out, numRows, numCols, callback1); From 66c2878ef9cbb75acc99c1a030c70d52f8d1a5c3 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Thu, 1 Oct 2026 14:29:42 -0700 Subject: [PATCH 04/27] Correct validation throwing function names --- quest/src/api/matrices.cpp | 4 ++-- quest/src/api/qureg.cpp | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/quest/src/api/matrices.cpp b/quest/src/api/matrices.cpp index a741cea40..10a889341 100644 --- a/quest/src/api/matrices.cpp +++ b/quest/src/api/matrices.cpp @@ -310,7 +310,7 @@ FullStateDiagMatr validateAndCreateCustomFullStateDiagMatr(int numQubits, int us out.cpuElems = cpu_allocArray(numElemsPerNode); out.gpuElems = (useGpuAccel)? gpu_allocArray(numElemsPerNode) : nullptr; - validateMatrixAllocs(out, __func__); + validateMatrixAllocs(out, caller); setInitialHeapFlags(out); return out; } @@ -743,7 +743,7 @@ extern "C" void setFullStateDiagMatrFromMultiDimLists(FullStateDiagMatr out, voi template void validateAndPrintMatrix(T matr, const char* caller) { validate_matrixFields(matr, caller); - validate_numReportedNewlinesAboveZero(__func__); // because trailing newline mandatory + validate_numReportedNewlinesAboveZero(caller); // because trailing newline mandatory // syncable matrices must be synced before reporting (though only CPU elems are printed) if constexpr (util_isHeapMatrixType()) diff --git a/quest/src/api/qureg.cpp b/quest/src/api/qureg.cpp index de726d7fa..a0a368b1a 100644 --- a/quest/src/api/qureg.cpp +++ b/quest/src/api/qureg.cpp @@ -166,7 +166,7 @@ Qureg validateAndCreateCustomQureg(int numQubits, int isDensMatr, int useDistrib // if any of the above mallocs failed, below validation will memory leak; so free first (but don't set to nullptr) freeAllMemoryIfAnyAllocsFailed(qureg); - validate_newQuregAllocs(qureg, __func__); + validate_newQuregAllocs(qureg, caller); // TODO / DEBUG: // TEMPORARILY FORBID HIP GPU ACCELERATION OF > 2^32 AMPS PER NODE, From 51c42bafcb90c079954288979594b5ccb682f8b8 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Thu, 1 Oct 2026 14:33:16 -0700 Subject: [PATCH 05/27] Include max nqubits in auto-deploy report Pretty much pointless since there's no new deployment that switches on only at the absolute maximum Qureg size before overflow, but eh, an aesthetic bug --- quest/src/api/environment.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/quest/src/api/environment.cpp b/quest/src/api/environment.cpp index e6fd2903a..37c523646 100644 --- a/quest/src/api/environment.cpp +++ b/quest/src/api/environment.cpp @@ -375,7 +375,7 @@ void printQuregAutoDeployments(bool isDensMatr) { // auto-deploy will still try to deploy there to (then subsequent validation will fail) int maxQubits = mem_getMaxNumQuregQubitsBeforeGlobalMemSizeofOverflow(isDensMatr, global_envPtr->numNodes); - for (int numQubits=1; numQubits Date: Thu, 1 Oct 2026 14:37:16 -0700 Subject: [PATCH 06/27] Patch C++ setDensityQuregAmps which erroneously passed (numRows,numRows) instead of (numRows,numCols), to the C API. This meant passing a rectangular (non-square) set of amplitudes either caused a seg-fault, or discarding some of the passed amplitudes. --- quest/src/api/initialisations.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/quest/src/api/initialisations.cpp b/quest/src/api/initialisations.cpp index 36f63910c..a8be7ae3e 100644 --- a/quest/src/api/initialisations.cpp +++ b/quest/src/api/initialisations.cpp @@ -293,7 +293,8 @@ void setDensityQuregAmps(Qureg qureg, qindex startRow, qindex startCol, vector0)? len : 0); // avoid seg-fault + size_t numCols = (len > 0)? amps[0].size() : 0; // avoid seg-fault + setDensityQuregAmps(qureg, startRow, startCol, ptrs.data(), len, numCols); } void setQuregAmps(Qureg qureg, qindex startInd, vector amps) { From baf63b2425d20f795969766a552fe741d8d76cb3 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Thu, 1 Oct 2026 14:39:14 -0700 Subject: [PATCH 07/27] Patch NUMA-aware destroyQureg destroyQureg previously erroneously used the non-NUMA dealloc --- quest/src/api/qureg.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/quest/src/api/qureg.cpp b/quest/src/api/qureg.cpp index a0a368b1a..57de85c0b 100644 --- a/quest/src/api/qureg.cpp +++ b/quest/src/api/qureg.cpp @@ -130,7 +130,7 @@ void freeAllMemoryIfAnyAllocsFailed(Qureg qureg) { return; // otherwise, free everything that was successfully allocated (freeing nullptr is legal) - cpu_deallocArray(qureg.cpuAmps); + cpu_deallocNumaArray(qureg.cpuAmps, qureg.numAmpsPerNode); cpu_deallocArray(qureg.cpuCommBuffer); // although we avoid calling GPU deallocation in non-GPU mode From 3b35c119ec8d0503ba0690e6e2a37f7ca7661c45 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Thu, 1 Oct 2026 14:50:04 -0700 Subject: [PATCH 08/27] Avoid overflows with safer types --- quest/src/api/trotterisation.cpp | 2 +- quest/src/core/bitwise.hpp | 4 ++-- quest/src/core/localiser.cpp | 2 +- quest/src/core/paulilogic.cpp | 6 +++--- quest/src/cpu/cpu_config.cpp | 2 +- quest/src/cpu/cpu_subroutines.cpp | 4 ++-- quest/src/gpu/gpu_kernels.cuh | 6 +++--- 7 files changed, 13 insertions(+), 13 deletions(-) diff --git a/quest/src/api/trotterisation.cpp b/quest/src/api/trotterisation.cpp index 6fd5781ba..56c09792b 100644 --- a/quest/src/api/trotterisation.cpp +++ b/quest/src/api/trotterisation.cpp @@ -36,7 +36,7 @@ void internal_applyFirstOrderTrotterRepetition( ) { // apply each sum term as a gadget, in forward or reverse order for (qindex i=0; i> groups; // group sum's terms into those with identical communication patterns - for (int i=0; i> src, qindex dim) { void cpu_copyPauliStrSum(PauliStrSum out, PauliStr* strings, qcomp* coeffs) { // serially copy data over to new heap memory - for (int i=0; i(blockIdx.x) * blockDim.x + threadIdx.x; } @@ -71,8 +71,8 @@ __host__ qindex getNumBlocks(qindex numThreads, int numThreadsPerBlock) { /// improve this with cudaOccupancyMaxPotentialBlockSize(), /// making it function specific - // CUDA ceil - return ceil(numThreads / static_cast(numThreadsPerBlock)); + // exact integer ceiling, without losing precision for large launches + return numThreads / numThreadsPerBlock + (numThreads % numThreadsPerBlock != 0); } From 0cd68f4c016e0af683c7c14ebd1f37fc8a089cf7 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Thu, 1 Oct 2026 14:50:49 -0700 Subject: [PATCH 09/27] Correct internal error typo --- quest/src/core/errors.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/quest/src/core/errors.cpp b/quest/src/core/errors.cpp index c78f1c727..5f98822c2 100644 --- a/quest/src/core/errors.cpp +++ b/quest/src/core/errors.cpp @@ -419,7 +419,7 @@ void assert_exponentMatchesTemplateParam(qcomp exponent, bool hasPower, bool use assert_exponentMatchesTemplateParam(exponent, hasPower); if (!hasPower && useRealPow) - raiseInternalError("A CPU or GPU subroutine received an invalid combination of HasPower=false and UseRealPow=false template parameters"); + raiseInternalError("A CPU or GPU subroutine received an invalid combination of HasPower=false and UseRealPow=true template parameters"); } From 878803439e676438bdf97a59647b76829f98c031 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Sun, 4 Oct 2026 11:39:59 -0700 Subject: [PATCH 10/27] Patch multi-GPU oversubscription check --- quest/src/comm/comm_routines.cpp | 18 +++++++++--------- quest/src/comm/comm_routines.hpp | 2 +- quest/src/gpu/gpu_config.cpp | 23 ++++++++--------------- 3 files changed, 18 insertions(+), 25 deletions(-) diff --git a/quest/src/comm/comm_routines.cpp b/quest/src/comm/comm_routines.cpp index 852316625..62bef2dfa 100644 --- a/quest/src/comm/comm_routines.cpp +++ b/quest/src/comm/comm_routines.cpp @@ -816,25 +816,25 @@ bool comm_isTrueOnRootNode(bool val) { */ -vector comm_gatherStringsToRoot(char* localChars, int maxNumLocalChars) { +vector comm_gatherStringsToRoot(const char* localChars, int numLocalChars) { #if QUEST_COMPILE_MPI // no need to validate array sizes and memory alloc successes; - // these are trivial O(#nodes)-size arrays containing <20 chars + // these are trivial O(#nodes)-size arrays containing few bytes int numNodes = comm_getNumNodes(); - // root makes one big contiguous array to receive all chars contiguously - vector allChars(maxNumLocalChars * numNodes); + // root makes one big contiguous array to receive all bytes contiguously + vector allChars(numLocalChars * numNodes); - // all nodes send root all their local chars + // all nodes send root their fixed-width byte sequence int recvRank = ROOT_RANK; - MPI_Gather(localChars, maxNumLocalChars, MPI_CHAR, allChars.data(), - maxNumLocalChars, MPI_CHAR, recvRank, comm_getMpiComm()); + MPI_Gather(localChars, numLocalChars, MPI_CHAR, allChars.data(), + numLocalChars, MPI_CHAR, recvRank, comm_getMpiComm()); - // divide allChars into stings, delimited by each node's terminal char + // preserve every byte, including embedded nulls, in fixed-width strings vector out(numNodes); for (int r=0; r comm_gatherStringsToRoot(char* localChars, int maxNumLocalChars); +vector comm_gatherStringsToRoot(const char* localChars, int numLocalChars); diff --git a/quest/src/gpu/gpu_config.cpp b/quest/src/gpu/gpu_config.cpp index 001cc62c0..5feff91a5 100644 --- a/quest/src/gpu/gpu_config.cpp +++ b/quest/src/gpu/gpu_config.cpp @@ -370,14 +370,11 @@ int gpu_getMaxNumThreadsPerBlock() { } -std::array getBoundGpuUuid() { +std::array getBoundGpuUuid() { #if QUEST_COMPILE_CUDA assert_gpuHasBeenBound(hasGpuBeenBound); - constexpr int numUuidChars = 16; - constexpr int numOutChars = numUuidChars + 1; // terminal char - - std::array out; + constexpr int numUuidBytes = 16; cudaUUID_t uuid; // ROCm v5's cudaDeviceProp doesn't have a uuid field @@ -391,12 +388,8 @@ std::array getBoundGpuUuid() { uuid = prop.uuid; #endif - // copy char[16] to out[17] (not human readable) - std::copy_n(uuid.bytes, numUuidChars, out.begin()); - - // include terminal char so that subsequent to-string - // operations will succeed without knowing string length - out[numOutChars-1] = '\0'; + std::array out; + std::copy_n(uuid.bytes, numUuidBytes, out.begin()); return out; #else @@ -437,11 +430,11 @@ bool gpu_areAnyNodesBoundToSameGpu() { if (!comm_isActive()) return false; - // obtain bound GPU's UUID; a unique identifier 16-char identifier - auto uuidStr = getBoundGpuUuid(); + // obtain the bound GPU's fixed-width binary UUID + auto uuid = getBoundGpuUuid(); - // we can repurpose string-to-root sending to collect all uuids - auto allUuids = comm_gatherStringsToRoot(uuidStr.data(), uuidStr.size()); + // gather every byte, including embedded nulls, and compare the identifiers + auto allUuids = comm_gatherStringsToRoot(uuid.data(), static_cast(uuid.size())); auto uniqueUuids = std::set(allUuids.begin(), allUuids.end()); // and assess whether they're all unique (non-root's bools are overwritten) From bd5692870fd373c325ce7b3e222e5263e077b90d Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Sun, 4 Oct 2026 11:43:53 -0700 Subject: [PATCH 11/27] Patch Hermiticity validation which previously did not validate matrix diagonals were real, and ergo returned true for non-Hermitian matrices! --- quest/src/core/utilities.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/quest/src/core/utilities.cpp b/quest/src/core/utilities.cpp index 71facdf64..9c3ab46bf 100644 --- a/quest/src/core/utilities.cpp +++ b/quest/src/core/utilities.cpp @@ -635,9 +635,9 @@ bool getHermiticity(T elems, qindex dim, qreal eps) { /// consider multithreading or GPU-accelerating this /// when caller is big and e.g. has GPU memory - // check adjoint(elems) == elems + // check adjoint(elems) == elems, including that each diagonal is real for (qindex r=0; r Date: Sun, 4 Oct 2026 12:06:29 -0700 Subject: [PATCH 12/27] Patch CPTP validation which previously numerically differed from the doc --- quest/src/core/utilities.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/quest/src/core/utilities.cpp b/quest/src/core/utilities.cpp index 9c3ab46bf..a5f1b3f46 100644 --- a/quest/src/core/utilities.cpp +++ b/quest/src/core/utilities.cpp @@ -857,6 +857,7 @@ bool util_isCPTP(KrausMap map, qreal eps) { /// use a multithreaded routine *(map.isApproxCPTP) = 1; + qreal epsSquared = eps * eps; // check whether each element satisfies Identity = sum dagger(m)*m for (qindex r=0; r eps) { + if (distSquared > epsSquared) { // by recording the result and returning immediately *(map.isApproxCPTP) = 0; From 6798e3a6808d751b0724f646b2f7e6871420597b Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Sun, 4 Oct 2026 12:36:31 -0700 Subject: [PATCH 13/27] Forbid zero-prob zero sample --- quest/src/core/randomiser.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/quest/src/core/randomiser.cpp b/quest/src/core/randomiser.cpp index 7b35a29fc..c3c47751f 100644 --- a/quest/src/core/randomiser.cpp +++ b/quest/src/core/randomiser.cpp @@ -135,7 +135,7 @@ int rand_getRandomSingleQubitOutcome(qreal probOfZero) { qreal sample = distrib(mainGenerator); // produces 0 with probOfZero probability - return sample > probOfZero; + return sample >= probOfZero; } From 701af9207a83a5779f8dc96107283799f3c573e9 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Sun, 4 Oct 2026 12:41:29 -0700 Subject: [PATCH 14/27] Do not oversubscribe threads in AnyTargDenseMatr --- quest/src/gpu/gpu_subroutines.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/quest/src/gpu/gpu_subroutines.cpp b/quest/src/gpu/gpu_subroutines.cpp index e9936c9c3..ac54957fd 100644 --- a/quest/src/gpu/gpu_subroutines.cpp +++ b/quest/src/gpu/gpu_subroutines.cpp @@ -498,9 +498,10 @@ void gpu_statevec_anyCtrlAnyTargDenseMatr_sub(Qureg qureg, ConstList64 ctrls, Co qindex numThreads = gpu_getMaxNumConcurrentThreads(); int numThreadsPerBlock = gpu_getNumThreadsPerBlock(); - // use strictly 2^# threads to maintain precondition of all kernels + // use at most a power-of-two number of concurrent threads to maintain + // the kernels' divisibility precondition without exceeding the hardware bound if (!isPowerOf2(numThreads)) - numThreads = util_getNextPowerOf2(numThreads); + numThreads = util_getNextPowerOf2(numThreads) / 2; // no point in dispatching more threads than batches if (numThreads > numBatches) From 17a5b10e19c42830006a9cb5e1dc69dfbcdd7eab Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Sun, 4 Oct 2026 12:44:40 -0700 Subject: [PATCH 15/27] Patch seed validation error msg --- quest/src/core/validation.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/quest/src/core/validation.cpp b/quest/src/core/validation.cpp index beeae12f1..0fd108283 100644 --- a/quest/src/core/validation.cpp +++ b/quest/src/core/validation.cpp @@ -1642,7 +1642,7 @@ void validate_randomSeeds(unsigned* seeds, int numSeeds, const char* caller) { } assertThat(!isNull, report::RANDOM_SEEDS_PTR_IS_NULL, caller); - assertThat(numRootSeeds > 0, report::INVALID_NUM_RANDOM_SEEDS, {{"${NUM_SEEDS}", numSeeds}}, caller); + assertThat(numRootSeeds > 0, report::INVALID_NUM_RANDOM_SEEDS, {{"${NUM_SEEDS}", numRootSeeds}}, caller); } void validate_newEpsilonValue(qreal eps, const char* caller) { From 4692841902b52f79eeb298480842f033e50890cb Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Sun, 4 Oct 2026 12:46:25 -0700 Subject: [PATCH 16/27] Patch setDensityQuregAmps dim validation --- quest/src/core/validation.cpp | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/quest/src/core/validation.cpp b/quest/src/core/validation.cpp index 0fd108283..a3cdff851 100644 --- a/quest/src/core/validation.cpp +++ b/quest/src/core/validation.cpp @@ -2585,21 +2585,19 @@ void validate_matrixRowsAllSameSize(vector> matrix, const char* ca size_t dim = matrix[0].size(); - size_t row=0; - bool allSame = true; - for (row=0; row Date: Sun, 4 Oct 2026 12:58:13 -0700 Subject: [PATCH 17/27] Extremely minor boundary changes --- quest/src/core/parser.cpp | 2 +- quest/src/core/utilities.cpp | 7 ++++++- quest/src/core/validation.cpp | 2 +- 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/quest/src/core/parser.cpp b/quest/src/core/parser.cpp index 9d9194a3f..86ca4f611 100644 --- a/quest/src/core/parser.cpp +++ b/quest/src/core/parser.cpp @@ -118,7 +118,7 @@ namespace regexes { bool isWhiteSpace(char ch) { // matches patterns:space and more (additionally; newlines) - return isspace(ch); + return std::isspace(static_cast(ch)); } bool isNotWhiteSpace(char ch) { diff --git a/quest/src/core/utilities.cpp b/quest/src/core/utilities.cpp index a5f1b3f46..14c4bb1b7 100644 --- a/quest/src/core/utilities.cpp +++ b/quest/src/core/utilities.cpp @@ -398,6 +398,11 @@ bool util_willProdOverflow(vector terms) { qindex max = std::numeric_limits::max(); qindex prod = 1; + // any zero factor makes the full product zero and cannot overflow + for (auto x : terms) + if (x == 0) + return false; + for (auto x : terms) { // division floors, so prod strictly exceed @@ -417,7 +422,7 @@ bool util_willSumOverflow(vector terms) { for (auto x : terms) { - if (sum >= max - x) + if (sum > max - x) return true; sum += x; diff --git a/quest/src/core/validation.cpp b/quest/src/core/validation.cpp index a3cdff851..2581a9892 100644 --- a/quest/src/core/validation.cpp +++ b/quest/src/core/validation.cpp @@ -3813,7 +3813,7 @@ void validate_newPauliStrSumParams(qindex numTerms, const char* caller) { size_t memPerTerm = sizeof(qcomp) + sizeof(PauliStr); size_t maxNumTerms = std::numeric_limits::max() / memPerTerm; tokenSubs vars = {{"${NUM_TERMS}", numTerms}, {"${NUM_BYTES_PER_TERM}", memPerTerm}, {"${MAX_NUM_TERMS}", maxNumTerms}}; - assertThat(numTerms < (qindex) maxNumTerms, report::NEW_PAULI_STR_SUM_MEM_WOULD_OVERFLOW, vars, caller); + assertThat(numTerms <= (qindex) maxNumTerms, report::NEW_PAULI_STR_SUM_MEM_WOULD_OVERFLOW, vars, caller); // attempt to fetch RAM, and simply return if we fail; if we unknowingly // didn't have enough RAM, then alloc validation will trigger later From 9a07f3f99b521a948db52754f75fd89253973e5e Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Sun, 4 Oct 2026 15:08:23 -0700 Subject: [PATCH 18/27] Correct validation error msgs --- quest/src/core/validation.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/quest/src/core/validation.cpp b/quest/src/core/validation.cpp index 2581a9892..d21988dc3 100644 --- a/quest/src/core/validation.cpp +++ b/quest/src/core/validation.cpp @@ -3404,7 +3404,7 @@ void validate_superOpDimMatchesTargs(SuperOp op, int numTargets, const char* cal if (!global_isValidationEnabled) return; - tokenSubs vars = {{"${OP_QUBITS}", op.numQubits}, {"{NUM_TARGS}", numTargets}}; + tokenSubs vars = {{"${OP_QUBITS}", op.numQubits}, {"${NUM_TARGS}", numTargets}}; assertThat(op.numQubits == numTargets, report::SUPER_OP_SIZE_MISMATCHES_NUM_TARGETS, vars, caller); } @@ -4405,7 +4405,7 @@ void validate_measurementOutcomesFitInGpuMem(Qureg qureg, int numQubits, const c size_t memNeeded = sizeof(qreal) * numOutcomes; tokenSubs vars = { - {"{NUM_QUBITS}", numQubits}, + {"${NUM_QUBITS}", numQubits}, {"${NUM_OUTCOMES}", numOutcomes}, {"${MEM_NEEDED}", memNeeded}, {"${MEM_AVAIL}", memAvail} From 610efe72fbcffd6982f8db54926565d3a7f767c8 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Sun, 4 Oct 2026 15:09:53 -0700 Subject: [PATCH 19/27] Patch maximally-distributed cuQuantum calcTotalProb --- quest/src/gpu/gpu_subroutines.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/quest/src/gpu/gpu_subroutines.cpp b/quest/src/gpu/gpu_subroutines.cpp index ac54957fd..0cb41c38d 100644 --- a/quest/src/gpu/gpu_subroutines.cpp +++ b/quest/src/gpu/gpu_subroutines.cpp @@ -1507,7 +1507,11 @@ INSTANTIATE_FUNC_OPTIMISED_FOR_NUM_TARGS( void, gpu_densmatr_partialTrace_sub, ( qreal gpu_statevec_calcTotalProb_sub(Qureg qureg) { #if QUEST_COMPILE_CUQUANTUM - return cuquantum_statevec_calcTotalProb_sub(qureg); + // cuStateVec's Z-basis routine requires a target qubit, which does not exist + // when a maximally distributed statevector stores one amplitude per node + return (qureg.logNumAmpsPerNode > 0)? + cuquantum_statevec_calcTotalProb_sub(qureg): + thrust_statevec_calcTotalProb_sub(qureg); #elif QUEST_COMPILE_CUDA return thrust_statevec_calcTotalProb_sub(qureg); From 99df89ddfff80dc8a7be865cf5bbb1c4f8a33e20 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Sun, 4 Oct 2026 15:20:37 -0700 Subject: [PATCH 20/27] Patch two's complement matrix initialisation --- quest/src/core/fastmath.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/quest/src/core/fastmath.hpp b/quest/src/core/fastmath.hpp index 79884659f..da4ef8d57 100644 --- a/quest/src/core/fastmath.hpp +++ b/quest/src/core/fastmath.hpp @@ -96,7 +96,7 @@ INLINE void fast_getSubQuregValues(qindex basisStateIndex, int* numQubitsPerSubQ if (areSigned) for (int v=0; v Date: Sun, 4 Oct 2026 15:39:50 -0700 Subject: [PATCH 21/27] Handle maximally distributed dephasing edgecase --- quest/src/core/errors.cpp | 5 +++++ quest/src/core/errors.hpp | 2 ++ quest/src/gpu/gpu_cuquantum.cuh | 5 +++++ 3 files changed, 12 insertions(+) diff --git a/quest/src/core/errors.cpp b/quest/src/core/errors.cpp index 5f98822c2..306f58ea2 100644 --- a/quest/src/core/errors.cpp +++ b/quest/src/core/errors.cpp @@ -756,6 +756,11 @@ void error_cuQuantumTempCpuAllocFailed() { raiseInternalError("Attempted allocation of temporary host-memory for a cuQuantum routine failed."); } +void error_cuQuantumIncompatibleWithMaximallyDistributedDephasing() { + + raiseInternalError("The current extreme distribution scheme (a single density matrix column per node) is incompatible with our optimised implementation of one-qubit dephasing using cuQuantum."); +} + /* diff --git a/quest/src/core/errors.hpp b/quest/src/core/errors.hpp index 7db4e69f5..d2c95a1f8 100644 --- a/quest/src/core/errors.hpp +++ b/quest/src/core/errors.hpp @@ -307,6 +307,8 @@ void error_cuQuantumInitOrFinalizedButNotCompiled(); void error_cuQuantumTempCpuAllocFailed(); +void error_cuQuantumIncompatibleWithMaximallyDistributedDephasing(); + /* diff --git a/quest/src/gpu/gpu_cuquantum.cuh b/quest/src/gpu/gpu_cuquantum.cuh index 6323f549f..44e43a97f 100644 --- a/quest/src/gpu/gpu_cuquantum.cuh +++ b/quest/src/gpu/gpu_cuquantum.cuh @@ -284,6 +284,11 @@ void cuquantum_densmatr_oneQubitDephasing_subB(Qureg qureg, int ketQubit, qreal // note it must still be a suffix qubit since cuQuantum does not know qureg is distributed int targ = qureg.logNumAmpsPerNode - 1; // leftmost suffix bra qubit + // only an extreme distribution (one density-matrix column per node) makes + // the chosen target equal the control, which cuStateVec forbids + if (targ == ketQubit) + error_cuQuantumIncompatibleWithMaximallyDistributedDephasing(); + bool conj = false; auto ctrls = lists_getList64({ketQubit}); auto states = lists_getList64({!braBit}); From fb2fe7f962fd4d1e625092a8d17bce9ae1f3482e Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Mon, 5 Oct 2026 00:17:51 -0700 Subject: [PATCH 22/27] Patch NUMA array allocation Avoid imprecise floating-point page counts, support multiword node masks, and skip binding empty page ranges --- quest/src/cpu/cpu_config.cpp | 33 +++++++++++++++++++++++---------- 1 file changed, 23 insertions(+), 10 deletions(-) diff --git a/quest/src/cpu/cpu_config.cpp b/quest/src/cpu/cpu_config.cpp index 139904f2a..34b234b89 100644 --- a/quest/src/cpu/cpu_config.cpp +++ b/quest/src/cpu/cpu_config.cpp @@ -136,7 +136,7 @@ int cpu_getCurrentNumThreads() { qindex getNumPagesToContainArray(long pageLen, qindex arrLen) { // round up to the nearest page - return static_cast(std::ceil(arrLen / (qreal) pageLen)); + return arrLen / pageLen + (arrLen % pageLen != 0); } @@ -223,6 +223,11 @@ qcomp* cpu_allocNumaArray(qindex length) { // use integer type for safe address arithmetic below uintptr_t offsetAddr = reinterpret_cast(rawAddr); + // mbind() accepts a bit array, which must span every configured NUMA node + size_t bitsPerMaskWord = 8 * sizeof(unsigned long); + size_t numMaskWords = (numNodes + bitsPerMaskWord - 1) / bitsPerMaskWord; + vector nodeMask(numMaskWords, 0); + for (int node=0, shift=numNodes; node < numNodes; ++node) { // decide number of pages to bind to NUMA node @@ -230,16 +235,24 @@ qcomp* cpu_allocNumaArray(qindex length) { qindex numPagesInNode = baseNumPagesPerNode + (shift <= 0); qindex numBytesInNode = numPagesInNode * pageSize; // validation prevents overflow - // bind those pages from the offset address to the node (identified by mask) - unsigned long nodeMask = 1UL << node; - unsigned long numBitsInMask = 8 * sizeof(nodeMask); - void* nodeAddr = reinterpret_cast(offsetAddr); - long success = mbind(nodeAddr, numBytesInNode, MPOL_BIND, &nodeMask, numBitsInMask, 0); + // skip nodes assigned no pages (which occurs when there are fewer pages than nodes) + if (numPagesInNode > 0) { + + // bind those pages from the offset address to the node (identified by mask) + size_t maskWord = node / bitsPerMaskWord; + unsigned long maskBit = 1UL << (node % bitsPerMaskWord); + nodeMask[maskWord] = maskBit; + + void* nodeAddr = reinterpret_cast(offsetAddr); + long success = mbind(nodeAddr, numBytesInNode, MPOL_BIND, nodeMask.data(), numMaskWords * bitsPerMaskWord, 0); + + nodeMask[maskWord] = 0; - // treat bind failure as internal error (even though it can result from insufficient kernel mem), - // rather than permitting silent fallback to non-NUMA awareness which might be astonishingly slow - if (success == -1) - error_numaBindingFailed(); + // treat bind failure as internal error (even though it can result from insufficient kernel mem), + // rather than permitting silent fallback to non-NUMA awareness which might be astonishingly slow + if (success == -1) + error_numaBindingFailed(); + } // prepare next node's address offsetAddr += numPagesInNode * pageSize; From 0715a68c61d288908bde1139b3863626ca8adc37 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Mon, 5 Oct 2026 14:42:10 -0700 Subject: [PATCH 23/27] Improve memory unit rounding --- quest/src/core/printer.cpp | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/quest/src/core/printer.cpp b/quest/src/core/printer.cpp index 863317c8f..a7207dcb2 100644 --- a/quest/src/core/printer.cpp +++ b/quest/src/core/printer.cpp @@ -473,18 +473,22 @@ void print_oneFewerNewlines() { string printer_getMemoryWithUnitStr(size_t numBytes) { // retain "bytes" instead of "B" since more recognisable - vector units = {"bytes", "KiB", "MiB", "GiB", "TiB", "PiB", "EiB", "ZiB", "YiB"}; + vector units = {"bytes", "KiB", "MiB", "GiB", "TiB", "PiB", "EiB"}; - // unit sizes (in bytes) may overflow; but then so too will numBytes + // unit sizes in bytes vector sizes(units.size()); for (size_t i=0; i 1 - int ind = 1; - while (numBytes > sizes[ind]) + // find the biggest unit no larger than numBytes + size_t ind = 0; + while (ind + 1 < sizes.size() && numBytes >= sizes[ind + 1]) ind++; - ind--; + + // at four significant figures, values at least halfway from 1023 to 1024 + // round to 1024 of the current unit, which is more cleanly shown as 1 of the next + if (ind + 1 < sizes.size() && numBytes >= sizes[ind + 1] - sizes[ind] / 2) + return "1 " + units[ind + 1]; // express numBytes in terms of new unit, forcefully rounding to 4 sig-figs max qreal frac = numBytes / static_cast(sizes[ind]); @@ -974,7 +978,7 @@ void populateManyColumnQcompmatr(qcompmatr &matr, qindex startRow, qindex startC // note this achieves node consensus; every node populates their matr, which is much more // demanding than only the root note, however we expect this is not a big deal; live printing // in distributed settings is exceptionally rare and will print small, tractable subsets - localiser_densmatr_getAmps(ptrs.data(), qureg, startRow, startCol, matr.size(), matr[0].capacity()); + localiser_densmatr_getAmps(ptrs.data(), qureg, startRow, startCol, matr.size(), matr[0].size()); } From d35635522f610ca1f5409e35fe6ba5518a6f26ed Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Mon, 5 Oct 2026 20:34:43 -0700 Subject: [PATCH 24/27] Patch sufficient memory validations Previous erroneous logic sometimes overestimated distributed capacity by 1-2 qubits, permitting validation to succeed when insufficient memory was available to allocate distributed structures, like Qureg and SuperOp --- quest/src/core/memory.cpp | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/quest/src/core/memory.cpp b/quest/src/core/memory.cpp index 7f11494a1..f864e1e00 100644 --- a/quest/src/core/memory.cpp +++ b/quest/src/core/memory.cpp @@ -48,17 +48,23 @@ if (hasBuffer && numNodes > 1) maxLocalNumAmps /= 2; // floors + // avoid log2(0) when even one amplitude cannot fit + if (maxLocalNumAmps == 0) + return -1; + + // distribution increases the global amplitude exponent before converting it to logical + // qubits; this order matters for density matrices, whose amplitude exponent is doubled + int maxGlobalAmpQubits = std::floor(std::log2(maxLocalNumAmps)) + logBase2(numNodes); + // density matrices require square more memory, so halve (flooring) the number of qubits - int maxLocalNumQubits = std::floor(std::log2(maxLocalNumAmps)); + int maxGlobalNumQubits = maxGlobalAmpQubits; if (isDensMatr) - maxLocalNumQubits /= 2; // floors + maxGlobalNumQubits /= 2; // floors // superoperators require square more memory still, so halve (flooring) the number of qubits if (isSuperOp) - maxLocalNumQubits /= 2; // floors + maxGlobalNumQubits /= 2; // floors - // doubling nodes permits 1 additional qubit - int maxGlobalNumQubits = maxLocalNumQubits + logBase2(numNodes); return maxGlobalNumQubits; } From 4c160b75087887698496b2d6ff4b9c537b5b0d2a Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Sun, 11 Oct 2026 15:44:15 -0400 Subject: [PATCH 25/27] Patch KrausMap memory validation which previously referenced a qindex overflow, rather than a size_t Commit also includes yet another ridiculous edge-case protection --- quest/src/core/memory.cpp | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/quest/src/core/memory.cpp b/quest/src/core/memory.cpp index f864e1e00..9ff65b05a 100644 --- a/quest/src/core/memory.cpp +++ b/quest/src/core/memory.cpp @@ -126,7 +126,7 @@ int mem_getMaxNumSuperOpQubitsBeforeIndexOverflow() { qindex mem_getMaxNumKrausMapMatricesBeforeIndexOverflow(int numQubits) { - qindex numElemPerMatrix = powerOf2(2 * numQubits); + qindex numElemPerMatrix = powerOf2(2 * numQubits); // caller prior validates no overflow qindex maxNumTotalElems = std::numeric_limits::max(); qindex maxNumMatrices = maxNumTotalElems / numElemPerMatrix; // floors return maxNumMatrices; @@ -202,9 +202,10 @@ int mem_getMaxNumSuperOpQubitsBeforeGlobalMemSizeofOverflow() { qindex mem_getMaxNumKrausMapMatricesBeforeLocalMemSizeofOverflow(int numQubits) { - qindex numMatrWithMaxTotalElems = mem_getMaxNumKrausMapMatricesBeforeIndexOverflow(numQubits); - qindex numMatrWithMaxTotalMem = numMatrWithMaxTotalElems / sizeof(qcomp); // floors - return numMatrWithMaxTotalMem; + qindex numElemsPerMatrix = powerOf2(2 * numQubits); // caller prior validates no overflow + size_t maxNumElems = std::numeric_limits::max() / sizeof(qcomp); // floors + size_t maxNumMatrices = maxNumElems / static_cast(numElemsPerMatrix); // floors + return static_cast(maxNumMatrices); } @@ -374,6 +375,11 @@ bool mem_canMatrixFitInMemory(int numQubits, bool isDense, int numNodes, qindex // work out the maximum "local" qubits that can fit in memory qindex maxLocalNumElems = memBytesPerNode / sizeof(qcomp); // floors + + // max distribution (useful for tests) has 1 amp per node, = 0 local qubits! + if (maxLocalNumElems == 0) + return false; + int maxLocalNumQubits = std::floor(std::log2(maxLocalNumElems)); // dense matrices (as opposed to diagonals) require square more memory From f39f0bc7e8c4c03b066f4ee7ae69007a0ae45180 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Sun, 11 Oct 2026 19:46:59 -0400 Subject: [PATCH 26/27] Patch MPI request leak and defend overflows --- quest/src/comm/comm_routines.cpp | 79 +++++++++++++++++++------------- quest/src/core/errors.cpp | 7 +++ quest/src/core/errors.hpp | 3 ++ 3 files changed, 57 insertions(+), 32 deletions(-) diff --git a/quest/src/comm/comm_routines.cpp b/quest/src/comm/comm_routines.cpp index 62bef2dfa..6ec3a3f1c 100644 --- a/quest/src/comm/comm_routines.cpp +++ b/quest/src/comm/comm_routines.cpp @@ -172,18 +172,23 @@ std::array dividePow2PayloadIntoMessages(qindex numAmps) { assert_commPayloadIsPowerOf2(numAmps); // use single message if possible - if (numAmps < MAX_MESSAGE_LENGTH) + if (numAmps < MAX_MESSAGE_LENGTH) { + assert_commMessageParamFitsIntoMpiInt(numAmps); return {numAmps, 1}; + } - // else, payload divides evenly between max-size messages (always fits in int) + // else, payload divides evenly between max-size messages qindex numMessages = numAmps / MAX_MESSAGE_LENGTH; // which we must be able to uniquely tag if (numMessages > getMaxNumMessages()) error_commNumMessagesExceedTagMax(); - // outputs always fit in 'int' but we force them to be 'qindex' since - // caller will multiply them with ints and could easily overflow + // ensure all outputs can safely narrow at the MPI boundary + // (they're kept qindex to avoid overflows in algebra) + assert_commMessageParamFitsIntoMpiInt(MAX_MESSAGE_LENGTH); + assert_commMessageParamFitsIntoMpiInt(numMessages); + return {MAX_MESSAGE_LENGTH, numMessages}; } @@ -191,19 +196,26 @@ std::array dividePow2PayloadIntoMessages(qindex numAmps) { std::array dividePayloadIntoMessages(qindex numAmps) { // use single message if possible - if (numAmps < MAX_MESSAGE_LENGTH) + if (numAmps < MAX_MESSAGE_LENGTH) { + assert_commMessageParamFitsIntoMpiInt(numAmps); return {numAmps, 1, 0}; + } // else, use as many max-size messages as possible, and one smaller msg qindex numMaxSizeMsgs = numAmps / MAX_MESSAGE_LENGTH; // floors qindex remainingMsgSize = numAmps - numMaxSizeMsgs * MAX_MESSAGE_LENGTH; // all of which we must be able to uniquely tag - if (numMaxSizeMsgs + 1 > getMaxNumMessages()) + qindex numMessages = numMaxSizeMsgs + (remainingMsgSize > 0); + if (numMessages > getMaxNumMessages()) error_commNumMessagesExceedTagMax(); - // outputs always fit in 'int' but we force them to be 'qindex' since - // caller will multiply them with ints and could easily overflow + // ensure all outputs can safely narrow at the MPI boundary + // (they're kept qindex to avoid overflows in algebra) + assert_commMessageParamFitsIntoMpiInt(MAX_MESSAGE_LENGTH); + assert_commMessageParamFitsIntoMpiInt(numMaxSizeMsgs); + assert_commMessageParamFitsIntoMpiInt(remainingMsgSize); + return {MAX_MESSAGE_LENGTH, numMaxSizeMsgs, remainingMsgSize}; } @@ -221,20 +233,21 @@ void exchangeArrays(qcomp* send, qcomp* recv, qindex numElems, int pairRank) { // each message is asynchronously dispatched with a final wait, as per arxiv.org/abs/2308.07402 - // we will send payload in multiple asynch messages (create two requests per msg for subsequent synch) + // we will send payload in multiple asynch messages auto [messageSize, numMessages] = dividePow2PayloadIntoMessages(numElems); - vector requests(2*numMessages, MPI_REQUEST_NULL); + vector recvRequests(numMessages, MPI_REQUEST_NULL); + vector sendRequests(numMessages, MPI_REQUEST_NULL); // asynchronously exchange the messages (effecting MPI_Isendrecv), using unique tags // so that messages are permitted to arrive out-of-order (supporting UCX adaptive-routing) for (qindex m=0; m(m); // gauranteed int, but m*messageSize needs qindex - MPI_Irecv(&recv[m*messageSize], messageSize, MPI_QCOMP, pairRank, tag, mpiComm, &requests[2*m]); - MPI_Isend(&send[m*messageSize], messageSize, MPI_QCOMP, pairRank, tag, mpiComm, &requests[2*m+1]); + MPI_Irecv(&recv[m * messageSize], messageSize, MPI_QCOMP, pairRank, /*tag=*/m, mpiComm, &recvRequests[m]); + MPI_Isend(&send[m * messageSize], messageSize, MPI_QCOMP, pairRank, /*tag=*/m, mpiComm, &sendRequests[m]); } // wait for all exchanges to complete (MPI will automatically free the request memory) - MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE); + MPI_Waitall(numMessages, recvRequests.data(), MPI_STATUSES_IGNORE); + MPI_Waitall(numMessages, sendRequests.data(), MPI_STATUSES_IGNORE); #else error_commButEnvNotDistributed(); @@ -253,16 +266,17 @@ void asynchSendArray(qcomp* send, qindex numElems, int pairRank) { MPI_Comm mpiComm = comm_getMpiComm(); - // we will not track nor wait for the asynch send; instead, the caller will later comm_sync() - MPI_Request nullReq = MPI_REQUEST_NULL; + // we will not track nor wait for the asynch sends; the receiver completes every + // matching receive before the caller's later comm_sync(), so req begins NULL + MPI_Request detachedReq = MPI_REQUEST_NULL; // divide the data into multiple messages auto [messageSize, numMessages] = dividePow2PayloadIntoMessages(numElems); - // asynchronously send the uniquely-tagged messages + // asynchronously send the uniquely-tagged messages, detaching requests so MPI can later release for (qindex m=0; m(m); // gauranteed int, but m*messageSize needs qindex - MPI_Isend(&send[m*messageSize], messageSize, MPI_QCOMP, pairRank, tag, mpiComm, &nullReq); + MPI_Isend(&send[m*messageSize], messageSize, MPI_QCOMP, pairRank, /*tag=*/m, mpiComm, &detachedReq); + MPI_Request_free(&detachedReq); } #else @@ -283,13 +297,11 @@ void receiveArray(qcomp* dest, qindex numElems, int pairRank) { vector requests(numMessages, MPI_REQUEST_NULL); // listen to receive each uniquely-tagged message asynchronously (as per arxiv.org/abs/2308.07402) - for (qindex m=0; m(m); // gauranteed int, but m*messageSize needs qindex - MPI_Irecv(&dest[m*messageSize], messageSize, MPI_QCOMP, pairRank, tag, mpiComm, &requests[m]); - } + for (qindex m=0; m requests; + qindex numRequests = 0; // each node broadcasts their partition (in-turn, but each is asynch)... for (int sendRank=0; sendRank 0); + assert_commMessageParamFitsIntoMpiInt(numRequests); + // these involve big asynch messages (could be 'numSend', or the max message size) - for (int m=0; m(m); + for (qindex m=0; m #include +#include using std::string; @@ -236,6 +237,12 @@ void assert_commPayloadIsPowerOf2(qindex numAmps) { raiseInternalError("A communication function was given a payload which was unexpectedly not a power of 2, as implied by preconditions."); } +void assert_commMessageParamFitsIntoMpiInt(qindex value) { + + if (value > std::numeric_limits::max()) + raiseInternalError("A function attempted to communicate with a message parameter which exceeded the maximum int accepted by MPI."); +} + void assert_commQuregIsDistributed(Qureg qureg) { if (!qureg.isDistributed) diff --git a/quest/src/core/errors.hpp b/quest/src/core/errors.hpp index d2c95a1f8..753a6d35b 100644 --- a/quest/src/core/errors.hpp +++ b/quest/src/core/errors.hpp @@ -21,6 +21,7 @@ #include "quest/include/qureg.h" #include "quest/include/matrices.h" +#include #include using std::string; @@ -111,6 +112,8 @@ void assert_commBoundsAreValid(Qureg qureg, qindex sendInd, qindex recvInd, qind void assert_commPayloadIsPowerOf2(qindex numAmps); +void assert_commMessageParamFitsIntoMpiInt(qindex value); + void assert_commQuregIsDistributed(Qureg qureg); void assert_commFullStateDiagMatrIsDistributed(FullStateDiagMatr matr); From 4e1ea7a45814dc83a8d952f14aa8b58aa2188d17 Mon Sep 17 00:00:00 2001 From: Tyson Jones Date: Sun, 11 Oct 2026 20:13:01 -0400 Subject: [PATCH 27/27] Patch MPI calcProbsOfAllMultiQubitOutcomes which failed to update the full array on each node when the number of local probabilities exceeded max-int (about 32 qubits worth) --- quest/src/comm/comm_routines.cpp | 57 +++++++++++++++++++++++--------- 1 file changed, 41 insertions(+), 16 deletions(-) diff --git a/quest/src/comm/comm_routines.cpp b/quest/src/comm/comm_routines.cpp index 6ec3a3f1c..900ece98c 100644 --- a/quest/src/comm/comm_routines.cpp +++ b/quest/src/comm/comm_routines.cpp @@ -168,17 +168,17 @@ int getMaxNumMessages() { } -std::array dividePow2PayloadIntoMessages(qindex numAmps) { - assert_commPayloadIsPowerOf2(numAmps); +std::array dividePow2PayloadIntoMessages(qindex messageLength) { + assert_commPayloadIsPowerOf2(messageLength); // use single message if possible - if (numAmps < MAX_MESSAGE_LENGTH) { - assert_commMessageParamFitsIntoMpiInt(numAmps); - return {numAmps, 1}; + if (messageLength < MAX_MESSAGE_LENGTH) { + assert_commMessageParamFitsIntoMpiInt(messageLength); + return {messageLength, 1}; } // else, payload divides evenly between max-size messages - qindex numMessages = numAmps / MAX_MESSAGE_LENGTH; + qindex numMessages = messageLength / MAX_MESSAGE_LENGTH; // which we must be able to uniquely tag if (numMessages > getMaxNumMessages()) @@ -193,17 +193,17 @@ std::array dividePow2PayloadIntoMessages(qindex numAmps) { } -std::array dividePayloadIntoMessages(qindex numAmps) { +std::array dividePayloadIntoMessages(qindex messageLength) { // use single message if possible - if (numAmps < MAX_MESSAGE_LENGTH) { - assert_commMessageParamFitsIntoMpiInt(numAmps); - return {numAmps, 1, 0}; + if (messageLength < MAX_MESSAGE_LENGTH) { + assert_commMessageParamFitsIntoMpiInt(messageLength); + return {messageLength, 1, 0}; } // else, use as many max-size messages as possible, and one smaller msg - qindex numMaxSizeMsgs = numAmps / MAX_MESSAGE_LENGTH; // floors - qindex remainingMsgSize = numAmps - numMaxSizeMsgs * MAX_MESSAGE_LENGTH; + qindex numMaxSizeMsgs = messageLength / MAX_MESSAGE_LENGTH; // floors + qindex remainingMsgSize = messageLength - numMaxSizeMsgs * MAX_MESSAGE_LENGTH; // all of which we must be able to uniquely tag qindex numMessages = numMaxSizeMsgs + (remainingMsgSize > 0); @@ -676,13 +676,23 @@ void comm_sendAmpsToRoot(int sendRank, qcomp* send, qcomp* recv, qindex numAmps) #if QUEST_COMPILE_MPI MPI_Comm mpiComm = comm_getMpiComm(); - - // only the sender and root nodes need to continue - int recvRank = ROOT_RANK; int myRank = comm_getRank(); + int recvRank = ROOT_RANK; + + // only the two sender and recv nodes need to continue if (myRank != sendRank && myRank != recvRank) return; + // avoid issuing an unmatched self-send when root already owns the payload + if (sendRank == recvRank) { + + // ensure recvRank has its 'send' overwritten with 'recv' + if (myRank == recvRank && send != recv) + cpu_copyArray(recv, send, numAmps); + + return; + } + // create an MPI_Request for every asynch MPI call auto [messageSize, numMessages] = dividePow2PayloadIntoMessages(numAmps); vector requests(numMessages, MPI_REQUEST_NULL); @@ -706,6 +716,9 @@ void comm_sendAmpsToRoot(int sendRank, qcomp* send, qcomp* recv, qindex numAmps) void comm_broadcastIntsFromRoot(int* arr, qindex length) { #if QUEST_COMPILE_MPI + // Bcast can only handle up to int-max; we presently only ever pass length=1 + assert_commMessageParamFitsIntoMpiInt(length); + int sendRank = ROOT_RANK; MPI_Bcast(arr, length, MPI_INT, sendRank, comm_getMpiComm()); @@ -718,6 +731,9 @@ void comm_broadcastIntsFromRoot(int* arr, qindex length) { void comm_broadcastUnsignedsFromRoot(unsigned* arr, qindex length) { #if QUEST_COMPILE_MPI + // Bcast can only handle up to int-max; we presently only ever pass length=1 + assert_commMessageParamFitsIntoMpiInt(length); + int sendRank = ROOT_RANK; MPI_Bcast(arr, length, MPI_UNSIGNED, sendRank, comm_getMpiComm()); @@ -768,7 +784,16 @@ void comm_reduceReal(qreal* localReal) { void comm_reduceReals(qreal* localReals, qindex numLocalReals) { #if QUEST_COMPILE_MPI - MPI_Allreduce(MPI_IN_PLACE, localReals, numLocalReals, MPI_QREAL, MPI_SUM, comm_getMpiComm()); + // since numLocalReals can exceed max-int, we must divide into max-length messages + auto [bigMsgSize, numBigMsgs, remMsgSize] = dividePayloadIntoMessages(numLocalReals); + + // send max-length messages + for (qindex m=0; m 0) + MPI_Allreduce(MPI_IN_PLACE, &localReals[numBigMsgs*bigMsgSize], remMsgSize, MPI_QREAL, MPI_SUM, comm_getMpiComm()); #else error_commButEnvNotDistributed();