Macros | Functions
QuEST_gpu.cu File Reference
#include "QuEST.h"
#include "QuEST_precision.h"
#include "QuEST_internal.h"
#include "mt19937ar.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

Go to the source code of this file.

Macros

#define DEBUG   0
 
#define REDUCE_SHARED_SIZE   512
 

Functions

__global__ void copySharedReduceBlock (qreal *arrayIn, qreal *reducedArray, int length)
 
void copyStateFromGPU (Qureg qureg)
 In GPU mode, this copies the state-vector (or density matrix) from GPU memory (qureg.deviceStateVec) to RAM (qureg.stateVec), where it can be accessed/modified by the user. More...
 
void copyStateToGPU (Qureg qureg)
 In GPU mode, this copies the state-vector (or density matrix) from RAM (qureg.stateVec) to VRAM / GPU-memory (qureg.deviceStateVec), which is the version operated upon by other calls to the API. More...
 
QuESTEnv createQuESTEnv (void)
 Create the QuEST execution environment. More...
 
qreal densmatr_calcFidelity (Qureg qureg, Qureg pureState)
 
__global__ void densmatr_calcFidelityKernel (Qureg dens, Qureg vec, long long int dim, qreal *reducedArray)
 computes one term of (vec^*T) dens * vec More...
 
qreal densmatr_calcHilbertSchmidtDistance (Qureg a, Qureg b)
 
__global__ void densmatr_calcHilbertSchmidtDistanceSquaredKernel (qreal *aRe, qreal *aIm, qreal *bRe, qreal *bIm, long long int numAmpsToSum, qreal *reducedArray)
 
qreal densmatr_calcInnerProduct (Qureg a, Qureg b)
 
__global__ void densmatr_calcInnerProductKernel (Qureg a, Qureg b, long long int numTermsToSum, qreal *reducedArray)
 computes Tr(conjTrans(a) b) = sum of (a_ij^* b_ij), which is a real number More...
 
qreal densmatr_calcProbOfOutcome (Qureg qureg, const int measureQubit, int outcome)
 
qreal densmatr_calcPurity (Qureg qureg)
 Computes the trace of the density matrix squared. More...
 
__global__ void densmatr_calcPurityKernel (qreal *vecReal, qreal *vecImag, long long int numAmpsToSum, qreal *reducedArray)
 
qreal densmatr_calcTotalProb (Qureg qureg)
 
void densmatr_collapseToKnownProbOutcome (Qureg qureg, const int measureQubit, int outcome, qreal outcomeProb)
 This involves finding |...i...><...j...| states and killing those where i!=j. More...
 
__global__ void densmatr_collapseToKnownProbOutcomeKernel (qreal outcomeProb, qreal *vecReal, qreal *vecImag, long long int numBasesToVisit, long long int part1, long long int part2, long long int part3, long long int rowBit, long long int colBit, long long int desired, long long int undesired)
 Maps thread ID to a |..0..><..0..| state and then locates |0><1|, |1><0| and |1><1|. More...
 
qreal densmatr_findProbabilityOfZero (Qureg qureg, const int measureQubit)
 
__global__ void densmatr_findProbabilityOfZeroKernel (Qureg qureg, const int measureQubit, qreal *reducedArray)
 
void densmatr_initClassicalState (Qureg qureg, long long int stateInd)
 
__global__ void densmatr_initClassicalStateKernel (long long int densityNumElems, qreal *densityReal, qreal *densityImag, long long int densityInd)
 
void densmatr_initPlusState (Qureg qureg)
 
__global__ void densmatr_initPlusStateKernel (long long int stateVecSize, qreal probFactor, qreal *stateVecReal, qreal *stateVecImag)
 
void densmatr_initPureState (Qureg targetQureg, Qureg copyQureg)
 
__global__ void densmatr_initPureStateKernel (long long int numPureAmps, qreal *targetVecReal, qreal *targetVecImag, qreal *copyVecReal, qreal *copyVecImag)
 
void densmatr_mixDamping (Qureg qureg, const int targetQubit, qreal damping)
 
__global__ void densmatr_mixDampingKernel (qreal damping, qreal *vecReal, qreal *vecImag, long long int numAmpsToVisit, long long int part1, long long int part2, long long int part3, long long int bothBits)
 Works like mixDephasing but modifies every other element, and elements are averaged in pairs. More...
 
void densmatr_mixDensityMatrix (Qureg combineQureg, qreal otherProb, Qureg otherQureg)
 
__global__ void densmatr_mixDensityMatrixKernel (Qureg combineQureg, qreal otherProb, Qureg otherQureg, long long int numAmpsToVisit)
 
void densmatr_mixDephasing (Qureg qureg, const int targetQubit, qreal dephase)
 
__global__ void densmatr_mixDephasingKernel (qreal fac, qreal *vecReal, qreal *vecImag, long long int numAmpsToVisit, long long int part1, long long int part2, long long int part3, long long int colBit, long long int rowBit)
 Called once for every 4 amplitudes in density matrix Works by establishing the |..0..><..0..| state (for its given index) then visiting |..1..><..0..| and |..0..><..1..|. More...
 
void densmatr_mixDepolarising (Qureg qureg, const int targetQubit, qreal depolLevel)
 
__global__ void densmatr_mixDepolarisingKernel (qreal depolLevel, qreal *vecReal, qreal *vecImag, long long int numAmpsToVisit, long long int part1, long long int part2, long long int part3, long long int bothBits)
 Works like mixDephasing but modifies every other element, and elements are averaged in pairs. More...
 
void densmatr_mixTwoQubitDephasing (Qureg qureg, int qubit1, int qubit2, qreal dephase)
 
__global__ void densmatr_mixTwoQubitDephasingKernel (qreal fac, qreal *vecReal, qreal *vecImag, long long int numBackgroundStates, long long int numAmpsToVisit, long long int part1, long long int part2, long long int part3, long long int part4, long long int part5, long long int colBit1, long long int rowBit1, long long int colBit2, long long int rowBit2)
 Called 12 times for every 16 amplitudes in density matrix Each sums from the |..0..0..><..0..0..| index to visit either |..0..0..><..0..1..|, |..0..0..><..1..0..|, |..0..0..><..1..1..|, |..0..1..><..0..0..| etc and so on to |..1..1..><..1..0|. More...
 
void densmatr_mixTwoQubitDepolarising (Qureg qureg, int qubit1, int qubit2, qreal depolLevel)
 
__global__ void densmatr_mixTwoQubitDepolarisingKernel (qreal depolLevel, qreal *vecReal, qreal *vecImag, long long int numAmpsToVisit, long long int part1, long long int part2, long long int part3, long long int part4, long long int part5, long long int rowCol1, long long int rowCol2)
 Called once for every 16 amplitudes. More...
 
void densmatr_oneQubitDegradeOffDiagonal (Qureg qureg, const int targetQubit, qreal dephFac)
 
void destroyQuESTEnv (QuESTEnv env)
 Destroy the QuEST environment. More...
 
__forceinline__ __device__ int extractBit (int locationOfBitFromRight, long long int theEncodedNumber)
 
__forceinline__ __device__ long long int flipBit (long long int number, int bitInd)
 
__forceinline__ __device__ int getBitMaskParity (long long int mask)
 
void getEnvironmentString (QuESTEnv env, Qureg qureg, char str[200])
 Sets str to a string containing the number of qubits in qureg, and the hardware facilities used (e.g. More...
 
int getNumReductionLevels (long long int numValuesToReduce, int numReducedPerLevel)
 
int GPUExists (void)
 
__forceinline__ __device__ long long int insertTwoZeroBits (long long int number, int bit1, int bit2)
 
__forceinline__ __device__ long long int insertZeroBit (long long int number, int index)
 
__forceinline__ __device__ long long int insertZeroBits (long long int number, int *inds, int numInds)
 
__device__ __host__ unsigned int log2Int (unsigned int x)
 
__device__ void reduceBlock (qreal *arrayIn, qreal *reducedArray, int length)
 
void reportQuESTEnv (QuESTEnv env)
 Report information about the QuEST environment. More...
 
void seedQuESTDefault ()
 Seed the Mersenne Twister used for random number generation in the QuEST environment with an example defualt seed. More...
 
Complex statevec_calcInnerProduct (Qureg bra, Qureg ket)
 Terrible code which unnecessarily individually computes and sums the real and imaginary components of the inner product, so as to not have to worry about keeping the sums separated during reduction. More...
 
__global__ void statevec_calcInnerProductKernel (int getRealComp, qreal *vecReal1, qreal *vecImag1, qreal *vecReal2, qreal *vecImag2, long long int numTermsToSum, qreal *reducedArray)
 computes either a real or imag term in the inner product More...
 
qreal statevec_calcProbOfOutcome (Qureg qureg, const int measureQubit, int outcome)
 
qreal statevec_calcTotalProb (Qureg qureg)
 
void statevec_cloneQureg (Qureg targetQureg, Qureg copyQureg)
 works for both statevectors and density matrices More...
 
void statevec_collapseToKnownProbOutcome (Qureg qureg, const int measureQubit, int outcome, qreal outcomeProb)
 
__global__ void statevec_collapseToKnownProbOutcomeKernel (Qureg qureg, int measureQubit, int outcome, qreal totalProbability)
 
void statevec_compactUnitary (Qureg qureg, const int targetQubit, Complex alpha, Complex beta)
 
__global__ void statevec_compactUnitaryKernel (Qureg qureg, const int rotQubit, Complex alpha, Complex beta)
 
int statevec_compareStates (Qureg mq1, Qureg mq2, qreal precision)
 
void statevec_controlledCompactUnitary (Qureg qureg, const int controlQubit, const int targetQubit, Complex alpha, Complex beta)
 
__global__ void statevec_controlledCompactUnitaryKernel (Qureg qureg, const int controlQubit, const int targetQubit, Complex alpha, Complex beta)
 
void statevec_controlledNot (Qureg qureg, const int controlQubit, const int targetQubit)
 
__global__ void statevec_controlledNotKernel (Qureg qureg, const int controlQubit, const int targetQubit)
 
void statevec_controlledPauliY (Qureg qureg, const int controlQubit, const int targetQubit)
 
void statevec_controlledPauliYConj (Qureg qureg, const int controlQubit, const int targetQubit)
 
__global__ void statevec_controlledPauliYKernel (Qureg qureg, const int controlQubit, const int targetQubit, const int conjFac)
 
void statevec_controlledPhaseFlip (Qureg qureg, const int idQubit1, const int idQubit2)
 
__global__ void statevec_controlledPhaseFlipKernel (Qureg qureg, const int idQubit1, const int idQubit2)
 
void statevec_controlledPhaseShift (Qureg qureg, const int idQubit1, const int idQubit2, qreal angle)
 
__global__ void statevec_controlledPhaseShiftKernel (Qureg qureg, const int idQubit1, const int idQubit2, qreal cosAngle, qreal sinAngle)
 
void statevec_controlledUnitary (Qureg qureg, const int controlQubit, const int targetQubit, ComplexMatrix2 u)
 
__global__ void statevec_controlledUnitaryKernel (Qureg qureg, const int controlQubit, const int targetQubit, ArgMatrix2 u)
 
void statevec_createQureg (Qureg *qureg, int numQubits, QuESTEnv env)
 
void statevec_destroyQureg (Qureg qureg, QuESTEnv env)
 
qreal statevec_findProbabilityOfZero (Qureg qureg, const int measureQubit)
 
__global__ void statevec_findProbabilityOfZeroKernel (Qureg qureg, const int measureQubit, qreal *reducedArray)
 
qreal statevec_getImagAmp (Qureg qureg, long long int index)
 
qreal statevec_getRealAmp (Qureg qureg, long long int index)
 
void statevec_hadamard (Qureg qureg, const int targetQubit)
 
__global__ void statevec_hadamardKernel (Qureg qureg, const int targetQubit)
 
void statevec_initBlankState (Qureg qureg)
 
__global__ void statevec_initBlankStateKernel (long long int stateVecSize, qreal *stateVecReal, qreal *stateVecImag)
 
void statevec_initClassicalState (Qureg qureg, long long int stateInd)
 
__global__ void statevec_initClassicalStateKernel (long long int stateVecSize, qreal *stateVecReal, qreal *stateVecImag, long long int stateInd)
 
void statevec_initDebugState (Qureg qureg)
 Initialise the state vector of probability amplitudes to an (unphysical) state with each component of each probability amplitude a unique floating point value. More...
 
__global__ void statevec_initDebugStateKernel (long long int stateVecSize, qreal *stateVecReal, qreal *stateVecImag)
 
void statevec_initPlusState (Qureg qureg)
 
__global__ void statevec_initPlusStateKernel (long long int stateVecSize, qreal *stateVecReal, qreal *stateVecImag)
 
int statevec_initStateFromSingleFile (Qureg *qureg, char filename[200], QuESTEnv env)
 
void statevec_initStateOfSingleQubit (Qureg *qureg, int qubitId, int outcome)
 Initialise the state vector of probability amplitudes such that one qubit is set to 'outcome' and all other qubits are in an equal superposition of zero and one. More...
 
__global__ void statevec_initStateOfSingleQubitKernel (long long int stateVecSize, qreal *stateVecReal, qreal *stateVecImag, int qubitId, int outcome)
 
void statevec_initZeroState (Qureg qureg)
 
__global__ void statevec_initZeroStateKernel (long long int stateVecSize, qreal *stateVecReal, qreal *stateVecImag)
 
void statevec_multiControlledMultiQubitUnitary (Qureg qureg, long long int ctrlMask, int *targs, const int numTargs, ComplexMatrixN u)
 This calls swapQubitAmps only when it would involve a distributed communication; if the qubit chunks already fit in the node, it operates the unitary direct. More...
 
__global__ void statevec_multiControlledMultiQubitUnitaryKernel (Qureg qureg, long long int ctrlMask, int *targs, int numTargs, qreal *uRe, qreal *uIm, long long int *ampInds, qreal *reAmps, qreal *imAmps, long long int numTargAmps)
 
void statevec_multiControlledPhaseFlip (Qureg qureg, int *controlQubits, int numControlQubits)
 
__global__ void statevec_multiControlledPhaseFlipKernel (Qureg qureg, long long int mask)
 
void statevec_multiControlledPhaseShift (Qureg qureg, int *controlQubits, int numControlQubits, qreal angle)
 
__global__ void statevec_multiControlledPhaseShiftKernel (Qureg qureg, long long int mask, qreal cosAngle, qreal sinAngle)
 
void statevec_multiControlledTwoQubitUnitary (Qureg qureg, long long int ctrlMask, const int q1, const int q2, ComplexMatrix4 u)
 This calls swapQubitAmps only when it would involve a distributed communication; if the qubit chunks already fit in the node, it operates the unitary direct. More...
 
__global__ void statevec_multiControlledTwoQubitUnitaryKernel (Qureg qureg, long long int ctrlMask, const int q1, const int q2, ArgMatrix4 u)
 
void statevec_multiControlledUnitary (Qureg qureg, long long int ctrlQubitsMask, long long int ctrlFlipMask, const int targetQubit, ComplexMatrix2 u)
 
__global__ void statevec_multiControlledUnitaryKernel (Qureg qureg, long long int ctrlQubitsMask, long long int ctrlFlipMask, const int targetQubit, ArgMatrix2 u)
 
void statevec_multiRotateZ (Qureg qureg, long long int mask, qreal angle)
 
__global__ void statevec_multiRotateZKernel (Qureg qureg, long long int mask, qreal cosAngle, qreal sinAngle)
 
void statevec_pauliX (Qureg qureg, const int targetQubit)
 
__global__ void statevec_pauliXKernel (Qureg qureg, const int targetQubit)
 
void statevec_pauliY (Qureg qureg, const int targetQubit)
 
void statevec_pauliYConj (Qureg qureg, const int targetQubit)
 
__global__ void statevec_pauliYKernel (Qureg qureg, const int targetQubit, const int conjFac)
 
void statevec_phaseShiftByTerm (Qureg qureg, const int targetQubit, Complex term)
 
__global__ void statevec_phaseShiftByTermKernel (Qureg qureg, const int targetQubit, qreal cosAngle, qreal sinAngle)
 
void statevec_reportStateToScreen (Qureg qureg, QuESTEnv env, int reportRank)
 Print the current state vector of probability amplitudes for a set of qubits to standard out. More...
 
void statevec_setAmps (Qureg qureg, long long int startInd, qreal *reals, qreal *imags, long long int numAmps)
 
void statevec_setWeightedQureg (Complex fac1, Qureg qureg1, Complex fac2, Qureg qureg2, Complex facOut, Qureg out)
 
__global__ void statevec_setWeightedQuregKernel (Complex fac1, Qureg qureg1, Complex fac2, Qureg qureg2, Complex facOut, Qureg out)
 
void statevec_swapQubitAmps (Qureg qureg, int qb1, int qb2)
 
__global__ void statevec_swapQubitAmpsKernel (Qureg qureg, int qb1, int qb2)
 
void statevec_unitary (Qureg qureg, const int targetQubit, ComplexMatrix2 u)
 
__global__ void statevec_unitaryKernel (Qureg qureg, const int targetQubit, ArgMatrix2 u)
 
void swapDouble (qreal **a, qreal **b)
 
void syncQuESTEnv (QuESTEnv env)
 Guarantees that all code up to the given point has been executed on all nodes (if running in distributed mode) More...
 
int syncQuESTSuccess (int successCode)
 Performs a logical AND on all successCodes held by all processes. More...
 

Detailed Description

An implementation of the backend in ../QuEST_internal.h for a GPU environment.

Author
Ania Brown
Tyson Jones

Definition in file QuEST_gpu.cu.

Macro Definition Documentation

◆ DEBUG

#define DEBUG   0

Definition at line 20 of file QuEST_gpu.cu.

◆ REDUCE_SHARED_SIZE

#define REDUCE_SHARED_SIZE   512

Definition at line 19 of file QuEST_gpu.cu.

Function Documentation

◆ copySharedReduceBlock()

__global__ void copySharedReduceBlock ( qreal *  arrayIn,
qreal *  reducedArray,
int  length 
)

Definition at line 1751 of file QuEST_gpu.cu.

1751  {
1752  extern __shared__ qreal tempReductionArray[];
1753  int blockOffset = blockIdx.x*length;
1754  tempReductionArray[threadIdx.x*2] = arrayIn[blockOffset + threadIdx.x*2];
1755  tempReductionArray[threadIdx.x*2+1] = arrayIn[blockOffset + threadIdx.x*2+1];
1756  __syncthreads();
1757  reduceBlock(tempReductionArray, reducedArray, length);
1758 }

References qreal, and reduceBlock().

Referenced by densmatr_calcFidelity(), densmatr_calcHilbertSchmidtDistance(), densmatr_calcInnerProduct(), densmatr_calcPurity(), densmatr_findProbabilityOfZero(), statevec_calcInnerProduct(), and statevec_findProbabilityOfZero().

◆ densmatr_calcFidelity()

qreal densmatr_calcFidelity ( Qureg  qureg,
Qureg  pureState 
)

Definition at line 2194 of file QuEST_gpu.cu.

2194  {
2195 
2196  // we're summing the square of every term in the density matrix
2197  long long int densityDim = 1LL << qureg.numQubitsRepresented;
2198  long long int numValuesToReduce = densityDim;
2199 
2200  int valuesPerCUDABlock, numCUDABlocks, sharedMemSize;
2201  int maxReducedPerLevel = REDUCE_SHARED_SIZE;
2202  int firstTime = 1;
2203 
2204  while (numValuesToReduce > 1) {
2205 
2206  // need less than one CUDA-BLOCK to reduce
2207  if (numValuesToReduce < maxReducedPerLevel) {
2208  valuesPerCUDABlock = numValuesToReduce;
2209  numCUDABlocks = 1;
2210  }
2211  // otherwise use only full CUDA-BLOCKS
2212  else {
2213  valuesPerCUDABlock = maxReducedPerLevel; // constrained by shared memory
2214  numCUDABlocks = ceil((qreal)numValuesToReduce/valuesPerCUDABlock);
2215  }
2216  // dictates size of reduction array
2217  sharedMemSize = valuesPerCUDABlock*sizeof(qreal);
2218 
2219  // spawn threads to sum the probs in each block
2220  // store the reduction in the pureState array
2221  if (firstTime) {
2222  densmatr_calcFidelityKernel<<<numCUDABlocks, valuesPerCUDABlock, sharedMemSize>>>(
2223  qureg, pureState, densityDim, pureState.firstLevelReduction);
2224  firstTime = 0;
2225 
2226  // sum the block probs
2227  } else {
2228  cudaDeviceSynchronize();
2229  copySharedReduceBlock<<<numCUDABlocks, valuesPerCUDABlock/2, sharedMemSize>>>(
2230  pureState.firstLevelReduction,
2231  pureState.secondLevelReduction, valuesPerCUDABlock);
2232  cudaDeviceSynchronize();
2233  swapDouble(&(pureState.firstLevelReduction), &(pureState.secondLevelReduction));
2234  }
2235 
2236  numValuesToReduce = numValuesToReduce/maxReducedPerLevel;
2237  }
2238 
2239  qreal fidelity;
2240  cudaMemcpy(&fidelity, pureState.firstLevelReduction, sizeof(qreal), cudaMemcpyDeviceToHost);
2241  return fidelity;
2242 }

References copySharedReduceBlock(), Qureg::firstLevelReduction, Qureg::numQubitsRepresented, qreal, REDUCE_SHARED_SIZE, Qureg::secondLevelReduction, and swapDouble().

Referenced by calcFidelity().

◆ densmatr_calcFidelityKernel()

__global__ void densmatr_calcFidelityKernel ( Qureg  dens,
Qureg  vec,
long long int  dim,
qreal *  reducedArray 
)

computes one term of (vec^*T) dens * vec

Definition at line 2156 of file QuEST_gpu.cu.

2156  {
2157 
2158  // figure out which density matrix row to consider
2159  long long int col;
2160  long long int row = blockIdx.x*blockDim.x + threadIdx.x;
2161  if (row >= dim) return;
2162 
2163  qreal* densReal = dens.deviceStateVec.real;
2164  qreal* densImag = dens.deviceStateVec.imag;
2165  qreal* vecReal = vec.deviceStateVec.real;
2166  qreal* vecImag = vec.deviceStateVec.imag;
2167 
2168  // compute the row-th element of the product dens*vec
2169  qreal prodReal = 0;
2170  qreal prodImag = 0;
2171  for (col=0LL; col < dim; col++) {
2172  qreal densElemReal = densReal[dim*col + row];
2173  qreal densElemImag = densImag[dim*col + row];
2174 
2175  prodReal += densElemReal*vecReal[col] - densElemImag*vecImag[col];
2176  prodImag += densElemReal*vecImag[col] + densElemImag*vecReal[col];
2177  }
2178 
2179  // multiply with row-th elem of (vec^*)
2180  qreal termReal = prodImag*vecImag[row] + prodReal*vecReal[row];
2181 
2182  // imag of every term should be zero, because each is a valid fidelity calc of an eigenstate
2183  //qreal termImag = prodImag*vecReal[row] - prodReal*vecImag[row];
2184 
2185  extern __shared__ qreal tempReductionArray[];
2186  tempReductionArray[threadIdx.x] = termReal;
2187  __syncthreads();
2188 
2189  // every second thread reduces
2190  if (threadIdx.x<blockDim.x/2)
2191  reduceBlock(tempReductionArray, reducedArray, blockDim.x);
2192 }

References Qureg::deviceStateVec, qreal, and reduceBlock().

◆ densmatr_calcHilbertSchmidtDistance()

qreal densmatr_calcHilbertSchmidtDistance ( Qureg  a,
Qureg  b 
)

Definition at line 2268 of file QuEST_gpu.cu.

2268  {
2269 
2270  // we're summing the square of every term in (a-b)
2271  long long int numValuesToReduce = a.numAmpsPerChunk;
2272 
2273  int valuesPerCUDABlock, numCUDABlocks, sharedMemSize;
2274  int maxReducedPerLevel = REDUCE_SHARED_SIZE;
2275  int firstTime = 1;
2276 
2277  while (numValuesToReduce > 1) {
2278 
2279  // need less than one CUDA-BLOCK to reduce
2280  if (numValuesToReduce < maxReducedPerLevel) {
2281  valuesPerCUDABlock = numValuesToReduce;
2282  numCUDABlocks = 1;
2283  }
2284  // otherwise use only full CUDA-BLOCKS
2285  else {
2286  valuesPerCUDABlock = maxReducedPerLevel; // constrained by shared memory
2287  numCUDABlocks = ceil((qreal)numValuesToReduce/valuesPerCUDABlock);
2288  }
2289  // dictates size of reduction array
2290  sharedMemSize = valuesPerCUDABlock*sizeof(qreal);
2291 
2292  // spawn threads to sum the probs in each block (store reduction temp values in a's reduction array)
2293  if (firstTime) {
2294  densmatr_calcHilbertSchmidtDistanceSquaredKernel<<<numCUDABlocks, valuesPerCUDABlock, sharedMemSize>>>(
2295  a.deviceStateVec.real, a.deviceStateVec.imag,
2296  b.deviceStateVec.real, b.deviceStateVec.imag,
2297  numValuesToReduce, a.firstLevelReduction);
2298  firstTime = 0;
2299 
2300  // sum the block probs
2301  } else {
2302  cudaDeviceSynchronize();
2303  copySharedReduceBlock<<<numCUDABlocks, valuesPerCUDABlock/2, sharedMemSize>>>(
2304  a.firstLevelReduction,
2305  a.secondLevelReduction, valuesPerCUDABlock);
2306  cudaDeviceSynchronize();
2307  swapDouble(&(a.firstLevelReduction), &(a.secondLevelReduction));
2308  }
2309 
2310  numValuesToReduce = numValuesToReduce/maxReducedPerLevel;
2311  }
2312 
2313  qreal trace;
2314  cudaMemcpy(&trace, a.firstLevelReduction, sizeof(qreal), cudaMemcpyDeviceToHost);
2315 
2316  qreal sqrtTrace = sqrt(trace);
2317  return sqrtTrace;
2318 }

References copySharedReduceBlock(), Qureg::deviceStateVec, Qureg::firstLevelReduction, Qureg::numAmpsPerChunk, qreal, REDUCE_SHARED_SIZE, Qureg::secondLevelReduction, and swapDouble().

Referenced by calcHilbertSchmidtDistance().

◆ densmatr_calcHilbertSchmidtDistanceSquaredKernel()

__global__ void densmatr_calcHilbertSchmidtDistanceSquaredKernel ( qreal *  aRe,
qreal *  aIm,
qreal *  bRe,
qreal *  bIm,
long long int  numAmpsToSum,
qreal *  reducedArray 
)

Definition at line 2244 of file QuEST_gpu.cu.

2247  {
2248  // figure out which density matrix term this thread is assigned
2249  long long int index = blockIdx.x*blockDim.x + threadIdx.x;
2250  if (index >= numAmpsToSum) return;
2251 
2252  // compute this thread's sum term
2253  qreal difRe = aRe[index] - bRe[index];
2254  qreal difIm = aIm[index] - bIm[index];
2255  qreal term = difRe*difRe + difIm*difIm;
2256 
2257  // array of each thread's collected term, to be summed
2258  extern __shared__ qreal tempReductionArray[];
2259  tempReductionArray[threadIdx.x] = term;
2260  __syncthreads();
2261 
2262  // every second thread reduces
2263  if (threadIdx.x<blockDim.x/2)
2264  reduceBlock(tempReductionArray, reducedArray, blockDim.x);
2265 }

References qreal, and reduceBlock().

◆ densmatr_calcInnerProduct()

qreal densmatr_calcInnerProduct ( Qureg  a,
Qureg  b 
)

Definition at line 1988 of file QuEST_gpu.cu.

1988  {
1989 
1990  // we're summing the square of every term in the density matrix
1991  long long int numValuesToReduce = a.numAmpsTotal;
1992 
1993  int valuesPerCUDABlock, numCUDABlocks, sharedMemSize;
1994  int maxReducedPerLevel = REDUCE_SHARED_SIZE;
1995  int firstTime = 1;
1996 
1997  while (numValuesToReduce > 1) {
1998 
1999  // need less than one CUDA-BLOCK to reduce
2000  if (numValuesToReduce < maxReducedPerLevel) {
2001  valuesPerCUDABlock = numValuesToReduce;
2002  numCUDABlocks = 1;
2003  }
2004  // otherwise use only full CUDA-BLOCKS
2005  else {
2006  valuesPerCUDABlock = maxReducedPerLevel; // constrained by shared memory
2007  numCUDABlocks = ceil((qreal)numValuesToReduce/valuesPerCUDABlock);
2008  }
2009  // dictates size of reduction array
2010  sharedMemSize = valuesPerCUDABlock*sizeof(qreal);
2011 
2012  // spawn threads to sum the terms in each block
2013  // arbitrarily store the reduction in the b qureg's array
2014  if (firstTime) {
2015  densmatr_calcInnerProductKernel<<<numCUDABlocks, valuesPerCUDABlock, sharedMemSize>>>(
2016  a, b, a.numAmpsTotal, b.firstLevelReduction);
2017  firstTime = 0;
2018  }
2019  // sum the block terms
2020  else {
2021  cudaDeviceSynchronize();
2022  copySharedReduceBlock<<<numCUDABlocks, valuesPerCUDABlock/2, sharedMemSize>>>(
2023  b.firstLevelReduction,
2024  b.secondLevelReduction, valuesPerCUDABlock);
2025  cudaDeviceSynchronize();
2026  swapDouble(&(b.firstLevelReduction), &(b.secondLevelReduction));
2027  }
2028 
2029  numValuesToReduce = numValuesToReduce/maxReducedPerLevel;
2030  }
2031 
2032  qreal innerprod;
2033  cudaMemcpy(&innerprod, b.firstLevelReduction, sizeof(qreal), cudaMemcpyDeviceToHost);
2034  return innerprod;
2035 }

References copySharedReduceBlock(), Qureg::firstLevelReduction, Qureg::numAmpsTotal, qreal, REDUCE_SHARED_SIZE, Qureg::secondLevelReduction, and swapDouble().

Referenced by calcDensityInnerProduct().

◆ densmatr_calcInnerProductKernel()

__global__ void densmatr_calcInnerProductKernel ( Qureg  a,
Qureg  b,
long long int  numTermsToSum,
qreal *  reducedArray 
)

computes Tr(conjTrans(a) b) = sum of (a_ij^* b_ij), which is a real number

Definition at line 1967 of file QuEST_gpu.cu.

1969  {
1970  long long int index = blockIdx.x*blockDim.x + threadIdx.x;
1971  if (index >= numTermsToSum) return;
1972 
1973  // Re{ conj(a) b } = Re{ (aRe - i aIm)(bRe + i bIm) } = aRe bRe + aIm bIm
1974  qreal prod = (
1975  a.deviceStateVec.real[index]*b.deviceStateVec.real[index]
1976  + a.deviceStateVec.imag[index]*b.deviceStateVec.imag[index]);
1977 
1978  // array of each thread's collected sum term, to be summed
1979  extern __shared__ qreal tempReductionArray[];
1980  tempReductionArray[threadIdx.x] = prod;
1981  __syncthreads();
1982 
1983  // every second thread reduces
1984  if (threadIdx.x<blockDim.x/2)
1985  reduceBlock(tempReductionArray, reducedArray, blockDim.x);
1986 }

References Qureg::deviceStateVec, qreal, and reduceBlock().

◆ densmatr_calcProbOfOutcome()

qreal densmatr_calcProbOfOutcome ( Qureg  qureg,
const int  measureQubit,
int  outcome 
)

Definition at line 1958 of file QuEST_gpu.cu.

1959 {
1960  qreal outcomeProb = densmatr_findProbabilityOfZero(qureg, measureQubit);
1961  if (outcome==1)
1962  outcomeProb = 1.0 - outcomeProb;
1963  return outcomeProb;
1964 }

References densmatr_findProbabilityOfZero(), and qreal.

Referenced by calcProbOfOutcome(), collapseToOutcome(), and densmatr_measureWithStats().

◆ densmatr_calcPurity()

qreal densmatr_calcPurity ( Qureg  qureg)

Computes the trace of the density matrix squared.

Definition at line 2339 of file QuEST_gpu.cu.

2339  {
2340 
2341  // we're summing the square of every term in the density matrix
2342  long long int numValuesToReduce = qureg.numAmpsPerChunk;
2343 
2344  int valuesPerCUDABlock, numCUDABlocks, sharedMemSize;
2345  int maxReducedPerLevel = REDUCE_SHARED_SIZE;
2346  int firstTime = 1;
2347 
2348  while (numValuesToReduce > 1) {
2349 
2350  // need less than one CUDA-BLOCK to reduce
2351  if (numValuesToReduce < maxReducedPerLevel) {
2352  valuesPerCUDABlock = numValuesToReduce;
2353  numCUDABlocks = 1;
2354  }
2355  // otherwise use only full CUDA-BLOCKS
2356  else {
2357  valuesPerCUDABlock = maxReducedPerLevel; // constrained by shared memory
2358  numCUDABlocks = ceil((qreal)numValuesToReduce/valuesPerCUDABlock);
2359  }
2360  // dictates size of reduction array
2361  sharedMemSize = valuesPerCUDABlock*sizeof(qreal);
2362 
2363  // spawn threads to sum the probs in each block
2364  if (firstTime) {
2365  densmatr_calcPurityKernel<<<numCUDABlocks, valuesPerCUDABlock, sharedMemSize>>>(
2366  qureg.deviceStateVec.real, qureg.deviceStateVec.imag,
2367  numValuesToReduce, qureg.firstLevelReduction);
2368  firstTime = 0;
2369 
2370  // sum the block probs
2371  } else {
2372  cudaDeviceSynchronize();
2373  copySharedReduceBlock<<<numCUDABlocks, valuesPerCUDABlock/2, sharedMemSize>>>(
2374  qureg.firstLevelReduction,
2375  qureg.secondLevelReduction, valuesPerCUDABlock);
2376  cudaDeviceSynchronize();
2377  swapDouble(&(qureg.firstLevelReduction), &(qureg.secondLevelReduction));
2378  }
2379 
2380  numValuesToReduce = numValuesToReduce/maxReducedPerLevel;
2381  }
2382 
2383  qreal traceDensSquared;
2384  cudaMemcpy(&traceDensSquared, qureg.firstLevelReduction, sizeof(qreal), cudaMemcpyDeviceToHost);
2385  return traceDensSquared;
2386 }

References copySharedReduceBlock(), Qureg::deviceStateVec, Qureg::firstLevelReduction, Qureg::numAmpsPerChunk, qreal, REDUCE_SHARED_SIZE, Qureg::secondLevelReduction, and swapDouble().

Referenced by calcPurity().

◆ densmatr_calcPurityKernel()

__global__ void densmatr_calcPurityKernel ( qreal *  vecReal,
qreal *  vecImag,
long long int  numAmpsToSum,
qreal *  reducedArray 
)

Definition at line 2320 of file QuEST_gpu.cu.

2320  {
2321 
2322  // figure out which density matrix term this thread is assigned
2323  long long int index = blockIdx.x*blockDim.x + threadIdx.x;
2324  if (index >= numAmpsToSum) return;
2325 
2326  qreal term = vecReal[index]*vecReal[index] + vecImag[index]*vecImag[index];
2327 
2328  // array of each thread's collected probability, to be summed
2329  extern __shared__ qreal tempReductionArray[];
2330  tempReductionArray[threadIdx.x] = term;
2331  __syncthreads();
2332 
2333  // every second thread reduces
2334  if (threadIdx.x<blockDim.x/2)
2335  reduceBlock(tempReductionArray, reducedArray, blockDim.x);
2336 }

References qreal, and reduceBlock().

◆ densmatr_calcTotalProb()

qreal densmatr_calcTotalProb ( Qureg  qureg)

Definition at line 1476 of file QuEST_gpu.cu.

1476  {
1477 
1478  // computes the trace using Kahan summation
1479  qreal pTotal=0;
1480  qreal y, t, c;
1481  c = 0;
1482 
1483  long long int numCols = 1LL << qureg.numQubitsRepresented;
1484  long long diagIndex;
1485 
1486  copyStateFromGPU(qureg);
1487 
1488  for (int col=0; col< numCols; col++) {
1489  diagIndex = col*(numCols + 1);
1490  y = qureg.stateVec.real[diagIndex] - c;
1491  t = pTotal + y;
1492  c = ( t - pTotal ) - y; // brackets are important
1493  pTotal = t;
1494  }
1495 
1496  return pTotal;
1497 }

References copyStateFromGPU(), Qureg::numQubitsRepresented, qreal, and Qureg::stateVec.

Referenced by calcTotalProb(), and statevec_calcExpecPauliProd().

◆ densmatr_collapseToKnownProbOutcome()

void densmatr_collapseToKnownProbOutcome ( Qureg  qureg,
const int  measureQubit,
int  outcome,
qreal  outcomeProb 
)

This involves finding |...i...><...j...| states and killing those where i!=j.

Renorms (/prob) every | * outcome * >< * outcome * | state, setting all others to zero.

Definition at line 2480 of file QuEST_gpu.cu.

2480  {
2481 
2482  int rowQubit = measureQubit + qureg.numQubitsRepresented;
2483 
2484  int colBit = 1LL << measureQubit;
2485  int rowBit = 1LL << rowQubit;
2486 
2487  long long int numBasesToVisit = qureg.numAmpsPerChunk/4;
2488  long long int part1 = colBit -1;
2489  long long int part2 = (rowBit >> 1) - colBit;
2490  long long int part3 = numBasesToVisit - (rowBit >> 1);
2491 
2492  long long int desired, undesired;
2493  if (outcome == 0) {
2494  desired = 0;
2495  undesired = colBit | rowBit;
2496  } else {
2497  desired = colBit | rowBit;
2498  undesired = 0;
2499  }
2500 
2501  int threadsPerCUDABlock, CUDABlocks;
2502  threadsPerCUDABlock = 128;
2503  CUDABlocks = ceil(numBasesToVisit / (qreal) threadsPerCUDABlock);
2504  densmatr_collapseToKnownProbOutcomeKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
2505  outcomeProb, qureg.deviceStateVec.real, qureg.deviceStateVec.imag, numBasesToVisit,
2506  part1, part2, part3, rowBit, colBit, desired, undesired);
2507 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, and qreal.

Referenced by collapseToOutcome(), and densmatr_measureWithStats().

◆ densmatr_collapseToKnownProbOutcomeKernel()

__global__ void densmatr_collapseToKnownProbOutcomeKernel ( qreal  outcomeProb,
qreal *  vecReal,
qreal *  vecImag,
long long int  numBasesToVisit,
long long int  part1,
long long int  part2,
long long int  part3,
long long int  rowBit,
long long int  colBit,
long long int  desired,
long long int  undesired 
)

Maps thread ID to a |..0..><..0..| state and then locates |0><1|, |1><0| and |1><1|.

Definition at line 2454 of file QuEST_gpu.cu.

2458 {
2459  long long int scanInd = blockIdx.x*blockDim.x + threadIdx.x;
2460  if (scanInd >= numBasesToVisit) return;
2461 
2462  long long int base = (scanInd&part1) + ((scanInd&part2)<<1) + ((scanInd&part3)<<2);
2463 
2464  // renormalise desired outcome
2465  vecReal[base + desired] /= outcomeProb;
2466  vecImag[base + desired] /= outcomeProb;
2467 
2468  // kill undesired outcome
2469  vecReal[base + undesired] = 0;
2470  vecImag[base + undesired] = 0;
2471 
2472  // kill |..0..><..1..| states
2473  vecReal[base + colBit] = 0;
2474  vecImag[base + colBit] = 0;
2475  vecReal[base + rowBit] = 0;
2476  vecImag[base + rowBit] = 0;
2477 }

◆ densmatr_findProbabilityOfZero()

qreal densmatr_findProbabilityOfZero ( Qureg  qureg,
const int  measureQubit 
)

Definition at line 1864 of file QuEST_gpu.cu.

1865 {
1866  long long int densityDim = 1LL << qureg.numQubitsRepresented;
1867  long long int numValuesToReduce = densityDim >> 1; // half of the diagonal has measureQubit=0
1868 
1869  int valuesPerCUDABlock, numCUDABlocks, sharedMemSize;
1870  int maxReducedPerLevel = REDUCE_SHARED_SIZE;
1871  int firstTime = 1;
1872 
1873  while (numValuesToReduce > 1) {
1874 
1875  // need less than one CUDA-BLOCK to reduce
1876  if (numValuesToReduce < maxReducedPerLevel) {
1877  valuesPerCUDABlock = numValuesToReduce;
1878  numCUDABlocks = 1;
1879  }
1880  // otherwise use only full CUDA-BLOCKS
1881  else {
1882  valuesPerCUDABlock = maxReducedPerLevel; // constrained by shared memory
1883  numCUDABlocks = ceil((qreal)numValuesToReduce/valuesPerCUDABlock);
1884  }
1885 
1886  sharedMemSize = valuesPerCUDABlock*sizeof(qreal);
1887 
1888  // spawn threads to sum the probs in each block
1889  if (firstTime) {
1890  densmatr_findProbabilityOfZeroKernel<<<numCUDABlocks, valuesPerCUDABlock, sharedMemSize>>>(
1891  qureg, measureQubit, qureg.firstLevelReduction);
1892  firstTime = 0;
1893 
1894  // sum the block probs
1895  } else {
1896  cudaDeviceSynchronize();
1897  copySharedReduceBlock<<<numCUDABlocks, valuesPerCUDABlock/2, sharedMemSize>>>(
1898  qureg.firstLevelReduction,
1899  qureg.secondLevelReduction, valuesPerCUDABlock);
1900  cudaDeviceSynchronize();
1901  swapDouble(&(qureg.firstLevelReduction), &(qureg.secondLevelReduction));
1902  }
1903 
1904  numValuesToReduce = numValuesToReduce/maxReducedPerLevel;
1905  }
1906 
1907  qreal zeroProb;
1908  cudaMemcpy(&zeroProb, qureg.firstLevelReduction, sizeof(qreal), cudaMemcpyDeviceToHost);
1909  return zeroProb;
1910 }

References copySharedReduceBlock(), Qureg::firstLevelReduction, Qureg::numQubitsRepresented, qreal, REDUCE_SHARED_SIZE, Qureg::secondLevelReduction, and swapDouble().

Referenced by densmatr_calcProbOfOutcome().

◆ densmatr_findProbabilityOfZeroKernel()

__global__ void densmatr_findProbabilityOfZeroKernel ( Qureg  qureg,
const int  measureQubit,
qreal *  reducedArray 
)

Definition at line 1760 of file QuEST_gpu.cu.

1762  {
1763  // run by each thread
1764  // use of block here refers to contiguous amplitudes where measureQubit = 0,
1765  // (then =1) and NOT the CUDA block, which is the partitioning of CUDA threads
1766 
1767  long long int densityDim = 1LL << qureg.numQubitsRepresented;
1768  long long int numTasks = densityDim >> 1;
1769  long long int sizeHalfBlock = 1LL << (measureQubit);
1770  long long int sizeBlock = 2LL * sizeHalfBlock;
1771 
1772  long long int thisBlock; // which block this thread is processing
1773  long long int thisTask; // which part of the block this thread is processing
1774  long long int basisIndex; // index of this thread's computational basis state
1775  long long int densityIndex; // " " index of |basis><basis| in the flat density matrix
1776 
1777  // array of each thread's collected probability, to be summed
1778  extern __shared__ qreal tempReductionArray[];
1779 
1780  // figure out which density matrix prob that this thread is assigned
1781  thisTask = blockIdx.x*blockDim.x + threadIdx.x;
1782  if (thisTask>=numTasks) return;
1783  thisBlock = thisTask / sizeHalfBlock;
1784  basisIndex = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
1785  densityIndex = (densityDim + 1) * basisIndex;
1786 
1787  // record the probability in the CUDA-BLOCK-wide array
1788  qreal prob = qureg.deviceStateVec.real[densityIndex]; // im[densityIndex] assumed ~ 0
1789  tempReductionArray[threadIdx.x] = prob;
1790 
1791  // sum the probs collected by this CUDA-BLOCK's threads into a per-CUDA-BLOCK array
1792  __syncthreads();
1793  if (threadIdx.x<blockDim.x/2){
1794  reduceBlock(tempReductionArray, reducedArray, blockDim.x);
1795  }
1796 }

References Qureg::deviceStateVec, Qureg::numQubitsRepresented, qreal, and reduceBlock().

◆ densmatr_initClassicalState()

void densmatr_initClassicalState ( Qureg  qureg,
long long int  stateInd 
)

Definition at line 258 of file QuEST_gpu.cu.

259 {
260  int threadsPerCUDABlock, CUDABlocks;
261  threadsPerCUDABlock = 128;
262  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk)/threadsPerCUDABlock);
263 
264  // index of the desired state in the flat density matrix
265  long long int densityDim = 1LL << qureg.numQubitsRepresented;
266  long long int densityInd = (densityDim + 1)*stateInd;
267 
268  // identical to pure version
269  densmatr_initClassicalStateKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
270  qureg.numAmpsPerChunk,
271  qureg.deviceStateVec.real,
272  qureg.deviceStateVec.imag, densityInd);
273 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, and qreal.

Referenced by initClassicalState().

◆ densmatr_initClassicalStateKernel()

__global__ void densmatr_initClassicalStateKernel ( long long int  densityNumElems,
qreal *  densityReal,
qreal *  densityImag,
long long int  densityInd 
)

Definition at line 239 of file QuEST_gpu.cu.

243 {
244  // initialise the state to all zeros
245  long long int index = blockIdx.x*blockDim.x + threadIdx.x;
246  if (index >= densityNumElems) return;
247 
248  densityReal[index] = 0.0;
249  densityImag[index] = 0.0;
250 
251  if (index==densityInd){
252  // classical state has probability 1
253  densityReal[densityInd] = 1.0;
254  densityImag[densityInd] = 0.0;
255  }
256 }

◆ densmatr_initPlusState()

void densmatr_initPlusState ( Qureg  qureg)

Definition at line 226 of file QuEST_gpu.cu.

227 {
228  qreal probFactor = 1.0/((qreal) (1LL << qureg.numQubitsRepresented));
229  int threadsPerCUDABlock, CUDABlocks;
230  threadsPerCUDABlock = 128;
231  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk)/threadsPerCUDABlock);
232  densmatr_initPlusStateKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
233  qureg.numAmpsPerChunk,
234  probFactor,
235  qureg.deviceStateVec.real,
236  qureg.deviceStateVec.imag);
237 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, and qreal.

Referenced by initPlusState().

◆ densmatr_initPlusStateKernel()

__global__ void densmatr_initPlusStateKernel ( long long int  stateVecSize,
qreal  probFactor,
qreal *  stateVecReal,
qreal *  stateVecImag 
)

Definition at line 216 of file QuEST_gpu.cu.

216  {
217  long long int index;
218 
219  index = blockIdx.x*blockDim.x + threadIdx.x;
220  if (index>=stateVecSize) return;
221 
222  stateVecReal[index] = probFactor;
223  stateVecImag[index] = 0.0;
224 }

◆ densmatr_initPureState()

void densmatr_initPureState ( Qureg  targetQureg,
Qureg  copyQureg 
)

Definition at line 205 of file QuEST_gpu.cu.

206 {
207  int threadsPerCUDABlock, CUDABlocks;
208  threadsPerCUDABlock = 128;
209  CUDABlocks = ceil((qreal)(copyQureg.numAmpsPerChunk)/threadsPerCUDABlock);
210  densmatr_initPureStateKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
211  copyQureg.numAmpsPerChunk,
212  targetQureg.deviceStateVec.real, targetQureg.deviceStateVec.imag,
213  copyQureg.deviceStateVec.real, copyQureg.deviceStateVec.imag);
214 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

Referenced by initPureState().

◆ densmatr_initPureStateKernel()

__global__ void densmatr_initPureStateKernel ( long long int  numPureAmps,
qreal *  targetVecReal,
qreal *  targetVecImag,
qreal *  copyVecReal,
qreal *  copyVecImag 
)

Definition at line 186 of file QuEST_gpu.cu.

190 {
191  // this is a particular index of the pure copyQureg
192  long long int index = blockIdx.x*blockDim.x + threadIdx.x;
193  if (index>=numPureAmps) return;
194 
195  qreal realRow = copyVecReal[index];
196  qreal imagRow = copyVecImag[index];
197  for (long long int col=0; col < numPureAmps; col++) {
198  qreal realCol = copyVecReal[col];
199  qreal imagCol = - copyVecImag[col]; // minus for conjugation
200  targetVecReal[col*numPureAmps + index] = realRow*realCol - imagRow*imagCol;
201  targetVecImag[col*numPureAmps + index] = realRow*imagCol + imagRow*realCol;
202  }
203 }

References qreal.

◆ densmatr_mixDamping()

void densmatr_mixDamping ( Qureg  qureg,
const int  targetQubit,
qreal  damping 
)

Definition at line 2723 of file QuEST_gpu.cu.

2723  {
2724 
2725  if (damping == 0)
2726  return;
2727 
2728  qreal dephase = sqrt(1-damping);
2729  densmatr_oneQubitDegradeOffDiagonal(qureg, targetQubit, dephase);
2730 
2731  long long int numAmpsToVisit = qureg.numAmpsPerChunk/4;
2732  int rowQubit = targetQubit + qureg.numQubitsRepresented;
2733 
2734  long long int colBit = 1LL << targetQubit;
2735  long long int rowBit = 1LL << rowQubit;
2736  long long int bothBits = colBit | rowBit;
2737 
2738  long long int part1 = colBit - 1;
2739  long long int part2 = (rowBit >> 1) - colBit;
2740  long long int part3 = numAmpsToVisit - (rowBit >> 1);
2741 
2742  int threadsPerCUDABlock, CUDABlocks;
2743  threadsPerCUDABlock = 128;
2744  CUDABlocks = ceil(numAmpsToVisit / (qreal) threadsPerCUDABlock);
2745  densmatr_mixDampingKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
2746  damping, qureg.deviceStateVec.real, qureg.deviceStateVec.imag, numAmpsToVisit,
2747  part1, part2, part3, bothBits);
2748 }

References densmatr_oneQubitDegradeOffDiagonal(), Qureg::deviceStateVec, Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, and qreal.

Referenced by mixDamping().

◆ densmatr_mixDampingKernel()

__global__ void densmatr_mixDampingKernel ( qreal  damping,
qreal *  vecReal,
qreal *  vecImag,
long long int  numAmpsToVisit,
long long int  part1,
long long int  part2,
long long int  part3,
long long int  bothBits 
)

Works like mixDephasing but modifies every other element, and elements are averaged in pairs.

Definition at line 2676 of file QuEST_gpu.cu.

2680 {
2681  long long int scanInd = blockIdx.x*blockDim.x + threadIdx.x;
2682  if (scanInd >= numAmpsToVisit) return;
2683 
2684  long long int baseInd = (scanInd&part1) + ((scanInd&part2)<<1) + ((scanInd&part3)<<2);
2685  long long int targetInd = baseInd + bothBits;
2686 
2687  qreal realAvDepol = damping * ( vecReal[targetInd]);
2688  qreal imagAvDepol = damping * ( vecImag[targetInd]);
2689 
2690  vecReal[targetInd] *= 1 - damping;
2691  vecImag[targetInd] *= 1 - damping;
2692 
2693  vecReal[baseInd] += realAvDepol;
2694  vecImag[baseInd] += imagAvDepol;
2695 }

References qreal.

◆ densmatr_mixDensityMatrix()

void densmatr_mixDensityMatrix ( Qureg  combineQureg,
qreal  otherProb,
Qureg  otherQureg 
)

Definition at line 2521 of file QuEST_gpu.cu.

2521  {
2522 
2523  long long int numAmpsToVisit = combineQureg.numAmpsPerChunk;
2524 
2525  int threadsPerCUDABlock, CUDABlocks;
2526  threadsPerCUDABlock = 128;
2527  CUDABlocks = ceil(numAmpsToVisit / (qreal) threadsPerCUDABlock);
2528  densmatr_mixDensityMatrixKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
2529  combineQureg, otherProb, otherQureg, numAmpsToVisit
2530  );
2531 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by mixDensityMatrix().

◆ densmatr_mixDensityMatrixKernel()

__global__ void densmatr_mixDensityMatrixKernel ( Qureg  combineQureg,
qreal  otherProb,
Qureg  otherQureg,
long long int  numAmpsToVisit 
)

Definition at line 2509 of file QuEST_gpu.cu.

2509  {
2510 
2511  long long int ampInd = blockIdx.x*blockDim.x + threadIdx.x;
2512  if (ampInd >= numAmpsToVisit) return;
2513 
2514  combineQureg.deviceStateVec.real[ampInd] *= 1-otherProb;
2515  combineQureg.deviceStateVec.imag[ampInd] *= 1-otherProb;
2516 
2517  combineQureg.deviceStateVec.real[ampInd] += otherProb*otherQureg.deviceStateVec.real[ampInd];
2518  combineQureg.deviceStateVec.imag[ampInd] += otherProb*otherQureg.deviceStateVec.imag[ampInd];
2519 }

References Qureg::deviceStateVec.

◆ densmatr_mixDephasing()

void densmatr_mixDephasing ( Qureg  qureg,
const int  targetQubit,
qreal  dephase 
)

Definition at line 2574 of file QuEST_gpu.cu.

2574  {
2575 
2576  if (dephase == 0)
2577  return;
2578 
2579  qreal dephFac = 1 - dephase;
2580  densmatr_oneQubitDegradeOffDiagonal(qureg, targetQubit, dephFac);
2581 }

References densmatr_oneQubitDegradeOffDiagonal(), and qreal.

Referenced by densmatr_mixDepolarising(), and mixDephasing().

◆ densmatr_mixDephasingKernel()

__global__ void densmatr_mixDephasingKernel ( qreal  fac,
qreal *  vecReal,
qreal *  vecImag,
long long int  numAmpsToVisit,
long long int  part1,
long long int  part2,
long long int  part3,
long long int  colBit,
long long int  rowBit 
)

Called once for every 4 amplitudes in density matrix Works by establishing the |..0..><..0..| state (for its given index) then visiting |..1..><..0..| and |..0..><..1..|.

Labels |part1 X pa><rt2 NOT(X) part3| From the brain of Simon Benjamin

Definition at line 2538 of file QuEST_gpu.cu.

2542 {
2543  long long int scanInd = blockIdx.x*blockDim.x + threadIdx.x;
2544  if (scanInd >= numAmpsToVisit) return;
2545 
2546  long long int ampInd = (scanInd&part1) + ((scanInd&part2)<<1) + ((scanInd&part3)<<2);
2547  vecReal[ampInd + colBit] *= fac;
2548  vecImag[ampInd + colBit] *= fac;
2549  vecReal[ampInd + rowBit] *= fac;
2550  vecImag[ampInd + rowBit] *= fac;
2551 }

◆ densmatr_mixDepolarising()

void densmatr_mixDepolarising ( Qureg  qureg,
const int  targetQubit,
qreal  depolLevel 
)

Definition at line 2697 of file QuEST_gpu.cu.

2697  {
2698 
2699  if (depolLevel == 0)
2700  return;
2701 
2702  densmatr_mixDephasing(qureg, targetQubit, depolLevel);
2703 
2704  long long int numAmpsToVisit = qureg.numAmpsPerChunk/4;
2705  int rowQubit = targetQubit + qureg.numQubitsRepresented;
2706 
2707  long long int colBit = 1LL << targetQubit;
2708  long long int rowBit = 1LL << rowQubit;
2709  long long int bothBits = colBit | rowBit;
2710 
2711  long long int part1 = colBit - 1;
2712  long long int part2 = (rowBit >> 1) - colBit;
2713  long long int part3 = numAmpsToVisit - (rowBit >> 1);
2714 
2715  int threadsPerCUDABlock, CUDABlocks;
2716  threadsPerCUDABlock = 128;
2717  CUDABlocks = ceil(numAmpsToVisit / (qreal) threadsPerCUDABlock);
2718  densmatr_mixDepolarisingKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
2719  depolLevel, qureg.deviceStateVec.real, qureg.deviceStateVec.imag, numAmpsToVisit,
2720  part1, part2, part3, bothBits);
2721 }

References densmatr_mixDephasing(), Qureg::deviceStateVec, Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, and qreal.

Referenced by mixDepolarising().

◆ densmatr_mixDepolarisingKernel()

__global__ void densmatr_mixDepolarisingKernel ( qreal  depolLevel,
qreal *  vecReal,
qreal *  vecImag,
long long int  numAmpsToVisit,
long long int  part1,
long long int  part2,
long long int  part3,
long long int  bothBits 
)

Works like mixDephasing but modifies every other element, and elements are averaged in pairs.

Definition at line 2650 of file QuEST_gpu.cu.

2654 {
2655  long long int scanInd = blockIdx.x*blockDim.x + threadIdx.x;
2656  if (scanInd >= numAmpsToVisit) return;
2657 
2658  long long int baseInd = (scanInd&part1) + ((scanInd&part2)<<1) + ((scanInd&part3)<<2);
2659  long long int targetInd = baseInd + bothBits;
2660 
2661  qreal realAvDepol = depolLevel * 0.5 * (vecReal[baseInd] + vecReal[targetInd]);
2662  qreal imagAvDepol = depolLevel * 0.5 * (vecImag[baseInd] + vecImag[targetInd]);
2663 
2664  vecReal[baseInd] *= 1 - depolLevel;
2665  vecImag[baseInd] *= 1 - depolLevel;
2666  vecReal[targetInd] *= 1 - depolLevel;
2667  vecImag[targetInd] *= 1 - depolLevel;
2668 
2669  vecReal[baseInd] += realAvDepol;
2670  vecImag[baseInd] += imagAvDepol;
2671  vecReal[targetInd] += realAvDepol;
2672  vecImag[targetInd] += imagAvDepol;
2673 }

References qreal.

◆ densmatr_mixTwoQubitDephasing()

void densmatr_mixTwoQubitDephasing ( Qureg  qureg,
int  qubit1,
int  qubit2,
qreal  dephase 
)

Definition at line 2613 of file QuEST_gpu.cu.

2613  {
2614 
2615  if (dephase == 0)
2616  return;
2617 
2618  // assumes qubit2 > qubit1
2619 
2620  int rowQubit1 = qubit1 + qureg.numQubitsRepresented;
2621  int rowQubit2 = qubit2 + qureg.numQubitsRepresented;
2622 
2623  long long int colBit1 = 1LL << qubit1;
2624  long long int rowBit1 = 1LL << rowQubit1;
2625  long long int colBit2 = 1LL << qubit2;
2626  long long int rowBit2 = 1LL << rowQubit2;
2627 
2628  long long int part1 = colBit1 - 1;
2629  long long int part2 = (colBit2 >> 1) - colBit1;
2630  long long int part3 = (rowBit1 >> 2) - (colBit2 >> 1);
2631  long long int part4 = (rowBit2 >> 3) - (rowBit1 >> 2);
2632  long long int part5 = (qureg.numAmpsPerChunk/16) - (rowBit2 >> 3);
2633  qreal dephFac = 1 - dephase;
2634 
2635  // refers to states |a 0 b 0 c><d 0 e 0 f| (target qubits are fixed)
2636  long long int numBackgroundStates = qureg.numAmpsPerChunk/16;
2637 
2638  // 12 of these states experience dephasing
2639  long long int numAmpsToVisit = 12 * numBackgroundStates;
2640 
2641  int threadsPerCUDABlock, CUDABlocks;
2642  threadsPerCUDABlock = 128;
2643  CUDABlocks = ceil(numAmpsToVisit / (qreal) threadsPerCUDABlock);
2644  densmatr_mixTwoQubitDephasingKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
2645  dephFac, qureg.deviceStateVec.real, qureg.deviceStateVec.imag, numBackgroundStates, numAmpsToVisit,
2646  part1, part2, part3, part4, part5, colBit1, rowBit1, colBit2, rowBit2);
2647 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, and qreal.

Referenced by densmatr_mixTwoQubitDepolarising(), and mixTwoQubitDephasing().

◆ densmatr_mixTwoQubitDephasingKernel()

__global__ void densmatr_mixTwoQubitDephasingKernel ( qreal  fac,
qreal *  vecReal,
qreal *  vecImag,
long long int  numBackgroundStates,
long long int  numAmpsToVisit,
long long int  part1,
long long int  part2,
long long int  part3,
long long int  part4,
long long int  part5,
long long int  colBit1,
long long int  rowBit1,
long long int  colBit2,
long long int  rowBit2 
)

Called 12 times for every 16 amplitudes in density matrix Each sums from the |..0..0..><..0..0..| index to visit either |..0..0..><..0..1..|, |..0..0..><..1..0..|, |..0..0..><..1..1..|, |..0..1..><..0..0..| etc and so on to |..1..1..><..1..0|.

Labels |part1 0 part2 0 par><t3 0 part4 0 part5|. From the brain of Simon Benjamin

Definition at line 2589 of file QuEST_gpu.cu.

2593 {
2594  long long int outerInd = blockIdx.x*blockDim.x + threadIdx.x;
2595  if (outerInd >= numAmpsToVisit) return;
2596 
2597  // sets meta in 1...14 excluding 5, 10, creating bit string DCBA for |..D..C..><..B..A|
2598  int meta = 1 + (outerInd/numBackgroundStates);
2599  if (meta > 4) meta++;
2600  if (meta > 9) meta++;
2601 
2602  long long int shift = rowBit2*((meta>>3)%2) + rowBit1*((meta>>2)%2) + colBit2*((meta>>1)%2) + colBit1*(meta%2);
2603  long long int scanInd = outerInd % numBackgroundStates;
2604  long long int stateInd = (
2605  shift +
2606  (scanInd&part1) + ((scanInd&part2)<<1) + ((scanInd&part3)<<2) + ((scanInd&part4)<<3) + ((scanInd&part5)<<4));
2607 
2608  vecReal[stateInd] *= fac;
2609  vecImag[stateInd] *= fac;
2610 }

◆ densmatr_mixTwoQubitDepolarising()

void densmatr_mixTwoQubitDepolarising ( Qureg  qureg,
int  qubit1,
int  qubit2,
qreal  depolLevel 
)

Definition at line 2783 of file QuEST_gpu.cu.

2783  {
2784 
2785  if (depolLevel == 0)
2786  return;
2787 
2788  // assumes qubit2 > qubit1
2789 
2790  densmatr_mixTwoQubitDephasing(qureg, qubit1, qubit2, depolLevel);
2791 
2792  int rowQubit1 = qubit1 + qureg.numQubitsRepresented;
2793  int rowQubit2 = qubit2 + qureg.numQubitsRepresented;
2794 
2795  long long int colBit1 = 1LL << qubit1;
2796  long long int rowBit1 = 1LL << rowQubit1;
2797  long long int colBit2 = 1LL << qubit2;
2798  long long int rowBit2 = 1LL << rowQubit2;
2799 
2800  long long int rowCol1 = colBit1 | rowBit1;
2801  long long int rowCol2 = colBit2 | rowBit2;
2802 
2803  long long int numAmpsToVisit = qureg.numAmpsPerChunk/16;
2804  long long int part1 = colBit1 - 1;
2805  long long int part2 = (colBit2 >> 1) - colBit1;
2806  long long int part3 = (rowBit1 >> 2) - (colBit2 >> 1);
2807  long long int part4 = (rowBit2 >> 3) - (rowBit1 >> 2);
2808  long long int part5 = numAmpsToVisit - (rowBit2 >> 3);
2809 
2810  int threadsPerCUDABlock, CUDABlocks;
2811  threadsPerCUDABlock = 128;
2812  CUDABlocks = ceil(numAmpsToVisit / (qreal) threadsPerCUDABlock);
2813  densmatr_mixTwoQubitDepolarisingKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
2814  depolLevel, qureg.deviceStateVec.real, qureg.deviceStateVec.imag, numAmpsToVisit,
2815  part1, part2, part3, part4, part5, rowCol1, rowCol2);
2816 }

References densmatr_mixTwoQubitDephasing(), Qureg::deviceStateVec, Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, and qreal.

Referenced by mixTwoQubitDepolarising().

◆ densmatr_mixTwoQubitDepolarisingKernel()

__global__ void densmatr_mixTwoQubitDepolarisingKernel ( qreal  depolLevel,
qreal *  vecReal,
qreal *  vecImag,
long long int  numAmpsToVisit,
long long int  part1,
long long int  part2,
long long int  part3,
long long int  part4,
long long int  part5,
long long int  rowCol1,
long long int  rowCol2 
)

Called once for every 16 amplitudes.

Definition at line 2751 of file QuEST_gpu.cu.

2756 {
2757  long long int scanInd = blockIdx.x*blockDim.x + threadIdx.x;
2758  if (scanInd >= numAmpsToVisit) return;
2759 
2760  // index of |..0..0..><..0..0|
2761  long long int ind00 = (scanInd&part1) + ((scanInd&part2)<<1) + ((scanInd&part3)<<2) + ((scanInd&part4)<<3) + ((scanInd&part5)<<4);
2762  long long int ind01 = ind00 + rowCol1;
2763  long long int ind10 = ind00 + rowCol2;
2764  long long int ind11 = ind00 + rowCol1 + rowCol2;
2765 
2766  qreal realAvDepol = depolLevel * 0.25 * (
2767  vecReal[ind00] + vecReal[ind01] + vecReal[ind10] + vecReal[ind11]);
2768  qreal imagAvDepol = depolLevel * 0.25 * (
2769  vecImag[ind00] + vecImag[ind01] + vecImag[ind10] + vecImag[ind11]);
2770 
2771  qreal retain = 1 - depolLevel;
2772  vecReal[ind00] *= retain; vecImag[ind00] *= retain;
2773  vecReal[ind01] *= retain; vecImag[ind01] *= retain;
2774  vecReal[ind10] *= retain; vecImag[ind10] *= retain;
2775  vecReal[ind11] *= retain; vecImag[ind11] *= retain;
2776 
2777  vecReal[ind00] += realAvDepol; vecImag[ind00] += imagAvDepol;
2778  vecReal[ind01] += realAvDepol; vecImag[ind01] += imagAvDepol;
2779  vecReal[ind10] += realAvDepol; vecImag[ind10] += imagAvDepol;
2780  vecReal[ind11] += realAvDepol; vecImag[ind11] += imagAvDepol;
2781 }

References qreal.

◆ densmatr_oneQubitDegradeOffDiagonal()

void densmatr_oneQubitDegradeOffDiagonal ( Qureg  qureg,
const int  targetQubit,
qreal  dephFac 
)

Definition at line 2554 of file QuEST_gpu.cu.

2554  {
2555 
2556  long long int numAmpsToVisit = qureg.numAmpsPerChunk/4;
2557 
2558  int rowQubit = targetQubit + qureg.numQubitsRepresented;
2559  long long int colBit = 1LL << targetQubit;
2560  long long int rowBit = 1LL << rowQubit;
2561 
2562  long long int part1 = colBit - 1;
2563  long long int part2 = (rowBit >> 1) - colBit;
2564  long long int part3 = numAmpsToVisit - (rowBit >> 1);
2565 
2566  int threadsPerCUDABlock, CUDABlocks;
2567  threadsPerCUDABlock = 128;
2568  CUDABlocks = ceil(numAmpsToVisit / (qreal) threadsPerCUDABlock);
2569  densmatr_mixDephasingKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
2570  dephFac, qureg.deviceStateVec.real, qureg.deviceStateVec.imag, numAmpsToVisit,
2571  part1, part2, part3, colBit, rowBit);
2572 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, and qreal.

Referenced by densmatr_mixDamping(), and densmatr_mixDephasing().

◆ extractBit()

__forceinline__ __device__ int extractBit ( int  locationOfBitFromRight,
long long int  theEncodedNumber 
)

Definition at line 82 of file QuEST_gpu.cu.

82  {
83  return (theEncodedNumber & ( 1LL << locationOfBitFromRight )) >> locationOfBitFromRight;
84 }

Referenced by compressPairVectorForSingleQubitDepolarise(), compressPairVectorForTwoQubitDepolarise(), densmatr_collapseToKnownProbOutcome(), densmatr_findProbabilityOfZeroLocal(), densmatr_mixDampingDistributed(), densmatr_mixDepolarisingDistributed(), densmatr_mixTwoQubitDepolarisingDistributed(), densmatr_mixTwoQubitDepolarisingQ1LocalQ2DistributedPart3(), getGlobalIndOfOddParityInChunk(), statevec_controlledCompactUnitaryDistributed(), statevec_controlledCompactUnitaryKernel(), statevec_controlledCompactUnitaryLocal(), statevec_controlledNotDistributed(), statevec_controlledNotKernel(), statevec_controlledNotLocal(), statevec_controlledPauliYDistributed(), statevec_controlledPauliYKernel(), statevec_controlledPauliYLocal(), statevec_controlledPhaseFlip(), statevec_controlledPhaseFlipKernel(), statevec_controlledPhaseShift(), statevec_controlledPhaseShiftKernel(), statevec_controlledUnitaryDistributed(), statevec_controlledUnitaryKernel(), statevec_controlledUnitaryLocal(), statevec_initStateOfSingleQubit(), statevec_initStateOfSingleQubitKernel(), statevec_multiControlledMultiQubitUnitaryKernel(), statevec_multiControlledMultiQubitUnitaryLocal(), and statevec_phaseShiftByTerm().

◆ flipBit()

__forceinline__ __device__ long long int flipBit ( long long int  number,
int  bitInd 
)

Definition at line 95 of file QuEST_gpu.cu.

95  {
96  return (number ^ (1LL << bitInd));
97 }

Referenced by getGlobalIndOfOddParityInChunk(), statevec_multiControlledMultiQubitUnitary(), statevec_multiControlledMultiQubitUnitaryKernel(), statevec_multiControlledMultiQubitUnitaryLocal(), statevec_multiControlledTwoQubitUnitary(), statevec_multiControlledTwoQubitUnitaryKernel(), statevec_multiControlledTwoQubitUnitaryLocal(), statevec_swapQubitAmps(), statevec_swapQubitAmpsDistributed(), statevec_swapQubitAmpsKernel(), and statevec_swapQubitAmpsLocal().

◆ getBitMaskParity()

__forceinline__ __device__ int getBitMaskParity ( long long int  mask)

Definition at line 86 of file QuEST_gpu.cu.

86  {
87  int parity = 0;
88  while (mask) {
89  parity = !parity;
90  mask = mask & (mask-1);
91  }
92  return parity;
93 }

Referenced by statevec_multiRotateZKernel().

◆ getNumReductionLevels()

int getNumReductionLevels ( long long int  numValuesToReduce,
int  numReducedPerLevel 
)

Definition at line 1848 of file QuEST_gpu.cu.

1848  {
1849  int levels=0;
1850  while (numValuesToReduce){
1851  numValuesToReduce = numValuesToReduce/numReducedPerLevel;
1852  levels++;
1853  }
1854  return levels;
1855 }

◆ GPUExists()

int GPUExists ( void  )

Definition at line 336 of file QuEST_gpu.cu.

336  {
337  int deviceCount, device;
338  int gpuDeviceCount = 0;
339  struct cudaDeviceProp properties;
340  cudaError_t cudaResultCode = cudaGetDeviceCount(&deviceCount);
341  if (cudaResultCode != cudaSuccess) deviceCount = 0;
342  /* machines with no GPUs can still report one emulation device */
343  for (device = 0; device < deviceCount; ++device) {
344  cudaGetDeviceProperties(&properties, device);
345  if (properties.major != 9999) { /* 9999 means emulation only */
346  ++gpuDeviceCount;
347  }
348  }
349  if (gpuDeviceCount) return 1;
350  else return 0;
351 }

Referenced by createQuESTEnv().

◆ insertTwoZeroBits()

__forceinline__ __device__ long long int insertTwoZeroBits ( long long int  number,
int  bit1,
int  bit2 
)

Definition at line 106 of file QuEST_gpu.cu.

106  {
107  int small = (bit1 < bit2)? bit1 : bit2;
108  int big = (bit1 < bit2)? bit2 : bit1;
109  return insertZeroBit(insertZeroBit(number, small), big);
110 }

References insertZeroBit().

Referenced by statevec_multiControlledTwoQubitUnitaryKernel(), statevec_multiControlledTwoQubitUnitaryLocal(), statevec_swapQubitAmpsKernel(), and statevec_swapQubitAmpsLocal().

◆ insertZeroBit()

__forceinline__ __device__ long long int insertZeroBit ( long long int  number,
int  index 
)

Definition at line 99 of file QuEST_gpu.cu.

99  {
100  long long int left, right;
101  left = (number >> index) << index;
102  right = number - left;
103  return (left << 1) ^ right;
104 }

Referenced by insertTwoZeroBits(), insertZeroBits(), and statevec_multiControlledMultiQubitUnitaryLocal().

◆ insertZeroBits()

__forceinline__ __device__ long long int insertZeroBits ( long long int  number,
int *  inds,
int  numInds 
)

Definition at line 112 of file QuEST_gpu.cu.

112  {
113  /* inserted bit inds must strictly increase, so that their final indices are correct.
114  * in-lieu of sorting (avoided since no C++ variable-size arrays, and since we're already
115  * memory bottle-necked so overhead eats this slowdown), we find the next-smallest index each
116  * at each insert. recall every element of inds (a positive or zero number) is unique.
117  * This function won't appear in the CPU code, which can use C99 variable-size arrays and
118  * ought to make a sorted array before threading
119  */
120  int curMin = inds[0];
121  int prevMin = -1;
122  for (int n=0; n < numInds; n++) {
123 
124  // find next min
125  for (int t=0; t < numInds; t++)
126  if (inds[t]>prevMin && inds[t]<curMin)
127  curMin = inds[t];
128 
129  number = insertZeroBit(number, curMin);
130 
131  // set curMin to an arbitrary non-visited elem
132  prevMin = curMin;
133  for (int t=0; t < numInds; t++)
134  if (inds[t] > curMin) {
135  curMin = inds[t];
136  break;
137  }
138  }
139  return number;
140 }

References insertZeroBit().

Referenced by statevec_multiControlledMultiQubitUnitaryKernel().

◆ log2Int()

__device__ __host__ unsigned int log2Int ( unsigned int  x)

Definition at line 1725 of file QuEST_gpu.cu.

1726 {
1727  unsigned int ans = 0 ;
1728  while( x>>=1 ) ans++;
1729  return ans ;
1730 }

Referenced by reduceBlock().

◆ reduceBlock()

__device__ void reduceBlock ( qreal *  arrayIn,
qreal *  reducedArray,
int  length 
)

Definition at line 1732 of file QuEST_gpu.cu.

1732  {
1733  int i, l, r;
1734  int threadMax, maxDepth;
1735  threadMax = length/2;
1736  maxDepth = log2Int(length/2);
1737 
1738  for (i=0; i<maxDepth+1; i++){
1739  if (threadIdx.x<threadMax){
1740  l = threadIdx.x;
1741  r = l + threadMax;
1742  arrayIn[l] = arrayIn[r] + arrayIn[l];
1743  }
1744  threadMax = threadMax >> 1;
1745  __syncthreads(); // optimise -- use warp shuffle instead
1746  }
1747 
1748  if (threadIdx.x==0) reducedArray[blockIdx.x] = arrayIn[0];
1749 }

References log2Int().

Referenced by copySharedReduceBlock(), densmatr_calcFidelityKernel(), densmatr_calcHilbertSchmidtDistanceSquaredKernel(), densmatr_calcInnerProductKernel(), densmatr_calcPurityKernel(), densmatr_findProbabilityOfZeroKernel(), statevec_calcInnerProductKernel(), and statevec_findProbabilityOfZeroKernel().

◆ statevec_calcInnerProduct()

Complex statevec_calcInnerProduct ( Qureg  bra,
Qureg  ket 
)

Terrible code which unnecessarily individually computes and sums the real and imaginary components of the inner product, so as to not have to worry about keeping the sums separated during reduction.

Truly disgusting, probably doubles runtime, please fix. @TODO could even do the kernel twice, storing real in bra.reduc and imag in ket.reduc?

Definition at line 2068 of file QuEST_gpu.cu.

2068  {
2069 
2070  qreal innerProdReal, innerProdImag;
2071 
2072  int getRealComp;
2073  long long int numValuesToReduce;
2074  int valuesPerCUDABlock, numCUDABlocks, sharedMemSize;
2075  int maxReducedPerLevel;
2076  int firstTime;
2077 
2078  // compute real component of inner product
2079  getRealComp = 1;
2080  numValuesToReduce = bra.numAmpsPerChunk;
2081  maxReducedPerLevel = REDUCE_SHARED_SIZE;
2082  firstTime = 1;
2083  while (numValuesToReduce > 1) {
2084  if (numValuesToReduce < maxReducedPerLevel) {
2085  valuesPerCUDABlock = numValuesToReduce;
2086  numCUDABlocks = 1;
2087  }
2088  else {
2089  valuesPerCUDABlock = maxReducedPerLevel;
2090  numCUDABlocks = ceil((qreal)numValuesToReduce/valuesPerCUDABlock);
2091  }
2092  sharedMemSize = valuesPerCUDABlock*sizeof(qreal);
2093  if (firstTime) {
2094  statevec_calcInnerProductKernel<<<numCUDABlocks, valuesPerCUDABlock, sharedMemSize>>>(
2095  getRealComp,
2096  bra.deviceStateVec.real, bra.deviceStateVec.imag,
2097  ket.deviceStateVec.real, ket.deviceStateVec.imag,
2098  numValuesToReduce,
2099  bra.firstLevelReduction);
2100  firstTime = 0;
2101  } else {
2102  cudaDeviceSynchronize();
2103  copySharedReduceBlock<<<numCUDABlocks, valuesPerCUDABlock/2, sharedMemSize>>>(
2104  bra.firstLevelReduction,
2105  bra.secondLevelReduction, valuesPerCUDABlock);
2106  cudaDeviceSynchronize();
2107  swapDouble(&(bra.firstLevelReduction), &(bra.secondLevelReduction));
2108  }
2109  numValuesToReduce = numValuesToReduce/maxReducedPerLevel;
2110  }
2111  cudaMemcpy(&innerProdReal, bra.firstLevelReduction, sizeof(qreal), cudaMemcpyDeviceToHost);
2112 
2113  // compute imag component of inner product
2114  getRealComp = 0;
2115  numValuesToReduce = bra.numAmpsPerChunk;
2116  maxReducedPerLevel = REDUCE_SHARED_SIZE;
2117  firstTime = 1;
2118  while (numValuesToReduce > 1) {
2119  if (numValuesToReduce < maxReducedPerLevel) {
2120  valuesPerCUDABlock = numValuesToReduce;
2121  numCUDABlocks = 1;
2122  }
2123  else {
2124  valuesPerCUDABlock = maxReducedPerLevel;
2125  numCUDABlocks = ceil((qreal)numValuesToReduce/valuesPerCUDABlock);
2126  }
2127  sharedMemSize = valuesPerCUDABlock*sizeof(qreal);
2128  if (firstTime) {
2129  statevec_calcInnerProductKernel<<<numCUDABlocks, valuesPerCUDABlock, sharedMemSize>>>(
2130  getRealComp,
2131  bra.deviceStateVec.real, bra.deviceStateVec.imag,
2132  ket.deviceStateVec.real, ket.deviceStateVec.imag,
2133  numValuesToReduce,
2134  bra.firstLevelReduction);
2135  firstTime = 0;
2136  } else {
2137  cudaDeviceSynchronize();
2138  copySharedReduceBlock<<<numCUDABlocks, valuesPerCUDABlock/2, sharedMemSize>>>(
2139  bra.firstLevelReduction,
2140  bra.secondLevelReduction, valuesPerCUDABlock);
2141  cudaDeviceSynchronize();
2142  swapDouble(&(bra.firstLevelReduction), &(bra.secondLevelReduction));
2143  }
2144  numValuesToReduce = numValuesToReduce/maxReducedPerLevel;
2145  }
2146  cudaMemcpy(&innerProdImag, bra.firstLevelReduction, sizeof(qreal), cudaMemcpyDeviceToHost);
2147 
2148  // return complex
2149  Complex innerProd;
2150  innerProd.real = innerProdReal;
2151  innerProd.imag = innerProdImag;
2152  return innerProd;
2153 }

References copySharedReduceBlock(), Qureg::deviceStateVec, Qureg::firstLevelReduction, Complex::imag, Qureg::numAmpsPerChunk, qreal, Complex::real, REDUCE_SHARED_SIZE, Qureg::secondLevelReduction, and swapDouble().

Referenced by calcInnerProduct(), statevec_calcExpecPauliProd(), and statevec_calcFidelity().

◆ statevec_calcInnerProductKernel()

__global__ void statevec_calcInnerProductKernel ( int  getRealComp,
qreal *  vecReal1,
qreal *  vecImag1,
qreal *  vecReal2,
qreal *  vecImag2,
long long int  numTermsToSum,
qreal *  reducedArray 
)

computes either a real or imag term in the inner product

Definition at line 2038 of file QuEST_gpu.cu.

2042 {
2043  long long int index = blockIdx.x*blockDim.x + threadIdx.x;
2044  if (index >= numTermsToSum) return;
2045 
2046  // choose whether to calculate the real or imaginary term of the inner product
2047  qreal innerProdTerm;
2048  if (getRealComp)
2049  innerProdTerm = vecReal1[index]*vecReal2[index] + vecImag1[index]*vecImag2[index];
2050  else
2051  innerProdTerm = vecReal1[index]*vecImag2[index] - vecImag1[index]*vecReal2[index];
2052 
2053  // array of each thread's collected sum term, to be summed
2054  extern __shared__ qreal tempReductionArray[];
2055  tempReductionArray[threadIdx.x] = innerProdTerm;
2056  __syncthreads();
2057 
2058  // every second thread reduces
2059  if (threadIdx.x<blockDim.x/2)
2060  reduceBlock(tempReductionArray, reducedArray, blockDim.x);
2061 }

References qreal, and reduceBlock().

◆ statevec_calcProbOfOutcome()

qreal statevec_calcProbOfOutcome ( Qureg  qureg,
const int  measureQubit,
int  outcome 
)

Definition at line 1950 of file QuEST_gpu.cu.

1951 {
1952  qreal outcomeProb = statevec_findProbabilityOfZero(qureg, measureQubit);
1953  if (outcome==1)
1954  outcomeProb = 1.0 - outcomeProb;
1955  return outcomeProb;
1956 }

References qreal, and statevec_findProbabilityOfZero().

Referenced by calcProbOfOutcome(), collapseToOutcome(), and statevec_measureWithStats().

◆ statevec_calcTotalProb()

qreal statevec_calcTotalProb ( Qureg  qureg)

Definition at line 1499 of file QuEST_gpu.cu.

1499  {
1500  /* IJB - implemented using Kahan summation for greater accuracy at a slight floating
1501  point operation overhead. For more details see https://en.wikipedia.org/wiki/Kahan_summation_algorithm */
1502  /* Don't change the bracketing in this routine! */
1503  qreal pTotal=0;
1504  qreal y, t, c;
1505  long long int index;
1506  long long int numAmpsPerRank = qureg.numAmpsPerChunk;
1507 
1508  copyStateFromGPU(qureg);
1509 
1510  c = 0.0;
1511  for (index=0; index<numAmpsPerRank; index++){
1512  /* Perform pTotal+=qureg.stateVec.real[index]*qureg.stateVec.real[index]; by Kahan */
1513  // pTotal+=qureg.stateVec.real[index]*qureg.stateVec.real[index];
1514  y = qureg.stateVec.real[index]*qureg.stateVec.real[index] - c;
1515  t = pTotal + y;
1516  c = ( t - pTotal ) - y;
1517  pTotal = t;
1518 
1519  /* Perform pTotal+=qureg.stateVec.imag[index]*qureg.stateVec.imag[index]; by Kahan */
1520  //pTotal+=qureg.stateVec.imag[index]*qureg.stateVec.imag[index];
1521  y = qureg.stateVec.imag[index]*qureg.stateVec.imag[index] - c;
1522  t = pTotal + y;
1523  c = ( t - pTotal ) - y;
1524  pTotal = t;
1525 
1526 
1527  }
1528  return pTotal;
1529 }

References copyStateFromGPU(), Qureg::numAmpsPerChunk, qreal, and Qureg::stateVec.

Referenced by calcTotalProb().

◆ statevec_cloneQureg()

void statevec_cloneQureg ( Qureg  targetQureg,
Qureg  copyQureg 
)

works for both statevectors and density matrices

Definition at line 170 of file QuEST_gpu.cu.

170  {
171 
172  // copy copyQureg's GPU statevec to targetQureg's GPU statevec
173  cudaDeviceSynchronize();
174  cudaMemcpy(
175  targetQureg.deviceStateVec.real,
176  copyQureg.deviceStateVec.real,
177  targetQureg.numAmpsPerChunk*sizeof(*(targetQureg.deviceStateVec.real)),
178  cudaMemcpyDeviceToDevice);
179  cudaMemcpy(
180  targetQureg.deviceStateVec.imag,
181  copyQureg.deviceStateVec.imag,
182  targetQureg.numAmpsPerChunk*sizeof(*(targetQureg.deviceStateVec.imag)),
183  cudaMemcpyDeviceToDevice);
184 }

References Qureg::deviceStateVec, and Qureg::numAmpsPerChunk.

Referenced by cloneQureg(), createCloneQureg(), initPureState(), and statevec_calcExpecPauliProd().

◆ statevec_collapseToKnownProbOutcome()

void statevec_collapseToKnownProbOutcome ( Qureg  qureg,
const int  measureQubit,
int  outcome,
qreal  outcomeProb 
)

Definition at line 2445 of file QuEST_gpu.cu.

2446 {
2447  int threadsPerCUDABlock, CUDABlocks;
2448  threadsPerCUDABlock = 128;
2449  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk>>1)/threadsPerCUDABlock);
2450  statevec_collapseToKnownProbOutcomeKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, measureQubit, outcome, outcomeProb);
2451 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by collapseToOutcome(), and statevec_measureWithStats().

◆ statevec_collapseToKnownProbOutcomeKernel()

__global__ void statevec_collapseToKnownProbOutcomeKernel ( Qureg  qureg,
int  measureQubit,
int  outcome,
qreal  totalProbability 
)

Definition at line 2388 of file QuEST_gpu.cu.

2389 {
2390  // ----- sizes
2391  long long int sizeBlock, // size of blocks
2392  sizeHalfBlock; // size of blocks halved
2393  // ----- indices
2394  long long int thisBlock, // current block
2395  index; // current index for first half block
2396  // ----- measured probability
2397  qreal renorm; // probability (returned) value
2398  // ----- temp variables
2399  long long int thisTask; // task based approach for expose loop with small granularity
2400  // (good for shared memory parallelism)
2401  long long int numTasks=qureg.numAmpsPerChunk>>1;
2402 
2403  // ---------------------------------------------------------------- //
2404  // dimensions //
2405  // ---------------------------------------------------------------- //
2406  sizeHalfBlock = 1LL << (measureQubit); // number of state vector elements to sum,
2407  // and then the number to skip
2408  sizeBlock = 2LL * sizeHalfBlock; // size of blocks (pairs of measure and skip entries)
2409 
2410  // ---------------------------------------------------------------- //
2411  // find probability //
2412  // ---------------------------------------------------------------- //
2413 
2414  //
2415  // --- task-based shared-memory parallel implementation
2416  //
2417  renorm=1/sqrt(totalProbability);
2418  qreal *stateVecReal = qureg.deviceStateVec.real;
2419  qreal *stateVecImag = qureg.deviceStateVec.imag;
2420 
2421  thisTask = blockIdx.x*blockDim.x + threadIdx.x;
2422  if (thisTask>=numTasks) return;
2423  thisBlock = thisTask / sizeHalfBlock;
2424  index = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2425 
2426  if (outcome==0){
2427  stateVecReal[index]=stateVecReal[index]*renorm;
2428  stateVecImag[index]=stateVecImag[index]*renorm;
2429 
2430  stateVecReal[index+sizeHalfBlock]=0;
2431  stateVecImag[index+sizeHalfBlock]=0;
2432  } else if (outcome==1){
2433  stateVecReal[index]=0;
2434  stateVecImag[index]=0;
2435 
2436  stateVecReal[index+sizeHalfBlock]=stateVecReal[index+sizeHalfBlock]*renorm;
2437  stateVecImag[index+sizeHalfBlock]=stateVecImag[index+sizeHalfBlock]*renorm;
2438  }
2439 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

◆ statevec_compactUnitary()

void statevec_compactUnitary ( Qureg  qureg,
const int  targetQubit,
Complex  alpha,
Complex  beta 
)

Definition at line 722 of file QuEST_gpu.cu.

723 {
724  int threadsPerCUDABlock, CUDABlocks;
725  threadsPerCUDABlock = 128;
726  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk>>1)/threadsPerCUDABlock);
727  statevec_compactUnitaryKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, targetQubit, alpha, beta);
728 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by compactUnitary(), statevec_multiRotatePauli(), statevec_rotateAroundAxis(), and statevec_rotateAroundAxisConj().

◆ statevec_compactUnitaryKernel()

__global__ void statevec_compactUnitaryKernel ( Qureg  qureg,
const int  rotQubit,
Complex  alpha,
Complex  beta 
)

fix – no necessary for GPU version

Definition at line 667 of file QuEST_gpu.cu.

667  {
668  // ----- sizes
669  long long int sizeBlock, // size of blocks
670  sizeHalfBlock; // size of blocks halved
671  // ----- indices
672  long long int thisBlock, // current block
673  indexUp,indexLo; // current index and corresponding index in lower half block
674 
675  // ----- temp variables
676  qreal stateRealUp,stateRealLo, // storage for previous state values
677  stateImagUp,stateImagLo; // (used in updates)
678  // ----- temp variables
679  long long int thisTask; // task based approach for expose loop with small granularity
680  const long long int numTasks=qureg.numAmpsPerChunk>>1;
681 
682  sizeHalfBlock = 1LL << rotQubit; // size of blocks halved
683  sizeBlock = 2LL * sizeHalfBlock; // size of blocks
684 
685  // ---------------------------------------------------------------- //
686  // rotate //
687  // ---------------------------------------------------------------- //
688 
690  qreal *stateVecReal = qureg.deviceStateVec.real;
691  qreal *stateVecImag = qureg.deviceStateVec.imag;
692  qreal alphaImag=alpha.imag, alphaReal=alpha.real;
693  qreal betaImag=beta.imag, betaReal=beta.real;
694 
695  thisTask = blockIdx.x*blockDim.x + threadIdx.x;
696  if (thisTask>=numTasks) return;
697 
698  thisBlock = thisTask / sizeHalfBlock;
699  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
700  indexLo = indexUp + sizeHalfBlock;
701 
702  // store current state vector values in temp variables
703  stateRealUp = stateVecReal[indexUp];
704  stateImagUp = stateVecImag[indexUp];
705 
706  stateRealLo = stateVecReal[indexLo];
707  stateImagLo = stateVecImag[indexLo];
708 
709  // state[indexUp] = alpha * state[indexUp] - conj(beta) * state[indexLo]
710  stateVecReal[indexUp] = alphaReal*stateRealUp - alphaImag*stateImagUp
711  - betaReal*stateRealLo - betaImag*stateImagLo;
712  stateVecImag[indexUp] = alphaReal*stateImagUp + alphaImag*stateRealUp
713  - betaReal*stateImagLo + betaImag*stateRealLo;
714 
715  // state[indexLo] = beta * state[indexUp] + conj(alpha) * state[indexLo]
716  stateVecReal[indexLo] = betaReal*stateRealUp - betaImag*stateImagUp
717  + alphaReal*stateRealLo + alphaImag*stateImagLo;
718  stateVecImag[indexLo] = betaReal*stateImagUp + betaImag*stateRealUp
719  + alphaReal*stateImagLo - alphaImag*stateRealLo;
720 }

References Qureg::deviceStateVec, Complex::imag, Qureg::numAmpsPerChunk, qreal, and Complex::real.

◆ statevec_compareStates()

int statevec_compareStates ( Qureg  mq1,
Qureg  mq2,
qreal  precision 
)

Definition at line 649 of file QuEST_gpu.cu.

649  {
650  qreal diff;
651  int chunkSize = mq1.numAmpsPerChunk;
652 
653  copyStateFromGPU(mq1);
654  copyStateFromGPU(mq2);
655 
656  for (int i=0; i<chunkSize; i++){
657  diff = mq1.stateVec.real[i] - mq2.stateVec.real[i];
658  if (diff<0) diff *= -1;
659  if (diff>precision) return 0;
660  diff = mq1.stateVec.imag[i] - mq2.stateVec.imag[i];
661  if (diff<0) diff *= -1;
662  if (diff>precision) return 0;
663  }
664  return 1;
665 }

References copyStateFromGPU(), Qureg::numAmpsPerChunk, qreal, and Qureg::stateVec.

Referenced by compareStates().

◆ statevec_controlledCompactUnitary()

void statevec_controlledCompactUnitary ( Qureg  qureg,
const int  controlQubit,
const int  targetQubit,
Complex  alpha,
Complex  beta 
)

Definition at line 789 of file QuEST_gpu.cu.

790 {
791  int threadsPerCUDABlock, CUDABlocks;
792  threadsPerCUDABlock = 128;
793  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk>>1)/threadsPerCUDABlock);
794  statevec_controlledCompactUnitaryKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, controlQubit, targetQubit, alpha, beta);
795 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by controlledCompactUnitary(), statevec_controlledRotateAroundAxis(), and statevec_controlledRotateAroundAxisConj().

◆ statevec_controlledCompactUnitaryKernel()

__global__ void statevec_controlledCompactUnitaryKernel ( Qureg  qureg,
const int  controlQubit,
const int  targetQubit,
Complex  alpha,
Complex  beta 
)

fix – no necessary for GPU version

Definition at line 730 of file QuEST_gpu.cu.

730  {
731  // ----- sizes
732  long long int sizeBlock, // size of blocks
733  sizeHalfBlock; // size of blocks halved
734  // ----- indices
735  long long int thisBlock, // current block
736  indexUp,indexLo; // current index and corresponding index in lower half block
737 
738  // ----- temp variables
739  qreal stateRealUp,stateRealLo, // storage for previous state values
740  stateImagUp,stateImagLo; // (used in updates)
741  // ----- temp variables
742  long long int thisTask; // task based approach for expose loop with small granularity
743  const long long int numTasks=qureg.numAmpsPerChunk>>1;
744  int controlBit;
745 
746  sizeHalfBlock = 1LL << targetQubit; // size of blocks halved
747  sizeBlock = 2LL * sizeHalfBlock; // size of blocks
748 
749  // ---------------------------------------------------------------- //
750  // rotate //
751  // ---------------------------------------------------------------- //
752 
754  qreal *stateVecReal = qureg.deviceStateVec.real;
755  qreal *stateVecImag = qureg.deviceStateVec.imag;
756  qreal alphaImag=alpha.imag, alphaReal=alpha.real;
757  qreal betaImag=beta.imag, betaReal=beta.real;
758 
759  thisTask = blockIdx.x*blockDim.x + threadIdx.x;
760  if (thisTask>=numTasks) return;
761 
762  thisBlock = thisTask / sizeHalfBlock;
763  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
764  indexLo = indexUp + sizeHalfBlock;
765 
766  controlBit = extractBit(controlQubit, indexUp);
767  if (controlBit){
768  // store current state vector values in temp variables
769  stateRealUp = stateVecReal[indexUp];
770  stateImagUp = stateVecImag[indexUp];
771 
772  stateRealLo = stateVecReal[indexLo];
773  stateImagLo = stateVecImag[indexLo];
774 
775  // state[indexUp] = alpha * state[indexUp] - conj(beta) * state[indexLo]
776  stateVecReal[indexUp] = alphaReal*stateRealUp - alphaImag*stateImagUp
777  - betaReal*stateRealLo - betaImag*stateImagLo;
778  stateVecImag[indexUp] = alphaReal*stateImagUp + alphaImag*stateRealUp
779  - betaReal*stateImagLo + betaImag*stateRealLo;
780 
781  // state[indexLo] = beta * state[indexUp] + conj(alpha) * state[indexLo]
782  stateVecReal[indexLo] = betaReal*stateRealUp - betaImag*stateImagUp
783  + alphaReal*stateRealLo + alphaImag*stateImagLo;
784  stateVecImag[indexLo] = betaReal*stateImagUp + betaImag*stateRealUp
785  + alphaReal*stateImagLo - alphaImag*stateRealLo;
786  }
787 }

References Qureg::deviceStateVec, extractBit(), Complex::imag, Qureg::numAmpsPerChunk, qreal, and Complex::real.

◆ statevec_controlledNot()

void statevec_controlledNot ( Qureg  qureg,
const int  controlQubit,
const int  targetQubit 
)

Definition at line 1717 of file QuEST_gpu.cu.

1718 {
1719  int threadsPerCUDABlock, CUDABlocks;
1720  threadsPerCUDABlock = 128;
1721  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk)/threadsPerCUDABlock);
1722  statevec_controlledNotKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, controlQubit, targetQubit);
1723 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by controlledNot().

◆ statevec_controlledNotKernel()

__global__ void statevec_controlledNotKernel ( Qureg  qureg,
const int  controlQubit,
const int  targetQubit 
)

Definition at line 1678 of file QuEST_gpu.cu.

1679 {
1680  long long int index;
1681  long long int sizeBlock, // size of blocks
1682  sizeHalfBlock; // size of blocks halved
1683  long long int stateVecSize;
1684  int controlBit;
1685 
1686  // ----- temp variables
1687  qreal stateRealUp, // storage for previous state values
1688  stateImagUp; // (used in updates)
1689  long long int thisBlock, // current block
1690  indexUp,indexLo; // current index and corresponding index in lower half block
1691  sizeHalfBlock = 1LL << targetQubit; // size of blocks halved
1692  sizeBlock = 2LL * sizeHalfBlock; // size of blocks
1693 
1694  stateVecSize = qureg.numAmpsPerChunk;
1695  qreal *stateVecReal = qureg.deviceStateVec.real;
1696  qreal *stateVecImag = qureg.deviceStateVec.imag;
1697 
1698  index = blockIdx.x*blockDim.x + threadIdx.x;
1699  if (index>=(stateVecSize>>1)) return;
1700  thisBlock = index / sizeHalfBlock;
1701  indexUp = thisBlock*sizeBlock + index%sizeHalfBlock;
1702  indexLo = indexUp + sizeHalfBlock;
1703 
1704  controlBit = extractBit(controlQubit, indexUp);
1705  if (controlBit){
1706  stateRealUp = stateVecReal[indexUp];
1707  stateImagUp = stateVecImag[indexUp];
1708 
1709  stateVecReal[indexUp] = stateVecReal[indexLo];
1710  stateVecImag[indexUp] = stateVecImag[indexLo];
1711 
1712  stateVecReal[indexLo] = stateRealUp;
1713  stateVecImag[indexLo] = stateImagUp;
1714  }
1715 }

References Qureg::deviceStateVec, extractBit(), Qureg::numAmpsPerChunk, and qreal.

◆ statevec_controlledPauliY()

void statevec_controlledPauliY ( Qureg  qureg,
const int  controlQubit,
const int  targetQubit 
)

Definition at line 1323 of file QuEST_gpu.cu.

1324 {
1325  int conjFactor = 1;
1326  int threadsPerCUDABlock, CUDABlocks;
1327  threadsPerCUDABlock = 128;
1328  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk)/threadsPerCUDABlock);
1329  statevec_controlledPauliYKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, controlQubit, targetQubit, conjFactor);
1330 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by controlledPauliY().

◆ statevec_controlledPauliYConj()

void statevec_controlledPauliYConj ( Qureg  qureg,
const int  controlQubit,
const int  targetQubit 
)

Definition at line 1332 of file QuEST_gpu.cu.

1333 {
1334  int conjFactor = -1;
1335  int threadsPerCUDABlock, CUDABlocks;
1336  threadsPerCUDABlock = 128;
1337  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk)/threadsPerCUDABlock);
1338  statevec_controlledPauliYKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, controlQubit, targetQubit, conjFactor);
1339 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by controlledPauliY().

◆ statevec_controlledPauliYKernel()

__global__ void statevec_controlledPauliYKernel ( Qureg  qureg,
const int  controlQubit,
const int  targetQubit,
const int  conjFac 
)

Definition at line 1287 of file QuEST_gpu.cu.

1288 {
1289  long long int index;
1290  long long int sizeBlock, sizeHalfBlock;
1291  long long int stateVecSize;
1292  int controlBit;
1293 
1294  qreal stateRealUp, stateImagUp;
1295  long long int thisBlock, indexUp, indexLo;
1296  sizeHalfBlock = 1LL << targetQubit;
1297  sizeBlock = 2LL * sizeHalfBlock;
1298 
1299  stateVecSize = qureg.numAmpsPerChunk;
1300  qreal *stateVecReal = qureg.deviceStateVec.real;
1301  qreal *stateVecImag = qureg.deviceStateVec.imag;
1302 
1303  index = blockIdx.x*blockDim.x + threadIdx.x;
1304  if (index>=(stateVecSize>>1)) return;
1305  thisBlock = index / sizeHalfBlock;
1306  indexUp = thisBlock*sizeBlock + index%sizeHalfBlock;
1307  indexLo = indexUp + sizeHalfBlock;
1308 
1309  controlBit = extractBit(controlQubit, indexUp);
1310  if (controlBit){
1311 
1312  stateRealUp = stateVecReal[indexUp];
1313  stateImagUp = stateVecImag[indexUp];
1314 
1315  // update under +-{{0, -i}, {i, 0}}
1316  stateVecReal[indexUp] = conjFac * stateVecImag[indexLo];
1317  stateVecImag[indexUp] = conjFac * -stateVecReal[indexLo];
1318  stateVecReal[indexLo] = conjFac * -stateImagUp;
1319  stateVecImag[indexLo] = conjFac * stateRealUp;
1320  }
1321 }

References Qureg::deviceStateVec, extractBit(), Qureg::numAmpsPerChunk, and qreal.

◆ statevec_controlledPhaseFlip()

void statevec_controlledPhaseFlip ( Qureg  qureg,
const int  idQubit1,
const int  idQubit2 
)

Definition at line 1552 of file QuEST_gpu.cu.

1553 {
1554  int threadsPerCUDABlock, CUDABlocks;
1555  threadsPerCUDABlock = 128;
1556  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk)/threadsPerCUDABlock);
1557  statevec_controlledPhaseFlipKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, idQubit1, idQubit2);
1558 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by controlledPhaseFlip().

◆ statevec_controlledPhaseFlipKernel()

__global__ void statevec_controlledPhaseFlipKernel ( Qureg  qureg,
const int  idQubit1,
const int  idQubit2 
)

Definition at line 1531 of file QuEST_gpu.cu.

1532 {
1533  long long int index;
1534  long long int stateVecSize;
1535  int bit1, bit2;
1536 
1537  stateVecSize = qureg.numAmpsPerChunk;
1538  qreal *stateVecReal = qureg.deviceStateVec.real;
1539  qreal *stateVecImag = qureg.deviceStateVec.imag;
1540 
1541  index = blockIdx.x*blockDim.x + threadIdx.x;
1542  if (index>=stateVecSize) return;
1543 
1544  bit1 = extractBit (idQubit1, index);
1545  bit2 = extractBit (idQubit2, index);
1546  if (bit1 && bit2) {
1547  stateVecReal [index] = - stateVecReal [index];
1548  stateVecImag [index] = - stateVecImag [index];
1549  }
1550 }

References Qureg::deviceStateVec, extractBit(), Qureg::numAmpsPerChunk, and qreal.

◆ statevec_controlledPhaseShift()

void statevec_controlledPhaseShift ( Qureg  qureg,
const int  idQubit1,
const int  idQubit2,
qreal  angle 
)

Definition at line 1405 of file QuEST_gpu.cu.

1406 {
1407  qreal cosAngle = cos(angle);
1408  qreal sinAngle = sin(angle);
1409 
1410  int threadsPerCUDABlock, CUDABlocks;
1411  threadsPerCUDABlock = 128;
1412  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk)/threadsPerCUDABlock);
1413  statevec_controlledPhaseShiftKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, idQubit1, idQubit2, cosAngle, sinAngle);
1414 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by controlledPhaseShift().

◆ statevec_controlledPhaseShiftKernel()

__global__ void statevec_controlledPhaseShiftKernel ( Qureg  qureg,
const int  idQubit1,
const int  idQubit2,
qreal  cosAngle,
qreal  sinAngle 
)

Definition at line 1380 of file QuEST_gpu.cu.

1381 {
1382  long long int index;
1383  long long int stateVecSize;
1384  int bit1, bit2;
1385  qreal stateRealLo, stateImagLo;
1386 
1387  stateVecSize = qureg.numAmpsPerChunk;
1388  qreal *stateVecReal = qureg.deviceStateVec.real;
1389  qreal *stateVecImag = qureg.deviceStateVec.imag;
1390 
1391  index = blockIdx.x*blockDim.x + threadIdx.x;
1392  if (index>=stateVecSize) return;
1393 
1394  bit1 = extractBit (idQubit1, index);
1395  bit2 = extractBit (idQubit2, index);
1396  if (bit1 && bit2) {
1397  stateRealLo = stateVecReal[index];
1398  stateImagLo = stateVecImag[index];
1399 
1400  stateVecReal[index] = cosAngle*stateRealLo - sinAngle*stateImagLo;
1401  stateVecImag[index] = sinAngle*stateRealLo + cosAngle*stateImagLo;
1402  }
1403 }

References Qureg::deviceStateVec, extractBit(), Qureg::numAmpsPerChunk, and qreal.

◆ statevec_controlledUnitary()

void statevec_controlledUnitary ( Qureg  qureg,
const int  controlQubit,
const int  targetQubit,
ComplexMatrix2  u 
)

Definition at line 1115 of file QuEST_gpu.cu.

1116 {
1117  int threadsPerCUDABlock, CUDABlocks;
1118  threadsPerCUDABlock = 128;
1119  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk>>1)/threadsPerCUDABlock);
1120  statevec_controlledUnitaryKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, controlQubit, targetQubit, argifyMatrix2(u));
1121 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by controlledUnitary().

◆ statevec_controlledUnitaryKernel()

__global__ void statevec_controlledUnitaryKernel ( Qureg  qureg,
const int  controlQubit,
const int  targetQubit,
ArgMatrix2  u 
)

fix – no necessary for GPU version

Definition at line 1057 of file QuEST_gpu.cu.

1057  {
1058  // ----- sizes
1059  long long int sizeBlock, // size of blocks
1060  sizeHalfBlock; // size of blocks halved
1061  // ----- indices
1062  long long int thisBlock, // current block
1063  indexUp,indexLo; // current index and corresponding index in lower half block
1064 
1065  // ----- temp variables
1066  qreal stateRealUp,stateRealLo, // storage for previous state values
1067  stateImagUp,stateImagLo; // (used in updates)
1068  // ----- temp variables
1069  long long int thisTask; // task based approach for expose loop with small granularity
1070  const long long int numTasks=qureg.numAmpsPerChunk>>1;
1071 
1072  int controlBit;
1073 
1074  sizeHalfBlock = 1LL << targetQubit; // size of blocks halved
1075  sizeBlock = 2LL * sizeHalfBlock; // size of blocks
1076 
1077  // ---------------------------------------------------------------- //
1078  // rotate //
1079  // ---------------------------------------------------------------- //
1080 
1082  qreal *stateVecReal = qureg.deviceStateVec.real;
1083  qreal *stateVecImag = qureg.deviceStateVec.imag;
1084 
1085  thisTask = blockIdx.x*blockDim.x + threadIdx.x;
1086  if (thisTask>=numTasks) return;
1087 
1088  thisBlock = thisTask / sizeHalfBlock;
1089  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
1090  indexLo = indexUp + sizeHalfBlock;
1091 
1092  // store current state vector values in temp variables
1093  stateRealUp = stateVecReal[indexUp];
1094  stateImagUp = stateVecImag[indexUp];
1095 
1096  stateRealLo = stateVecReal[indexLo];
1097  stateImagLo = stateVecImag[indexLo];
1098 
1099  controlBit = extractBit(controlQubit, indexUp);
1100  if (controlBit){
1101  // state[indexUp] = u00 * state[indexUp] + u01 * state[indexLo]
1102  stateVecReal[indexUp] = u.r0c0.real*stateRealUp - u.r0c0.imag*stateImagUp
1103  + u.r0c1.real*stateRealLo - u.r0c1.imag*stateImagLo;
1104  stateVecImag[indexUp] = u.r0c0.real*stateImagUp + u.r0c0.imag*stateRealUp
1105  + u.r0c1.real*stateImagLo + u.r0c1.imag*stateRealLo;
1106 
1107  // state[indexLo] = u10 * state[indexUp] + u11 * state[indexLo]
1108  stateVecReal[indexLo] = u.r1c0.real*stateRealUp - u.r1c0.imag*stateImagUp
1109  + u.r1c1.real*stateRealLo - u.r1c1.imag*stateImagLo;
1110  stateVecImag[indexLo] = u.r1c0.real*stateImagUp + u.r1c0.imag*stateRealUp
1111  + u.r1c1.real*stateImagLo + u.r1c1.imag*stateRealLo;
1112  }
1113 }

References Qureg::deviceStateVec, extractBit(), Qureg::numAmpsPerChunk, and qreal.

◆ statevec_createQureg()

void statevec_createQureg ( Qureg *  qureg,
int  numQubits,
QuESTEnv  env 
)

Definition at line 275 of file QuEST_gpu.cu.

276 {
277  // allocate CPU memory
278  long long int numAmps = 1L << numQubits;
279  long long int numAmpsPerRank = numAmps/env.numRanks;
280  qureg->stateVec.real = (qreal*) malloc(numAmpsPerRank * sizeof(qureg->stateVec.real));
281  qureg->stateVec.imag = (qreal*) malloc(numAmpsPerRank * sizeof(qureg->stateVec.imag));
282  if (env.numRanks>1){
283  qureg->pairStateVec.real = (qreal*) malloc(numAmpsPerRank * sizeof(qureg->pairStateVec.real));
284  qureg->pairStateVec.imag = (qreal*) malloc(numAmpsPerRank * sizeof(qureg->pairStateVec.imag));
285  }
286 
287  // check cpu memory allocation was successful
288  if ( (!(qureg->stateVec.real) || !(qureg->stateVec.imag))
289  && numAmpsPerRank ) {
290  printf("Could not allocate memory!\n");
291  exit (EXIT_FAILURE);
292  }
293  if ( env.numRanks>1 && (!(qureg->pairStateVec.real) || !(qureg->pairStateVec.imag))
294  && numAmpsPerRank ) {
295  printf("Could not allocate memory!\n");
296  exit (EXIT_FAILURE);
297  }
298 
299  qureg->numQubitsInStateVec = numQubits;
300  qureg->numAmpsPerChunk = numAmpsPerRank;
301  qureg->numAmpsTotal = numAmps;
302  qureg->chunkId = env.rank;
303  qureg->numChunks = env.numRanks;
304  qureg->isDensityMatrix = 0;
305 
306  // allocate GPU memory
307  cudaMalloc(&(qureg->deviceStateVec.real), qureg->numAmpsPerChunk*sizeof(*(qureg->deviceStateVec.real)));
308  cudaMalloc(&(qureg->deviceStateVec.imag), qureg->numAmpsPerChunk*sizeof(*(qureg->deviceStateVec.imag)));
309  cudaMalloc(&(qureg->firstLevelReduction), ceil(qureg->numAmpsPerChunk/(qreal)REDUCE_SHARED_SIZE)*sizeof(qreal));
310  cudaMalloc(&(qureg->secondLevelReduction), ceil(qureg->numAmpsPerChunk/(qreal)(REDUCE_SHARED_SIZE*REDUCE_SHARED_SIZE))*
311  sizeof(qreal));
312 
313  // check gpu memory allocation was successful
314  if (!(qureg->deviceStateVec.real) || !(qureg->deviceStateVec.imag)){
315  printf("Could not allocate memory on GPU!\n");
316  exit (EXIT_FAILURE);
317  }
318 
319 }

References Qureg::chunkId, Qureg::deviceStateVec, Qureg::firstLevelReduction, Qureg::isDensityMatrix, Qureg::numAmpsPerChunk, Qureg::numAmpsTotal, Qureg::numChunks, Qureg::numQubitsInStateVec, QuESTEnv::numRanks, Qureg::pairStateVec, qreal, QuESTEnv::rank, REDUCE_SHARED_SIZE, Qureg::secondLevelReduction, and Qureg::stateVec.

Referenced by createCloneQureg(), createDensityQureg(), and createQureg().

◆ statevec_destroyQureg()

void statevec_destroyQureg ( Qureg  qureg,
QuESTEnv  env 
)

Definition at line 321 of file QuEST_gpu.cu.

322 {
323  // Free CPU memory
324  free(qureg.stateVec.real);
325  free(qureg.stateVec.imag);
326  if (env.numRanks>1){
327  free(qureg.pairStateVec.real);
328  free(qureg.pairStateVec.imag);
329  }
330 
331  // Free GPU memory
332  cudaFree(qureg.deviceStateVec.real);
333  cudaFree(qureg.deviceStateVec.imag);
334 }

References Qureg::deviceStateVec, QuESTEnv::numRanks, Qureg::pairStateVec, and Qureg::stateVec.

Referenced by destroyQureg().

◆ statevec_findProbabilityOfZero()

qreal statevec_findProbabilityOfZero ( Qureg  qureg,
const int  measureQubit 
)

Definition at line 1912 of file QuEST_gpu.cu.

1913 {
1914  long long int numValuesToReduce = qureg.numAmpsPerChunk>>1;
1915  int valuesPerCUDABlock, numCUDABlocks, sharedMemSize;
1916  qreal stateProb=0;
1917  int firstTime=1;
1918  int maxReducedPerLevel = REDUCE_SHARED_SIZE;
1919 
1920  while(numValuesToReduce>1){
1921  if (numValuesToReduce<maxReducedPerLevel){
1922  // Need less than one CUDA block to reduce values
1923  valuesPerCUDABlock = numValuesToReduce;
1924  numCUDABlocks = 1;
1925  } else {
1926  // Use full CUDA blocks, with block size constrained by shared mem usage
1927  valuesPerCUDABlock = maxReducedPerLevel;
1928  numCUDABlocks = ceil((qreal)numValuesToReduce/valuesPerCUDABlock);
1929  }
1930  sharedMemSize = valuesPerCUDABlock*sizeof(qreal);
1931 
1932  if (firstTime){
1933  statevec_findProbabilityOfZeroKernel<<<numCUDABlocks, valuesPerCUDABlock, sharedMemSize>>>(
1934  qureg, measureQubit, qureg.firstLevelReduction);
1935  firstTime=0;
1936  } else {
1937  cudaDeviceSynchronize();
1938  copySharedReduceBlock<<<numCUDABlocks, valuesPerCUDABlock/2, sharedMemSize>>>(
1939  qureg.firstLevelReduction,
1940  qureg.secondLevelReduction, valuesPerCUDABlock);
1941  cudaDeviceSynchronize();
1942  swapDouble(&(qureg.firstLevelReduction), &(qureg.secondLevelReduction));
1943  }
1944  numValuesToReduce = numValuesToReduce/maxReducedPerLevel;
1945  }
1946  cudaMemcpy(&stateProb, qureg.firstLevelReduction, sizeof(qreal), cudaMemcpyDeviceToHost);
1947  return stateProb;
1948 }

References copySharedReduceBlock(), Qureg::firstLevelReduction, Qureg::numAmpsPerChunk, qreal, REDUCE_SHARED_SIZE, Qureg::secondLevelReduction, and swapDouble().

Referenced by statevec_calcProbOfOutcome().

◆ statevec_findProbabilityOfZeroKernel()

__global__ void statevec_findProbabilityOfZeroKernel ( Qureg  qureg,
const int  measureQubit,
qreal *  reducedArray 
)

Definition at line 1798 of file QuEST_gpu.cu.

1800  {
1801  // ----- sizes
1802  long long int sizeBlock, // size of blocks
1803  sizeHalfBlock; // size of blocks halved
1804  // ----- indices
1805  long long int thisBlock, // current block
1806  index; // current index for first half block
1807  // ----- temp variables
1808  long long int thisTask; // task based approach for expose loop with small granularity
1809  long long int numTasks=qureg.numAmpsPerChunk>>1;
1810  // (good for shared memory parallelism)
1811 
1812  extern __shared__ qreal tempReductionArray[];
1813 
1814  // ---------------------------------------------------------------- //
1815  // dimensions //
1816  // ---------------------------------------------------------------- //
1817  sizeHalfBlock = 1LL << (measureQubit); // number of state vector elements to sum,
1818  // and then the number to skip
1819  sizeBlock = 2LL * sizeHalfBlock; // size of blocks (pairs of measure and skip entries)
1820 
1821  // ---------------------------------------------------------------- //
1822  // find probability //
1823  // ---------------------------------------------------------------- //
1824 
1825  //
1826  // --- task-based shared-memory parallel implementation
1827  //
1828 
1829  qreal *stateVecReal = qureg.deviceStateVec.real;
1830  qreal *stateVecImag = qureg.deviceStateVec.imag;
1831 
1832  thisTask = blockIdx.x*blockDim.x + threadIdx.x;
1833  if (thisTask>=numTasks) return;
1834 
1835  thisBlock = thisTask / sizeHalfBlock;
1836  index = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
1837  qreal realVal, imagVal;
1838  realVal = stateVecReal[index];
1839  imagVal = stateVecImag[index];
1840  tempReductionArray[threadIdx.x] = realVal*realVal + imagVal*imagVal;
1841  __syncthreads();
1842 
1843  if (threadIdx.x<blockDim.x/2){
1844  reduceBlock(tempReductionArray, reducedArray, blockDim.x);
1845  }
1846 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, qreal, and reduceBlock().

◆ statevec_getImagAmp()

qreal statevec_getImagAmp ( Qureg  qureg,
long long int  index 
)

Definition at line 454 of file QuEST_gpu.cu.

454  {
455  qreal el=0;
456  cudaMemcpy(&el, &(qureg.deviceStateVec.imag[index]),
457  sizeof(*(qureg.deviceStateVec.imag)), cudaMemcpyDeviceToHost);
458  return el;
459 }

References Qureg::deviceStateVec, and qreal.

Referenced by getAmp(), getDensityAmp(), getImagAmp(), and statevec_getProbAmp().

◆ statevec_getRealAmp()

qreal statevec_getRealAmp ( Qureg  qureg,
long long int  index 
)

Definition at line 447 of file QuEST_gpu.cu.

447  {
448  qreal el=0;
449  cudaMemcpy(&el, &(qureg.deviceStateVec.real[index]),
450  sizeof(*(qureg.deviceStateVec.real)), cudaMemcpyDeviceToHost);
451  return el;
452 }

References Qureg::deviceStateVec, and qreal.

Referenced by getAmp(), getDensityAmp(), getRealAmp(), and statevec_getProbAmp().

◆ statevec_hadamard()

void statevec_hadamard ( Qureg  qureg,
const int  targetQubit 
)

Definition at line 1670 of file QuEST_gpu.cu.

1671 {
1672  int threadsPerCUDABlock, CUDABlocks;
1673  threadsPerCUDABlock = 128;
1674  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk>>1)/threadsPerCUDABlock);
1675  statevec_hadamardKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, targetQubit);
1676 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by hadamard().

◆ statevec_hadamardKernel()

__global__ void statevec_hadamardKernel ( Qureg  qureg,
const int  targetQubit 
)

fix – no necessary for GPU version

Definition at line 1621 of file QuEST_gpu.cu.

1621  {
1622  // ----- sizes
1623  long long int sizeBlock, // size of blocks
1624  sizeHalfBlock; // size of blocks halved
1625  // ----- indices
1626  long long int thisBlock, // current block
1627  indexUp,indexLo; // current index and corresponding index in lower half block
1628 
1629  // ----- temp variables
1630  qreal stateRealUp,stateRealLo, // storage for previous state values
1631  stateImagUp,stateImagLo; // (used in updates)
1632  // ----- temp variables
1633  long long int thisTask; // task based approach for expose loop with small granularity
1634  const long long int numTasks=qureg.numAmpsPerChunk>>1;
1635 
1636  sizeHalfBlock = 1LL << targetQubit; // size of blocks halved
1637  sizeBlock = 2LL * sizeHalfBlock; // size of blocks
1638 
1639  // ---------------------------------------------------------------- //
1640  // rotate //
1641  // ---------------------------------------------------------------- //
1642 
1644  qreal *stateVecReal = qureg.deviceStateVec.real;
1645  qreal *stateVecImag = qureg.deviceStateVec.imag;
1646 
1647  qreal recRoot2 = 1.0/sqrt(2.0);
1648 
1649  thisTask = blockIdx.x*blockDim.x + threadIdx.x;
1650  if (thisTask>=numTasks) return;
1651 
1652  thisBlock = thisTask / sizeHalfBlock;
1653  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
1654  indexLo = indexUp + sizeHalfBlock;
1655 
1656  // store current state vector values in temp variables
1657  stateRealUp = stateVecReal[indexUp];
1658  stateImagUp = stateVecImag[indexUp];
1659 
1660  stateRealLo = stateVecReal[indexLo];
1661  stateImagLo = stateVecImag[indexLo];
1662 
1663  stateVecReal[indexUp] = recRoot2*(stateRealUp + stateRealLo);
1664  stateVecImag[indexUp] = recRoot2*(stateImagUp + stateImagLo);
1665 
1666  stateVecReal[indexLo] = recRoot2*(stateRealUp - stateRealLo);
1667  stateVecImag[indexLo] = recRoot2*(stateImagUp - stateImagLo);
1668 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

◆ statevec_initBlankState()

void statevec_initBlankState ( Qureg  qureg)

Definition at line 471 of file QuEST_gpu.cu.

472 {
473  int threadsPerCUDABlock, CUDABlocks;
474  threadsPerCUDABlock = 128;
475  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk)/threadsPerCUDABlock);
476  statevec_initBlankStateKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
477  qureg.numAmpsPerChunk,
478  qureg.deviceStateVec.real,
479  qureg.deviceStateVec.imag);
480 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

Referenced by initBlankState(), and statevec_applyPauliSum().

◆ statevec_initBlankStateKernel()

__global__ void statevec_initBlankStateKernel ( long long int  stateVecSize,
qreal *  stateVecReal,
qreal *  stateVecImag 
)

Definition at line 461 of file QuEST_gpu.cu.

461  {
462  long long int index;
463 
464  // initialise the statevector to be all-zeros
465  index = blockIdx.x*blockDim.x + threadIdx.x;
466  if (index>=stateVecSize) return;
467  stateVecReal[index] = 0.0;
468  stateVecImag[index] = 0.0;
469 }

◆ statevec_initClassicalState()

void statevec_initClassicalState ( Qureg  qureg,
long long int  stateInd 
)

Definition at line 546 of file QuEST_gpu.cu.

547 {
548  int threadsPerCUDABlock, CUDABlocks;
549  threadsPerCUDABlock = 128;
550  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk)/threadsPerCUDABlock);
551  statevec_initClassicalStateKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
552  qureg.numAmpsPerChunk,
553  qureg.deviceStateVec.real,
554  qureg.deviceStateVec.imag, stateInd);
555 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

Referenced by initClassicalState().

◆ statevec_initClassicalStateKernel()

__global__ void statevec_initClassicalStateKernel ( long long int  stateVecSize,
qreal *  stateVecReal,
qreal *  stateVecImag,
long long int  stateInd 
)

Definition at line 531 of file QuEST_gpu.cu.

531  {
532  long long int index;
533 
534  // initialise the state to |stateInd>
535  index = blockIdx.x*blockDim.x + threadIdx.x;
536  if (index>=stateVecSize) return;
537  stateVecReal[index] = 0.0;
538  stateVecImag[index] = 0.0;
539 
540  if (index==stateInd){
541  // classical state has probability 1
542  stateVecReal[stateInd] = 1.0;
543  stateVecImag[stateInd] = 0.0;
544  }
545 }

◆ statevec_initDebugState()

void statevec_initDebugState ( Qureg  qureg)

Initialise the state vector of probability amplitudes to an (unphysical) state with each component of each probability amplitude a unique floating point value.

For debugging processes

Parameters
[in,out]quregobject representing the set of qubits to be initialised

Definition at line 567 of file QuEST_gpu.cu.

568 {
569  int threadsPerCUDABlock, CUDABlocks;
570  threadsPerCUDABlock = 128;
571  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk)/threadsPerCUDABlock);
572  statevec_initDebugStateKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
573  qureg.numAmpsPerChunk,
574  qureg.deviceStateVec.real,
575  qureg.deviceStateVec.imag);
576 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

Referenced by initDebugState().

◆ statevec_initDebugStateKernel()

__global__ void statevec_initDebugStateKernel ( long long int  stateVecSize,
qreal *  stateVecReal,
qreal *  stateVecImag 
)

Definition at line 557 of file QuEST_gpu.cu.

557  {
558  long long int index;
559 
560  index = blockIdx.x*blockDim.x + threadIdx.x;
561  if (index>=stateVecSize) return;
562 
563  stateVecReal[index] = (index*2.0)/10.0;
564  stateVecImag[index] = (index*2.0+1.0)/10.0;
565 }

◆ statevec_initPlusState()

void statevec_initPlusState ( Qureg  qureg)

Definition at line 520 of file QuEST_gpu.cu.

521 {
522  int threadsPerCUDABlock, CUDABlocks;
523  threadsPerCUDABlock = 128;
524  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk)/threadsPerCUDABlock);
525  statevec_initPlusStateKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
526  qureg.numAmpsPerChunk,
527  qureg.deviceStateVec.real,
528  qureg.deviceStateVec.imag);
529 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

Referenced by initPlusState().

◆ statevec_initPlusStateKernel()

__global__ void statevec_initPlusStateKernel ( long long int  stateVecSize,
qreal *  stateVecReal,
qreal *  stateVecImag 
)

Definition at line 509 of file QuEST_gpu.cu.

509  {
510  long long int index;
511 
512  index = blockIdx.x*blockDim.x + threadIdx.x;
513  if (index>=stateVecSize) return;
514 
515  qreal normFactor = 1.0/sqrt((qreal)stateVecSize);
516  stateVecReal[index] = normFactor;
517  stateVecImag[index] = 0.0;
518 }

References qreal.

◆ statevec_initStateFromSingleFile()

int statevec_initStateFromSingleFile ( Qureg *  qureg,
char  filename[200],
QuESTEnv  env 
)

Definition at line 605 of file QuEST_gpu.cu.

605  {
606  long long int chunkSize, stateVecSize;
607  long long int indexInChunk, totalIndex;
608 
609  chunkSize = qureg->numAmpsPerChunk;
610  stateVecSize = chunkSize*qureg->numChunks;
611 
612  qreal *stateVecReal = qureg->stateVec.real;
613  qreal *stateVecImag = qureg->stateVec.imag;
614 
615  FILE *fp;
616  char line[200];
617 
618  fp = fopen(filename, "r");
619  if (fp == NULL)
620  return 0;
621 
622  indexInChunk = 0; totalIndex = 0;
623  while (fgets(line, sizeof(char)*200, fp) != NULL && totalIndex<stateVecSize){
624  if (line[0]!='#'){
625  int chunkId = totalIndex/chunkSize;
626  if (chunkId==qureg->chunkId){
627  # if QuEST_PREC==1
628  sscanf(line, "%f, %f", &(stateVecReal[indexInChunk]),
629  &(stateVecImag[indexInChunk]));
630  # elif QuEST_PREC==2
631  sscanf(line, "%lf, %lf", &(stateVecReal[indexInChunk]),
632  &(stateVecImag[indexInChunk]));
633  # elif QuEST_PREC==4
634  sscanf(line, "%lf, %lf", &(stateVecReal[indexInChunk]),
635  &(stateVecImag[indexInChunk]));
636  # endif
637  indexInChunk += 1;
638  }
639  totalIndex += 1;
640  }
641  }
642  fclose(fp);
643  copyStateToGPU(*qureg);
644 
645  // indicate success
646  return 1;
647 }

References Qureg::chunkId, copyStateToGPU(), Qureg::numAmpsPerChunk, Qureg::numChunks, qreal, and Qureg::stateVec.

Referenced by initStateFromSingleFile().

◆ statevec_initStateOfSingleQubit()

void statevec_initStateOfSingleQubit ( Qureg *  qureg,
int  qubitId,
int  outcome 
)

Initialise the state vector of probability amplitudes such that one qubit is set to 'outcome' and all other qubits are in an equal superposition of zero and one.

Parameters
[in,out]quregobject representing the set of qubits to be initialised
[in]qubitIdid of qubit to set to state 'outcome'
[in]valueof qubit 'qubitId'

Definition at line 596 of file QuEST_gpu.cu.

597 {
598  int threadsPerCUDABlock, CUDABlocks;
599  threadsPerCUDABlock = 128;
600  CUDABlocks = ceil((qreal)(qureg->numAmpsPerChunk)/threadsPerCUDABlock);
601  statevec_initStateOfSingleQubitKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg->numAmpsPerChunk, qureg->deviceStateVec.real, qureg->deviceStateVec.imag, qubitId, outcome);
602 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

Referenced by initStateOfSingleQubit().

◆ statevec_initStateOfSingleQubitKernel()

__global__ void statevec_initStateOfSingleQubitKernel ( long long int  stateVecSize,
qreal *  stateVecReal,
qreal *  stateVecImag,
int  qubitId,
int  outcome 
)

Definition at line 578 of file QuEST_gpu.cu.

578  {
579  long long int index;
580  int bit;
581 
582  index = blockIdx.x*blockDim.x + threadIdx.x;
583  if (index>=stateVecSize) return;
584 
585  qreal normFactor = 1.0/sqrt((qreal)stateVecSize/2);
586  bit = extractBit(qubitId, index);
587  if (bit==outcome) {
588  stateVecReal[index] = normFactor;
589  stateVecImag[index] = 0.0;
590  } else {
591  stateVecReal[index] = 0.0;
592  stateVecImag[index] = 0.0;
593  }
594 }

References extractBit(), and qreal.

◆ statevec_initZeroState()

void statevec_initZeroState ( Qureg  qureg)

Definition at line 498 of file QuEST_gpu.cu.

499 {
500  int threadsPerCUDABlock, CUDABlocks;
501  threadsPerCUDABlock = 128;
502  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk)/threadsPerCUDABlock);
503  statevec_initZeroStateKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
504  qureg.numAmpsPerChunk,
505  qureg.deviceStateVec.real,
506  qureg.deviceStateVec.imag);
507 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

Referenced by initZeroState().

◆ statevec_initZeroStateKernel()

__global__ void statevec_initZeroStateKernel ( long long int  stateVecSize,
qreal *  stateVecReal,
qreal *  stateVecImag 
)

Definition at line 482 of file QuEST_gpu.cu.

482  {
483  long long int index;
484 
485  // initialise the state to |0000..0000>
486  index = blockIdx.x*blockDim.x + threadIdx.x;
487  if (index>=stateVecSize) return;
488  stateVecReal[index] = 0.0;
489  stateVecImag[index] = 0.0;
490 
491  if (index==0){
492  // zero state |0000..0000> has probability 1
493  stateVecReal[0] = 1.0;
494  stateVecImag[0] = 0.0;
495  }
496 }

◆ statevec_multiControlledMultiQubitUnitary()

void statevec_multiControlledMultiQubitUnitary ( Qureg  qureg,
long long int  ctrlMask,
int *  targs,
const int  numTargs,
ComplexMatrixN  u 
)

This calls swapQubitAmps only when it would involve a distributed communication; if the qubit chunks already fit in the node, it operates the unitary direct.

It is already gauranteed here that all target qubits can fit on each node (this is validated in the front-end)

@TODO: refactor so that the 'swap back' isn't performed; instead the qubit locations are updated.

Definition at line 917 of file QuEST_gpu.cu.

918 {
919  int threadsPerCUDABlock = 128;
920  int CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk>>numTargs)/threadsPerCUDABlock);
921 
922  // allocate device space for global {targs} (length: numTargs) and populate
923  int *d_targs;
924  size_t targMemSize = numTargs * sizeof *d_targs;
925  cudaMalloc(&d_targs, targMemSize);
926  cudaMemcpy(d_targs, targs, targMemSize, cudaMemcpyHostToDevice);
927 
928  // flatten out the u.real and u.imag lists
929  int uNumRows = (1 << u.numQubits);
930  qreal* uReFlat = (qreal*) malloc(uNumRows*uNumRows * sizeof *uReFlat);
931  qreal* uImFlat = (qreal*) malloc(uNumRows*uNumRows * sizeof *uImFlat);
932  long long int i = 0;
933  for (int r=0; r < uNumRows; r++)
934  for (int c=0; c < uNumRows; c++) {
935  uReFlat[i] = u.real[r][c];
936  uImFlat[i] = u.imag[r][c];
937  i++;
938  }
939 
940  // allocate device space for global u.real and u.imag (flatten by concatenating rows) and populate
941  qreal* d_uRe;
942  qreal* d_uIm;
943  size_t uMemSize = uNumRows*uNumRows * sizeof *d_uRe; // size of each of d_uRe and d_uIm
944  cudaMalloc(&d_uRe, uMemSize);
945  cudaMalloc(&d_uIm, uMemSize);
946  cudaMemcpy(d_uRe, uReFlat, uMemSize, cudaMemcpyHostToDevice);
947  cudaMemcpy(d_uIm, uImFlat, uMemSize, cudaMemcpyHostToDevice);
948 
949  // allocate device Wspace for thread-local {ampInds}, {reAmps}, {imAmps} (length: 1<<numTargs)
950  long long int *d_ampInds;
951  qreal *d_reAmps;
952  qreal *d_imAmps;
953  size_t gridSize = (size_t) threadsPerCUDABlock * CUDABlocks;
954  int numTargAmps = uNumRows;
955  cudaMalloc(&d_ampInds, numTargAmps*gridSize * sizeof *d_ampInds);
956  cudaMalloc(&d_reAmps, numTargAmps*gridSize * sizeof *d_reAmps);
957  cudaMalloc(&d_imAmps, numTargAmps*gridSize * sizeof *d_imAmps);
958 
959  // call kernel
960  statevec_multiControlledMultiQubitUnitaryKernel<<<CUDABlocks,threadsPerCUDABlock>>>(
961  qureg, ctrlMask, d_targs, numTargs, d_uRe, d_uIm, d_ampInds, d_reAmps, d_imAmps, numTargAmps);
962 
963  // free kernel memory
964  free(uReFlat);
965  free(uImFlat);
966  cudaFree(d_targs);
967  cudaFree(d_uRe);
968  cudaFree(d_uIm);
969  cudaFree(d_ampInds);
970  cudaFree(d_reAmps);
971  cudaFree(d_imAmps);
972 }

References ComplexMatrixN::imag, Qureg::numAmpsPerChunk, ComplexMatrixN::numQubits, qreal, and ComplexMatrixN::real.

Referenced by densmatr_applyMultiQubitKrausSuperoperator(), densmatr_applyTwoQubitKrausSuperoperator(), multiControlledMultiQubitUnitary(), statevec_controlledMultiQubitUnitary(), and statevec_multiQubitUnitary().

◆ statevec_multiControlledMultiQubitUnitaryKernel()

__global__ void statevec_multiControlledMultiQubitUnitaryKernel ( Qureg  qureg,
long long int  ctrlMask,
int *  targs,
int  numTargs,
qreal *  uRe,
qreal *  uIm,
long long int *  ampInds,
qreal *  reAmps,
qreal *  imAmps,
long long int  numTargAmps 
)

Definition at line 858 of file QuEST_gpu.cu.

861 {
862 
863  // decide the amplitudes this thread will modify
864  long long int thisTask = blockIdx.x*blockDim.x + threadIdx.x;
865  long long int numTasks = qureg.numAmpsPerChunk >> numTargs; // kernel called on every 1 in 2^numTargs amplitudes
866  if (thisTask>=numTasks) return;
867 
868  // find this task's start index (where all targs are 0)
869  long long int ind00 = insertZeroBits(thisTask, targs, numTargs);
870 
871  // this task only modifies amplitudes if control qubits are 1 for this state
872  if (ctrlMask && (ctrlMask&ind00) != ctrlMask)
873  return;
874 
875  qreal *reVec = qureg.deviceStateVec.real;
876  qreal *imVec = qureg.deviceStateVec.imag;
877 
878  /*
879  each thread needs:
880  long long int ampInds[numAmps];
881  qreal reAmps[numAmps];
882  qreal imAmps[numAmps];
883  but instead has access to shared arrays, with below stride and offset
884  */
885  size_t stride = gridDim.x*blockDim.x;
886  size_t offset = blockIdx.x*blockDim.x + threadIdx.x;
887 
888  // determine the indices and record values of target amps
889  long long int ind;
890  for (int i=0; i < numTargAmps; i++) {
891 
892  // get global index of current target qubit assignment
893  ind = ind00;
894  for (int t=0; t < numTargs; t++)
895  if (extractBit(t, i))
896  ind = flipBit(ind, targs[t]);
897 
898  ampInds[i*stride+offset] = ind;
899  reAmps [i*stride+offset] = reVec[ind];
900  imAmps [i*stride+offset] = imVec[ind];
901  }
902 
903  // update the amplitudes
904  for (int r=0; r < numTargAmps; r++) {
905  ind = ampInds[r*stride+offset];
906  reVec[ind] = 0;
907  imVec[ind] = 0;
908  for (int c=0; c < numTargAmps; c++) {
909  qreal uReElem = uRe[c + r*numTargAmps];
910  qreal uImElem = uIm[c + r*numTargAmps];
911  reVec[ind] += reAmps[c*stride+offset]*uReElem - imAmps[c*stride+offset]*uImElem;
912  imVec[ind] += reAmps[c*stride+offset]*uImElem + imAmps[c*stride+offset]*uReElem;
913  }
914  }
915 }

References Qureg::deviceStateVec, extractBit(), flipBit(), insertZeroBits(), Qureg::numAmpsPerChunk, and qreal.

◆ statevec_multiControlledPhaseFlip()

void statevec_multiControlledPhaseFlip ( Qureg  qureg,
int *  controlQubits,
int  numControlQubits 
)

Definition at line 1578 of file QuEST_gpu.cu.

1579 {
1580  int threadsPerCUDABlock, CUDABlocks;
1581  long long int mask = getQubitBitMask(controlQubits, numControlQubits);
1582  threadsPerCUDABlock = 128;
1583  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk)/threadsPerCUDABlock);
1584  statevec_multiControlledPhaseFlipKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, mask);
1585 }

References getQubitBitMask(), Qureg::numAmpsPerChunk, and qreal.

Referenced by multiControlledPhaseFlip().

◆ statevec_multiControlledPhaseFlipKernel()

__global__ void statevec_multiControlledPhaseFlipKernel ( Qureg  qureg,
long long int  mask 
)

Definition at line 1560 of file QuEST_gpu.cu.

1561 {
1562  long long int index;
1563  long long int stateVecSize;
1564 
1565  stateVecSize = qureg.numAmpsPerChunk;
1566  qreal *stateVecReal = qureg.deviceStateVec.real;
1567  qreal *stateVecImag = qureg.deviceStateVec.imag;
1568 
1569  index = blockIdx.x*blockDim.x + threadIdx.x;
1570  if (index>=stateVecSize) return;
1571 
1572  if (mask == (mask & index) ){
1573  stateVecReal [index] = - stateVecReal [index];
1574  stateVecImag [index] = - stateVecImag [index];
1575  }
1576 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

◆ statevec_multiControlledPhaseShift()

void statevec_multiControlledPhaseShift ( Qureg  qureg,
int *  controlQubits,
int  numControlQubits,
qreal  angle 
)

Definition at line 1436 of file QuEST_gpu.cu.

1437 {
1438  qreal cosAngle = cos(angle);
1439  qreal sinAngle = sin(angle);
1440 
1441  long long int mask = getQubitBitMask(controlQubits, numControlQubits);
1442 
1443  int threadsPerCUDABlock, CUDABlocks;
1444  threadsPerCUDABlock = 128;
1445  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk)/threadsPerCUDABlock);
1446  statevec_multiControlledPhaseShiftKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, mask, cosAngle, sinAngle);
1447 }

References getQubitBitMask(), Qureg::numAmpsPerChunk, and qreal.

Referenced by multiControlledPhaseShift().

◆ statevec_multiControlledPhaseShiftKernel()

__global__ void statevec_multiControlledPhaseShiftKernel ( Qureg  qureg,
long long int  mask,
qreal  cosAngle,
qreal  sinAngle 
)

Definition at line 1416 of file QuEST_gpu.cu.

1416  {
1417  qreal stateRealLo, stateImagLo;
1418  long long int index;
1419  long long int stateVecSize;
1420 
1421  stateVecSize = qureg.numAmpsPerChunk;
1422  qreal *stateVecReal = qureg.deviceStateVec.real;
1423  qreal *stateVecImag = qureg.deviceStateVec.imag;
1424 
1425  index = blockIdx.x*blockDim.x + threadIdx.x;
1426  if (index>=stateVecSize) return;
1427 
1428  if (mask == (mask & index) ){
1429  stateRealLo = stateVecReal[index];
1430  stateImagLo = stateVecImag[index];
1431  stateVecReal[index] = cosAngle*stateRealLo - sinAngle*stateImagLo;
1432  stateVecImag[index] = sinAngle*stateRealLo + cosAngle*stateImagLo;
1433  }
1434 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

◆ statevec_multiControlledTwoQubitUnitary()

void statevec_multiControlledTwoQubitUnitary ( Qureg  qureg,
long long int  ctrlMask,
const int  q1,
const int  q2,
ComplexMatrix4  u 
)

This calls swapQubitAmps only when it would involve a distributed communication; if the qubit chunks already fit in the node, it operates the unitary direct.

Note the order of q1 and q2 in the call to twoQubitUnitaryLocal is important.

@TODO: refactor so that the 'swap back' isn't performed; instead the qubit locations are updated. @TODO: the double swap (q1,q2 to 0,1) may be possible simultaneously by a bespoke swap routine.

Definition at line 1050 of file QuEST_gpu.cu.

1051 {
1052  int threadsPerCUDABlock = 128;
1053  int CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk>>2)/threadsPerCUDABlock); // one kernel eval for every 4 amplitudes
1054  statevec_multiControlledTwoQubitUnitaryKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, ctrlMask, q1, q2, argifyMatrix4(u));
1055 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by densmatr_applyKrausSuperoperator(), multiControlledTwoQubitUnitary(), statevec_controlledTwoQubitUnitary(), and statevec_twoQubitUnitary().

◆ statevec_multiControlledTwoQubitUnitaryKernel()

__global__ void statevec_multiControlledTwoQubitUnitaryKernel ( Qureg  qureg,
long long int  ctrlMask,
const int  q1,
const int  q2,
ArgMatrix4  u 
)

Definition at line 974 of file QuEST_gpu.cu.

974  {
975 
976  // decide the 4 amplitudes this thread will modify
977  long long int thisTask = blockIdx.x*blockDim.x + threadIdx.x;
978  long long int numTasks = qureg.numAmpsPerChunk >> 2; // kernel called on every 1 in 4 amplitudes
979  if (thisTask>=numTasks) return;
980 
981  qreal *reVec = qureg.deviceStateVec.real;
982  qreal *imVec = qureg.deviceStateVec.imag;
983 
984  // find indices of amplitudes to modify (treat q1 as the least significant bit)
985  long long int ind00, ind01, ind10, ind11;
986  ind00 = insertTwoZeroBits(thisTask, q1, q2);
987 
988  // modify only if control qubits are 1 for this state
989  if (ctrlMask && (ctrlMask&ind00) != ctrlMask)
990  return;
991 
992  ind01 = flipBit(ind00, q1);
993  ind10 = flipBit(ind00, q2);
994  ind11 = flipBit(ind01, q2);
995 
996  // extract statevec amplitudes
997  qreal re00, re01, re10, re11;
998  qreal im00, im01, im10, im11;
999  re00 = reVec[ind00]; im00 = imVec[ind00];
1000  re01 = reVec[ind01]; im01 = imVec[ind01];
1001  re10 = reVec[ind10]; im10 = imVec[ind10];
1002  re11 = reVec[ind11]; im11 = imVec[ind11];
1003 
1004  // apply u * {amp00, amp01, amp10, amp11}
1005  reVec[ind00] =
1006  u.r0c0.real*re00 - u.r0c0.imag*im00 +
1007  u.r0c1.real*re01 - u.r0c1.imag*im01 +
1008  u.r0c2.real*re10 - u.r0c2.imag*im10 +
1009  u.r0c3.real*re11 - u.r0c3.imag*im11;
1010  imVec[ind00] =
1011  u.r0c0.imag*re00 + u.r0c0.real*im00 +
1012  u.r0c1.imag*re01 + u.r0c1.real*im01 +
1013  u.r0c2.imag*re10 + u.r0c2.real*im10 +
1014  u.r0c3.imag*re11 + u.r0c3.real*im11;
1015 
1016  reVec[ind01] =
1017  u.r1c0.real*re00 - u.r1c0.imag*im00 +
1018  u.r1c1.real*re01 - u.r1c1.imag*im01 +
1019  u.r1c2.real*re10 - u.r1c2.imag*im10 +
1020  u.r1c3.real*re11 - u.r1c3.imag*im11;
1021  imVec[ind01] =
1022  u.r1c0.imag*re00 + u.r1c0.real*im00 +
1023  u.r1c1.imag*re01 + u.r1c1.real*im01 +
1024  u.r1c2.imag*re10 + u.r1c2.real*im10 +
1025  u.r1c3.imag*re11 + u.r1c3.real*im11;
1026 
1027  reVec[ind10] =
1028  u.r2c0.real*re00 - u.r2c0.imag*im00 +
1029  u.r2c1.real*re01 - u.r2c1.imag*im01 +
1030  u.r2c2.real*re10 - u.r2c2.imag*im10 +
1031  u.r2c3.real*re11 - u.r2c3.imag*im11;
1032  imVec[ind10] =
1033  u.r2c0.imag*re00 + u.r2c0.real*im00 +
1034  u.r2c1.imag*re01 + u.r2c1.real*im01 +
1035  u.r2c2.imag*re10 + u.r2c2.real*im10 +
1036  u.r2c3.imag*re11 + u.r2c3.real*im11;
1037 
1038  reVec[ind11] =
1039  u.r3c0.real*re00 - u.r3c0.imag*im00 +
1040  u.r3c1.real*re01 - u.r3c1.imag*im01 +
1041  u.r3c2.real*re10 - u.r3c2.imag*im10 +
1042  u.r3c3.real*re11 - u.r3c3.imag*im11;
1043  imVec[ind11] =
1044  u.r3c0.imag*re00 + u.r3c0.real*im00 +
1045  u.r3c1.imag*re01 + u.r3c1.real*im01 +
1046  u.r3c2.imag*re10 + u.r3c2.real*im10 +
1047  u.r3c3.imag*re11 + u.r3c3.real*im11;
1048 }

References Qureg::deviceStateVec, flipBit(), insertTwoZeroBits(), Qureg::numAmpsPerChunk, and qreal.

◆ statevec_multiControlledUnitary()

void statevec_multiControlledUnitary ( Qureg  qureg,
long long int  ctrlQubitsMask,
long long int  ctrlFlipMask,
const int  targetQubit,
ComplexMatrix2  u 
)

Definition at line 1183 of file QuEST_gpu.cu.

1187  {
1188  int threadsPerCUDABlock = 128;
1189  int CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk>>1)/threadsPerCUDABlock);
1190  statevec_multiControlledUnitaryKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
1191  qureg, ctrlQubitsMask, ctrlFlipMask, targetQubit, argifyMatrix2(u));
1192 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by multiControlledUnitary(), and multiStateControlledUnitary().

◆ statevec_multiControlledUnitaryKernel()

__global__ void statevec_multiControlledUnitaryKernel ( Qureg  qureg,
long long int  ctrlQubitsMask,
long long int  ctrlFlipMask,
const int  targetQubit,
ArgMatrix2  u 
)

fix – no necessary for GPU version

Definition at line 1123 of file QuEST_gpu.cu.

1127  {
1128  // ----- sizes
1129  long long int sizeBlock, // size of blocks
1130  sizeHalfBlock; // size of blocks halved
1131  // ----- indices
1132  long long int thisBlock, // current block
1133  indexUp,indexLo; // current index and corresponding index in lower half block
1134 
1135  // ----- temp variables
1136  qreal stateRealUp,stateRealLo, // storage for previous state values
1137  stateImagUp,stateImagLo; // (used in updates)
1138  // ----- temp variables
1139  long long int thisTask; // task based approach for expose loop with small granularity
1140  const long long int numTasks=qureg.numAmpsPerChunk>>1;
1141 
1142 
1143  sizeHalfBlock = 1LL << targetQubit; // size of blocks halved
1144  sizeBlock = 2LL * sizeHalfBlock; // size of blocks
1145 
1146  // ---------------------------------------------------------------- //
1147  // rotate //
1148  // ---------------------------------------------------------------- //
1149 
1151  qreal *stateVecReal = qureg.deviceStateVec.real;
1152  qreal *stateVecImag = qureg.deviceStateVec.imag;
1153 
1154  thisTask = blockIdx.x*blockDim.x + threadIdx.x;
1155  if (thisTask>=numTasks) return;
1156 
1157  thisBlock = thisTask / sizeHalfBlock;
1158  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
1159  indexLo = indexUp + sizeHalfBlock;
1160 
1161  if (ctrlQubitsMask == (ctrlQubitsMask & (indexUp ^ ctrlFlipMask))) {
1162  // store current state vector values in temp variables
1163  stateRealUp = stateVecReal[indexUp];
1164  stateImagUp = stateVecImag[indexUp];
1165 
1166  stateRealLo = stateVecReal[indexLo];
1167  stateImagLo = stateVecImag[indexLo];
1168 
1169  // state[indexUp] = u00 * state[indexUp] + u01 * state[indexLo]
1170  stateVecReal[indexUp] = u.r0c0.real*stateRealUp - u.r0c0.imag*stateImagUp
1171  + u.r0c1.real*stateRealLo - u.r0c1.imag*stateImagLo;
1172  stateVecImag[indexUp] = u.r0c0.real*stateImagUp + u.r0c0.imag*stateRealUp
1173  + u.r0c1.real*stateImagLo + u.r0c1.imag*stateRealLo;
1174 
1175  // state[indexLo] = u10 * state[indexUp] + u11 * state[indexLo]
1176  stateVecReal[indexLo] = u.r1c0.real*stateRealUp - u.r1c0.imag*stateImagUp
1177  + u.r1c1.real*stateRealLo - u.r1c1.imag*stateImagLo;
1178  stateVecImag[indexLo] = u.r1c0.real*stateImagUp + u.r1c0.imag*stateRealUp
1179  + u.r1c1.real*stateImagLo + u.r1c1.imag*stateRealLo;
1180  }
1181 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

◆ statevec_multiRotateZ()

void statevec_multiRotateZ ( Qureg  qureg,
long long int  mask,
qreal  angle 
)

Definition at line 1466 of file QuEST_gpu.cu.

1467 {
1468  qreal cosAngle = cos(angle/2.0);
1469  qreal sinAngle = sin(angle/2.0);
1470 
1471  int threadsPerCUDABlock, CUDABlocks;
1472  threadsPerCUDABlock = 128;
1473  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk)/threadsPerCUDABlock);
1474  statevec_multiRotateZKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, mask, cosAngle, sinAngle);
1475 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by multiRotateZ(), and statevec_multiRotatePauli().

◆ statevec_multiRotateZKernel()

__global__ void statevec_multiRotateZKernel ( Qureg  qureg,
long long int  mask,
qreal  cosAngle,
qreal  sinAngle 
)

Definition at line 1449 of file QuEST_gpu.cu.

1449  {
1450 
1451  long long int stateVecSize = qureg.numAmpsPerChunk;
1452  long long int index = blockIdx.x*blockDim.x + threadIdx.x;
1453  if (index>=stateVecSize) return;
1454 
1455  qreal *stateVecReal = qureg.deviceStateVec.real;
1456  qreal *stateVecImag = qureg.deviceStateVec.imag;
1457 
1458  int fac = getBitMaskParity(mask & index)? -1 : 1;
1459  qreal stateReal = stateVecReal[index];
1460  qreal stateImag = stateVecImag[index];
1461 
1462  stateVecReal[index] = cosAngle*stateReal + fac * sinAngle*stateImag;
1463  stateVecImag[index] = - fac * sinAngle*stateReal + cosAngle*stateImag;
1464 }

References Qureg::deviceStateVec, getBitMaskParity(), Qureg::numAmpsPerChunk, and qreal.

◆ statevec_pauliX()

void statevec_pauliX ( Qureg  qureg,
const int  targetQubit 
)

Definition at line 1238 of file QuEST_gpu.cu.

1239 {
1240  int threadsPerCUDABlock, CUDABlocks;
1241  threadsPerCUDABlock = 128;
1242  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk>>1)/threadsPerCUDABlock);
1243  statevec_pauliXKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, targetQubit);
1244 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by pauliX(), and statevec_applyPauliProd().

◆ statevec_pauliXKernel()

__global__ void statevec_pauliXKernel ( Qureg  qureg,
const int  targetQubit 
)

fix – no necessary for GPU version

Definition at line 1194 of file QuEST_gpu.cu.

1194  {
1195  // ----- sizes
1196  long long int sizeBlock, // size of blocks
1197  sizeHalfBlock; // size of blocks halved
1198  // ----- indices
1199  long long int thisBlock, // current block
1200  indexUp,indexLo; // current index and corresponding index in lower half block
1201 
1202  // ----- temp variables
1203  qreal stateRealUp, // storage for previous state values
1204  stateImagUp; // (used in updates)
1205  // ----- temp variables
1206  long long int thisTask; // task based approach for expose loop with small granularity
1207  const long long int numTasks=qureg.numAmpsPerChunk>>1;
1208 
1209  sizeHalfBlock = 1LL << targetQubit; // size of blocks halved
1210  sizeBlock = 2LL * sizeHalfBlock; // size of blocks
1211 
1212  // ---------------------------------------------------------------- //
1213  // rotate //
1214  // ---------------------------------------------------------------- //
1215 
1217  qreal *stateVecReal = qureg.deviceStateVec.real;
1218  qreal *stateVecImag = qureg.deviceStateVec.imag;
1219 
1220  thisTask = blockIdx.x*blockDim.x + threadIdx.x;
1221  if (thisTask>=numTasks) return;
1222 
1223  thisBlock = thisTask / sizeHalfBlock;
1224  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
1225  indexLo = indexUp + sizeHalfBlock;
1226 
1227  // store current state vector values in temp variables
1228  stateRealUp = stateVecReal[indexUp];
1229  stateImagUp = stateVecImag[indexUp];
1230 
1231  stateVecReal[indexUp] = stateVecReal[indexLo];
1232  stateVecImag[indexUp] = stateVecImag[indexLo];
1233 
1234  stateVecReal[indexLo] = stateRealUp;
1235  stateVecImag[indexLo] = stateImagUp;
1236 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

◆ statevec_pauliY()

void statevec_pauliY ( Qureg  qureg,
const int  targetQubit 
)

Definition at line 1271 of file QuEST_gpu.cu.

1272 {
1273  int threadsPerCUDABlock, CUDABlocks;
1274  threadsPerCUDABlock = 128;
1275  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk>>1)/threadsPerCUDABlock);
1276  statevec_pauliYKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, targetQubit, 1);
1277 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by pauliY(), and statevec_applyPauliProd().

◆ statevec_pauliYConj()

void statevec_pauliYConj ( Qureg  qureg,
const int  targetQubit 
)

Definition at line 1279 of file QuEST_gpu.cu.

1280 {
1281  int threadsPerCUDABlock, CUDABlocks;
1282  threadsPerCUDABlock = 128;
1283  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk>>1)/threadsPerCUDABlock);
1284  statevec_pauliYKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, targetQubit, -1);
1285 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by pauliY().

◆ statevec_pauliYKernel()

__global__ void statevec_pauliYKernel ( Qureg  qureg,
const int  targetQubit,
const int  conjFac 
)

Definition at line 1246 of file QuEST_gpu.cu.

1246  {
1247 
1248  long long int sizeHalfBlock = 1LL << targetQubit;
1249  long long int sizeBlock = 2LL * sizeHalfBlock;
1250  long long int numTasks = qureg.numAmpsPerChunk >> 1;
1251  long long int thisTask = blockIdx.x*blockDim.x + threadIdx.x;
1252  if (thisTask>=numTasks) return;
1253 
1254  long long int thisBlock = thisTask / sizeHalfBlock;
1255  long long int indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
1256  long long int indexLo = indexUp + sizeHalfBlock;
1257  qreal stateRealUp, stateImagUp;
1258 
1259  qreal *stateVecReal = qureg.deviceStateVec.real;
1260  qreal *stateVecImag = qureg.deviceStateVec.imag;
1261  stateRealUp = stateVecReal[indexUp];
1262  stateImagUp = stateVecImag[indexUp];
1263 
1264  // update under +-{{0, -i}, {i, 0}}
1265  stateVecReal[indexUp] = conjFac * stateVecImag[indexLo];
1266  stateVecImag[indexUp] = conjFac * -stateVecReal[indexLo];
1267  stateVecReal[indexLo] = conjFac * -stateImagUp;
1268  stateVecImag[indexLo] = conjFac * stateRealUp;
1269 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

◆ statevec_phaseShiftByTerm()

void statevec_phaseShiftByTerm ( Qureg  qureg,
const int  targetQubit,
Complex  term 
)

Definition at line 1369 of file QuEST_gpu.cu.

1370 {
1371  qreal cosAngle = term.real;
1372  qreal sinAngle = term.imag;
1373 
1374  int threadsPerCUDABlock, CUDABlocks;
1375  threadsPerCUDABlock = 128;
1376  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk>>1)/threadsPerCUDABlock);
1377  statevec_phaseShiftByTermKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, targetQubit, cosAngle, sinAngle);
1378 }

References Complex::imag, Qureg::numAmpsPerChunk, qreal, and Complex::real.

Referenced by statevec_pauliZ(), statevec_phaseShift(), statevec_sGate(), statevec_sGateConj(), statevec_tGate(), and statevec_tGateConj().

◆ statevec_phaseShiftByTermKernel()

__global__ void statevec_phaseShiftByTermKernel ( Qureg  qureg,
const int  targetQubit,
qreal  cosAngle,
qreal  sinAngle 
)

Definition at line 1341 of file QuEST_gpu.cu.

1341  {
1342 
1343  long long int sizeBlock, sizeHalfBlock;
1344  long long int thisBlock, indexUp,indexLo;
1345 
1346  qreal stateRealLo, stateImagLo;
1347  long long int thisTask;
1348  const long long int numTasks = qureg.numAmpsPerChunk >> 1;
1349 
1350  sizeHalfBlock = 1LL << targetQubit;
1351  sizeBlock = 2LL * sizeHalfBlock;
1352 
1353  qreal *stateVecReal = qureg.deviceStateVec.real;
1354  qreal *stateVecImag = qureg.deviceStateVec.imag;
1355 
1356  thisTask = blockIdx.x*blockDim.x + threadIdx.x;
1357  if (thisTask>=numTasks) return;
1358  thisBlock = thisTask / sizeHalfBlock;
1359  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
1360  indexLo = indexUp + sizeHalfBlock;
1361 
1362  stateRealLo = stateVecReal[indexLo];
1363  stateImagLo = stateVecImag[indexLo];
1364 
1365  stateVecReal[indexLo] = cosAngle*stateRealLo - sinAngle*stateImagLo;
1366  stateVecImag[indexLo] = sinAngle*stateRealLo + cosAngle*stateImagLo;
1367 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

◆ statevec_reportStateToScreen()

void statevec_reportStateToScreen ( Qureg  qureg,
QuESTEnv  env,
int  reportRank 
)

Print the current state vector of probability amplitudes for a set of qubits to standard out.

For debugging purposes. Each rank should print output serially. Only print output for systems <= 5 qubits

Definition at line 421 of file QuEST_gpu.cu.

421  {
422  long long int index;
423  int rank;
424  copyStateFromGPU(qureg);
425  if (qureg.numQubitsInStateVec<=5){
426  for (rank=0; rank<qureg.numChunks; rank++){
427  if (qureg.chunkId==rank){
428  if (reportRank) {
429  printf("Reporting state from rank %d [\n", qureg.chunkId);
430  //printf("\trank, index, real, imag\n");
431  printf("real, imag\n");
432  } else if (rank==0) {
433  printf("Reporting state [\n");
434  printf("real, imag\n");
435  }
436 
437  for(index=0; index<qureg.numAmpsPerChunk; index++){
438  printf(REAL_STRING_FORMAT ", " REAL_STRING_FORMAT "\n", qureg.stateVec.real[index], qureg.stateVec.imag[index]);
439  }
440  if (reportRank || rank==qureg.numChunks-1) printf("]\n");
441  }
442  syncQuESTEnv(env);
443  }
444  }
445 }

References Qureg::chunkId, copyStateFromGPU(), Qureg::numAmpsPerChunk, Qureg::numChunks, Qureg::numQubitsInStateVec, Qureg::stateVec, and syncQuESTEnv().

Referenced by reportStateToScreen().

◆ statevec_setAmps()

void statevec_setAmps ( Qureg  qureg,
long long int  startInd,
qreal *  reals,
qreal *  imags,
long long int  numAmps 
)

Definition at line 153 of file QuEST_gpu.cu.

153  {
154 
155  cudaDeviceSynchronize();
156  cudaMemcpy(
157  qureg.deviceStateVec.real + startInd,
158  reals,
159  numAmps * sizeof(*(qureg.deviceStateVec.real)),
160  cudaMemcpyHostToDevice);
161  cudaMemcpy(
162  qureg.deviceStateVec.imag + startInd,
163  imags,
164  numAmps * sizeof(*(qureg.deviceStateVec.real)),
165  cudaMemcpyHostToDevice);
166 }

References Qureg::deviceStateVec.

Referenced by initStateFromAmps(), setAmps(), and setDensityAmps().

◆ statevec_setWeightedQureg()

void statevec_setWeightedQureg ( Complex  fac1,
Qureg  qureg1,
Complex  fac2,
Qureg  qureg2,
Complex  facOut,
Qureg  out 
)

Definition at line 2850 of file QuEST_gpu.cu.

2850  {
2851 
2852  long long int numAmpsToVisit = qureg1.numAmpsPerChunk;
2853 
2854  int threadsPerCUDABlock, CUDABlocks;
2855  threadsPerCUDABlock = 128;
2856  CUDABlocks = ceil(numAmpsToVisit / (qreal) threadsPerCUDABlock);
2857  statevec_setWeightedQuregKernel<<<CUDABlocks, threadsPerCUDABlock>>>(
2858  fac1, qureg1, fac2, qureg2, facOut, out
2859  );
2860 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by setWeightedQureg(), and statevec_applyPauliSum().

◆ statevec_setWeightedQuregKernel()

__global__ void statevec_setWeightedQuregKernel ( Complex  fac1,
Qureg  qureg1,
Complex  fac2,
Qureg  qureg2,
Complex  facOut,
Qureg  out 
)

Definition at line 2818 of file QuEST_gpu.cu.

2818  {
2819 
2820  long long int ampInd = blockIdx.x*blockDim.x + threadIdx.x;
2821  long long int numAmpsToVisit = qureg1.numAmpsPerChunk;
2822  if (ampInd >= numAmpsToVisit) return;
2823 
2824  qreal *vecRe1 = qureg1.deviceStateVec.real;
2825  qreal *vecIm1 = qureg1.deviceStateVec.imag;
2826  qreal *vecRe2 = qureg2.deviceStateVec.real;
2827  qreal *vecIm2 = qureg2.deviceStateVec.imag;
2828  qreal *vecReOut = out.deviceStateVec.real;
2829  qreal *vecImOut = out.deviceStateVec.imag;
2830 
2831  qreal facRe1 = fac1.real;
2832  qreal facIm1 = fac1.imag;
2833  qreal facRe2 = fac2.real;
2834  qreal facIm2 = fac2.imag;
2835  qreal facReOut = facOut.real;
2836  qreal facImOut = facOut.imag;
2837 
2838  qreal re1,im1, re2,im2, reOut,imOut;
2839  long long int index = ampInd;
2840 
2841  re1 = vecRe1[index]; im1 = vecIm1[index];
2842  re2 = vecRe2[index]; im2 = vecIm2[index];
2843  reOut = vecReOut[index];
2844  imOut = vecImOut[index];
2845 
2846  vecReOut[index] = (facReOut*reOut - facImOut*imOut) + (facRe1*re1 - facIm1*im1) + (facRe2*re2 - facIm2*im2);
2847  vecImOut[index] = (facReOut*imOut + facImOut*reOut) + (facRe1*im1 + facIm1*re1) + (facRe2*im2 + facIm2*re2);
2848 }

References Qureg::deviceStateVec, Complex::imag, Qureg::numAmpsPerChunk, qreal, and Complex::real.

◆ statevec_swapQubitAmps()

void statevec_swapQubitAmps ( Qureg  qureg,
int  qb1,
int  qb2 
)

Definition at line 1613 of file QuEST_gpu.cu.

1614 {
1615  int threadsPerCUDABlock, CUDABlocks;
1616  threadsPerCUDABlock = 128;
1617  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk>>2)/threadsPerCUDABlock);
1618  statevec_swapQubitAmpsKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, qb1, qb2);
1619 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by swapGate().

◆ statevec_swapQubitAmpsKernel()

__global__ void statevec_swapQubitAmpsKernel ( Qureg  qureg,
int  qb1,
int  qb2 
)

Definition at line 1587 of file QuEST_gpu.cu.

1587  {
1588 
1589  qreal *reVec = qureg.deviceStateVec.real;
1590  qreal *imVec = qureg.deviceStateVec.imag;
1591 
1592  long long int numTasks = qureg.numAmpsPerChunk >> 2; // each iteration updates 2 amps and skips 2 amps
1593  long long int thisTask = blockIdx.x*blockDim.x + threadIdx.x;
1594  if (thisTask>=numTasks) return;
1595 
1596  long long int ind00, ind01, ind10;
1597  qreal re01, re10, im01, im10;
1598 
1599  // determine ind00 of |..0..0..>, |..0..1..> and |..1..0..>
1600  ind00 = insertTwoZeroBits(thisTask, qb1, qb2);
1601  ind01 = flipBit(ind00, qb1);
1602  ind10 = flipBit(ind00, qb2);
1603 
1604  // extract statevec amplitudes
1605  re01 = reVec[ind01]; im01 = imVec[ind01];
1606  re10 = reVec[ind10]; im10 = imVec[ind10];
1607 
1608  // swap 01 and 10 amps
1609  reVec[ind01] = re10; reVec[ind10] = re01;
1610  imVec[ind01] = im10; imVec[ind10] = im01;
1611 }

References Qureg::deviceStateVec, flipBit(), insertTwoZeroBits(), Qureg::numAmpsPerChunk, and qreal.

◆ statevec_unitary()

void statevec_unitary ( Qureg  qureg,
const int  targetQubit,
ComplexMatrix2  u 
)

Definition at line 850 of file QuEST_gpu.cu.

851 {
852  int threadsPerCUDABlock, CUDABlocks;
853  threadsPerCUDABlock = 128;
854  CUDABlocks = ceil((qreal)(qureg.numAmpsPerChunk>>1)/threadsPerCUDABlock);
855  statevec_unitaryKernel<<<CUDABlocks, threadsPerCUDABlock>>>(qureg, targetQubit, argifyMatrix2(u));
856 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by unitary().

◆ statevec_unitaryKernel()

__global__ void statevec_unitaryKernel ( Qureg  qureg,
const int  targetQubit,
ArgMatrix2  u 
)

fix – no necessary for GPU version

Definition at line 797 of file QuEST_gpu.cu.

797  {
798  // ----- sizes
799  long long int sizeBlock, // size of blocks
800  sizeHalfBlock; // size of blocks halved
801  // ----- indices
802  long long int thisBlock, // current block
803  indexUp,indexLo; // current index and corresponding index in lower half block
804 
805  // ----- temp variables
806  qreal stateRealUp,stateRealLo, // storage for previous state values
807  stateImagUp,stateImagLo; // (used in updates)
808  // ----- temp variables
809  long long int thisTask; // task based approach for expose loop with small granularity
810  const long long int numTasks=qureg.numAmpsPerChunk>>1;
811 
812  sizeHalfBlock = 1LL << targetQubit; // size of blocks halved
813  sizeBlock = 2LL * sizeHalfBlock; // size of blocks
814 
815  // ---------------------------------------------------------------- //
816  // rotate //
817  // ---------------------------------------------------------------- //
818 
820  qreal *stateVecReal = qureg.deviceStateVec.real;
821  qreal *stateVecImag = qureg.deviceStateVec.imag;
822 
823  thisTask = blockIdx.x*blockDim.x + threadIdx.x;
824  if (thisTask>=numTasks) return;
825 
826  thisBlock = thisTask / sizeHalfBlock;
827  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
828  indexLo = indexUp + sizeHalfBlock;
829 
830  // store current state vector values in temp variables
831  stateRealUp = stateVecReal[indexUp];
832  stateImagUp = stateVecImag[indexUp];
833 
834  stateRealLo = stateVecReal[indexLo];
835  stateImagLo = stateVecImag[indexLo];
836 
837  // state[indexUp] = u00 * state[indexUp] + u01 * state[indexLo]
838  stateVecReal[indexUp] = u.r0c0.real*stateRealUp - u.r0c0.imag*stateImagUp
839  + u.r0c1.real*stateRealLo - u.r0c1.imag*stateImagLo;
840  stateVecImag[indexUp] = u.r0c0.real*stateImagUp + u.r0c0.imag*stateRealUp
841  + u.r0c1.real*stateImagLo + u.r0c1.imag*stateRealLo;
842 
843  // state[indexLo] = u10 * state[indexUp] + u11 * state[indexLo]
844  stateVecReal[indexLo] = u.r1c0.real*stateRealUp - u.r1c0.imag*stateImagUp
845  + u.r1c1.real*stateRealLo - u.r1c1.imag*stateImagLo;
846  stateVecImag[indexLo] = u.r1c0.real*stateImagUp + u.r1c0.imag*stateRealUp
847  + u.r1c1.real*stateImagLo + u.r1c1.imag*stateRealLo;
848 }

References Qureg::deviceStateVec, Qureg::numAmpsPerChunk, and qreal.

◆ swapDouble()

void swapDouble ( qreal **  a,
qreal **  b 
)

Definition at line 1857 of file QuEST_gpu.cu.

1857  {
1858  qreal *temp;
1859  temp = *a;
1860  *a = *b;
1861  *b = temp;
1862 }

References qreal.

Referenced by densmatr_calcFidelity(), densmatr_calcHilbertSchmidtDistance(), densmatr_calcInnerProduct(), densmatr_calcPurity(), densmatr_findProbabilityOfZero(), statevec_calcInnerProduct(), and statevec_findProbabilityOfZero().

log2Int
__device__ __host__ unsigned int log2Int(unsigned int x)
Definition: QuEST_gpu.cu:1725
copyStateFromGPU
void copyStateFromGPU(Qureg qureg)
In GPU mode, this copies the state-vector (or density matrix) from GPU memory (qureg....
Definition: QuEST_gpu.cu:407
syncQuESTEnv
void syncQuESTEnv(QuESTEnv env)
Guarantees that all code up to the given point has been executed on all nodes (if running in distribu...
Definition: QuEST_gpu.cu:369
QuESTEnv::rank
int rank
Definition: QuEST.h:201
copySharedReduceBlock
__global__ void copySharedReduceBlock(qreal *arrayIn, qreal *reducedArray, int length)
Definition: QuEST_gpu.cu:1751
swapDouble
void swapDouble(qreal **a, qreal **b)
Definition: QuEST_gpu.cu:1857
Qureg::pairStateVec
ComplexArray pairStateVec
Temporary storage for a chunk of the state vector received from another process in the MPI version.
Definition: QuEST.h:181
Qureg::numChunks
int numChunks
Number of chunks the state vector is broken up into – the number of MPI processes used.
Definition: QuEST.h:176
insertZeroBits
__forceinline__ __device__ long long int insertZeroBits(long long int number, int *inds, int numInds)
Definition: QuEST_gpu.cu:112
densmatr_oneQubitDegradeOffDiagonal
void densmatr_oneQubitDegradeOffDiagonal(Qureg qureg, const int targetQubit, qreal dephFac)
Definition: QuEST_gpu.cu:2554
extractBit
__forceinline__ __device__ int extractBit(int locationOfBitFromRight, long long int theEncodedNumber)
Definition: QuEST_gpu.cu:82
Qureg::deviceStateVec
ComplexArray deviceStateVec
Storage for wavefunction amplitudes in the GPU version.
Definition: QuEST.h:184
qreal
#define qreal
Qureg::numQubitsInStateVec
int numQubitsInStateVec
Number of qubits in the state-vector - this is double the number represented for mixed states.
Definition: QuEST.h:167
densmatr_findProbabilityOfZero
qreal densmatr_findProbabilityOfZero(Qureg qureg, const int measureQubit)
Definition: QuEST_gpu.cu:1864
insertTwoZeroBits
__forceinline__ __device__ long long int insertTwoZeroBits(long long int number, int bit1, int bit2)
Definition: QuEST_gpu.cu:106
getQubitBitMask
long long int getQubitBitMask(int *qubits, const int numQubits)
Definition: QuEST_common.c:43
Qureg::chunkId
int chunkId
The position of the chunk of the state vector held by this process in the full state vector.
Definition: QuEST.h:174
getBitMaskParity
__forceinline__ __device__ int getBitMaskParity(long long int mask)
Definition: QuEST_gpu.cu:86
Qureg::numAmpsPerChunk
long long int numAmpsPerChunk
Number of probability amplitudes held in stateVec by this process In the non-MPI version,...
Definition: QuEST.h:170
copyStateToGPU
void copyStateToGPU(Qureg qureg)
In GPU mode, this copies the state-vector (or density matrix) from RAM (qureg.stateVec) to VRAM / GPU...
Definition: QuEST_gpu.cu:397
REDUCE_SHARED_SIZE
#define REDUCE_SHARED_SIZE
Definition: QuEST_gpu.cu:19
QuESTEnv::numRanks
int numRanks
Definition: QuEST.h:202
ComplexMatrixN::real
qreal ** real
Definition: QuEST.h:139
Qureg::secondLevelReduction
qreal * secondLevelReduction
Definition: QuEST.h:186
flipBit
__forceinline__ __device__ long long int flipBit(long long int number, int bitInd)
Definition: QuEST_gpu.cu:95
densmatr_mixTwoQubitDephasing
void densmatr_mixTwoQubitDephasing(Qureg qureg, int qubit1, int qubit2, qreal dephase)
Definition: QuEST_gpu.cu:2613
ComplexMatrixN::imag
qreal ** imag
Definition: QuEST.h:140
Qureg::stateVec
ComplexArray stateVec
Computational state amplitudes - a subset thereof in the MPI version.
Definition: QuEST.h:179
Qureg::isDensityMatrix
int isDensityMatrix
Whether this instance is a density-state representation.
Definition: QuEST.h:163
ComplexMatrixN::numQubits
int numQubits
Definition: QuEST.h:138
statevec_findProbabilityOfZero
qreal statevec_findProbabilityOfZero(Qureg qureg, const int measureQubit)
Definition: QuEST_gpu.cu:1912
Qureg::numQubitsRepresented
int numQubitsRepresented
The number of qubits represented in either the state-vector or density matrix.
Definition: QuEST.h:165
Qureg::numAmpsTotal
long long int numAmpsTotal
Total number of amplitudes, which are possibly distributed among machines.
Definition: QuEST.h:172
Complex::real
qreal real
Definition: QuEST.h:105
reduceBlock
__device__ void reduceBlock(qreal *arrayIn, qreal *reducedArray, int length)
Definition: QuEST_gpu.cu:1732
Complex::imag
qreal imag
Definition: QuEST.h:106
insertZeroBit
__forceinline__ __device__ long long int insertZeroBit(long long int number, int index)
Definition: QuEST_gpu.cu:99
Complex
Represents one complex number.
Definition: QuEST.h:103
Qureg::firstLevelReduction
qreal * firstLevelReduction
Storage for reduction of probabilities on GPU.
Definition: QuEST.h:186
densmatr_mixDephasing
void densmatr_mixDephasing(Qureg qureg, const int targetQubit, qreal dephase)
Definition: QuEST_gpu.cu:2574