diff --git a/quest/src/core/bitwise.hpp b/quest/src/core/bitwise.hpp
index 4d455c2d..841e1cd1 100644
--- a/quest/src/core/bitwise.hpp
+++ b/quest/src/core/bitwise.hpp
@@ -172,6 +172,15 @@ INLINE qindex insertBits(qindex number, int* bitIndices, int numIndices, int bit
     return number;
 }
 
+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, int* bitIndices, int numIndices, qindex bitsValue) {
     
@@ -184,6 +193,17 @@ INLINE qindex setBits(qindex number, int* bitIndices, int numIndices, qindex bit
     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;
+}
+
 
 INLINE qindex getValueOfBits(qindex number, int* bitIndices, int numIndices) {
 
@@ -196,6 +216,16 @@ INLINE qindex getValueOfBits(qindex number, int* bitIndices, int numIndices) {
     return value;
 }
 
+INLINE qindex getValueOfBits(qindex number, const int* bitIndices, int numIndices) {
+
+    // bits are arbitrarily ordered, which affects value
+    qindex value = 0;
+
+    for (int i=0; i<numIndices; i++)
+        value |= getBit(number, bitIndices[i]) << i;
+
+    return value;
+}
 
 
 /*
@@ -211,6 +241,12 @@ INLINE qindex insertBitsWithMaskedValues(qindex number, int* bitInds, int numBit
     return mask | insertBits(number, bitInds, numBits, 0);
 }
 
+INLINE qindex insertBitsWithMaskedValues(qindex number, const int* bitInds, int numBits, qindex mask) {
+
+    // bitInds must be sorted (increasing), and mask must be zero everywhere except bitInds
+    return mask | insertBits(number, bitInds, numBits, 0);
+}
+
 
 INLINE int getTwoBits(qindex number, int highInd, int lowInd) {
 
@@ -379,4 +415,4 @@ INLINE void getBitsFromInteger(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 30377850..fd00f0c5 100644
--- a/quest/src/gpu/gpu_kernels.cuh
+++ b/quest/src/gpu/gpu_kernels.cuh
@@ -43,7 +43,19 @@
  */
 
 
-const int NUM_THREADS_PER_BLOCK = 128;
+const int NUM_THREADS_PER_BLOCK =128;
+
+// __device__ __constant__ int ctrl_device[30]; 
+
+
+typedef struct {
+    int ctrl_device[64];
+} ctrl_device_t;
+
+struct QubitList_t {
+    int indices[64];
+    int length;
+};
 
 
 __forceinline__ __device__ qindex getThreadInd() {
@@ -90,15 +102,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 QubitList_t 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.length);
 
     // i = nth local index where qubits are active
-    qindex i = insertBitsWithMaskedValues(n, qubits, numBits, qubitStateMask);
+    qindex i = insertBitsWithMaskedValues(n, qubits.indices, numBits, qubitStateMask);
 
     // caller offsets buffer by sub-buffer send-index
     buffer[n] = amps[i];
@@ -107,8 +118,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
@@ -129,16 +139,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 QubitList_t 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.length);
     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.indices, numQubitBits, ctrlsAndTargsMask);
     qindex i10 = flipTwoBits(i01, targ2, targ1);
 
     // swap amps
@@ -151,15 +160,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 QubitList_t 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.length);
 
     // i = nth local index where ctrls are active
-    qindex i = insertBitsWithMaskedValues(n, ctrls, numCtrlBits, ctrlStateMask);
+    qindex i = insertBitsWithMaskedValues(n, ctrls.indices, numCtrlBits, ctrlStateMask);
 
     // caller offsets buffer if necessary
     amps[i] = buffer[n];
@@ -169,16 +178,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 QubitList_t 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.length);
     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.indices, numQubitBits, ctrlsAndTargMask);
 
     // caller offsets buffer if necessary
     amps[i] = buffer[n];
@@ -193,8 +202,8 @@ __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 ctrl_device_t ctrl,
+    int numCtrls, qindex ctrlStateMask, int targ, 
     gpu_qcomp m00, gpu_qcomp m01, gpu_qcomp m10, gpu_qcomp m11
 ) {
     GET_THREAD_IND(n, numThreads);
@@ -203,7 +212,7 @@ __global__ void kernel_statevec_anyCtrlOneTargDenseMatr_subA(
     SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
 
     // 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.ctrl_device, numCtrlBits + 1, ctrlStateMask);
     qindex i1 = flipBit(i0, targ);
 
     // note amps are strided by 2^targ
@@ -218,16 +227,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 QubitList_t 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.length);
 
     // i = nth local index where ctrl bits are active
-    qindex i = insertBitsWithMaskedValues(n, ctrls, numCtrlBits, ctrlStateMask);
+    qindex i = insertBitsWithMaskedValues(n, ctrls.indices, numCtrlBits, ctrlStateMask);
 
     // caller offsets buffer by receive-index
     amps[i] = fac0*amps[i] + fac1*buffer[n];
@@ -243,7 +252,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 QubitList_t 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,
@@ -252,10 +261,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.length);
 
     // 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.indices, numCtrlBits + 2, ctrlStateMask);
     qindex i01 = flipBit(i00, targ1);
     qindex i10 = flipBit(i00, targ2);
     qindex i11 = flipBit(i01, targ2);
@@ -292,7 +301,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 QubitList_t ctrlsAndTargs, qindex ctrlsAndTargsMask, __grid_constant__ const QubitList_t targs,
     gpu_qcomp* flatMatrElems
 ) {
     GET_THREAD_IND(n, numThreads);
@@ -309,18 +318,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.length);
     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.indices, 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.indices, NumTargs, k); // loop will be unrolled
 
         // write to thread-private cache at compile-time known index
         privateCache[k] = amps[i];
@@ -330,7 +339,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.indices, NumTargs, k); // loop will be unrolled
         amps[i] = getGpuQcomp(0, 0);
     
         // force unroll to ensure compile-time cache indices
@@ -360,14 +369,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 QubitList_t ctrlsAndTargs, qindex ctrlsAndTargsMask, 
+    __grid_constant__ const QubitList_t 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.length);
 
     // unlike all other kernels, each thread modifies multiple batches of amplitudes
     for (qindex b=0; b<numBatchesPerThread; b++) {
@@ -376,13 +384,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.indices, numCtrlBits + targs.length, 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.indices, targs.length, 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 +401,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.indices, targs.length, k); // loop may be unrolled
             amps[i] = getGpuQcomp(0, 0);
         
             for (qindex l=0; l<numTargAmps; l++) {
@@ -431,8 +439,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 ctrl_device_t ctrl,
+    int numCtrls, qindex ctrlStateMask, int targ, 
     gpu_qcomp m1, gpu_qcomp m2
 ) {
     GET_THREAD_IND(n, numThreads);
@@ -452,7 +460,7 @@ __global__ void kernel_statevec_anyCtrlOneTargDiagMatr_sub(
     SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
 
     // j = nth local index where ctrls are active (in the specified states)
-    qindex j = insertBitsWithMaskedValues(n, ctrls, numCtrlBits, ctrlStateMask);
+    qindex j = insertBitsWithMaskedValues(n, ctrl.ctrl_device, numCtrlBits, ctrlStateMask);
 
     // i = global index corresponding to j
     qindex i = concatenateBits(rank, j, logNumAmpsPerNode);
@@ -470,8 +478,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 ctrl_device_t ctrl,
+    int numCtrls, qindex ctrlStateMask, int targ1, int targ2,
     gpu_qcomp m1, gpu_qcomp m2, gpu_qcomp m3, gpu_qcomp m4
 ) {
     GET_THREAD_IND(n, numThreads);
@@ -491,7 +499,7 @@ __global__ void kernel_statevec_anyCtrlTwoTargDiagMatr_sub(
     SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, numCtrls);
 
     // j = nth local index where ctrls are active (in the specified states)
-    qindex j = insertBitsWithMaskedValues(n, ctrls, numCtrlBits, ctrlStateMask);
+    qindex j = insertBitsWithMaskedValues(n, ctrl.ctrl_device, numCtrlBits, ctrlStateMask);
 
     // i = global index corresponding to j
     qindex i = concatenateBits(rank, j, logNumAmpsPerNode);
@@ -511,8 +519,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 ctrl_device_t ctrl,
+    int numCtrls, qindex ctrlStateMask, __grid_constant__ const QubitList_t targs,
     gpu_qcomp* elems, gpu_qcomp exponent
 ) {
     GET_THREAD_IND(n, numThreads);
@@ -530,16 +538,16 @@ __global__ void kernel_statevec_anyCtrlAnyTargDiagMatr_sub(
 
     // 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, numTargBits, NumTargs, targs.length);
 
     // j = nth local index where ctrls are active (in the specified states)
-    qindex j = insertBitsWithMaskedValues(n, ctrls, numCtrlBits, ctrlStateMask);
+    qindex j = insertBitsWithMaskedValues(n, ctrl.ctrl_device, 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.indices, numTargBits);
 
     gpu_qcomp elem = elems[t];
 
@@ -607,15 +615,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 QubitList_t ctrlsAndTargs, qindex ctrlsAndTargsStateMask, 
+    __grid_constant__ const QubitList_t 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.length);
+    SET_VAR_AT_COMPILE_TIME(int, numTargBits, NumTargs, targsXY.length);
 
     // n = local index of amp sub-batch with common i0, v = value of target bits
     qindex numInnerIts = powerOf2(numTargBits) / 2;
@@ -623,10 +631,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.indices, 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.indices, numTargBits, v); // may be unrolled
     qindex iB = flipBits(iA, maskXY);
 
     // determine whether to multiply amps by +-1 or +-i
@@ -647,17 +655,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 QubitList_t 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.length);
 
     // i = nth local index where ctrl bits are in specified states
-    qindex i = insertBitsWithMaskedValues(n, ctrls, numCtrlBits, ctrlStateMask);
+    qindex i = insertBitsWithMaskedValues(n, ctrls.indices, numCtrlBits, ctrlStateMask);
 
     // j = buffer index of amp to be mixed with i
     qindex j = flipBits(n, bufferMaskXY);
@@ -682,16 +690,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 QubitList_t 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.length);
 
     // i = nth local index where ctrl bits are in specified states
-    qindex i = insertBitsWithMaskedValues(n, ctrls, numCtrlBits, ctrlStateMask);
+    qindex i = insertBitsWithMaskedValues(n, ctrls.indices, numCtrlBits, ctrlStateMask);
 
     // apply phase to amp depending on parity of targets in global index 
     int p = cudaGetBitMaskParity(i & targMask);
@@ -1130,12 +1138,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 QubitList_t ketTargs,
+    __grid_constant__ const QubitList_t pairTargs,
+    __grid_constant__ const QubitList_t 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.length);
 
     // may be inferred at compile-time
     int numAllTargs = 2 * numTargPairs;
@@ -1147,7 +1157,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.indices, numAllTargs, 0); // loop may be unrolled
 
     // each outQureg amp results from summing 2^targs inQureg amps
     gpu_qcomp outAmp = getGpuQcomp(0, 0);
@@ -1157,8 +1167,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.indices,  numTargPairs, j); // loops may be unrolled
+        i = setBits(i, pairTargs.indices, numTargPairs, j);
 
         outAmp += ampsIn[i];
     }
@@ -1177,7 +1187,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 QubitList_t qubits
 ) {
     GET_THREAD_IND(n, numThreads);
 
@@ -1188,7 +1198,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.length);
 
     qreal prob = norm(amps[n]);
 
@@ -1196,7 +1206,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.indices, numBits); // loop therein may be unrolled
 
     atomicAdd(&outProbs[j], prob);
 }
@@ -1207,12 +1217,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 QubitList_t 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.length);
 
     // i = index of nth local diagonal elem
     qindex i = fast_getQuregLocalIndexOfDiagonalAmp(n, firstDiagInd, numAmpsPerCol);
@@ -1222,7 +1232,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.indices, 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 cd473ee1..862a29c3 100644
--- a/quest/src/gpu/gpu_subroutines.cpp
+++ b/quest/src/gpu/gpu_subroutines.cpp
@@ -144,12 +144,16 @@ qindex gpu_statevec_packAmpsIntoBuffer(Qureg qureg, vector<int> qubits, vector<i
     qindex numBlocks = getNumBlocks(numThreads);
     qindex sendInd = getSubBufferSendInd(qureg);
 
-    devints sortedQubits = util_getSorted(qubits);
+    vector<int> sortedQubits = util_getSorted(qubits); // change
     qindex qubitStateMask  = util_getBitMask(qubits, qubitStates);
 
+    QubitList_t qubits_dev;
+    std::copy(sortedQubits.begin(), sortedQubits.end(), qubits_dev.indices);
+    qubits_dev.length = sortedQubits.size();
+
     kernel_statevec_packAmpsIntoBuffer <NumQubits> <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + sendInd, numThreads, 
-        getPtr(sortedQubits), qubits.size(), qubitStateMask
+        qubits_dev, qubitStateMask
     );
 
     // return the number of packed amps
@@ -210,12 +214,16 @@ void gpu_statevec_anyCtrlSwap_subA(Qureg qureg, vector<int> ctrls, vector<int> c
     qindex numThreads = qureg.numAmpsPerNode / powerOf2(2 + ctrls.size());
     qindex numBlocks = getNumBlocks(numThreads);
 
-    devints sortedQubits = util_getSorted(ctrls, {targ2, targ1});
+    vector<int> sortedQubits = util_getSorted(ctrls, {targ2, targ1}); // change for performance
     qindex qubitStateMask = util_getBitMask(ctrls, ctrlStates, {targ2, targ1}, {0, 1});
 
+    QubitList_t qubits_dev;
+    std::copy(sortedQubits.begin(), sortedQubits.end(), qubits_dev.indices);
+    qubits_dev.length = sortedQubits.size();
+
     kernel_statevec_anyCtrlSwap_subA <NumCtrls> <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
         getGpuQcompPtr(qureg.gpuAmps), numThreads, 
-        getPtr(sortedQubits), ctrls.size(), qubitStateMask, targ1, targ2
+        qubits_dev, qubitStateMask, targ1, targ2
     );
 
 #else
@@ -235,12 +243,16 @@ void gpu_statevec_anyCtrlSwap_subB(Qureg qureg, vector<int> ctrls, vector<int> c
     qindex numBlocks = getNumBlocks(numThreads);
     qindex recvInd = getBufferRecvInd();
 
-    devints sortedCtrls = util_getSorted(ctrls);
+    vector<int> sortedCtrls = util_getSorted(ctrls); // change for performance
     qindex ctrlStateMask = util_getBitMask(ctrls, ctrlStates);
 
+    QubitList_t Ctrls_dev;
+    std::copy(sortedCtrls.begin(), sortedCtrls.end(), Ctrls_dev.indices);
+    Ctrls_dev.length = sortedCtrls.size();
+
     kernel_statevec_anyCtrlSwap_subB <NumCtrls> <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + recvInd, numThreads, 
-        getPtr(sortedCtrls), ctrls.size(), ctrlStateMask
+        Ctrls_dev, ctrlStateMask
     );
 
 #else
@@ -260,13 +272,16 @@ void gpu_statevec_anyCtrlSwap_subC(Qureg qureg, vector<int> ctrls, vector<int> c
     qindex numBlocks = getNumBlocks(numThreads);
     qindex recvInd = getBufferRecvInd();
 
-    devints sortedQubits = util_getSorted(ctrls, {targ});
+    vector<int> sortedQubits = util_getSorted(ctrls, {targ}); // change for performance
     qindex qubitStateMask = util_getBitMask(ctrls, ctrlStates, {targ}, {targState});
 
+    QubitList_t Qubits_dev;
+    std::copy(sortedQubits.begin(), sortedQubits.end(), Qubits_dev.indices);
+    Qubits_dev.length = sortedQubits.size();
+
     kernel_statevec_anyCtrlSwap_subC <NumCtrls> <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + recvInd, numThreads, 
-        getPtr(sortedQubits), ctrls.size(), qubitStateMask
-    );
+        Qubits_dev, qubitStateMask);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -301,16 +316,27 @@ void gpu_statevec_anyCtrlOneTargDenseMatr_subA(Qureg qureg, vector<int> ctrls, v
     qindex numThreads = qureg.numAmpsPerNode / powerOf2(ctrls.size() + 1);
     qindex numBlocks = getNumBlocks(numThreads);
 
-    devints sortedQubits = util_getSorted(ctrls, {targ});
+    //devints sortedQubits = util_getSorted(ctrls, {targ});
+    
+    vector<int> sortedQubits = 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!
 
+    ctrl_device_t ctrl; // change for performance Standardise
+
+    std::copy(sortedQubits.begin(), sortedQubits.end(), ctrl.ctrl_device);
+
+
+    //int ctrl_device[sortedQubits.size()];
+
+    //cudaMemcpyToSymbol(ctrl_device, sortedQubits.data(), sortedQubits.size()*sizeof(int));
+
     kernel_statevec_anyCtrlOneTargDenseMatr_subA <NumCtrls> <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
-        getGpuQcompPtr(qureg.gpuAmps), numThreads, 
-        getPtr(sortedQubits), ctrls.size(), qubitStateMask, targ, 
-        m00, m01, m10, m11
-    );
+        getGpuQcompPtr(qureg.gpuAmps), numThreads, ctrl,
+        ctrls.size(), qubitStateMask, targ, 
+        m00, m01, m10, m11);
 
 #else
     error_gpuSimButGpuNotCompiled();
@@ -329,12 +355,16 @@ void gpu_statevec_anyCtrlOneTargDenseMatr_subB(Qureg qureg, vector<int> ctrls, v
     qindex numBlocks = getNumBlocks(numThreads);
     qindex recvInd = getBufferRecvInd();
 
-    devints sortedCtrls = util_getSorted(ctrls);
+    vector<int> sortedCtrls = util_getSorted(ctrls); // change for performance
     qindex ctrlStateMask = util_getBitMask(ctrls, ctrlStates);
 
+    QubitList_t Ctrls_dev;
+    std::copy(sortedCtrls.begin(), sortedCtrls.end(), Ctrls_dev.indices);
+    Ctrls_dev.length = sortedCtrls.size();
+
     kernel_statevec_anyCtrlOneTargDenseMatr_subB <NumCtrls> <<<numBlocks,NUM_THREADS_PER_BLOCK>>> (
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + recvInd, numThreads, 
-        getPtr(sortedCtrls), ctrls.size(), ctrlStateMask, 
+        Ctrls_dev, ctrlStateMask, 
         getGpuQcomp(fac0), getGpuQcomp(fac1)
     );
 
@@ -370,15 +400,19 @@ void gpu_statevec_anyCtrlTwoTargDenseMatr_sub(Qureg qureg, vector<int> ctrls, ve
     qindex numThreads = qureg.numAmpsPerNode / powerOf2(ctrls.size() + 2);
     qindex numBlocks = getNumBlocks(numThreads);
 
-    devints sortedQubits = util_getSorted(ctrls, {targ1,targ2});
+    vector<int> sortedQubits = util_getSorted(ctrls, {targ1,targ2}); // change for performance
     qindex qubitStateMask = util_getBitMask(ctrls, ctrlStates, {targ1,targ2}, {0,0});
 
+    QubitList_t qubits_dev;
+    std::copy(sortedQubits.begin(), sortedQubits.end(), qubits_dev.indices);
+    qubits_dev.length = sortedQubits.size();
+
     // unpack matrix elems which are more efficiently accessed by kernels as args than shared mem (... maybe...)
     auto m = getFlattenedGpuQcompMatrix<4>(matr.elems); // explicit template for MSVC, grr!
 
     kernel_statevec_anyCtrlTwoTargDenseMatr_sub <NumCtrls> <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
         getGpuQcompPtr(qureg.gpuAmps), numThreads, 
-        getPtr(sortedQubits), ctrls.size(), qubitStateMask, targ1, targ2,
+        qubits_dev, 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]
     );
@@ -432,16 +466,24 @@ void gpu_statevec_anyCtrlAnyTargDenseMatr_sub(Qureg qureg, vector<int> ctrls, ve
     // task each thread with processing more than a single batch
     qindex numBatches = qureg.numAmpsPerNode / powerOf2(ctrls.size() + targs.size());
 
-    devints deviceTargs = targs;
-    devints deviceQubits = util_getSorted(ctrls, targs);
+    vector<int> deviceTargs = targs;
+    vector<int> deviceQubits = util_getSorted(ctrls, targs); // change for performance
     qindex qubitStateMask = util_getBitMask(ctrls, ctrlStates, targs, vector<int>(targs.size(),0));
 
+    QubitList_t qubits_dev;
+    std::copy(deviceQubits.begin(), deviceQubits.end(), qubits_dev.indices);
+    qubits_dev.length = deviceQubits.size();
+
+    QubitList_t targs_dev;
+    std::copy(deviceTargs.begin(), deviceTargs.end(), targs_dev.indices);
+    targs_dev.length = deviceTargs.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
@@ -469,8 +511,8 @@ void gpu_statevec_anyCtrlAnyTargDenseMatr_sub(Qureg qureg, vector<int> ctrls, ve
             <NumCtrls, NumTargs, ApplyConj, ApplyTransp> 
             <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
                 ampsPtr, numThreads, 
-                qubitsPtr, nCtrls, qubitStateMask, 
-                targsPtr, matrPtr
+                qubits_dev, qubitStateMask, 
+                targs_dev, matrPtr
         );
 
     } else {
@@ -508,8 +550,8 @@ void gpu_statevec_anyCtrlAnyTargDenseMatr_sub(Qureg qureg, vector<int> ctrls, ve
             <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
                 getGpuQcompPtr(cache),
                 ampsPtr, numThreads, numBatchesPerThread, 
-                qubitsPtr, nCtrls, qubitStateMask, 
-                targsPtr, targs.size(), powerOf2(targs.size()), matrPtr
+                qubits_dev, qubitStateMask, 
+                targs_dev, powerOf2(targs.size()), matrPtr
         );
     }
 
@@ -568,13 +610,33 @@ void gpu_statevec_anyCtrlOneTargDiagMatr_sub(Qureg qureg, vector<int> ctrls, vec
     qindex numThreads = qureg.numAmpsPerNode / powerOf2(ctrls.size());
     qindex numBlocks = getNumBlocks(numThreads);
 
-    devints deviceCtrls = util_getSorted(ctrls);
+
+    // removed implicit thrust mem copy
+    vector<int> sortedCtrls = util_getSorted(ctrls); // change for performance and standarisation
+
+    ctrl_device_t ctrl;
+
+    std::copy(sortedCtrls.begin(), sortedCtrls.end(), ctrl.ctrl_device);
+
+
+    // Assume size of ctls is at most one per qubit so small enough for device contant memory
+    //int ctrl_device[ctrls.size()];
+
+    // cudaMemcpyToSymbol(ctrl_device, sortedCtrls.data(), ctrls.size()*sizeof(int));
+
+//    cudaMemcpyToSymbol (const char * 	symbol,
+// const void * 	src,
+// size_t 	count,
+// size_t 	offset = 0,
+// enum cudaMemcpyKind 	kind = cudaMemcpyHostToDevice	 
+// )	
+
     qindex ctrlStateMask = util_getBitMask(ctrls, ctrlStates);
     auto elems = getGpuQcompArray<2>(matr.elems); // explicit template for MSVC, grr!
 
     kernel_statevec_anyCtrlOneTargDiagMatr_sub <NumCtrls> <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
-        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, ctrl,
+        ctrls.size(), ctrlStateMask, targ, elems[0], elems[1]
     );
 
     // explicitly return to avoid runtime error below
@@ -636,13 +698,27 @@ void gpu_statevec_anyCtrlTwoTargDiagMatr_sub(Qureg qureg, vector<int> ctrls, vec
     qindex numThreads = qureg.numAmpsPerNode / powerOf2(ctrls.size());
     qindex numBlocks = getNumBlocks(numThreads);
 
-    devints deviceCtrls = util_getSorted(ctrls);
+    // devints deviceCtrls = util_getSorted(ctrls);
+
+    // removed implicit thrust mem copy
+    vector<int> sortedCtrls = util_getSorted(ctrls); // change for performance and standardisationn
+
+    ctrl_device_t ctrl;
+
+    std::copy(sortedCtrls.begin(), sortedCtrls.end(), ctrl.ctrl_device);
+
+    // Assume size of ctls is at most one per qubit so small enough for device contant memory
+    // int ctrl_device[ctrls.size()];
+
+    //cudaMemcpyToSymbol(ctrl_device, sortedCtrls.data(), ctrls.size()*sizeof(int));
+    
+    
     qindex ctrlStateMask = util_getBitMask(ctrls, ctrlStates);
     auto elems = getGpuQcompArray<4>(matr.elems); // explicit template for MSVC, grr!
 
     kernel_statevec_anyCtrlTwoTargDiagMatr_sub <NumCtrls> <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
-        getGpuQcompPtr(qureg.gpuAmps), numThreads, qureg.rank, qureg.logNumAmpsPerNode,
-        getPtr(deviceCtrls), ctrls.size(), ctrlStateMask, targ1, targ2,
+        getGpuQcompPtr(qureg.gpuAmps), numThreads, qureg.rank, qureg.logNumAmpsPerNode, ctrl,
+        ctrls.size(), ctrlStateMask, targ1, targ2,
         elems[0], elems[1], elems[2], elems[3]
     );
 
@@ -704,13 +780,31 @@ void gpu_statevec_anyCtrlAnyTargDiagMatr_sub(Qureg qureg, vector<int> ctrls, vec
     qindex numThreads = qureg.numAmpsPerNode / powerOf2(ctrls.size());
     qindex numBlocks = getNumBlocks(numThreads);
 
-    devints deviceTargs = targs;
-    devints deviceCtrls = util_getSorted(ctrls);
+    vector<int> deviceTargs = targs; // change for performance and standardise 
+    // devints deviceCtrls = util_getSorted(ctrls);
+    
+    QubitList_t targs_dev;
+    std::copy(deviceTargs.begin(), deviceTargs.end(), targs_dev.indices);
+    targs_dev.length = deviceTargs.size();
+
+    // removed implicit thrust mem copy
+    vector<int> sortedCtrls = util_getSorted(ctrls);
+
+    ctrl_device_t ctrl;
+
+    std::copy(sortedCtrls.begin(), sortedCtrls.end(), ctrl.ctrl_device);
+
+    // Assume size of ctls is at most one per qubit so small enough for device contant memory
+    //int ctrl_device[ctrls.size()];
+
+    //cudaMemcpyToSymbol(ctrl_device, sortedCtrls.data(), ctrls.size()*sizeof(int));
+    
+
     qindex ctrlStateMask = util_getBitMask(ctrls, ctrlStates);
 
     kernel_statevec_anyCtrlAnyTargDiagMatr_sub <NumCtrls, NumTargs, ApplyConj, HasPower> <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
-        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, ctrl,
+        ctrls.size(), ctrlStateMask, targs_dev, 
         getGpuQcompPtr(util_getGpuMemPtr(matr)), getGpuQcomp(exponent)
     );
 
@@ -811,10 +905,18 @@ void gpu_statevector_anyCtrlPauliTensorOrGadget_subA(Qureg qureg, vector<int> ct
     auto maskXY  = util_getBitMask(targsXY);
     auto maskYZ  = util_getBitMask(util_getConcatenated(y, z));
 
-    devints deviceTargs   = targsXY;
-    devints deviceQubits  = util_getSorted(ctrls, targsXY);
+    vector<int> deviceTargs   = targsXY; 
+    vector<int> deviceQubits  = util_getSorted(ctrls, targsXY); // change for performance
     qindex qubitStateMask = util_getBitMask(ctrls, ctrlStates, targsXY, vector<int>(targsXY.size(),0));
 
+    QubitList_t qubits_dev;
+    std::copy(deviceQubits.begin(), deviceQubits.end(), qubits_dev.indices);
+    qubits_dev.length = deviceQubits.size();
+
+    QubitList_t targs_dev;
+    std::copy(deviceTargs.begin(), deviceTargs.end(), targs_dev.indices);
+    targs_dev.length = deviceTargs.size();
+
     // unlike the analogous cpu routine, this function has only a single parallelisation
     // granularity; where every pair-of-amps is modified by an independent thread, despite
     // that many threads share a common i0 value (appearing in the kernel). This turns out
@@ -824,8 +926,8 @@ void gpu_statevector_anyCtrlPauliTensorOrGadget_subA(Qureg qureg, vector<int> ct
     qindex numBlocks = getNumBlocks(numThreads);
     kernel_statevector_anyCtrlPauliTensorOrGadget_subA <NumCtrls, NumTargs> <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
         getGpuQcompPtr(qureg.gpuAmps), numThreads,
-        getPtr(deviceQubits), ctrls.size(), qubitStateMask, 
-        getPtr(deviceTargs), deviceTargs.size(),
+        qubits_dev, qubitStateMask, 
+        targs_dev,
         maskXY, maskYZ, getGpuQcomp(powI), getGpuQcomp(ampFac), getGpuQcomp(pairAmpFac)
     );
 
@@ -850,12 +952,16 @@ void gpu_statevector_anyCtrlPauliTensorOrGadget_subB(Qureg qureg, vector<int> ct
     auto maskXY = util_getBitMask(util_getConcatenated(x, y));
     auto maskYZ = util_getBitMask(util_getConcatenated(y, z));
 
-    devints sortedCtrls = util_getSorted(ctrls);
+    vector<int> sortedCtrls = util_getSorted(ctrls); // change for performance
     qindex ctrlStateMask = util_getBitMask(ctrls, ctrlStates);
 
+    QubitList_t Ctrls_dev;
+    std::copy(sortedCtrls.begin(), sortedCtrls.end(), Ctrls_dev.indices);
+    Ctrls_dev.length = sortedCtrls.size();
+
     kernel_statevector_anyCtrlPauliTensorOrGadget_subB <NumCtrls> <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
         getGpuQcompPtr(qureg.gpuAmps), getGpuQcompPtr(qureg.gpuCommBuffer) + recvInd, numThreads, 
-        getPtr(sortedCtrls), ctrls.size(), ctrlStateMask,
+        Ctrls_dev, ctrlStateMask,
         maskXY, maskYZ, bufferMaskXY,
         getGpuQcomp(powI), getGpuQcomp(ampFac), getGpuQcomp(pairAmpFac)
     );
@@ -886,13 +992,17 @@ void gpu_statevector_anyCtrlAnyTargZOrPhaseGadget_sub(Qureg qureg, vector<int> c
     qindex numThreads = qureg.numAmpsPerNode / powerOf2(ctrls.size());
     qindex numBlocks = getNumBlocks(numThreads);
 
-    devints sortedCtrls = util_getSorted(ctrls);
+    vector<int> sortedCtrls = util_getSorted(ctrls); // change for performance
     qindex ctrlStateMask = util_getBitMask(ctrls, ctrlStates);
     qindex targMask = util_getBitMask(targs);
 
+    QubitList_t Ctrls_dev;
+    std::copy(sortedCtrls.begin(), sortedCtrls.end(), Ctrls_dev.indices);
+    Ctrls_dev.length = sortedCtrls.size();
+
     kernel_statevector_anyCtrlAnyTargZOrPhaseGadget_sub <NumCtrls> <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
         getGpuQcompPtr(qureg.gpuAmps), numThreads,
-        getPtr(sortedCtrls), ctrls.size(), ctrlStateMask, targMask,
+        Ctrls_dev, ctrlStateMask, targMask,
         getGpuQcomp(fac0), getGpuQcomp(fac1)
     );
 
@@ -926,8 +1036,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, NUM_THREADS_PER_BLOCK>>> (
         getGpuQcompPtr(outQureg.gpuAmps), numThreads,
@@ -1439,13 +1549,25 @@ void gpu_densmatr_partialTrace_sub(Qureg inQureg, Qureg outQureg, vector<int> ta
     qindex numThreads = outQureg.numAmpsPerNode;
     qindex numBlocks = getNumBlocks(numThreads);
 
-    devints devTargs = targs;
-    devints devPairTargs = pairTargs;
-    devints devAllTargs = util_getSorted(targs, pairTargs);
+    vector<int> devTargs = targs; // change for performance
+    vector<int> devPairTargs = pairTargs;
+    vector<int> devAllTargs = util_getSorted(targs, pairTargs);
+
+    QubitList_t targs_dev;
+    std::copy(devTargs.begin(), devTargs.end(), targs_dev.indices);
+    targs_dev.length = devTargs.size();
+
+    QubitList_t pairTargs_dev;
+    std::copy(devPairTargs.begin(), devPairTargs.end(), pairTargs_dev.indices);
+    pairTargs_dev.length = devPairTargs.size();
+
+    QubitList_t allTargs_dev;
+    std::copy(devAllTargs.begin(), devAllTargs.end(), allTargs_dev.indices);
+    allTargs_dev.length = devAllTargs.size();
 
     kernel_densmatr_partialTrace_sub <NumTargs> <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
         getGpuQcompPtr(inQureg.gpuAmps), getGpuQcompPtr(outQureg.gpuAmps), numThreads,
-        getPtr(devTargs), getPtr(devPairTargs), getPtr(devAllTargs), targs.size()
+        targs_dev, pairTargs_dev, allTargs_dev
     );
 
 #else
@@ -1560,12 +1682,16 @@ void gpu_statevec_calcProbsOfAllMultiQubitOutcomes_sub(qreal* outProbs, Qureg qu
     qindex numBlocks = getNumBlocks(numThreads);
 
     // allocate exponentially-big temporary memory (error if failed)
-    devints devQubits = qubits;
+    vector<int> devQubits = qubits; // change for performance
     devreals devProbs = getDeviceRealsVec(powerOf2(qubits.size())); // throws
 
+    QubitList_t qubits_dev;
+    std::copy(devQubits.begin(), devQubits.end(), qubits_dev.indices);
+    qubits_dev.length = devQubits.size();
+
     kernel_statevec_calcProbsOfAllMultiQubitOutcomes_sub<NumQubits> <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
         getPtr(devProbs), getGpuQcompPtr(qureg.gpuAmps), numThreads, 
-        qureg.rank, qureg.logNumAmpsPerNode, getPtr(devQubits), devQubits.size()
+        qureg.rank, qureg.logNumAmpsPerNode, qubits_dev
     );
 
     // overwrite outProbs with GPU memory
@@ -1597,14 +1723,18 @@ void gpu_densmatr_calcProbsOfAllMultiQubitOutcomes_sub(qreal* outProbs, Qureg qu
     qindex numAmpsPerCol = powerOf2(qureg.numQubits);
 
     // allocate exponentially-big temporary memory (error if failed)
-    devints devQubits = qubits;
+    vector<int> devQubits = qubits; // change for performance
     devreals devProbs = getDeviceRealsVec(powerOf2(qubits.size())); // throws
 
+    QubitList_t qubits_dev;
+    std::copy(devQubits.begin(), devQubits.end(), qubits_dev.indices);
+    qubits_dev.length = devQubits.size();
+
     kernel_densmatr_calcProbsOfAllMultiQubitOutcomes_sub<NumQubits> <<<numBlocks, NUM_THREADS_PER_BLOCK>>> (
         getPtr(devProbs), getGpuQcompPtr(qureg.gpuAmps), 
         numThreads, firstDiagInd, numAmpsPerCol,
         qureg.rank, qureg.logNumAmpsPerNode, 
-        getPtr(devQubits), devQubits.size()
+        qubits_dev
     );
 
     // overwrite outProbs with GPU memory
diff --git a/quest/src/gpu/gpu_thrust.cuh b/quest/src/gpu/gpu_thrust.cuh
index 07a65054..92366056 100644
--- a/quest/src/gpu/gpu_thrust.cuh
+++ b/quest/src/gpu/gpu_thrust.cuh
@@ -68,7 +68,7 @@
  */
 
 
-using devints = thrust::device_vector<int>;
+using devints = thrust::device_vector<int>; // remove for performance
 
 int* getPtr(devints& qubits) {
 
@@ -781,7 +781,7 @@ qreal thrust_densmatr_calcTotalProb_sub(Qureg qureg) {
 template <int NumQubits>
 qreal thrust_statevec_calcProbOfMultiQubitOutcome_sub(Qureg qureg, vector<int> qubits, vector<int> outcomes) {
 
-    devints sortedQubits = util_getSorted(qubits);
+    devints sortedQubits = util_getSorted(qubits); // change for performance
     qindex valueMask = util_getBitMask(qubits, outcomes);
 
     auto indFunctor = functor_insertBits<NumQubits>(getPtr(sortedQubits), valueMask, qubits.size());
@@ -803,7 +803,7 @@ qreal thrust_densmatr_calcProbOfMultiQubitOutcome_sub(Qureg qureg, vector<int> 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 = util_getSorted(qubits);
+    devints sortedQubits = util_getSorted(qubits); // change for performance
     qindex valueMask = util_getBitMask(qubits, outcomes);
 
     auto basisIndFunctor = functor_insertBits<NumQubits>(getPtr(sortedQubits), valueMask, qubits.size());
@@ -1007,7 +1007,7 @@ gpu_qcomp thrust_densmatr_calcExpecFullStateDiagMatr_sub(Qureg qureg, FullStateD
 template <int NumQubits>
 void thrust_statevec_multiQubitProjector_sub(Qureg qureg, vector<int> qubits, vector<int> outcomes, qreal renorm) {
 
-    devints devQubits = qubits;
+    devints devQubits = qubits; // change for performance
     qindex retainValue = getIntegerFromBits(outcomes.data(), outcomes.size());
     auto projFunctor = functor_projectStateVec<NumQubits>(
         getPtr(devQubits), qubits.size(), retainValue, renorm);
@@ -1023,7 +1023,7 @@ void thrust_statevec_multiQubitProjector_sub(Qureg qureg, vector<int> qubits, ve
 template <int NumQubits>
 void thrust_densmatr_multiQubitProjector_sub(Qureg qureg, vector<int> qubits, vector<int> outcomes, qreal renorm) {
 
-    devints devQubits = qubits;
+    devints devQubits = qubits; // change for performance
     qindex retainValue = getIntegerFromBits(outcomes.data(), outcomes.size());
     auto projFunctor = functor_projectDensMatr<NumQubits>(
         getPtr(devQubits), qubits.size(), qureg.rank, qureg.numQubits,
@@ -1078,4 +1078,4 @@ void thrust_statevec_initUnnormalisedUniformlyRandomPureStateAmps_sub(Qureg qure
 
 
 
-#endif // GPU_THRUST_HPP
\ No newline at end of file
+#endif // GPU_THRUST_HPP