Functions
QuEST_cpu_internal.h File Reference
#include "QuEST_precision.h"

Go to the source code of this file.

Functions

void densmatr_applyDiagonalOpLocal (Qureg qureg, DiagonalOp op)
 
Complex densmatr_calcExpecDiagonalOpLocal (Qureg qureg, DiagonalOp op)
 
qreal densmatr_calcFidelityLocal (Qureg qureg, Qureg pureState)
 computes a few dens-columns-worth of (vec^*T) dens * vec More...
 
qreal densmatr_calcHilbertSchmidtDistanceSquaredLocal (Qureg a, Qureg b)
 computes Tr((a-b) conjTrans(a-b)) = sum of abs values of (a-b) More...
 
qreal densmatr_calcInnerProductLocal (Qureg a, Qureg b)
 computes Tr(conjTrans(a) b) = sum of (a_ij^* b_ij) More...
 
void densmatr_calcProbOfAllOutcomesLocal (qreal *retProbs, Qureg qureg, int *qubits, int numQubits)
 
qreal densmatr_calcPurityLocal (Qureg qureg)
 
qreal densmatr_findProbabilityOfZeroLocal (Qureg qureg, int measureQubit)
 
void densmatr_initPureStateLocal (Qureg targetQureg, Qureg copyQureg)
 
void densmatr_mixDampingDistributed (Qureg qureg, int targetQubit, qreal damping)
 
void densmatr_mixDampingLocal (Qureg qureg, int targetQubit, qreal damping)
 
void densmatr_mixDepolarisingDistributed (Qureg qureg, int targetQubit, qreal depolLevel)
 
void densmatr_mixDepolarisingLocal (Qureg qureg, int targetQubit, qreal depolLevel)
 
void densmatr_mixTwoQubitDepolarisingDistributed (Qureg qureg, int targetQubit, int qubit2, qreal delta, qreal gamma)
 
void densmatr_mixTwoQubitDepolarisingLocal (Qureg qureg, int qubit1, int qubit2, qreal delta, qreal gamma)
 
void densmatr_mixTwoQubitDepolarisingLocalPart1 (Qureg qureg, int qubit1, int qubit2, qreal delta)
 
void densmatr_mixTwoQubitDepolarisingQ1LocalQ2DistributedPart3 (Qureg qureg, int targetQubit, int qubit2, qreal delta, qreal gamma)
 
static int extractBit (const int locationOfBitFromRight, const long long int theEncodedNumber)
 
static long long int flipBit (const long long int number, const int bitInd)
 
static long long int insertTwoZeroBits (const long long int number, const int bit1, const int bit2)
 
static long long int insertZeroBit (const long long int number, const int index)
 
static int isOddParity (const long long int number, const int qb1, const int qb2)
 
static int maskContainsBit (const long long int mask, const int bitInd)
 
Complex statevec_calcExpecDiagonalOpLocal (Qureg qureg, DiagonalOp op)
 
Complex statevec_calcInnerProductLocal (Qureg bra, Qureg ket)
 
void statevec_calcProbOfAllOutcomesLocal (qreal *retProbs, Qureg qureg, int *qubits, int numQubits)
 
void statevec_collapseToKnownProbOutcomeDistributedRenorm (Qureg qureg, int measureQubit, qreal totalProbability)
 Renormalise parts of the state vector where measureQubit=0 or 1, based on the total probability of that qubit being in state 0 or 1. More...
 
void statevec_collapseToKnownProbOutcomeLocal (Qureg qureg, int measureQubit, int outcome, qreal totalProbability)
 Update the state vector to be consistent with measuring measureQubit=0 if outcome=0 and measureQubit=1 if outcome=1. More...
 
void statevec_collapseToOutcomeDistributedSetZero (Qureg qureg)
 
void statevec_compactUnitaryDistributed (Qureg qureg, Complex rot1, Complex rot2, ComplexArray stateVecUp, ComplexArray stateVecLo, ComplexArray stateVecOut)
 Rotate a single qubit in the state vector of probability amplitudes, given two complex numbers alpha and beta, and a subset of the state vector with upper and lower block values stored seperately. More...
 
void statevec_compactUnitaryLocal (Qureg qureg, int targetQubit, Complex alpha, Complex beta)
 
void statevec_controlledCompactUnitaryDistributed (Qureg qureg, int controlQubit, Complex rot1, Complex rot2, ComplexArray stateVecUp, ComplexArray stateVecLo, ComplexArray stateVecOut)
 Rotate a single qubit in the state vector of probability amplitudes, given two complex numbers alpha and beta and a subset of the state vector with upper and lower block values stored seperately. More...
 
void statevec_controlledCompactUnitaryLocal (Qureg qureg, int controlQubit, int targetQubit, Complex alpha, Complex beta)
 
void statevec_controlledNotDistributed (Qureg qureg, int controlQubit, ComplexArray stateVecIn, ComplexArray stateVecOut)
 Rotate a single qubit by {{0,1},{1,0}. More...
 
void statevec_controlledNotLocal (Qureg qureg, int controlQubit, int targetQubit)
 
void statevec_controlledPauliYDistributed (Qureg qureg, int controlQubit, ComplexArray stateVecIn, ComplexArray stateVecOut, int conjFactor)
 
void statevec_controlledPauliYLocal (Qureg qureg, int controlQubit, int targetQubit, int conjFactor)
 
void statevec_controlledUnitaryDistributed (Qureg qureg, int controlQubit, Complex rot1, Complex rot2, ComplexArray stateVecUp, ComplexArray stateVecLo, ComplexArray stateVecOut)
 Rotate a single qubit in the state vector of probability amplitudes, given two complex numbers alpha and beta and a subset of the state vector with upper and lower block values stored seperately. More...
 
void statevec_controlledUnitaryLocal (Qureg qureg, int controlQubit, int targetQubit, ComplexMatrix2 u)
 
qreal statevec_findProbabilityOfZeroDistributed (Qureg qureg)
 Measure the probability of a specified qubit being in the zero state across all amplitudes held in this chunk. More...
 
qreal statevec_findProbabilityOfZeroLocal (Qureg qureg, int measureQubit)
 Measure the total probability of a specified qubit being in the zero state across all amplitudes in this chunk. More...
 
void statevec_hadamardDistributed (Qureg qureg, ComplexArray stateVecUp, ComplexArray stateVecLo, ComplexArray stateVecOut, int updateUpper)
 Rotate a single qubit by {{1,1},{1,-1}}/sqrt2. More...
 
void statevec_hadamardLocal (Qureg qureg, int targetQubit)
 
void statevec_multiControlledMultiQubitNotDistributed (Qureg qureg, int ctrlMask, int targMask, ComplexArray stateVecIn, ComplexArray stateVecOut)
 
void statevec_multiControlledMultiQubitNotLocal (Qureg qureg, int ctrlMask, int targMask)
 
void statevec_multiControlledMultiQubitUnitaryLocal (Qureg qureg, long long int ctrlMask, int *targs, int numTargs, ComplexMatrixN u)
 
void statevec_multiControlledTwoQubitUnitaryLocal (Qureg qureg, long long int ctrlMask, int q1, int q2, ComplexMatrix4 u)
 
void statevec_multiControlledUnitaryDistributed (Qureg qureg, int targetQubit, long long int ctrlQubitsMask, long long int ctrlFlipMask, Complex rot1, Complex rot2, ComplexArray stateVecUp, ComplexArray stateVecLo, ComplexArray stateVecOut)
 Apply a unitary operation to a single qubit in the state vector of probability amplitudes, given a subset of the state vector with upper and lower block values stored seperately. More...
 
void statevec_multiControlledUnitaryLocal (Qureg qureg, int targetQubit, long long int ctrlQubitsMask, long long int ctrlFlipMask, ComplexMatrix2 u)
 
void statevec_pauliXDistributed (Qureg qureg, ComplexArray stateVecIn, ComplexArray stateVecOut)
 Rotate a single qubit by {{0,1},{1,0}. More...
 
void statevec_pauliXLocal (Qureg qureg, int targetQubit)
 
void statevec_pauliYDistributed (Qureg qureg, ComplexArray stateVecIn, ComplexArray stateVecOut, int updateUpper, int conjFac)
 Rotate a single qubit by +-{{0,-i},{i,0}. More...
 
void statevec_pauliYLocal (Qureg qureg, int targetQubit, int conjFac)
 
void statevec_swapQubitAmpsDistributed (Qureg qureg, int pairRank, int qb1, int qb2)
 qureg.pairStateVec contains the entire set of amplitudes of the paired node which includes the set of all amplitudes which need to be swapped between |..0..1..> and |..1..0..> More...
 
void statevec_swapQubitAmpsLocal (Qureg qureg, int qb1, int qb2)
 It is ensured that all amplitudes needing to be swapped are on this node. More...
 
void statevec_unitaryDistributed (Qureg qureg, Complex rot1, Complex rot2, ComplexArray stateVecUp, ComplexArray stateVecLo, ComplexArray stateVecOut)
 Apply a unitary operation to a single qubit given a subset of the state vector with upper and lower block values stored seperately. More...
 
void statevec_unitaryLocal (Qureg qureg, int targetQubit, ComplexMatrix2 u)
 

Detailed Description

Internal functions used to implement the pure backend in ../QuEST_ops_pure.h. Do not call these functions directly. In general, qubits_cpu_local.c and qubits_cpu_mpi.c will implement the pure backend by choosing the correct function or combination of functions to use from those included here, which are defined in QuEST_cpu.c

Author
Ania Brown
Tyson Jones
Balint Koczor

Definition in file QuEST_cpu_internal.h.

Function Documentation

◆ densmatr_applyDiagonalOpLocal()

void densmatr_applyDiagonalOpLocal ( Qureg  qureg,
DiagonalOp  op 
)

Definition at line 4082 of file QuEST_cpu.c.

4082  {
4083 
4084  /* ALL values of op are pre-loaded into qureg.pairStateVector (on every node).
4085  * Furthermore, since it's gauranteed each node contains an integer number of
4086  * columns of qureg (because op upperlimits the number of nodes; 1 per element),
4087  * then we know iteration below begins at the 'top' of a column, and there is
4088  * no offset for op (pairStateVector)
4089  */
4090 
4091  long long int numAmps = qureg.numAmpsPerChunk;
4092  int opDim = (1 << op.numQubits);
4093 
4094  qreal* stateRe = qureg.stateVec.real;
4095  qreal* stateIm = qureg.stateVec.imag;
4096  qreal* opRe = qureg.pairStateVec.real;
4097  qreal* opIm = qureg.pairStateVec.imag;
4098 
4099  qreal a,b,c,d;
4100  long long int index;
4101 
4102 # ifdef _OPENMP
4103 # pragma omp parallel \
4104  shared (stateRe,stateIm, opRe,opIm, numAmps,opDim) \
4105  private (index, a,b,c,d)
4106 # endif
4107  {
4108 # ifdef _OPENMP
4109 # pragma omp for schedule (static)
4110 # endif
4111  for (index=0LL; index<numAmps; index++) {
4112  a = stateRe[index];
4113  b = stateIm[index];
4114  c = opRe[index % opDim];
4115  d = opIm[index % opDim];
4116 
4117  // (a + b i)(c + d i) = (a c - b d) + i (a d + b c)
4118  stateRe[index] = a*c - b*d;
4119  stateIm[index] = a*d + b*c;
4120  }
4121  }
4122 }

References Qureg::numAmpsPerChunk, DiagonalOp::numQubits, Qureg::pairStateVec, qreal, and Qureg::stateVec.

Referenced by densmatr_applyDiagonalOp().

◆ densmatr_calcExpecDiagonalOpLocal()

Complex densmatr_calcExpecDiagonalOpLocal ( Qureg  qureg,
DiagonalOp  op 
)

Definition at line 4167 of file QuEST_cpu.c.

4167  {
4168 
4169  /* since for every 1 element in \p op, there exists a column in \p qureg,
4170  * we know that the elements in \p op live on the same node as the
4171  * corresponding diagonal elements of \p qureg. This means, the problem is
4172  * embarrassingly parallelisable, and the code below works for both
4173  * serial and distributed modes.
4174  */
4175 
4176  // computes first local index containing a diagonal element
4177  long long int diagSpacing = 1LL + (1LL << qureg.numQubitsRepresented);
4178  long long int numPrevDiags = (qureg.chunkId>0)? 1+(qureg.chunkId*qureg.numAmpsPerChunk)/diagSpacing : 0;
4179  long long int globalIndNextDiag = diagSpacing * numPrevDiags;
4180  long long int localIndNextDiag = globalIndNextDiag % qureg.numAmpsPerChunk;
4181  long long int numAmps = qureg.numAmpsPerChunk;
4182 
4183  qreal* stateReal = qureg.stateVec.real;
4184  qreal* stateImag = qureg.stateVec.imag;
4185  qreal* opReal = op.real;
4186  qreal* opImag = op.imag;
4187 
4188  qreal expecRe = 0;
4189  qreal expecIm = 0;
4190 
4191  long long int stateInd;
4192  long long int opInd;
4193  qreal matRe, matIm, opRe, opIm;
4194 
4195  // visits every diagonal element with global index (2^n + 1)i for i in [0, 2^n-1]
4196 
4197 # ifdef _OPENMP
4198 # pragma omp parallel \
4199  shared (stateReal,stateImag, opReal,opImag, localIndNextDiag,diagSpacing,numAmps) \
4200  private (stateInd,opInd, matRe,matIm, opRe,opIm) \
4201  reduction ( +:expecRe, expecIm )
4202 # endif
4203  {
4204 # ifdef _OPENMP
4205 # pragma omp for schedule (static)
4206 # endif
4207  for (stateInd=localIndNextDiag; stateInd < numAmps; stateInd += diagSpacing) {
4208 
4209  matRe = stateReal[stateInd];
4210  matIm = stateImag[stateInd];
4211  opInd = (stateInd - localIndNextDiag) / diagSpacing;
4212  opRe = opReal[opInd];
4213  opIm = opImag[opInd];
4214 
4215  // (matRe + matIm i)(opRe + opIm i) =
4216  // (matRe opRe - matIm opIm) + i (matRe opIm + matIm opRe)
4217  expecRe += matRe * opRe - matIm * opIm;
4218  expecIm += matRe * opIm + matIm * opRe;
4219  }
4220  }
4221 
4222  Complex expecVal;
4223  expecVal.real = expecRe;
4224  expecVal.imag = expecIm;
4225  return expecVal;
4226 }

References Qureg::chunkId, Complex::imag, DiagonalOp::imag, Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, qreal, Complex::real, DiagonalOp::real, and Qureg::stateVec.

Referenced by densmatr_calcExpecDiagonalOp().

◆ densmatr_calcFidelityLocal()

qreal densmatr_calcFidelityLocal ( Qureg  qureg,
Qureg  pureState 
)

computes a few dens-columns-worth of (vec^*T) dens * vec

Definition at line 1001 of file QuEST_cpu.c.

1001  {
1002 
1003  /* Here, elements of pureState are not accessed (instead grabbed from qureg.pair).
1004  * We only consult the attributes.
1005  *
1006  * qureg is a density matrix, and pureState is a statevector.
1007  * Every node contains as many columns of qureg as amps by pureState.
1008  * (each node contains an integer, exponent-of-2 number of whole columns of qureg)
1009  * Ergo, this node contains columns:
1010  * qureg.chunkID * pureState.numAmpsPerChunk to
1011  * (qureg.chunkID + 1) * pureState.numAmpsPerChunk
1012  *
1013  * The first pureState.numAmpsTotal elements of qureg.pairStateVec are the
1014  * entire pureState state-vector
1015  */
1016 
1017  // unpack everything for OPENMP
1018  qreal* vecRe = qureg.pairStateVec.real;
1019  qreal* vecIm = qureg.pairStateVec.imag;
1020  qreal* densRe = qureg.stateVec.real;
1021  qreal* densIm = qureg.stateVec.imag;
1022 
1023  int row, col;
1024  int dim = (int) pureState.numAmpsTotal;
1025  int colsPerNode = (int) pureState.numAmpsPerChunk;
1026  // using only int, because density matrix has squared as many amps so its
1027  // iteration would be impossible if the pureStates numAmpsTotal didn't fit into int
1028 
1029  // starting GLOBAL column index of the qureg columns on this node
1030  int startCol = (int) (qureg.chunkId * pureState.numAmpsPerChunk);
1031 
1032  qreal densElemRe, densElemIm;
1033  qreal prefacRe, prefacIm;
1034  qreal rowSumRe, rowSumIm;
1035  qreal vecElemRe, vecElemIm;
1036 
1037  // quantity computed by this node
1038  qreal globalSumRe = 0; // imag-component is assumed zero
1039 
1040 # ifdef _OPENMP
1041 # pragma omp parallel \
1042  shared (vecRe,vecIm,densRe,densIm, dim,colsPerNode,startCol) \
1043  private (row,col, prefacRe,prefacIm, rowSumRe,rowSumIm, densElemRe,densElemIm, vecElemRe,vecElemIm) \
1044  reduction ( +:globalSumRe )
1045 # endif
1046  {
1047 # ifdef _OPENMP
1048 # pragma omp for schedule (static)
1049 # endif
1050  // indices of my GLOBAL row
1051  for (row=0; row < dim; row++) {
1052 
1053  // single element of conj(pureState)
1054  prefacRe = vecRe[row];
1055  prefacIm = - vecIm[row];
1056 
1057  rowSumRe = 0;
1058  rowSumIm = 0;
1059 
1060  // indices of my LOCAL column
1061  for (col=0; col < colsPerNode; col++) {
1062 
1063  // my local density element
1064  densElemRe = densRe[row + dim*col];
1065  densElemIm = densIm[row + dim*col];
1066 
1067  // state-vector element
1068  vecElemRe = vecRe[startCol + col];
1069  vecElemIm = vecIm[startCol + col];
1070 
1071  rowSumRe += densElemRe*vecElemRe - densElemIm*vecElemIm;
1072  rowSumIm += densElemRe*vecElemIm + densElemIm*vecElemRe;
1073  }
1074 
1075  globalSumRe += rowSumRe*prefacRe - rowSumIm*prefacIm;
1076  }
1077  }
1078 
1079  return globalSumRe;
1080 }

References Qureg::chunkId, Qureg::numAmpsPerChunk, Qureg::numAmpsTotal, Qureg::pairStateVec, qreal, and Qureg::stateVec.

Referenced by densmatr_calcFidelity().

◆ densmatr_calcHilbertSchmidtDistanceSquaredLocal()

qreal densmatr_calcHilbertSchmidtDistanceSquaredLocal ( Qureg  a,
Qureg  b 
)

computes Tr((a-b) conjTrans(a-b)) = sum of abs values of (a-b)

Definition at line 934 of file QuEST_cpu.c.

934  {
935 
936  long long int index;
937  long long int numAmps = a.numAmpsPerChunk;
938 
939  qreal *aRe = a.stateVec.real;
940  qreal *aIm = a.stateVec.imag;
941  qreal *bRe = b.stateVec.real;
942  qreal *bIm = b.stateVec.imag;
943 
944  qreal trace = 0;
945  qreal difRe, difIm;
946 
947 # ifdef _OPENMP
948 # pragma omp parallel \
949  shared (aRe,aIm, bRe,bIm, numAmps) \
950  private (index,difRe,difIm) \
951  reduction ( +:trace )
952 # endif
953  {
954 # ifdef _OPENMP
955 # pragma omp for schedule (static)
956 # endif
957  for (index=0LL; index<numAmps; index++) {
958 
959  difRe = aRe[index] - bRe[index];
960  difIm = aIm[index] - bIm[index];
961  trace += difRe*difRe + difIm*difIm;
962  }
963  }
964 
965  return trace;
966 }

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

Referenced by densmatr_calcHilbertSchmidtDistance().

◆ densmatr_calcInnerProductLocal()

qreal densmatr_calcInnerProductLocal ( Qureg  a,
Qureg  b 
)

computes Tr(conjTrans(a) b) = sum of (a_ij^* b_ij)

Definition at line 969 of file QuEST_cpu.c.

969  {
970 
971  long long int index;
972  long long int numAmps = a.numAmpsPerChunk;
973 
974  qreal *aRe = a.stateVec.real;
975  qreal *aIm = a.stateVec.imag;
976  qreal *bRe = b.stateVec.real;
977  qreal *bIm = b.stateVec.imag;
978 
979  qreal trace = 0;
980 
981 # ifdef _OPENMP
982 # pragma omp parallel \
983  shared (aRe,aIm, bRe,bIm, numAmps) \
984  private (index) \
985  reduction ( +:trace )
986 # endif
987  {
988 # ifdef _OPENMP
989 # pragma omp for schedule (static)
990 # endif
991  for (index=0LL; index<numAmps; index++) {
992  trace += aRe[index]*bRe[index] + aIm[index]*bIm[index];
993  }
994  }
995 
996  return trace;
997 }

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

Referenced by densmatr_calcInnerProduct().

◆ densmatr_calcProbOfAllOutcomesLocal()

void densmatr_calcProbOfAllOutcomesLocal ( qreal *  retProbs,
Qureg  qureg,
int *  qubits,
int  numQubits 
)

Definition at line 3616 of file QuEST_cpu.c.

3616  {
3617 
3618  // clear outcomeProbs (in parallel, in case it's large)
3619  long long int numOutcomeProbs = (1 << numQubits);
3620  long long int j;
3621 
3622 # ifdef _OPENMP
3623 # pragma omp parallel \
3624  default (none) \
3625  shared (numOutcomeProbs,outcomeProbs) \
3626  private (j)
3627 # endif
3628  {
3629 # ifdef _OPENMP
3630 # pragma omp for schedule (static)
3631 # endif
3632  for (j=0; j<numOutcomeProbs; j++)
3633  outcomeProbs[j] = 0;
3634  }
3635 
3636  // compute first local index containing a diagonal element
3637  long long int localNumAmps = qureg.numAmpsPerChunk;
3638  long long int densityDim = (1LL << qureg.numQubitsRepresented);
3639  long long int diagSpacing = 1LL + densityDim;
3640  long long int maxNumDiagsPerChunk = 1 + localNumAmps / diagSpacing;
3641  long long int numPrevDiags = (qureg.chunkId>0)? 1+(qureg.chunkId*localNumAmps)/diagSpacing : 0;
3642  long long int globalIndNextDiag = diagSpacing * numPrevDiags;
3643  long long int localIndNextDiag = globalIndNextDiag % localNumAmps;
3644 
3645  // computes how many diagonals are contained in this chunk
3646  long long int numDiagsInThisChunk = maxNumDiagsPerChunk;
3647  if (localIndNextDiag + (numDiagsInThisChunk-1)*diagSpacing >= localNumAmps)
3648  numDiagsInThisChunk -= 1;
3649 
3650  long long int visitedDiags; // number of visited diagonals in this chunk so far
3651  long long int basisStateInd; // current diagonal index being considered
3652  long long int index; // index in the local chunk
3653 
3654  int q;
3655  long long int outcomeInd;
3656  qreal *stateRe = qureg.stateVec.real;
3657 
3658 # ifdef _OPENMP
3659 # pragma omp parallel \
3660  shared (localIndNextDiag, numPrevDiags, diagSpacing, stateRe, numDiagsInThisChunk, qubits,numQubits, outcomeProbs) \
3661  private (visitedDiags, basisStateInd, index, q,outcomeInd)
3662 # endif
3663  {
3664 # ifdef _OPENMP
3665 # pragma omp for schedule (static)
3666 # endif
3667  // sums the diagonal elems of the density matrix where measureQubit=0
3668  for (visitedDiags = 0; visitedDiags < numDiagsInThisChunk; visitedDiags++) {
3669 
3670  basisStateInd = numPrevDiags + visitedDiags;
3671  index = localIndNextDiag + diagSpacing * visitedDiags;
3672 
3673  // determine outcome implied by basisStateInd
3674  outcomeInd = 0;
3675  for (q=0; q<numQubits; q++)
3676  outcomeInd += extractBit(qubits[q], basisStateInd) * (1LL << q);
3677 
3678  // atomicly update corresponding outcome array element
3679  # ifdef _OPENMP
3680  # pragma omp atomic
3681  # endif
3682  outcomeProbs[outcomeInd] += stateRe[index];
3683  }
3684  }
3685 }

References Qureg::chunkId, extractBit(), Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, qreal, and Qureg::stateVec.

Referenced by densmatr_calcProbOfAllOutcomes().

◆ densmatr_calcPurityLocal()

qreal densmatr_calcPurityLocal ( Qureg  qureg)

Definition at line 872 of file QuEST_cpu.c.

872  {
873 
874  /* sum of qureg^2, which is sum_i |qureg[i]|^2 */
875  long long int index;
876  long long int numAmps = qureg.numAmpsPerChunk;
877 
878  qreal trace = 0;
879  qreal *vecRe = qureg.stateVec.real;
880  qreal *vecIm = qureg.stateVec.imag;
881 
882 # ifdef _OPENMP
883 # pragma omp parallel \
884  shared (vecRe, vecIm, numAmps) \
885  private (index) \
886  reduction ( +:trace )
887 # endif
888  {
889 # ifdef _OPENMP
890 # pragma omp for schedule (static)
891 # endif
892  for (index=0LL; index<numAmps; index++) {
893 
894  trace += vecRe[index]*vecRe[index] + vecIm[index]*vecIm[index];
895  }
896  }
897 
898  return trace;
899 }

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

Referenced by densmatr_calcPurity().

◆ densmatr_findProbabilityOfZeroLocal()

qreal densmatr_findProbabilityOfZeroLocal ( Qureg  qureg,
int  measureQubit 
)

Definition at line 3402 of file QuEST_cpu.c.

3402  {
3403 
3404  // computes first local index containing a diagonal element
3405  long long int localNumAmps = qureg.numAmpsPerChunk;
3406  long long int densityDim = (1LL << qureg.numQubitsRepresented);
3407  long long int diagSpacing = 1LL + densityDim;
3408  long long int maxNumDiagsPerChunk = 1 + localNumAmps / diagSpacing;
3409  long long int numPrevDiags = (qureg.chunkId>0)? 1+(qureg.chunkId*localNumAmps)/diagSpacing : 0;
3410  long long int globalIndNextDiag = diagSpacing * numPrevDiags;
3411  long long int localIndNextDiag = globalIndNextDiag % localNumAmps;
3412 
3413  // computes how many diagonals are contained in this chunk
3414  long long int numDiagsInThisChunk = maxNumDiagsPerChunk;
3415  if (localIndNextDiag + (numDiagsInThisChunk-1)*diagSpacing >= localNumAmps)
3416  numDiagsInThisChunk -= 1;
3417 
3418  long long int visitedDiags; // number of visited diagonals in this chunk so far
3419  long long int basisStateInd; // current diagonal index being considered
3420  long long int index; // index in the local chunk
3421 
3422  qreal zeroProb = 0;
3423  qreal *stateVecReal = qureg.stateVec.real;
3424 
3425 # ifdef _OPENMP
3426 # pragma omp parallel \
3427  shared (localIndNextDiag, numPrevDiags, diagSpacing, stateVecReal, numDiagsInThisChunk) \
3428  private (visitedDiags, basisStateInd, index) \
3429  reduction ( +:zeroProb )
3430 # endif
3431  {
3432 # ifdef _OPENMP
3433 # pragma omp for schedule (static)
3434 # endif
3435  // sums the diagonal elems of the density matrix where measureQubit=0
3436  for (visitedDiags = 0; visitedDiags < numDiagsInThisChunk; visitedDiags++) {
3437 
3438  basisStateInd = numPrevDiags + visitedDiags;
3439  index = localIndNextDiag + diagSpacing * visitedDiags;
3440 
3441  if (extractBit(measureQubit, basisStateInd) == 0)
3442  zeroProb += stateVecReal[index]; // assume imag[diagonls] ~ 0
3443 
3444  }
3445  }
3446 
3447  return zeroProb;
3448 }

References Qureg::chunkId, extractBit(), Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, qreal, and Qureg::stateVec.

Referenced by densmatr_calcProbOfOutcome().

◆ densmatr_initPureStateLocal()

void densmatr_initPureStateLocal ( Qureg  targetQureg,
Qureg  copyQureg 
)

Definition at line 1195 of file QuEST_cpu.c.

1195  {
1196 
1197  /* copyQureg amps aren't explicitly used - they're accessed through targetQureg.pair,
1198  * which contains the full pure statevector.
1199  * targetQureg has as many columns on node as copyQureg has amps
1200  */
1201 
1202  long long int colOffset = targetQureg.chunkId * copyQureg.numAmpsPerChunk;
1203  long long int colsPerNode = copyQureg.numAmpsPerChunk;
1204  long long int rowsPerNode = copyQureg.numAmpsTotal;
1205 
1206  // unpack vars for OpenMP
1207  qreal* vecRe = targetQureg.pairStateVec.real;
1208  qreal* vecIm = targetQureg.pairStateVec.imag;
1209  qreal* densRe = targetQureg.stateVec.real;
1210  qreal* densIm = targetQureg.stateVec.imag;
1211 
1212  long long int col, row, index;
1213 
1214  // a_i conj(a_j) |i><j|
1215  qreal ketRe, ketIm, braRe, braIm;
1216 
1217 # ifdef _OPENMP
1218 # pragma omp parallel \
1219  default (none) \
1220  shared (colOffset, colsPerNode,rowsPerNode, vecRe,vecIm,densRe,densIm) \
1221  private (col,row, ketRe,ketIm,braRe,braIm, index)
1222 # endif
1223  {
1224 # ifdef _OPENMP
1225 # pragma omp for schedule (static)
1226 # endif
1227  // local column
1228  for (col=0; col < colsPerNode; col++) {
1229 
1230  // global row
1231  for (row=0; row < rowsPerNode; row++) {
1232 
1233  // get pure state amps
1234  ketRe = vecRe[row];
1235  ketIm = vecIm[row];
1236  braRe = vecRe[col + colOffset];
1237  braIm = - vecIm[col + colOffset]; // minus for conjugation
1238 
1239  // update density matrix
1240  index = row + col*rowsPerNode; // local ind
1241  densRe[index] = ketRe*braRe - ketIm*braIm;
1242  densIm[index] = ketRe*braIm + ketIm*braRe;
1243  }
1244  }
1245  }
1246 }

References Qureg::chunkId, Qureg::numAmpsPerChunk, Qureg::numAmpsTotal, Qureg::pairStateVec, qreal, and Qureg::stateVec.

Referenced by densmatr_initPureState().

◆ densmatr_mixDampingDistributed()

void densmatr_mixDampingDistributed ( Qureg  qureg,
int  targetQubit,
qreal  damping 
)

Definition at line 306 of file QuEST_cpu.c.

306  {
307  qreal retain=1-damping;
308  qreal dephase=sqrt(1-damping);
309 
310  // multiply the off-diagonal (|0><1| and |1><0|) terms by sqrt(1-damping)
311  densmatr_oneQubitDegradeOffDiagonal(qureg, targetQubit, dephase);
312 
313  // below, we modify the diagonals terms which require |1><1| to |0><0| communication
314 
315  long long int sizeInnerBlock, sizeInnerHalfBlock;
316  long long int sizeOuterColumn, sizeOuterHalfColumn;
317  long long int thisInnerBlock, // current block
318  thisOuterColumn, // current column in density matrix
319  thisIndex, // current index in (density matrix representation) state vector
320  thisIndexInOuterColumn,
321  thisIndexInInnerBlock;
322  int outerBit;
323  int stateBit;
324 
325  long long int thisTask;
326  long long int numTasks=qureg.numAmpsPerChunk>>1;
327 
328  // set dimensions
329  sizeInnerHalfBlock = 1LL << targetQubit;
330  sizeInnerBlock = 2LL * sizeInnerHalfBlock;
331  sizeOuterColumn = 1LL << qureg.numQubitsRepresented;
332  sizeOuterHalfColumn = sizeOuterColumn >> 1;
333 
334 # ifdef _OPENMP
335 # pragma omp parallel \
336  default (none) \
337  shared (sizeInnerBlock,sizeInnerHalfBlock,sizeOuterColumn,sizeOuterHalfColumn, \
338  qureg,damping, retain, dephase, numTasks,targetQubit) \
339  private (thisTask,thisInnerBlock,thisOuterColumn,thisIndex,thisIndexInOuterColumn, \
340  thisIndexInInnerBlock,outerBit, stateBit)
341 # endif
342  {
343 # ifdef _OPENMP
344 # pragma omp for schedule (static)
345 # endif
346  // thisTask iterates over half the elements in this process' chunk of the density matrix
347  // treat this as iterating over all columns, then iterating over half the values
348  // within one column.
349  // If this function has been called, this process' chunk contains half an
350  // outer block or less
351  for (thisTask=0; thisTask<numTasks; thisTask++) {
352  // we want to process all columns in the density matrix,
353  // updating the values for half of each column (one half of each inner block)
354  thisOuterColumn = thisTask / sizeOuterHalfColumn;
355  thisIndexInOuterColumn = thisTask&(sizeOuterHalfColumn-1); // thisTask % sizeOuterHalfColumn
356  thisInnerBlock = thisIndexInOuterColumn/sizeInnerHalfBlock;
357  // get index in state vector corresponding to upper inner block
358  thisIndexInInnerBlock = thisTask&(sizeInnerHalfBlock-1); // thisTask % sizeInnerHalfBlock
359  thisIndex = thisOuterColumn*sizeOuterColumn + thisInnerBlock*sizeInnerBlock
360  + thisIndexInInnerBlock;
361  // check if we are in the upper or lower half of an outer block
362  outerBit = extractBit(targetQubit, (thisIndex+qureg.numAmpsPerChunk*qureg.chunkId)>>qureg.numQubitsRepresented);
363  // if we are in the lower half of an outer block, shift to be in the lower half
364  // of the inner block as well (we want to dephase |0><0| and |1><1| only)
365  thisIndex += outerBit*(sizeInnerHalfBlock);
366 
367  // NOTE: at this point thisIndex should be the index of the element we want to
368  // dephase in the chunk of the state vector on this process, in the
369  // density matrix representation.
370  // thisTask is the index of the pair element in pairStateVec
371 
372  // Extract state bit, is 0 if thisIndex corresponds to a state with 0 in the target qubit
373  // and is 1 if thisIndex corresponds to a state with 1 in the target qubit
374  stateBit = extractBit(targetQubit, (thisIndex+qureg.numAmpsPerChunk*qureg.chunkId));
375 
376  // state[thisIndex] = (1-depolLevel)*state[thisIndex] + depolLevel*(state[thisIndex]
377  // + pair[thisTask])/2
378  if(stateBit == 0){
379  qureg.stateVec.real[thisIndex] = qureg.stateVec.real[thisIndex] +
380  damping*( qureg.pairStateVec.real[thisTask]);
381 
382  qureg.stateVec.imag[thisIndex] = qureg.stateVec.imag[thisIndex] +
383  damping*( qureg.pairStateVec.imag[thisTask]);
384  } else{
385  qureg.stateVec.real[thisIndex] = retain*qureg.stateVec.real[thisIndex];
386 
387  qureg.stateVec.imag[thisIndex] = retain*qureg.stateVec.imag[thisIndex];
388  }
389  }
390  }
391 }

References Qureg::chunkId, densmatr_oneQubitDegradeOffDiagonal(), extractBit(), Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, Qureg::pairStateVec, qreal, and Qureg::stateVec.

Referenced by densmatr_mixDamping().

◆ densmatr_mixDampingLocal()

void densmatr_mixDampingLocal ( Qureg  qureg,
int  targetQubit,
qreal  damping 
)

Definition at line 180 of file QuEST_cpu.c.

180  {
181  qreal retain=1-damping;
182  qreal dephase=sqrt(retain);
183 
184  long long int numTasks = qureg.numAmpsPerChunk;
185  long long int innerMask = 1LL << targetQubit;
186  long long int outerMask = 1LL << (targetQubit + (qureg.numQubitsRepresented));
187  long long int totMask = innerMask|outerMask;
188 
189  long long int thisTask;
190  long long int partner;
191  long long int thisPattern;
192 
193  //qreal realAv, imagAv;
194 
195 # ifdef _OPENMP
196 # pragma omp parallel \
197  default (none) \
198  shared (innerMask,outerMask,totMask,qureg,retain,damping,dephase,numTasks) \
199  private (thisTask,partner,thisPattern)
200 # endif
201  {
202 # ifdef _OPENMP
203 # pragma omp for schedule (static)
204 # endif
205  for (thisTask=0; thisTask<numTasks; thisTask++){
206  thisPattern = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMask;
207  if ((thisPattern==innerMask) || (thisPattern==outerMask)){
208  // do dephase
209  // the lines below will degrade the off-diagonal terms |..0..><..1..| and |..1..><..0..|
210  qureg.stateVec.real[thisTask] = dephase*qureg.stateVec.real[thisTask];
211  qureg.stateVec.imag[thisTask] = dephase*qureg.stateVec.imag[thisTask];
212  } else {
213  if ((thisTask&totMask)==0){ //this element relates to targetQubit in state 0
214  // do depolarise
215  partner = thisTask | totMask;
216  //realAv = (qureg.stateVec.real[thisTask] + qureg.stateVec.real[partner]) /2 ;
217  //imagAv = (qureg.stateVec.imag[thisTask] + qureg.stateVec.imag[partner]) /2 ;
218 
219  qureg.stateVec.real[thisTask] = qureg.stateVec.real[thisTask] + damping*qureg.stateVec.real[partner];
220  qureg.stateVec.imag[thisTask] = qureg.stateVec.imag[thisTask] + damping*qureg.stateVec.imag[partner];
221 
222  qureg.stateVec.real[partner] = retain*qureg.stateVec.real[partner];
223  qureg.stateVec.imag[partner] = retain*qureg.stateVec.imag[partner];
224  }
225  }
226  }
227  }
228 }

References Qureg::chunkId, Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, qreal, and Qureg::stateVec.

Referenced by densmatr_mixDamping().

◆ densmatr_mixDepolarisingDistributed()

void densmatr_mixDepolarisingDistributed ( Qureg  qureg,
int  targetQubit,
qreal  depolLevel 
)

Definition at line 230 of file QuEST_cpu.c.

230  {
231 
232  // first do dephase part.
233  // TODO -- this might be more efficient to do at the same time as the depolarise if we move to
234  // iterating over all elements in the state vector for the purpose of vectorisation
235  // TODO -- if we keep this split, move this function to densmatr_mixDepolarising()
236  densmatr_mixDephasing(qureg, targetQubit, depolLevel);
237 
238  long long int sizeInnerBlock, sizeInnerHalfBlock;
239  long long int sizeOuterColumn, sizeOuterHalfColumn;
240  long long int thisInnerBlock, // current block
241  thisOuterColumn, // current column in density matrix
242  thisIndex, // current index in (density matrix representation) state vector
243  thisIndexInOuterColumn,
244  thisIndexInInnerBlock;
245  int outerBit;
246 
247  long long int thisTask;
248  long long int numTasks=qureg.numAmpsPerChunk>>1;
249 
250  // set dimensions
251  sizeInnerHalfBlock = 1LL << targetQubit;
252  sizeInnerBlock = 2LL * sizeInnerHalfBlock;
253  sizeOuterColumn = 1LL << qureg.numQubitsRepresented;
254  sizeOuterHalfColumn = sizeOuterColumn >> 1;
255 
256 # ifdef _OPENMP
257 # pragma omp parallel \
258  default (none) \
259  shared (sizeInnerBlock,sizeInnerHalfBlock,sizeOuterColumn,sizeOuterHalfColumn, \
260  qureg,depolLevel,numTasks,targetQubit) \
261  private (thisTask,thisInnerBlock,thisOuterColumn,thisIndex,thisIndexInOuterColumn, \
262  thisIndexInInnerBlock,outerBit)
263 # endif
264  {
265 # ifdef _OPENMP
266 # pragma omp for schedule (static)
267 # endif
268  // thisTask iterates over half the elements in this process' chunk of the density matrix
269  // treat this as iterating over all columns, then iterating over half the values
270  // within one column.
271  // If this function has been called, this process' chunk contains half an
272  // outer block or less
273  for (thisTask=0; thisTask<numTasks; thisTask++) {
274  // we want to process all columns in the density matrix,
275  // updating the values for half of each column (one half of each inner block)
276  thisOuterColumn = thisTask / sizeOuterHalfColumn;
277  thisIndexInOuterColumn = thisTask&(sizeOuterHalfColumn-1); // thisTask % sizeOuterHalfColumn
278  thisInnerBlock = thisIndexInOuterColumn/sizeInnerHalfBlock;
279  // get index in state vector corresponding to upper inner block
280  thisIndexInInnerBlock = thisTask&(sizeInnerHalfBlock-1); // thisTask % sizeInnerHalfBlock
281  thisIndex = thisOuterColumn*sizeOuterColumn + thisInnerBlock*sizeInnerBlock
282  + thisIndexInInnerBlock;
283  // check if we are in the upper or lower half of an outer block
284  outerBit = extractBit(targetQubit, (thisIndex+qureg.numAmpsPerChunk*qureg.chunkId)>>qureg.numQubitsRepresented);
285  // if we are in the lower half of an outer block, shift to be in the lower half
286  // of the inner block as well (we want to dephase |0><0| and |1><1| only)
287  thisIndex += outerBit*(sizeInnerHalfBlock);
288 
289  // NOTE: at this point thisIndex should be the index of the element we want to
290  // dephase in the chunk of the state vector on this process, in the
291  // density matrix representation.
292  // thisTask is the index of the pair element in pairStateVec
293 
294 
295  // state[thisIndex] = (1-depolLevel)*state[thisIndex] + depolLevel*(state[thisIndex]
296  // + pair[thisTask])/2
297  qureg.stateVec.real[thisIndex] = (1-depolLevel)*qureg.stateVec.real[thisIndex] +
298  depolLevel*(qureg.stateVec.real[thisIndex] + qureg.pairStateVec.real[thisTask])/2;
299 
300  qureg.stateVec.imag[thisIndex] = (1-depolLevel)*qureg.stateVec.imag[thisIndex] +
301  depolLevel*(qureg.stateVec.imag[thisIndex] + qureg.pairStateVec.imag[thisTask])/2;
302  }
303  }
304 }

References Qureg::chunkId, densmatr_mixDephasing(), extractBit(), Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, Qureg::pairStateVec, and Qureg::stateVec.

Referenced by densmatr_mixDepolarising().

◆ densmatr_mixDepolarisingLocal()

void densmatr_mixDepolarisingLocal ( Qureg  qureg,
int  targetQubit,
qreal  depolLevel 
)

Definition at line 131 of file QuEST_cpu.c.

131  {
132  qreal retain=1-depolLevel;
133 
134  long long int numTasks = qureg.numAmpsPerChunk;
135  long long int innerMask = 1LL << targetQubit;
136  long long int outerMask = 1LL << (targetQubit + (qureg.numQubitsRepresented));
137  long long int totMask = innerMask|outerMask;
138 
139  long long int thisTask;
140  long long int partner;
141  long long int thisPattern;
142 
143  qreal realAv, imagAv;
144 
145 # ifdef _OPENMP
146 # pragma omp parallel \
147  default (none) \
148  shared (innerMask,outerMask,totMask,qureg,retain,depolLevel,numTasks) \
149  private (thisTask,partner,thisPattern,realAv,imagAv)
150 # endif
151  {
152 # ifdef _OPENMP
153 # pragma omp for schedule (static)
154 # endif
155  for (thisTask=0; thisTask<numTasks; thisTask++){
156  thisPattern = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMask;
157  if ((thisPattern==innerMask) || (thisPattern==outerMask)){
158  // do dephase
159  // the lines below will degrade the off-diagonal terms |..0..><..1..| and |..1..><..0..|
160  qureg.stateVec.real[thisTask] = retain*qureg.stateVec.real[thisTask];
161  qureg.stateVec.imag[thisTask] = retain*qureg.stateVec.imag[thisTask];
162  } else {
163  if ((thisTask&totMask)==0){ //this element relates to targetQubit in state 0
164  // do depolarise
165  partner = thisTask | totMask;
166  realAv = (qureg.stateVec.real[thisTask] + qureg.stateVec.real[partner]) /2 ;
167  imagAv = (qureg.stateVec.imag[thisTask] + qureg.stateVec.imag[partner]) /2 ;
168 
169  qureg.stateVec.real[thisTask] = retain*qureg.stateVec.real[thisTask] + depolLevel*realAv;
170  qureg.stateVec.imag[thisTask] = retain*qureg.stateVec.imag[thisTask] + depolLevel*imagAv;
171 
172  qureg.stateVec.real[partner] = retain*qureg.stateVec.real[partner] + depolLevel*realAv;
173  qureg.stateVec.imag[partner] = retain*qureg.stateVec.imag[partner] + depolLevel*imagAv;
174  }
175  }
176  }
177  }
178 }

References Qureg::chunkId, Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, qreal, and Qureg::stateVec.

Referenced by densmatr_mixDepolarising().

◆ densmatr_mixTwoQubitDepolarisingDistributed()

void densmatr_mixTwoQubitDepolarisingDistributed ( Qureg  qureg,
int  targetQubit,
int  qubit2,
qreal  delta,
qreal  gamma 
)

Definition at line 547 of file QuEST_cpu.c.

548  {
549 
550  long long int sizeInnerBlockQ1, sizeInnerHalfBlockQ1;
551  long long int sizeInnerBlockQ2, sizeInnerHalfBlockQ2, sizeInnerQuarterBlockQ2;
552  long long int sizeOuterColumn, sizeOuterQuarterColumn;
553  long long int thisInnerBlockQ2,
554  thisOuterColumn, // current column in density matrix
555  thisIndex, // current index in (density matrix representation) state vector
556  thisIndexInOuterColumn,
557  thisIndexInInnerBlockQ1,
558  thisIndexInInnerBlockQ2,
559  thisInnerBlockQ1InInnerBlockQ2;
560  int outerBitQ1, outerBitQ2;
561 
562  long long int thisTask;
563  long long int numTasks=qureg.numAmpsPerChunk>>2;
564 
565  // set dimensions
566  sizeInnerHalfBlockQ1 = 1LL << targetQubit;
567  sizeInnerHalfBlockQ2 = 1LL << qubit2;
568  sizeInnerQuarterBlockQ2 = sizeInnerHalfBlockQ2 >> 1;
569  sizeInnerBlockQ2 = sizeInnerHalfBlockQ2 << 1;
570  sizeInnerBlockQ1 = 2LL * sizeInnerHalfBlockQ1;
571  sizeOuterColumn = 1LL << qureg.numQubitsRepresented;
572  sizeOuterQuarterColumn = sizeOuterColumn >> 2;
573 
574 # ifdef _OPENMP
575 # pragma omp parallel \
576  default (none) \
577  shared (sizeInnerBlockQ1,sizeInnerHalfBlockQ1,sizeInnerBlockQ2,sizeInnerHalfBlockQ2,sizeInnerQuarterBlockQ2,\
578  sizeOuterColumn,sizeOuterQuarterColumn,qureg,delta,gamma,numTasks,targetQubit,qubit2) \
579  private (thisTask,thisInnerBlockQ2,thisInnerBlockQ1InInnerBlockQ2, \
580  thisOuterColumn,thisIndex,thisIndexInOuterColumn, \
581  thisIndexInInnerBlockQ1,thisIndexInInnerBlockQ2,outerBitQ1,outerBitQ2)
582 # endif
583  {
584 # ifdef _OPENMP
585 # pragma omp for schedule (static)
586 # endif
587  // thisTask iterates over half the elements in this process' chunk of the density matrix
588  // treat this as iterating over all columns, then iterating over half the values
589  // within one column.
590  // If this function has been called, this process' chunk contains half an
591  // outer block or less
592  for (thisTask=0; thisTask<numTasks; thisTask++) {
593  // we want to process all columns in the density matrix,
594  // updating the values for half of each column (one half of each inner block)
595  thisOuterColumn = thisTask / sizeOuterQuarterColumn;
596  // thisTask % sizeOuterQuarterColumn
597  thisIndexInOuterColumn = thisTask&(sizeOuterQuarterColumn-1);
598  thisInnerBlockQ2 = thisIndexInOuterColumn / sizeInnerQuarterBlockQ2;
599  // thisTask % sizeInnerQuarterBlockQ2;
600  thisIndexInInnerBlockQ2 = thisTask&(sizeInnerQuarterBlockQ2-1);
601  thisInnerBlockQ1InInnerBlockQ2 = thisIndexInInnerBlockQ2 / sizeInnerHalfBlockQ1;
602  // thisTask % sizeInnerHalfBlockQ1;
603  thisIndexInInnerBlockQ1 = thisTask&(sizeInnerHalfBlockQ1-1);
604 
605  // get index in state vector corresponding to upper inner block
606  thisIndex = thisOuterColumn*sizeOuterColumn + thisInnerBlockQ2*sizeInnerBlockQ2
607  + thisInnerBlockQ1InInnerBlockQ2*sizeInnerBlockQ1 + thisIndexInInnerBlockQ1;
608 
609  // check if we are in the upper or lower half of an outer block for Q1
610  outerBitQ1 = extractBit(targetQubit, (thisIndex+qureg.numAmpsPerChunk*qureg.chunkId)>>qureg.numQubitsRepresented);
611  // if we are in the lower half of an outer block, shift to be in the lower half
612  // of the inner block as well (we want to dephase |0><0| and |1><1| only)
613  thisIndex += outerBitQ1*(sizeInnerHalfBlockQ1);
614 
615  // check if we are in the upper or lower half of an outer block for Q2
616  outerBitQ2 = extractBit(qubit2, (thisIndex+qureg.numAmpsPerChunk*qureg.chunkId)>>qureg.numQubitsRepresented);
617  // if we are in the lower half of an outer block, shift to be in the lower half
618  // of the inner block as well (we want to dephase |0><0| and |1><1| only)
619  thisIndex += outerBitQ2*(sizeInnerQuarterBlockQ2<<1);
620 
621  // NOTE: at this point thisIndex should be the index of the element we want to
622  // dephase in the chunk of the state vector on this process, in the
623  // density matrix representation.
624  // thisTask is the index of the pair element in pairStateVec
625 
626 
627  // state[thisIndex] = (1-depolLevel)*state[thisIndex] + depolLevel*(state[thisIndex]
628  // + pair[thisTask])/2
629  // NOTE: must set gamma=1 if using this function for steps 1 or 2
630  qureg.stateVec.real[thisIndex] = gamma*(qureg.stateVec.real[thisIndex] +
631  delta*qureg.pairStateVec.real[thisTask]);
632  qureg.stateVec.imag[thisIndex] = gamma*(qureg.stateVec.imag[thisIndex] +
633  delta*qureg.pairStateVec.imag[thisTask]);
634  }
635  }
636 }

References Qureg::chunkId, extractBit(), Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, Qureg::pairStateVec, and Qureg::stateVec.

Referenced by densmatr_mixTwoQubitDepolarising().

◆ densmatr_mixTwoQubitDepolarisingLocal()

void densmatr_mixTwoQubitDepolarisingLocal ( Qureg  qureg,
int  qubit1,
int  qubit2,
qreal  delta,
qreal  gamma 
)

Definition at line 393 of file QuEST_cpu.c.

393  {
394  long long int numTasks = qureg.numAmpsPerChunk;
395  long long int innerMaskQubit1 = 1LL << qubit1;
396  long long int outerMaskQubit1= 1LL << (qubit1 + qureg.numQubitsRepresented);
397  long long int totMaskQubit1 = innerMaskQubit1 | outerMaskQubit1;
398  long long int innerMaskQubit2 = 1LL << qubit2;
399  long long int outerMaskQubit2 = 1LL << (qubit2 + qureg.numQubitsRepresented);
400  long long int totMaskQubit2 = innerMaskQubit2 | outerMaskQubit2;
401 
402  long long int thisTask;
403  long long int partner;
404  long long int thisPatternQubit1, thisPatternQubit2;
405 
406  qreal real00, imag00;
407 
408 # ifdef _OPENMP
409 # pragma omp parallel \
410  default (none) \
411  shared (totMaskQubit1,totMaskQubit2,qureg,delta,gamma,numTasks) \
412  private (thisTask,partner,thisPatternQubit1,thisPatternQubit2,real00,imag00)
413 # endif
414  {
415 # ifdef _OPENMP
416 # pragma omp for schedule (static)
417 # endif
418  //--------------------------------------- STEP ONE ---------------------
419  for (thisTask=0; thisTask<numTasks; thisTask++){
420  thisPatternQubit1 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit1;
421  thisPatternQubit2 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit2;
422  if ((thisPatternQubit1==0) && ((thisPatternQubit2==0)
423  || (thisPatternQubit2==totMaskQubit2))){
424  //this element of form |...X...0...><...X...0...| for X either 0 or 1.
425  partner = thisTask | totMaskQubit1;
426  real00 = qureg.stateVec.real[thisTask];
427  imag00 = qureg.stateVec.imag[thisTask];
428 
429  qureg.stateVec.real[thisTask] = qureg.stateVec.real[thisTask]
430  + delta*qureg.stateVec.real[partner];
431  qureg.stateVec.imag[thisTask] = qureg.stateVec.imag[thisTask]
432  + delta*qureg.stateVec.imag[partner];
433 
434  qureg.stateVec.real[partner] = qureg.stateVec.real[partner] + delta*real00;
435  qureg.stateVec.imag[partner] = qureg.stateVec.imag[partner] + delta*imag00;
436 
437  }
438  }
439 # ifdef _OPENMP
440 # pragma omp for schedule (static)
441 # endif
442  //--------------------------------------- STEP TWO ---------------------
443  for (thisTask=0; thisTask<numTasks; thisTask++){
444  thisPatternQubit1 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit1;
445  thisPatternQubit2 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit2;
446  if ((thisPatternQubit2==0) && ((thisPatternQubit1==0)
447  || (thisPatternQubit1==totMaskQubit1))){
448  //this element of form |...0...X...><...0...X...| for X either 0 or 1.
449  partner = thisTask | totMaskQubit2;
450  real00 = qureg.stateVec.real[thisTask];
451  imag00 = qureg.stateVec.imag[thisTask];
452 
453  qureg.stateVec.real[thisTask] = qureg.stateVec.real[thisTask]
454  + delta*qureg.stateVec.real[partner];
455  qureg.stateVec.imag[thisTask] = qureg.stateVec.imag[thisTask]
456  + delta*qureg.stateVec.imag[partner];
457 
458  qureg.stateVec.real[partner] = qureg.stateVec.real[partner] + delta*real00;
459  qureg.stateVec.imag[partner] = qureg.stateVec.imag[partner] + delta*imag00;
460 
461  }
462  }
463 
464 # ifdef _OPENMP
465 # pragma omp for schedule (static)
466 # endif
467  //--------------------------------------- STEP THREE ---------------------
468  for (thisTask=0; thisTask<numTasks; thisTask++){
469  thisPatternQubit1 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit1;
470  thisPatternQubit2 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit2;
471  if ((thisPatternQubit2==0) && ((thisPatternQubit1==0)
472  || (thisPatternQubit1==totMaskQubit1))){
473  //this element of form |...0...X...><...0...X...| for X either 0 or 1.
474  partner = thisTask | totMaskQubit2;
475  partner = partner ^ totMaskQubit1;
476  real00 = qureg.stateVec.real[thisTask];
477  imag00 = qureg.stateVec.imag[thisTask];
478 
479  qureg.stateVec.real[thisTask] = gamma * (qureg.stateVec.real[thisTask]
480  + delta*qureg.stateVec.real[partner]);
481  qureg.stateVec.imag[thisTask] = gamma * (qureg.stateVec.imag[thisTask]
482  + delta*qureg.stateVec.imag[partner]);
483 
484  qureg.stateVec.real[partner] = gamma * (qureg.stateVec.real[partner]
485  + delta*real00);
486  qureg.stateVec.imag[partner] = gamma * (qureg.stateVec.imag[partner]
487  + delta*imag00);
488 
489  }
490  }
491  }
492 }

References Qureg::chunkId, Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, qreal, and Qureg::stateVec.

Referenced by densmatr_mixTwoQubitDepolarising().

◆ densmatr_mixTwoQubitDepolarisingLocalPart1()

void densmatr_mixTwoQubitDepolarisingLocalPart1 ( Qureg  qureg,
int  qubit1,
int  qubit2,
qreal  delta 
)

Definition at line 494 of file QuEST_cpu.c.

494  {
495  long long int numTasks = qureg.numAmpsPerChunk;
496  long long int innerMaskQubit1 = 1LL << qubit1;
497  long long int outerMaskQubit1= 1LL << (qubit1 + qureg.numQubitsRepresented);
498  long long int totMaskQubit1 = innerMaskQubit1 | outerMaskQubit1;
499  long long int innerMaskQubit2 = 1LL << qubit2;
500  long long int outerMaskQubit2 = 1LL << (qubit2 + qureg.numQubitsRepresented);
501  long long int totMaskQubit2 = innerMaskQubit2 | outerMaskQubit2;
502  // correct for being in a particular chunk
503  //totMaskQubit2 = totMaskQubit2&(qureg.numAmpsPerChunk-1); // totMaskQubit2 % numAmpsPerChunk
504 
505 
506  long long int thisTask;
507  long long int partner;
508  long long int thisPatternQubit1, thisPatternQubit2;
509 
510  qreal real00, imag00;
511 
512 # ifdef _OPENMP
513 # pragma omp parallel \
514  default (none) \
515  shared (totMaskQubit1,totMaskQubit2,qureg,delta,numTasks) \
516  private (thisTask,partner,thisPatternQubit1,thisPatternQubit2,real00,imag00)
517 # endif
518  {
519 
520 # ifdef _OPENMP
521 # pragma omp for schedule (static)
522 # endif
523  //--------------------------------------- STEP ONE ---------------------
524  for (thisTask=0; thisTask<numTasks; thisTask ++){
525  thisPatternQubit1 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit1;
526  thisPatternQubit2 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit2;
527  if ((thisPatternQubit1==0) && ((thisPatternQubit2==0)
528  || (thisPatternQubit2==totMaskQubit2))){
529  //this element of form |...X...0...><...X...0...| for X either 0 or 1.
530  partner = thisTask | totMaskQubit1;
531  real00 = qureg.stateVec.real[thisTask];
532  imag00 = qureg.stateVec.imag[thisTask];
533 
534  qureg.stateVec.real[thisTask] = qureg.stateVec.real[thisTask]
535  + delta*qureg.stateVec.real[partner];
536  qureg.stateVec.imag[thisTask] = qureg.stateVec.imag[thisTask]
537  + delta*qureg.stateVec.imag[partner];
538 
539  qureg.stateVec.real[partner] = qureg.stateVec.real[partner] + delta*real00;
540  qureg.stateVec.imag[partner] = qureg.stateVec.imag[partner] + delta*imag00;
541 
542  }
543  }
544  }
545 }

References Qureg::chunkId, Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, qreal, and Qureg::stateVec.

Referenced by densmatr_mixTwoQubitDepolarising().

◆ densmatr_mixTwoQubitDepolarisingQ1LocalQ2DistributedPart3()

void densmatr_mixTwoQubitDepolarisingQ1LocalQ2DistributedPart3 ( Qureg  qureg,
int  targetQubit,
int  qubit2,
qreal  delta,
qreal  gamma 
)

Definition at line 638 of file QuEST_cpu.c.

639  {
640 
641  long long int sizeInnerBlockQ1, sizeInnerHalfBlockQ1;
642  long long int sizeInnerBlockQ2, sizeInnerHalfBlockQ2, sizeInnerQuarterBlockQ2;
643  long long int sizeOuterColumn, sizeOuterQuarterColumn;
644  long long int thisInnerBlockQ2,
645  thisOuterColumn, // current column in density matrix
646  thisIndex, // current index in (density matrix representation) state vector
647  thisIndexInPairVector,
648  thisIndexInOuterColumn,
649  thisIndexInInnerBlockQ1,
650  thisIndexInInnerBlockQ2,
651  thisInnerBlockQ1InInnerBlockQ2;
652  int outerBitQ1, outerBitQ2;
653 
654  long long int thisTask;
655  long long int numTasks=qureg.numAmpsPerChunk>>2;
656 
657  // set dimensions
658  sizeInnerHalfBlockQ1 = 1LL << targetQubit;
659  sizeInnerHalfBlockQ2 = 1LL << qubit2;
660  sizeInnerQuarterBlockQ2 = sizeInnerHalfBlockQ2 >> 1;
661  sizeInnerBlockQ2 = sizeInnerHalfBlockQ2 << 1;
662  sizeInnerBlockQ1 = 2LL * sizeInnerHalfBlockQ1;
663  sizeOuterColumn = 1LL << qureg.numQubitsRepresented;
664  sizeOuterQuarterColumn = sizeOuterColumn >> 2;
665 
666 //# if 0
667 # ifdef _OPENMP
668 # pragma omp parallel \
669  default (none) \
670  shared (sizeInnerBlockQ1,sizeInnerHalfBlockQ1,sizeInnerBlockQ2,sizeInnerHalfBlockQ2,sizeInnerQuarterBlockQ2,\
671  sizeOuterColumn,sizeOuterQuarterColumn,qureg,delta,gamma, numTasks,targetQubit,qubit2) \
672  private (thisTask,thisInnerBlockQ2,thisInnerBlockQ1InInnerBlockQ2, \
673  thisOuterColumn,thisIndex,thisIndexInPairVector,thisIndexInOuterColumn, \
674  thisIndexInInnerBlockQ1,thisIndexInInnerBlockQ2,outerBitQ1,outerBitQ2)
675 # endif
676  {
677 # ifdef _OPENMP
678 # pragma omp for schedule (static)
679 # endif
680 //# endif
681  // thisTask iterates over half the elements in this process' chunk of the density matrix
682  // treat this as iterating over all columns, then iterating over half the values
683  // within one column.
684  // If this function has been called, this process' chunk contains half an
685  // outer block or less
686  for (thisTask=0; thisTask<numTasks; thisTask++) {
687  // we want to process all columns in the density matrix,
688  // updating the values for half of each column (one half of each inner block)
689  thisOuterColumn = thisTask / sizeOuterQuarterColumn;
690  // thisTask % sizeOuterQuarterColumn
691  thisIndexInOuterColumn = thisTask&(sizeOuterQuarterColumn-1);
692  thisInnerBlockQ2 = thisIndexInOuterColumn / sizeInnerQuarterBlockQ2;
693  // thisTask % sizeInnerQuarterBlockQ2;
694  thisIndexInInnerBlockQ2 = thisTask&(sizeInnerQuarterBlockQ2-1);
695  thisInnerBlockQ1InInnerBlockQ2 = thisIndexInInnerBlockQ2 / sizeInnerHalfBlockQ1;
696  // thisTask % sizeInnerHalfBlockQ1;
697  thisIndexInInnerBlockQ1 = thisTask&(sizeInnerHalfBlockQ1-1);
698 
699  // get index in state vector corresponding to upper inner block
700  thisIndex = thisOuterColumn*sizeOuterColumn + thisInnerBlockQ2*sizeInnerBlockQ2
701  + thisInnerBlockQ1InInnerBlockQ2*sizeInnerBlockQ1 + thisIndexInInnerBlockQ1;
702 
703  // check if we are in the upper or lower half of an outer block for Q1
704  outerBitQ1 = extractBit(targetQubit, (thisIndex+qureg.numAmpsPerChunk*qureg.chunkId)>>qureg.numQubitsRepresented);
705  // if we are in the lower half of an outer block, shift to be in the lower half
706  // of the inner block as well (we want to dephase |0><0| and |1><1| only)
707  thisIndex += outerBitQ1*(sizeInnerHalfBlockQ1);
708 
709  // For part 3 we need to match elements such that (my Q1 != pair Q1) AND (my Q2 != pair Q2)
710  // Find correct index in pairStateVector
711  thisIndexInPairVector = thisTask + (1-outerBitQ1)*sizeInnerHalfBlockQ1*sizeOuterQuarterColumn -
712  outerBitQ1*sizeInnerHalfBlockQ1*sizeOuterQuarterColumn;
713 
714  // check if we are in the upper or lower half of an outer block for Q2
715  outerBitQ2 = extractBit(qubit2, (thisIndex+qureg.numAmpsPerChunk*qureg.chunkId)>>qureg.numQubitsRepresented);
716  // if we are in the lower half of an outer block, shift to be in the lower half
717  // of the inner block as well (we want to dephase |0><0| and |1><1| only)
718  thisIndex += outerBitQ2*(sizeInnerQuarterBlockQ2<<1);
719 
720 
721  // NOTE: at this point thisIndex should be the index of the element we want to
722  // dephase in the chunk of the state vector on this process, in the
723  // density matrix representation.
724 
725 
726  // state[thisIndex] = (1-depolLevel)*state[thisIndex] + depolLevel*(state[thisIndex]
727  // + pair[thisIndexInPairVector])/2
728  qureg.stateVec.real[thisIndex] = gamma*(qureg.stateVec.real[thisIndex] +
729  delta*qureg.pairStateVec.real[thisIndexInPairVector]);
730 
731  qureg.stateVec.imag[thisIndex] = gamma*(qureg.stateVec.imag[thisIndex] +
732  delta*qureg.pairStateVec.imag[thisIndexInPairVector]);
733  }
734  }
735 
736 }

References Qureg::chunkId, extractBit(), Qureg::numAmpsPerChunk, Qureg::numQubitsRepresented, Qureg::pairStateVec, and Qureg::stateVec.

Referenced by densmatr_mixTwoQubitDepolarising().

◆ extractBit()

static int extractBit ( const int  locationOfBitFromRight,
const long long int  theEncodedNumber 
)
inlinestatic

Definition at line 26 of file QuEST_cpu_internal.h.

26  {
27  return (theEncodedNumber & ( 1LL << locationOfBitFromRight )) >> locationOfBitFromRight;
28 }

Referenced by isOddParity().

◆ flipBit()

static long long int flipBit ( const long long int  number,
const int  bitInd 
)
inlinestatic

Definition at line 30 of file QuEST_cpu_internal.h.

30  {
31  return (number ^ (1LL << bitInd));
32 }

◆ insertTwoZeroBits()

static long long int insertTwoZeroBits ( const long long int  number,
const int  bit1,
const int  bit2 
)
inlinestatic

Definition at line 49 of file QuEST_cpu_internal.h.

49  {
50  int small = (bit1 < bit2)? bit1 : bit2;
51  int big = (bit1 < bit2)? bit2 : bit1;
52  return insertZeroBit(insertZeroBit(number, small), big);
53 }

References insertZeroBit().

◆ insertZeroBit()

static long long int insertZeroBit ( const long long int  number,
const int  index 
)
inlinestatic

Definition at line 42 of file QuEST_cpu_internal.h.

42  {
43  long long int left, right;
44  left = (number >> index) << index;
45  right = number - left;
46  return (left << 1) ^ right;
47 }

Referenced by insertTwoZeroBits().

◆ isOddParity()

static int isOddParity ( const long long int  number,
const int  qb1,
const int  qb2 
)
inlinestatic

Definition at line 38 of file QuEST_cpu_internal.h.

38  {
39  return extractBit(qb1, number) != extractBit(qb2, number);
40 }

References extractBit().

Referenced by statevec_swapQubitAmpsDistributed().

◆ maskContainsBit()

static int maskContainsBit ( const long long int  mask,
const int  bitInd 
)
inlinestatic

Definition at line 34 of file QuEST_cpu_internal.h.

34  {
35  return mask & (1LL << bitInd);
36 }

Referenced by statevec_multiControlledMultiQubitUnitary(), and statevec_multiControlledTwoQubitUnitary().

◆ statevec_calcExpecDiagonalOpLocal()

Complex statevec_calcExpecDiagonalOpLocal ( Qureg  qureg,
DiagonalOp  op 
)

Definition at line 4124 of file QuEST_cpu.c.

4124  {
4125 
4126  qreal expecRe = 0;
4127  qreal expecIm = 0;
4128 
4129  long long int index;
4130  long long int numAmps = qureg.numAmpsPerChunk;
4131  qreal *stateReal = qureg.stateVec.real;
4132  qreal *stateImag = qureg.stateVec.imag;
4133  qreal *opReal = op.real;
4134  qreal *opImag = op.imag;
4135 
4136  qreal vecRe,vecIm,vecAbs, opRe, opIm;
4137 
4138 # ifdef _OPENMP
4139 # pragma omp parallel \
4140  shared (stateReal, stateImag, opReal, opImag, numAmps) \
4141  private (index, vecRe,vecIm,vecAbs, opRe,opIm) \
4142  reduction ( +:expecRe, expecIm )
4143 # endif
4144  {
4145 # ifdef _OPENMP
4146 # pragma omp for schedule (static)
4147 # endif
4148  for (index=0; index < numAmps; index++) {
4149  vecRe = stateReal[index];
4150  vecIm = stateImag[index];
4151  opRe = opReal[index];
4152  opIm = opImag[index];
4153 
4154  // abs(vec)^2 op
4155  vecAbs = vecRe*vecRe + vecIm*vecIm;
4156  expecRe += vecAbs*opRe;
4157  expecIm += vecAbs*opIm;
4158  }
4159  }
4160 
4161  Complex innerProd;
4162  innerProd.real = expecRe;
4163  innerProd.imag = expecIm;
4164  return innerProd;
4165 }

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

Referenced by statevec_calcExpecDiagonalOp().

◆ statevec_calcInnerProductLocal()

Complex statevec_calcInnerProductLocal ( Qureg  bra,
Qureg  ket 
)

Definition at line 1082 of file QuEST_cpu.c.

1082  {
1083 
1084  qreal innerProdReal = 0;
1085  qreal innerProdImag = 0;
1086 
1087  long long int index;
1088  long long int numAmps = bra.numAmpsPerChunk;
1089  qreal *braVecReal = bra.stateVec.real;
1090  qreal *braVecImag = bra.stateVec.imag;
1091  qreal *ketVecReal = ket.stateVec.real;
1092  qreal *ketVecImag = ket.stateVec.imag;
1093 
1094  qreal braRe, braIm, ketRe, ketIm;
1095 
1096 # ifdef _OPENMP
1097 # pragma omp parallel \
1098  shared (braVecReal, braVecImag, ketVecReal, ketVecImag, numAmps) \
1099  private (index, braRe, braIm, ketRe, ketIm) \
1100  reduction ( +:innerProdReal, innerProdImag )
1101 # endif
1102  {
1103 # ifdef _OPENMP
1104 # pragma omp for schedule (static)
1105 # endif
1106  for (index=0; index < numAmps; index++) {
1107  braRe = braVecReal[index];
1108  braIm = braVecImag[index];
1109  ketRe = ketVecReal[index];
1110  ketIm = ketVecImag[index];
1111 
1112  // conj(bra_i) * ket_i
1113  innerProdReal += braRe*ketRe + braIm*ketIm;
1114  innerProdImag += braRe*ketIm - braIm*ketRe;
1115  }
1116  }
1117 
1118  Complex innerProd;
1119  innerProd.real = innerProdReal;
1120  innerProd.imag = innerProdImag;
1121  return innerProd;
1122 }

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

Referenced by statevec_calcInnerProduct().

◆ statevec_calcProbOfAllOutcomesLocal()

void statevec_calcProbOfAllOutcomesLocal ( qreal *  retProbs,
Qureg  qureg,
int *  qubits,
int  numQubits 
)

Definition at line 3549 of file QuEST_cpu.c.

3549  {
3550 
3551  /* Below, we manually reduce amplitudes into outcomeProbs by using atomic update.
3552  * This maintains OpenMP 3.1 compatibility. An alternative is to use array reduction
3553  * (requires OpenMP 4.5, limits #qubits since outcomeProbs must be a local stack array)
3554  * or a dynamic list of omp locks (duplicates memory cost of outcomeProbs).
3555  * Using locks was always slower than the method below. Using reduction was only
3556  * faster for very few threads, or very few outcomeProbs.
3557  * Finally, we exclude the 'update' clause after 'atomic' to maintain MSVC compatibility
3558  */
3559 
3560  long long int numOutcomeProbs = (1 << numQubits);
3561  long long int j;
3562 
3563  // clear outcomeProbs (in parallel, in case it's large)
3564 # ifdef _OPENMP
3565 # pragma omp parallel \
3566  default (none) \
3567  shared (numOutcomeProbs,outcomeProbs) \
3568  private (j)
3569 # endif
3570  {
3571 # ifdef _OPENMP
3572 # pragma omp for schedule (static)
3573 # endif
3574  for (j=0; j<numOutcomeProbs; j++)
3575  outcomeProbs[j] = 0;
3576  }
3577 
3578  long long int numTasks = qureg.numAmpsPerChunk;
3579  long long int offset = qureg.chunkId*qureg.numAmpsPerChunk;
3580  qreal* stateRe = qureg.stateVec.real;
3581  qreal* stateIm = qureg.stateVec.imag;
3582 
3583  long long int i;
3584  long long int outcomeInd;
3585  int q;
3586  qreal prob;
3587 
3588 # ifdef _OPENMP
3589 # pragma omp parallel \
3590  shared (numTasks,offset, qubits,numQubits, stateRe,stateIm, outcomeProbs) \
3591  private (i, q, outcomeInd, prob)
3592 # endif
3593  {
3594 # ifdef _OPENMP
3595 # pragma omp for schedule (static)
3596 # endif
3597  // every amplitude contributes to a single element of retProbs
3598  for (i=0; i<numTasks; i++) {
3599 
3600  // determine index informed by qubits outcome
3601  outcomeInd = 0;
3602  for (q=0; q<numQubits; q++)
3603  outcomeInd += extractBit(qubits[q], i + offset) * (1LL << q);
3604 
3605  prob = stateRe[i]*stateRe[i] + stateIm[i]*stateIm[i];
3606 
3607  // atomicly update corresponding outcome array element
3608  # ifdef _OPENMP
3609  # pragma omp atomic
3610  # endif
3611  outcomeProbs[outcomeInd] += prob;
3612  }
3613  }
3614 }

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

Referenced by statevec_calcProbOfAllOutcomes().

◆ statevec_collapseToKnownProbOutcomeDistributedRenorm()

void statevec_collapseToKnownProbOutcomeDistributedRenorm ( Qureg  qureg,
int  measureQubit,
qreal  totalProbability 
)

Renormalise parts of the state vector where measureQubit=0 or 1, based on the total probability of that qubit being in state 0 or 1.

Measure in Zero performs an irreversible change to the state vector: it updates the vector according to the event that the value 'outcome' has been measured on the qubit indicated by measureQubit (where this label starts from 0, of course). It achieves this by setting all inconsistent amplitudes to 0 and then renormalising based on the total probability of measuring measureQubit=0 if outcome=0 and measureQubit=1 if outcome=1. In the distributed version, one block (with measureQubit=0 in the first half of the block and measureQubit=1 in the second half of the block) is spread over multiple chunks, meaning that each chunks performs only renormalisation or only setting amplitudes to 0. This function handles the renormalisation.

Parameters
[in,out]quregobject representing the set of qubits
[in]measureQubitqubit to measure
[in]totalProbabilityprobability of qubit measureQubit being zero

Definition at line 3849 of file QuEST_cpu.c.

3850 {
3851  // ----- temp variables
3852  long long int thisTask;
3853  long long int numTasks=qureg.numAmpsPerChunk;
3854 
3855  qreal renorm=1/sqrt(totalProbability);
3856 
3857  qreal *stateVecReal = qureg.stateVec.real;
3858  qreal *stateVecImag = qureg.stateVec.imag;
3859 
3860 # ifdef _OPENMP
3861 # pragma omp parallel \
3862  shared (numTasks,stateVecReal,stateVecImag) \
3863  private (thisTask)
3864 # endif
3865  {
3866 # ifdef _OPENMP
3867 # pragma omp for schedule (static)
3868 # endif
3869  for (thisTask=0; thisTask<numTasks; thisTask++) {
3870  stateVecReal[thisTask] = stateVecReal[thisTask]*renorm;
3871  stateVecImag[thisTask] = stateVecImag[thisTask]*renorm;
3872  }
3873  }
3874 }

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

Referenced by statevec_collapseToKnownProbOutcome().

◆ statevec_collapseToKnownProbOutcomeLocal()

void statevec_collapseToKnownProbOutcomeLocal ( Qureg  qureg,
int  measureQubit,
int  outcome,
qreal  totalProbability 
)

Update the state vector to be consistent with measuring measureQubit=0 if outcome=0 and measureQubit=1 if outcome=1.

Performs an irreversible change to the state vector: it updates the vector according to the event that an outcome have been measured on the qubit indicated by measureQubit (where this label starts from 0, of course). It achieves this by setting all inconsistent amplitudes to 0 and then renormalising based on the total probability of measuring measureQubit=0 or 1 according to the value of outcome. In the local version, one or more blocks (with measureQubit=0 in the first half of the block and measureQubit=1 in the second half of the block) fit entirely into one chunk.

Parameters
[in,out]quregobject representing the set of qubits
[in]measureQubitqubit to measure
[in]totalProbabilityprobability of qubit measureQubit being either zero or one
[in]outcometo measure the probability of and set the state to – either zero or one

Definition at line 3767 of file QuEST_cpu.c.

3768 {
3769  // ----- sizes
3770  long long int sizeBlock, // size of blocks
3771  sizeHalfBlock; // size of blocks halved
3772  // ----- indices
3773  long long int thisBlock, // current block
3774  index; // current index for first half block
3775  // ----- measured probability
3776  qreal renorm; // probability (returned) value
3777  // ----- temp variables
3778  long long int thisTask; // task based approach for expose loop with small granularity
3779  // (good for shared memory parallelism)
3780  long long int numTasks=qureg.numAmpsPerChunk>>1;
3781 
3782  // ---------------------------------------------------------------- //
3783  // dimensions //
3784  // ---------------------------------------------------------------- //
3785  sizeHalfBlock = 1LL << (measureQubit); // number of state vector elements to sum,
3786  // and then the number to skip
3787  sizeBlock = 2LL * sizeHalfBlock; // size of blocks (pairs of measure and skip entries)
3788 
3789  renorm=1/sqrt(totalProbability);
3790  qreal *stateVecReal = qureg.stateVec.real;
3791  qreal *stateVecImag = qureg.stateVec.imag;
3792 
3793 
3794 # ifdef _OPENMP
3795 # pragma omp parallel \
3796  default (none) \
3797  shared (numTasks,sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag,renorm,outcome) \
3798  private (thisTask,thisBlock,index)
3799 # endif
3800  {
3801  if (outcome==0){
3802  // measure qubit is 0
3803 # ifdef _OPENMP
3804 # pragma omp for schedule (static)
3805 # endif
3806  for (thisTask=0; thisTask<numTasks; thisTask++) {
3807  thisBlock = thisTask / sizeHalfBlock;
3808  index = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
3809  stateVecReal[index]=stateVecReal[index]*renorm;
3810  stateVecImag[index]=stateVecImag[index]*renorm;
3811 
3812  stateVecReal[index+sizeHalfBlock]=0;
3813  stateVecImag[index+sizeHalfBlock]=0;
3814  }
3815  } else {
3816  // measure qubit is 1
3817 # ifdef _OPENMP
3818 # pragma omp for schedule (static)
3819 # endif
3820  for (thisTask=0; thisTask<numTasks; thisTask++) {
3821  thisBlock = thisTask / sizeHalfBlock;
3822  index = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
3823  stateVecReal[index]=0;
3824  stateVecImag[index]=0;
3825 
3826  stateVecReal[index+sizeHalfBlock]=stateVecReal[index+sizeHalfBlock]*renorm;
3827  stateVecImag[index+sizeHalfBlock]=stateVecImag[index+sizeHalfBlock]*renorm;
3828  }
3829  }
3830  }
3831 
3832 }

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

Referenced by statevec_collapseToKnownProbOutcome().

◆ statevec_collapseToOutcomeDistributedSetZero()

void statevec_collapseToOutcomeDistributedSetZero ( Qureg  qureg)

Definition at line 3887 of file QuEST_cpu.c.

3888 {
3889  // ----- temp variables
3890  long long int thisTask;
3891  long long int numTasks=qureg.numAmpsPerChunk;
3892 
3893  // ---------------------------------------------------------------- //
3894  // find probability //
3895  // ---------------------------------------------------------------- //
3896 
3897  qreal *stateVecReal = qureg.stateVec.real;
3898  qreal *stateVecImag = qureg.stateVec.imag;
3899 
3900 # ifdef _OPENMP
3901 # pragma omp parallel \
3902  shared (numTasks,stateVecReal,stateVecImag) \
3903  private (thisTask)
3904 # endif
3905  {
3906 # ifdef _OPENMP
3907 # pragma omp for schedule (static)
3908 # endif
3909  for (thisTask=0; thisTask<numTasks; thisTask++) {
3910  stateVecReal[thisTask] = 0;
3911  stateVecImag[thisTask] = 0;
3912  }
3913  }
3914 }

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

Referenced by statevec_collapseToKnownProbOutcome().

◆ statevec_compactUnitaryDistributed()

void statevec_compactUnitaryDistributed ( Qureg  qureg,
Complex  rot1,
Complex  rot2,
ComplexArray  stateVecUp,
ComplexArray  stateVecLo,
ComplexArray  stateVecOut 
)

Rotate a single qubit in the state vector of probability amplitudes, given two complex numbers alpha and beta, and a subset of the state vector with upper and lower block values stored seperately.

Parameters
[in,out]quregobject representing the set of qubits
[in]rot1rotation angle
[in]rot2rotation angle
[in]stateVecUpprobability amplitudes in upper half of a block
[in]stateVecLoprobability amplitudes in lower half of a block
[out]stateVecOutarray section to update (will correspond to either the lower or upper half of a block)

Definition at line 2095 of file QuEST_cpu.c.

2100 {
2101 
2102  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2103  long long int thisTask;
2104  long long int numTasks=qureg.numAmpsPerChunk;
2105 
2106  qreal rot1Real=rot1.real, rot1Imag=rot1.imag;
2107  qreal rot2Real=rot2.real, rot2Imag=rot2.imag;
2108  qreal *stateVecRealUp=stateVecUp.real, *stateVecImagUp=stateVecUp.imag;
2109  qreal *stateVecRealLo=stateVecLo.real, *stateVecImagLo=stateVecLo.imag;
2110  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2111 
2112 # ifdef _OPENMP
2113 # pragma omp parallel \
2114  default (none) \
2115  shared (stateVecRealUp,stateVecImagUp,stateVecRealLo,stateVecImagLo,stateVecRealOut,stateVecImagOut, \
2116  rot1Real,rot1Imag, rot2Real,rot2Imag,numTasks) \
2117  private (thisTask,stateRealUp,stateImagUp,stateRealLo,stateImagLo)
2118 # endif
2119  {
2120 # ifdef _OPENMP
2121 # pragma omp for schedule (static)
2122 # endif
2123  for (thisTask=0; thisTask<numTasks; thisTask++) {
2124  // store current state vector values in temp variables
2125  stateRealUp = stateVecRealUp[thisTask];
2126  stateImagUp = stateVecImagUp[thisTask];
2127 
2128  stateRealLo = stateVecRealLo[thisTask];
2129  stateImagLo = stateVecImagLo[thisTask];
2130 
2131  // state[indexUp] = alpha * state[indexUp] - conj(beta) * state[indexLo]
2132  stateVecRealOut[thisTask] = rot1Real*stateRealUp - rot1Imag*stateImagUp + rot2Real*stateRealLo + rot2Imag*stateImagLo;
2133  stateVecImagOut[thisTask] = rot1Real*stateImagUp + rot1Imag*stateRealUp + rot2Real*stateImagLo - rot2Imag*stateRealLo;
2134  }
2135  }
2136 }

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

Referenced by statevec_compactUnitary().

◆ statevec_compactUnitaryLocal()

void statevec_compactUnitaryLocal ( Qureg  qureg,
int  targetQubit,
Complex  alpha,
Complex  beta 
)

Definition at line 1754 of file QuEST_cpu.c.

1755 {
1756  long long int sizeBlock, sizeHalfBlock;
1757  long long int thisBlock, // current block
1758  indexUp,indexLo; // current index and corresponding index in lower half block
1759 
1760  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
1761  long long int thisTask;
1762  long long int numTasks=qureg.numAmpsPerChunk>>1;
1763 
1764  // set dimensions
1765  sizeHalfBlock = 1LL << targetQubit;
1766  sizeBlock = 2LL * sizeHalfBlock;
1767 
1768  // Can't use qureg.stateVec as a private OMP var
1769  qreal *stateVecReal = qureg.stateVec.real;
1770  qreal *stateVecImag = qureg.stateVec.imag;
1771  qreal alphaImag=alpha.imag, alphaReal=alpha.real;
1772  qreal betaImag=beta.imag, betaReal=beta.real;
1773 
1774 # ifdef _OPENMP
1775 # pragma omp parallel \
1776  default (none) \
1777  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, alphaReal,alphaImag, betaReal,betaImag, numTasks) \
1778  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp,stateRealLo,stateImagLo)
1779 # endif
1780  {
1781 # ifdef _OPENMP
1782 # pragma omp for schedule (static)
1783 # endif
1784  for (thisTask=0; thisTask<numTasks; thisTask++) {
1785 
1786  thisBlock = thisTask / sizeHalfBlock;
1787  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
1788  indexLo = indexUp + sizeHalfBlock;
1789 
1790  // store current state vector values in temp variables
1791  stateRealUp = stateVecReal[indexUp];
1792  stateImagUp = stateVecImag[indexUp];
1793 
1794  stateRealLo = stateVecReal[indexLo];
1795  stateImagLo = stateVecImag[indexLo];
1796 
1797  // state[indexUp] = alpha * state[indexUp] - conj(beta) * state[indexLo]
1798  stateVecReal[indexUp] = alphaReal*stateRealUp - alphaImag*stateImagUp
1799  - betaReal*stateRealLo - betaImag*stateImagLo;
1800  stateVecImag[indexUp] = alphaReal*stateImagUp + alphaImag*stateRealUp
1801  - betaReal*stateImagLo + betaImag*stateRealLo;
1802 
1803  // state[indexLo] = beta * state[indexUp] + conj(alpha) * state[indexLo]
1804  stateVecReal[indexLo] = betaReal*stateRealUp - betaImag*stateImagUp
1805  + alphaReal*stateRealLo + alphaImag*stateImagLo;
1806  stateVecImag[indexLo] = betaReal*stateImagUp + betaImag*stateRealUp
1807  + alphaReal*stateImagLo - alphaImag*stateRealLo;
1808  }
1809  }
1810 
1811 }

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

Referenced by statevec_compactUnitary().

◆ statevec_controlledCompactUnitaryDistributed()

void statevec_controlledCompactUnitaryDistributed ( Qureg  qureg,
int  controlQubit,
Complex  rot1,
Complex  rot2,
ComplexArray  stateVecUp,
ComplexArray  stateVecLo,
ComplexArray  stateVecOut 
)

Rotate a single qubit in the state vector of probability amplitudes, given two complex numbers alpha and beta and a subset of the state vector with upper and lower block values stored seperately.

Only perform the rotation where the control qubit is one.

Parameters
[in,out]quregobject representing the set of qubits
[in]controlQubitqubit to determine whether or not to perform a rotation
[in]rot1rotation angle
[in]rot2rotation angle
[in]stateVecUpprobability amplitudes in upper half of a block
[in]stateVecLoprobability amplitudes in lower half of a block
[out]stateVecOutarray section to update (will correspond to either the lower or upper half of a block)

Definition at line 2414 of file QuEST_cpu.c.

2419 {
2420 
2421  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2422  long long int thisTask;
2423  long long int numTasks=qureg.numAmpsPerChunk;
2424  long long int chunkSize=qureg.numAmpsPerChunk;
2425  long long int chunkId=qureg.chunkId;
2426 
2427  int controlBit;
2428 
2429  qreal rot1Real=rot1.real, rot1Imag=rot1.imag;
2430  qreal rot2Real=rot2.real, rot2Imag=rot2.imag;
2431  qreal *stateVecRealUp=stateVecUp.real, *stateVecImagUp=stateVecUp.imag;
2432  qreal *stateVecRealLo=stateVecLo.real, *stateVecImagLo=stateVecLo.imag;
2433  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2434 
2435 # ifdef _OPENMP
2436 # pragma omp parallel \
2437  default (none) \
2438  shared (stateVecRealUp,stateVecImagUp,stateVecRealLo,stateVecImagLo,stateVecRealOut,stateVecImagOut, \
2439  rot1Real,rot1Imag, rot2Real,rot2Imag,numTasks,chunkId,chunkSize,controlQubit) \
2440  private (thisTask,stateRealUp,stateImagUp,stateRealLo,stateImagLo,controlBit)
2441 # endif
2442  {
2443 # ifdef _OPENMP
2444 # pragma omp for schedule (static)
2445 # endif
2446  for (thisTask=0; thisTask<numTasks; thisTask++) {
2447  controlBit = extractBit (controlQubit, thisTask+chunkId*chunkSize);
2448  if (controlBit){
2449  // store current state vector values in temp variables
2450  stateRealUp = stateVecRealUp[thisTask];
2451  stateImagUp = stateVecImagUp[thisTask];
2452 
2453  stateRealLo = stateVecRealLo[thisTask];
2454  stateImagLo = stateVecImagLo[thisTask];
2455 
2456  // state[indexUp] = alpha * state[indexUp] - conj(beta) * state[indexLo]
2457  stateVecRealOut[thisTask] = rot1Real*stateRealUp - rot1Imag*stateImagUp + rot2Real*stateRealLo + rot2Imag*stateImagLo;
2458  stateVecImagOut[thisTask] = rot1Real*stateImagUp + rot1Imag*stateRealUp + rot2Real*stateImagLo - rot2Imag*stateRealLo;
2459  }
2460  }
2461  }
2462 }

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

Referenced by statevec_controlledCompactUnitary().

◆ statevec_controlledCompactUnitaryLocal()

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

Definition at line 2196 of file QuEST_cpu.c.

2198 {
2199  long long int sizeBlock, sizeHalfBlock;
2200  long long int thisBlock, // current block
2201  indexUp,indexLo; // current index and corresponding index in lower half block
2202 
2203  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2204  long long int thisTask;
2205  long long int numTasks=qureg.numAmpsPerChunk>>1;
2206  long long int chunkSize=qureg.numAmpsPerChunk;
2207  long long int chunkId=qureg.chunkId;
2208 
2209  int controlBit;
2210 
2211  // set dimensions
2212  sizeHalfBlock = 1LL << targetQubit;
2213  sizeBlock = 2LL * sizeHalfBlock;
2214 
2215  // Can't use qureg.stateVec as a private OMP var
2216  qreal *stateVecReal = qureg.stateVec.real;
2217  qreal *stateVecImag = qureg.stateVec.imag;
2218  qreal alphaImag=alpha.imag, alphaReal=alpha.real;
2219  qreal betaImag=beta.imag, betaReal=beta.real;
2220 
2221 # ifdef _OPENMP
2222 # pragma omp parallel \
2223  default (none) \
2224  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, alphaReal,alphaImag, betaReal,betaImag, \
2225  numTasks,chunkId,chunkSize,controlQubit) \
2226  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp,stateRealLo,stateImagLo,controlBit)
2227 # endif
2228  {
2229 # ifdef _OPENMP
2230 # pragma omp for schedule (static)
2231 # endif
2232  for (thisTask=0; thisTask<numTasks; thisTask++) {
2233 
2234  thisBlock = thisTask / sizeHalfBlock;
2235  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2236  indexLo = indexUp + sizeHalfBlock;
2237 
2238  controlBit = extractBit (controlQubit, indexUp+chunkId*chunkSize);
2239  if (controlBit){
2240  // store current state vector values in temp variables
2241  stateRealUp = stateVecReal[indexUp];
2242  stateImagUp = stateVecImag[indexUp];
2243 
2244  stateRealLo = stateVecReal[indexLo];
2245  stateImagLo = stateVecImag[indexLo];
2246 
2247  // state[indexUp] = alpha * state[indexUp] - conj(beta) * state[indexLo]
2248  stateVecReal[indexUp] = alphaReal*stateRealUp - alphaImag*stateImagUp
2249  - betaReal*stateRealLo - betaImag*stateImagLo;
2250  stateVecImag[indexUp] = alphaReal*stateImagUp + alphaImag*stateRealUp
2251  - betaReal*stateImagLo + betaImag*stateRealLo;
2252 
2253  // state[indexLo] = beta * state[indexUp] + conj(alpha) * state[indexLo]
2254  stateVecReal[indexLo] = betaReal*stateRealUp - betaImag*stateImagUp
2255  + alphaReal*stateRealLo + alphaImag*stateImagLo;
2256  stateVecImag[indexLo] = betaReal*stateImagUp + betaImag*stateRealUp
2257  + alphaReal*stateImagLo - alphaImag*stateRealLo;
2258  }
2259  }
2260  }
2261 
2262 }

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

Referenced by statevec_controlledCompactUnitary().

◆ statevec_controlledNotDistributed()

void statevec_controlledNotDistributed ( Qureg  qureg,
int  controlQubit,
ComplexArray  stateVecIn,
ComplexArray  stateVecOut 
)

Rotate a single qubit by {{0,1},{1,0}.

Operate on a subset of the state vector with upper and lower block values stored seperately. This rotation is just swapping upper and lower values, and stateVecIn must already be the correct section for this chunk. Only perform the rotation for elements where controlQubit is one.

Parameters
[in,out]quregobject representing the set of qubits
[in]controlQubitthe control qubit
[in]stateVecInprobability amplitudes in lower or upper half of a block depending on chunkId
[out]stateVecOutarray section to update (will correspond to either the lower or upper half of a block)

Definition at line 2742 of file QuEST_cpu.c.

2745 {
2746 
2747  long long int thisTask;
2748  long long int numTasks=qureg.numAmpsPerChunk;
2749  long long int chunkSize=qureg.numAmpsPerChunk;
2750  long long int chunkId=qureg.chunkId;
2751 
2752  int controlBit;
2753 
2754  qreal *stateVecRealIn=stateVecIn.real, *stateVecImagIn=stateVecIn.imag;
2755  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2756 
2757 # ifdef _OPENMP
2758 # pragma omp parallel \
2759  default (none) \
2760  shared (stateVecRealIn,stateVecImagIn,stateVecRealOut,stateVecImagOut, \
2761  numTasks,chunkId,chunkSize,controlQubit) \
2762  private (thisTask,controlBit)
2763 # endif
2764  {
2765 # ifdef _OPENMP
2766 # pragma omp for schedule (static)
2767 # endif
2768  for (thisTask=0; thisTask<numTasks; thisTask++) {
2769  controlBit = extractBit (controlQubit, thisTask+chunkId*chunkSize);
2770  if (controlBit){
2771  stateVecRealOut[thisTask] = stateVecRealIn[thisTask];
2772  stateVecImagOut[thisTask] = stateVecImagIn[thisTask];
2773  }
2774  }
2775  }
2776 }

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

Referenced by statevec_controlledNot().

◆ statevec_controlledNotLocal()

void statevec_controlledNotLocal ( Qureg  qureg,
int  controlQubit,
int  targetQubit 
)

Definition at line 2679 of file QuEST_cpu.c.

2680 {
2681  long long int sizeBlock, sizeHalfBlock;
2682  long long int thisBlock, // current block
2683  indexUp,indexLo; // current index and corresponding index in lower half block
2684 
2685  qreal stateRealUp,stateImagUp;
2686  long long int thisTask;
2687  long long int numTasks=qureg.numAmpsPerChunk>>1;
2688  long long int chunkSize=qureg.numAmpsPerChunk;
2689  long long int chunkId=qureg.chunkId;
2690 
2691  int controlBit;
2692 
2693  // set dimensions
2694  sizeHalfBlock = 1LL << targetQubit;
2695  sizeBlock = 2LL * sizeHalfBlock;
2696 
2697  // Can't use qureg.stateVec as a private OMP var
2698  qreal *stateVecReal = qureg.stateVec.real;
2699  qreal *stateVecImag = qureg.stateVec.imag;
2700 
2701 # ifdef _OPENMP
2702 # pragma omp parallel \
2703  default (none) \
2704  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag,numTasks,chunkId,chunkSize,controlQubit) \
2705  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp,controlBit)
2706 # endif
2707  {
2708 # ifdef _OPENMP
2709 # pragma omp for schedule (static)
2710 # endif
2711  for (thisTask=0; thisTask<numTasks; thisTask++) {
2712  thisBlock = thisTask / sizeHalfBlock;
2713  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2714  indexLo = indexUp + sizeHalfBlock;
2715 
2716  controlBit = extractBit(controlQubit, indexUp+chunkId*chunkSize);
2717  if (controlBit){
2718  stateRealUp = stateVecReal[indexUp];
2719  stateImagUp = stateVecImag[indexUp];
2720 
2721  stateVecReal[indexUp] = stateVecReal[indexLo];
2722  stateVecImag[indexUp] = stateVecImag[indexLo];
2723 
2724  stateVecReal[indexLo] = stateRealUp;
2725  stateVecImag[indexLo] = stateImagUp;
2726  }
2727  }
2728  }
2729 }

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

Referenced by statevec_controlledNot().

◆ statevec_controlledPauliYDistributed()

void statevec_controlledPauliYDistributed ( Qureg  qureg,
int  controlQubit,
ComplexArray  stateVecIn,
ComplexArray  stateVecOut,
int  conjFactor 
)

Definition at line 3036 of file QuEST_cpu.c.

3039 {
3040 
3041  long long int thisTask;
3042  long long int numTasks=qureg.numAmpsPerChunk;
3043  long long int chunkSize=qureg.numAmpsPerChunk;
3044  long long int chunkId=qureg.chunkId;
3045 
3046  int controlBit;
3047 
3048  qreal *stateVecRealIn=stateVecIn.real, *stateVecImagIn=stateVecIn.imag;
3049  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
3050 
3051 # ifdef _OPENMP
3052 # pragma omp parallel \
3053  default (none) \
3054  shared (stateVecRealIn,stateVecImagIn,stateVecRealOut,stateVecImagOut, \
3055  numTasks,chunkId,chunkSize,controlQubit,conjFac) \
3056  private (thisTask,controlBit)
3057 # endif
3058  {
3059 # ifdef _OPENMP
3060 # pragma omp for schedule (static)
3061 # endif
3062  for (thisTask=0; thisTask<numTasks; thisTask++) {
3063  controlBit = extractBit (controlQubit, thisTask+chunkId*chunkSize);
3064  if (controlBit){
3065  stateVecRealOut[thisTask] = conjFac * stateVecImagIn[thisTask];
3066  stateVecImagOut[thisTask] = conjFac * -stateVecRealIn[thisTask];
3067  }
3068  }
3069  }
3070 }

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

Referenced by statevec_controlledPauliY(), and statevec_controlledPauliYConj().

◆ statevec_controlledPauliYLocal()

void statevec_controlledPauliYLocal ( Qureg  qureg,
int  controlQubit,
int  targetQubit,
int  conjFactor 
)

Definition at line 2982 of file QuEST_cpu.c.

2983 {
2984  long long int sizeBlock, sizeHalfBlock;
2985  long long int thisBlock, // current block
2986  indexUp,indexLo; // current index and corresponding index in lower half block
2987 
2988  qreal stateRealUp,stateImagUp;
2989  long long int thisTask;
2990  long long int numTasks=qureg.numAmpsPerChunk>>1;
2991  long long int chunkSize=qureg.numAmpsPerChunk;
2992  long long int chunkId=qureg.chunkId;
2993 
2994  int controlBit;
2995 
2996  // set dimensions
2997  sizeHalfBlock = 1LL << targetQubit;
2998  sizeBlock = 2LL * sizeHalfBlock;
2999 
3000  // Can't use qureg.stateVec as a private OMP var
3001  qreal *stateVecReal = qureg.stateVec.real;
3002  qreal *stateVecImag = qureg.stateVec.imag;
3003 
3004 # ifdef _OPENMP
3005 # pragma omp parallel \
3006  default (none) \
3007  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, numTasks,chunkId, \
3008  chunkSize,controlQubit,conjFac) \
3009  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp,controlBit)
3010 # endif
3011  {
3012 # ifdef _OPENMP
3013 # pragma omp for schedule (static)
3014 # endif
3015  for (thisTask=0; thisTask<numTasks; thisTask++) {
3016  thisBlock = thisTask / sizeHalfBlock;
3017  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
3018  indexLo = indexUp + sizeHalfBlock;
3019 
3020  controlBit = extractBit(controlQubit, indexUp+chunkId*chunkSize);
3021  if (controlBit){
3022  stateRealUp = stateVecReal[indexUp];
3023  stateImagUp = stateVecImag[indexUp];
3024 
3025  // update under +-{{0, -i}, {i, 0}}
3026  stateVecReal[indexUp] = conjFac * stateVecImag[indexLo];
3027  stateVecImag[indexUp] = conjFac * -stateVecReal[indexLo];
3028  stateVecReal[indexLo] = conjFac * -stateImagUp;
3029  stateVecImag[indexLo] = conjFac * stateRealUp;
3030  }
3031  }
3032  }
3033 }

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

Referenced by statevec_controlledPauliY(), and statevec_controlledPauliYConj().

◆ statevec_controlledUnitaryDistributed()

void statevec_controlledUnitaryDistributed ( Qureg  qureg,
int  controlQubit,
Complex  rot1,
Complex  rot2,
ComplexArray  stateVecUp,
ComplexArray  stateVecLo,
ComplexArray  stateVecOut 
)

Rotate a single qubit in the state vector of probability amplitudes, given two complex numbers alpha and beta and a subset of the state vector with upper and lower block values stored seperately.

Only perform the rotation where the control qubit is one.

Parameters
[in,out]quregobject representing the set of qubits
[in]controlQubitqubit to determine whether or not to perform a rotation
[in]rot1rotation angle
[in]rot2rotation angle
[in]stateVecUpprobability amplitudes in upper half of a block
[in]stateVecLoprobability amplitudes in lower half of a block
[out]stateVecOutarray section to update (will correspond to either the lower or upper half of a block)

Definition at line 2476 of file QuEST_cpu.c.

2481 {
2482 
2483  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2484  long long int thisTask;
2485  long long int numTasks=qureg.numAmpsPerChunk;
2486  long long int chunkSize=qureg.numAmpsPerChunk;
2487  long long int chunkId=qureg.chunkId;
2488 
2489  int controlBit;
2490 
2491  qreal rot1Real=rot1.real, rot1Imag=rot1.imag;
2492  qreal rot2Real=rot2.real, rot2Imag=rot2.imag;
2493  qreal *stateVecRealUp=stateVecUp.real, *stateVecImagUp=stateVecUp.imag;
2494  qreal *stateVecRealLo=stateVecLo.real, *stateVecImagLo=stateVecLo.imag;
2495  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2496 
2497 # ifdef _OPENMP
2498 # pragma omp parallel \
2499  default (none) \
2500  shared (stateVecRealUp,stateVecImagUp,stateVecRealLo,stateVecImagLo,stateVecRealOut,stateVecImagOut, \
2501  rot1Real,rot1Imag, rot2Real,rot2Imag, numTasks,chunkId,chunkSize,controlQubit) \
2502  private (thisTask,stateRealUp,stateImagUp,stateRealLo,stateImagLo,controlBit)
2503 # endif
2504  {
2505 # ifdef _OPENMP
2506 # pragma omp for schedule (static)
2507 # endif
2508  for (thisTask=0; thisTask<numTasks; thisTask++) {
2509  controlBit = extractBit (controlQubit, thisTask+chunkId*chunkSize);
2510  if (controlBit){
2511  // store current state vector values in temp variables
2512  stateRealUp = stateVecRealUp[thisTask];
2513  stateImagUp = stateVecImagUp[thisTask];
2514 
2515  stateRealLo = stateVecRealLo[thisTask];
2516  stateImagLo = stateVecImagLo[thisTask];
2517 
2518  stateVecRealOut[thisTask] = rot1Real*stateRealUp - rot1Imag*stateImagUp
2519  + rot2Real*stateRealLo - rot2Imag*stateImagLo;
2520  stateVecImagOut[thisTask] = rot1Real*stateImagUp + rot1Imag*stateRealUp
2521  + rot2Real*stateImagLo + rot2Imag*stateRealLo;
2522  }
2523  }
2524  }
2525 }

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

Referenced by statevec_controlledUnitary().

◆ statevec_controlledUnitaryLocal()

void statevec_controlledUnitaryLocal ( Qureg  qureg,
int  controlQubit,
int  targetQubit,
ComplexMatrix2  u 
)

Definition at line 2336 of file QuEST_cpu.c.

2338 {
2339  long long int sizeBlock, sizeHalfBlock;
2340  long long int thisBlock, // current block
2341  indexUp,indexLo; // current index and corresponding index in lower half block
2342 
2343  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2344  long long int thisTask;
2345  long long int numTasks=qureg.numAmpsPerChunk>>1;
2346  long long int chunkSize=qureg.numAmpsPerChunk;
2347  long long int chunkId=qureg.chunkId;
2348 
2349  int controlBit;
2350 
2351  // set dimensions
2352  sizeHalfBlock = 1LL << targetQubit;
2353  sizeBlock = 2LL * sizeHalfBlock;
2354 
2355  // Can't use qureg.stateVec as a private OMP var
2356  qreal *stateVecReal = qureg.stateVec.real;
2357  qreal *stateVecImag = qureg.stateVec.imag;
2358 
2359 # ifdef _OPENMP
2360 # pragma omp parallel \
2361  default (none) \
2362  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, u,numTasks,chunkId,chunkSize,controlQubit) \
2363  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp,stateRealLo,stateImagLo,controlBit)
2364 # endif
2365  {
2366 # ifdef _OPENMP
2367 # pragma omp for schedule (static)
2368 # endif
2369  for (thisTask=0; thisTask<numTasks; thisTask++) {
2370 
2371  thisBlock = thisTask / sizeHalfBlock;
2372  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2373  indexLo = indexUp + sizeHalfBlock;
2374 
2375  controlBit = extractBit (controlQubit, indexUp+chunkId*chunkSize);
2376  if (controlBit){
2377  // store current state vector values in temp variables
2378  stateRealUp = stateVecReal[indexUp];
2379  stateImagUp = stateVecImag[indexUp];
2380 
2381  stateRealLo = stateVecReal[indexLo];
2382  stateImagLo = stateVecImag[indexLo];
2383 
2384 
2385  // state[indexUp] = u00 * state[indexUp] + u01 * state[indexLo]
2386  stateVecReal[indexUp] = u.real[0][0]*stateRealUp - u.imag[0][0]*stateImagUp
2387  + u.real[0][1]*stateRealLo - u.imag[0][1]*stateImagLo;
2388  stateVecImag[indexUp] = u.real[0][0]*stateImagUp + u.imag[0][0]*stateRealUp
2389  + u.real[0][1]*stateImagLo + u.imag[0][1]*stateRealLo;
2390 
2391  // state[indexLo] = u10 * state[indexUp] + u11 * state[indexLo]
2392  stateVecReal[indexLo] = u.real[1][0]*stateRealUp - u.imag[1][0]*stateImagUp
2393  + u.real[1][1]*stateRealLo - u.imag[1][1]*stateImagLo;
2394  stateVecImag[indexLo] = u.real[1][0]*stateImagUp + u.imag[1][0]*stateRealUp
2395  + u.real[1][1]*stateImagLo + u.imag[1][1]*stateRealLo;
2396  }
2397  }
2398  }
2399 
2400 }

References Qureg::chunkId, extractBit(), ComplexMatrix2::imag, Qureg::numAmpsPerChunk, qreal, ComplexMatrix2::real, and Qureg::stateVec.

Referenced by statevec_controlledUnitary().

◆ statevec_findProbabilityOfZeroDistributed()

qreal statevec_findProbabilityOfZeroDistributed ( Qureg  qureg)

Measure the probability of a specified qubit being in the zero state across all amplitudes held in this chunk.

Size of regions to skip is a multiple of chunkSize. The results are communicated and aggregated by the caller

Parameters
[in]quregobject representing the set of qubits
Returns
probability of qubit measureQubit being zero

Definition at line 3513 of file QuEST_cpu.c.

3513  {
3514  // ----- measured probability
3515  qreal totalProbability; // probability (returned) value
3516  // ----- temp variables
3517  long long int thisTask; // task based approach for expose loop with small granularity
3518  long long int numTasks=qureg.numAmpsPerChunk;
3519 
3520  // ---------------------------------------------------------------- //
3521  // find probability //
3522  // ---------------------------------------------------------------- //
3523 
3524  // initialise returned value
3525  totalProbability = 0.0;
3526 
3527  qreal *stateVecReal = qureg.stateVec.real;
3528  qreal *stateVecImag = qureg.stateVec.imag;
3529 
3530 # ifdef _OPENMP
3531 # pragma omp parallel \
3532  shared (numTasks,stateVecReal,stateVecImag) \
3533  private (thisTask) \
3534  reduction ( +:totalProbability )
3535 # endif
3536  {
3537 # ifdef _OPENMP
3538 # pragma omp for schedule (static)
3539 # endif
3540  for (thisTask=0; thisTask<numTasks; thisTask++) {
3541  totalProbability += stateVecReal[thisTask]*stateVecReal[thisTask]
3542  + stateVecImag[thisTask]*stateVecImag[thisTask];
3543  }
3544  }
3545 
3546  return totalProbability;
3547 }

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

Referenced by statevec_calcProbOfOutcome().

◆ statevec_findProbabilityOfZeroLocal()

qreal statevec_findProbabilityOfZeroLocal ( Qureg  qureg,
int  measureQubit 
)

Measure the total probability of a specified qubit being in the zero state across all amplitudes in this chunk.

Size of regions to skip is less than the size of one chunk.

Parameters
[in]quregobject representing the set of qubits
[in]measureQubitqubit to measure
Returns
probability of qubit measureQubit being zero

Definition at line 3457 of file QuEST_cpu.c.

3459 {
3460  // ----- sizes
3461  long long int sizeBlock, // size of blocks
3462  sizeHalfBlock; // size of blocks halved
3463  // ----- indices
3464  long long int thisBlock, // current block
3465  index; // current index for first half block
3466  // ----- measured probability
3467  qreal totalProbability; // probability (returned) value
3468  // ----- temp variables
3469  long long int thisTask;
3470  long long int numTasks=qureg.numAmpsPerChunk>>1;
3471 
3472  // ---------------------------------------------------------------- //
3473  // dimensions //
3474  // ---------------------------------------------------------------- //
3475  sizeHalfBlock = 1LL << (measureQubit); // number of state vector elements to sum,
3476  // and then the number to skip
3477  sizeBlock = 2LL * sizeHalfBlock; // size of blocks (pairs of measure and skip entries)
3478 
3479  // initialise returned value
3480  totalProbability = 0.0;
3481 
3482  qreal *stateVecReal = qureg.stateVec.real;
3483  qreal *stateVecImag = qureg.stateVec.imag;
3484 
3485 # ifdef _OPENMP
3486 # pragma omp parallel \
3487  shared (numTasks,sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag) \
3488  private (thisTask,thisBlock,index) \
3489  reduction ( +:totalProbability )
3490 # endif
3491  {
3492 # ifdef _OPENMP
3493 # pragma omp for schedule (static)
3494 # endif
3495  for (thisTask=0; thisTask<numTasks; thisTask++) {
3496  thisBlock = thisTask / sizeHalfBlock;
3497  index = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
3498 
3499  totalProbability += stateVecReal[index]*stateVecReal[index]
3500  + stateVecImag[index]*stateVecImag[index];
3501  }
3502  }
3503  return totalProbability;
3504 }

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

Referenced by statevec_calcProbOfOutcome().

◆ statevec_hadamardDistributed()

void statevec_hadamardDistributed ( Qureg  qureg,
ComplexArray  stateVecUp,
ComplexArray  stateVecLo,
ComplexArray  stateVecOut,
int  updateUpper 
)

Rotate a single qubit by {{1,1},{1,-1}}/sqrt2.

Operate on a subset of the state vector with upper and lower block values stored seperately. This rotation is just swapping upper and lower values, and stateVecIn must already be the correct section for this chunk

Parameters
[in,out]quregobject representing the set of qubits
[in]stateVecUpprobability amplitudes in upper half of a block depending on chunkId
[in]stateVecLoprobability amplitudes in lower half of a block depending on chunkId
[out]stateVecOutarray section to update (will correspond to either the lower or upper half of a block)
[in]updateUpperflag, 1: updating upper values, 0: updating lower values in block

Definition at line 3139 of file QuEST_cpu.c.

3144 {
3145 
3146  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
3147  long long int thisTask;
3148  long long int numTasks=qureg.numAmpsPerChunk;
3149 
3150  int sign;
3151  if (updateUpper) sign=1;
3152  else sign=-1;
3153 
3154  qreal recRoot2 = 1.0/sqrt(2);
3155 
3156  qreal *stateVecRealUp=stateVecUp.real, *stateVecImagUp=stateVecUp.imag;
3157  qreal *stateVecRealLo=stateVecLo.real, *stateVecImagLo=stateVecLo.imag;
3158  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
3159 
3160 # ifdef _OPENMP
3161 # pragma omp parallel \
3162  default (none) \
3163  shared (stateVecRealUp,stateVecImagUp,stateVecRealLo,stateVecImagLo,stateVecRealOut,stateVecImagOut, \
3164  recRoot2, sign, numTasks) \
3165  private (thisTask,stateRealUp,stateImagUp,stateRealLo,stateImagLo)
3166 # endif
3167  {
3168 # ifdef _OPENMP
3169 # pragma omp for schedule (static)
3170 # endif
3171  for (thisTask=0; thisTask<numTasks; thisTask++) {
3172  // store current state vector values in temp variables
3173  stateRealUp = stateVecRealUp[thisTask];
3174  stateImagUp = stateVecImagUp[thisTask];
3175 
3176  stateRealLo = stateVecRealLo[thisTask];
3177  stateImagLo = stateVecImagLo[thisTask];
3178 
3179  stateVecRealOut[thisTask] = recRoot2*(stateRealUp + sign*stateRealLo);
3180  stateVecImagOut[thisTask] = recRoot2*(stateImagUp + sign*stateImagLo);
3181  }
3182  }
3183 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by statevec_hadamard().

◆ statevec_hadamardLocal()

void statevec_hadamardLocal ( Qureg  qureg,
int  targetQubit 
)

Definition at line 3078 of file QuEST_cpu.c.

3079 {
3080  long long int sizeBlock, sizeHalfBlock;
3081  long long int thisBlock, // current block
3082  indexUp,indexLo; // current index and corresponding index in lower half block
3083 
3084  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
3085  long long int thisTask;
3086  long long int numTasks=qureg.numAmpsPerChunk>>1;
3087 
3088  // set dimensions
3089  sizeHalfBlock = 1LL << targetQubit;
3090  sizeBlock = 2LL * sizeHalfBlock;
3091 
3092  // Can't use qureg.stateVec as a private OMP var
3093  qreal *stateVecReal = qureg.stateVec.real;
3094  qreal *stateVecImag = qureg.stateVec.imag;
3095 
3096  qreal recRoot2 = 1.0/sqrt(2);
3097 
3098 # ifdef _OPENMP
3099 # pragma omp parallel \
3100  default (none) \
3101  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, recRoot2, numTasks) \
3102  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp,stateRealLo,stateImagLo)
3103 # endif
3104  {
3105 # ifdef _OPENMP
3106 # pragma omp for schedule (static)
3107 # endif
3108  for (thisTask=0; thisTask<numTasks; thisTask++) {
3109  thisBlock = thisTask / sizeHalfBlock;
3110  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
3111  indexLo = indexUp + sizeHalfBlock;
3112 
3113  stateRealUp = stateVecReal[indexUp];
3114  stateImagUp = stateVecImag[indexUp];
3115 
3116  stateRealLo = stateVecReal[indexLo];
3117  stateImagLo = stateVecImag[indexLo];
3118 
3119  stateVecReal[indexUp] = recRoot2*(stateRealUp + stateRealLo);
3120  stateVecImag[indexUp] = recRoot2*(stateImagUp + stateImagLo);
3121 
3122  stateVecReal[indexLo] = recRoot2*(stateRealUp - stateRealLo);
3123  stateVecImag[indexLo] = recRoot2*(stateImagUp - stateImagLo);
3124  }
3125  }
3126 }

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

Referenced by statevec_hadamard().

◆ statevec_multiControlledMultiQubitNotDistributed()

void statevec_multiControlledMultiQubitNotDistributed ( Qureg  qureg,
int  ctrlMask,
int  targMask,
ComplexArray  stateVecIn,
ComplexArray  stateVecOut 
)

Definition at line 2834 of file QuEST_cpu.c.

2838  {
2839  long long int numAmps = qureg.numAmpsPerChunk;
2840  long long int globalOffset = qureg.chunkId * numAmps;
2841 
2842  /* stateVecOut is qureg's local state-vector partition, which we modify.
2843  * stateVecIn is the pair node's state-vector partition, in an order which
2844  * does not necessarily correlate to stateVecOut's
2845  */
2846  qreal* inReal = stateVecIn.real;
2847  qreal* inImag = stateVecIn.imag;
2848  qreal* outReal = stateVecOut.real;
2849  qreal* outImag = stateVecOut.imag;
2850 
2851  long long int outInd, outIndGlobal, inInd, inIndGlobal;
2852 
2853 # ifdef _OPENMP
2854 # pragma omp parallel \
2855  default (none) \
2856  shared (inReal,inImag,outReal,outImag, numAmps,globalOffset, ctrlMask,targMask) \
2857  private (outInd,outIndGlobal, inInd,inIndGlobal)
2858 # endif
2859  {
2860 # ifdef _OPENMP
2861 # pragma omp for schedule (static)
2862 # endif
2863  for (outInd = 0; outInd < numAmps; outInd++) {
2864 
2865  // modify amplitude only if control qubits are 1 for this state
2866  outIndGlobal = outInd + globalOffset;
2867  if (ctrlMask && ((ctrlMask & outIndGlobal) != ctrlMask))
2868  continue;
2869  /* it is a premature optimisation to remove this seemingly wasteful abort above,
2870  * because the maximum skipped amplitudes is 1/2 that stored in the node
2871  * (since this function is not called if all amps should be skipped)
2872  */
2873 
2874  /* unlike statevec_controlledNotDistributed(), we cannot assume stateVecOut
2875  * maps contiguously/parallel into stateVecIn; we must map each amplitude, bit-wise.
2876  * However, the arithmetic doesn't necessitate knowing the rank of stateVecIn
2877  */
2878  inIndGlobal = outIndGlobal ^ targMask;
2879  inInd = inIndGlobal % numAmps; // = inIndGlobal - pairRank * numAmps
2880 
2881  outReal[outInd] = inReal[inInd];
2882  outImag[outInd] = inImag[inInd];
2883  }
2884  }
2885 }

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

Referenced by statevec_multiControlledMultiQubitNot().

◆ statevec_multiControlledMultiQubitNotLocal()

void statevec_multiControlledMultiQubitNotLocal ( Qureg  qureg,
int  ctrlMask,
int  targMask 
)

Definition at line 2778 of file QuEST_cpu.c.

2778  {
2779  long long int numAmps = qureg.numAmpsPerChunk;
2780  qreal* stateRe = qureg.stateVec.real;
2781  qreal* stateIm = qureg.stateVec.imag;
2782 
2783  long long int globalOffset = qureg.chunkId * numAmps;
2784 
2785  // each amplitude is swapped with a 'mate' amplitude
2786  long long int ampInd, mateInd, globalInd;
2787  qreal mateRe, mateIm;
2788 
2789 # ifdef _OPENMP
2790 # pragma omp parallel \
2791  default (none) \
2792  shared (stateRe,stateIm, numAmps, ctrlMask,targMask, globalOffset) \
2793  private (ampInd, mateInd,mateRe,mateIm, globalInd)
2794 # endif
2795  {
2796 # ifdef _OPENMP
2797 # pragma omp for schedule (static)
2798 # endif
2799  for (ampInd = 0; ampInd < numAmps; ampInd++) {
2800 
2801  /* it may be a premature optimisation to remove the seemingly wasteful continues below,
2802  * because the maximum skipped amplitudes is 1/2 that stored in the node
2803  * (e.g. since this function is not called if all amps should be skipped via controls),
2804  * and since we're memory-bandwidth bottlenecked.
2805  */
2806 
2807  // although amps are local, we may still be running in distributed mode,
2808  // and hence need to consult the global index to determine the values of
2809  // the control qubits
2810  globalInd = ampInd + globalOffset;
2811 
2812  // modify amplitude only if control qubits are 1 for this state
2813  if (ctrlMask && ((ctrlMask & globalInd) != ctrlMask))
2814  continue;
2815 
2816  mateInd = ampInd ^ targMask;
2817 
2818  // if the mate is behind, it was already processed
2819  if (mateInd < ampInd)
2820  continue;
2821 
2822  mateRe = stateRe[mateInd];
2823  mateIm = stateIm[mateInd];
2824 
2825  // swap amp with mate
2826  stateRe[mateInd] = stateRe[ampInd];
2827  stateIm[mateInd] = stateIm[ampInd];
2828  stateRe[ampInd] = mateRe;
2829  stateIm[ampInd] = mateIm;
2830  }
2831  }
2832 }

References Qureg::chunkId, Qureg::numAmpsPerChunk, qreal, and Qureg::stateVec.

Referenced by statevec_multiControlledMultiQubitNot().

◆ statevec_multiControlledMultiQubitUnitaryLocal()

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

Definition at line 1912 of file QuEST_cpu.c.

1913 {
1914  // can't use qureg.stateVec as a private OMP var
1915  qreal *reVec = qureg.stateVec.real;
1916  qreal *imVec = qureg.stateVec.imag;
1917 
1918  long long int numTasks = qureg.numAmpsPerChunk >> numTargs; // kernel called on every 1 in 2^numTargs amplitudes
1919  long long int numTargAmps = 1 << u.numQubits; // num amps to be modified by each task
1920 
1921  // the global (between all nodes) index of this node's start index
1922  long long int globalIndStart = qureg.chunkId*qureg.numAmpsPerChunk;
1923 
1924  long long int thisTask;
1925  long long int thisInd00; // this thread's index of |..0..0..> (target qubits = 0)
1926  long long int thisGlobalInd00; // the global (between all nodes) index of this thread's |..0..0..> state
1927  long long int ind; // each thread's iteration of amplitudes to modify
1928  int i, t, r, c; // each thread's iteration of amps and targets
1929  qreal reElem, imElem; // each thread's iteration of u elements
1930 
1931  // each thread/task will record and modify numTargAmps amplitudes, privately
1932  // (of course, tasks eliminated by the ctrlMask won't edit their allocation)
1933  //
1934  // If we're NOT on windows, we can fortunately use the stack directly
1935  #ifndef _WIN32
1936  long long int ampInds[numTargAmps];
1937  qreal reAmps[numTargAmps];
1938  qreal imAmps[numTargAmps];
1939 
1940  int sortedTargs[numTargs];
1941  // on Windows, with no VLA, we can use _malloca to allocate on stack (must free)
1942  #else
1943  long long int* ampInds;
1944  qreal* reAmps;
1945  qreal* imAmps;
1946  int* sortedTargs = (int*) _malloca(numTargs * sizeof *sortedTargs);
1947  #endif
1948 
1949  // we need a sorted targets list to find thisInd00 for each task.
1950  // we can't modify targets, because the user-ordering of targets matters in u
1951  for (int t=0; t < numTargs; t++)
1952  sortedTargs[t] = targs[t];
1953  qsort(sortedTargs, numTargs, sizeof(int), qsortComp);
1954 
1955 # ifdef _OPENMP
1956 # pragma omp parallel \
1957  default (none) \
1958  shared (reVec,imVec, numTasks,numTargAmps,globalIndStart, ctrlMask,targs,sortedTargs,u,numTargs) \
1959  private (thisTask,thisInd00,thisGlobalInd00,ind,i,t,r,c,reElem,imElem, ampInds,reAmps,imAmps)
1960 # endif
1961  {
1962  // when manually allocating array memory (on Windows), this must be done in each thread
1963  // separately and is not performed automatically by declaring a var as omp-private
1964  # ifdef _WIN32
1965  ampInds = (long long int*) _malloca(numTargAmps * sizeof *ampInds);
1966  reAmps = (qreal*) _malloca(numTargAmps * sizeof *reAmps);
1967  imAmps = (qreal*) _malloca(numTargAmps * sizeof *imAmps);
1968  # endif
1969 # ifdef _OPENMP
1970 # pragma omp for schedule (static)
1971 # endif
1972  for (thisTask=0; thisTask<numTasks; thisTask++) {
1973 
1974  // find this task's start index (where all targs are 0)
1975  thisInd00 = thisTask;
1976  for (t=0; t < numTargs; t++)
1977  thisInd00 = insertZeroBit(thisInd00, sortedTargs[t]);
1978 
1979  // this task only modifies amplitudes if control qubits are 1 for this state
1980  thisGlobalInd00 = thisInd00 + globalIndStart;
1981  if (ctrlMask && ((ctrlMask & thisGlobalInd00) != ctrlMask))
1982  continue;
1983 
1984  // determine the indices and record values of this tasks's target amps
1985  for (i=0; i < numTargAmps; i++) {
1986 
1987  // get statevec index of current target qubit assignment
1988  ind = thisInd00;
1989  for (t=0; t < numTargs; t++)
1990  if (extractBit(t, i))
1991  ind = flipBit(ind, targs[t]);
1992 
1993  // update this tasks's private arrays
1994  ampInds[i] = ind;
1995  reAmps [i] = reVec[ind];
1996  imAmps [i] = imVec[ind];
1997  }
1998 
1999  // modify this tasks's target amplitudes
2000  for (r=0; r < numTargAmps; r++) {
2001  ind = ampInds[r];
2002  reVec[ind] = 0;
2003  imVec[ind] = 0;
2004 
2005  for (c=0; c < numTargAmps; c++) {
2006  reElem = u.real[r][c];
2007  imElem = u.imag[r][c];
2008  reVec[ind] += reAmps[c]*reElem - imAmps[c]*imElem;
2009  imVec[ind] += reAmps[c]*imElem + imAmps[c]*reElem;
2010  }
2011  }
2012  }
2013  // on Windows, we must explicitly free the stack structures
2014  #ifdef _WIN32
2015  _freea(ampInds);
2016  _freea(reAmps);
2017  _freea(imAmps);
2018  #endif
2019  }
2020 
2021  #ifdef _WIN32
2022  _freea(sortedTargs);
2023  #endif
2024 }

References Qureg::chunkId, extractBit(), flipBit(), ComplexMatrixN::imag, insertZeroBit(), Qureg::numAmpsPerChunk, ComplexMatrixN::numQubits, qreal, qsortComp(), ComplexMatrixN::real, and Qureg::stateVec.

Referenced by statevec_multiControlledMultiQubitUnitary().

◆ statevec_multiControlledTwoQubitUnitaryLocal()

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

Definition at line 1813 of file QuEST_cpu.c.

1813  {
1814 
1815  // can't use qureg.stateVec as a private OMP var
1816  qreal *reVec = qureg.stateVec.real;
1817  qreal *imVec = qureg.stateVec.imag;
1818 
1819  // the global (between all nodes) index of this node's start index
1820  long long int globalIndStart = qureg.chunkId*qureg.numAmpsPerChunk;
1821 
1822  long long int numTasks = qureg.numAmpsPerChunk >> 2; // each iteration updates 4 amplitudes
1823  long long int thisTask;
1824  long long int thisGlobalInd00;
1825  long long int ind00, ind01, ind10, ind11;
1826  qreal re00, re01, re10, re11;
1827  qreal im00, im01, im10, im11;
1828 
1829 # ifdef _OPENMP
1830 # pragma omp parallel \
1831  default (none) \
1832  shared (reVec,imVec,globalIndStart,numTasks,ctrlMask,u,q2,q1) \
1833  private (thisTask, thisGlobalInd00, ind00,ind01,ind10,ind11, re00,re01,re10,re11, im00,im01,im10,im11)
1834 # endif
1835  {
1836 # ifdef _OPENMP
1837 # pragma omp for schedule (static)
1838 # endif
1839  for (thisTask=0; thisTask<numTasks; thisTask++) {
1840 
1841  // determine ind00 of |..0..0..>
1842  ind00 = insertTwoZeroBits(thisTask, q1, q2);
1843 
1844  // skip amplitude if controls aren't in 1 state (overloaded for speed)
1845  thisGlobalInd00 = ind00 + globalIndStart;
1846  if (ctrlMask && ((ctrlMask & thisGlobalInd00) != ctrlMask))
1847  continue;
1848 
1849  // inds of |..0..1..>, |..1..0..> and |..1..1..>
1850  ind01 = flipBit(ind00, q1);
1851  ind10 = flipBit(ind00, q2);
1852  ind11 = flipBit(ind01, q2);
1853 
1854  // extract statevec amplitudes
1855  re00 = reVec[ind00]; im00 = imVec[ind00];
1856  re01 = reVec[ind01]; im01 = imVec[ind01];
1857  re10 = reVec[ind10]; im10 = imVec[ind10];
1858  re11 = reVec[ind11]; im11 = imVec[ind11];
1859 
1860  // apply u * {amp00, amp01, amp10, amp11}
1861  reVec[ind00] =
1862  u.real[0][0]*re00 - u.imag[0][0]*im00 +
1863  u.real[0][1]*re01 - u.imag[0][1]*im01 +
1864  u.real[0][2]*re10 - u.imag[0][2]*im10 +
1865  u.real[0][3]*re11 - u.imag[0][3]*im11;
1866  imVec[ind00] =
1867  u.imag[0][0]*re00 + u.real[0][0]*im00 +
1868  u.imag[0][1]*re01 + u.real[0][1]*im01 +
1869  u.imag[0][2]*re10 + u.real[0][2]*im10 +
1870  u.imag[0][3]*re11 + u.real[0][3]*im11;
1871 
1872  reVec[ind01] =
1873  u.real[1][0]*re00 - u.imag[1][0]*im00 +
1874  u.real[1][1]*re01 - u.imag[1][1]*im01 +
1875  u.real[1][2]*re10 - u.imag[1][2]*im10 +
1876  u.real[1][3]*re11 - u.imag[1][3]*im11;
1877  imVec[ind01] =
1878  u.imag[1][0]*re00 + u.real[1][0]*im00 +
1879  u.imag[1][1]*re01 + u.real[1][1]*im01 +
1880  u.imag[1][2]*re10 + u.real[1][2]*im10 +
1881  u.imag[1][3]*re11 + u.real[1][3]*im11;
1882 
1883  reVec[ind10] =
1884  u.real[2][0]*re00 - u.imag[2][0]*im00 +
1885  u.real[2][1]*re01 - u.imag[2][1]*im01 +
1886  u.real[2][2]*re10 - u.imag[2][2]*im10 +
1887  u.real[2][3]*re11 - u.imag[2][3]*im11;
1888  imVec[ind10] =
1889  u.imag[2][0]*re00 + u.real[2][0]*im00 +
1890  u.imag[2][1]*re01 + u.real[2][1]*im01 +
1891  u.imag[2][2]*re10 + u.real[2][2]*im10 +
1892  u.imag[2][3]*re11 + u.real[2][3]*im11;
1893 
1894  reVec[ind11] =
1895  u.real[3][0]*re00 - u.imag[3][0]*im00 +
1896  u.real[3][1]*re01 - u.imag[3][1]*im01 +
1897  u.real[3][2]*re10 - u.imag[3][2]*im10 +
1898  u.real[3][3]*re11 - u.imag[3][3]*im11;
1899  imVec[ind11] =
1900  u.imag[3][0]*re00 + u.real[3][0]*im00 +
1901  u.imag[3][1]*re01 + u.real[3][1]*im01 +
1902  u.imag[3][2]*re10 + u.real[3][2]*im10 +
1903  u.imag[3][3]*re11 + u.real[3][3]*im11;
1904  }
1905  }
1906 }

References Qureg::chunkId, flipBit(), ComplexMatrix4::imag, insertTwoZeroBits(), Qureg::numAmpsPerChunk, qreal, ComplexMatrix4::real, and Qureg::stateVec.

Referenced by statevec_multiControlledTwoQubitUnitary().

◆ statevec_multiControlledUnitaryDistributed()

void statevec_multiControlledUnitaryDistributed ( Qureg  qureg,
int  targetQubit,
long long int  ctrlQubitsMask,
long long int  ctrlFlipMask,
Complex  rot1,
Complex  rot2,
ComplexArray  stateVecUp,
ComplexArray  stateVecLo,
ComplexArray  stateVecOut 
)

Apply a unitary operation to a single qubit in the state vector of probability amplitudes, given a subset of the state vector with upper and lower block values stored seperately.

Only perform the rotation where all the control qubits are 1.

Parameters
[in,out]quregobject representing the set of qubits
[in]targetQubitqubit to rotate
[in]ctrlQubitsMaska bit mask indicating whether each qubit is a control (1) or not (0)
[in]ctrlFlipMaska bit mask indicating whether each qubit (only controls are relevant) should be flipped when checking the control condition
[in]rot1rotation angle
[in]rot2rotation angle
[in]stateVecUpprobability amplitudes in upper half of a block
[in]stateVecLoprobability amplitudes in lower half of a block
[out]stateVecOutarray section to update (will correspond to either the lower or upper half of a block)

Definition at line 2542 of file QuEST_cpu.c.

2550 {
2551 
2552  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2553  long long int thisTask;
2554  long long int numTasks=qureg.numAmpsPerChunk;
2555  long long int chunkSize=qureg.numAmpsPerChunk;
2556  long long int chunkId=qureg.chunkId;
2557 
2558  qreal rot1Real=rot1.real, rot1Imag=rot1.imag;
2559  qreal rot2Real=rot2.real, rot2Imag=rot2.imag;
2560  qreal *stateVecRealUp=stateVecUp.real, *stateVecImagUp=stateVecUp.imag;
2561  qreal *stateVecRealLo=stateVecLo.real, *stateVecImagLo=stateVecLo.imag;
2562  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2563 
2564 # ifdef _OPENMP
2565 # pragma omp parallel \
2566  default (none) \
2567  shared (stateVecRealUp,stateVecImagUp,stateVecRealLo,stateVecImagLo,stateVecRealOut,stateVecImagOut, \
2568  rot1Real,rot1Imag, rot2Real,rot2Imag, ctrlQubitsMask,ctrlFlipMask, numTasks,chunkId,chunkSize) \
2569  private (thisTask,stateRealUp,stateImagUp,stateRealLo,stateImagLo)
2570 # endif
2571  {
2572 # ifdef _OPENMP
2573 # pragma omp for schedule (static)
2574 # endif
2575  for (thisTask=0; thisTask<numTasks; thisTask++) {
2576  if (ctrlQubitsMask == (ctrlQubitsMask & ((thisTask+chunkId*chunkSize) ^ ctrlFlipMask))) {
2577  // store current state vector values in temp variables
2578  stateRealUp = stateVecRealUp[thisTask];
2579  stateImagUp = stateVecImagUp[thisTask];
2580 
2581  stateRealLo = stateVecRealLo[thisTask];
2582  stateImagLo = stateVecImagLo[thisTask];
2583 
2584  stateVecRealOut[thisTask] = rot1Real*stateRealUp - rot1Imag*stateImagUp
2585  + rot2Real*stateRealLo - rot2Imag*stateImagLo;
2586  stateVecImagOut[thisTask] = rot1Real*stateImagUp + rot1Imag*stateRealUp
2587  + rot2Real*stateImagLo + rot2Imag*stateRealLo;
2588  }
2589  }
2590  }
2591 }

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

Referenced by statevec_multiControlledUnitary().

◆ statevec_multiControlledUnitaryLocal()

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

Definition at line 2268 of file QuEST_cpu.c.

2272 {
2273  long long int sizeBlock, sizeHalfBlock;
2274  long long int thisBlock, // current block
2275  indexUp,indexLo; // current index and corresponding index in lower half block
2276 
2277  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2278  long long int thisTask;
2279  long long int numTasks=qureg.numAmpsPerChunk>>1;
2280  long long int chunkSize=qureg.numAmpsPerChunk;
2281  long long int chunkId=qureg.chunkId;
2282 
2283  // set dimensions
2284  sizeHalfBlock = 1LL << targetQubit;
2285  sizeBlock = 2LL * sizeHalfBlock;
2286 
2287  // Can't use qureg.stateVec as a private OMP var
2288  qreal *stateVecReal = qureg.stateVec.real;
2289  qreal *stateVecImag = qureg.stateVec.imag;
2290 
2291 # ifdef _OPENMP
2292 # pragma omp parallel \
2293  default (none) \
2294  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, u, ctrlQubitsMask,ctrlFlipMask, \
2295  numTasks,chunkId,chunkSize) \
2296  private (thisTask,thisBlock, indexUp,indexLo, stateRealUp,stateImagUp,stateRealLo,stateImagLo)
2297 # endif
2298  {
2299 # ifdef _OPENMP
2300 # pragma omp for schedule (static)
2301 # endif
2302  for (thisTask=0; thisTask<numTasks; thisTask++) {
2303 
2304  thisBlock = thisTask / sizeHalfBlock;
2305  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2306  indexLo = indexUp + sizeHalfBlock;
2307 
2308 
2309  // take the basis index, flip the designated (XOR) 'control' bits, AND with the controls.
2310  // if this equals the control mask, the control qubits have the desired values in the basis index
2311  if (ctrlQubitsMask == (ctrlQubitsMask & ((indexUp+chunkId*chunkSize) ^ ctrlFlipMask))) {
2312  // store current state vector values in temp variables
2313  stateRealUp = stateVecReal[indexUp];
2314  stateImagUp = stateVecImag[indexUp];
2315 
2316  stateRealLo = stateVecReal[indexLo];
2317  stateImagLo = stateVecImag[indexLo];
2318 
2319  // state[indexUp] = u00 * state[indexUp] + u01 * state[indexLo]
2320  stateVecReal[indexUp] = u.real[0][0]*stateRealUp - u.imag[0][0]*stateImagUp
2321  + u.real[0][1]*stateRealLo - u.imag[0][1]*stateImagLo;
2322  stateVecImag[indexUp] = u.real[0][0]*stateImagUp + u.imag[0][0]*stateRealUp
2323  + u.real[0][1]*stateImagLo + u.imag[0][1]*stateRealLo;
2324 
2325  // state[indexLo] = u10 * state[indexUp] + u11 * state[indexLo]
2326  stateVecReal[indexLo] = u.real[1][0]*stateRealUp - u.imag[1][0]*stateImagUp
2327  + u.real[1][1]*stateRealLo - u.imag[1][1]*stateImagLo;
2328  stateVecImag[indexLo] = u.real[1][0]*stateImagUp + u.imag[1][0]*stateRealUp
2329  + u.real[1][1]*stateImagLo + u.imag[1][1]*stateRealLo;
2330  }
2331  }
2332  }
2333 
2334 }

References Qureg::chunkId, ComplexMatrix2::imag, Qureg::numAmpsPerChunk, qreal, ComplexMatrix2::real, and Qureg::stateVec.

Referenced by statevec_multiControlledUnitary().

◆ statevec_pauliXDistributed()

void statevec_pauliXDistributed ( Qureg  qureg,
ComplexArray  stateVecIn,
ComplexArray  stateVecOut 
)

Rotate a single qubit by {{0,1},{1,0}.

Operate on a subset of the state vector with upper and lower block values stored seperately. This rotation is just swapping upper and lower values, and stateVecIn must already be the correct section for this chunk

Remarks
Qubits are zero-based and the
the first qubit is the rightmost
Parameters
[in,out]quregobject representing the set of qubits
[in]stateVecInprobability amplitudes in lower or upper half of a block depending on chunkId
[out]stateVecOutarray section to update (will correspond to either the lower or upper half of a block)

Definition at line 2651 of file QuEST_cpu.c.

2654 {
2655 
2656  long long int thisTask;
2657  long long int numTasks=qureg.numAmpsPerChunk;
2658 
2659  qreal *stateVecRealIn=stateVecIn.real, *stateVecImagIn=stateVecIn.imag;
2660  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2661 
2662 # ifdef _OPENMP
2663 # pragma omp parallel \
2664  default (none) \
2665  shared (stateVecRealIn,stateVecImagIn,stateVecRealOut,stateVecImagOut,numTasks) \
2666  private (thisTask)
2667 # endif
2668  {
2669 # ifdef _OPENMP
2670 # pragma omp for schedule (static)
2671 # endif
2672  for (thisTask=0; thisTask<numTasks; thisTask++) {
2673  stateVecRealOut[thisTask] = stateVecRealIn[thisTask];
2674  stateVecImagOut[thisTask] = stateVecImagIn[thisTask];
2675  }
2676  }
2677 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by statevec_pauliX().

◆ statevec_pauliXLocal()

void statevec_pauliXLocal ( Qureg  qureg,
int  targetQubit 
)

Definition at line 2593 of file QuEST_cpu.c.

2594 {
2595  long long int sizeBlock, sizeHalfBlock;
2596  long long int thisBlock, // current block
2597  indexUp,indexLo; // current index and corresponding index in lower half block
2598 
2599  qreal stateRealUp,stateImagUp;
2600  long long int thisTask;
2601  long long int numTasks=qureg.numAmpsPerChunk>>1;
2602 
2603  // set dimensions
2604  sizeHalfBlock = 1LL << targetQubit;
2605  sizeBlock = 2LL * sizeHalfBlock;
2606 
2607  // Can't use qureg.stateVec as a private OMP var
2608  qreal *stateVecReal = qureg.stateVec.real;
2609  qreal *stateVecImag = qureg.stateVec.imag;
2610 
2611 # ifdef _OPENMP
2612 # pragma omp parallel \
2613  default (none) \
2614  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, numTasks) \
2615  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp)
2616 # endif
2617  {
2618 # ifdef _OPENMP
2619 # pragma omp for schedule (static)
2620 # endif
2621  for (thisTask=0; thisTask<numTasks; thisTask++) {
2622  thisBlock = thisTask / sizeHalfBlock;
2623  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2624  indexLo = indexUp + sizeHalfBlock;
2625 
2626  stateRealUp = stateVecReal[indexUp];
2627  stateImagUp = stateVecImag[indexUp];
2628 
2629  stateVecReal[indexUp] = stateVecReal[indexLo];
2630  stateVecImag[indexUp] = stateVecImag[indexLo];
2631 
2632  stateVecReal[indexLo] = stateRealUp;
2633  stateVecImag[indexLo] = stateImagUp;
2634  }
2635  }
2636 
2637 }

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

Referenced by statevec_pauliX().

◆ statevec_pauliYDistributed()

void statevec_pauliYDistributed ( Qureg  qureg,
ComplexArray  stateVecIn,
ComplexArray  stateVecOut,
int  updateUpper,
int  conjFac 
)

Rotate a single qubit by +-{{0,-i},{i,0}.

Operate on a subset of the state vector with upper and lower block values stored seperately. This rotation is just swapping upper and lower values, and stateVecIn must already be the correct section for this chunk

Remarks
Qubits are zero-based and the
the first qubit is the rightmost
Parameters
[in,out]quregobject representing the set of qubits
[in]stateVecInprobability amplitudes in lower or upper half of a block depending on chunkId
[out]stateVecOutarray section to update (will correspond to either the lower or upper half of a block)
[in]updateUpperflag, 1: updating upper values, 0: updating lower values in block
[in]conjFac1: apply conj(pauliY), 0: apply pauliY

Definition at line 2945 of file QuEST_cpu.c.

2949 {
2950 
2951  long long int thisTask;
2952  long long int numTasks=qureg.numAmpsPerChunk;
2953 
2954  qreal *stateVecRealIn=stateVecIn.real, *stateVecImagIn=stateVecIn.imag;
2955  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2956 
2957  int realSign=1, imagSign=1;
2958  if (updateUpper) imagSign=-1;
2959  else realSign = -1;
2960 
2961 # ifdef _OPENMP
2962 # pragma omp parallel \
2963  default (none) \
2964  shared (stateVecRealIn,stateVecImagIn,stateVecRealOut,stateVecImagOut, \
2965  realSign,imagSign, numTasks,conjFac) \
2966  private (thisTask)
2967 # endif
2968  {
2969 # ifdef _OPENMP
2970 # pragma omp for schedule (static)
2971 # endif
2972  for (thisTask=0; thisTask<numTasks; thisTask++) {
2973  stateVecRealOut[thisTask] = conjFac * realSign * stateVecImagIn[thisTask];
2974  stateVecImagOut[thisTask] = conjFac * imagSign * stateVecRealIn[thisTask];
2975  }
2976  }
2977 }

References Qureg::numAmpsPerChunk, and qreal.

Referenced by statevec_pauliY(), and statevec_pauliYConj().

◆ statevec_pauliYLocal()

void statevec_pauliYLocal ( Qureg  qureg,
int  targetQubit,
int  conjFac 
)

Definition at line 2887 of file QuEST_cpu.c.

2888 {
2889  long long int sizeBlock, sizeHalfBlock;
2890  long long int thisBlock, // current block
2891  indexUp,indexLo; // current index and corresponding index in lower half block
2892 
2893  qreal stateRealUp,stateImagUp;
2894  long long int thisTask;
2895  long long int numTasks=qureg.numAmpsPerChunk>>1;
2896 
2897  // set dimensions
2898  sizeHalfBlock = 1LL << targetQubit;
2899  sizeBlock = 2LL * sizeHalfBlock;
2900 
2901  // Can't use qureg.stateVec as a private OMP var
2902  qreal *stateVecReal = qureg.stateVec.real;
2903  qreal *stateVecImag = qureg.stateVec.imag;
2904 
2905 # ifdef _OPENMP
2906 # pragma omp parallel \
2907  default (none) \
2908  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, numTasks,conjFac) \
2909  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp)
2910 # endif
2911  {
2912 # ifdef _OPENMP
2913 # pragma omp for schedule (static)
2914 # endif
2915  for (thisTask=0; thisTask<numTasks; thisTask++) {
2916  thisBlock = thisTask / sizeHalfBlock;
2917  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2918  indexLo = indexUp + sizeHalfBlock;
2919 
2920  stateRealUp = stateVecReal[indexUp];
2921  stateImagUp = stateVecImag[indexUp];
2922 
2923  stateVecReal[indexUp] = conjFac * stateVecImag[indexLo];
2924  stateVecImag[indexUp] = conjFac * -stateVecReal[indexLo];
2925  stateVecReal[indexLo] = conjFac * -stateImagUp;
2926  stateVecImag[indexLo] = conjFac * stateRealUp;
2927  }
2928  }
2929 }

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

Referenced by statevec_pauliY(), and statevec_pauliYConj().

◆ statevec_swapQubitAmpsDistributed()

void statevec_swapQubitAmpsDistributed ( Qureg  qureg,
int  pairRank,
int  qb1,
int  qb2 
)

qureg.pairStateVec contains the entire set of amplitudes of the paired node which includes the set of all amplitudes which need to be swapped between |..0..1..> and |..1..0..>

Definition at line 3965 of file QuEST_cpu.c.

3965  {
3966 
3967  // can't use qureg.stateVec as a private OMP var
3968  qreal *reVec = qureg.stateVec.real;
3969  qreal *imVec = qureg.stateVec.imag;
3970  qreal *rePairVec = qureg.pairStateVec.real;
3971  qreal *imPairVec = qureg.pairStateVec.imag;
3972 
3973  long long int numLocalAmps = qureg.numAmpsPerChunk;
3974  long long int globalStartInd = qureg.chunkId * numLocalAmps;
3975  long long int pairGlobalStartInd = pairRank * numLocalAmps;
3976 
3977  long long int localInd, globalInd;
3978  long long int pairLocalInd, pairGlobalInd;
3979 
3980 # ifdef _OPENMP
3981 # pragma omp parallel \
3982  default (none) \
3983  shared (reVec,imVec,rePairVec,imPairVec,numLocalAmps,globalStartInd,pairGlobalStartInd,qb1,qb2) \
3984  private (localInd,globalInd, pairLocalInd,pairGlobalInd)
3985 # endif
3986  {
3987 # ifdef _OPENMP
3988 # pragma omp for schedule (static)
3989 # endif
3990  for (localInd=0; localInd < numLocalAmps; localInd++) {
3991 
3992  globalInd = globalStartInd + localInd;
3993  if (isOddParity(globalInd, qb1, qb2)) {
3994 
3995  pairGlobalInd = flipBit(flipBit(globalInd, qb1), qb2);
3996  pairLocalInd = pairGlobalInd - pairGlobalStartInd;
3997 
3998  reVec[localInd] = rePairVec[pairLocalInd];
3999  imVec[localInd] = imPairVec[pairLocalInd];
4000  }
4001  }
4002  }
4003 }

References Qureg::chunkId, flipBit(), isOddParity(), Qureg::numAmpsPerChunk, Qureg::pairStateVec, qreal, and Qureg::stateVec.

Referenced by statevec_swapQubitAmps().

◆ statevec_swapQubitAmpsLocal()

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

It is ensured that all amplitudes needing to be swapped are on this node.

This means that amplitudes for |a 0..0..> to |a 1..1..> all exist on this node and each node has a different bit-string prefix "a". The prefix 'a' (and ergo, the chunkID) don't enter the calculations for the offset of |a 0..1..> and |a 1..0..> from |a 0..0..> and ergo are not featured below.

Definition at line 3922 of file QuEST_cpu.c.

3922  {
3923 
3924  // can't use qureg.stateVec as a private OMP var
3925  qreal *reVec = qureg.stateVec.real;
3926  qreal *imVec = qureg.stateVec.imag;
3927 
3928  long long int numTasks = qureg.numAmpsPerChunk >> 2; // each iteration updates 2 amps and skips 2 amps
3929  long long int thisTask;
3930  long long int ind00, ind01, ind10;
3931  qreal re01, re10;
3932  qreal im01, im10;
3933 
3934 # ifdef _OPENMP
3935 # pragma omp parallel \
3936  default (none) \
3937  shared (reVec,imVec,numTasks,qb1,qb2) \
3938  private (thisTask, ind00,ind01,ind10, re01,re10, im01,im10)
3939 # endif
3940  {
3941 # ifdef _OPENMP
3942 # pragma omp for schedule (static)
3943 # endif
3944  for (thisTask=0; thisTask<numTasks; thisTask++) {
3945  // determine ind00 of |..0..0..>, |..0..1..> and |..1..0..>
3946  ind00 = insertTwoZeroBits(thisTask, qb1, qb2);
3947  ind01 = flipBit(ind00, qb1);
3948  ind10 = flipBit(ind00, qb2);
3949 
3950  // extract statevec amplitudes
3951  re01 = reVec[ind01]; im01 = imVec[ind01];
3952  re10 = reVec[ind10]; im10 = imVec[ind10];
3953 
3954  // swap 01 and 10 amps
3955  reVec[ind01] = re10; reVec[ind10] = re01;
3956  imVec[ind01] = im10; imVec[ind10] = im01;
3957  }
3958  }
3959 }

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

Referenced by statevec_swapQubitAmps().

◆ statevec_unitaryDistributed()

void statevec_unitaryDistributed ( Qureg  qureg,
Complex  rot1,
Complex  rot2,
ComplexArray  stateVecUp,
ComplexArray  stateVecLo,
ComplexArray  stateVecOut 
)

Apply a unitary operation to a single qubit given a subset of the state vector with upper and lower block values stored seperately.

Remarks
Qubits are zero-based and the first qubit is the rightmost
Parameters
[in,out]quregobject representing the set of qubits
[in]rot1
[in]rot2
[in]stateVecUpprobability amplitudes in upper half of a block
[in]stateVecLoprobability amplitudes in lower half of a block
[out]stateVecOutarray section to update (will correspond to either the lower or upper half of a block)

Definition at line 2151 of file QuEST_cpu.c.

2156 {
2157 
2158  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2159  long long int thisTask;
2160  long long int numTasks=qureg.numAmpsPerChunk;
2161 
2162  qreal rot1Real=rot1.real, rot1Imag=rot1.imag;
2163  qreal rot2Real=rot2.real, rot2Imag=rot2.imag;
2164  qreal *stateVecRealUp=stateVecUp.real, *stateVecImagUp=stateVecUp.imag;
2165  qreal *stateVecRealLo=stateVecLo.real, *stateVecImagLo=stateVecLo.imag;
2166  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2167 
2168 
2169 # ifdef _OPENMP
2170 # pragma omp parallel \
2171  default (none) \
2172  shared (stateVecRealUp,stateVecImagUp,stateVecRealLo,stateVecImagLo,stateVecRealOut,stateVecImagOut, \
2173  rot1Real, rot1Imag, rot2Real, rot2Imag,numTasks) \
2174  private (thisTask,stateRealUp,stateImagUp,stateRealLo,stateImagLo)
2175 # endif
2176  {
2177 # ifdef _OPENMP
2178 # pragma omp for schedule (static)
2179 # endif
2180  for (thisTask=0; thisTask<numTasks; thisTask++) {
2181  // store current state vector values in temp variables
2182  stateRealUp = stateVecRealUp[thisTask];
2183  stateImagUp = stateVecImagUp[thisTask];
2184 
2185  stateRealLo = stateVecRealLo[thisTask];
2186  stateImagLo = stateVecImagLo[thisTask];
2187 
2188  stateVecRealOut[thisTask] = rot1Real*stateRealUp - rot1Imag*stateImagUp
2189  + rot2Real*stateRealLo - rot2Imag*stateImagLo;
2190  stateVecImagOut[thisTask] = rot1Real*stateImagUp + rot1Imag*stateRealUp
2191  + rot2Real*stateImagLo + rot2Imag*stateRealLo;
2192  }
2193  }
2194 }

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

Referenced by statevec_unitary().

◆ statevec_unitaryLocal()

void statevec_unitaryLocal ( Qureg  qureg,
int  targetQubit,
ComplexMatrix2  u 
)

Definition at line 2026 of file QuEST_cpu.c.

2027 {
2028  long long int sizeBlock, sizeHalfBlock;
2029  long long int thisBlock, // current block
2030  indexUp,indexLo; // current index and corresponding index in lower half block
2031 
2032  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2033  long long int thisTask;
2034  long long int numTasks=qureg.numAmpsPerChunk>>1;
2035 
2036  // set dimensions
2037  sizeHalfBlock = 1LL << targetQubit;
2038  sizeBlock = 2LL * sizeHalfBlock;
2039 
2040  // Can't use qureg.stateVec as a private OMP var
2041  qreal *stateVecReal = qureg.stateVec.real;
2042  qreal *stateVecImag = qureg.stateVec.imag;
2043 
2044 # ifdef _OPENMP
2045 # pragma omp parallel \
2046  default (none) \
2047  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, u,numTasks) \
2048  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp,stateRealLo,stateImagLo)
2049 # endif
2050  {
2051 # ifdef _OPENMP
2052 # pragma omp for schedule (static)
2053 # endif
2054  for (thisTask=0; thisTask<numTasks; thisTask++) {
2055 
2056  thisBlock = thisTask / sizeHalfBlock;
2057  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2058  indexLo = indexUp + sizeHalfBlock;
2059 
2060  // store current state vector values in temp variables
2061  stateRealUp = stateVecReal[indexUp];
2062  stateImagUp = stateVecImag[indexUp];
2063 
2064  stateRealLo = stateVecReal[indexLo];
2065  stateImagLo = stateVecImag[indexLo];
2066 
2067 
2068  // state[indexUp] = u00 * state[indexUp] + u01 * state[indexLo]
2069  stateVecReal[indexUp] = u.real[0][0]*stateRealUp - u.imag[0][0]*stateImagUp
2070  + u.real[0][1]*stateRealLo - u.imag[0][1]*stateImagLo;
2071  stateVecImag[indexUp] = u.real[0][0]*stateImagUp + u.imag[0][0]*stateRealUp
2072  + u.real[0][1]*stateImagLo + u.imag[0][1]*stateRealLo;
2073 
2074  // state[indexLo] = u10 * state[indexUp] + u11 * state[indexLo]
2075  stateVecReal[indexLo] = u.real[1][0]*stateRealUp - u.imag[1][0]*stateImagUp
2076  + u.real[1][1]*stateRealLo - u.imag[1][1]*stateImagLo;
2077  stateVecImag[indexLo] = u.real[1][0]*stateImagUp + u.imag[1][0]*stateRealUp
2078  + u.real[1][1]*stateImagLo + u.imag[1][1]*stateRealLo;
2079 
2080  }
2081  }
2082 }

References ComplexMatrix2::imag, Qureg::numAmpsPerChunk, qreal, ComplexMatrix2::real, and Qureg::stateVec.

Referenced by statevec_unitary().

qsortComp
int qsortComp(const void *a, const void *b)
Definition: QuEST_cpu.c:1908
ComplexMatrix4::real
qreal real[4][4]
Definition: QuEST.h:177
Qureg::pairStateVec
ComplexArray pairStateVec
Temporary storage for a chunk of the state vector received from another process in the MPI version.
Definition: QuEST.h:343
densmatr_mixDephasing
void densmatr_mixDephasing(Qureg qureg, int targetQubit, qreal dephase)
Definition: QuEST_cpu.c:85
qreal
#define qreal
flipBit
__forceinline__ __device__ long long int flipBit(const long long int number, const int bitInd)
Definition: QuEST_gpu.cu:95
extractBit
static int extractBit(const int locationOfBitFromRight, const long long int theEncodedNumber)
Definition: QuEST_cpu_internal.h:26
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:336
ComplexMatrix2::imag
qreal imag[2][2]
Definition: QuEST.h:140
insertZeroBit
__forceinline__ __device__ long long int insertZeroBit(const long long int number, const int index)
Definition: QuEST_gpu.cu:99
DiagonalOp::imag
qreal * imag
The imaginary values of the 2^numQubits complex elements.
Definition: QuEST.h:310
Qureg::numAmpsPerChunk
long long int numAmpsPerChunk
Number of probability amplitudes held in stateVec by this process In the non-MPI version,...
Definition: QuEST.h:332
densmatr_oneQubitDegradeOffDiagonal
void densmatr_oneQubitDegradeOffDiagonal(Qureg qureg, int targetQubit, qreal retain)
Definition: QuEST_cpu.c:54
ComplexMatrix4::imag
qreal imag[4][4]
Definition: QuEST.h:178
DiagonalOp::numQubits
int numQubits
The number of qubits this operator can act on (informing its size)
Definition: QuEST.h:300
ComplexMatrixN::real
qreal ** real
Definition: QuEST.h:189
extractBit
__forceinline__ __device__ int extractBit(const int locationOfBitFromRight, const long long int theEncodedNumber)
Definition: QuEST_gpu.cu:82
insertZeroBit
static long long int insertZeroBit(const long long int number, const int index)
Definition: QuEST_cpu_internal.h:42
insertTwoZeroBits
__forceinline__ __device__ long long int insertTwoZeroBits(const long long int number, const int bit1, const int bit2)
Definition: QuEST_gpu.cu:106
ComplexMatrixN::imag
qreal ** imag
Definition: QuEST.h:190
Qureg::stateVec
ComplexArray stateVec
Computational state amplitudes - a subset thereof in the MPI version.
Definition: QuEST.h:341
ComplexMatrix2::real
qreal real[2][2]
Definition: QuEST.h:139
ComplexMatrixN::numQubits
int numQubits
Definition: QuEST.h:188
isOddParity
static int isOddParity(const long long int number, const int qb1, const int qb2)
Definition: QuEST_cpu_internal.h:38
Qureg::numQubitsRepresented
int numQubitsRepresented
The number of qubits represented in either the state-vector or density matrix.
Definition: QuEST.h:327
Qureg::numAmpsTotal
long long int numAmpsTotal
Total number of amplitudes, which are possibly distributed among machines.
Definition: QuEST.h:334
DiagonalOp::real
qreal * real
The real values of the 2^numQubits complex elements.
Definition: QuEST.h:308
Complex::real
qreal real
Definition: QuEST.h:105
Complex::imag
qreal imag
Definition: QuEST.h:106
Complex
Represents one complex number.
Definition: QuEST.h:103