diff --git a/quest/src/core/bitwise.hpp b/quest/src/core/bitwise.hpp
index f5266afa4..17f384af3 100644
--- a/quest/src/core/bitwise.hpp
+++ b/quest/src/core/bitwise.hpp
@@ -164,23 +164,22 @@ INLINE int getBitMaskParity(qindex mask) {
 
 
 INLINE qindex insertBits(qindex number, const int* bitIndices, int numIndices, int bitValue) {
-    
+
     // bitIndices must be strictly increasing
     for (int i=0; i<numIndices; i++)
         number = insertBit(number, bitIndices[i], bitValue);
-        
+
     return number;
 }
 
-
 INLINE qindex setBits(qindex number, const int* bitIndices, int numIndices, qindex bitsValue) {
-    
+
     // bitIndices are arbitrarily ordered, which does not affect number
     for (int i=0; i<numIndices; i++) {
         int bit = getBit(bitsValue, i);
         number = setBit(number, bitIndices[i], bit);
     }
-    
+
     return number;
 }
 
@@ -196,8 +195,6 @@ INLINE qindex getValueOfBits(qindex number, const int* bitIndices, int numIndice
     return value;
 }
 
-
-
 /*
  * PERFORMANCE-CRITICAL CONVENIENCE FUNCTIONS
  *
@@ -211,7 +208,6 @@ INLINE qindex insertBitsWithMaskedValues(qindex number, const int* bitInds, int
     return mask | insertBits(number, bitInds, numBits, 0);
 }
 
-
 INLINE int getTwoBits(qindex number, int highInd, int lowInd) {
 
     int b1 = getBit(number, lowInd);
@@ -379,4 +375,4 @@ INLINE void setToBitsOfInteger(int* bits, qindex number, int numBits) {
 
 
 
-#endif // BITWISE_HPP
\ No newline at end of file
+#endif // BITWISE_HPP
diff --git a/quest/src/gpu/gpu_kernels.cuh b/quest/src/gpu/gpu_kernels.cuh
index b6954f701..314847430 100644
--- a/quest/src/gpu/gpu_kernels.cuh
+++ b/quest/src/gpu/gpu_kernels.cuh
@@ -86,15 +86,14 @@ __forceinline__ __device__ int cudaGetBitMaskParity(qindex mask) {
 template <int NumCtrls>
 __global__ void kernel_statevec_packAmpsIntoBuffer(
     gpu_qcomp* amps, gpu_qcomp* buffer, qindex numThreads, 
-    int* qubits, int numQubits, qindex qubitStateMask
-) {
+    __grid_constant__ const List64 qubits, qindex qubitStateMask) {
     GET_THREAD_IND(n, numThreads);
 
     // use template param to compile-time unroll loop in insertBits()
-    SET_VAR_AT_COMPILE_TIME(int, numBits, NumCtrls, numQubits);
+    SET_VAR_AT_COMPILE_TIME(int, numBits, NumCtrls, qubits.size());
 
     // i = nth local index where qubits are active
-    qindex i = insertBitsWithMaskedValues(n, qubits, numBits, qubitStateMask);
+    qindex i = insertBitsWithMaskedValues(n, qubits.data(), numBits, qubitStateMask);
 
     // caller offsets buffer by sub-buffer send-index
     buffer[n] = amps[i];
@@ -103,8 +102,7 @@ __global__ void kernel_statevec_packAmpsIntoBuffer(
 
 __global__ void kernel_statevec_packPairSummedAmpsIntoBuffer(
     gpu_qcomp* amps, gpu_qcomp* buffer, qindex numThreads, 
-    int qubit1, int qubit2, int qubit3, int bit2
-) {
+    int qubit1, int qubit2, int qubit3, int bit2) {
     GET_THREAD_IND(n, numThreads);
 
     // i000 = nth local index where all qubits are 0
@@ -125,16 +123,15 @@ __global__ void kernel_statevec_packPairSummedAmpsIntoBuffer(
 template <int NumCtrls> 
 __global__ void kernel_statevec_anyCtrlSwap_subA(
     gpu_qcomp* amps, qindex numThreads, 
-    int* ctrlsAndTargs, int numCtrls, qindex ctrlsAndTargsMask, int targ1, int targ2
-) {
+    __grid_constant__ const List64 ctrlsAndTargs, qindex ctrlsAndTargsMask, int targ1, int targ2) {
     GET_THREAD_IND(n, numThreads);
 
     // use template param to compile-time unroll loop in insertBits()
-    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
+    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrlsAndTargs.size());
     int numQubitBits = 2 + numCtrlBits;
 
     // i01 = nth local index where ctrls are active, targ2=0 and targ1=1
-    qindex i01 = insertBitsWithMaskedValues(n, ctrlsAndTargs, numQubitBits, ctrlsAndTargsMask);
+    qindex i01 = insertBitsWithMaskedValues(n, ctrlsAndTargs.data(), numQubitBits, ctrlsAndTargsMask);
     qindex i10 = flipTwoBits(i01, targ2, targ1);
 
     // swap amps
@@ -147,15 +144,15 @@ __global__ void kernel_statevec_anyCtrlSwap_subA(
 template <int NumCtrls> 
 __global__ void kernel_statevec_anyCtrlSwap_subB(
     gpu_qcomp* amps, gpu_qcomp* buffer, qindex numThreads, 
-    int* ctrls, int numCtrls, qindex ctrlStateMask
+    __grid_constant__ const List64 ctrls, qindex ctrlStateMask
 ) {
     GET_THREAD_IND(n, numThreads);
 
     // use template param to compile-time unroll loop in insertBits()
-    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
+    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrls.size());
 
     // i = nth local index where ctrls are active
-    qindex i = insertBitsWithMaskedValues(n, ctrls, numCtrlBits, ctrlStateMask);
+    qindex i = insertBitsWithMaskedValues(n, ctrls.data(), numCtrlBits, ctrlStateMask);
 
     // caller offsets buffer if necessary
     amps[i] = buffer[n];
@@ -165,16 +162,16 @@ __global__ void kernel_statevec_anyCtrlSwap_subB(
 template <int NumCtrls> 
 __global__ void kernel_statevec_anyCtrlSwap_subC(
     gpu_qcomp* amps, gpu_qcomp* buffer, qindex numThreads, 
-    int* ctrlsAndTarg, int numCtrls, qindex ctrlsAndTargMask
+    __grid_constant__ const List64 ctrlsAndTarg, qindex ctrlsAndTargMask
 ) {
     GET_THREAD_IND(n, numThreads);
 
     // use template param to compile-time unroll loop in insertBits()
-    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
+    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrlsAndTarg.size());
     int numQubitBits = numCtrlBits + 1;
 
     // i = nth local index where ctrls and targ are in specified states
-    qindex i = insertBitsWithMaskedValues(n, ctrlsAndTarg, numQubitBits, ctrlsAndTargMask);
+    qindex i = insertBitsWithMaskedValues(n, ctrlsAndTarg.data(), numQubitBits, ctrlsAndTargMask);
 
     // caller offsets buffer if necessary
     amps[i] = buffer[n];
@@ -189,17 +186,17 @@ __global__ void kernel_statevec_anyCtrlSwap_subC(
 
 template <int NumCtrls>
 __global__ void kernel_statevec_anyCtrlOneTargDenseMatr_subA(
-    gpu_qcomp* amps, qindex numThreads, 
-    int* ctrlsAndTarg, int numCtrls, qindex ctrlStateMask, int targ, 
+    gpu_qcomp* amps, qindex numThreads, __grid_constant__ const List64 ctrl,
+    qindex ctrlStateMask, int targ, 
     gpu_qcomp m00, gpu_qcomp m01, gpu_qcomp m10, gpu_qcomp m11
 ) {
     GET_THREAD_IND(n, numThreads);
 
     // use template param to compile-time unroll loop in insertBits()
-    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
+    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrl.size());
 
     // i0 = nth local index where ctrls are active and targ is 0
-    qindex i0 = insertBitsWithMaskedValues(n, ctrlsAndTarg, numCtrlBits + 1, ctrlStateMask);
+    qindex i0 = insertBitsWithMaskedValues(n, ctrl.data(), numCtrlBits + 1, ctrlStateMask);
     qindex i1 = flipBit(i0, targ);
 
     // note amps are strided by 2^targ
@@ -214,16 +211,16 @@ __global__ void kernel_statevec_anyCtrlOneTargDenseMatr_subA(
 template <int NumCtrls>
 __global__ void kernel_statevec_anyCtrlOneTargDenseMatr_subB(
     gpu_qcomp* amps, gpu_qcomp* buffer, qindex numThreads, 
-    int* ctrls, int numCtrls, qindex ctrlStateMask,
+    __grid_constant__ const List64 ctrls, qindex ctrlStateMask,
     gpu_qcomp fac0, gpu_qcomp fac1
 ) {
     GET_THREAD_IND(n, numThreads);
 
     // use template param to compile-time unroll loop in insertBits()
-    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
+    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrls.size());
 
     // i = nth local index where ctrl bits are active
-    qindex i = insertBitsWithMaskedValues(n, ctrls, numCtrlBits, ctrlStateMask);
+    qindex i = insertBitsWithMaskedValues(n, ctrls.data(), numCtrlBits, ctrlStateMask);
 
     // caller offsets buffer by receive-index
     amps[i] = fac0*amps[i] + fac1*buffer[n];
@@ -239,7 +236,7 @@ __global__ void kernel_statevec_anyCtrlOneTargDenseMatr_subB(
 template <int NumCtrls>
 __global__ void kernel_statevec_anyCtrlTwoTargDenseMatr_sub(
     gpu_qcomp* amps, qindex numThreads, 
-    int* ctrlsAndTarg, int numCtrls, qindex ctrlStateMask, int targ1, int targ2,
+    __grid_constant__ const List64 ctrlsAndTarg, qindex ctrlStateMask, int targ1, int targ2,
     gpu_qcomp m00, gpu_qcomp m01, gpu_qcomp m02, gpu_qcomp m03,
     gpu_qcomp m10, gpu_qcomp m11, gpu_qcomp m12, gpu_qcomp m13,
     gpu_qcomp m20, gpu_qcomp m21, gpu_qcomp m22, gpu_qcomp m23,
@@ -248,10 +245,10 @@ __global__ void kernel_statevec_anyCtrlTwoTargDenseMatr_sub(
     GET_THREAD_IND(n, numThreads);
 
     // use template param to compile-time unroll loop in insertBits()
-    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
+    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrlsAndTarg.size());
 
     // i00 = nth local index where ctrls are active and both targs are 0
-    qindex i00 = insertBitsWithMaskedValues(n, ctrlsAndTarg, numCtrlBits + 2, ctrlStateMask);
+    qindex i00 = insertBitsWithMaskedValues(n, ctrlsAndTarg.data(), numCtrlBits + 2, ctrlStateMask);
     qindex i01 = flipBit(i00, targ1);
     qindex i10 = flipBit(i00, targ2);
     qindex i11 = flipBit(i01, targ2);
@@ -288,7 +285,7 @@ __forceinline__ __device__ qindex getThreadsNthGlobalArrInd(qindex n, qindex thr
 template <int NumCtrls, int NumTargs, bool ApplyConj, bool ApplyTransp>
 __global__ void kernel_statevec_anyCtrlFewTargDenseMatr(
     gpu_qcomp* amps, qindex numThreads,
-    int* ctrlsAndTargs, int numCtrls, qindex ctrlsAndTargsMask, int* targs,
+    __grid_constant__ const List64 ctrlsAndTargs, qindex ctrlsAndTargsMask, __grid_constant__ const List64 targs,
     gpu_qcomp* flatMatrElems
 ) {
     GET_THREAD_IND(n, numThreads);
@@ -309,18 +306,18 @@ __global__ void kernel_statevec_anyCtrlFewTargDenseMatr(
     REGISTER gpu_qcomp privateCache[1 << NumTargs];
 
     // we know NumTargs <= 5, though NumCtrls is permitted anything (including -1)
-    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
+    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrlsAndTargs.size());
     constexpr qindex numTargAmps = (1 << NumTargs); // explicit, in lieu of powerOf2
 
     // i0 = nth local index where ctrls are active and targs are all zero
-    qindex i0 = insertBitsWithMaskedValues(n, ctrlsAndTargs, numCtrlBits + NumTargs, ctrlsAndTargsMask); // loop may be unrolled
+    qindex i0 = insertBitsWithMaskedValues(n, ctrlsAndTargs.data(), numCtrlBits + NumTargs, ctrlsAndTargsMask); // loop may be unrolled
 
     // populate cache (force unroll to ensure compile-time cache indices)
     #pragma unroll  
     for (qindex k=0; k<numTargAmps; k++) {
 
         // i = nth local index where ctrls are active and targs form value k
-        qindex i = setBits(i0, targs, NumTargs, k); // loop will be unrolled
+        qindex i = setBits(i0, targs.data(), NumTargs, k); // loop will be unrolled
 
         // write to thread-private cache at compile-time known index
         privateCache[k] = amps[i];
@@ -330,7 +327,7 @@ __global__ void kernel_statevec_anyCtrlFewTargDenseMatr(
     for (qindex k=0; k<numTargAmps; k++) {
 
         // i = nth local index where ctrls are active and targs form value k
-        qindex i = setBits(i0, targs, NumTargs, k); // loop will be unrolled
+        qindex i = setBits(i0, targs.data(), NumTargs, k); // loop will be unrolled
         amps[i] = getGpuQcomp(0, 0);
     
         // force unroll to ensure compile-time cache indices
@@ -360,14 +357,13 @@ template <int NumCtrls, bool ApplyConj, bool ApplyTransp>
 __global__ void kernel_statevec_anyCtrlManyTargDenseMatr(
     gpu_qcomp* globalCache,
     gpu_qcomp* amps, qindex numThreads, qindex numBatchesPerThread,
-    int* ctrlsAndTargs, int numCtrls, qindex ctrlsAndTargsMask, 
-    int* targs, int numTargBits, qindex numTargAmps,
-    gpu_qcomp* flatMatrElems
-) {
+    __grid_constant__ const List64 ctrlsAndTargs, qindex ctrlsAndTargsMask, 
+    __grid_constant__ const List64 targs, qindex numTargAmps,
+    gpu_qcomp* flatMatrElems) {
     GET_THREAD_IND(t, numThreads);
 
     // NumCtrls might be compile-time known, but numTargBits>5 is always unknown/runtime
-    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
+    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrlsAndTargs.size());
 
     // unlike all other kernels, each thread modifies multiple batches of amplitudes
     for (qindex b=0; b<numBatchesPerThread; b++) {
@@ -376,13 +372,13 @@ __global__ void kernel_statevec_anyCtrlManyTargDenseMatr(
         qindex n = t + b * numThreads;
 
         // i0 = nth local index where ctrls are active and targs are all zero
-        qindex i0 = insertBitsWithMaskedValues(n, ctrlsAndTargs, numCtrlBits + numTargBits, ctrlsAndTargsMask);
+        qindex i0 = insertBitsWithMaskedValues(n, ctrlsAndTargs.data(), numCtrlBits + targs.size(), ctrlsAndTargsMask);
 
         // collect and cache all to-be-modified amps (loop might be unrolled)        
         for (qindex k=0; k<numTargAmps; k++) {
 
             // i = nth local index where ctrls are active and targs form value k
-            qindex i = setBits(i0, targs, numTargBits, k); // loop may be unrolled
+            qindex i = setBits(i0, targs.data(), targs.size(), k); // loop may be unrolled
 
             // j = index of k-th element of thread's private cache partition
             qindex j = getThreadsNthGlobalArrInd(k, t, numThreads);
@@ -393,7 +389,7 @@ __global__ void kernel_statevec_anyCtrlManyTargDenseMatr(
         for (qindex k=0; k<numTargAmps; k++) {
 
             // i = nth local index where ctrls are active and targs form value k
-            qindex i = setBits(i0, targs, numTargBits, k); // loop may be unrolled
+            qindex i = setBits(i0, targs.data(), targs.size(), k); // loop may be unrolled
             amps[i] = getGpuQcomp(0, 0);
         
             for (qindex l=0; l<numTargAmps; l++) {
@@ -431,8 +427,8 @@ __global__ void kernel_statevec_anyCtrlManyTargDenseMatr(
 
 template <int NumCtrls>
 __global__ void kernel_statevec_anyCtrlOneTargDiagMatr_sub(
-    gpu_qcomp* amps, qindex numThreads, int rank, qindex logNumAmpsPerNode,
-    int* ctrls, int numCtrls, qindex ctrlStateMask, int targ, 
+    gpu_qcomp* amps, qindex numThreads, int rank, qindex logNumAmpsPerNode, __grid_constant__ const List64 ctrl,
+    qindex ctrlStateMask, int targ, 
     gpu_qcomp m1, gpu_qcomp m2
 ) {
     GET_THREAD_IND(n, numThreads);
@@ -449,10 +445,10 @@ __global__ void kernel_statevec_anyCtrlOneTargDiagMatr_sub(
     /// We should verify this!
 
     // use template params to compile-time unroll loops in insertBits()
-    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
+    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrl.size());
 
     // j = nth local index where ctrls are active (in the specified states)
-    qindex j = insertBitsWithMaskedValues(n, ctrls, numCtrlBits, ctrlStateMask);
+    qindex j = insertBitsWithMaskedValues(n, ctrl.data(), numCtrlBits, ctrlStateMask);
 
     // i = global index corresponding to j
     qindex i = concatenateBits(rank, j, logNumAmpsPerNode);
@@ -470,8 +466,8 @@ __global__ void kernel_statevec_anyCtrlOneTargDiagMatr_sub(
 
 template <int NumCtrls>
 __global__ void kernel_statevec_anyCtrlTwoTargDiagMatr_sub(
-    gpu_qcomp* amps, qindex numThreads, int rank, qindex logNumAmpsPerNode,
-    int* ctrls, int numCtrls, qindex ctrlStateMask, int targ1, int targ2,
+    gpu_qcomp* amps, qindex numThreads, int rank, qindex logNumAmpsPerNode, __grid_constant__ const List64 ctrl,
+    qindex ctrlStateMask, int targ1, int targ2,
     gpu_qcomp m1, gpu_qcomp m2, gpu_qcomp m3, gpu_qcomp m4
 ) {
     GET_THREAD_IND(n, numThreads);
@@ -488,10 +484,10 @@ __global__ void kernel_statevec_anyCtrlTwoTargDiagMatr_sub(
     /// We should verify this!
 
     // use template params to compile-time unroll loops in insertBits()
-    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
+    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrl.size());
 
     // j = nth local index where ctrls are active (in the specified states)
-    qindex j = insertBitsWithMaskedValues(n, ctrls, numCtrlBits, ctrlStateMask);
+    qindex j = insertBitsWithMaskedValues(n, ctrl.data(), numCtrlBits, ctrlStateMask);
 
     // i = global index corresponding to j
     qindex i = concatenateBits(rank, j, logNumAmpsPerNode);
@@ -511,8 +507,8 @@ __global__ void kernel_statevec_anyCtrlTwoTargDiagMatr_sub(
 
 template <int NumCtrls, int NumTargs, bool ApplyConj, bool HasPower>
 __global__ void kernel_statevec_anyCtrlAnyTargDiagMatr_sub(
-    gpu_qcomp* amps, qindex numThreads, int rank, qindex logNumAmpsPerNode,
-    int* ctrls, int numCtrls, qindex ctrlStateMask, int* targs, int numTargs,
+    gpu_qcomp* amps, qindex numThreads, int rank, qindex logNumAmpsPerNode, __grid_constant__ const List64 ctrl,
+    qindex ctrlStateMask, __grid_constant__ const List64 targs,
     gpu_qcomp* elems, gpu_qcomp exponent
 ) {
     GET_THREAD_IND(n, numThreads);
@@ -529,17 +525,17 @@ __global__ void kernel_statevec_anyCtrlAnyTargDiagMatr_sub(
     /// We should verify this!
 
     // use template params to compile-time unroll loops in insertBits() and getValueOfBits()
-    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
-    SET_VAR_AT_COMPILE_TIME(int, numTargBits, NumTargs, numTargs);
+    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrl.size());
+    SET_VAR_AT_COMPILE_TIME(int, numTargBits, NumTargs, targs.size());
 
     // j = nth local index where ctrls are active (in the specified states)
-    qindex j = insertBitsWithMaskedValues(n, ctrls, numCtrlBits, ctrlStateMask);
+    qindex j = insertBitsWithMaskedValues(n, ctrl.data(), numCtrlBits, ctrlStateMask);
 
     // i = global index corresponding to j
     qindex i = concatenateBits(rank, j, logNumAmpsPerNode);
 
     // t = value of targeted bits, which may be in the prefix substate
-    qindex t = getValueOfBits(i, targs, numTargBits);
+    qindex t = getValueOfBits(i, targs.data(), numTargBits);
 
     gpu_qcomp elem = elems[t];
 
@@ -607,15 +603,15 @@ __global__ void kernel_densmatr_allTargDiagMatr_sub(
 template <int NumCtrls, int NumTargs> 
 __global__ void kernel_statevector_anyCtrlPauliTensorOrGadget_subA(
     gpu_qcomp* amps, qindex numThreads,
-    int* ctrlsAndTargs, int numCtrls, qindex ctrlsAndTargsStateMask, 
-    int* targsXY, int numXY, qindex maskXY, qindex maskYZ, 
+    __grid_constant__ const List64 ctrlsAndTargs, qindex ctrlsAndTargsStateMask, 
+    __grid_constant__ const List64 targsXY, qindex maskXY, qindex maskYZ, 
     gpu_qcomp powI, gpu_qcomp ampFac, gpu_qcomp pairAmpFac
 ) {
     GET_THREAD_IND(t, numThreads);
 
     // use template params to compile-time unroll loops in insertBits() and setBits()
-    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
-    SET_VAR_AT_COMPILE_TIME(int, numTargBits, NumTargs, numXY);
+    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrlsAndTargs.size());
+    SET_VAR_AT_COMPILE_TIME(int, numTargBits, NumTargs, targsXY.size());
 
     // n = local index of amp sub-batch with common i0, v = value of target bits
     qindex numInnerIts = powerOf2(numTargBits) / 2;
@@ -623,10 +619,10 @@ __global__ void kernel_statevector_anyCtrlPauliTensorOrGadget_subA(
     qindex v = t % numInnerIts;
 
     // i0 = nth local index where ctrls are active and targs are all zero (loop therein may be unrolled)
-    qindex i0 = insertBitsWithMaskedValues(n, ctrlsAndTargs, numCtrlBits + numTargBits, ctrlsAndTargsStateMask);
+    qindex i0 = insertBitsWithMaskedValues(n, ctrlsAndTargs.data(), numCtrlBits + numTargBits, ctrlsAndTargsStateMask);
 
     // iA = nth local index where targs have value v, iB = (last - nth) such index
-    qindex iA = setBits(i0, targsXY, numTargBits, v); // may be unrolled
+    qindex iA = setBits(i0, targsXY.data(), numTargBits, v); // may be unrolled
     qindex iB = flipBits(iA, maskXY);
 
     // determine whether to multiply amps by +-1 or +-i
@@ -647,17 +643,17 @@ __global__ void kernel_statevector_anyCtrlPauliTensorOrGadget_subA(
 template <int NumCtrls>
 __global__ void kernel_statevector_anyCtrlPauliTensorOrGadget_subB(
     gpu_qcomp* amps, gpu_qcomp* buffer, qindex numThreads,
-    int* ctrls, int numCtrls, qindex ctrlStateMask,
+    __grid_constant__ const List64 ctrls, qindex ctrlStateMask,
     qindex maskXY, qindex maskYZ, qindex bufferMaskXY,
     gpu_qcomp powI, gpu_qcomp thisAmpFac, gpu_qcomp otherAmpFac
 ) {
     GET_THREAD_IND(n, numThreads);
 
     // use template param to compile-time unroll loop in insertBits()
-    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
+    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrls.size());
 
     // i = nth local index where ctrl bits are in specified states
-    qindex i = insertBitsWithMaskedValues(n, ctrls, numCtrlBits, ctrlStateMask);
+    qindex i = insertBitsWithMaskedValues(n, ctrls.data(), numCtrlBits, ctrlStateMask);
 
     // j = buffer index of amp to be mixed with i
     qindex j = flipBits(n, bufferMaskXY);
@@ -682,16 +678,16 @@ __global__ void kernel_statevector_anyCtrlPauliTensorOrGadget_subB(
 template <int NumCtrls>
 __global__ void kernel_statevector_anyCtrlAnyTargZOrPhaseGadget_sub(
     gpu_qcomp* amps, qindex numThreads,
-    int* ctrls, int numCtrls, qindex ctrlStateMask, qindex targMask,
+    __grid_constant__ const List64 ctrls, qindex ctrlStateMask, qindex targMask,
     gpu_qcomp fac0, gpu_qcomp fac1
 ) {
     GET_THREAD_IND(n, numThreads);
 
     // use template param to compile-time unroll loop in insertBits()
-    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
+    SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrls.size());
 
     // i = nth local index where ctrl bits are in specified states
-    qindex i = insertBitsWithMaskedValues(n, ctrls, numCtrlBits, ctrlStateMask);
+    qindex i = insertBitsWithMaskedValues(n, ctrls.data(), numCtrlBits, ctrlStateMask);
 
     // apply phase to amp depending on parity of targets in global index 
     int p = cudaGetBitMaskParity(i & targMask);
@@ -1130,12 +1126,14 @@ __global__ void kernel_densmatr_oneQubitDamping_subD(
 template <int NumTargs>
 __global__ void kernel_densmatr_partialTrace_sub(
     gpu_qcomp* ampsIn, gpu_qcomp* ampsOut, qindex numThreads,
-    int* ketTargs, int* pairTargs, int* allTargs, int numKetTargs
+    __grid_constant__ const List64 ketTargs,
+    __grid_constant__ const List64 pairTargs,
+    __grid_constant__ const List64 allTargs
 ) {
     GET_THREAD_IND(n, numThreads);
 
     // use template param to compile-time unroll below loops
-    SET_VAR_AT_COMPILE_TIME(int, numTargPairs, NumTargs, numKetTargs);
+    SET_VAR_AT_COMPILE_TIME(int, numTargPairs, NumTargs, ketTargs.size());
 
     // may be inferred at compile-time
     int numAllTargs = 2 * numTargPairs;
@@ -1147,7 +1145,7 @@ __global__ void kernel_densmatr_partialTrace_sub(
     /// should change the parallelisation axis in this scenario, or preclude it with validation!
 
     // k = nth local index of inQureg where all targs and pairs are zero
-    qindex k = insertBits(n, allTargs, numAllTargs, 0); // loop may be unrolled
+    qindex k = insertBits(n, allTargs.data(), numAllTargs, 0); // loop may be unrolled
 
     // each outQureg amp results from summing 2^targs inQureg amps
     gpu_qcomp outAmp = getGpuQcomp(0, 0);
@@ -1157,8 +1155,8 @@ __global__ void kernel_densmatr_partialTrace_sub(
 
         // i = nth local index of inQureg where targs=j and pairTargs=j
         qindex i = k;
-        i = setBits(i, ketTargs,  numTargPairs, j); // loops may be unrolled
-        i = setBits(i, pairTargs, numTargPairs, j);
+        i = setBits(i, ketTargs.data(),  numTargPairs, j); // loops may be unrolled
+        i = setBits(i, pairTargs.data(), numTargPairs, j);
 
         outAmp += ampsIn[i];
     }
@@ -1177,7 +1175,7 @@ template <int NumQubits>
 __global__ void kernel_statevec_calcProbsOfAllMultiQubitOutcomes_sub(
     qreal* outProbs, gpu_qcomp* amps, qindex numThreads, 
     int rank, qindex logNumAmpsPerNode,
-    int* qubits, int numQubits
+    __grid_constant__ const List64 qubits
 ) {
     GET_THREAD_IND(n, numThreads);
 
@@ -1188,7 +1186,7 @@ __global__ void kernel_statevec_calcProbsOfAllMultiQubitOutcomes_sub(
     /// whether this is worthwhile and faster!
 
     // use template param to compile-time unroll below loops
-    SET_VAR_AT_COMPILE_TIME(int, numBits, NumQubits, numQubits);
+    SET_VAR_AT_COMPILE_TIME(int, numBits, NumQubits, qubits.size());
 
     qreal prob = norm(amps[n]);
 
@@ -1196,7 +1194,7 @@ __global__ void kernel_statevec_calcProbsOfAllMultiQubitOutcomes_sub(
     qindex i = concatenateBits(rank, n, logNumAmpsPerNode);
 
     // j = outcome index corresponding to prob
-    qindex j = getValueOfBits(i, qubits, numBits); // loop therein may be unrolled
+    qindex j = getValueOfBits(i, qubits.data(), numBits); // loop therein may be unrolled
 
     atomicAdd(&outProbs[j], prob);
 }
@@ -1207,12 +1205,12 @@ __global__ void kernel_densmatr_calcProbsOfAllMultiQubitOutcomes_sub(
     qreal* outProbs, gpu_qcomp* amps, qindex numThreads, 
     qindex firstDiagInd, qindex numAmpsPerCol,
     int rank, qindex logNumAmpsPerNode,
-    int* qubits, int numQubits
+    __grid_constant__ const List64 qubits
 ) {
     GET_THREAD_IND(n, numThreads);
 
     // use template param to compile-time unroll loop in insertBits()
-    SET_VAR_AT_COMPILE_TIME(int, numBits, NumQubits, numQubits);
+    SET_VAR_AT_COMPILE_TIME(int, numBits, NumQubits, qubits.size());
 
     // i = index of nth local diagonal elem
     qindex i = fast_getQuregLocalIndexOfDiagonalAmp(n, firstDiagInd, numAmpsPerCol);
@@ -1222,7 +1220,7 @@ __global__ void kernel_densmatr_calcProbsOfAllMultiQubitOutcomes_sub(
     qindex j = concatenateBits(rank, i, logNumAmpsPerNode);
 
     // k = outcome index corresponding to 
-    qindex k = getValueOfBits(j, qubits, numBits); // loop therein may be unrolled
+    qindex k = getValueOfBits(j, qubits.data(), numBits); // loop therein may be unrolled
 
     atomicAdd(&outProbs[k], prob);
 }
diff --git a/quest/src/gpu/gpu_subroutines.cpp b/quest/src/gpu/gpu_subroutines.cpp
index 9b8e819b5..9e3d2780d 100644
--- a/quest/src/gpu/gpu_subroutines.cpp
+++ b/quest/src/gpu/gpu_subroutines.cpp
@@ -144,13 +144,12 @@ qindex gpu_statevec_packAmpsIntoBuffer(Qureg qureg, ConstList64 qubits, ConstLis
     qindex numBlocks = getNumBlocks(numThreads, numThreadsPerBlock);
     qindex sendInd = getSubBufferSendInd(qureg);
 
-    devints sortedQubits = getDevInts(util_getSorted(qubits));
+    List64 sortedQubits = util_getSorted(qubits);
     qindex qubitStateMask  = util_getBitMask(qubits, qubitStates);
 
     kernel_statevec_packAmpsIntoBuffer <NumQubits> <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + sendInd, numThreads, 
-        getPtr(sortedQubits), qubits.size(), qubitStateMask
-    );
+        sortedQubits, qubitStateMask);
 
     // return the number of packed amps
     return numThreads;
@@ -212,12 +211,12 @@ void gpu_statevec_anyCtrlSwap_subA(Qureg qureg, ConstList64 ctrls, ConstList64 c
     int numThreadsPerBlock = gpu_getNumThreadsPerBlock();
     qindex numBlocks = getNumBlocks(numThreads, numThreadsPerBlock);
 
-    devints sortedQubits = getDevInts(util_getSorted(ctrls, {targ2, targ1}));
+    List64 sortedQubits = util_getSorted(ctrls, {targ2, targ1});
     qindex qubitStateMask = util_getBitMask(ctrls, ctrlStates, {targ2, targ1}, {0, 1});
 
     kernel_statevec_anyCtrlSwap_subA <NumCtrls> <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(qureg.gpuAmps), numThreads, 
-        getPtr(sortedQubits), ctrls.size(), qubitStateMask, targ1, targ2
+        sortedQubits, qubitStateMask, targ1, targ2
     );
 
 #else
@@ -238,12 +237,12 @@ void gpu_statevec_anyCtrlSwap_subB(Qureg qureg, ConstList64 ctrls, ConstList64 c
     qindex numBlocks = getNumBlocks(numThreads, numThreadsPerBlock);
     qindex recvInd = getBufferRecvInd();
 
-    devints sortedCtrls = getDevInts(util_getSorted(ctrls));
+    List64 sortedCtrls = util_getSorted(ctrls);
     qindex ctrlStateMask = util_getBitMask(ctrls, ctrlStates);
 
     kernel_statevec_anyCtrlSwap_subB <NumCtrls> <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + recvInd, numThreads, 
-        getPtr(sortedCtrls), ctrls.size(), ctrlStateMask
+        sortedCtrls, ctrlStateMask
     );
 
 #else
@@ -264,13 +263,12 @@ void gpu_statevec_anyCtrlSwap_subC(Qureg qureg, ConstList64 ctrls, ConstList64 c
     qindex numBlocks = getNumBlocks(numThreads, numThreadsPerBlock);
     qindex recvInd = getBufferRecvInd();
 
-    devints sortedQubits = getDevInts(util_getSorted(ctrls, {targ}));
+    List64 sortedQubits = util_getSorted(ctrls, {targ});
     qindex qubitStateMask = util_getBitMask(ctrls, ctrlStates, {targ}, {targState});
 
     kernel_statevec_anyCtrlSwap_subC <NumCtrls> <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + recvInd, numThreads, 
-        getPtr(sortedQubits), ctrls.size(), qubitStateMask
-    );
+        sortedQubits, qubitStateMask);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -306,17 +304,17 @@ void gpu_statevec_anyCtrlOneTargDenseMatr_subA(Qureg qureg, ConstList64 ctrls, C
     qindex numThreads = qureg.numAmpsPerNode / powerOf2(ctrls.size() + 1);
     int numThreadsPerBlock = gpu_getNumThreadsPerBlock();
     qindex numBlocks = getNumBlocks(numThreads, numThreadsPerBlock);
+    
+    List64 sortedQubits = util_getSorted(ctrls, {targ}); 
 
-    devints sortedQubits = getDevInts(util_getSorted(ctrls, {targ}));
     qindex qubitStateMask = util_getBitMask(ctrls, ctrlStates, {targ}, {0});
 
-    auto [m00, m01, m10, m11] = getFlattenedGpuQcompMatrix<2>(matr.elems); // explicit template for MSVC, grr!
+    auto [m00, m01, m10, m11] = getFlattenedGpuQcompMatrix<2>(matr.elems); // explicit template for MSVC, grrr!
 
     kernel_statevec_anyCtrlOneTargDenseMatr_subA <NumCtrls> <<<numBlocks, numThreadsPerBlock>>> (
-        getGpuQcompPtr(qureg.gpuAmps), numThreads, 
-        getPtr(sortedQubits), ctrls.size(), qubitStateMask, targ, 
-        m00, m01, m10, m11
-    );
+        getGpuQcompPtr(qureg.gpuAmps), numThreads, sortedQubits,
+        qubitStateMask, targ, 
+        m00, m01, m10, m11);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -336,12 +334,12 @@ void gpu_statevec_anyCtrlOneTargDenseMatr_subB(Qureg qureg, ConstList64 ctrls, C
     qindex numBlocks = getNumBlocks(numThreads, numThreadsPerBlock);
     qindex recvInd = getBufferRecvInd();
 
-    devints sortedCtrls = getDevInts(util_getSorted(ctrls));
+    List64 sortedCtrls = util_getSorted(ctrls);
     qindex ctrlStateMask = util_getBitMask(ctrls, ctrlStates);
 
     kernel_statevec_anyCtrlOneTargDenseMatr_subB <NumCtrls> <<<numBlocks,numThreadsPerBlock>>> (
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + recvInd, numThreads, 
-        getPtr(sortedCtrls), ctrls.size(), ctrlStateMask, 
+        sortedCtrls, ctrlStateMask, 
         getGpuQcomp(fac0), getGpuQcomp(fac1)
     );
 
@@ -379,7 +377,7 @@ void gpu_statevec_anyCtrlTwoTargDenseMatr_sub(Qureg qureg, ConstList64 ctrls, Co
     int numThreadsPerBlock = gpu_getNumThreadsPerBlock();
     qindex numBlocks = getNumBlocks(numThreads, numThreadsPerBlock);
 
-    devints sortedQubits = getDevInts(util_getSorted(ctrls, {targ1,targ2}));
+    List64 sortedQubits = util_getSorted(ctrls, {targ1,targ2});
     qindex qubitStateMask = util_getBitMask(ctrls, ctrlStates, {targ1,targ2}, {0,0});
 
     // unpack matrix elems which are more efficiently accessed by kernels as args than shared mem (... maybe...)
@@ -387,7 +385,7 @@ void gpu_statevec_anyCtrlTwoTargDenseMatr_sub(Qureg qureg, ConstList64 ctrls, Co
 
     kernel_statevec_anyCtrlTwoTargDenseMatr_sub <NumCtrls> <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(qureg.gpuAmps), numThreads, 
-        getPtr(sortedQubits), ctrls.size(), qubitStateMask, targ1, targ2,
+        sortedQubits, qubitStateMask, targ1, targ2,
         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]
     );
@@ -441,16 +439,15 @@ void gpu_statevec_anyCtrlAnyTargDenseMatr_sub(Qureg qureg, ConstList64 ctrls, Co
     // task each thread with processing more than a single batch
     qindex numBatches = qureg.numAmpsPerNode / powerOf2(ctrls.size() + targs.size());
 
-    devints deviceTargs = getDevInts(targs);
-    devints deviceQubits = getDevInts(util_getSorted(ctrls, targs));
+    List64 sortedQubits = util_getSorted(ctrls, targs);
     qindex qubitStateMask = util_getBitMask(ctrls, ctrlStates, targs, util_getConstantList(0,targs.size()));
 
     // unpacking args (to better distinguish below signatures)
     auto ampsPtr   = getGpuQcompPtr(qureg.gpuAmps);
     auto matrPtr   = getGpuQcompPtr(matr.gpuElemsFlat);
-    auto qubitsPtr = getPtr(deviceQubits);
-    auto targsPtr  = getPtr(deviceTargs);
-    auto nCtrls    = ctrls.size();
+    // auto qubitsPtr = getPtr(deviceQubits);
+    // auto targsPtr  = getPtr(deviceTargs);
+    // auto nCtrls    = ctrls.size();
 
     // this function updates amplitudes in batches of 2^NumTargs, where each is
     // determined by distinct mixtures of the existing 2^NumTargs values, which
@@ -482,8 +479,8 @@ void gpu_statevec_anyCtrlAnyTargDenseMatr_sub(Qureg qureg, ConstList64 ctrls, Co
             <NumCtrls, NumTargs, ApplyConj, ApplyTransp> 
             <<<numBlocks, numThreadsPerBlock>>> (
                 ampsPtr, numThreads, 
-                qubitsPtr, nCtrls, qubitStateMask, 
-                targsPtr, matrPtr
+                sortedQubits, qubitStateMask, 
+                targs, matrPtr
         );
 
     } else {
@@ -522,8 +519,8 @@ void gpu_statevec_anyCtrlAnyTargDenseMatr_sub(Qureg qureg, ConstList64 ctrls, Co
             <<<numBlocks, numThreadsPerBlock>>> (
                 getGpuQcompPtr(cache),
                 ampsPtr, numThreads, numBatchesPerThread, 
-                qubitsPtr, nCtrls, qubitStateMask, 
-                targsPtr, targs.size(), powerOf2(targs.size()), matrPtr
+                sortedQubits, qubitStateMask, 
+                targs, powerOf2(targs.size()), matrPtr
         );
     }
 
@@ -584,13 +581,15 @@ void gpu_statevec_anyCtrlOneTargDiagMatr_sub(Qureg qureg, ConstList64 ctrls, Con
     int numThreadsPerBlock = gpu_getNumThreadsPerBlock();
     qindex numBlocks = getNumBlocks(numThreads, numThreadsPerBlock);
 
-    devints deviceCtrls = getDevInts(util_getSorted(ctrls));
+    // removed implicit thrust mem copy
+    List64 sortedCtrls = util_getSorted(ctrls);
+
     qindex ctrlStateMask = util_getBitMask(ctrls, ctrlStates);
     auto elems = getGpuQcompArray<2>(matr.elems); // explicit template for MSVC, grr!
 
     kernel_statevec_anyCtrlOneTargDiagMatr_sub <NumCtrls> <<<numBlocks, numThreadsPerBlock>>> (
-        getGpuQcompPtr(qureg.gpuAmps), numThreads, qureg.rank, qureg.logNumAmpsPerNode,
-        getPtr(deviceCtrls), ctrls.size(), ctrlStateMask, targ, elems[0], elems[1]
+        getGpuQcompPtr(qureg.gpuAmps), numThreads, qureg.rank, qureg.logNumAmpsPerNode, sortedCtrls,
+       	ctrlStateMask, targ, elems[0], elems[1]
     );
 
     // explicitly return to avoid runtime error below
@@ -655,13 +654,14 @@ void gpu_statevec_anyCtrlTwoTargDiagMatr_sub(Qureg qureg, ConstList64 ctrls, Con
     int numThreadsPerBlock = gpu_getNumThreadsPerBlock();
     qindex numBlocks = getNumBlocks(numThreads, numThreadsPerBlock);
 
-    devints deviceCtrls = getDevInts(util_getSorted(ctrls));
+    List64 sortedCtrls = util_getSorted(ctrls);
+    
     qindex ctrlStateMask = util_getBitMask(ctrls, ctrlStates);
     auto elems = getGpuQcompArray<4>(matr.elems); // explicit template for MSVC, grr!
 
     kernel_statevec_anyCtrlTwoTargDiagMatr_sub <NumCtrls> <<<numBlocks, numThreadsPerBlock>>> (
-        getGpuQcompPtr(qureg.gpuAmps), numThreads, qureg.rank, qureg.logNumAmpsPerNode,
-        getPtr(deviceCtrls), ctrls.size(), ctrlStateMask, targ1, targ2,
+        getGpuQcompPtr(qureg.gpuAmps), numThreads, qureg.rank, qureg.logNumAmpsPerNode, sortedCtrls,
+        ctrlStateMask, targ1, targ2,
         elems[0], elems[1], elems[2], elems[3]
     );
 
@@ -724,13 +724,12 @@ void gpu_statevec_anyCtrlAnyTargDiagMatr_sub(Qureg qureg, ConstList64 ctrls, Con
     int numThreadsPerBlock = gpu_getNumThreadsPerBlock();
     qindex numBlocks = getNumBlocks(numThreads, numThreadsPerBlock);
 
-    devints deviceTargs = getDevInts(targs);
-    devints deviceCtrls = getDevInts(util_getSorted(ctrls));
+    List64 sortedCtrls = util_getSorted(ctrls);
     qindex ctrlStateMask = util_getBitMask(ctrls, ctrlStates);
 
     kernel_statevec_anyCtrlAnyTargDiagMatr_sub <NumCtrls, NumTargs, ApplyConj, HasPower> <<<numBlocks, numThreadsPerBlock>>> (
-        getGpuQcompPtr(qureg.gpuAmps), numThreads, qureg.rank, qureg.logNumAmpsPerNode,
-        getPtr(deviceCtrls), ctrls.size(), ctrlStateMask, getPtr(deviceTargs), targs.size(), 
+        getGpuQcompPtr(qureg.gpuAmps), numThreads, qureg.rank, qureg.logNumAmpsPerNode, sortedCtrls,
+        ctrlStateMask, targs, 
         getGpuQcompPtr(util_getGpuMemPtr(matr)), getGpuQcomp(exponent)
     );
 
@@ -831,9 +830,8 @@ void gpu_statevector_anyCtrlPauliTensorOrGadget_subA(Qureg qureg, ConstList64 ct
     auto targsXY = util_getConcatenated(x, y);
     auto maskXY  = util_getBitMask(targsXY);
     auto maskYZ  = util_getBitMask(util_getConcatenated(y, z));
-
-    devints deviceTargs   = getDevInts(targsXY);
-    devints deviceQubits  = getDevInts(util_getSorted(ctrls, targsXY));
+     
+    List64 sortedQubits  = util_getSorted(ctrls, targsXY);
     qindex qubitStateMask = util_getBitMask(ctrls, ctrlStates, targsXY, util_getConstantList(0,targsXY.size()));
 
     // unlike the analogous cpu routine, this function has only a single parallelisation
@@ -846,9 +844,9 @@ void gpu_statevector_anyCtrlPauliTensorOrGadget_subA(Qureg qureg, ConstList64 ct
     qindex numBlocks = getNumBlocks(numThreads, numThreadsPerBlock);
     kernel_statevector_anyCtrlPauliTensorOrGadget_subA <NumCtrls, NumTargs> <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(qureg.gpuAmps), numThreads,
-        getPtr(deviceQubits), ctrls.size(), qubitStateMask, 
-        getPtr(deviceTargs), deviceTargs.size(),
-        maskXY, maskYZ, getGpuQcomp(powI), getGpuQcomp(ampFac), getGpuQcomp(pairAmpFac)
+        sortedQubits, qubitStateMask, 
+        targsXY, maskXY, maskYZ,
+       	getGpuQcomp(powI), getGpuQcomp(ampFac), getGpuQcomp(pairAmpFac)
     );
 
 #else
@@ -873,12 +871,12 @@ void gpu_statevector_anyCtrlPauliTensorOrGadget_subB(Qureg qureg, ConstList64 ct
     auto maskXY = util_getBitMask(util_getConcatenated(x, y));
     auto maskYZ = util_getBitMask(util_getConcatenated(y, z));
 
-    devints sortedCtrls = getDevInts(util_getSorted(ctrls));
+    List64 sortedCtrls = util_getSorted(ctrls);
     qindex ctrlStateMask = util_getBitMask(ctrls, ctrlStates);
 
     kernel_statevector_anyCtrlPauliTensorOrGadget_subB <NumCtrls> <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + recvInd, numThreads, 
-        getPtr(sortedCtrls), ctrls.size(), ctrlStateMask,
+        sortedCtrls, ctrlStateMask,
         maskXY, maskYZ, bufferMaskXY,
         getGpuQcomp(powI), getGpuQcomp(ampFac), getGpuQcomp(pairAmpFac)
     );
@@ -910,13 +908,13 @@ void gpu_statevector_anyCtrlAnyTargZOrPhaseGadget_sub(Qureg qureg, ConstList64 c
     int numThreadsPerBlock = gpu_getNumThreadsPerBlock();
     qindex numBlocks = getNumBlocks(numThreads, numThreadsPerBlock);
 
-    devints sortedCtrls = getDevInts(util_getSorted(ctrls));
+    List64 sortedCtrls = util_getSorted(ctrls); // change for performance
     qindex ctrlStateMask = util_getBitMask(ctrls, ctrlStates);
     qindex targMask = util_getBitMask(targs);
 
     kernel_statevector_anyCtrlAnyTargZOrPhaseGadget_sub <NumCtrls> <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(qureg.gpuAmps), numThreads,
-        getPtr(sortedCtrls), ctrls.size(), ctrlStateMask, targMask,
+        sortedCtrls, ctrlStateMask, targMask,
         getGpuQcomp(fac0), getGpuQcomp(fac1)
     );
 
@@ -951,8 +949,8 @@ void gpu_statevec_setQuregToWeightedSum_sub(Qureg outQureg, vector<qcomp> coeffs
         ptrs.push_back(getGpuQcompPtr(qureg.gpuAmps));
     
     // copy coeff and qureg lists into GPU memory
-    devgpuqcompptrs devQuregAmps = ptrs;
-    devcomps devCoeffs = coeffs;
+    devgpuqcompptrs devQuregAmps = ptrs; // review performance
+    devcomps devCoeffs = coeffs;  // review performance
 
     kernel_statevec_setQuregToWeightedSum_sub <NumQuregs> <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(outQureg.gpuAmps), numThreads,
@@ -1484,13 +1482,11 @@ void gpu_densmatr_partialTrace_sub(Qureg inQureg, Qureg outQureg, ConstList64 ta
     int numThreadsPerBlock = gpu_getNumThreadsPerBlock();
     qindex numBlocks = getNumBlocks(numThreads, numThreadsPerBlock);
 
-    devints devTargs = getDevInts(targs);
-    devints devPairTargs = getDevInts(pairTargs);
-    devints devAllTargs = getDevInts(util_getSorted(targs, pairTargs));
+    List64 allTargs = util_getSorted(targs, pairTargs);
 
     kernel_densmatr_partialTrace_sub <NumTargs> <<<numBlocks, numThreadsPerBlock>>> (
         getGpuQcompPtr(inQureg.gpuAmps), getGpuQcompPtr(outQureg.gpuAmps), numThreads,
-        getPtr(devTargs), getPtr(devPairTargs), getPtr(devAllTargs), targs.size()
+        targs, pairTargs, allTargs
     );
 
 #else
@@ -1606,12 +1602,11 @@ void gpu_statevec_calcProbsOfAllMultiQubitOutcomes_sub(qreal* outProbs, Qureg qu
     qindex numBlocks = getNumBlocks(numThreads, numThreadsPerBlock);
 
     // allocate exponentially-big temporary memory (error if failed)
-    devints devQubits = getDevInts(qubits);
     devreals devProbs = getDeviceRealsVec(powerOf2(qubits.size())); // throws
 
     kernel_statevec_calcProbsOfAllMultiQubitOutcomes_sub<NumQubits> <<<numBlocks, numThreadsPerBlock>>> (
-        getPtr(devProbs), getGpuQcompPtr(qureg.gpuAmps), numThreads, 
-        qureg.rank, qureg.logNumAmpsPerNode, getPtr(devQubits), devQubits.size()
+       getPtr(devProbs), getGpuQcompPtr(qureg.gpuAmps), numThreads, 
+        qureg.rank, qureg.logNumAmpsPerNode, qubits
     );
 
     // overwrite outProbs with GPU memory
@@ -1644,14 +1639,13 @@ void gpu_densmatr_calcProbsOfAllMultiQubitOutcomes_sub(qreal* outProbs, Qureg qu
     qindex numAmpsPerCol = powerOf2(qureg.numQubits);
 
     // allocate exponentially-big temporary memory (error if failed)
-    devints devQubits = getDevInts(qubits);
     devreals devProbs = getDeviceRealsVec(powerOf2(qubits.size())); // throws
 
     kernel_densmatr_calcProbsOfAllMultiQubitOutcomes_sub<NumQubits> <<<numBlocks, numThreadsPerBlock>>> (
-        getPtr(devProbs), getGpuQcompPtr(qureg.gpuAmps), 
+       getPtr(devProbs), getGpuQcompPtr(qureg.gpuAmps), 
         numThreads, firstDiagInd, numAmpsPerCol,
         qureg.rank, qureg.logNumAmpsPerNode, 
-        getPtr(devQubits), devQubits.size()
+        qubits
     );
 
     // overwrite outProbs with GPU memory
diff --git a/quest/src/gpu/gpu_thrust.cuh b/quest/src/gpu/gpu_thrust.cuh
index 864cca5f8..5fb6063e2 100644
--- a/quest/src/gpu/gpu_thrust.cuh
+++ b/quest/src/gpu/gpu_thrust.cuh
@@ -69,7 +69,7 @@
  */
 
 
-using devints = thrust::device_vector<int>;
+using devints = thrust::device_vector<int>; // remove for performance
 
 devints getDevInts(ConstList64 h_list) {
 
@@ -814,6 +814,7 @@ qreal thrust_densmatr_calcProbOfMultiQubitOutcome_sub(Qureg qureg, ConstList64 q
 
     // cannot move these into functor_insertBits constructor, since the memory
     // would dangle - and we cannot bind deviceints as an attribute - it's host-only!
+
     devints sortedQubits = getDevInts(util_getSorted(qubits));
     qindex valueMask = util_getBitMask(qubits, outcomes);
 
@@ -1089,4 +1090,4 @@ void thrust_statevec_initUnnormalisedUniformlyRandomPureStateAmps_sub(Qureg qure
 
 
 
-#endif // GPU_THRUST_HPP
\ No newline at end of file
+#endif // GPU_THRUST_HPP
diff --git a/tests/unit/experimental.cpp b/tests/unit/experimental.cpp
index 943645831..a8497eb8f 100644
--- a/tests/unit/experimental.cpp
+++ b/tests/unit/experimental.cpp
@@ -84,7 +84,7 @@ TEST_CASE( "setQuESTNumGpuThreadsPerBlock", TEST_CATEGORY ) {
 
         SECTION( "Exceeds device maximum" ) {
 
-            int badNumTPB = 999999; // exceeds expected 1024 max
+            int badNumTPB = 102400; // exceeds expected 1024 max
 
             // Cannot be tested (since validation not imposed) when GPU is not actively used
             if (getQuESTEnv().isGpuAccelerated)