Silent wrong answers on GPU from 32 qubits upward: kernel launch failures are never checked

[ipasichnyk]

On an AMD GPU, any QuEST kernel that needs 2^32 or more threads fails to launch with hipErrorInvalidConfiguration. QuEST does not check kernel launch errors, so the operation silently does nothing, the API returns success, and the caller gets a wrong answer with no diagnostic. This begins at 32 qubits for uncontrolled diagonal operations and 33 qubits for one-target dense gates.

Environment

• QuEST aa19534 (v4.3-release), unmodified
• AMD Instinct MI355X (gfx950), 288 GiB, single GPU, no MPI
• ROCm 7.2.1, HIP 7.2.53211, amdclang++, Release, QUEST_FLOAT_PRECISION=2

Reproducer

#include <cstdlib>
#include <cmath>
#include <iostream>
#include "quest.h"

// declared directly to avoid HIP headers on the host compile line
extern "C" int hipGetLastError(void);
extern "C" const char* hipGetErrorString(int);

static int failures = 0;

static void check(const char* what, qreal got, qreal want) {
    int ok = std::abs(got - want) < 1e-9;
    failures += !ok;
    std::cout << (ok ? "  pass  " : "  FAIL  ") << what
              << " = " << got << "  (expected " << want << ")" << std::endl;
}

int main(int argc, char* argv[]) {
    int n = (argc > 1) ? std::atoi(argv[1]) : 32;

    initQuESTEnv();
    Qureg q = createQureg(n);

    qindex numAmps = (qindex) 1 << n;
    qindex half    = (qindex) 1 << (n - 1);

    // [A] uncontrolled diagonal -> 2^n threads
    std::cout << "[A] applyPauliZ (" << numAmps << " threads)" << std::endl;
    initPlusState(q);
    applyPauliZ(q, 0);

    int err = hipGetLastError();
    std::cout << "  returned normally; hipGetLastError() = "
              << hipGetErrorString(err) << std::endl;

    qreal amp = 1 / std::sqrt((qreal) numAmps);
    check("amp[1]", std::real(getQuregAmp(q, 1)), -amp);   // Z negates odd indices

    // [B] one-target dense -> 2^(n-1) threads
    std::cout << "[B] applyHadamard (" << half << " threads)" << std::endl;
    initZeroState(q);
    applyHadamard(q, n - 1);

    qreal invRoot2 = 0.70710678118654752;
    check("amp[0]",    std::real(getQuregAmp(q, 0)),    invRoot2);
    check("amp[half]", std::real(getQuregAmp(q, half)), invRoot2);

    std::cout << (failures ? "RESULT: FAIL" : "RESULT: PASS") << std::endl;

    destroyQureg(q);
    finalizeQuESTEnv();
    return failures ? EXIT_FAILURE : EXIT_SUCCESS;
}

Built with -DUSER_SOURCE_NAMES=repro.cpp -DUSER_OUTPUT_EXE_NAME=repro, plus -DCMAKE_EXE_LINKER_FLAGS="-L$ROCM_PATH/lib -lamdhip64" for the two extern "C" declarations.

Observed

$ ./repro 31
[A] applyPauliZ (2147483648 threads)
  returned normally; hipGetLastError() = no error
  pass  amp[1] = -2.15792e-05  (expected -2.15792e-05)
[B] applyHadamard (1073741824 threads)
  pass  amp[0] = 0.707107  (expected 0.707107)
  pass  amp[half] = 0.707107  (expected 0.707107)
RESULT: PASS

$ ./repro 32
[A] applyPauliZ (4294967296 threads)
  returned normally; hipGetLastError() = invalid configuration argument
  FAIL  amp[1] = 1.52588e-05  (expected -1.52588e-05)
[B] applyHadamard (2147483648 threads)
  pass  amp[0] = 0.707107  (expected 0.707107)
  pass  amp[half] = 0.707107  (expected 0.707107)
RESULT: FAIL

$ ./repro 33
[A] applyPauliZ (8589934592 threads)
  returned normally; hipGetLastError() = invalid configuration argument
  FAIL  amp[1] = 1.07896e-05  (expected -1.07896e-05)
[B] applyHadamard (4294967296 threads)
  FAIL  amp[0] = 1  (expected 0.707107)
  FAIL  amp[half] = 0  (expected 0.707107)
RESULT: FAIL

At 32 qubits applyPauliZ is a complete no-op yet returns normally. At 33 qubits applyHadamard also becomes a no-op: amp[0] = 1, amp[half] = 0 means the register is untouched.

Cause

ROCm limits a grid to 2^32 - 1 total work-items; rocminfo reports Grid Max Size: 4294967295(0xffffffff). QuEST launches one thread per iteration:
gpu_statevector_anyCtrlAnyTargZOrPhaseGadget_sub: numThreads = qureg.numAmpsPerNode / powerOf2(ctrls.size()), so 2^n uncontrolled
gpu_statevec_anyCtrlOneTargDenseMatr_subA: numThreads = qureg.numAmpsPerNode / powerOf2(ctrls.size() + 1), so 2^(n-1) uncontrolled
At 32 qubits the first is exactly 2^32 and the launch is rejected. The failure is silent because kernel launches are not error-checked. CUDA_CHECK is applied consistently to API calls, but no <<<...>>> site is followed by a cudaGetLastError() check.
Thrust-based operations are unaffected: initPlusState and initZeroState are correct at 2^33 amplitudes, since Thrust chunks large ranges internally. Only the hand-written kernels fail.

[otbrown]

Thanks for this @ipasichnyk! Aside from correcting the error-handling, it would also be nice if the AMD GPU implementation ran correctly above 2^32 elements per GPU.

I'll have a chat with the others so we can take a view on whether that's a v4.3.0 problem, or a v4.3.1 problem.

[TysonRayJones]

This results from QuEST's design of "do the minimum amount of work per thread". It's a very naive approach, but beautifully simple and makes for very easy reading of kernels:

__global__ void kernel_densmatr_oneQubitDepolarising_subA(
    gpu_qcomp* amps, qindex numThreads, 
    int ketQubit, int braQubit, qreal facAA, qreal facBB, qreal facAB
) {
    GET_THREAD_IND(n, numThreads);

    // i00 = nth local index where both qubits are 0
    qindex i00 = insertTwoBits(n, braQubit, 0, ketQubit, 0);
    qindex i01 = flipBit(i00, ketQubit);
    qindex i10 = flipBit(i00, braQubit);
    qindex i11 = flipBit(i01, braQubit);

    // modify 4 amps, mixing a pair, and scaling the other
    gpu_qcomp amp00 = amps[i00];
    amps[i00] = (facAA * amp00) + (facBB * amps[i11]);
    amps[i01] *= facAB;
    amps[i10] *= facAB;
    amps[i11] = (facAA * amps[i11]) + (facBB * amp00);
}

Alas it appears there is no way to retain this simplicity, if we indeed wish to support AMD processing of >2^32 amplitudes per node. We will have to introduce looping within a kernel, the associated boilerplate, and tedium of selecting related hyperparameters 😢 That's a substantial refactor that cannot be done this release (or so I'd expect).

What we can do in the near term is:

1. Add the missing kernel-launch checks - I have no idea how I neglected them! Sincerely sorry for the enormous gap in my touted "defensive design"! 😓
2. Add a temporary internal error inside the create() functions when trying to alloc more than 2^32-sized structures, when using an AMD GPU

These should be sufficiently isolated from any logic that they shouldn't introduce any bugs, or jeopardise the existing v4.3 tests. I'll make the PR, and @otbrown can decide if worth squeezing into this release.

Thanks very much for the thorough diagnostics, and apologies for the bug!

[TysonRayJones]

Addressed respectively in #816 and #817. I've manually triggered/tested all internal error handling paths (CUDA, cuStateVec, Thrust, and AMD bug guard). In my opinion they're safe to merge without retesting, but I leave it up to @otbrown !