diff --git a/quest/src/core/errors.cpp b/quest/src/core/errors.cpp
index 807cad105..abe1bce92 100644
--- a/quest/src/core/errors.cpp
+++ b/quest/src/core/errors.cpp
@@ -694,11 +694,13 @@ void assert_gpuHasBeenBound(bool isBound) {
  * CUDA ERRORS
  */
 
-void error_cudaCallFailed(const char* msg, const char* func, const char* caller, const char* file, int line) {
+void internal_cudaLibCallFailed(const char* libname, const char* msg, const char* func, const char* caller, const char* file, int line) {
 
     // using operator overloads to cast const char[] literals to std::string, to concat with const char*.
     string err = "";
-    err += "A CUDA (or cuQuantum) API function (\"";
+    err += "A ";
+    err += libname;
+    err += " API function (\"";
     err += func;
     err += "\", called by \"";
     err += caller;
@@ -712,21 +714,39 @@ void error_cudaCallFailed(const char* msg, const char* func, const char* caller,
     raiseInternalError(err);
 }
 
+void error_cudaCallFailed(const char* msg, const char* func, const char* caller, const char* file, int line) {
+
+    internal_cudaLibCallFailed("CUDA", msg, func, caller, file, line);
+}
+
 void error_cudaEncounteredIrrecoverableError() {
 
     raiseInternalError("The CUDA API encountered an irrecoverable \"sticky\" error which was attemptedly cleared as if it were non-sticky.");
 }
 
+void error_cudaKernelLaunchFailed(const char* caller, const char* cudaErrMsg) {
+
+    string err = "";
+    err += "A CUDA kernel invoked within '";
+    err += caller;
+    err += "' failed to launch - or a prior kernel called from elsewhere asynchronously failed -";
+    err += " with CUDA error message: \"";
+    err += cudaErrMsg;
+    err += "\". ";
+    raiseInternalError(err);
+}
+
+// Looking for assert_lastKernelLaunchSucceeded(const char*)? It's in gpu_config :^)
+
 
 
 /*
  * THRUST ERRORS
  */
 
+void error_thrustCallFailed(const char* msg, const char* func, const char* caller, const char* file, int line) {
 
-void error_thrustTempGpuAllocFailed() {
-
-    raiseInternalError("Thrust failed to allocate temporary GPU memory.");
+    internal_cudaLibCallFailed("Thrust", msg, func, caller, file, line);
 }
 
 
@@ -735,6 +755,11 @@ void error_thrustTempGpuAllocFailed() {
  * CUQUANTUM ERRORS
  */
 
+void error_cuQuantumCallFailed(const char* msg, const char* func, const char* caller, const char* file, int line) {
+
+    internal_cudaLibCallFailed("cuQuantum (specifically cuStateVec)", msg, func, caller, file, line);
+}
+
 void error_cuQuantumInitOrFinalizedButNotCompiled() {
 
     raiseInternalError("Attempted to initialise or finalise cuQuantum, but cuQuantum was not compiled.");
diff --git a/quest/src/core/errors.hpp b/quest/src/core/errors.hpp
index f91f890b0..a1f20a60f 100644
--- a/quest/src/core/errors.hpp
+++ b/quest/src/core/errors.hpp
@@ -281,13 +281,17 @@ void error_cudaCallFailed(const char* msg, const char* func, const char* caller,
 
 void error_cudaEncounteredIrrecoverableError();
 
+void error_cudaKernelLaunchFailed(const char* caller, const char* cudaErrMsg);
+
+// Looking for assert_lastKernelLaunchSucceeded(const char*)? It's in gpu_config :^)
+
 
 
 /*
  * THRUST ERRORS
  */
 
-void error_thrustTempGpuAllocFailed();
+void error_thrustCallFailed(const char* msg, const char* func, const char* caller, const char* file, int line);
 
 
 
@@ -295,6 +299,8 @@ void error_thrustTempGpuAllocFailed();
  * CUQUANTUM ERRORS
  */
 
+void error_cuQuantumCallFailed(const char* msg, const char* func, const char* caller, const char* file, int line);
+
 void error_cuQuantumInitOrFinalizedButNotCompiled();
 
 void error_cuQuantumTempCpuAllocFailed();
diff --git a/quest/src/gpu/gpu_config.cpp b/quest/src/gpu/gpu_config.cpp
index 001cc62c0..7f63b831e 100644
--- a/quest/src/gpu/gpu_config.cpp
+++ b/quest/src/gpu/gpu_config.cpp
@@ -91,6 +91,17 @@ void clearPossibleCudaError() {
         error_cudaEncounteredIrrecoverableError();
 }
 
+void assert_lastKernelLaunchSucceeded(const char* caller) {
+
+    // note that we are only checking the kernel LAUNCH succeeded; it remains
+    // possible for the kernel to subsequently fail, which would only be detected
+    // at a subsequent cudaGetLastError(), or a cudaDeviceSynchronize().
+
+    cudaError_t status = cudaGetLastError();
+    if (status != cudaSuccess)
+        error_cudaKernelLaunchFailed(caller, cudaGetErrorString(status));
+}
+
 #endif
 
 
@@ -360,7 +371,7 @@ int gpu_getMaxNumThreadsPerBlock() {
 #if QUEST_COMPILE_CUDA
 
     cudaDeviceProp prop;
-    cudaGetDeviceProperties(&prop, getBoundGpuId());
+    CUDA_CHECK( cudaGetDeviceProperties(&prop, getBoundGpuId()) );
     return prop.maxThreadsPerBlock; // HIP compatible
 
 #else
diff --git a/quest/src/gpu/gpu_config.hpp b/quest/src/gpu/gpu_config.hpp
index 98cb9c8a3..ff2b83edd 100644
--- a/quest/src/gpu/gpu_config.hpp
+++ b/quest/src/gpu/gpu_config.hpp
@@ -40,6 +40,8 @@ constexpr int gpu_HIP_WARP_SIZE = 64;
 
 void assertCudaCallSucceeded(int code, const char* call, const char* caller, const char* file, int line);
 
+void assert_lastKernelLaunchSucceeded(const char* caller);
+
 #endif 
 
 
diff --git a/quest/src/gpu/gpu_cuquantum.cuh b/quest/src/gpu/gpu_cuquantum.cuh
index 6323f549f..61b1102f9 100644
--- a/quest/src/gpu/gpu_cuquantum.cuh
+++ b/quest/src/gpu/gpu_cuquantum.cuh
@@ -44,6 +44,7 @@
 
 #include "quest/include/precision.h"
 
+#include "quest/src/core/errors.hpp"
 #include "quest/src/core/lists.hpp"
 #include "quest/src/core/utilities.hpp"
 #include "quest/src/gpu/gpu_config.hpp"
@@ -77,6 +78,21 @@ using std::vector;
 
 
 
+/*
+ * CUSTATEVEC ERROR HANDLING
+ */
+
+inline void assertCuStateVecCallSucceeded(custatevecStatus_t status, const char* call, const char* caller, const char* file, int line) {
+
+    if (status != CUSTATEVEC_STATUS_SUCCESS)
+        error_cuQuantumCallFailed(custatevecGetErrorString(status), call, caller, file, line);
+}
+
+#define CUSV_CHECK(cmd) \
+    assertCuStateVecCallSucceeded((cmd), #cmd, __func__, __FILE__, __LINE__)
+
+
+
 /*
  * ENVIRONMENT MANAGEMENT
  */
@@ -143,7 +159,7 @@ void gpu_initCuQuantum() {
     // prior validation prevent it (disabled by an environment variable)
 
     // create new stream and cuQuantum handle, binding to global config
-    CUDA_CHECK( custatevecCreate(&config.handle) );
+    CUSV_CHECK( custatevecCreate(&config.handle) );
     CUDA_CHECK( cudaStreamCreate(&config.stream) );
 
     // get and configure existing memory pool (for later automatic alloc/dealloc of gate matrices)
@@ -157,15 +173,15 @@ void gpu_initCuQuantum() {
     strcpy(config.memhandler.name, "mempool");
 
     // bind memory handler and stream to cuQuantum handle
-    CUDA_CHECK( custatevecSetDeviceMemHandler(config.handle, &config.memhandler) );
-    CUDA_CHECK( custatevecSetStream(config.handle, config.stream) );
+    CUSV_CHECK( custatevecSetDeviceMemHandler(config.handle, &config.memhandler) );
+    CUSV_CHECK( custatevecSetStream(config.handle, config.stream) );
 }
 
 
 void gpu_finalizeCuQuantum() {
 
     CUDA_CHECK( cudaStreamDestroy(config.stream) );
-    CUDA_CHECK( custatevecDestroy(config.handle) );
+    CUSV_CHECK( custatevecDestroy(config.handle) );
 }
 
 
@@ -181,7 +197,7 @@ void cuquantum_statevec_anyCtrlSwap_subA(Qureg qureg, ConstList64 ctrls, ConstLi
     int2 targPairs[] = {{targ1, targ2}};;
     int numTargPairs = 1;
 
-    CUDA_CHECK( custatevecSwapIndexBits(
+    CUSV_CHECK( custatevecSwapIndexBits(
         config.handle,
         getGpuQcompPtr(qureg.gpuAmps), CUQUANTUM_QCOMP, qureg.logNumAmpsPerNode,
         targPairs, numTargPairs,
@@ -209,7 +225,7 @@ void cuquantum_statevec_anyCtrlAnyTargDenseMatrix_subA(Qureg qureg, ConstList64
     void* work = nullptr;
     size_t workSize = 0;
 
-    CUDA_CHECK( custatevecApplyMatrix(
+    CUSV_CHECK( custatevecApplyMatrix(
         config.handle, 
         getGpuQcompPtr(qureg.gpuAmps), CUQUANTUM_QCOMP, qureg.logNumAmpsPerNode, 
         flatMatrElems, CUQUANTUM_QCOMP, CUSTATEVEC_MATRIX_LAYOUT_ROW, applyAdj, 
@@ -238,7 +254,7 @@ void cuquantum_statevec_anyCtrlAnyTargDiagMatr_sub(Qureg qureg, ConstList64 ctrl
     void* work = nullptr;
     size_t workSize = 0;
 
-    CUDA_CHECK( custatevecApplyGeneralizedPermutationMatrix(
+    CUSV_CHECK( custatevecApplyGeneralizedPermutationMatrix(
         config.handle,
         getGpuQcompPtr(qureg.gpuAmps), CUQUANTUM_QCOMP, qureg.logNumAmpsPerNode,
         perm, flatMatrElems, CUQUANTUM_QCOMP, adj, 
@@ -335,7 +351,7 @@ qreal cuquantum_statevec_calcTotalProb_sub(Qureg qureg) {
     int qubit = qureg.logNumAmpsPerNode - 1;
     int numQubits = 1;
 
-    CUDA_CHECK( custatevecAbs2SumOnZBasis(
+    CUSV_CHECK( custatevecAbs2SumOnZBasis(
         config.handle,
         getGpuQcompPtr(qureg.gpuAmps), CUQUANTUM_QCOMP, qureg.logNumAmpsPerNode,
         &prob0, &prob1, &qubit, numQubits ) );
@@ -350,7 +366,7 @@ qreal cuquantum_statevec_calcProbOfMultiQubitOutcome_sub(Qureg qureg, ConstList6
     // cuQuantum probabilities are always double
     double prob;
 
-    CUDA_CHECK( custatevecAbs2SumArray(
+    CUSV_CHECK( custatevecAbs2SumArray(
         config.handle,
         getGpuQcompPtr(qureg.gpuAmps), CUQUANTUM_QCOMP, qureg.logNumAmpsPerNode,
         &prob, nullptr, 0, outcomes.data(), qubits.data(), qubits.size()) );
@@ -371,7 +387,7 @@ void cuquantum_statevec_calcProbsOfAllMultiQubitOutcomes_sub(qreal* outProbs, Qu
         double* outPtr = tmpProbs.data();
     #endif
 
-    CUDA_CHECK( custatevecAbs2SumArray(
+    CUSV_CHECK( custatevecAbs2SumArray(
         config.handle,
         getGpuQcompPtr(qureg.gpuAmps), CUQUANTUM_QCOMP, qureg.logNumAmpsPerNode,
         outPtr, qubits.data(), qubits.size(), nullptr, nullptr, 0) );
@@ -413,7 +429,7 @@ qreal cuquantum_statevec_calcExpecPauliStr_subA(Qureg qureg, ConstList64 x, Cons
     // cuStateVec output is always double
     double value = 0;
 
-    CUDA_CHECK( custatevecComputeExpectationsOnPauliBasis(
+    CUSV_CHECK( custatevecComputeExpectationsOnPauliBasis(
         config.handle,
         getGpuQcompPtr(qureg.gpuAmps), CUQUANTUM_QCOMP, qureg.logNumAmpsPerNode,
         &value, termPaulis, numTerms, termTargets, numPaulisPerTerm) );
@@ -437,11 +453,11 @@ qreal cuquantum_statevec_calcExpecAnyTargZ_sub(Qureg qureg, ConstList64 targs) {
 
 void cuquantum_statevec_multiQubitProjector_sub(Qureg qureg, ConstList64 qubits, ConstList64 outcomes, qreal prob) {
 
-    CUDA_CHECK( custatevecCollapseByBitString(
+    CUSV_CHECK( custatevecCollapseByBitString(
         config.handle,
         getGpuQcompPtr(qureg.gpuAmps), CUQUANTUM_QCOMP, qureg.logNumAmpsPerNode,
         outcomes.data(), qubits.data(), qubits.size(), prob) );
 }
 
 
-#endif // GPU_CUQUANTUM_HPP
\ No newline at end of file
+#endif // GPU_CUQUANTUM_HPP
diff --git a/quest/src/gpu/gpu_subroutines.cpp b/quest/src/gpu/gpu_subroutines.cpp
index e9936c9c3..e367a121f 100644
--- a/quest/src/gpu/gpu_subroutines.cpp
+++ b/quest/src/gpu/gpu_subroutines.cpp
@@ -152,6 +152,7 @@ qindex gpu_statevec_packAmpsIntoBuffer(Qureg qureg, ConstList64 qubits, ConstLis
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + sendInd, numThreads, 
         sortedQubits, qubitStateMask
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
     // return the number of packed amps
     return numThreads;
@@ -178,6 +179,7 @@ qindex gpu_statevec_packPairSummedAmpsIntoBuffer(Qureg qureg, int qubit1, int qu
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + sendInd, numThreads, 
         qubit1, qubit2, qubit3, bit2
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
     // return the number of packed amps
     return numThreads;
@@ -220,6 +222,7 @@ void gpu_statevec_anyCtrlSwap_subA(Qureg qureg, ConstList64 ctrls, ConstList64 c
         getGpuQcompPtr(qureg.gpuAmps), numThreads, 
         sortedQubits, qubitStateMask, targ1, targ2
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -246,6 +249,7 @@ void gpu_statevec_anyCtrlSwap_subB(Qureg qureg, ConstList64 ctrls, ConstList64 c
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + recvInd, numThreads, 
         sortedCtrls, ctrlStateMask
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -272,6 +276,7 @@ void gpu_statevec_anyCtrlSwap_subC(Qureg qureg, ConstList64 ctrls, ConstList64 c
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + recvInd, numThreads, 
         sortedQubits, qubitStateMask
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -318,6 +323,7 @@ void gpu_statevec_anyCtrlOneTargDenseMatr_subA(Qureg qureg, ConstList64 ctrls, C
         qubitStateMask, targ, 
         m00, m01, m10, m11
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -345,6 +351,7 @@ void gpu_statevec_anyCtrlOneTargDenseMatr_subB(Qureg qureg, ConstList64 ctrls, C
         sortedCtrls, ctrlStateMask, 
         getGpuQcomp(fac0), getGpuQcomp(fac1)
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -392,6 +399,7 @@ void gpu_statevec_anyCtrlTwoTargDenseMatr_sub(Qureg qureg, ConstList64 ctrls, Co
         m[0], m[1], m[2],  m[3],  m[4],  m[5],  m[6],  m[7],
         m[8], m[9], m[10], m[11], m[12], m[13], m[14], m[15]
     );
+    assert_lastKernelLaunchSucceeded(__func__);
  
 #else
     error_gpuSimButGpuNotCompiled();
@@ -482,6 +490,7 @@ void gpu_statevec_anyCtrlAnyTargDenseMatr_sub(Qureg qureg, ConstList64 ctrls, Co
                 sortedQubits, qubitStateMask, 
                 targs, matrPtr
         );
+        assert_lastKernelLaunchSucceeded(__func__);
 
     } else {
 
@@ -522,6 +531,7 @@ void gpu_statevec_anyCtrlAnyTargDenseMatr_sub(Qureg qureg, ConstList64 ctrls, Co
                 sortedQubits, qubitStateMask, 
                 targs, powerOf2(targs.size()), matrPtr
         );
+        assert_lastKernelLaunchSucceeded(__func__);
     }
 
 #else
@@ -589,6 +599,7 @@ void gpu_statevec_anyCtrlOneTargDiagMatr_sub(Qureg qureg, ConstList64 ctrls, Con
         getGpuQcompPtr(qureg.gpuAmps), numThreads, qureg.rank, qureg.logNumAmpsPerNode, sortedCtrls,
        	ctrlStateMask, targ, elems[0], elems[1]
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
     // explicitly return to avoid runtime error below
     return;
@@ -661,6 +672,7 @@ void gpu_statevec_anyCtrlTwoTargDiagMatr_sub(Qureg qureg, ConstList64 ctrls, Con
         ctrlStateMask, targ1, targ2,
         elems[0], elems[1], elems[2], elems[3]
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
     // explicitly return to avoid runtime error below
     return;
@@ -729,6 +741,7 @@ void gpu_statevec_anyCtrlAnyTargDiagMatr_sub(Qureg qureg, ConstList64 ctrls, Con
         ctrlStateMask, targs, 
         getGpuQcompPtr(util_getGpuMemPtr(matr)), getGpuQcomp(exponent)
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
     // must return to avoid runtime error below
     return;
@@ -784,6 +797,7 @@ void gpu_densmatr_allTargDiagMatr_sub(Qureg qureg, FullStateDiagMatr matr, qcomp
             getGpuQcompPtr(qureg.gpuAmps), numThreads, qureg.rank, qureg.logNumAmpsPerNode,
             getGpuQcompPtr(util_getGpuMemPtr(matr)), matr.numElems, getGpuQcomp(exponent)
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -846,6 +860,7 @@ void gpu_statevector_anyCtrlPauliTensorOrGadget_subA(Qureg qureg, ConstList64 ct
         targsXY, maskXY, maskYZ,
        	getGpuQcomp(powI), getGpuQcomp(ampFac), getGpuQcomp(pairAmpFac)
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -878,6 +893,7 @@ void gpu_statevector_anyCtrlPauliTensorOrGadget_subB(Qureg qureg, ConstList64 ct
         maskXY, maskYZ, bufferMaskXY,
         getGpuQcomp(powI), getGpuQcomp(ampFac), getGpuQcomp(pairAmpFac)
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -915,6 +931,7 @@ void gpu_statevector_anyCtrlAnyTargZOrPhaseGadget_sub(Qureg qureg, ConstList64 c
         sortedCtrls, ctrlStateMask, targMask,
         getGpuQcomp(fac0), getGpuQcomp(fac1)
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -948,13 +965,16 @@ void gpu_statevec_setQuregToWeightedSum_sub(Qureg outQureg, vector<qcomp> coeffs
     
     // copy coeff and qureg lists into GPU memory, allocating new device memory
     // which will be a visible overhead when the Qureg are small. But eh!
-    devgpuqcompptrs devQuregAmps = ptrs;
-    devcomps devCoeffs = coeffs;
+    devgpuqcompptrs devQuregAmps;
+    devcomps devCoeffs;
+    THRUST_CHECK( devQuregAmps.assign(ptrs.begin(), ptrs.end()) );
+    THRUST_CHECK( devCoeffs.assign(coeffs.begin(), coeffs.end()) );
 
     kernel_statevec_setQuregToWeightedSum_sub <NumQuregs> <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(outQureg.gpuAmps), numThreads,
         getPtr(devCoeffs), getPtr(devQuregAmps), inQuregs.size()
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -986,6 +1006,7 @@ void gpu_densmatr_mixQureg_subB(qreal outProb, Qureg outQureg, qreal inProb, Qur
         outProb, getGpuQcompPtr(outQureg.gpuAmps), inProb, getGpuQcompPtr(inQureg.gpuAmps),
         numThreads, inQureg.numAmps
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1005,6 +1026,7 @@ void gpu_densmatr_mixQureg_subC(qreal outProb, Qureg outQureg, qreal inProb) {
         outProb, getGpuQcompPtr(outQureg.gpuAmps), inProb, getGpuQcompPtr(outQureg.gpuCommBuffer),
         numThreads, outQureg.rank, powerOf2(outQureg.numQubits), outQureg.logNumAmpsPerNode        
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1040,6 +1062,7 @@ void gpu_densmatr_oneQubitDephasing_subA(Qureg qureg, int ketQubit, qreal prob)
     kernel_densmatr_oneQubitDephasing_subA <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(qureg.gpuAmps), numThreads, ketQubit, braQubit, fac
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1067,6 +1090,7 @@ void gpu_densmatr_oneQubitDephasing_subB(Qureg qureg, int ketQubit, qreal prob)
     kernel_densmatr_oneQubitDephasing_subB <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(qureg.gpuAmps), numThreads, ketQubit, braBit, fac
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1115,6 +1139,7 @@ void gpu_densmatr_twoQubitDephasing_subB(Qureg qureg, int ketQubitA, int ketQubi
         getGpuQcompPtr(qureg.gpuAmps), numThreads, qureg.rank, qureg.logNumAmpsPerNode, // numAmps, not numCols
         ketQubitA, ketQubitB, braQubitA, braQubitB, term
     );
+    assert_lastKernelLaunchSucceeded(__func__);
     
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1142,6 +1167,7 @@ void gpu_densmatr_oneQubitDepolarising_subA(Qureg qureg, int ketQubit, qreal pro
     kernel_densmatr_oneQubitDepolarising_subA <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(qureg.gpuAmps), numThreads, ketQubit, braQubit, factors.c1, factors.c2, factors.c3
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1165,6 +1191,7 @@ void gpu_densmatr_oneQubitDepolarising_subB(Qureg qureg, int ketQubit, qreal pro
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + recvInd, numThreads, 
         ketQubit, braBit, factors.c1, factors.c2, factors.c3
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1194,6 +1221,7 @@ void gpu_densmatr_twoQubitDepolarising_subA(Qureg qureg, int ketQb1, int ketQb2,
         getGpuQcompPtr(qureg.gpuAmps), numThreads,
         ketQb1, ketQb2, braQb1, braQb2, c3
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1220,6 +1248,7 @@ void gpu_densmatr_twoQubitDepolarising_subB(Qureg qureg, int ketQb1, int ketQb2,
         getGpuQcompPtr(qureg.gpuAmps), numThreads,
         ketQb1, ketQb2, braQb1, braQb2, altc1, factors.c2
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1243,6 +1272,7 @@ void gpu_densmatr_twoQubitDepolarising_subC(Qureg qureg, int ketQb1, int ketQb2,
         getGpuQcompPtr(qureg.gpuAmps), numThreads,
         ketQb1, ketQb2, braQb1, braBit2, c3
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1267,6 +1297,7 @@ void gpu_densmatr_twoQubitDepolarising_subD(Qureg qureg, int ketQb1, int ketQb2,
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + offset, numThreads,
         ketQb1, ketQb2, braQb1, braBit2, factors.c1, factors.c2
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1293,6 +1324,7 @@ void gpu_densmatr_twoQubitDepolarising_subE(Qureg qureg, int ketQb1, int ketQb2,
         getGpuQcompPtr(qureg.gpuAmps), numThreads,
         ketQb1, ketQb2, braBit1, braBit2, fac0, fac1
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1317,6 +1349,7 @@ void gpu_densmatr_twoQubitDepolarising_subF(Qureg qureg, int ketQb1, int ketQb2,
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + offset, numThreads,
         ketQb1, ketQb2, braBit1, braBit2, c2
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1345,6 +1378,7 @@ void gpu_densmatr_oneQubitPauliChannel_subA(Qureg qureg, int ketQubit, qreal pI,
         getGpuQcompPtr(qureg.gpuAmps), numThreads, ketQubit, braQubit, 
         factors.c1, factors.c2, factors.c3, factors.c4
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1368,6 +1402,7 @@ void gpu_densmatr_oneQubitPauliChannel_subB(Qureg qureg, int ketQubit, qreal pI,
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + recvInd, numThreads, 
         ketQubit, braBit, factors.c1, factors.c2, factors.c3, factors.c4
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1396,6 +1431,7 @@ void gpu_densmatr_oneQubitDamping_subA(Qureg qureg, int ketQubit, qreal prob) {
         getGpuQcompPtr(qureg.gpuAmps), numThreads,
         ketQubit, braQubit, prob, factors.c1, factors.c2
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1416,6 +1452,7 @@ void gpu_densmatr_oneQubitDamping_subB(Qureg qureg, int qubit, qreal prob) {
     kernel_densmatr_oneQubitDamping_subB <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(qureg.gpuAmps), numThreads, qubit, c2
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1437,6 +1474,7 @@ void gpu_densmatr_oneQubitDamping_subC(Qureg qureg, int ketQubit, qreal prob) {
     kernel_densmatr_oneQubitDamping_subC <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(qureg.gpuAmps), numThreads, ketQubit, braBit, c1
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1457,6 +1495,7 @@ void gpu_densmatr_oneQubitDamping_subD(Qureg qureg, int qubit, qreal prob) {
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + recvInd, numThreads, 
         qubit, prob
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1487,6 +1526,7 @@ void gpu_densmatr_partialTrace_sub(Qureg inQureg, Qureg outQureg, ConstList64 ta
         getGpuQcompPtr(inQureg.gpuAmps), getGpuQcompPtr(outQureg.gpuAmps), numThreads,
         targs, pairTargs, allTargs
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -1601,12 +1641,13 @@ void gpu_statevec_calcProbsOfAllMultiQubitOutcomes_sub(qreal* outProbs, Qureg qu
     qindex numBlocks = getNumBlocks(numThreads, numThreadsPerBlock);
 
     // allocate exponentially-big temporary memory (error if failed)
-    devreals devProbs = getDeviceRealsVec(powerOf2(qubits.size())); // throws
+    devreals devProbs = getDeviceRealsVec(powerOf2(qubits.size()));
 
     kernel_statevec_calcProbsOfAllMultiQubitOutcomes_sub<NumQubits> <<<numBlocks, numThreadsPerBlock>>> (
        getPtr(devProbs), getGpuQcompPtr(qureg.gpuAmps), numThreads, 
         qureg.rank, qureg.logNumAmpsPerNode, qubits
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
     // overwrite outProbs with GPU memory
     copyFromDeviceVec(devProbs, outProbs);
@@ -1638,7 +1679,7 @@ void gpu_densmatr_calcProbsOfAllMultiQubitOutcomes_sub(qreal* outProbs, Qureg qu
     qindex numAmpsPerCol = powerOf2(qureg.numQubits);
 
     // allocate exponentially-big temporary memory (error if failed)
-    devreals devProbs = getDeviceRealsVec(powerOf2(qubits.size())); // throws
+    devreals devProbs = getDeviceRealsVec(powerOf2(qubits.size()));
 
     kernel_densmatr_calcProbsOfAllMultiQubitOutcomes_sub<NumQubits> <<<numBlocks, numThreadsPerBlock>>> (
         getPtr(devProbs), getGpuQcompPtr(qureg.gpuAmps), 
@@ -1646,6 +1687,7 @@ void gpu_densmatr_calcProbsOfAllMultiQubitOutcomes_sub(qreal* outProbs, Qureg qu
         qureg.rank, qureg.logNumAmpsPerNode, 
         qubits
     );
+    assert_lastKernelLaunchSucceeded(__func__);
 
     // overwrite outProbs with GPU memory
     copyFromDeviceVec(devProbs, outProbs);
diff --git a/quest/src/gpu/gpu_thrust.cuh b/quest/src/gpu/gpu_thrust.cuh
index ae226545c..63974a90a 100644
--- a/quest/src/gpu/gpu_thrust.cuh
+++ b/quest/src/gpu/gpu_thrust.cuh
@@ -58,6 +58,25 @@
 #include <thrust/iterator/transform_iterator.h>
 #include <thrust/system/system_error.h>
 
+#include <exception>
+
+
+
+/*
+ * THRUST ERROR HANDLING
+ */
+
+#define THRUST_CHECK(cmd)                                                                                       \
+    do {                                                                                                        \
+        try {                                                                                                   \
+            cmd;                                                                                                \
+        } catch (const std::exception& err) {                                                                   \
+            error_thrustCallFailed(err.what(), #cmd, __func__, __FILE__, __LINE__);                             \
+        } catch (...) {                                                                                         \
+            error_thrustCallFailed("Unknown non-standard exception", #cmd, __func__, __FILE__, __LINE__);       \
+        }                                                                                                       \
+    } while (0)
+
 
 
 /*
@@ -82,20 +101,15 @@ qreal* getPtr(devreals& reals) {
 
 void copyFromDeviceVec(devreals& reals, qreal* out) {
 
-    thrust::copy(reals.begin(), reals.end(), out);
+    THRUST_CHECK( thrust::copy(reals.begin(), reals.end(), out) );
 }
 
 devreals getDeviceRealsVec(qindex dim) {
 
     devreals out;
 
-    try  {
-        out.resize(dim);
-        thrust::fill(out.begin(), out.end(), 0.);
-
-    } catch (thrust::system_error &e) {
-        error_thrustTempGpuAllocFailed();
-    }
+    THRUST_CHECK( out.resize(dim) );
+    THRUST_CHECK( thrust::fill(out.begin(), out.end(), 0.) );
 
     return out;
 }
@@ -626,8 +640,10 @@ void thrust_fullstatediagmatr_setElemsToPauliStrSum(FullStateDiagMatr out, Pauli
 
     // copy 'in' lists into GPU memory, which is not a big deal even when 'in'
     // is very large, because we only do this during FullStateDiagMatr initialisation
-    thrust::device_vector<qcomp> devCoeffs(in.coeffs, in.coeffs + in.numTerms);
-    thrust::device_vector<PauliStr> devStrings(in.strings, in.strings + in.numTerms);
+    thrust::device_vector<qcomp> devCoeffs;
+    thrust::device_vector<PauliStr> devStrings;
+    THRUST_CHECK( devCoeffs.assign(in.coeffs, in.coeffs + in.numTerms) );
+    THRUST_CHECK( devStrings.assign(in.strings, in.strings + in.numTerms) );
     
     // obtain raw pointers which can be passed to fastmath.hpp routines
     gpu_qcomp* devCoeffsPtr = getGpuQcompPtr(thrust::raw_pointer_cast(devCoeffs.data()));
@@ -643,7 +659,7 @@ void thrust_fullstatediagmatr_setElemsToPauliStrSum(FullStateDiagMatr out, Pauli
 
     auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + out.numElemsPerNode;
-    thrust::for_each(indIter, endIter, functor);
+    THRUST_CHECK( thrust::for_each(indIter, endIter, functor) );
 }
 
 
@@ -656,7 +672,7 @@ void thrust_fullstatediagmatr_setElemsToPauliStrSum(FullStateDiagMatr out, Pauli
 void thrust_setElemsToConjugate(gpu_qcomp* matrElemsPtr, qindex matrElemsLen) {
 
     auto ptr = getStartPtr(matrElemsPtr);
-    thrust::transform(ptr, ptr + matrElemsLen, ptr, functor_getAmpConj());
+    THRUST_CHECK( thrust::transform(ptr, ptr + matrElemsLen, ptr, functor_getAmpConj()) );
 }
 
 
@@ -672,8 +688,10 @@ void thrust_densmatr_setAmpsToPauliStrSum_sub(Qureg qureg, PauliStrSum sum) {
 
     // copy sum lists into GPU memory, which is not a big deal even when sum
     // is very large, because we only do this during Qureg initialisation (infrequent)
-    thrust::device_vector<qcomp> devCoeffs(sum.coeffs, sum.coeffs + sum.numTerms);
-    thrust::device_vector<PauliStr> devStrings(sum.strings, sum.strings + sum.numTerms);
+    thrust::device_vector<qcomp> devCoeffs;
+    thrust::device_vector<PauliStr> devStrings;
+    THRUST_CHECK( devCoeffs.assign(sum.coeffs, sum.coeffs + sum.numTerms) );
+    THRUST_CHECK( devStrings.assign(sum.strings, sum.strings + sum.numTerms) );
     
     // obtain raw pointers which can be passed to fastmath.hpp routines
     gpu_qcomp* devCoeffsPtr = getGpuQcompPtr(thrust::raw_pointer_cast(devCoeffs.data()));
@@ -686,26 +704,26 @@ void thrust_densmatr_setAmpsToPauliStrSum_sub(Qureg qureg, PauliStrSum sum) {
 
     auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + qureg.numAmpsPerNode;
-    thrust::for_each(indIter, endIter, functor);
+    THRUST_CHECK( thrust::for_each(indIter, endIter, functor) );
 }
 
 
 void thrust_densmatr_mixQureg_subA(qreal outProb, Qureg outQureg, qreal inProb, Qureg inQureg) {
 
-    thrust::transform(
+    THRUST_CHECK( thrust::transform(
         getStartPtr(outQureg), getEndPtr(outQureg), 
         getStartPtr(inQureg),  getStartPtr(outQureg), // 4th arg is output pointer
-        functor_mixAmps(outProb, inProb));
+        functor_mixAmps(outProb, inProb)) );
 }
 
 
 template <bool HasPower>
 void thrust_statevec_allTargDiagMatr_sub(Qureg qureg, FullStateDiagMatr matr, gpu_qcomp exponent) {
 
-    thrust::transform(
+    THRUST_CHECK( thrust::transform(
         getStartPtr(qureg), getEndPtr(qureg), 
         getStartPtr(matr),  getStartPtr(qureg), // 4th arg is output pointer
-        functor_multiplyElemPowerWithAmpOrNorm<HasPower,false,false>(exponent));
+        functor_multiplyElemPowerWithAmpOrNorm<HasPower,false,false>(exponent)) );
 }
 
 
@@ -729,9 +747,10 @@ qreal thrust_statevec_calcTotalProb_sub(Qureg qureg) {
     // literal causes a silent Thrust error. Grr...
     qreal init = 0.0;
 
-    qreal prob = thrust::transform_reduce(
+    qreal prob = 0;
+    THRUST_CHECK( prob = thrust::transform_reduce(
         getStartPtr(qureg), getEndPtr(qureg), 
-        functor_getAmpNorm(), init, thrust::plus<qreal>());
+        functor_getAmpNorm(), init, thrust::plus<qreal>()) );
 
     return prob;
 }
@@ -753,7 +772,8 @@ qreal thrust_densmatr_calcTotalProb_sub(Qureg qureg) {
     auto probIter= thrust::make_transform_iterator(ampIter, functor_getAmpReal());
 
     qindex numIts = powerOf2(qureg.logNumColsPerNode);
-    qreal prob = thrust::reduce(probIter, probIter + numIts);
+    qreal prob = 0;
+    THRUST_CHECK( prob = thrust::reduce(probIter, probIter + numIts) );
     return prob;
 }
 
@@ -774,7 +794,8 @@ qreal thrust_statevec_calcProbOfMultiQubitOutcome_sub(Qureg qureg, ConstList64 q
     auto probIter = thrust::make_transform_iterator(ampIter, probFunctor);
 
     qindex numIts = qureg.numAmpsPerNode / powerOf2(qubits.size());
-    qreal prob = thrust::reduce(probIter, probIter + numIts);
+    qreal prob = 0;
+    THRUST_CHECK( prob = thrust::reduce(probIter, probIter + numIts) );
     return prob;
 }
 
@@ -797,7 +818,8 @@ qreal thrust_densmatr_calcProbOfMultiQubitOutcome_sub(Qureg qureg, ConstList64 q
     auto probIter= thrust::make_transform_iterator(ampIter, probFunctor);
 
     qindex numIts = powerOf2(qureg.logNumColsPerNode - qubits.size());
-    qreal prob = thrust::reduce(probIter, probIter + numIts);
+    qreal prob = 0;
+    THRUST_CHECK( prob = thrust::reduce(probIter, probIter + numIts) );
     return prob;
 }
 
@@ -812,9 +834,10 @@ gpu_qcomp thrust_statevec_calcInnerProduct_sub(Qureg quregA, Qureg quregB) {
 
     gpu_qcomp init = getGpuQcomp(0, 0);
 
-    gpu_qcomp prod = thrust::inner_product(
+    gpu_qcomp prod = getGpuQcomp(0, 0);
+    THRUST_CHECK( prod = thrust::inner_product(
         getStartPtr(quregA), getEndPtr(quregA), getStartPtr(quregB), 
-        init, thrust::plus<gpu_qcomp>(), functor_getAmpConjProd());
+        init, thrust::plus<gpu_qcomp>(), functor_getAmpConjProd()) );
 
     return prod;
 }
@@ -824,9 +847,10 @@ qreal thrust_densmatr_calcHilbertSchmidtDistance_sub(Qureg quregA, Qureg quregB)
 
     qreal init = 0;
 
-    qreal dist = thrust::inner_product(
+    qreal dist = 0;
+    THRUST_CHECK( dist = thrust::inner_product(
         getStartPtr(quregA), getEndPtr(quregA), getStartPtr(quregB), 
-        init, thrust::plus<qreal>(), functor_getNormOfAmpDif());
+        init, thrust::plus<qreal>(), functor_getNormOfAmpDif()) );
 
     return dist;
 }
@@ -844,9 +868,10 @@ gpu_qcomp thrust_densmatr_calcFidelityWithPureState_sub(Qureg rho, Qureg psi) {
     qindex numIts = rho.numAmpsPerNode;
 
     gpu_qcomp init = getGpuQcomp(0, 0);
-    gpu_qcomp fid = thrust::transform_reduce(
+    gpu_qcomp fid = getGpuQcomp(0, 0);
+    THRUST_CHECK( fid = thrust::transform_reduce(
         indIter, indIter + numIts, 
-        functor, init, thrust::plus<gpu_qcomp>());
+        functor, init, thrust::plus<gpu_qcomp>()) );
 
     return fid;
 }
@@ -867,9 +892,11 @@ qreal thrust_statevec_calcExpecAnyTargZ_sub(Qureg qureg, ConstList64 targs) {
     auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + qureg.numAmpsPerNode;
 
-    return thrust::inner_product(
+    qreal value = 0;
+    THRUST_CHECK( value = thrust::inner_product(
         indIter, endIter, getStartPtr(qureg), 
-        init, thrust::plus<qreal>(), functor);
+        init, thrust::plus<qreal>(), functor) );
+    return value;
 }
 
 
@@ -884,7 +911,9 @@ gpu_qcomp thrust_densmatr_calcExpecAnyTargZ_sub(Qureg qureg, ConstList64 targs)
     auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + powerOf2(qureg.logNumColsPerNode);
 
-    return thrust::transform_reduce(indIter, endIter, functor, init, thrust::plus<gpu_qcomp>());
+    gpu_qcomp value = getGpuQcomp(0, 0);
+    THRUST_CHECK( value = thrust::transform_reduce(indIter, endIter, functor, init, thrust::plus<gpu_qcomp>()) );
+    return value;
 }
 
 
@@ -899,7 +928,8 @@ gpu_qcomp thrust_statevec_calcExpecPauliStr_subA(Qureg qureg, ConstList64 x, Con
     auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + qureg.numAmpsPerNode;
 
-    gpu_qcomp value = thrust::transform_reduce(indIter, endIter, functor, init, thrust::plus<gpu_qcomp>());
+    gpu_qcomp value = getGpuQcomp(0, 0);
+    THRUST_CHECK( value = thrust::transform_reduce(indIter, endIter, functor, init, thrust::plus<gpu_qcomp>()) );
 
     return value * getGpuQcomp(util_getPowerOfI(y.size()));
 }
@@ -917,7 +947,8 @@ gpu_qcomp thrust_statevec_calcExpecPauliStr_subB(Qureg qureg, ConstList64 x, Con
     auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + qureg.numAmpsPerNode;
 
-    gpu_qcomp value = thrust::transform_reduce(indIter, endIter, functor, init, thrust::plus<gpu_qcomp>());
+    gpu_qcomp value = getGpuQcomp(0, 0);
+    THRUST_CHECK( value = thrust::transform_reduce(indIter, endIter, functor, init, thrust::plus<gpu_qcomp>()) );
 
     return value * getGpuQcomp(util_getPowerOfI(y.size()));
 }
@@ -935,7 +966,8 @@ gpu_qcomp thrust_densmatr_calcExpecPauliStr_sub(Qureg qureg, ConstList64 x, Cons
     auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + powerOf2(qureg.logNumColsPerNode);
 
-    gpu_qcomp value = thrust::transform_reduce(indIter, endIter, functor, init, thrust::plus<gpu_qcomp>());
+    gpu_qcomp value = getGpuQcomp(0, 0);
+    THRUST_CHECK( value = thrust::transform_reduce(indIter, endIter, functor, init, thrust::plus<gpu_qcomp>()) );
 
     return value * getGpuQcomp(util_getPowerOfI(y.size()));
 }
@@ -953,9 +985,10 @@ gpu_qcomp thrust_statevec_calcExpecFullStateDiagMatr_sub(Qureg qureg, FullStateD
     gpu_qcomp init = getGpuQcomp(0, 0);
     auto functor = functor_multiplyElemPowerWithAmpOrNorm<HasPower,UseRealPow,true>(expo);
 
-    gpu_qcomp value = thrust::inner_product(
+    gpu_qcomp value = getGpuQcomp(0, 0);
+    THRUST_CHECK( value = thrust::inner_product(
         getStartPtr(qureg), getEndPtr(qureg), getStartPtr(matr), 
-        init, thrust::plus<gpu_qcomp>(), functor);
+        init, thrust::plus<gpu_qcomp>(), functor) );
 
     return value;
 }
@@ -974,7 +1007,9 @@ gpu_qcomp thrust_densmatr_calcExpecFullStateDiagMatr_sub(Qureg qureg, FullStateD
     auto indIter = thrust::make_counting_iterator(QINDEX_ZERO);
     auto endIter = indIter + powerOf2(qureg.logNumColsPerNode);
 
-    return thrust::transform_reduce(indIter, endIter, functor, init, thrust::plus<gpu_qcomp>());
+    gpu_qcomp value = getGpuQcomp(0, 0);
+    THRUST_CHECK( value = thrust::transform_reduce(indIter, endIter, functor, init, thrust::plus<gpu_qcomp>()) );
+    return value;
 }
 
 
@@ -995,7 +1030,7 @@ void thrust_statevec_multiQubitProjector_sub(Qureg qureg, ConstList64 qubits, Co
     auto ampIter = getStartPtr(qureg);
 
     qindex numIts = qureg.numAmpsPerNode;
-    thrust::transform(indIter, indIter + numIts, ampIter, ampIter, projFunctor); // 4th arg gets modified
+    THRUST_CHECK( thrust::transform(indIter, indIter + numIts, ampIter, ampIter, projFunctor) ); // 4th arg gets modified
 }
 
 
@@ -1011,7 +1046,7 @@ void thrust_densmatr_multiQubitProjector_sub(Qureg qureg, ConstList64 qubits, Co
     auto ampIter = getStartPtr(qureg);
 
     qindex numIts = qureg.numAmpsPerNode;
-    thrust::transform(indIter, indIter + numIts, ampIter, ampIter, projFunctor); // 4th arg gets modified
+    THRUST_CHECK( thrust::transform(indIter, indIter + numIts, ampIter, ampIter, projFunctor) ); // 4th arg gets modified
 }
 
 
@@ -1023,7 +1058,7 @@ void thrust_densmatr_multiQubitProjector_sub(Qureg qureg, ConstList64 qubits, Co
 
 void thrust_statevec_initUniformState(Qureg qureg, gpu_qcomp amp) {
 
-    thrust::fill(getStartPtr(qureg), getEndPtr(qureg), amp);
+    THRUST_CHECK( thrust::fill(getStartPtr(qureg), getEndPtr(qureg), amp) );
 }
 
 
@@ -1037,7 +1072,7 @@ void thrust_statevec_initDebugState_sub(Qureg qureg) {
     qindex n = util_getGlobalIndexOfFirstLocalAmp(qureg);
     gpu_qcomp init = getGpuQcomp(2*n/10., (2*n+1)/10.);
 
-    thrust::sequence(getStartPtr(qureg), getEndPtr(qureg), init, step);
+    THRUST_CHECK( thrust::sequence(getStartPtr(qureg), getEndPtr(qureg), init, step) );
 }
 
 
@@ -1051,7 +1086,7 @@ void thrust_statevec_initUnnormalisedUniformlyRandomPureStateAmps_sub(Qureg qure
     auto ampIter = getStartPtr(qureg);
 
     qindex numIts = qureg.numAmpsPerNode;
-    thrust::transform(indIter, indIter + numIts, ampIter, functor); // 3rd arg gets modified
+    THRUST_CHECK( thrust::transform(indIter, indIter + numIts, ampIter, functor) ); // 3rd arg gets modified
 }