Investigate potential single-CPU performance regression in v4

[TysonRayJones]

In v4, the core algorithm of practically all QuEST routines fundamentally changed as per this manuscript. As a comparison, compare this v3.6 routine:

QuEST/QuEST/src/CPU/QuEST_cpu.c

Lines 2196 to 2262 in cc19264

	void statevec_multiControlledUnitaryLocal(
	Qureg qureg, int targetQubit,
	long long int ctrlQubitsMask, long long int ctrlFlipMask,
	ComplexMatrix2 u)
	{
	long long int sizeBlock, sizeHalfBlock;
	long long int thisBlock, // current block
	indexUp,indexLo; // current index and corresponding index in lower half block
	qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
	long long int thisTask;
	long long int numTasks=qureg.numAmpsPerChunk>>1;
	long long int chunkSize=qureg.numAmpsPerChunk;
	long long int chunkId=qureg.chunkId;
	// set dimensions
	sizeHalfBlock = 1LL << targetQubit;
	sizeBlock = 2LL * sizeHalfBlock;
	// Can't use qureg.stateVec as a private OMP var
	qreal *stateVecReal = qureg.stateVec.real;
	qreal *stateVecImag = qureg.stateVec.imag;
	# ifdef _OPENMP
	# pragma omp parallel \
	default (none) \
	shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, u, ctrlQubitsMask,ctrlFlipMask, \
	numTasks,chunkId,chunkSize) \
	private (thisTask,thisBlock, indexUp,indexLo, stateRealUp,stateImagUp,stateRealLo,stateImagLo)
	# endif
	{
	# ifdef _OPENMP
	# pragma omp for schedule (static)
	# endif
	for (thisTask=0; thisTask<numTasks; thisTask++) {
	thisBlock = thisTask / sizeHalfBlock;
	indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
	indexLo = indexUp + sizeHalfBlock;
	// take the basis index, flip the designated (XOR) 'control' bits, AND with the controls.
	// if this equals the control mask, the control qubits have the desired values in the basis index
	if (ctrlQubitsMask == (ctrlQubitsMask & ((indexUp+chunkId*chunkSize) ^ ctrlFlipMask))) {
	// store current state vector values in temp variables
	stateRealUp = stateVecReal[indexUp];
	stateImagUp = stateVecImag[indexUp];
	stateRealLo = stateVecReal[indexLo];
	stateImagLo = stateVecImag[indexLo];
	// state[indexUp] = u00 * state[indexUp] + u01 * state[indexLo]
	stateVecReal[indexUp] = u.real[0][0]*stateRealUp - u.imag[0][0]*stateImagUp
	+ u.real[0][1]*stateRealLo - u.imag[0][1]*stateImagLo;
	stateVecImag[indexUp] = u.real[0][0]*stateImagUp + u.imag[0][0]*stateRealUp
	+ u.real[0][1]*stateImagLo + u.imag[0][1]*stateRealLo;
	// state[indexLo] = u10 * state[indexUp] + u11 * state[indexLo]
	stateVecReal[indexLo] = u.real[1][0]*stateRealUp - u.imag[1][0]*stateImagUp
	+ u.real[1][1]*stateRealLo - u.imag[1][1]*stateImagLo;
	stateVecImag[indexLo] = u.real[1][0]*stateImagUp + u.imag[1][0]*stateRealUp
	+ u.real[1][1]*stateImagLo + u.imag[1][1]*stateRealLo;
	}
	}
	}
	}

to its v4 equivalent:

QuEST/quest/src/cpu/cpu_subroutines.cpp

Lines 361 to 390 in 4d44ec8

	template <int NumCtrls>
	void cpu_statevec_anyCtrlOneTargDenseMatr_subA(Qureg qureg, vector<int> ctrls, vector<int> ctrlStates, int targ, CompMatr1 matr) {
	assert_numCtrlsMatchesNumCtrlStatesAndTemplateParam(ctrls.size(), ctrlStates.size(), NumCtrls);
	// each control qubit halves the needed iterations, and each iteration modifies two amplitudes
	qindex numIts = qureg.numAmpsPerNode / powerOf2(ctrls.size() + 1);
	auto sortedQubits = util_getSorted(ctrls, {targ});
	auto qubitStateMask = util_getBitMask(ctrls, ctrlStates, {targ}, {0});
	// use template param to compile-time unroll loop in insertBits()
	SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrls.size());
	int numQubitBits = numCtrlBits + 1;
	#pragma omp parallel for if(qureg.isMultithreaded)
	for (qindex n=0; n<numIts; n++) {
	// i0 = nth local index where ctrl bits are in specified states and targ is 0
	qindex i0 = insertBitsWithMaskedValues(n, sortedQubits.data(), numQubitBits, qubitStateMask);
	qindex i1 = flipBit(i0, targ);
	// note the two amplitudes are likely strided and not adjacent (separated by 2^t)
	qcomp amp0 = qureg.cpuAmps[i0];
	qcomp amp1 = qureg.cpuAmps[i1];
	qureg.cpuAmps[i0] = matr.elems[0][0]*amp0 + matr.elems[0][1]*amp1;
	qureg.cpuAmps[i1] = matr.elems[1][0]*amp0 + matr.elems[1][1]*amp1;
	}
	}

The v4 optimisations are focused upon avoiding redundant enumeration, avoiding branching within hot loops, and replacing integer arithmetic (which informs updated amplitude indices) with bitwise operations. The code was also cleaned to remove default and superfluous OpenMP boilerplate. As a consequence, the order of accessed amplitudes in the state has substantially changed; from conditional albeit contiguous to unconditional but potentially unpredictable (by a compiler or branch-predictor).

These changes are a benefit to GPU simulation which is detrimentally affected by warp-level looped branching. Alas, despite isolated testing during development, some changed functions may have regressed CPU performance as per #631 (and potentially #616). This might be remedied in multithreaded multi-NUMA-node settings by distributing over nodes, or impending memory-pinning. But in any case, the potential damage to single-CPU performance remains.

TODO: get to the bottom of this! A quick isolated test shows v4 bitwise non-branched enumeration still beats v3's logic (at least on my Apple M1 Max). Figure out which functions have regressed and why!

[TysonRayJones]

Update

A ~5x slowdown seen in serial functions (on my 2017 Macbook) like applyRotateX and applyRotateY - and here's why!

In v3, functions rotateX/Y wrapped the compactUnitary function, passing coefficients $\alpha$, $\beta$, which saw the simulation of matrix:

$$ e^{- \theta/2 i \sigma} \equiv \begin{pmatrix} \alpha & -\beta^* \\ \beta & \alpha^* \end{pmatrix}. $$

The logic resembles

// params of index algebra
long long int numTasks= qureg.numAmpsPerChunk / 2;
qindex sizeHalfBlock = 1LL << targetQubit; 
qindex sizeBlock = 2LL * sizeHalfBlock; 

for (thisTask=0; thisTask<numTasks; thisTask++) {

 // get indices of pair of amplitudes to modify
 qindex thisBlock = thisTask / sizeHalfBlock;
 qindex indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
 qindex indexLo = indexUp + sizeHalfBlock;

 // code which effects...
 state[indexUp] = alpha * state[indexUp] - conj(beta) * state[indexLo];
 state[indexLo] = beta * state[indexUp] + conj(alpha) * state[indexLo];
}

This performs near enough to simulating an arbitrary single-qubit matrix, whereby amplitude pairs are linearly combined via fixed, pre-determined coefficients.

In contrast, v4's applyRotateX/Y invoke a single-qubit edgecase of the multi-qubit applyPauliGadget function. That function uses a totally distinct logic which avoids invoking any matrix representation of the operator, and instead leverages that

$$ \exp(\theta i \hat{\sigma}) \equiv \cos(\theta) \hat{I} + i \sin(\theta) \hat{\sigma} $$

The logic to determine the factor on the pair amplitude, as informed by $\hat{\sigma}$, is somewhat complicated and determined by the parity of several quantities, as explained on page 24 here. Further complicating QuEST v4's code is that it supports rotation around a Pauli string of any size, with any number of control qubits, and where those control qubits are in any arbitrary basis state (rather than all one). Extensive use of templating and compile-time unrolling/inlining ensures (or tries to ensure) that this generality does not worsen the performance versus a bespoke function for the non-controlled single-qubit case when simulated via the same logic.

QuEST v4's full code to do this, containing the many compile-time tricks and optimisations, is:

template <int NumTargs>
INLINE void applyPauliUponAmpPair(
 Qureg qureg, qindex v, qindex i0, int* indXY, int numXY, 
 qindex maskXY, qindex maskYZ, qcomp ampFac, qcomp pairAmpFac
) {
 // this is a subroutine of cpu_statevector_anyCtrlPauliTensorOrGadget_subA() below
 // called in a hot-loop (hence it is here inlined) which exists because the caller
 // chooses one of two possible OpenMP parallelisation granularities

 // remind compiler when NumTargs is compile-time to unroll loop in setBits()
 SET_VAR_AT_COMPILE_TIME(int, numTargBits, NumTargs, numXY);

 // iA = nth local index where targs have value m, iB = (last - nth) such index
 qindex iA = setBits(i0, indXY, numTargBits, v);
 qindex iB = flipBits(iA, maskXY);

 // sign of amps due to Y and Z (excludes Y's i factor which is inside pairAmpFac)
 int signA = fast_getPlusOrMinusMaskedBitParity(iA, maskYZ);
 int signB = fast_getPlusOrMinusMaskedBitParity(iB, maskYZ);

 // mix or swap scaled amp pair (where pairAmpFac includes Y's i factor)
 qcomp ampA = qureg.cpuAmps[iA];
 qcomp ampB = qureg.cpuAmps[iB];
 qureg.cpuAmps[iA] = (ampFac * ampA) + (pairAmpFac * signB * ampB);
 qureg.cpuAmps[iB] = (ampFac * ampB) + (pairAmpFac * signA * ampA);
}


template <int NumCtrls, int NumTargs>
void cpu_statevector_anyCtrlPauliTensorOrGadget_subA(
 Qureg qureg, vector<int> ctrls, vector<int> ctrlStates, 
 vector<int> x, vector<int> y, vector<int> z, qcomp ampFac, qcomp pairAmpFac
) {
 assert_numCtrlsMatchesNumCtrlStatesAndTemplateParam(ctrls.size(), ctrlStates.size(), NumCtrls);
 assert_numTargsMatchesTemplateParam(x.size() + y.size(), NumTargs);
 
 // only X and Y count as targets
 vector<int> sortedTargsXY = util_getSorted(util_getConcatenated(x, y));

 // prepare a mask which yields ctrls in specified state, and X-Y targs in all-zero
 auto sortedQubits = util_getSorted(ctrls, sortedTargsXY);
 auto qubitStateMask = util_getBitMask(ctrls, ctrlStates, sortedTargsXY, vector<int>(sortedTargsXY.size(),0));

 // prepare masks for extracting Pauli parities
 auto maskXY = util_getBitMask(sortedTargsXY);
 auto maskYZ = util_getBitMask(util_getConcatenated(y, z));

 // we will scale pairAmp below by i^numY, so that each amp need only choose the +-1 sign
 pairAmpFac *= util_getPowerOfI(y.size());

 // use template params to compile-time unroll loops in insertBits() and inner-loop below
 SET_VAR_AT_COMPILE_TIME(int, numCtrlBits, NumCtrls, ctrls.size());
 SET_VAR_AT_COMPILE_TIME(int, numTargBits, NumTargs, sortedTargsXY.size());
 int numQubitBits = numCtrlBits + numTargBits;

 // each outer iteration handles all assignments of the target qubits, and each ctrl halves the outer iterations
 qindex numOuterIts = qureg.numAmpsPerNode / powerOf2(numCtrlBits + numTargBits);

 // each inner iteration modifies 2 amplitudes (may be compile-time sized) 
 qindex numInnerIts = powerOf2(numTargBits) / 2; // divides evenly

 // we must choose whether to parallelise the outer or inner loop, since
 // their relative sizes change exponentially depending on numTargAmps.
 // For example, when x+y=O(1), the outer loop is O(2^N) and is parallelised 
 // while the inner loop is O(1). In contrast, when x+y=O(N), the outer loop 
 // is O(1) and left serial while the inner loop is O(2^N) and parallelised.
 // While outer-loop parallelisation can be simply decided with a runtime 
 // if() clause in the pragma, the inner-loop parallelisation must be chosen
 // at compile-time to avoid a huge per-out-loop slowdown when x+y~1. 
 // We will opt to parallelise the outer loop, leaving the inner loop serial,
 // whenever each thread has at least 1 iteration for itself. And of course
 // we serialise both inner and outer loops when qureg multithreading is off.

 if (!qureg.isMultithreaded || numOuterIts >= cpu_getAvailableNumThreads()) {
 
 // parallel
 #pragma omp parallel for if(qureg.isMultithreaded)
 for (qindex n=0; n<numOuterIts; n++) {

 // i0 = nth local index where ctrls are active and targs are all zero
 qindex i0 = insertBitsWithMaskedValues(n, sortedQubits.data(), numQubitBits, qubitStateMask);

 // serial
 for (qindex v=0; v<numInnerIts; v++)
 applyPauliUponAmpPair<NumTargs>(
 qureg, v, i0, sortedTargsXY.data(), numTargBits, maskXY, maskYZ, ampFac, pairAmpFac);
 }

 } else {

 // serial
 for (qindex n=0; n<numOuterIts; n++) {

 // i0 = nth local index where ctrls are active and targs are all zero
 qindex i0 = insertBitsWithMaskedValues(n, sortedQubits.data(), numQubitBits, qubitStateMask);

 // parallel
 #pragma omp parallel for
 for (qindex v=0; v<numInnerIts; v++)
 applyPauliUponAmpPair<NumTargs>(
 qureg, v, i0, sortedTargsXY.data(), numTargBits, maskXY, maskYZ, ampFac, pairAmpFac);
 }
 }
}

Even with all compile-time optimisations working correctly, v4's logic does not reduce to v3's. Instead, v4 must still evaluate the pair amplitude coefficient dynamically, twice invoking fast_getPlusOrMinusMaskedBitParity() in the hot loop. The result is determined only by present Y and Z operators, so when the PauliStr contains only X (like for applyRotateX), the result is known a priori to be one. Worse, that parity function in-turn calls __builtin_parityll() which has vastly varying performance between CPUs! In essence, QuESTv4 does more work per-amplitude than v3's less general method.

This leads to a performance regression only seen in non-controlled single-target rotations. A quick fix is to redirect such functions (applyRotateX and applyRotateY) to wrap applyCompMatr, passing dense 2x2 matrices. This may also be worthwhile to do for any number of control qubits too, though that should be tested. It is expectedly not worthwhile for more qubits (instantiating dense 2^n by 2^n matrices) because of the distributed issues.

A separate optimisation to the anyCtrlPauliTensorOrGadget subroutines is to pass another template parameter which specifies whether any Y or Z are present; if not, we compile-time avoid invoking fast_getPlusOrMinusMaskedBitParity(). Alas this final optimisation is unexpectedly insufficient to restore v3 performance of applyRotateX on my 2017 Macbook. Generally, it seems that the compile-time inlining and unrolling is not working (clang 12.0.0) despite working fine (from memory) on other tested platforms, like an M4 Macbook. That is, the compile-time optimisations which allow us to avoid code duplication are not known to older compilers, crippling performance.

We could make the bitwise functions themselves accept a template parameter, forcing optimisations like loop-unrolling. It seems such an inelegant shame however! 😢

[TysonRayJones]

See PR #682 for an improvement of the one-qubit Pauli performance and relevant testing scripts

[TysonRayJones]

Further testing revealed a remaining Clang-specific performance regression; ~3x when single-threaded and up to 50x when multithreaded! The issue is Clang's slow std::complex arithmetic overloads which include branching for nan and inf handling, disrupting compiler optimisations (like inlining and vectorising). They also cause issues in LLVM's OpenMP implementation.

Remedying this through code requires replacing use of std::complex arithmetic operators with manual calculation (ew!), and unpacking Qureg struct fields before their reference in OpenMP (ew!).

Thankfully we can instead use additional compiler flags, like -ffast-math, though this is dangerous (e.g. invokes associativity and breaks Kahan summation). I instead found on an ARM Mac M4 that...

-ffinite-math-only -fno-signed-zeros -ffp-contract=fast

fully recovers the speedup of manual complex arithmetic, solving the performance regression. The final flag permits contraction of a*b+c into op(a,b,c) when op is a CPU intrinsic. While it can worsen precision in complex arithmetic, -ffp-contract=on is reportedly default since Clang 14.0 (even when worsening performance because op is not an intrinsic?!) and is anyway the default in GCC anyway.

For safety, we should invoke these flags only when the compiler is Clang/LLVM, the CMAKE_BUILD_TYPE is "Release", and only upon the file cpu_subroutines.cpp which contains all use of std::complex arithmetic in hot-loops.

[TysonRayJones]

Testing of PR #684 is being done through comparison to v3 which performed manual complex arithmetic.

v3

#include "QuEST.h"

#include <chrono>
#include <vector>
#include <random>
#include <iostream>
#include <algorithm>

using std::vector;

static std::mt19937 RNG;


QuESTEnv env;


vector<int> getRandomCtrls(int numCtrls, int maxCtrlIndExcl, int excludeTarget) {

	vector<int> ctrls;
	for (int i=0; i<maxCtrlIndExcl; i++)
		if (i != excludeTarget)
			ctrls.push_back(i);

	// shuffle controls only when non-distributed (lazy)
	if (env.numRanks <= 1)
		std::shuffle(ctrls.begin(), ctrls.end(), RNG);

	return vector<int>(ctrls.begin(), ctrls.begin() + numCtrls);
}


template <typename T>
double getFuncDuration(int nQb, T func) {

	double out = 0;

	for (int t=0; t<nQb; t++) {

		auto ctrls = getRandomCtrls(nQb-1, nQb, t);

		// warmup
		for(int r=0; r<5; r++)
			func(ctrls.data(), t);

		auto start = std::chrono::high_resolution_clock::now();

		func(ctrls.data(), t);
		syncQuESTEnv(env);

		auto end = std::chrono::high_resolution_clock::now();
		out += (end - start).count();
	}

	return out;
}


template <typename T1>
void profileFunc(std::string label, int nQb, T1 func1) {

	std::cout << label << "\t" << getFuncDuration(nQb, func1) << std::endl;
}


int main() {

	env = createQuESTEnv();

	// suppress non-root output
	if (env.rank != 0)
		std::cout.setstate(std::ios_base::failbit);

 // prepare random-ctrl-list RNG
 std::random_device cspnrg;
 unsigned seed = cspnrg();
 RNG.seed(seed);

	// operation args (lazily fixed)
	int st[] = {0,1,0,1};
	qreal a = 0.123;
 ComplexMatrix2 m = {
 .real = {{1,2},{3,4}},
 .imag = {{5,6},{7,8}}
 };

	
	for (int nQb : {5, 10, 15, 20, 25, 30}) {

		std::cout << "\n[" << nQb << " qubits]" << std::endl;

		// beware; will use all available accelerations
		Qureg q = createQureg(nQb, env); // createForcedQureg(nQb);


 profileFunc("+", nQb, 
			[&](int* c, int t) { initPlusState(q); });

 profileFunc("u", nQb, 
			[&](int* c, int t) { applyMatrix2(q, t, m); });


		profileFunc("X", nQb, 
			[&](int* c, int t) { pauliX(q, t); });

 profileFunc("cX", nQb, 
			[&](int* c, int t) { multiControlledMultiQubitNot(q, c, 1, &t, 1); });

 profileFunc("ccccX", nQb, 
			[&](int* c, int t) { multiControlledMultiQubitNot(q, c, 4, &t, 1); });


		profileFunc("Y", nQb, 
			[&](int* c, int t) { pauliY(q, t); });

		profileFunc("cY", nQb, 
			[&](int* c, int t) { controlledPauliY(q, c[0], t); });


		profileFunc("Z", nQb, 
			[&](int* c, int t) { pauliZ(q, t); });

		profileFunc("cZ", nQb, 
			[&](int* c, int t) { 
 int all[] = {c[0], t};
 multiControlledPhaseFlip(q, all, 2); });

		profileFunc("ccccZ", nQb, 
			[&](int* c, int t) { 
 int all[] = {c[0], c[1], c[2], c[3], t};
 multiControlledPhaseFlip(q, all, 5); });


		profileFunc("Rx", nQb, 
			[&](int* c, int t) { rotateX(q, t, a); });

		profileFunc("cRx", nQb, 
			[&](int* c, int t) { controlledRotateX(q, c[0], t, a); });


		profileFunc("Ry", nQb, 
			[&](int* c, int t) { rotateY(q, t, a); });

		profileFunc("cRy", nQb, 
			[&](int* c, int t) { controlledRotateY(q, c[0], t, a); });


		profileFunc("Rz", nQb, 
			[&](int* c, int t) { rotateZ(q, t, a); });

		profileFunc("cRz", nQb, 
			[&](int* c, int t) { controlledRotateZ(q, c[0], t, a); });


		destroyQureg(q, env);
	}
	
	destroyQuESTEnv(env);
	return 0;
}

v4

#include "quest.h"

#include <chrono>
#include <vector>
#include <random>
#include <iostream>
#include <algorithm>

using std::vector;

static std::mt19937 RNG;


QuESTEnv env;


vector<int> getRandomCtrls(int numCtrls, int maxCtrlIndExcl, int excludeTarget) {

	vector<int> ctrls;
	for (int i=0; i<maxCtrlIndExcl; i++)
		if (i != excludeTarget)
			ctrls.push_back(i);

	// shuffle controls only when non-distributed (lazy)
	if (!env.isDistributed)
		std::shuffle(ctrls.begin(), ctrls.end(), RNG);

	return vector<int>(ctrls.begin(), ctrls.begin() + numCtrls);
}


template <typename T>
double getFuncDuration(int nQb, T func) {

	double out = 0;

	for (int t=0; t<nQb; t++) {

		auto ctrls = getRandomCtrls(nQb-1, nQb, t);

		// warmup
		for(int r=0; r<5; r++)
			func(ctrls.data(), t);

		auto start = std::chrono::high_resolution_clock::now();

		func(ctrls.data(), t);
		syncQuESTEnv();

		auto end = std::chrono::high_resolution_clock::now();
		out += (end - start).count();
	}

	return out;
}


template <typename T1>
void profileFunc(std::string label, int nQb, T1 func1) {

	std::cout << label << "\t" << getFuncDuration(nQb, func1) << std::endl;
}


int main() {
	initQuESTEnv();
 env = getQuESTEnv();

 setValidationOff();

	// suppress non-root output
	if (env.rank != 0)
		std::cout.setstate(std::ios_base::failbit);

				// // prepare matrix alternatives to Pauli functions
				// CompMatr1 x = getCompMatr1({{0,1}, {1,0}});
				// CompMatr1 y = getCompMatr1({{0,-1i}, {1i,0}});
				// DiagMatr1 z = getDiagMatr1({1, -1});

				// qreal a = 0.123;
				// qreal c = std::cos(-a/2);
				// qreal s = std::sin(-a/2);
				// qcomp v = std::exp(qcomp(0, -a/2));

				// CompMatr1 rx = getCompMatr1({
				// {qcomp(c,0), qcomp(0,s)},
				// {qcomp(0,s), qcomp(c,0)}
				// });
				// CompMatr1 ry = getCompMatr1({
				// { c, s},
				// {-s, c}
				// });
				// DiagMatr1 rz = getDiagMatr1(
				// 	{v, std::conj(v)}
				// );

 // prepare random-ctrl-list RNG
 std::random_device cspnrg;
 unsigned seed = cspnrg();
 RNG.seed(seed);

	// operation args (lazily fixed)
	int st[] = {0,1,0,1};
	qreal a = 0.123;
	CompMatr1 u = getCompMatr1({
		{1+2_i,3+4_i},
		{5+6_i,7+8_i}
	});


	for (int nQb : {5, 10, 15, 20, 25, 30}) {

		std::cout << "\n[" << nQb << " qubits]" << std::endl;

		// beware; will use all available accelerations
		Qureg q = createForcedQureg(nQb);


 profileFunc("+", nQb, 
			[&](int* c, int t) { initPlusState(q); });

 profileFunc("u", nQb, 
			[&](int* c, int t) { applyCompMatr1(q, t, u); });


		profileFunc("X", nQb, 
			[&](int* c, int t) { applyPauliX(q, t); });

		profileFunc("cX", nQb, 
			[&](int* c, int t) { applyControlledPauliX(q, c[0], t); });

		profileFunc("ccccX", nQb, 
			[&](int* c, int t) { applyMultiControlledPauliX(q, c, 4, t); });


		profileFunc("Y", nQb, 
			[&](int* c, int t) { applyPauliY(q, t); });

		profileFunc("cY", nQb, 
			[&](int* c, int t) { applyControlledPauliY(q, c[0], t); });


		profileFunc("Z", nQb, 
			[&](int* c, int t) { applyPauliZ(q, t); });

		profileFunc("cZ", nQb, 
			[&](int* c, int t) { applyControlledPauliZ(q, c[0], t); });

		profileFunc("ccccZ", nQb, 
			[&](int* c, int t) { applyMultiControlledPauliZ(q, c, 4, t); });


		profileFunc("Rx", nQb, 
			[&](int* c, int t) { applyRotateX(q, t, a); });

		profileFunc("cRx", nQb, 
			[&](int* c, int t) { applyControlledRotateX(q, c[0], t, a); });


		profileFunc("Ry", nQb, 
			[&](int* c, int t) { applyRotateY(q, t, a); });

		profileFunc("cRy", nQb, 
			[&](int* c, int t) { applyControlledRotateY(q, c[0], t, a); });


		profileFunc("Rz", nQb, 
			[&](int* c, int t) { applyRotateZ(q, t, a); });

		profileFunc("cRz", nQb, 
			[&](int* c, int t) { applyControlledRotateZ(q, c[0], t, a); });


		destroyQureg(q);
	}
	
	finalizeQuESTEnv();
	return 0;
}

[TysonRayJones]

Patched (although perhaps unsatisfyingly) in 60f9b3a