QuEST_cpu.c
Go to the documentation of this file.
1 // Distributed under MIT licence. See https://github.com/QuEST-Kit/QuEST/blob/master/LICENCE.txt for details
2 
13 # include "QuEST.h"
14 # include "QuEST_internal.h"
15 # include "QuEST_precision.h"
16 # include "mt19937ar.h"
17 
18 # include "QuEST_cpu_internal.h"
19 
20 # include <math.h>
21 # include <stdio.h>
22 # include <stdlib.h>
23 # include <stdint.h>
24 # include <assert.h>
25 
26 # ifdef _OPENMP
27 # include <omp.h>
28 # endif
29 
30 
31 
32 /*
33  * overloads for consistent API with GPU
34  */
35 
36 void copyStateToGPU(Qureg qureg) {
37 }
38 
39 void copyStateFromGPU(Qureg qureg) {
40 }
41 
42 
43 
44 /*
45  * state vector and density matrix operations
46  */
47 
48 void densmatr_oneQubitDegradeOffDiagonal(Qureg qureg, const int targetQubit, qreal retain){
49  const long long int numTasks = qureg.numAmpsPerChunk;
50  long long int innerMask = 1LL << targetQubit;
51  long long int outerMask = 1LL << (targetQubit + (qureg.numQubitsRepresented));
52 
53  long long int thisTask;
54  long long int thisPattern;
55  long long int totMask = innerMask|outerMask;
56 
57  # ifdef _OPENMP
58 # pragma omp parallel \
59  default (none) \
60  shared (innerMask,outerMask,totMask,qureg,retain) \
61  private (thisTask,thisPattern)
62 # endif
63  {
64 # ifdef _OPENMP
65 # pragma omp for schedule (static)
66 # endif
67  for (thisTask=0; thisTask<numTasks; thisTask++){
68  thisPattern = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMask;
69  if ((thisPattern==innerMask) || (thisPattern==outerMask)){
70  // do dephase
71  // the lines below will degrade the off-diagonal terms |..0..><..1..| and |..1..><..0..|
72  qureg.stateVec.real[thisTask] = retain*qureg.stateVec.real[thisTask];
73  qureg.stateVec.imag[thisTask] = retain*qureg.stateVec.imag[thisTask];
74  }
75  }
76  }
77 }
78 
79 void densmatr_mixDephasing(Qureg qureg, const int targetQubit, qreal dephase) {
80  qreal retain=1-dephase;
81  densmatr_oneQubitDegradeOffDiagonal(qureg, targetQubit, retain);
82 }
83 
84 void densmatr_mixTwoQubitDephasing(Qureg qureg, const int qubit1, const int qubit2, qreal dephase) {
85  qreal retain=1-dephase;
86 
87  const long long int numTasks = qureg.numAmpsPerChunk;
88  long long int innerMaskQubit1 = 1LL << qubit1;
89  long long int outerMaskQubit1 = 1LL << (qubit1 + (qureg.numQubitsRepresented));
90  long long int innerMaskQubit2 = 1LL << qubit2;
91  long long int outerMaskQubit2 = 1LL << (qubit2 + (qureg.numQubitsRepresented));
92  long long int totMaskQubit1 = innerMaskQubit1|outerMaskQubit1;
93  long long int totMaskQubit2 = innerMaskQubit2|outerMaskQubit2;
94 
95  long long int thisTask;
96  long long int thisPatternQubit1, thisPatternQubit2;
97 
98 # ifdef _OPENMP
99 # pragma omp parallel \
100  default (none) \
101  shared (innerMaskQubit1,outerMaskQubit1,totMaskQubit1,innerMaskQubit2,outerMaskQubit2, \
102  totMaskQubit2,qureg,retain) \
103  private (thisTask,thisPatternQubit1,thisPatternQubit2)
104 # endif
105  {
106 # ifdef _OPENMP
107 # pragma omp for schedule (static)
108 # endif
109  for (thisTask=0; thisTask<numTasks; thisTask++){
110  thisPatternQubit1 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit1;
111  thisPatternQubit2 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit2;
112 
113  // any mismatch |...0...><...1...| etc
114  if ( (thisPatternQubit1==innerMaskQubit1) || (thisPatternQubit1==outerMaskQubit1) ||
115  (thisPatternQubit2==innerMaskQubit2) || (thisPatternQubit2==outerMaskQubit2) ){
116  // do dephase
117  // the lines below will degrade the off-diagonal terms |..0..><..1..| and |..1..><..0..|
118  qureg.stateVec.real[thisTask] = retain*qureg.stateVec.real[thisTask];
119  qureg.stateVec.imag[thisTask] = retain*qureg.stateVec.imag[thisTask];
120  }
121  }
122  }
123 }
124 
125 void densmatr_mixDepolarisingLocal(Qureg qureg, const int targetQubit, qreal depolLevel) {
126  qreal retain=1-depolLevel;
127 
128  const long long int numTasks = qureg.numAmpsPerChunk;
129  long long int innerMask = 1LL << targetQubit;
130  long long int outerMask = 1LL << (targetQubit + (qureg.numQubitsRepresented));
131  long long int totMask = innerMask|outerMask;
132 
133  long long int thisTask;
134  long long int partner;
135  long long int thisPattern;
136 
137  qreal realAv, imagAv;
138 
139 # ifdef _OPENMP
140 # pragma omp parallel \
141  default (none) \
142  shared (innerMask,outerMask,totMask,qureg,retain,depolLevel) \
143  private (thisTask,partner,thisPattern,realAv,imagAv)
144 # endif
145  {
146 # ifdef _OPENMP
147 # pragma omp for schedule (static)
148 # endif
149  for (thisTask=0; thisTask<numTasks; thisTask++){
150  thisPattern = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMask;
151  if ((thisPattern==innerMask) || (thisPattern==outerMask)){
152  // do dephase
153  // the lines below will degrade the off-diagonal terms |..0..><..1..| and |..1..><..0..|
154  qureg.stateVec.real[thisTask] = retain*qureg.stateVec.real[thisTask];
155  qureg.stateVec.imag[thisTask] = retain*qureg.stateVec.imag[thisTask];
156  } else {
157  if ((thisTask&totMask)==0){ //this element relates to targetQubit in state 0
158  // do depolarise
159  partner = thisTask | totMask;
160  realAv = (qureg.stateVec.real[thisTask] + qureg.stateVec.real[partner]) /2 ;
161  imagAv = (qureg.stateVec.imag[thisTask] + qureg.stateVec.imag[partner]) /2 ;
162 
163  qureg.stateVec.real[thisTask] = retain*qureg.stateVec.real[thisTask] + depolLevel*realAv;
164  qureg.stateVec.imag[thisTask] = retain*qureg.stateVec.imag[thisTask] + depolLevel*imagAv;
165 
166  qureg.stateVec.real[partner] = retain*qureg.stateVec.real[partner] + depolLevel*realAv;
167  qureg.stateVec.imag[partner] = retain*qureg.stateVec.imag[partner] + depolLevel*imagAv;
168  }
169  }
170  }
171  }
172 }
173 
174 void densmatr_mixDampingLocal(Qureg qureg, const int targetQubit, qreal damping) {
175  qreal retain=1-damping;
176  qreal dephase=sqrt(retain);
177 
178  const long long int numTasks = qureg.numAmpsPerChunk;
179  long long int innerMask = 1LL << targetQubit;
180  long long int outerMask = 1LL << (targetQubit + (qureg.numQubitsRepresented));
181  long long int totMask = innerMask|outerMask;
182 
183  long long int thisTask;
184  long long int partner;
185  long long int thisPattern;
186 
187  //qreal realAv, imagAv;
188 
189 # ifdef _OPENMP
190 # pragma omp parallel \
191  default (none) \
192  shared (innerMask,outerMask,totMask,qureg,retain,damping,dephase) \
193  private (thisTask,partner,thisPattern)
194 # endif
195  {
196 # ifdef _OPENMP
197 # pragma omp for schedule (static)
198 # endif
199  for (thisTask=0; thisTask<numTasks; thisTask++){
200  thisPattern = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMask;
201  if ((thisPattern==innerMask) || (thisPattern==outerMask)){
202  // do dephase
203  // the lines below will degrade the off-diagonal terms |..0..><..1..| and |..1..><..0..|
204  qureg.stateVec.real[thisTask] = dephase*qureg.stateVec.real[thisTask];
205  qureg.stateVec.imag[thisTask] = dephase*qureg.stateVec.imag[thisTask];
206  } else {
207  if ((thisTask&totMask)==0){ //this element relates to targetQubit in state 0
208  // do depolarise
209  partner = thisTask | totMask;
210  //realAv = (qureg.stateVec.real[thisTask] + qureg.stateVec.real[partner]) /2 ;
211  //imagAv = (qureg.stateVec.imag[thisTask] + qureg.stateVec.imag[partner]) /2 ;
212 
213  qureg.stateVec.real[thisTask] = qureg.stateVec.real[thisTask] + damping*qureg.stateVec.real[partner];
214  qureg.stateVec.imag[thisTask] = qureg.stateVec.imag[thisTask] + damping*qureg.stateVec.imag[partner];
215 
216  qureg.stateVec.real[partner] = retain*qureg.stateVec.real[partner];
217  qureg.stateVec.imag[partner] = retain*qureg.stateVec.imag[partner];
218  }
219  }
220  }
221  }
222 }
223 
224 void densmatr_mixDepolarisingDistributed(Qureg qureg, const int targetQubit, qreal depolLevel) {
225 
226  // first do dephase part.
227  // TODO -- this might be more efficient to do at the same time as the depolarise if we move to
228  // iterating over all elements in the state vector for the purpose of vectorisation
229  // TODO -- if we keep this split, move this function to densmatr_mixDepolarising()
230  densmatr_mixDephasing(qureg, targetQubit, depolLevel);
231 
232  long long int sizeInnerBlock, sizeInnerHalfBlock;
233  long long int sizeOuterColumn, sizeOuterHalfColumn;
234  long long int thisInnerBlock, // current block
235  thisOuterColumn, // current column in density matrix
236  thisIndex, // current index in (density matrix representation) state vector
237  thisIndexInOuterColumn,
238  thisIndexInInnerBlock;
239  int outerBit;
240 
241  long long int thisTask;
242  const long long int numTasks=qureg.numAmpsPerChunk>>1;
243 
244  // set dimensions
245  sizeInnerHalfBlock = 1LL << targetQubit;
246  sizeInnerBlock = 2LL * sizeInnerHalfBlock;
247  sizeOuterColumn = 1LL << qureg.numQubitsRepresented;
248  sizeOuterHalfColumn = sizeOuterColumn >> 1;
249 
250 # ifdef _OPENMP
251 # pragma omp parallel \
252  default (none) \
253  shared (sizeInnerBlock,sizeInnerHalfBlock,sizeOuterColumn,sizeOuterHalfColumn,qureg,depolLevel) \
254  private (thisTask,thisInnerBlock,thisOuterColumn,thisIndex,thisIndexInOuterColumn, \
255  thisIndexInInnerBlock,outerBit)
256 # endif
257  {
258 # ifdef _OPENMP
259 # pragma omp for schedule (static)
260 # endif
261  // thisTask iterates over half the elements in this process' chunk of the density matrix
262  // treat this as iterating over all columns, then iterating over half the values
263  // within one column.
264  // If this function has been called, this process' chunk contains half an
265  // outer block or less
266  for (thisTask=0; thisTask<numTasks; thisTask++) {
267  // we want to process all columns in the density matrix,
268  // updating the values for half of each column (one half of each inner block)
269  thisOuterColumn = thisTask / sizeOuterHalfColumn;
270  thisIndexInOuterColumn = thisTask&(sizeOuterHalfColumn-1); // thisTask % sizeOuterHalfColumn
271  thisInnerBlock = thisIndexInOuterColumn/sizeInnerHalfBlock;
272  // get index in state vector corresponding to upper inner block
273  thisIndexInInnerBlock = thisTask&(sizeInnerHalfBlock-1); // thisTask % sizeInnerHalfBlock
274  thisIndex = thisOuterColumn*sizeOuterColumn + thisInnerBlock*sizeInnerBlock
275  + thisIndexInInnerBlock;
276  // check if we are in the upper or lower half of an outer block
277  outerBit = extractBit(targetQubit, (thisIndex+qureg.numAmpsPerChunk*qureg.chunkId)>>qureg.numQubitsRepresented);
278  // if we are in the lower half of an outer block, shift to be in the lower half
279  // of the inner block as well (we want to dephase |0><0| and |1><1| only)
280  thisIndex += outerBit*(sizeInnerHalfBlock);
281 
282  // NOTE: at this point thisIndex should be the index of the element we want to
283  // dephase in the chunk of the state vector on this process, in the
284  // density matrix representation.
285  // thisTask is the index of the pair element in pairStateVec
286 
287 
288  // state[thisIndex] = (1-depolLevel)*state[thisIndex] + depolLevel*(state[thisIndex]
289  // + pair[thisTask])/2
290  qureg.stateVec.real[thisIndex] = (1-depolLevel)*qureg.stateVec.real[thisIndex] +
291  depolLevel*(qureg.stateVec.real[thisIndex] + qureg.pairStateVec.real[thisTask])/2;
292 
293  qureg.stateVec.imag[thisIndex] = (1-depolLevel)*qureg.stateVec.imag[thisIndex] +
294  depolLevel*(qureg.stateVec.imag[thisIndex] + qureg.pairStateVec.imag[thisTask])/2;
295  }
296  }
297 }
298 
299 void densmatr_mixDampingDistributed(Qureg qureg, const int targetQubit, qreal damping) {
300  qreal retain=1-damping;
301  qreal dephase=sqrt(1-damping);
302 
303  // multiply the off-diagonal (|0><1| and |1><0|) terms by sqrt(1-damping)
304  densmatr_oneQubitDegradeOffDiagonal(qureg, targetQubit, dephase);
305 
306  // below, we modify the diagonals terms which require |1><1| to |0><0| communication
307 
308  long long int sizeInnerBlock, sizeInnerHalfBlock;
309  long long int sizeOuterColumn, sizeOuterHalfColumn;
310  long long int thisInnerBlock, // current block
311  thisOuterColumn, // current column in density matrix
312  thisIndex, // current index in (density matrix representation) state vector
313  thisIndexInOuterColumn,
314  thisIndexInInnerBlock;
315  int outerBit;
316  int stateBit;
317 
318  long long int thisTask;
319  const long long int numTasks=qureg.numAmpsPerChunk>>1;
320 
321  // set dimensions
322  sizeInnerHalfBlock = 1LL << targetQubit;
323  sizeInnerBlock = 2LL * sizeInnerHalfBlock;
324  sizeOuterColumn = 1LL << qureg.numQubitsRepresented;
325  sizeOuterHalfColumn = sizeOuterColumn >> 1;
326 
327 # ifdef _OPENMP
328 # pragma omp parallel \
329  default (none) \
330  shared (sizeInnerBlock,sizeInnerHalfBlock,sizeOuterColumn,sizeOuterHalfColumn,qureg,damping, retain, dephase) \
331  private (thisTask,thisInnerBlock,thisOuterColumn,thisIndex,thisIndexInOuterColumn, \
332  thisIndexInInnerBlock,outerBit, stateBit)
333 # endif
334  {
335 # ifdef _OPENMP
336 # pragma omp for schedule (static)
337 # endif
338  // thisTask iterates over half the elements in this process' chunk of the density matrix
339  // treat this as iterating over all columns, then iterating over half the values
340  // within one column.
341  // If this function has been called, this process' chunk contains half an
342  // outer block or less
343  for (thisTask=0; thisTask<numTasks; thisTask++) {
344  // we want to process all columns in the density matrix,
345  // updating the values for half of each column (one half of each inner block)
346  thisOuterColumn = thisTask / sizeOuterHalfColumn;
347  thisIndexInOuterColumn = thisTask&(sizeOuterHalfColumn-1); // thisTask % sizeOuterHalfColumn
348  thisInnerBlock = thisIndexInOuterColumn/sizeInnerHalfBlock;
349  // get index in state vector corresponding to upper inner block
350  thisIndexInInnerBlock = thisTask&(sizeInnerHalfBlock-1); // thisTask % sizeInnerHalfBlock
351  thisIndex = thisOuterColumn*sizeOuterColumn + thisInnerBlock*sizeInnerBlock
352  + thisIndexInInnerBlock;
353  // check if we are in the upper or lower half of an outer block
354  outerBit = extractBit(targetQubit, (thisIndex+qureg.numAmpsPerChunk*qureg.chunkId)>>qureg.numQubitsRepresented);
355  // if we are in the lower half of an outer block, shift to be in the lower half
356  // of the inner block as well (we want to dephase |0><0| and |1><1| only)
357  thisIndex += outerBit*(sizeInnerHalfBlock);
358 
359  // NOTE: at this point thisIndex should be the index of the element we want to
360  // dephase in the chunk of the state vector on this process, in the
361  // density matrix representation.
362  // thisTask is the index of the pair element in pairStateVec
363 
364  // Extract state bit, is 0 if thisIndex corresponds to a state with 0 in the target qubit
365  // and is 1 if thisIndex corresponds to a state with 1 in the target qubit
366  stateBit = extractBit(targetQubit, (thisIndex+qureg.numAmpsPerChunk*qureg.chunkId));
367 
368  // state[thisIndex] = (1-depolLevel)*state[thisIndex] + depolLevel*(state[thisIndex]
369  // + pair[thisTask])/2
370  if(stateBit == 0){
371  qureg.stateVec.real[thisIndex] = qureg.stateVec.real[thisIndex] +
372  damping*( qureg.pairStateVec.real[thisTask]);
373 
374  qureg.stateVec.imag[thisIndex] = qureg.stateVec.imag[thisIndex] +
375  damping*( qureg.pairStateVec.imag[thisTask]);
376  } else{
377  qureg.stateVec.real[thisIndex] = retain*qureg.stateVec.real[thisIndex];
378 
379  qureg.stateVec.imag[thisIndex] = retain*qureg.stateVec.imag[thisIndex];
380  }
381  }
382  }
383 }
384 
385 void densmatr_mixTwoQubitDepolarisingLocal(Qureg qureg, int qubit1, int qubit2, qreal delta, qreal gamma) {
386  const long long int numTasks = qureg.numAmpsPerChunk;
387  long long int innerMaskQubit1 = 1LL << qubit1;
388  long long int outerMaskQubit1= 1LL << (qubit1 + qureg.numQubitsRepresented);
389  long long int totMaskQubit1 = innerMaskQubit1 | outerMaskQubit1;
390  long long int innerMaskQubit2 = 1LL << qubit2;
391  long long int outerMaskQubit2 = 1LL << (qubit2 + qureg.numQubitsRepresented);
392  long long int totMaskQubit2 = innerMaskQubit2 | outerMaskQubit2;
393 
394  long long int thisTask;
395  long long int partner;
396  long long int thisPatternQubit1, thisPatternQubit2;
397 
398  qreal real00, imag00;
399 
400 # ifdef _OPENMP
401 # pragma omp parallel \
402  default (none) \
403  shared (totMaskQubit1,totMaskQubit2,qureg,delta,gamma) \
404  private (thisTask,partner,thisPatternQubit1,thisPatternQubit2,real00,imag00)
405 # endif
406  {
407 
408 # ifdef _OPENMP
409 # pragma omp for schedule (static)
410 # endif
411  //--------------------------------------- STEP ONE ---------------------
412  for (thisTask=0; thisTask<numTasks; thisTask++){
413  thisPatternQubit1 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit1;
414  thisPatternQubit2 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit2;
415  if ((thisPatternQubit1==0) && ((thisPatternQubit2==0)
416  || (thisPatternQubit2==totMaskQubit2))){
417  //this element of form |...X...0...><...X...0...| for X either 0 or 1.
418  partner = thisTask | totMaskQubit1;
419  real00 = qureg.stateVec.real[thisTask];
420  imag00 = qureg.stateVec.imag[thisTask];
421 
422  qureg.stateVec.real[thisTask] = qureg.stateVec.real[thisTask]
423  + delta*qureg.stateVec.real[partner];
424  qureg.stateVec.imag[thisTask] = qureg.stateVec.imag[thisTask]
425  + delta*qureg.stateVec.imag[partner];
426 
427  qureg.stateVec.real[partner] = qureg.stateVec.real[partner] + delta*real00;
428  qureg.stateVec.imag[partner] = qureg.stateVec.imag[partner] + delta*imag00;
429 
430  }
431  }
432 # ifdef _OPENMP
433 # pragma omp for schedule (static)
434 # endif
435  //--------------------------------------- STEP TWO ---------------------
436  for (thisTask=0; thisTask<numTasks; thisTask++){
437  thisPatternQubit1 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit1;
438  thisPatternQubit2 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit2;
439  if ((thisPatternQubit2==0) && ((thisPatternQubit1==0)
440  || (thisPatternQubit1==totMaskQubit1))){
441  //this element of form |...0...X...><...0...X...| for X either 0 or 1.
442  partner = thisTask | totMaskQubit2;
443  real00 = qureg.stateVec.real[thisTask];
444  imag00 = qureg.stateVec.imag[thisTask];
445 
446  qureg.stateVec.real[thisTask] = qureg.stateVec.real[thisTask]
447  + delta*qureg.stateVec.real[partner];
448  qureg.stateVec.imag[thisTask] = qureg.stateVec.imag[thisTask]
449  + delta*qureg.stateVec.imag[partner];
450 
451  qureg.stateVec.real[partner] = qureg.stateVec.real[partner] + delta*real00;
452  qureg.stateVec.imag[partner] = qureg.stateVec.imag[partner] + delta*imag00;
453 
454  }
455  }
456 
457 # ifdef _OPENMP
458 # pragma omp for schedule (static)
459 # endif
460  //--------------------------------------- STEP THREE ---------------------
461  for (thisTask=0; thisTask<numTasks; thisTask++){
462  thisPatternQubit1 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit1;
463  thisPatternQubit2 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit2;
464  if ((thisPatternQubit2==0) && ((thisPatternQubit1==0)
465  || (thisPatternQubit1==totMaskQubit1))){
466  //this element of form |...0...X...><...0...X...| for X either 0 or 1.
467  partner = thisTask | totMaskQubit2;
468  partner = partner ^ totMaskQubit1;
469  real00 = qureg.stateVec.real[thisTask];
470  imag00 = qureg.stateVec.imag[thisTask];
471 
472  qureg.stateVec.real[thisTask] = gamma * (qureg.stateVec.real[thisTask]
473  + delta*qureg.stateVec.real[partner]);
474  qureg.stateVec.imag[thisTask] = gamma * (qureg.stateVec.imag[thisTask]
475  + delta*qureg.stateVec.imag[partner]);
476 
477  qureg.stateVec.real[partner] = gamma * (qureg.stateVec.real[partner]
478  + delta*real00);
479  qureg.stateVec.imag[partner] = gamma * (qureg.stateVec.imag[partner]
480  + delta*imag00);
481 
482  }
483  }
484  }
485 }
486 
487 void densmatr_mixTwoQubitDepolarisingLocalPart1(Qureg qureg, int qubit1, int qubit2, qreal delta) {
488  const long long int numTasks = qureg.numAmpsPerChunk;
489  long long int innerMaskQubit1 = 1LL << qubit1;
490  long long int outerMaskQubit1= 1LL << (qubit1 + qureg.numQubitsRepresented);
491  long long int totMaskQubit1 = innerMaskQubit1 | outerMaskQubit1;
492  long long int innerMaskQubit2 = 1LL << qubit2;
493  long long int outerMaskQubit2 = 1LL << (qubit2 + qureg.numQubitsRepresented);
494  long long int totMaskQubit2 = innerMaskQubit2 | outerMaskQubit2;
495  // correct for being in a particular chunk
496  //totMaskQubit2 = totMaskQubit2&(qureg.numAmpsPerChunk-1); // totMaskQubit2 % numAmpsPerChunk
497 
498 
499  long long int thisTask;
500  long long int partner;
501  long long int thisPatternQubit1, thisPatternQubit2;
502 
503  qreal real00, imag00;
504 
505 # ifdef _OPENMP
506 # pragma omp parallel \
507  default (none) \
508  shared (totMaskQubit1,totMaskQubit2,qureg,delta) \
509  private (thisTask,partner,thisPatternQubit1,thisPatternQubit2,real00,imag00)
510 # endif
511  {
512 
513 # ifdef _OPENMP
514 # pragma omp for schedule (static)
515 # endif
516  //--------------------------------------- STEP ONE ---------------------
517  for (thisTask=0; thisTask<numTasks; thisTask ++){
518  thisPatternQubit1 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit1;
519  thisPatternQubit2 = (thisTask+qureg.numAmpsPerChunk*qureg.chunkId)&totMaskQubit2;
520  if ((thisPatternQubit1==0) && ((thisPatternQubit2==0)
521  || (thisPatternQubit2==totMaskQubit2))){
522  //this element of form |...X...0...><...X...0...| for X either 0 or 1.
523  partner = thisTask | totMaskQubit1;
524  real00 = qureg.stateVec.real[thisTask];
525  imag00 = qureg.stateVec.imag[thisTask];
526 
527  qureg.stateVec.real[thisTask] = qureg.stateVec.real[thisTask]
528  + delta*qureg.stateVec.real[partner];
529  qureg.stateVec.imag[thisTask] = qureg.stateVec.imag[thisTask]
530  + delta*qureg.stateVec.imag[partner];
531 
532  qureg.stateVec.real[partner] = qureg.stateVec.real[partner] + delta*real00;
533  qureg.stateVec.imag[partner] = qureg.stateVec.imag[partner] + delta*imag00;
534 
535  }
536  }
537  }
538 }
539 
540 void densmatr_mixTwoQubitDepolarisingDistributed(Qureg qureg, const int targetQubit,
541  const int qubit2, qreal delta, qreal gamma) {
542 
543  long long int sizeInnerBlockQ1, sizeInnerHalfBlockQ1;
544  long long int sizeInnerBlockQ2, sizeInnerHalfBlockQ2, sizeInnerQuarterBlockQ2;
545  long long int sizeOuterColumn, sizeOuterQuarterColumn;
546  long long int thisInnerBlockQ2,
547  thisOuterColumn, // current column in density matrix
548  thisIndex, // current index in (density matrix representation) state vector
549  thisIndexInOuterColumn,
550  thisIndexInInnerBlockQ1,
551  thisIndexInInnerBlockQ2,
552  thisInnerBlockQ1InInnerBlockQ2;
553  int outerBitQ1, outerBitQ2;
554 
555  long long int thisTask;
556  const long long int numTasks=qureg.numAmpsPerChunk>>2;
557 
558  // set dimensions
559  sizeInnerHalfBlockQ1 = 1LL << targetQubit;
560  sizeInnerHalfBlockQ2 = 1LL << qubit2;
561  sizeInnerQuarterBlockQ2 = sizeInnerHalfBlockQ2 >> 1;
562  sizeInnerBlockQ2 = sizeInnerHalfBlockQ2 << 1;
563  sizeInnerBlockQ1 = 2LL * sizeInnerHalfBlockQ1;
564  sizeOuterColumn = 1LL << qureg.numQubitsRepresented;
565  sizeOuterQuarterColumn = sizeOuterColumn >> 2;
566 
567 # ifdef _OPENMP
568 # pragma omp parallel \
569  default (none) \
570  shared (sizeInnerBlockQ1,sizeInnerHalfBlockQ1,sizeInnerBlockQ2,sizeInnerHalfBlockQ2,sizeInnerQuarterBlockQ2,\
571  sizeOuterColumn,sizeOuterQuarterColumn,qureg,delta,gamma) \
572  private (thisTask,thisInnerBlockQ2,thisInnerBlockQ1InInnerBlockQ2, \
573  thisOuterColumn,thisIndex,thisIndexInOuterColumn, \
574  thisIndexInInnerBlockQ1,thisIndexInInnerBlockQ2,outerBitQ1,outerBitQ2)
575 # endif
576  {
577 # ifdef _OPENMP
578 # pragma omp for schedule (static)
579 # endif
580  // thisTask iterates over half the elements in this process' chunk of the density matrix
581  // treat this as iterating over all columns, then iterating over half the values
582  // within one column.
583  // If this function has been called, this process' chunk contains half an
584  // outer block or less
585  for (thisTask=0; thisTask<numTasks; thisTask++) {
586  // we want to process all columns in the density matrix,
587  // updating the values for half of each column (one half of each inner block)
588  thisOuterColumn = thisTask / sizeOuterQuarterColumn;
589  // thisTask % sizeOuterQuarterColumn
590  thisIndexInOuterColumn = thisTask&(sizeOuterQuarterColumn-1);
591  thisInnerBlockQ2 = thisIndexInOuterColumn / sizeInnerQuarterBlockQ2;
592  // thisTask % sizeInnerQuarterBlockQ2;
593  thisIndexInInnerBlockQ2 = thisTask&(sizeInnerQuarterBlockQ2-1);
594  thisInnerBlockQ1InInnerBlockQ2 = thisIndexInInnerBlockQ2 / sizeInnerHalfBlockQ1;
595  // thisTask % sizeInnerHalfBlockQ1;
596  thisIndexInInnerBlockQ1 = thisTask&(sizeInnerHalfBlockQ1-1);
597 
598  // get index in state vector corresponding to upper inner block
599  thisIndex = thisOuterColumn*sizeOuterColumn + thisInnerBlockQ2*sizeInnerBlockQ2
600  + thisInnerBlockQ1InInnerBlockQ2*sizeInnerBlockQ1 + thisIndexInInnerBlockQ1;
601 
602  // check if we are in the upper or lower half of an outer block for Q1
603  outerBitQ1 = extractBit(targetQubit, (thisIndex+qureg.numAmpsPerChunk*qureg.chunkId)>>qureg.numQubitsRepresented);
604  // if we are in the lower half of an outer block, shift to be in the lower half
605  // of the inner block as well (we want to dephase |0><0| and |1><1| only)
606  thisIndex += outerBitQ1*(sizeInnerHalfBlockQ1);
607 
608  // check if we are in the upper or lower half of an outer block for Q2
609  outerBitQ2 = extractBit(qubit2, (thisIndex+qureg.numAmpsPerChunk*qureg.chunkId)>>qureg.numQubitsRepresented);
610  // if we are in the lower half of an outer block, shift to be in the lower half
611  // of the inner block as well (we want to dephase |0><0| and |1><1| only)
612  thisIndex += outerBitQ2*(sizeInnerQuarterBlockQ2<<1);
613 
614  // NOTE: at this point thisIndex should be the index of the element we want to
615  // dephase in the chunk of the state vector on this process, in the
616  // density matrix representation.
617  // thisTask is the index of the pair element in pairStateVec
618 
619 
620  // state[thisIndex] = (1-depolLevel)*state[thisIndex] + depolLevel*(state[thisIndex]
621  // + pair[thisTask])/2
622  // NOTE: must set gamma=1 if using this function for steps 1 or 2
623  qureg.stateVec.real[thisIndex] = gamma*(qureg.stateVec.real[thisIndex] +
624  delta*qureg.pairStateVec.real[thisTask]);
625  qureg.stateVec.imag[thisIndex] = gamma*(qureg.stateVec.imag[thisIndex] +
626  delta*qureg.pairStateVec.imag[thisTask]);
627  }
628  }
629 }
630 
631 void densmatr_mixTwoQubitDepolarisingQ1LocalQ2DistributedPart3(Qureg qureg, const int targetQubit,
632  const int qubit2, qreal delta, qreal gamma) {
633 
634  long long int sizeInnerBlockQ1, sizeInnerHalfBlockQ1;
635  long long int sizeInnerBlockQ2, sizeInnerHalfBlockQ2, sizeInnerQuarterBlockQ2;
636  long long int sizeOuterColumn, sizeOuterQuarterColumn;
637  long long int thisInnerBlockQ2,
638  thisOuterColumn, // current column in density matrix
639  thisIndex, // current index in (density matrix representation) state vector
640  thisIndexInPairVector,
641  thisIndexInOuterColumn,
642  thisIndexInInnerBlockQ1,
643  thisIndexInInnerBlockQ2,
644  thisInnerBlockQ1InInnerBlockQ2;
645  int outerBitQ1, outerBitQ2;
646 
647  long long int thisTask;
648  const long long int numTasks=qureg.numAmpsPerChunk>>2;
649 
650  // set dimensions
651  sizeInnerHalfBlockQ1 = 1LL << targetQubit;
652  sizeInnerHalfBlockQ2 = 1LL << qubit2;
653  sizeInnerQuarterBlockQ2 = sizeInnerHalfBlockQ2 >> 1;
654  sizeInnerBlockQ2 = sizeInnerHalfBlockQ2 << 1;
655  sizeInnerBlockQ1 = 2LL * sizeInnerHalfBlockQ1;
656  sizeOuterColumn = 1LL << qureg.numQubitsRepresented;
657  sizeOuterQuarterColumn = sizeOuterColumn >> 2;
658 
659 //# if 0
660 # ifdef _OPENMP
661 # pragma omp parallel \
662  default (none) \
663  shared (sizeInnerBlockQ1,sizeInnerHalfBlockQ1,sizeInnerBlockQ2,sizeInnerHalfBlockQ2,sizeInnerQuarterBlockQ2,\
664  sizeOuterColumn,sizeOuterQuarterColumn,qureg,delta,gamma) \
665  private (thisTask,thisInnerBlockQ2,thisInnerBlockQ1InInnerBlockQ2, \
666  thisOuterColumn,thisIndex,thisIndexInPairVector,thisIndexInOuterColumn, \
667  thisIndexInInnerBlockQ1,thisIndexInInnerBlockQ2,outerBitQ1,outerBitQ2)
668 # endif
669  {
670 # ifdef _OPENMP
671 # pragma omp for schedule (static)
672 # endif
673 //# endif
674  // thisTask iterates over half the elements in this process' chunk of the density matrix
675  // treat this as iterating over all columns, then iterating over half the values
676  // within one column.
677  // If this function has been called, this process' chunk contains half an
678  // outer block or less
679  for (thisTask=0; thisTask<numTasks; thisTask++) {
680  // we want to process all columns in the density matrix,
681  // updating the values for half of each column (one half of each inner block)
682  thisOuterColumn = thisTask / sizeOuterQuarterColumn;
683  // thisTask % sizeOuterQuarterColumn
684  thisIndexInOuterColumn = thisTask&(sizeOuterQuarterColumn-1);
685  thisInnerBlockQ2 = thisIndexInOuterColumn / sizeInnerQuarterBlockQ2;
686  // thisTask % sizeInnerQuarterBlockQ2;
687  thisIndexInInnerBlockQ2 = thisTask&(sizeInnerQuarterBlockQ2-1);
688  thisInnerBlockQ1InInnerBlockQ2 = thisIndexInInnerBlockQ2 / sizeInnerHalfBlockQ1;
689  // thisTask % sizeInnerHalfBlockQ1;
690  thisIndexInInnerBlockQ1 = thisTask&(sizeInnerHalfBlockQ1-1);
691 
692  // get index in state vector corresponding to upper inner block
693  thisIndex = thisOuterColumn*sizeOuterColumn + thisInnerBlockQ2*sizeInnerBlockQ2
694  + thisInnerBlockQ1InInnerBlockQ2*sizeInnerBlockQ1 + thisIndexInInnerBlockQ1;
695 
696  // check if we are in the upper or lower half of an outer block for Q1
697  outerBitQ1 = extractBit(targetQubit, (thisIndex+qureg.numAmpsPerChunk*qureg.chunkId)>>qureg.numQubitsRepresented);
698  // if we are in the lower half of an outer block, shift to be in the lower half
699  // of the inner block as well (we want to dephase |0><0| and |1><1| only)
700  thisIndex += outerBitQ1*(sizeInnerHalfBlockQ1);
701 
702  // For part 3 we need to match elements such that (my Q1 != pair Q1) AND (my Q2 != pair Q2)
703  // Find correct index in pairStateVector
704  thisIndexInPairVector = thisTask + (1-outerBitQ1)*sizeInnerHalfBlockQ1*sizeOuterQuarterColumn -
705  outerBitQ1*sizeInnerHalfBlockQ1*sizeOuterQuarterColumn;
706 
707  // check if we are in the upper or lower half of an outer block for Q2
708  outerBitQ2 = extractBit(qubit2, (thisIndex+qureg.numAmpsPerChunk*qureg.chunkId)>>qureg.numQubitsRepresented);
709  // if we are in the lower half of an outer block, shift to be in the lower half
710  // of the inner block as well (we want to dephase |0><0| and |1><1| only)
711  thisIndex += outerBitQ2*(sizeInnerQuarterBlockQ2<<1);
712 
713 
714  // NOTE: at this point thisIndex should be the index of the element we want to
715  // dephase in the chunk of the state vector on this process, in the
716  // density matrix representation.
717 
718 
719  // state[thisIndex] = (1-depolLevel)*state[thisIndex] + depolLevel*(state[thisIndex]
720  // + pair[thisIndexInPairVector])/2
721  qureg.stateVec.real[thisIndex] = gamma*(qureg.stateVec.real[thisIndex] +
722  delta*qureg.pairStateVec.real[thisIndexInPairVector]);
723 
724  qureg.stateVec.imag[thisIndex] = gamma*(qureg.stateVec.imag[thisIndex] +
725  delta*qureg.pairStateVec.imag[thisIndexInPairVector]);
726  }
727  }
728 
729 }
730 
731 
732 /* Without nested parallelisation, only the outer most loops which call below are parallelised */
733 void zeroSomeAmps(Qureg qureg, long long int startInd, long long int numAmps) {
734  long long int i;
735 # ifdef _OPENMP
736 # pragma omp parallel for schedule (static)
737 # endif
738  for (i=startInd; i < startInd+numAmps; i++) {
739  qureg.stateVec.real[i] = 0;
740  qureg.stateVec.imag[i] = 0;
741  }
742 }
743 void normaliseSomeAmps(Qureg qureg, qreal norm, long long int startInd, long long int numAmps) {
744  long long int i;
745 # ifdef _OPENMP
746 # pragma omp parallel for schedule (static)
747 # endif
748  for (i=startInd; i < startInd+numAmps; i++) {
749  qureg.stateVec.real[i] /= norm;
750  qureg.stateVec.imag[i] /= norm;
751  }
752 }
753 void alternateNormZeroingSomeAmpBlocks(
754  Qureg qureg, qreal norm, int normFirst,
755  long long int startAmpInd, long long int numAmps, long long int blockSize
756 ) {
757  long long int numDubBlocks = numAmps / (2*blockSize);
758  long long int blockStartInd;
759 
760  if (normFirst) {
761  long long int dubBlockInd;
762 # ifdef _OPENMP
763 # pragma omp parallel for schedule (static) private (blockStartInd)
764 # endif
765  for (dubBlockInd=0; dubBlockInd < numDubBlocks; dubBlockInd++) {
766  blockStartInd = startAmpInd + dubBlockInd*2*blockSize;
767  normaliseSomeAmps(qureg, norm, blockStartInd, blockSize); // |0><0|
768  zeroSomeAmps( qureg, blockStartInd + blockSize, blockSize);
769  }
770  } else {
771  long long int dubBlockInd;
772 # ifdef _OPENMP
773 # pragma omp parallel for schedule (static) private (blockStartInd)
774 # endif
775  for (dubBlockInd=0; dubBlockInd < numDubBlocks; dubBlockInd++) {
776  blockStartInd = startAmpInd + dubBlockInd*2*blockSize;
777  zeroSomeAmps( qureg, blockStartInd, blockSize);
778  normaliseSomeAmps(qureg, norm, blockStartInd + blockSize, blockSize); // |1><1|
779  }
780  }
781 }
782 
784 void densmatr_collapseToKnownProbOutcome(Qureg qureg, const int measureQubit, int outcome, qreal totalStateProb) {
785 
786  // only (global) indices (as bit sequence): '* outcome *(n+q) outcome *q are spared
787  // where n = measureQubit, q = qureg.numQubitsRepresented.
788  // We can thus step in blocks of 2^q+n, killing every second, and inside the others,
789  // stepping in sub-blocks of 2^q, killing every second.
790  // When outcome=1, we offset the start of these blocks by their size.
791  long long int innerBlockSize = (1LL << measureQubit);
792  long long int outerBlockSize = (1LL << (measureQubit + qureg.numQubitsRepresented));
793 
794  // Because there are 2^a number of nodes(/chunks), each node will contain 2^b number of blocks,
795  // or each block will span 2^c number of nodes. Similarly for the innerblocks.
796  long long int locNumAmps = qureg.numAmpsPerChunk;
797  long long int globalStartInd = qureg.chunkId * locNumAmps;
798  int innerBit = extractBit(measureQubit, globalStartInd);
799  int outerBit = extractBit(measureQubit + qureg.numQubitsRepresented, globalStartInd);
800 
801  // If this chunk's amps are entirely inside an outer block
802  if (locNumAmps <= outerBlockSize) {
803 
804  // if this is an undesired outer block, kill all elems
805  if (outerBit != outcome)
806  return zeroSomeAmps(qureg, 0, qureg.numAmpsPerChunk);
807 
808  // othwerwise, if this is a desired outer block, and also entirely an inner block
809  if (locNumAmps <= innerBlockSize) {
810 
811  // and that inner block is undesired, kill all elems
812  if (innerBit != outcome)
813  return zeroSomeAmps(qureg, 0, qureg.numAmpsPerChunk);
814  // otherwise normalise all elems
815  else
816  return normaliseSomeAmps(qureg, totalStateProb, 0, qureg.numAmpsPerChunk);
817  }
818 
819  // otherwise this is a desired outer block which contains 2^a inner blocks; kill/renorm every second inner block
820  return alternateNormZeroingSomeAmpBlocks(
821  qureg, totalStateProb, innerBit==outcome, 0, qureg.numAmpsPerChunk, innerBlockSize);
822  }
823 
824  // Otherwise, this chunk's amps contain multiple outer blocks (and hence multiple inner blocks)
825  long long int numOuterDoubleBlocks = locNumAmps / (2*outerBlockSize);
826  long long int firstBlockInd;
827 
828  // alternate norming* and zeroing the outer blocks (with order based on the desired outcome)
829  // These loops aren't parallelised, since they could have 1 or 2 iterations and will prevent
830  // inner parallelisation
831  if (outerBit == outcome) {
832 
833  for (long long int outerDubBlockInd = 0; outerDubBlockInd < numOuterDoubleBlocks; outerDubBlockInd++) {
834  firstBlockInd = outerDubBlockInd*2*outerBlockSize;
835 
836  // *norm only the desired inner blocks in the desired outer block
837  alternateNormZeroingSomeAmpBlocks(
838  qureg, totalStateProb, innerBit==outcome,
839  firstBlockInd, outerBlockSize, innerBlockSize);
840 
841  // zero the undesired outer block
842  zeroSomeAmps(qureg, firstBlockInd + outerBlockSize, outerBlockSize);
843  }
844 
845  } else {
846 
847  for (long long int outerDubBlockInd = 0; outerDubBlockInd < numOuterDoubleBlocks; outerDubBlockInd++) {
848  firstBlockInd = outerDubBlockInd*2*outerBlockSize;
849 
850  // same thing but undesired outer blocks come first
851  zeroSomeAmps(qureg, firstBlockInd, outerBlockSize);
852  alternateNormZeroingSomeAmpBlocks(
853  qureg, totalStateProb, innerBit==outcome,
854  firstBlockInd + outerBlockSize, outerBlockSize, innerBlockSize);
855  }
856  }
857 
858 }
859 
860 qreal densmatr_calcPurityLocal(Qureg qureg) {
861 
862  /* sum of qureg^2, which is sum_i |qureg[i]|^2 */
863  long long int index;
864  long long int numAmps = qureg.numAmpsPerChunk;
865 
866  qreal trace = 0;
867  qreal *vecRe = qureg.stateVec.real;
868  qreal *vecIm = qureg.stateVec.imag;
869 
870 # ifdef _OPENMP
871 # pragma omp parallel \
872  shared (vecRe, vecIm, numAmps) \
873  private (index) \
874  reduction ( +:trace )
875 # endif
876  {
877 # ifdef _OPENMP
878 # pragma omp for schedule (static)
879 # endif
880  for (index=0LL; index<numAmps; index++) {
881 
882  trace += vecRe[index]*vecRe[index] + vecIm[index]*vecIm[index];
883  }
884  }
885 
886  return trace;
887 }
888 
889 void densmatr_mixDensityMatrix(Qureg combineQureg, qreal otherProb, Qureg otherQureg) {
890 
891  /* corresponding amplitudes live on the same node (same dimensions) */
892 
893  // unpack vars for OpenMP
894  qreal* combineVecRe = combineQureg.stateVec.real;
895  qreal* combineVecIm = combineQureg.stateVec.imag;
896  qreal* otherVecRe = otherQureg.stateVec.real;
897  qreal* otherVecIm = otherQureg.stateVec.imag;
898  long long int numAmps = combineQureg.numAmpsPerChunk;
899  long long int index;
900 
901 # ifdef _OPENMP
902 # pragma omp parallel \
903  default (none) \
904  shared (combineVecRe,combineVecIm,otherVecRe,otherVecIm, otherProb, numAmps) \
905  private (index)
906 # endif
907  {
908 # ifdef _OPENMP
909 # pragma omp for schedule (static)
910 # endif
911  for (index=0; index < numAmps; index++) {
912  combineVecRe[index] *= 1-otherProb;
913  combineVecIm[index] *= 1-otherProb;
914 
915  combineVecRe[index] += otherProb * otherVecRe[index];
916  combineVecIm[index] += otherProb * otherVecIm[index];
917  }
918  }
919 }
920 
922 qreal densmatr_calcHilbertSchmidtDistanceSquaredLocal(Qureg a, Qureg b) {
923 
924  long long int index;
925  long long int numAmps = a.numAmpsPerChunk;
926 
927  qreal *aRe = a.stateVec.real;
928  qreal *aIm = a.stateVec.imag;
929  qreal *bRe = b.stateVec.real;
930  qreal *bIm = b.stateVec.imag;
931 
932  qreal trace = 0;
933  qreal difRe, difIm;
934 
935 # ifdef _OPENMP
936 # pragma omp parallel \
937  shared (aRe,aIm, bRe,bIm, numAmps) \
938  private (index,difRe,difIm) \
939  reduction ( +:trace )
940 # endif
941  {
942 # ifdef _OPENMP
943 # pragma omp for schedule (static)
944 # endif
945  for (index=0LL; index<numAmps; index++) {
946 
947  difRe = aRe[index] - bRe[index];
948  difIm = aIm[index] - bIm[index];
949  trace += difRe*difRe + difIm*difIm;
950  }
951  }
952 
953  return trace;
954 }
955 
957 qreal densmatr_calcInnerProductLocal(Qureg a, Qureg b) {
958 
959  long long int index;
960  long long int numAmps = a.numAmpsPerChunk;
961 
962  qreal *aRe = a.stateVec.real;
963  qreal *aIm = a.stateVec.imag;
964  qreal *bRe = b.stateVec.real;
965  qreal *bIm = b.stateVec.imag;
966 
967  qreal trace = 0;
968 
969 # ifdef _OPENMP
970 # pragma omp parallel \
971  shared (aRe,aIm, bRe,bIm, numAmps) \
972  private (index) \
973  reduction ( +:trace )
974 # endif
975  {
976 # ifdef _OPENMP
977 # pragma omp for schedule (static)
978 # endif
979  for (index=0LL; index<numAmps; index++) {
980  trace += aRe[index]*bRe[index] + aIm[index]*bIm[index];
981  }
982  }
983 
984  return trace;
985 }
986 
987 
989 qreal densmatr_calcFidelityLocal(Qureg qureg, Qureg pureState) {
990 
991  /* Here, elements of pureState are not accessed (instead grabbed from qureg.pair).
992  * We only consult the attributes.
993  *
994  * qureg is a density matrix, and pureState is a statevector.
995  * Every node contains as many columns of qureg as amps by pureState.
996  * (each node contains an integer, exponent-of-2 number of whole columns of qureg)
997  * Ergo, this node contains columns:
998  * qureg.chunkID * pureState.numAmpsPerChunk to
999  * (qureg.chunkID + 1) * pureState.numAmpsPerChunk
1000  *
1001  * The first pureState.numAmpsTotal elements of qureg.pairStateVec are the
1002  * entire pureState state-vector
1003  */
1004 
1005  // unpack everything for OPENMP
1006  qreal* vecRe = qureg.pairStateVec.real;
1007  qreal* vecIm = qureg.pairStateVec.imag;
1008  qreal* densRe = qureg.stateVec.real;
1009  qreal* densIm = qureg.stateVec.imag;
1010 
1011  int row, col;
1012  int dim = (int) pureState.numAmpsTotal;
1013  int colsPerNode = (int) pureState.numAmpsPerChunk;
1014  // using only int, because density matrix has squared as many amps so its
1015  // iteration would be impossible if the pureStates numAmpsTotal didn't fit into int
1016 
1017  // starting GLOBAL column index of the qureg columns on this node
1018  int startCol = (int) (qureg.chunkId * pureState.numAmpsPerChunk);
1019 
1020  qreal densElemRe, densElemIm;
1021  qreal prefacRe, prefacIm;
1022  qreal rowSumRe, rowSumIm;
1023  qreal vecElemRe, vecElemIm;
1024 
1025  // quantity computed by this node
1026  qreal globalSumRe = 0; // imag-component is assumed zero
1027 
1028 # ifdef _OPENMP
1029 # pragma omp parallel \
1030  shared (vecRe,vecIm,densRe,densIm, dim,colsPerNode,startCol) \
1031  private (row,col, prefacRe,prefacIm, rowSumRe,rowSumIm, densElemRe,densElemIm, vecElemRe,vecElemIm) \
1032  reduction ( +:globalSumRe )
1033 # endif
1034  {
1035 # ifdef _OPENMP
1036 # pragma omp for schedule (static)
1037 # endif
1038  // indices of my GLOBAL row
1039  for (row=0; row < dim; row++) {
1040 
1041  // single element of conj(pureState)
1042  prefacRe = vecRe[row];
1043  prefacIm = - vecIm[row];
1044 
1045  rowSumRe = 0;
1046  rowSumIm = 0;
1047 
1048  // indices of my LOCAL column
1049  for (col=0; col < colsPerNode; col++) {
1050 
1051  // my local density element
1052  densElemRe = densRe[row + dim*col];
1053  densElemIm = densIm[row + dim*col];
1054 
1055  // state-vector element
1056  vecElemRe = vecRe[startCol + col];
1057  vecElemIm = vecIm[startCol + col];
1058 
1059  rowSumRe += densElemRe*vecElemRe - densElemIm*vecElemIm;
1060  rowSumIm += densElemRe*vecElemIm + densElemIm*vecElemRe;
1061  }
1062 
1063  globalSumRe += rowSumRe*prefacRe - rowSumIm*prefacIm;
1064  }
1065  }
1066 
1067  return globalSumRe;
1068 }
1069 
1070 Complex statevec_calcInnerProductLocal(Qureg bra, Qureg ket) {
1071 
1072  qreal innerProdReal = 0;
1073  qreal innerProdImag = 0;
1074 
1075  long long int index;
1076  long long int numAmps = bra.numAmpsPerChunk;
1077  qreal *braVecReal = bra.stateVec.real;
1078  qreal *braVecImag = bra.stateVec.imag;
1079  qreal *ketVecReal = ket.stateVec.real;
1080  qreal *ketVecImag = ket.stateVec.imag;
1081 
1082  qreal braRe, braIm, ketRe, ketIm;
1083 
1084 # ifdef _OPENMP
1085 # pragma omp parallel \
1086  shared (braVecReal, braVecImag, ketVecReal, ketVecImag, numAmps) \
1087  private (index, braRe, braIm, ketRe, ketIm) \
1088  reduction ( +:innerProdReal, innerProdImag )
1089 # endif
1090  {
1091 # ifdef _OPENMP
1092 # pragma omp for schedule (static)
1093 # endif
1094  for (index=0; index < numAmps; index++) {
1095  braRe = braVecReal[index];
1096  braIm = braVecImag[index];
1097  ketRe = ketVecReal[index];
1098  ketIm = ketVecImag[index];
1099 
1100  // conj(bra_i) * ket_i
1101  innerProdReal += braRe*ketRe + braIm*ketIm;
1102  innerProdImag += braRe*ketIm - braIm*ketRe;
1103  }
1104  }
1105 
1106  Complex innerProd;
1107  innerProd.real = innerProdReal;
1108  innerProd.imag = innerProdImag;
1109  return innerProd;
1110 }
1111 
1112 
1113 
1114 void densmatr_initClassicalState (Qureg qureg, long long int stateInd)
1115 {
1116  // dimension of the state vector
1117  long long int densityNumElems = qureg.numAmpsPerChunk;
1118 
1119  // Can't use qureg->stateVec as a private OMP var
1120  qreal *densityReal = qureg.stateVec.real;
1121  qreal *densityImag = qureg.stateVec.imag;
1122 
1123  // initialise the state to all zeros
1124  long long int index;
1125 # ifdef _OPENMP
1126 # pragma omp parallel \
1127  default (none) \
1128  shared (densityNumElems, densityReal, densityImag) \
1129  private (index)
1130 # endif
1131  {
1132 # ifdef _OPENMP
1133 # pragma omp for schedule (static)
1134 # endif
1135  for (index=0; index<densityNumElems; index++) {
1136  densityReal[index] = 0.0;
1137  densityImag[index] = 0.0;
1138  }
1139  }
1140 
1141  // index of the single density matrix elem to set non-zero
1142  long long int densityDim = 1LL << qureg.numQubitsRepresented;
1143  long long int densityInd = (densityDim + 1)*stateInd;
1144 
1145  // give the specified classical state prob 1
1146  if (qureg.chunkId == densityInd / densityNumElems){
1147  densityReal[densityInd % densityNumElems] = 1.0;
1148  densityImag[densityInd % densityNumElems] = 0.0;
1149  }
1150 }
1151 
1152 
1153 void densmatr_initPlusState (Qureg qureg)
1154 {
1155  // |+><+| = sum_i 1/sqrt(2^N) |i> 1/sqrt(2^N) <j| = sum_ij 1/2^N |i><j|
1156  long long int dim = (1LL << qureg.numQubitsRepresented);
1157  qreal probFactor = 1.0/((qreal) dim);
1158 
1159  // Can't use qureg->stateVec as a private OMP var
1160  qreal *densityReal = qureg.stateVec.real;
1161  qreal *densityImag = qureg.stateVec.imag;
1162 
1163  long long int index;
1164  long long int chunkSize = qureg.numAmpsPerChunk;
1165  // initialise the state to |+++..+++> = 1/normFactor {1, 1, 1, ...}
1166 # ifdef _OPENMP
1167 # pragma omp parallel \
1168  default (none) \
1169  shared (chunkSize, densityReal, densityImag, probFactor) \
1170  private (index)
1171 # endif
1172  {
1173 # ifdef _OPENMP
1174 # pragma omp for schedule (static)
1175 # endif
1176  for (index=0; index<chunkSize; index++) {
1177  densityReal[index] = probFactor;
1178  densityImag[index] = 0.0;
1179  }
1180  }
1181 }
1182 
1183 void densmatr_initPureStateLocal(Qureg targetQureg, Qureg copyQureg) {
1184 
1185  /* copyQureg amps aren't explicitly used - they're accessed through targetQureg.pair,
1186  * which contains the full pure statevector.
1187  * targetQureg has as many columns on node as copyQureg has amps
1188  */
1189 
1190  long long int colOffset = targetQureg.chunkId * copyQureg.numAmpsPerChunk;
1191  long long int colsPerNode = copyQureg.numAmpsPerChunk;
1192  long long int rowsPerNode = copyQureg.numAmpsTotal;
1193 
1194  // unpack vars for OpenMP
1195  qreal* vecRe = targetQureg.pairStateVec.real;
1196  qreal* vecIm = targetQureg.pairStateVec.imag;
1197  qreal* densRe = targetQureg.stateVec.real;
1198  qreal* densIm = targetQureg.stateVec.imag;
1199 
1200  long long int col, row, index;
1201 
1202  // a_i conj(a_j) |i><j|
1203  qreal ketRe, ketIm, braRe, braIm;
1204 
1205 # ifdef _OPENMP
1206 # pragma omp parallel \
1207  default (none) \
1208  shared (colOffset, colsPerNode,rowsPerNode, vecRe,vecIm,densRe,densIm) \
1209  private (col,row, ketRe,ketIm,braRe,braIm, index)
1210 # endif
1211  {
1212 # ifdef _OPENMP
1213 # pragma omp for schedule (static)
1214 # endif
1215  // local column
1216  for (col=0; col < colsPerNode; col++) {
1217 
1218  // global row
1219  for (row=0; row < rowsPerNode; row++) {
1220 
1221  // get pure state amps
1222  ketRe = vecRe[row];
1223  ketIm = vecIm[row];
1224  braRe = vecRe[col + colOffset];
1225  braIm = - vecIm[col + colOffset]; // minus for conjugation
1226 
1227  // update density matrix
1228  index = row + col*rowsPerNode; // local ind
1229  densRe[index] = ketRe*braRe - ketIm*braIm;
1230  densIm[index] = ketRe*braIm + ketIm*braRe;
1231  }
1232  }
1233  }
1234 }
1235 
1236 void statevec_setAmps(Qureg qureg, long long int startInd, qreal* reals, qreal* imags, long long int numAmps) {
1237 
1238  /* this is actually distributed, since the user's code runs on every node */
1239 
1240  // local start/end indices of the given amplitudes, assuming they fit in this chunk
1241  // these may be negative or above qureg.numAmpsPerChunk
1242  long long int localStartInd = startInd - qureg.chunkId*qureg.numAmpsPerChunk;
1243  long long int localEndInd = localStartInd + numAmps; // exclusive
1244 
1245  // add this to a local index to get corresponding elem in reals & imags
1246  long long int offset = qureg.chunkId*qureg.numAmpsPerChunk - startInd;
1247 
1248  // restrict these indices to fit into this chunk
1249  if (localStartInd < 0)
1250  localStartInd = 0;
1251  if (localEndInd > qureg.numAmpsPerChunk)
1252  localEndInd = qureg.numAmpsPerChunk;
1253  // they may now be out of order = no iterations
1254 
1255  // unpacking OpenMP vars
1256  long long int index;
1257  qreal* vecRe = qureg.stateVec.real;
1258  qreal* vecIm = qureg.stateVec.imag;
1259 
1260 # ifdef _OPENMP
1261 # pragma omp parallel \
1262  default (none) \
1263  shared (localStartInd,localEndInd, vecRe,vecIm, reals,imags, offset) \
1264  private (index)
1265 # endif
1266  {
1267 # ifdef _OPENMP
1268 # pragma omp for schedule (static)
1269 # endif
1270  // iterate these local inds - this might involve no iterations
1271  for (index=localStartInd; index < localEndInd; index++) {
1272  vecRe[index] = reals[index + offset];
1273  vecIm[index] = imags[index + offset];
1274  }
1275  }
1276 }
1277 
1278 void statevec_createQureg(Qureg *qureg, int numQubits, QuESTEnv env)
1279 {
1280  long long int numAmps = 1LL << numQubits;
1281  long long int numAmpsPerRank = numAmps/env.numRanks;
1282 
1283  if (numAmpsPerRank > SIZE_MAX) {
1284  printf("Could not allocate memory (cannot fit numAmps into size_t)!");
1285  exit (EXIT_FAILURE);
1286  }
1287 
1288  size_t arrSize = (size_t) (numAmpsPerRank * sizeof(*(qureg->stateVec.real)));
1289  qureg->stateVec.real = malloc(arrSize);
1290  qureg->stateVec.imag = malloc(arrSize);
1291  if (env.numRanks>1){
1292  qureg->pairStateVec.real = malloc(arrSize);
1293  qureg->pairStateVec.imag = malloc(arrSize);
1294  }
1295 
1296  if ( (!(qureg->stateVec.real) || !(qureg->stateVec.imag))
1297  && numAmpsPerRank ) {
1298  printf("Could not allocate memory!");
1299  exit (EXIT_FAILURE);
1300  }
1301 
1302  if ( env.numRanks>1 && (!(qureg->pairStateVec.real) || !(qureg->pairStateVec.imag))
1303  && numAmpsPerRank ) {
1304  printf("Could not allocate memory!");
1305  exit (EXIT_FAILURE);
1306  }
1307 
1308  qureg->numQubitsInStateVec = numQubits;
1309  qureg->numAmpsTotal = numAmps;
1310  qureg->numAmpsPerChunk = numAmpsPerRank;
1311  qureg->chunkId = env.rank;
1312  qureg->numChunks = env.numRanks;
1313  qureg->isDensityMatrix = 0;
1314 }
1315 
1316 void statevec_destroyQureg(Qureg qureg, QuESTEnv env){
1317 
1318  qureg.numQubitsInStateVec = 0;
1319  qureg.numAmpsTotal = 0;
1320  qureg.numAmpsPerChunk = 0;
1321 
1322  free(qureg.stateVec.real);
1323  free(qureg.stateVec.imag);
1324  if (env.numRanks>1){
1325  free(qureg.pairStateVec.real);
1326  free(qureg.pairStateVec.imag);
1327  }
1328  qureg.stateVec.real = NULL;
1329  qureg.stateVec.imag = NULL;
1330  qureg.pairStateVec.real = NULL;
1331  qureg.pairStateVec.imag = NULL;
1332 }
1333 
1334 void statevec_reportStateToScreen(Qureg qureg, QuESTEnv env, int reportRank){
1335  long long int index;
1336  int rank;
1337  if (qureg.numQubitsInStateVec<=5){
1338  for (rank=0; rank<qureg.numChunks; rank++){
1339  if (qureg.chunkId==rank){
1340  if (reportRank) {
1341  printf("Reporting state from rank %d [\n", qureg.chunkId);
1342  printf("real, imag\n");
1343  } else if (rank==0) {
1344  printf("Reporting state [\n");
1345  printf("real, imag\n");
1346  }
1347 
1348  for(index=0; index<qureg.numAmpsPerChunk; index++){
1349  //printf(REAL_STRING_FORMAT ", " REAL_STRING_FORMAT "\n", qureg.pairStateVec.real[index], qureg.pairStateVec.imag[index]);
1350  printf(REAL_STRING_FORMAT ", " REAL_STRING_FORMAT "\n", qureg.stateVec.real[index], qureg.stateVec.imag[index]);
1351  }
1352  if (reportRank || rank==qureg.numChunks-1) printf("]\n");
1353  }
1354  syncQuESTEnv(env);
1355  }
1356  } else printf("Error: reportStateToScreen will not print output for systems of more than 5 qubits.\n");
1357 }
1358 void statevec_getEnvironmentString(QuESTEnv env, Qureg qureg, char str[200]){
1359  int numThreads=1;
1360 # ifdef _OPENMP
1361  numThreads=omp_get_max_threads();
1362 # endif
1363  sprintf(str, "%dqubits_CPU_%dranksx%dthreads", qureg.numQubitsInStateVec, env.numRanks, numThreads);
1364 }
1365 
1366 void statevec_initBlankState (Qureg qureg)
1367 {
1368  long long int stateVecSize;
1369  long long int index;
1370 
1371  // dimension of the state vector
1372  stateVecSize = qureg.numAmpsPerChunk;
1373 
1374  // Can't use qureg->stateVec as a private OMP var
1375  qreal *stateVecReal = qureg.stateVec.real;
1376  qreal *stateVecImag = qureg.stateVec.imag;
1377 
1378  // initialise the state-vector to all-zeroes
1379 # ifdef _OPENMP
1380 # pragma omp parallel \
1381  default (none) \
1382  shared (stateVecSize, stateVecReal, stateVecImag) \
1383  private (index)
1384 # endif
1385  {
1386 # ifdef _OPENMP
1387 # pragma omp for schedule (static)
1388 # endif
1389  for (index=0; index<stateVecSize; index++) {
1390  stateVecReal[index] = 0.0;
1391  stateVecImag[index] = 0.0;
1392  }
1393  }
1394 }
1395 
1396 void statevec_initZeroState (Qureg qureg)
1397 {
1398  statevec_initBlankState(qureg);
1399  if (qureg.chunkId==0){
1400  // zero state |0000..0000> has probability 1
1401  qureg.stateVec.real[0] = 1.0;
1402  qureg.stateVec.imag[0] = 0.0;
1403  }
1404 }
1405 
1406 void statevec_initPlusState (Qureg qureg)
1407 {
1408  long long int chunkSize, stateVecSize;
1409  long long int index;
1410 
1411  // dimension of the state vector
1412  chunkSize = qureg.numAmpsPerChunk;
1413  stateVecSize = chunkSize*qureg.numChunks;
1414  qreal normFactor = 1.0/sqrt((qreal)stateVecSize);
1415 
1416  // Can't use qureg->stateVec as a private OMP var
1417  qreal *stateVecReal = qureg.stateVec.real;
1418  qreal *stateVecImag = qureg.stateVec.imag;
1419 
1420  // initialise the state to |+++..+++> = 1/normFactor {1, 1, 1, ...}
1421 # ifdef _OPENMP
1422 # pragma omp parallel \
1423  default (none) \
1424  shared (chunkSize, stateVecReal, stateVecImag, normFactor) \
1425  private (index)
1426 # endif
1427  {
1428 # ifdef _OPENMP
1429 # pragma omp for schedule (static)
1430 # endif
1431  for (index=0; index<chunkSize; index++) {
1432  stateVecReal[index] = normFactor;
1433  stateVecImag[index] = 0.0;
1434  }
1435  }
1436 }
1437 
1438 void statevec_initClassicalState (Qureg qureg, long long int stateInd)
1439 {
1440  long long int stateVecSize;
1441  long long int index;
1442 
1443  // dimension of the state vector
1444  stateVecSize = qureg.numAmpsPerChunk;
1445 
1446  // Can't use qureg->stateVec as a private OMP var
1447  qreal *stateVecReal = qureg.stateVec.real;
1448  qreal *stateVecImag = qureg.stateVec.imag;
1449 
1450  // initialise the state to vector to all zeros
1451 # ifdef _OPENMP
1452 # pragma omp parallel \
1453  default (none) \
1454  shared (stateVecSize, stateVecReal, stateVecImag) \
1455  private (index)
1456 # endif
1457  {
1458 # ifdef _OPENMP
1459 # pragma omp for schedule (static)
1460 # endif
1461  for (index=0; index<stateVecSize; index++) {
1462  stateVecReal[index] = 0.0;
1463  stateVecImag[index] = 0.0;
1464  }
1465  }
1466 
1467  // give the specified classical state prob 1
1468  if (qureg.chunkId == stateInd/stateVecSize){
1469  stateVecReal[stateInd % stateVecSize] = 1.0;
1470  stateVecImag[stateInd % stateVecSize] = 0.0;
1471  }
1472 }
1473 
1474 void statevec_cloneQureg(Qureg targetQureg, Qureg copyQureg) {
1475 
1476  // registers are equal sized, so nodes hold the same state-vector partitions
1477  long long int stateVecSize;
1478  long long int index;
1479 
1480  // dimension of the state vector
1481  stateVecSize = targetQureg.numAmpsPerChunk;
1482 
1483  // Can't use qureg->stateVec as a private OMP var
1484  qreal *targetStateVecReal = targetQureg.stateVec.real;
1485  qreal *targetStateVecImag = targetQureg.stateVec.imag;
1486  qreal *copyStateVecReal = copyQureg.stateVec.real;
1487  qreal *copyStateVecImag = copyQureg.stateVec.imag;
1488 
1489  // initialise the state to |0000..0000>
1490 # ifdef _OPENMP
1491 # pragma omp parallel \
1492  default (none) \
1493  shared (stateVecSize, targetStateVecReal, targetStateVecImag, copyStateVecReal, copyStateVecImag) \
1494  private (index)
1495 # endif
1496  {
1497 # ifdef _OPENMP
1498 # pragma omp for schedule (static)
1499 # endif
1500  for (index=0; index<stateVecSize; index++) {
1501  targetStateVecReal[index] = copyStateVecReal[index];
1502  targetStateVecImag[index] = copyStateVecImag[index];
1503  }
1504  }
1505 }
1506 
1513 void statevec_initStateOfSingleQubit(Qureg *qureg, int qubitId, int outcome)
1514 {
1515  long long int chunkSize, stateVecSize;
1516  long long int index;
1517  int bit;
1518  const long long int chunkId=qureg->chunkId;
1519 
1520  // dimension of the state vector
1521  chunkSize = qureg->numAmpsPerChunk;
1522  stateVecSize = chunkSize*qureg->numChunks;
1523  qreal normFactor = 1.0/sqrt((qreal)stateVecSize/2.0);
1524 
1525  // Can't use qureg->stateVec as a private OMP var
1526  qreal *stateVecReal = qureg->stateVec.real;
1527  qreal *stateVecImag = qureg->stateVec.imag;
1528 
1529  // initialise the state to |0000..0000>
1530 # ifdef _OPENMP
1531 # pragma omp parallel \
1532  default (none) \
1533  shared (chunkSize, stateVecReal, stateVecImag, normFactor, qubitId, outcome) \
1534  private (index, bit)
1535 # endif
1536  {
1537 # ifdef _OPENMP
1538 # pragma omp for schedule (static)
1539 # endif
1540  for (index=0; index<chunkSize; index++) {
1541  bit = extractBit(qubitId, index+chunkId*chunkSize);
1542  if (bit==outcome) {
1543  stateVecReal[index] = normFactor;
1544  stateVecImag[index] = 0.0;
1545  } else {
1546  stateVecReal[index] = 0.0;
1547  stateVecImag[index] = 0.0;
1548  }
1549  }
1550  }
1551 }
1552 
1553 
1559 void statevec_initDebugState (Qureg qureg)
1560 {
1561  long long int chunkSize;
1562  long long int index;
1563  long long int indexOffset;
1564 
1565  // dimension of the state vector
1566  chunkSize = qureg.numAmpsPerChunk;
1567 
1568  // Can't use qureg->stateVec as a private OMP var
1569  qreal *stateVecReal = qureg.stateVec.real;
1570  qreal *stateVecImag = qureg.stateVec.imag;
1571 
1572  indexOffset = chunkSize * qureg.chunkId;
1573 
1574  // initialise the state to |0000..0000>
1575 # ifdef _OPENMP
1576 # pragma omp parallel \
1577  default (none) \
1578  shared (chunkSize, stateVecReal, stateVecImag, indexOffset) \
1579  private (index)
1580 # endif
1581  {
1582 # ifdef _OPENMP
1583 # pragma omp for schedule (static)
1584 # endif
1585  for (index=0; index<chunkSize; index++) {
1586  stateVecReal[index] = ((indexOffset + index)*2.0)/10.0;
1587  stateVecImag[index] = ((indexOffset + index)*2.0+1.0)/10.0;
1588  }
1589  }
1590 }
1591 
1592 // returns 1 if successful, else 0
1593 int statevec_initStateFromSingleFile(Qureg *qureg, char filename[200], QuESTEnv env){
1594  long long int chunkSize, stateVecSize;
1595  long long int indexInChunk, totalIndex;
1596 
1597  chunkSize = qureg->numAmpsPerChunk;
1598  stateVecSize = chunkSize*qureg->numChunks;
1599 
1600  qreal *stateVecReal = qureg->stateVec.real;
1601  qreal *stateVecImag = qureg->stateVec.imag;
1602 
1603  FILE *fp;
1604  char line[200];
1605 
1606  for (int rank=0; rank<(qureg->numChunks); rank++){
1607  if (rank==qureg->chunkId){
1608  fp = fopen(filename, "r");
1609 
1610  // indicate file open failure
1611  if (fp == NULL)
1612  return 0;
1613 
1614  indexInChunk = 0; totalIndex = 0;
1615  while (fgets(line, sizeof(char)*200, fp) != NULL && totalIndex<stateVecSize){
1616  if (line[0]!='#'){
1617  int chunkId = (int) (totalIndex/chunkSize);
1618  if (chunkId==qureg->chunkId){
1619  # if QuEST_PREC==1
1620  sscanf(line, "%f, %f", &(stateVecReal[indexInChunk]),
1621  &(stateVecImag[indexInChunk]));
1622  # elif QuEST_PREC==2
1623  sscanf(line, "%lf, %lf", &(stateVecReal[indexInChunk]),
1624  &(stateVecImag[indexInChunk]));
1625  # elif QuEST_PREC==4
1626  sscanf(line, "%Lf, %Lf", &(stateVecReal[indexInChunk]),
1627  &(stateVecImag[indexInChunk]));
1628  # endif
1629  indexInChunk += 1;
1630  }
1631  totalIndex += 1;
1632  }
1633  }
1634  fclose(fp);
1635  }
1636  syncQuESTEnv(env);
1637  }
1638 
1639  // indicate success
1640  return 1;
1641 }
1642 
1643 int statevec_compareStates(Qureg mq1, Qureg mq2, qreal precision){
1644  qreal diff;
1645  long long int chunkSize = mq1.numAmpsPerChunk;
1646 
1647  for (long long int i=0; i<chunkSize; i++){
1648  diff = absReal(mq1.stateVec.real[i] - mq2.stateVec.real[i]);
1649  if (diff>precision) return 0;
1650  diff = absReal(mq1.stateVec.imag[i] - mq2.stateVec.imag[i]);
1651  if (diff>precision) return 0;
1652  }
1653  return 1;
1654 }
1655 
1656 void statevec_compactUnitaryLocal (Qureg qureg, const int targetQubit, Complex alpha, Complex beta)
1657 {
1658  long long int sizeBlock, sizeHalfBlock;
1659  long long int thisBlock, // current block
1660  indexUp,indexLo; // current index and corresponding index in lower half block
1661 
1662  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
1663  long long int thisTask;
1664  const long long int numTasks=qureg.numAmpsPerChunk>>1;
1665 
1666  // set dimensions
1667  sizeHalfBlock = 1LL << targetQubit;
1668  sizeBlock = 2LL * sizeHalfBlock;
1669 
1670  // Can't use qureg.stateVec as a private OMP var
1671  qreal *stateVecReal = qureg.stateVec.real;
1672  qreal *stateVecImag = qureg.stateVec.imag;
1673  qreal alphaImag=alpha.imag, alphaReal=alpha.real;
1674  qreal betaImag=beta.imag, betaReal=beta.real;
1675 
1676 # ifdef _OPENMP
1677 # pragma omp parallel \
1678  default (none) \
1679  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, alphaReal,alphaImag, betaReal,betaImag) \
1680  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp,stateRealLo,stateImagLo)
1681 # endif
1682  {
1683 # ifdef _OPENMP
1684 # pragma omp for schedule (static)
1685 # endif
1686  for (thisTask=0; thisTask<numTasks; thisTask++) {
1687 
1688  thisBlock = thisTask / sizeHalfBlock;
1689  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
1690  indexLo = indexUp + sizeHalfBlock;
1691 
1692  // store current state vector values in temp variables
1693  stateRealUp = stateVecReal[indexUp];
1694  stateImagUp = stateVecImag[indexUp];
1695 
1696  stateRealLo = stateVecReal[indexLo];
1697  stateImagLo = stateVecImag[indexLo];
1698 
1699  // state[indexUp] = alpha * state[indexUp] - conj(beta) * state[indexLo]
1700  stateVecReal[indexUp] = alphaReal*stateRealUp - alphaImag*stateImagUp
1701  - betaReal*stateRealLo - betaImag*stateImagLo;
1702  stateVecImag[indexUp] = alphaReal*stateImagUp + alphaImag*stateRealUp
1703  - betaReal*stateImagLo + betaImag*stateRealLo;
1704 
1705  // state[indexLo] = beta * state[indexUp] + conj(alpha) * state[indexLo]
1706  stateVecReal[indexLo] = betaReal*stateRealUp - betaImag*stateImagUp
1707  + alphaReal*stateRealLo + alphaImag*stateImagLo;
1708  stateVecImag[indexLo] = betaReal*stateImagUp + betaImag*stateRealUp
1709  + alphaReal*stateImagLo - alphaImag*stateRealLo;
1710  }
1711  }
1712 
1713 }
1714 
1715 void statevec_multiControlledTwoQubitUnitaryLocal(Qureg qureg, long long int ctrlMask, const int q1, const int q2, ComplexMatrix4 u) {
1716 
1717  // can't use qureg.stateVec as a private OMP var
1718  qreal *reVec = qureg.stateVec.real;
1719  qreal *imVec = qureg.stateVec.imag;
1720 
1721  // the global (between all nodes) index of this node's start index
1722  long long int globalIndStart = qureg.chunkId*qureg.numAmpsPerChunk;
1723 
1724  long long int numTasks = qureg.numAmpsPerChunk >> 2; // each iteration updates 4 amplitudes
1725  long long int thisTask;
1726  long long int thisGlobalInd00;
1727  long long int ind00, ind01, ind10, ind11;
1728  qreal re00, re01, re10, re11;
1729  qreal im00, im01, im10, im11;
1730 
1731 # ifdef _OPENMP
1732 # pragma omp parallel \
1733  default (none) \
1734  shared (reVec,imVec,globalIndStart,numTasks,ctrlMask,u) \
1735  private (thisTask, thisGlobalInd00, ind00,ind01,ind10,ind11, re00,re01,re10,re11, im00,im01,im10,im11)
1736 # endif
1737  {
1738 # ifdef _OPENMP
1739 # pragma omp for schedule (static)
1740 # endif
1741  for (thisTask=0; thisTask<numTasks; thisTask++) {
1742 
1743  // determine ind00 of |..0..0..>
1744  ind00 = insertTwoZeroBits(thisTask, q1, q2);
1745 
1746  // skip amplitude if controls aren't in 1 state (overloaded for speed)
1747  thisGlobalInd00 = ind00 + globalIndStart;
1748  if (ctrlMask && ((ctrlMask & thisGlobalInd00) != ctrlMask))
1749  continue;
1750 
1751  // inds of |..0..1..>, |..1..0..> and |..1..1..>
1752  ind01 = flipBit(ind00, q1);
1753  ind10 = flipBit(ind00, q2);
1754  ind11 = flipBit(ind01, q2);
1755 
1756  // extract statevec amplitudes
1757  re00 = reVec[ind00]; im00 = imVec[ind00];
1758  re01 = reVec[ind01]; im01 = imVec[ind01];
1759  re10 = reVec[ind10]; im10 = imVec[ind10];
1760  re11 = reVec[ind11]; im11 = imVec[ind11];
1761 
1762  // apply u * {amp00, amp01, amp10, amp11}
1763  reVec[ind00] =
1764  u.real[0][0]*re00 - u.imag[0][0]*im00 +
1765  u.real[0][1]*re01 - u.imag[0][1]*im01 +
1766  u.real[0][2]*re10 - u.imag[0][2]*im10 +
1767  u.real[0][3]*re11 - u.imag[0][3]*im11;
1768  imVec[ind00] =
1769  u.imag[0][0]*re00 + u.real[0][0]*im00 +
1770  u.imag[0][1]*re01 + u.real[0][1]*im01 +
1771  u.imag[0][2]*re10 + u.real[0][2]*im10 +
1772  u.imag[0][3]*re11 + u.real[0][3]*im11;
1773 
1774  reVec[ind01] =
1775  u.real[1][0]*re00 - u.imag[1][0]*im00 +
1776  u.real[1][1]*re01 - u.imag[1][1]*im01 +
1777  u.real[1][2]*re10 - u.imag[1][2]*im10 +
1778  u.real[1][3]*re11 - u.imag[1][3]*im11;
1779  imVec[ind01] =
1780  u.imag[1][0]*re00 + u.real[1][0]*im00 +
1781  u.imag[1][1]*re01 + u.real[1][1]*im01 +
1782  u.imag[1][2]*re10 + u.real[1][2]*im10 +
1783  u.imag[1][3]*re11 + u.real[1][3]*im11;
1784 
1785  reVec[ind10] =
1786  u.real[2][0]*re00 - u.imag[2][0]*im00 +
1787  u.real[2][1]*re01 - u.imag[2][1]*im01 +
1788  u.real[2][2]*re10 - u.imag[2][2]*im10 +
1789  u.real[2][3]*re11 - u.imag[2][3]*im11;
1790  imVec[ind10] =
1791  u.imag[2][0]*re00 + u.real[2][0]*im00 +
1792  u.imag[2][1]*re01 + u.real[2][1]*im01 +
1793  u.imag[2][2]*re10 + u.real[2][2]*im10 +
1794  u.imag[2][3]*re11 + u.real[2][3]*im11;
1795 
1796  reVec[ind11] =
1797  u.real[3][0]*re00 - u.imag[3][0]*im00 +
1798  u.real[3][1]*re01 - u.imag[3][1]*im01 +
1799  u.real[3][2]*re10 - u.imag[3][2]*im10 +
1800  u.real[3][3]*re11 - u.imag[3][3]*im11;
1801  imVec[ind11] =
1802  u.imag[3][0]*re00 + u.real[3][0]*im00 +
1803  u.imag[3][1]*re01 + u.real[3][1]*im01 +
1804  u.imag[3][2]*re10 + u.real[3][2]*im10 +
1805  u.imag[3][3]*re11 + u.real[3][3]*im11;
1806  }
1807  }
1808 }
1809 
1810 int qsortComp(const void *a, const void *b) {
1811  return *(int*)a - *(int*)b;
1812 }
1813 
1814 void statevec_multiControlledMultiQubitUnitaryLocal(Qureg qureg, long long int ctrlMask, int* targs, const int numTargs, ComplexMatrixN u)
1815 {
1816  // can't use qureg.stateVec as a private OMP var
1817  qreal *reVec = qureg.stateVec.real;
1818  qreal *imVec = qureg.stateVec.imag;
1819 
1820  long long int numTasks = qureg.numAmpsPerChunk >> numTargs; // kernel called on every 1 in 2^numTargs amplitudes
1821  long long int numTargAmps = 1 << u.numQubits; // num amps to be modified by each task
1822 
1823  // the global (between all nodes) index of this node's start index
1824  long long int globalIndStart = qureg.chunkId*qureg.numAmpsPerChunk;
1825 
1826  long long int thisTask;
1827  long long int thisInd00; // this thread's index of |..0..0..> (target qubits = 0)
1828  long long int thisGlobalInd00; // the global (between all nodes) index of this thread's |..0..0..> state
1829  long long int ind; // each thread's iteration of amplitudes to modify
1830  int i, t, r, c; // each thread's iteration of amps and targets
1831  qreal reElem, imElem; // each thread's iteration of u elements
1832 
1833  // each thread/task will record and modify numTargAmps amplitudes, privately
1834  // (of course, tasks eliminated by the ctrlMask won't edit their allocation)
1835  long long int ampInds[numTargAmps];
1836  qreal reAmps[numTargAmps];
1837  qreal imAmps[numTargAmps];
1838 
1839  // we need a sorted targets list to find thisInd00 for each task.
1840  // we can't modify targets, because the user-ordering of targets matters in u
1841  int sortedTargs[numTargs];
1842  for (int t=0; t < numTargs; t++)
1843  sortedTargs[t] = targs[t];
1844  qsort(sortedTargs, numTargs, sizeof(int), qsortComp);
1845 
1846 # ifdef _OPENMP
1847 # pragma omp parallel \
1848  default (none) \
1849  shared (reVec,imVec, numTasks,numTargAmps,globalIndStart, ctrlMask,targs,sortedTargs,u) \
1850  private (thisTask,thisInd00,thisGlobalInd00,ind,i,t,r,c,reElem,imElem, ampInds,reAmps,imAmps)
1851 # endif
1852  {
1853 # ifdef _OPENMP
1854 # pragma omp for schedule (static)
1855 # endif
1856  for (thisTask=0; thisTask<numTasks; thisTask++) {
1857 
1858  // find this task's start index (where all targs are 0)
1859  thisInd00 = thisTask;
1860  for (t=0; t < numTargs; t++)
1861  thisInd00 = insertZeroBit(thisInd00, sortedTargs[t]);
1862 
1863  // this task only modifies amplitudes if control qubits are 1 for this state
1864  thisGlobalInd00 = thisInd00 + globalIndStart;
1865  if (ctrlMask && ((ctrlMask & thisGlobalInd00) != ctrlMask))
1866  continue;
1867 
1868  // determine the indices and record values of this tasks's target amps
1869  for (i=0; i < numTargAmps; i++) {
1870 
1871  // get statevec index of current target qubit assignment
1872  ind = thisInd00;
1873  for (t=0; t < numTargs; t++)
1874  if (extractBit(t, i))
1875  ind = flipBit(ind, targs[t]);
1876 
1877  // update this tasks's private arrays
1878  ampInds[i] = ind;
1879  reAmps [i] = reVec[ind];
1880  imAmps [i] = imVec[ind];
1881  }
1882 
1883  // modify this tasks's target amplitudes
1884  for (r=0; r < numTargAmps; r++) {
1885  ind = ampInds[r];
1886  reVec[ind] = 0;
1887  imVec[ind] = 0;
1888 
1889  for (c=0; c < numTargAmps; c++) {
1890  reElem = u.real[r][c];
1891  imElem = u.imag[r][c];
1892  reVec[ind] += reAmps[c]*reElem - imAmps[c]*imElem;
1893  imVec[ind] += reAmps[c]*imElem + imAmps[c]*reElem;
1894  }
1895  }
1896  }
1897  }
1898 }
1899 
1900 void statevec_unitaryLocal(Qureg qureg, const int targetQubit, ComplexMatrix2 u)
1901 {
1902  long long int sizeBlock, sizeHalfBlock;
1903  long long int thisBlock, // current block
1904  indexUp,indexLo; // current index and corresponding index in lower half block
1905 
1906  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
1907  long long int thisTask;
1908  const long long int numTasks=qureg.numAmpsPerChunk>>1;
1909 
1910  // set dimensions
1911  sizeHalfBlock = 1LL << targetQubit;
1912  sizeBlock = 2LL * sizeHalfBlock;
1913 
1914  // Can't use qureg.stateVec as a private OMP var
1915  qreal *stateVecReal = qureg.stateVec.real;
1916  qreal *stateVecImag = qureg.stateVec.imag;
1917 
1918 # ifdef _OPENMP
1919 # pragma omp parallel \
1920  default (none) \
1921  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, u) \
1922  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp,stateRealLo,stateImagLo)
1923 # endif
1924  {
1925 # ifdef _OPENMP
1926 # pragma omp for schedule (static)
1927 # endif
1928  for (thisTask=0; thisTask<numTasks; thisTask++) {
1929 
1930  thisBlock = thisTask / sizeHalfBlock;
1931  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
1932  indexLo = indexUp + sizeHalfBlock;
1933 
1934  // store current state vector values in temp variables
1935  stateRealUp = stateVecReal[indexUp];
1936  stateImagUp = stateVecImag[indexUp];
1937 
1938  stateRealLo = stateVecReal[indexLo];
1939  stateImagLo = stateVecImag[indexLo];
1940 
1941 
1942  // state[indexUp] = u00 * state[indexUp] + u01 * state[indexLo]
1943  stateVecReal[indexUp] = u.real[0][0]*stateRealUp - u.imag[0][0]*stateImagUp
1944  + u.real[0][1]*stateRealLo - u.imag[0][1]*stateImagLo;
1945  stateVecImag[indexUp] = u.real[0][0]*stateImagUp + u.imag[0][0]*stateRealUp
1946  + u.real[0][1]*stateImagLo + u.imag[0][1]*stateRealLo;
1947 
1948  // state[indexLo] = u10 * state[indexUp] + u11 * state[indexLo]
1949  stateVecReal[indexLo] = u.real[1][0]*stateRealUp - u.imag[1][0]*stateImagUp
1950  + u.real[1][1]*stateRealLo - u.imag[1][1]*stateImagLo;
1951  stateVecImag[indexLo] = u.real[1][0]*stateImagUp + u.imag[1][0]*stateRealUp
1952  + u.real[1][1]*stateImagLo + u.imag[1][1]*stateRealLo;
1953 
1954  }
1955  }
1956 }
1957 
1969 void statevec_compactUnitaryDistributed (Qureg qureg,
1970  Complex rot1, Complex rot2,
1971  ComplexArray stateVecUp,
1972  ComplexArray stateVecLo,
1973  ComplexArray stateVecOut)
1974 {
1975 
1976  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
1977  long long int thisTask;
1978  const long long int numTasks=qureg.numAmpsPerChunk;
1979 
1980  qreal rot1Real=rot1.real, rot1Imag=rot1.imag;
1981  qreal rot2Real=rot2.real, rot2Imag=rot2.imag;
1982  qreal *stateVecRealUp=stateVecUp.real, *stateVecImagUp=stateVecUp.imag;
1983  qreal *stateVecRealLo=stateVecLo.real, *stateVecImagLo=stateVecLo.imag;
1984  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
1985 
1986 # ifdef _OPENMP
1987 # pragma omp parallel \
1988  default (none) \
1989  shared (stateVecRealUp,stateVecImagUp,stateVecRealLo,stateVecImagLo,stateVecRealOut,stateVecImagOut, \
1990  rot1Real,rot1Imag, rot2Real,rot2Imag) \
1991  private (thisTask,stateRealUp,stateImagUp,stateRealLo,stateImagLo)
1992 # endif
1993  {
1994 # ifdef _OPENMP
1995 # pragma omp for schedule (static)
1996 # endif
1997  for (thisTask=0; thisTask<numTasks; thisTask++) {
1998  // store current state vector values in temp variables
1999  stateRealUp = stateVecRealUp[thisTask];
2000  stateImagUp = stateVecImagUp[thisTask];
2001 
2002  stateRealLo = stateVecRealLo[thisTask];
2003  stateImagLo = stateVecImagLo[thisTask];
2004 
2005  // state[indexUp] = alpha * state[indexUp] - conj(beta) * state[indexLo]
2006  stateVecRealOut[thisTask] = rot1Real*stateRealUp - rot1Imag*stateImagUp + rot2Real*stateRealLo + rot2Imag*stateImagLo;
2007  stateVecImagOut[thisTask] = rot1Real*stateImagUp + rot1Imag*stateRealUp + rot2Real*stateImagLo - rot2Imag*stateRealLo;
2008  }
2009  }
2010 }
2011 
2024 void statevec_unitaryDistributed (Qureg qureg,
2025  Complex rot1, Complex rot2,
2026  ComplexArray stateVecUp,
2027  ComplexArray stateVecLo,
2028  ComplexArray stateVecOut)
2029 {
2030 
2031  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2032  long long int thisTask;
2033  const long long int numTasks=qureg.numAmpsPerChunk;
2034 
2035  qreal rot1Real=rot1.real, rot1Imag=rot1.imag;
2036  qreal rot2Real=rot2.real, rot2Imag=rot2.imag;
2037  qreal *stateVecRealUp=stateVecUp.real, *stateVecImagUp=stateVecUp.imag;
2038  qreal *stateVecRealLo=stateVecLo.real, *stateVecImagLo=stateVecLo.imag;
2039  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2040 
2041 
2042 # ifdef _OPENMP
2043 # pragma omp parallel \
2044  default (none) \
2045  shared (stateVecRealUp,stateVecImagUp,stateVecRealLo,stateVecImagLo,stateVecRealOut,stateVecImagOut, \
2046  rot1Real, rot1Imag, rot2Real, rot2Imag) \
2047  private (thisTask,stateRealUp,stateImagUp,stateRealLo,stateImagLo)
2048 # endif
2049  {
2050 # ifdef _OPENMP
2051 # pragma omp for schedule (static)
2052 # endif
2053  for (thisTask=0; thisTask<numTasks; thisTask++) {
2054  // store current state vector values in temp variables
2055  stateRealUp = stateVecRealUp[thisTask];
2056  stateImagUp = stateVecImagUp[thisTask];
2057 
2058  stateRealLo = stateVecRealLo[thisTask];
2059  stateImagLo = stateVecImagLo[thisTask];
2060 
2061  stateVecRealOut[thisTask] = rot1Real*stateRealUp - rot1Imag*stateImagUp
2062  + rot2Real*stateRealLo - rot2Imag*stateImagLo;
2063  stateVecImagOut[thisTask] = rot1Real*stateImagUp + rot1Imag*stateRealUp
2064  + rot2Real*stateImagLo + rot2Imag*stateRealLo;
2065  }
2066  }
2067 }
2068 
2069 void statevec_controlledCompactUnitaryLocal (Qureg qureg, const int controlQubit, const int targetQubit,
2070  Complex alpha, Complex beta)
2071 {
2072  long long int sizeBlock, sizeHalfBlock;
2073  long long int thisBlock, // current block
2074  indexUp,indexLo; // current index and corresponding index in lower half block
2075 
2076  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2077  long long int thisTask;
2078  const long long int numTasks=qureg.numAmpsPerChunk>>1;
2079  const long long int chunkSize=qureg.numAmpsPerChunk;
2080  const long long int chunkId=qureg.chunkId;
2081 
2082  int controlBit;
2083 
2084  // set dimensions
2085  sizeHalfBlock = 1LL << targetQubit;
2086  sizeBlock = 2LL * sizeHalfBlock;
2087 
2088  // Can't use qureg.stateVec as a private OMP var
2089  qreal *stateVecReal = qureg.stateVec.real;
2090  qreal *stateVecImag = qureg.stateVec.imag;
2091  qreal alphaImag=alpha.imag, alphaReal=alpha.real;
2092  qreal betaImag=beta.imag, betaReal=beta.real;
2093 
2094 # ifdef _OPENMP
2095 # pragma omp parallel \
2096  default (none) \
2097  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, alphaReal,alphaImag, betaReal,betaImag) \
2098  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp,stateRealLo,stateImagLo,controlBit)
2099 # endif
2100  {
2101 # ifdef _OPENMP
2102 # pragma omp for schedule (static)
2103 # endif
2104  for (thisTask=0; thisTask<numTasks; thisTask++) {
2105 
2106  thisBlock = thisTask / sizeHalfBlock;
2107  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2108  indexLo = indexUp + sizeHalfBlock;
2109 
2110  controlBit = extractBit (controlQubit, indexUp+chunkId*chunkSize);
2111  if (controlBit){
2112  // store current state vector values in temp variables
2113  stateRealUp = stateVecReal[indexUp];
2114  stateImagUp = stateVecImag[indexUp];
2115 
2116  stateRealLo = stateVecReal[indexLo];
2117  stateImagLo = stateVecImag[indexLo];
2118 
2119  // state[indexUp] = alpha * state[indexUp] - conj(beta) * state[indexLo]
2120  stateVecReal[indexUp] = alphaReal*stateRealUp - alphaImag*stateImagUp
2121  - betaReal*stateRealLo - betaImag*stateImagLo;
2122  stateVecImag[indexUp] = alphaReal*stateImagUp + alphaImag*stateRealUp
2123  - betaReal*stateImagLo + betaImag*stateRealLo;
2124 
2125  // state[indexLo] = beta * state[indexUp] + conj(alpha) * state[indexLo]
2126  stateVecReal[indexLo] = betaReal*stateRealUp - betaImag*stateImagUp
2127  + alphaReal*stateRealLo + alphaImag*stateImagLo;
2128  stateVecImag[indexLo] = betaReal*stateImagUp + betaImag*stateRealUp
2129  + alphaReal*stateImagLo - alphaImag*stateRealLo;
2130  }
2131  }
2132  }
2133 
2134 }
2135 
2136 /* ctrlQubitsMask is a bit mask indicating which qubits are control Qubits
2137  * ctrlFlipMask is a bit mask indicating which control qubits should be 'flipped'
2138  * in the condition, i.e. they should have value 0 when the unitary is applied
2139  */
2140 void statevec_multiControlledUnitaryLocal(
2141  Qureg qureg, const int targetQubit,
2142  long long int ctrlQubitsMask, long long int ctrlFlipMask,
2143  ComplexMatrix2 u)
2144 {
2145  long long int sizeBlock, sizeHalfBlock;
2146  long long int thisBlock, // current block
2147  indexUp,indexLo; // current index and corresponding index in lower half block
2148 
2149  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2150  long long int thisTask;
2151  const long long int numTasks=qureg.numAmpsPerChunk>>1;
2152  const long long int chunkSize=qureg.numAmpsPerChunk;
2153  const long long int chunkId=qureg.chunkId;
2154 
2155  // set dimensions
2156  sizeHalfBlock = 1LL << targetQubit;
2157  sizeBlock = 2LL * sizeHalfBlock;
2158 
2159  // Can't use qureg.stateVec as a private OMP var
2160  qreal *stateVecReal = qureg.stateVec.real;
2161  qreal *stateVecImag = qureg.stateVec.imag;
2162 
2163 # ifdef _OPENMP
2164 # pragma omp parallel \
2165  default (none) \
2166  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, u, ctrlQubitsMask,ctrlFlipMask) \
2167  private (thisTask,thisBlock, indexUp,indexLo, stateRealUp,stateImagUp,stateRealLo,stateImagLo)
2168 # endif
2169  {
2170 # ifdef _OPENMP
2171 # pragma omp for schedule (static)
2172 # endif
2173  for (thisTask=0; thisTask<numTasks; thisTask++) {
2174 
2175  thisBlock = thisTask / sizeHalfBlock;
2176  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2177  indexLo = indexUp + sizeHalfBlock;
2178 
2179 
2180  // take the basis index, flip the designated (XOR) 'control' bits, AND with the controls.
2181  // if this equals the control mask, the control qubits have the desired values in the basis index
2182  if (ctrlQubitsMask == (ctrlQubitsMask & ((indexUp+chunkId*chunkSize) ^ ctrlFlipMask))) {
2183  // store current state vector values in temp variables
2184  stateRealUp = stateVecReal[indexUp];
2185  stateImagUp = stateVecImag[indexUp];
2186 
2187  stateRealLo = stateVecReal[indexLo];
2188  stateImagLo = stateVecImag[indexLo];
2189 
2190  // state[indexUp] = u00 * state[indexUp] + u01 * state[indexLo]
2191  stateVecReal[indexUp] = u.real[0][0]*stateRealUp - u.imag[0][0]*stateImagUp
2192  + u.real[0][1]*stateRealLo - u.imag[0][1]*stateImagLo;
2193  stateVecImag[indexUp] = u.real[0][0]*stateImagUp + u.imag[0][0]*stateRealUp
2194  + u.real[0][1]*stateImagLo + u.imag[0][1]*stateRealLo;
2195 
2196  // state[indexLo] = u10 * state[indexUp] + u11 * state[indexLo]
2197  stateVecReal[indexLo] = u.real[1][0]*stateRealUp - u.imag[1][0]*stateImagUp
2198  + u.real[1][1]*stateRealLo - u.imag[1][1]*stateImagLo;
2199  stateVecImag[indexLo] = u.real[1][0]*stateImagUp + u.imag[1][0]*stateRealUp
2200  + u.real[1][1]*stateImagLo + u.imag[1][1]*stateRealLo;
2201  }
2202  }
2203  }
2204 
2205 }
2206 
2207 void statevec_controlledUnitaryLocal(Qureg qureg, const int controlQubit, const int targetQubit,
2208  ComplexMatrix2 u)
2209 {
2210  long long int sizeBlock, sizeHalfBlock;
2211  long long int thisBlock, // current block
2212  indexUp,indexLo; // current index and corresponding index in lower half block
2213 
2214  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2215  long long int thisTask;
2216  const long long int numTasks=qureg.numAmpsPerChunk>>1;
2217  const long long int chunkSize=qureg.numAmpsPerChunk;
2218  const long long int chunkId=qureg.chunkId;
2219 
2220  int controlBit;
2221 
2222  // set dimensions
2223  sizeHalfBlock = 1LL << targetQubit;
2224  sizeBlock = 2LL * sizeHalfBlock;
2225 
2226  // Can't use qureg.stateVec as a private OMP var
2227  qreal *stateVecReal = qureg.stateVec.real;
2228  qreal *stateVecImag = qureg.stateVec.imag;
2229 
2230 # ifdef _OPENMP
2231 # pragma omp parallel \
2232  default (none) \
2233  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, u) \
2234  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp,stateRealLo,stateImagLo,controlBit)
2235 # endif
2236  {
2237 # ifdef _OPENMP
2238 # pragma omp for schedule (static)
2239 # endif
2240  for (thisTask=0; thisTask<numTasks; thisTask++) {
2241 
2242  thisBlock = thisTask / sizeHalfBlock;
2243  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2244  indexLo = indexUp + sizeHalfBlock;
2245 
2246  controlBit = extractBit (controlQubit, indexUp+chunkId*chunkSize);
2247  if (controlBit){
2248  // store current state vector values in temp variables
2249  stateRealUp = stateVecReal[indexUp];
2250  stateImagUp = stateVecImag[indexUp];
2251 
2252  stateRealLo = stateVecReal[indexLo];
2253  stateImagLo = stateVecImag[indexLo];
2254 
2255 
2256  // state[indexUp] = u00 * state[indexUp] + u01 * state[indexLo]
2257  stateVecReal[indexUp] = u.real[0][0]*stateRealUp - u.imag[0][0]*stateImagUp
2258  + u.real[0][1]*stateRealLo - u.imag[0][1]*stateImagLo;
2259  stateVecImag[indexUp] = u.real[0][0]*stateImagUp + u.imag[0][0]*stateRealUp
2260  + u.real[0][1]*stateImagLo + u.imag[0][1]*stateRealLo;
2261 
2262  // state[indexLo] = u10 * state[indexUp] + u11 * state[indexLo]
2263  stateVecReal[indexLo] = u.real[1][0]*stateRealUp - u.imag[1][0]*stateImagUp
2264  + u.real[1][1]*stateRealLo - u.imag[1][1]*stateImagLo;
2265  stateVecImag[indexLo] = u.real[1][0]*stateImagUp + u.imag[1][0]*stateRealUp
2266  + u.real[1][1]*stateImagLo + u.imag[1][1]*stateRealLo;
2267  }
2268  }
2269  }
2270 
2271 }
2272 
2285 void statevec_controlledCompactUnitaryDistributed (Qureg qureg, const int controlQubit,
2286  Complex rot1, Complex rot2,
2287  ComplexArray stateVecUp,
2288  ComplexArray stateVecLo,
2289  ComplexArray stateVecOut)
2290 {
2291 
2292  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2293  long long int thisTask;
2294  const long long int numTasks=qureg.numAmpsPerChunk;
2295  const long long int chunkSize=qureg.numAmpsPerChunk;
2296  const long long int chunkId=qureg.chunkId;
2297 
2298  int controlBit;
2299 
2300  qreal rot1Real=rot1.real, rot1Imag=rot1.imag;
2301  qreal rot2Real=rot2.real, rot2Imag=rot2.imag;
2302  qreal *stateVecRealUp=stateVecUp.real, *stateVecImagUp=stateVecUp.imag;
2303  qreal *stateVecRealLo=stateVecLo.real, *stateVecImagLo=stateVecLo.imag;
2304  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2305 
2306 # ifdef _OPENMP
2307 # pragma omp parallel \
2308  default (none) \
2309  shared (stateVecRealUp,stateVecImagUp,stateVecRealLo,stateVecImagLo,stateVecRealOut,stateVecImagOut, \
2310  rot1Real,rot1Imag, rot2Real,rot2Imag) \
2311  private (thisTask,stateRealUp,stateImagUp,stateRealLo,stateImagLo,controlBit)
2312 # endif
2313  {
2314 # ifdef _OPENMP
2315 # pragma omp for schedule (static)
2316 # endif
2317  for (thisTask=0; thisTask<numTasks; thisTask++) {
2318  controlBit = extractBit (controlQubit, thisTask+chunkId*chunkSize);
2319  if (controlBit){
2320  // store current state vector values in temp variables
2321  stateRealUp = stateVecRealUp[thisTask];
2322  stateImagUp = stateVecImagUp[thisTask];
2323 
2324  stateRealLo = stateVecRealLo[thisTask];
2325  stateImagLo = stateVecImagLo[thisTask];
2326 
2327  // state[indexUp] = alpha * state[indexUp] - conj(beta) * state[indexLo]
2328  stateVecRealOut[thisTask] = rot1Real*stateRealUp - rot1Imag*stateImagUp + rot2Real*stateRealLo + rot2Imag*stateImagLo;
2329  stateVecImagOut[thisTask] = rot1Real*stateImagUp + rot1Imag*stateRealUp + rot2Real*stateImagLo - rot2Imag*stateRealLo;
2330  }
2331  }
2332  }
2333 }
2334 
2347 void statevec_controlledUnitaryDistributed (Qureg qureg, const int controlQubit,
2348  Complex rot1, Complex rot2,
2349  ComplexArray stateVecUp,
2350  ComplexArray stateVecLo,
2351  ComplexArray stateVecOut)
2352 {
2353 
2354  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2355  long long int thisTask;
2356  const long long int numTasks=qureg.numAmpsPerChunk;
2357  const long long int chunkSize=qureg.numAmpsPerChunk;
2358  const long long int chunkId=qureg.chunkId;
2359 
2360  int controlBit;
2361 
2362  qreal rot1Real=rot1.real, rot1Imag=rot1.imag;
2363  qreal rot2Real=rot2.real, rot2Imag=rot2.imag;
2364  qreal *stateVecRealUp=stateVecUp.real, *stateVecImagUp=stateVecUp.imag;
2365  qreal *stateVecRealLo=stateVecLo.real, *stateVecImagLo=stateVecLo.imag;
2366  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2367 
2368 # ifdef _OPENMP
2369 # pragma omp parallel \
2370  default (none) \
2371  shared (stateVecRealUp,stateVecImagUp,stateVecRealLo,stateVecImagLo,stateVecRealOut,stateVecImagOut, \
2372  rot1Real,rot1Imag, rot2Real,rot2Imag) \
2373  private (thisTask,stateRealUp,stateImagUp,stateRealLo,stateImagLo,controlBit)
2374 # endif
2375  {
2376 # ifdef _OPENMP
2377 # pragma omp for schedule (static)
2378 # endif
2379  for (thisTask=0; thisTask<numTasks; thisTask++) {
2380  controlBit = extractBit (controlQubit, thisTask+chunkId*chunkSize);
2381  if (controlBit){
2382  // store current state vector values in temp variables
2383  stateRealUp = stateVecRealUp[thisTask];
2384  stateImagUp = stateVecImagUp[thisTask];
2385 
2386  stateRealLo = stateVecRealLo[thisTask];
2387  stateImagLo = stateVecImagLo[thisTask];
2388 
2389  stateVecRealOut[thisTask] = rot1Real*stateRealUp - rot1Imag*stateImagUp
2390  + rot2Real*stateRealLo - rot2Imag*stateImagLo;
2391  stateVecImagOut[thisTask] = rot1Real*stateImagUp + rot1Imag*stateRealUp
2392  + rot2Real*stateImagLo + rot2Imag*stateRealLo;
2393  }
2394  }
2395  }
2396 }
2397 
2413 void statevec_multiControlledUnitaryDistributed (
2414  Qureg qureg,
2415  const int targetQubit,
2416  long long int ctrlQubitsMask, long long int ctrlFlipMask,
2417  Complex rot1, Complex rot2,
2418  ComplexArray stateVecUp,
2419  ComplexArray stateVecLo,
2420  ComplexArray stateVecOut)
2421 {
2422 
2423  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2424  long long int thisTask;
2425  const long long int numTasks=qureg.numAmpsPerChunk;
2426  const long long int chunkSize=qureg.numAmpsPerChunk;
2427  const long long int chunkId=qureg.chunkId;
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, ctrlQubitsMask,ctrlFlipMask) \
2440  private (thisTask,stateRealUp,stateImagUp,stateRealLo,stateImagLo)
2441 # endif
2442  {
2443 # ifdef _OPENMP
2444 # pragma omp for schedule (static)
2445 # endif
2446  for (thisTask=0; thisTask<numTasks; thisTask++) {
2447  if (ctrlQubitsMask == (ctrlQubitsMask & ((thisTask+chunkId*chunkSize) ^ ctrlFlipMask))) {
2448  // store current state vector values in temp variables
2449  stateRealUp = stateVecRealUp[thisTask];
2450  stateImagUp = stateVecImagUp[thisTask];
2451 
2452  stateRealLo = stateVecRealLo[thisTask];
2453  stateImagLo = stateVecImagLo[thisTask];
2454 
2455  stateVecRealOut[thisTask] = rot1Real*stateRealUp - rot1Imag*stateImagUp
2456  + rot2Real*stateRealLo - rot2Imag*stateImagLo;
2457  stateVecImagOut[thisTask] = rot1Real*stateImagUp + rot1Imag*stateRealUp
2458  + rot2Real*stateImagLo + rot2Imag*stateRealLo;
2459  }
2460  }
2461  }
2462 }
2463 
2464 void statevec_pauliXLocal(Qureg qureg, const int targetQubit)
2465 {
2466  long long int sizeBlock, sizeHalfBlock;
2467  long long int thisBlock, // current block
2468  indexUp,indexLo; // current index and corresponding index in lower half block
2469 
2470  qreal stateRealUp,stateImagUp;
2471  long long int thisTask;
2472  const long long int numTasks=qureg.numAmpsPerChunk>>1;
2473 
2474  // set dimensions
2475  sizeHalfBlock = 1LL << targetQubit;
2476  sizeBlock = 2LL * sizeHalfBlock;
2477 
2478  // Can't use qureg.stateVec as a private OMP var
2479  qreal *stateVecReal = qureg.stateVec.real;
2480  qreal *stateVecImag = qureg.stateVec.imag;
2481 
2482 # ifdef _OPENMP
2483 # pragma omp parallel \
2484  default (none) \
2485  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag) \
2486  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp)
2487 # endif
2488  {
2489 # ifdef _OPENMP
2490 # pragma omp for schedule (static)
2491 # endif
2492  for (thisTask=0; thisTask<numTasks; thisTask++) {
2493  thisBlock = thisTask / sizeHalfBlock;
2494  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2495  indexLo = indexUp + sizeHalfBlock;
2496 
2497  stateRealUp = stateVecReal[indexUp];
2498  stateImagUp = stateVecImag[indexUp];
2499 
2500  stateVecReal[indexUp] = stateVecReal[indexLo];
2501  stateVecImag[indexUp] = stateVecImag[indexLo];
2502 
2503  stateVecReal[indexLo] = stateRealUp;
2504  stateVecImag[indexLo] = stateImagUp;
2505  }
2506  }
2507 
2508 }
2509 
2522 void statevec_pauliXDistributed (Qureg qureg,
2523  ComplexArray stateVecIn,
2524  ComplexArray stateVecOut)
2525 {
2526 
2527  long long int thisTask;
2528  const long long int numTasks=qureg.numAmpsPerChunk;
2529 
2530  qreal *stateVecRealIn=stateVecIn.real, *stateVecImagIn=stateVecIn.imag;
2531  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2532 
2533 # ifdef _OPENMP
2534 # pragma omp parallel \
2535  default (none) \
2536  shared (stateVecRealIn,stateVecImagIn,stateVecRealOut,stateVecImagOut) \
2537  private (thisTask)
2538 # endif
2539  {
2540 # ifdef _OPENMP
2541 # pragma omp for schedule (static)
2542 # endif
2543  for (thisTask=0; thisTask<numTasks; thisTask++) {
2544  stateVecRealOut[thisTask] = stateVecRealIn[thisTask];
2545  stateVecImagOut[thisTask] = stateVecImagIn[thisTask];
2546  }
2547  }
2548 }
2549 
2550 void statevec_controlledNotLocal(Qureg qureg, const int controlQubit, const int targetQubit)
2551 {
2552  long long int sizeBlock, sizeHalfBlock;
2553  long long int thisBlock, // current block
2554  indexUp,indexLo; // current index and corresponding index in lower half block
2555 
2556  qreal stateRealUp,stateImagUp;
2557  long long int thisTask;
2558  const long long int numTasks=qureg.numAmpsPerChunk>>1;
2559  const long long int chunkSize=qureg.numAmpsPerChunk;
2560  const long long int chunkId=qureg.chunkId;
2561 
2562  int controlBit;
2563 
2564  // set dimensions
2565  sizeHalfBlock = 1LL << targetQubit;
2566  sizeBlock = 2LL * sizeHalfBlock;
2567 
2568  // Can't use qureg.stateVec as a private OMP var
2569  qreal *stateVecReal = qureg.stateVec.real;
2570  qreal *stateVecImag = qureg.stateVec.imag;
2571 
2572 # ifdef _OPENMP
2573 # pragma omp parallel \
2574  default (none) \
2575  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag) \
2576  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp,controlBit)
2577 # endif
2578  {
2579 # ifdef _OPENMP
2580 # pragma omp for schedule (static)
2581 # endif
2582  for (thisTask=0; thisTask<numTasks; thisTask++) {
2583  thisBlock = thisTask / sizeHalfBlock;
2584  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2585  indexLo = indexUp + sizeHalfBlock;
2586 
2587  controlBit = extractBit(controlQubit, indexUp+chunkId*chunkSize);
2588  if (controlBit){
2589  stateRealUp = stateVecReal[indexUp];
2590  stateImagUp = stateVecImag[indexUp];
2591 
2592  stateVecReal[indexUp] = stateVecReal[indexLo];
2593  stateVecImag[indexUp] = stateVecImag[indexLo];
2594 
2595  stateVecReal[indexLo] = stateRealUp;
2596  stateVecImag[indexLo] = stateImagUp;
2597  }
2598  }
2599  }
2600 }
2601 
2612 void statevec_controlledNotDistributed (Qureg qureg, const int controlQubit,
2613  ComplexArray stateVecIn,
2614  ComplexArray stateVecOut)
2615 {
2616 
2617  long long int thisTask;
2618  const long long int numTasks=qureg.numAmpsPerChunk;
2619  const long long int chunkSize=qureg.numAmpsPerChunk;
2620  const long long int chunkId=qureg.chunkId;
2621 
2622  int controlBit;
2623 
2624  qreal *stateVecRealIn=stateVecIn.real, *stateVecImagIn=stateVecIn.imag;
2625  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2626 
2627 # ifdef _OPENMP
2628 # pragma omp parallel \
2629  default (none) \
2630  shared (stateVecRealIn,stateVecImagIn,stateVecRealOut,stateVecImagOut) \
2631  private (thisTask,controlBit)
2632 # endif
2633  {
2634 # ifdef _OPENMP
2635 # pragma omp for schedule (static)
2636 # endif
2637  for (thisTask=0; thisTask<numTasks; thisTask++) {
2638  controlBit = extractBit (controlQubit, thisTask+chunkId*chunkSize);
2639  if (controlBit){
2640  stateVecRealOut[thisTask] = stateVecRealIn[thisTask];
2641  stateVecImagOut[thisTask] = stateVecImagIn[thisTask];
2642  }
2643  }
2644  }
2645 }
2646 
2647 void statevec_pauliYLocal(Qureg qureg, const int targetQubit, const int conjFac)
2648 {
2649  long long int sizeBlock, sizeHalfBlock;
2650  long long int thisBlock, // current block
2651  indexUp,indexLo; // current index and corresponding index in lower half block
2652 
2653  qreal stateRealUp,stateImagUp;
2654  long long int thisTask;
2655  const long long int numTasks=qureg.numAmpsPerChunk>>1;
2656 
2657  // set dimensions
2658  sizeHalfBlock = 1LL << targetQubit;
2659  sizeBlock = 2LL * sizeHalfBlock;
2660 
2661  // Can't use qureg.stateVec as a private OMP var
2662  qreal *stateVecReal = qureg.stateVec.real;
2663  qreal *stateVecImag = qureg.stateVec.imag;
2664 
2665 # ifdef _OPENMP
2666 # pragma omp parallel \
2667  default (none) \
2668  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag) \
2669  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp)
2670 # endif
2671  {
2672 # ifdef _OPENMP
2673 # pragma omp for schedule (static)
2674 # endif
2675  for (thisTask=0; thisTask<numTasks; thisTask++) {
2676  thisBlock = thisTask / sizeHalfBlock;
2677  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2678  indexLo = indexUp + sizeHalfBlock;
2679 
2680  stateRealUp = stateVecReal[indexUp];
2681  stateImagUp = stateVecImag[indexUp];
2682 
2683  stateVecReal[indexUp] = conjFac * stateVecImag[indexLo];
2684  stateVecImag[indexUp] = conjFac * -stateVecReal[indexLo];
2685  stateVecReal[indexLo] = conjFac * -stateImagUp;
2686  stateVecImag[indexLo] = conjFac * stateRealUp;
2687  }
2688  }
2689 }
2690 
2704 void statevec_pauliYDistributed(Qureg qureg,
2705  ComplexArray stateVecIn,
2706  ComplexArray stateVecOut,
2707  int updateUpper, const int conjFac)
2708 {
2709 
2710  long long int thisTask;
2711  const long long int numTasks=qureg.numAmpsPerChunk;
2712 
2713  qreal *stateVecRealIn=stateVecIn.real, *stateVecImagIn=stateVecIn.imag;
2714  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2715 
2716  int realSign=1, imagSign=1;
2717  if (updateUpper) imagSign=-1;
2718  else realSign = -1;
2719 
2720 # ifdef _OPENMP
2721 # pragma omp parallel \
2722  default (none) \
2723  shared (stateVecRealIn,stateVecImagIn,stateVecRealOut,stateVecImagOut,realSign,imagSign) \
2724  private (thisTask)
2725 # endif
2726  {
2727 # ifdef _OPENMP
2728 # pragma omp for schedule (static)
2729 # endif
2730  for (thisTask=0; thisTask<numTasks; thisTask++) {
2731  stateVecRealOut[thisTask] = conjFac * realSign * stateVecImagIn[thisTask];
2732  stateVecImagOut[thisTask] = conjFac * imagSign * stateVecRealIn[thisTask];
2733  }
2734  }
2735 }
2736 
2737 
2738 
2739 
2740 void statevec_controlledPauliYLocal(Qureg qureg, const int controlQubit, const int targetQubit, const int conjFac)
2741 {
2742  long long int sizeBlock, sizeHalfBlock;
2743  long long int thisBlock, // current block
2744  indexUp,indexLo; // current index and corresponding index in lower half block
2745 
2746  qreal stateRealUp,stateImagUp;
2747  long long int thisTask;
2748  const long long int numTasks=qureg.numAmpsPerChunk>>1;
2749  const long long int chunkSize=qureg.numAmpsPerChunk;
2750  const long long int chunkId=qureg.chunkId;
2751 
2752  int controlBit;
2753 
2754  // set dimensions
2755  sizeHalfBlock = 1LL << targetQubit;
2756  sizeBlock = 2LL * sizeHalfBlock;
2757 
2758  // Can't use qureg.stateVec as a private OMP var
2759  qreal *stateVecReal = qureg.stateVec.real;
2760  qreal *stateVecImag = qureg.stateVec.imag;
2761 
2762 # ifdef _OPENMP
2763 # pragma omp parallel \
2764  default (none) \
2765  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag) \
2766  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp,controlBit)
2767 # endif
2768  {
2769 # ifdef _OPENMP
2770 # pragma omp for schedule (static)
2771 # endif
2772  for (thisTask=0; thisTask<numTasks; thisTask++) {
2773  thisBlock = thisTask / sizeHalfBlock;
2774  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2775  indexLo = indexUp + sizeHalfBlock;
2776 
2777  controlBit = extractBit(controlQubit, indexUp+chunkId*chunkSize);
2778  if (controlBit){
2779  stateRealUp = stateVecReal[indexUp];
2780  stateImagUp = stateVecImag[indexUp];
2781 
2782  // update under +-{{0, -i}, {i, 0}}
2783  stateVecReal[indexUp] = conjFac * stateVecImag[indexLo];
2784  stateVecImag[indexUp] = conjFac * -stateVecReal[indexLo];
2785  stateVecReal[indexLo] = conjFac * -stateImagUp;
2786  stateVecImag[indexLo] = conjFac * stateRealUp;
2787  }
2788  }
2789  }
2790 }
2791 
2792 
2793 void statevec_controlledPauliYDistributed (Qureg qureg, const int controlQubit,
2794  ComplexArray stateVecIn,
2795  ComplexArray stateVecOut, const int conjFac)
2796 {
2797 
2798  long long int thisTask;
2799  const long long int numTasks=qureg.numAmpsPerChunk;
2800  const long long int chunkSize=qureg.numAmpsPerChunk;
2801  const long long int chunkId=qureg.chunkId;
2802 
2803  int controlBit;
2804 
2805  qreal *stateVecRealIn=stateVecIn.real, *stateVecImagIn=stateVecIn.imag;
2806  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2807 
2808 # ifdef _OPENMP
2809 # pragma omp parallel \
2810  default (none) \
2811  shared (stateVecRealIn,stateVecImagIn,stateVecRealOut,stateVecImagOut) \
2812  private (thisTask,controlBit)
2813 # endif
2814  {
2815 # ifdef _OPENMP
2816 # pragma omp for schedule (static)
2817 # endif
2818  for (thisTask=0; thisTask<numTasks; thisTask++) {
2819  controlBit = extractBit (controlQubit, thisTask+chunkId*chunkSize);
2820  if (controlBit){
2821  stateVecRealOut[thisTask] = conjFac * stateVecImagIn[thisTask];
2822  stateVecImagOut[thisTask] = conjFac * -stateVecRealIn[thisTask];
2823  }
2824  }
2825  }
2826 }
2827 
2828 
2829 
2830 
2831 
2832 
2833 
2834 void statevec_hadamardLocal(Qureg qureg, const int targetQubit)
2835 {
2836  long long int sizeBlock, sizeHalfBlock;
2837  long long int thisBlock, // current block
2838  indexUp,indexLo; // current index and corresponding index in lower half block
2839 
2840  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2841  long long int thisTask;
2842  const long long int numTasks=qureg.numAmpsPerChunk>>1;
2843 
2844  // set dimensions
2845  sizeHalfBlock = 1LL << targetQubit;
2846  sizeBlock = 2LL * sizeHalfBlock;
2847 
2848  // Can't use qureg.stateVec as a private OMP var
2849  qreal *stateVecReal = qureg.stateVec.real;
2850  qreal *stateVecImag = qureg.stateVec.imag;
2851 
2852  qreal recRoot2 = 1.0/sqrt(2);
2853 
2854 # ifdef _OPENMP
2855 # pragma omp parallel \
2856  default (none) \
2857  shared (sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag, recRoot2) \
2858  private (thisTask,thisBlock ,indexUp,indexLo, stateRealUp,stateImagUp,stateRealLo,stateImagLo)
2859 # endif
2860  {
2861 # ifdef _OPENMP
2862 # pragma omp for schedule (static)
2863 # endif
2864  for (thisTask=0; thisTask<numTasks; thisTask++) {
2865  thisBlock = thisTask / sizeHalfBlock;
2866  indexUp = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
2867  indexLo = indexUp + sizeHalfBlock;
2868 
2869  stateRealUp = stateVecReal[indexUp];
2870  stateImagUp = stateVecImag[indexUp];
2871 
2872  stateRealLo = stateVecReal[indexLo];
2873  stateImagLo = stateVecImag[indexLo];
2874 
2875  stateVecReal[indexUp] = recRoot2*(stateRealUp + stateRealLo);
2876  stateVecImag[indexUp] = recRoot2*(stateImagUp + stateImagLo);
2877 
2878  stateVecReal[indexLo] = recRoot2*(stateRealUp - stateRealLo);
2879  stateVecImag[indexLo] = recRoot2*(stateImagUp - stateImagLo);
2880  }
2881  }
2882 }
2883 
2894 void statevec_hadamardDistributed(Qureg qureg,
2895  ComplexArray stateVecUp,
2896  ComplexArray stateVecLo,
2897  ComplexArray stateVecOut,
2898  int updateUpper)
2899 {
2900 
2901  qreal stateRealUp,stateRealLo,stateImagUp,stateImagLo;
2902  long long int thisTask;
2903  const long long int numTasks=qureg.numAmpsPerChunk;
2904 
2905  int sign;
2906  if (updateUpper) sign=1;
2907  else sign=-1;
2908 
2909  qreal recRoot2 = 1.0/sqrt(2);
2910 
2911  qreal *stateVecRealUp=stateVecUp.real, *stateVecImagUp=stateVecUp.imag;
2912  qreal *stateVecRealLo=stateVecLo.real, *stateVecImagLo=stateVecLo.imag;
2913  qreal *stateVecRealOut=stateVecOut.real, *stateVecImagOut=stateVecOut.imag;
2914 
2915 # ifdef _OPENMP
2916 # pragma omp parallel \
2917  default (none) \
2918  shared (stateVecRealUp,stateVecImagUp,stateVecRealLo,stateVecImagLo,stateVecRealOut,stateVecImagOut, \
2919  recRoot2, sign) \
2920  private (thisTask,stateRealUp,stateImagUp,stateRealLo,stateImagLo)
2921 # endif
2922  {
2923 # ifdef _OPENMP
2924 # pragma omp for schedule (static)
2925 # endif
2926  for (thisTask=0; thisTask<numTasks; thisTask++) {
2927  // store current state vector values in temp variables
2928  stateRealUp = stateVecRealUp[thisTask];
2929  stateImagUp = stateVecImagUp[thisTask];
2930 
2931  stateRealLo = stateVecRealLo[thisTask];
2932  stateImagLo = stateVecImagLo[thisTask];
2933 
2934  stateVecRealOut[thisTask] = recRoot2*(stateRealUp + sign*stateRealLo);
2935  stateVecImagOut[thisTask] = recRoot2*(stateImagUp + sign*stateImagLo);
2936  }
2937  }
2938 }
2939 
2940 void statevec_phaseShiftByTerm (Qureg qureg, const int targetQubit, Complex term)
2941 {
2942  long long int index;
2943  long long int stateVecSize;
2944  int targetBit;
2945 
2946  const long long int chunkSize=qureg.numAmpsPerChunk;
2947  const long long int chunkId=qureg.chunkId;
2948 
2949  // dimension of the state vector
2950  stateVecSize = qureg.numAmpsPerChunk;
2951  qreal *stateVecReal = qureg.stateVec.real;
2952  qreal *stateVecImag = qureg.stateVec.imag;
2953 
2954  qreal stateRealLo, stateImagLo;
2955  const qreal cosAngle = term.real;
2956  const qreal sinAngle = term.imag;
2957 
2958 # ifdef _OPENMP
2959 # pragma omp parallel for \
2960  default (none) \
2961  shared (stateVecSize, stateVecReal,stateVecImag ) \
2962  private (index,targetBit,stateRealLo,stateImagLo) \
2963  schedule (static)
2964 # endif
2965  for (index=0; index<stateVecSize; index++) {
2966 
2967  // update the coeff of the |1> state of the target qubit
2968  targetBit = extractBit (targetQubit, index+chunkId*chunkSize);
2969  if (targetBit) {
2970 
2971  stateRealLo = stateVecReal[index];
2972  stateImagLo = stateVecImag[index];
2973 
2974  stateVecReal[index] = cosAngle*stateRealLo - sinAngle*stateImagLo;
2975  stateVecImag[index] = sinAngle*stateRealLo + cosAngle*stateImagLo;
2976  }
2977  }
2978 }
2979 
2980 void statevec_controlledPhaseShift (Qureg qureg, const int idQubit1, const int idQubit2, qreal angle)
2981 {
2982  long long int index;
2983  long long int stateVecSize;
2984  int bit1, bit2;
2985 
2986  const long long int chunkSize=qureg.numAmpsPerChunk;
2987  const long long int chunkId=qureg.chunkId;
2988 
2989  // dimension of the state vector
2990  stateVecSize = qureg.numAmpsPerChunk;
2991  qreal *stateVecReal = qureg.stateVec.real;
2992  qreal *stateVecImag = qureg.stateVec.imag;
2993 
2994  qreal stateRealLo, stateImagLo;
2995  const qreal cosAngle = cos(angle);
2996  const qreal sinAngle = sin(angle);
2997 
2998 # ifdef _OPENMP
2999 # pragma omp parallel for \
3000  default (none) \
3001  shared (stateVecSize, stateVecReal,stateVecImag ) \
3002  private (index,bit1,bit2,stateRealLo,stateImagLo) \
3003  schedule (static)
3004 # endif
3005  for (index=0; index<stateVecSize; index++) {
3006  bit1 = extractBit (idQubit1, index+chunkId*chunkSize);
3007  bit2 = extractBit (idQubit2, index+chunkId*chunkSize);
3008  if (bit1 && bit2) {
3009 
3010  stateRealLo = stateVecReal[index];
3011  stateImagLo = stateVecImag[index];
3012 
3013  stateVecReal[index] = cosAngle*stateRealLo - sinAngle*stateImagLo;
3014  stateVecImag[index] = sinAngle*stateRealLo + cosAngle*stateImagLo;
3015  }
3016  }
3017 }
3018 
3019 void statevec_multiControlledPhaseShift(Qureg qureg, int *controlQubits, int numControlQubits, qreal angle)
3020 {
3021  long long int index;
3022  long long int stateVecSize;
3023 
3024  const long long int chunkSize=qureg.numAmpsPerChunk;
3025  const long long int chunkId=qureg.chunkId;
3026 
3027  long long int mask = getQubitBitMask(controlQubits, numControlQubits);
3028 
3029  stateVecSize = qureg.numAmpsPerChunk;
3030  qreal *stateVecReal = qureg.stateVec.real;
3031  qreal *stateVecImag = qureg.stateVec.imag;
3032 
3033  qreal stateRealLo, stateImagLo;
3034  const qreal cosAngle = cos(angle);
3035  const qreal sinAngle = sin(angle);
3036 
3037 # ifdef _OPENMP
3038 # pragma omp parallel \
3039  default (none) \
3040  shared (stateVecSize, stateVecReal, stateVecImag, mask) \
3041  private (index, stateRealLo, stateImagLo)
3042 # endif
3043  {
3044 # ifdef _OPENMP
3045 # pragma omp for schedule (static)
3046 # endif
3047  for (index=0; index<stateVecSize; index++) {
3048  if (mask == (mask & (index+chunkId*chunkSize)) ){
3049 
3050  stateRealLo = stateVecReal[index];
3051  stateImagLo = stateVecImag[index];
3052 
3053  stateVecReal[index] = cosAngle*stateRealLo - sinAngle*stateImagLo;
3054  stateVecImag[index] = sinAngle*stateRealLo + cosAngle*stateImagLo;
3055  }
3056  }
3057  }
3058 }
3059 
3060 int getBitMaskParity(long long int mask) {
3061  int parity = 0;
3062  while (mask) {
3063  parity = !parity;
3064  mask = mask & (mask-1);
3065  }
3066  return parity;
3067 }
3068 
3069 void statevec_multiRotateZ(Qureg qureg, long long int mask, qreal angle)
3070 {
3071  long long int index;
3072  long long int stateVecSize;
3073 
3074  const long long int chunkSize=qureg.numAmpsPerChunk;
3075  const long long int chunkId=qureg.chunkId;
3076 
3077  stateVecSize = qureg.numAmpsPerChunk;
3078  qreal *stateVecReal = qureg.stateVec.real;
3079  qreal *stateVecImag = qureg.stateVec.imag;
3080 
3081  qreal stateReal, stateImag;
3082  const qreal cosAngle = cos(angle/2.0);
3083  const qreal sinAngle = sin(angle/2.0);
3084 
3085  // = +-1, to flip sinAngle based on target qubit parity, to effect
3086  // exp(-angle/2 i fac_j)|j>
3087  int fac;
3088 
3089 # ifdef _OPENMP
3090 # pragma omp parallel \
3091  default (none) \
3092  shared (stateVecSize, stateVecReal, stateVecImag, mask) \
3093  private (index, fac, stateReal, stateImag)
3094 # endif
3095  {
3096 # ifdef _OPENMP
3097 # pragma omp for schedule (static)
3098 # endif
3099  for (index=0; index<stateVecSize; index++) {
3100  stateReal = stateVecReal[index];
3101  stateImag = stateVecImag[index];
3102 
3103  // odd-parity target qubits get fac_j = -1
3104  fac = getBitMaskParity(mask & (index+chunkId*chunkSize))? -1 : 1;
3105  stateVecReal[index] = cosAngle*stateReal + fac * sinAngle*stateImag;
3106  stateVecImag[index] = - fac * sinAngle*stateReal + cosAngle*stateImag;
3107  }
3108  }
3109 }
3110 
3111 qreal densmatr_findProbabilityOfZeroLocal(Qureg qureg, const int measureQubit) {
3112 
3113  // computes first local index containing a diagonal element
3114  long long int localNumAmps = qureg.numAmpsPerChunk;
3115  long long int densityDim = (1LL << qureg.numQubitsRepresented);
3116  long long int diagSpacing = 1LL + densityDim;
3117  long long int maxNumDiagsPerChunk = 1 + localNumAmps / diagSpacing;
3118  long long int numPrevDiags = (qureg.chunkId>0)? 1+(qureg.chunkId*localNumAmps)/diagSpacing : 0;
3119  long long int globalIndNextDiag = diagSpacing * numPrevDiags;
3120  long long int localIndNextDiag = globalIndNextDiag % localNumAmps;
3121 
3122  // computes how many diagonals are contained in this chunk
3123  long long int numDiagsInThisChunk = maxNumDiagsPerChunk;
3124  if (localIndNextDiag + (numDiagsInThisChunk-1)*diagSpacing >= localNumAmps)
3125  numDiagsInThisChunk -= 1;
3126 
3127  long long int visitedDiags; // number of visited diagonals in this chunk so far
3128  long long int basisStateInd; // current diagonal index being considered
3129  long long int index; // index in the local chunk
3130 
3131  qreal zeroProb = 0;
3132  qreal *stateVecReal = qureg.stateVec.real;
3133 
3134 # ifdef _OPENMP
3135 # pragma omp parallel \
3136  shared (localIndNextDiag, numPrevDiags, diagSpacing, stateVecReal, numDiagsInThisChunk) \
3137  private (visitedDiags, basisStateInd, index) \
3138  reduction ( +:zeroProb )
3139 # endif
3140  {
3141 # ifdef _OPENMP
3142 # pragma omp for schedule (static)
3143 # endif
3144  // sums the diagonal elems of the density matrix where measureQubit=0
3145  for (visitedDiags = 0; visitedDiags < numDiagsInThisChunk; visitedDiags++) {
3146 
3147  basisStateInd = numPrevDiags + visitedDiags;
3148  index = localIndNextDiag + diagSpacing * visitedDiags;
3149 
3150  if (extractBit(measureQubit, basisStateInd) == 0)
3151  zeroProb += stateVecReal[index]; // assume imag[diagonls] ~ 0
3152 
3153  }
3154  }
3155 
3156  return zeroProb;
3157 }
3158 
3166 qreal statevec_findProbabilityOfZeroLocal (Qureg qureg,
3167  const int measureQubit)
3168 {
3169  // ----- sizes
3170  long long int sizeBlock, // size of blocks
3171  sizeHalfBlock; // size of blocks halved
3172  // ----- indices
3173  long long int thisBlock, // current block
3174  index; // current index for first half block
3175  // ----- measured probability
3176  qreal totalProbability; // probability (returned) value
3177  // ----- temp variables
3178  long long int thisTask;
3179  long long int numTasks=qureg.numAmpsPerChunk>>1;
3180 
3181  // ---------------------------------------------------------------- //
3182  // dimensions //
3183  // ---------------------------------------------------------------- //
3184  sizeHalfBlock = 1LL << (measureQubit); // number of state vector elements to sum,
3185  // and then the number to skip
3186  sizeBlock = 2LL * sizeHalfBlock; // size of blocks (pairs of measure and skip entries)
3187 
3188  // initialise returned value
3189  totalProbability = 0.0;
3190 
3191  qreal *stateVecReal = qureg.stateVec.real;
3192  qreal *stateVecImag = qureg.stateVec.imag;
3193 
3194 # ifdef _OPENMP
3195 # pragma omp parallel \
3196  shared (numTasks,sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag) \
3197  private (thisTask,thisBlock,index) \
3198  reduction ( +:totalProbability )
3199 # endif
3200  {
3201 # ifdef _OPENMP
3202 # pragma omp for schedule (static)
3203 # endif
3204  for (thisTask=0; thisTask<numTasks; thisTask++) {
3205  thisBlock = thisTask / sizeHalfBlock;
3206  index = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
3207 
3208  totalProbability += stateVecReal[index]*stateVecReal[index]
3209  + stateVecImag[index]*stateVecImag[index];
3210  }
3211  }
3212  return totalProbability;
3213 }
3214 
3222 qreal statevec_findProbabilityOfZeroDistributed (Qureg qureg) {
3223  // ----- measured probability
3224  qreal totalProbability; // probability (returned) value
3225  // ----- temp variables
3226  long long int thisTask; // task based approach for expose loop with small granularity
3227  long long int numTasks=qureg.numAmpsPerChunk;
3228 
3229  // ---------------------------------------------------------------- //
3230  // find probability //
3231  // ---------------------------------------------------------------- //
3232 
3233  // initialise returned value
3234  totalProbability = 0.0;
3235 
3236  qreal *stateVecReal = qureg.stateVec.real;
3237  qreal *stateVecImag = qureg.stateVec.imag;
3238 
3239 # ifdef _OPENMP
3240 # pragma omp parallel \
3241  shared (numTasks,stateVecReal,stateVecImag) \
3242  private (thisTask) \
3243  reduction ( +:totalProbability )
3244 # endif
3245  {
3246 # ifdef _OPENMP
3247 # pragma omp for schedule (static)
3248 # endif
3249  for (thisTask=0; thisTask<numTasks; thisTask++) {
3250  totalProbability += stateVecReal[thisTask]*stateVecReal[thisTask]
3251  + stateVecImag[thisTask]*stateVecImag[thisTask];
3252  }
3253  }
3254 
3255  return totalProbability;
3256 }
3257 
3258 
3259 
3260 void statevec_controlledPhaseFlip (Qureg qureg, const int idQubit1, const int idQubit2)
3261 {
3262  long long int index;
3263  long long int stateVecSize;
3264  int bit1, bit2;
3265 
3266  const long long int chunkSize=qureg.numAmpsPerChunk;
3267  const long long int chunkId=qureg.chunkId;
3268 
3269  // dimension of the state vector
3270  stateVecSize = qureg.numAmpsPerChunk;
3271  qreal *stateVecReal = qureg.stateVec.real;
3272  qreal *stateVecImag = qureg.stateVec.imag;
3273 
3274 # ifdef _OPENMP
3275 # pragma omp parallel for \
3276  default (none) \
3277  shared (stateVecSize, stateVecReal,stateVecImag ) \
3278  private (index,bit1,bit2) \
3279  schedule (static)
3280 # endif
3281  for (index=0; index<stateVecSize; index++) {
3282  bit1 = extractBit (idQubit1, index+chunkId*chunkSize);
3283  bit2 = extractBit (idQubit2, index+chunkId*chunkSize);
3284  if (bit1 && bit2) {
3285  stateVecReal [index] = - stateVecReal [index];
3286  stateVecImag [index] = - stateVecImag [index];
3287  }
3288  }
3289 }
3290 
3291 void statevec_multiControlledPhaseFlip(Qureg qureg, int *controlQubits, int numControlQubits)
3292 {
3293  long long int index;
3294  long long int stateVecSize;
3295 
3296  const long long int chunkSize=qureg.numAmpsPerChunk;
3297  const long long int chunkId=qureg.chunkId;
3298 
3299  long long int mask = getQubitBitMask(controlQubits, numControlQubits);
3300 
3301  stateVecSize = qureg.numAmpsPerChunk;
3302  qreal *stateVecReal = qureg.stateVec.real;
3303  qreal *stateVecImag = qureg.stateVec.imag;
3304 
3305 # ifdef _OPENMP
3306 # pragma omp parallel \
3307  default (none) \
3308  shared (stateVecSize, stateVecReal,stateVecImag, mask ) \
3309  private (index)
3310 # endif
3311  {
3312 # ifdef _OPENMP
3313 # pragma omp for schedule (static)
3314 # endif
3315  for (index=0; index<stateVecSize; index++) {
3316  if (mask == (mask & (index+chunkId*chunkSize)) ){
3317  stateVecReal [index] = - stateVecReal [index];
3318  stateVecImag [index] = - stateVecImag [index];
3319  }
3320  }
3321  }
3322 }
3323 
3340 void statevec_collapseToKnownProbOutcomeLocal(Qureg qureg, int measureQubit, int outcome, qreal totalProbability)
3341 {
3342  // ----- sizes
3343  long long int sizeBlock, // size of blocks
3344  sizeHalfBlock; // size of blocks halved
3345  // ----- indices
3346  long long int thisBlock, // current block
3347  index; // current index for first half block
3348  // ----- measured probability
3349  qreal renorm; // probability (returned) value
3350  // ----- temp variables
3351  long long int thisTask; // task based approach for expose loop with small granularity
3352  // (good for shared memory parallelism)
3353  long long int numTasks=qureg.numAmpsPerChunk>>1;
3354 
3355  // ---------------------------------------------------------------- //
3356  // dimensions //
3357  // ---------------------------------------------------------------- //
3358  sizeHalfBlock = 1LL << (measureQubit); // number of state vector elements to sum,
3359  // and then the number to skip
3360  sizeBlock = 2LL * sizeHalfBlock; // size of blocks (pairs of measure and skip entries)
3361 
3362  renorm=1/sqrt(totalProbability);
3363  qreal *stateVecReal = qureg.stateVec.real;
3364  qreal *stateVecImag = qureg.stateVec.imag;
3365 
3366 
3367 # ifdef _OPENMP
3368 # pragma omp parallel \
3369  default (none) \
3370  shared (numTasks,sizeBlock,sizeHalfBlock, stateVecReal,stateVecImag,renorm,outcome) \
3371  private (thisTask,thisBlock,index)
3372 # endif
3373  {
3374  if (outcome==0){
3375  // measure qubit is 0
3376 # ifdef _OPENMP
3377 # pragma omp for schedule (static)
3378 # endif
3379  for (thisTask=0; thisTask<numTasks; thisTask++) {
3380  thisBlock = thisTask / sizeHalfBlock;
3381  index = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
3382  stateVecReal[index]=stateVecReal[index]*renorm;
3383  stateVecImag[index]=stateVecImag[index]*renorm;
3384 
3385  stateVecReal[index+sizeHalfBlock]=0;
3386  stateVecImag[index+sizeHalfBlock]=0;
3387  }
3388  } else {
3389  // measure qubit is 1
3390 # ifdef _OPENMP
3391 # pragma omp for schedule (static)
3392 # endif
3393  for (thisTask=0; thisTask<numTasks; thisTask++) {
3394  thisBlock = thisTask / sizeHalfBlock;
3395  index = thisBlock*sizeBlock + thisTask%sizeHalfBlock;
3396  stateVecReal[index]=0;
3397  stateVecImag[index]=0;
3398 
3399  stateVecReal[index+sizeHalfBlock]=stateVecReal[index+sizeHalfBlock]*renorm;
3400  stateVecImag[index+sizeHalfBlock]=stateVecImag[index+sizeHalfBlock]*renorm;
3401  }
3402  }
3403  }
3404 
3405 }
3406 
3422 void statevec_collapseToKnownProbOutcomeDistributedRenorm (Qureg qureg, const int measureQubit, const qreal totalProbability)
3423 {
3424  // ----- temp variables
3425  long long int thisTask;
3426  long long int numTasks=qureg.numAmpsPerChunk;
3427 
3428  qreal renorm=1/sqrt(totalProbability);
3429 
3430  qreal *stateVecReal = qureg.stateVec.real;
3431  qreal *stateVecImag = qureg.stateVec.imag;
3432 
3433 # ifdef _OPENMP
3434 # pragma omp parallel \
3435  shared (numTasks,stateVecReal,stateVecImag) \
3436  private (thisTask)
3437 # endif
3438  {
3439 # ifdef _OPENMP
3440 # pragma omp for schedule (static)
3441 # endif
3442  for (thisTask=0; thisTask<numTasks; thisTask++) {
3443  stateVecReal[thisTask] = stateVecReal[thisTask]*renorm;
3444  stateVecImag[thisTask] = stateVecImag[thisTask]*renorm;
3445  }
3446  }
3447 }
3448 
3461 void statevec_collapseToOutcomeDistributedSetZero(Qureg qureg)
3462 {
3463  // ----- temp variables
3464  long long int thisTask;
3465  long long int numTasks=qureg.numAmpsPerChunk;
3466 
3467  // ---------------------------------------------------------------- //
3468  // find probability //
3469  // ---------------------------------------------------------------- //
3470 
3471  qreal *stateVecReal = qureg.stateVec.real;
3472  qreal *stateVecImag = qureg.stateVec.imag;
3473 
3474 # ifdef _OPENMP
3475 # pragma omp parallel \
3476  shared (numTasks,stateVecReal,stateVecImag) \
3477  private (thisTask)
3478 # endif
3479  {
3480 # ifdef _OPENMP
3481 # pragma omp for schedule (static)
3482 # endif
3483  for (thisTask=0; thisTask<numTasks; thisTask++) {
3484  stateVecReal[thisTask] = 0;
3485  stateVecImag[thisTask] = 0;
3486  }
3487  }
3488 }
3489 
3496 void statevec_swapQubitAmpsLocal(Qureg qureg, int qb1, int qb2) {
3497 
3498  // can't use qureg.stateVec as a private OMP var
3499  qreal *reVec = qureg.stateVec.real;
3500  qreal *imVec = qureg.stateVec.imag;
3501 
3502  long long int numTasks = qureg.numAmpsPerChunk >> 2; // each iteration updates 2 amps and skips 2 amps
3503  long long int thisTask;
3504  long long int ind00, ind01, ind10;
3505  qreal re01, re10;
3506  qreal im01, im10;
3507 
3508 # ifdef _OPENMP
3509 # pragma omp parallel \
3510  default (none) \
3511  shared (reVec,imVec,numTasks,qb1,qb2) \
3512  private (thisTask, ind00,ind01,ind10, re01,re10, im01,im10)
3513 # endif
3514  {
3515 # ifdef _OPENMP
3516 # pragma omp for schedule (static)
3517 # endif
3518  for (thisTask=0; thisTask<numTasks; thisTask++) {
3519  // determine ind00 of |..0..0..>, |..0..1..> and |..1..0..>
3520  ind00 = insertTwoZeroBits(thisTask, qb1, qb2);
3521  ind01 = flipBit(ind00, qb1);
3522  ind10 = flipBit(ind00, qb2);
3523 
3524  // extract statevec amplitudes
3525  re01 = reVec[ind01]; im01 = imVec[ind01];
3526  re10 = reVec[ind10]; im10 = imVec[ind10];
3527 
3528  // swap 01 and 10 amps
3529  reVec[ind01] = re10; reVec[ind10] = re01;
3530  imVec[ind01] = im10; imVec[ind10] = im01;
3531  }
3532  }
3533 }
3534 
3539 void statevec_swapQubitAmpsDistributed(Qureg qureg, int pairRank, int qb1, int qb2) {
3540 
3541  // can't use qureg.stateVec as a private OMP var
3542  qreal *reVec = qureg.stateVec.real;
3543  qreal *imVec = qureg.stateVec.imag;
3544  qreal *rePairVec = qureg.pairStateVec.real;
3545  qreal *imPairVec = qureg.pairStateVec.imag;
3546 
3547  long long int numLocalAmps = qureg.numAmpsPerChunk;
3548  long long int globalStartInd = qureg.chunkId * numLocalAmps;
3549  long long int pairGlobalStartInd = pairRank * numLocalAmps;
3550 
3551  long long int localInd, globalInd;
3552  long long int pairLocalInd, pairGlobalInd;
3553 
3554 # ifdef _OPENMP
3555 # pragma omp parallel \
3556  default (none) \
3557  shared (reVec,imVec,rePairVec,imPairVec,numLocalAmps,globalStartInd,pairGlobalStartInd,qb1,qb2) \
3558  private (localInd,globalInd, pairLocalInd,pairGlobalInd)
3559 # endif
3560  {
3561 # ifdef _OPENMP
3562 # pragma omp for schedule (static)
3563 # endif
3564  for (localInd=0; localInd < numLocalAmps; localInd++) {
3565 
3566  globalInd = globalStartInd + localInd;
3567  if (isOddParity(globalInd, qb1, qb2)) {
3568 
3569  pairGlobalInd = flipBit(flipBit(globalInd, qb1), qb2);
3570  pairLocalInd = pairGlobalInd - pairGlobalStartInd;
3571 
3572  reVec[localInd] = rePairVec[pairLocalInd];
3573  imVec[localInd] = imPairVec[pairLocalInd];
3574  }
3575  }
3576  }
3577 }
3578 
3579 void statevec_setWeightedQureg(Complex fac1, Qureg qureg1, Complex fac2, Qureg qureg2, Complex facOut, Qureg out) {
3580 
3581  long long int numAmps = qureg1.numAmpsPerChunk;
3582 
3583  qreal *vecRe1 = qureg1.stateVec.real;
3584  qreal *vecIm1 = qureg1.stateVec.imag;
3585  qreal *vecRe2 = qureg2.stateVec.real;
3586  qreal *vecIm2 = qureg2.stateVec.imag;
3587  qreal *vecReOut = out.stateVec.real;
3588  qreal *vecImOut = out.stateVec.imag;
3589 
3590  qreal facRe1 = fac1.real;
3591  qreal facIm1 = fac1.imag;
3592  qreal facRe2 = fac2.real;
3593  qreal facIm2 = fac2.imag;
3594  qreal facReOut = facOut.real;
3595  qreal facImOut = facOut.imag;
3596 
3597  qreal re1,im1, re2,im2, reOut,imOut;
3598  long long int index;
3599 
3600 # ifdef _OPENMP
3601 # pragma omp parallel \
3602  shared (vecRe1,vecIm1, vecRe2,vecIm2, vecReOut,vecImOut, facRe1,facIm1,facRe2,facIm2, numAmps) \
3603  private (index, re1,im1, re2,im2, reOut,imOut)
3604 # endif
3605  {
3606 # ifdef _OPENMP
3607 # pragma omp for schedule (static)
3608 # endif
3609  for (index=0LL; index<numAmps; index++) {
3610  re1 = vecRe1[index]; im1 = vecIm1[index];
3611  re2 = vecRe2[index]; im2 = vecIm2[index];
3612  reOut = vecReOut[index];
3613  imOut = vecImOut[index];
3614 
3615  vecReOut[index] = (facReOut*reOut - facImOut*imOut) + (facRe1*re1 - facIm1*im1) + (facRe2*re2 - facIm2*im2);
3616  vecImOut[index] = (facReOut*imOut + facImOut*reOut) + (facRe1*im1 + facIm1*re1) + (facRe2*im2 + facIm2*re2);
3617  }
3618  }
3619 }
3620 
statevec_phaseShiftByTerm
void statevec_phaseShiftByTerm(Qureg qureg, const int targetQubit, Complex term)
Definition: QuEST_cpu.c:2940
qsortComp
int qsortComp(const void *a, const void *b)
Definition: QuEST_cpu.c:1810
copyStateFromGPU
void copyStateFromGPU(Qureg qureg)
In GPU mode, this copies the state-vector (or density matrix) from GPU memory (qureg....
Definition: QuEST_cpu.c:39
ComplexMatrix4::real
qreal real[4][4]
Definition: QuEST.h:127
syncQuESTEnv
void syncQuESTEnv(QuESTEnv env)
Guarantees that all code up to the given point has been executed on all nodes (if running in distribu...
Definition: QuEST_cpu_distributed.c:162
statevec_hadamardLocal
void statevec_hadamardLocal(Qureg qureg, const int targetQubit)
Definition: QuEST_cpu.c:2834
densmatr_calcHilbertSchmidtDistanceSquaredLocal
qreal densmatr_calcHilbertSchmidtDistanceSquaredLocal(Qureg a, Qureg b)
computes Tr((a-b) conjTrans(a-b)) = sum of abs values of (a-b)
Definition: QuEST_cpu.c:922
densmatr_mixTwoQubitDepolarisingDistributed
void densmatr_mixTwoQubitDepolarisingDistributed(Qureg qureg, const int targetQubit, const int qubit2, qreal delta, qreal gamma)
Definition: QuEST_cpu.c:540
QuESTEnv::rank
int rank
Definition: QuEST.h:201
statevec_controlledNotDistributed
void statevec_controlledNotDistributed(Qureg qureg, const int controlQubit, ComplexArray stateVecIn, ComplexArray stateVecOut)
Rotate a single qubit by {{0,1},{1,0}.
Definition: QuEST_cpu.c:2612
statevec_controlledPhaseShift
void statevec_controlledPhaseShift(Qureg qureg, const int idQubit1, const int idQubit2, qreal angle)
Definition: QuEST_cpu.c:2980
densmatr_initClassicalState
void densmatr_initClassicalState(Qureg qureg, long long int stateInd)
Definition: QuEST_cpu.c:1114
statevec_controlledPhaseFlip
void statevec_controlledPhaseFlip(Qureg qureg, const int idQubit1, const int idQubit2)
Definition: QuEST_cpu.c:3260
Qureg::pairStateVec
ComplexArray pairStateVec
Temporary storage for a chunk of the state vector received from another process in the MPI version.
Definition: QuEST.h:181
statevec_initDebugState
void statevec_initDebugState(Qureg qureg)
Initialise the state vector of probability amplitudes to an (unphysical) state with each component of...
Definition: QuEST_cpu.c:1559
QuEST_internal.h
statevec_findProbabilityOfZeroDistributed
qreal statevec_findProbabilityOfZeroDistributed(Qureg qureg)
Measure the probability of a specified qubit being in the zero state across all amplitudes held in th...
Definition: QuEST_cpu.c:3222
statevec_multiControlledTwoQubitUnitaryLocal
void statevec_multiControlledTwoQubitUnitaryLocal(Qureg qureg, long long int ctrlMask, const int q1, const int q2, ComplexMatrix4 u)
Definition: QuEST_cpu.c:1715
densmatr_mixTwoQubitDepolarisingQ1LocalQ2DistributedPart3
void densmatr_mixTwoQubitDepolarisingQ1LocalQ2DistributedPart3(Qureg qureg, const int targetQubit, const int qubit2, qreal delta, qreal gamma)
Definition: QuEST_cpu.c:631
Qureg::numChunks
int numChunks
Number of chunks the state vector is broken up into – the number of MPI processes used.
Definition: QuEST.h:176
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.
Definition: QuEST_cpu.c:3496
getBitMaskParity
int getBitMaskParity(long long int mask)
Definition: QuEST_cpu.c:3060
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=...
Definition: QuEST_cpu.c:3340
statevec_cloneQureg
void statevec_cloneQureg(Qureg targetQureg, Qureg copyQureg)
Definition: QuEST_cpu.c:1474
extractBit
__forceinline__ __device__ int extractBit(int locationOfBitFromRight, long long int theEncodedNumber)
Definition: QuEST_gpu.cu:82
statevec_findProbabilityOfZeroLocal
qreal statevec_findProbabilityOfZeroLocal(Qureg qureg, const int measureQubit)
Measure the total probability of a specified qubit being in the zero state across all amplitudes in t...
Definition: QuEST_cpu.c:3166
statevec_setAmps
void statevec_setAmps(Qureg qureg, long long int startInd, qreal *reals, qreal *imags, long long int numAmps)
Definition: QuEST_cpu.c:1236
densmatr_calcPurityLocal
qreal densmatr_calcPurityLocal(Qureg qureg)
Definition: QuEST_cpu.c:860
ComplexMatrix4
Represents a 4x4 matrix of complex numbers.
Definition: QuEST.h:125
QuESTEnv
Information about the environment the program is running in.
Definition: QuEST.h:199
statevec_initBlankState
void statevec_initBlankState(Qureg qureg)
Definition: QuEST_cpu.c:1366
ComplexMatrixN
Represents a general 2^N by 2^N matrix of complex numbers.
Definition: QuEST.h:136
statevec_initClassicalState
void statevec_initClassicalState(Qureg qureg, long long int stateInd)
Definition: QuEST_cpu.c:1438
statevec_pauliXDistributed
void statevec_pauliXDistributed(Qureg qureg, ComplexArray stateVecIn, ComplexArray stateVecOut)
Rotate a single qubit by {{0,1},{1,0}.
Definition: QuEST_cpu.c:2522
qreal
#define qreal
statevec_multiControlledPhaseShift
void statevec_multiControlledPhaseShift(Qureg qureg, int *controlQubits, int numControlQubits, qreal angle)
Definition: QuEST_cpu.c:3019
Qureg::numQubitsInStateVec
int numQubitsInStateVec
Number of qubits in the state-vector - this is double the number represented for mixed states.
Definition: QuEST.h:167
statevec_controlledNotLocal
void statevec_controlledNotLocal(Qureg qureg, const int controlQubit, const int targetQubit)
Definition: QuEST_cpu.c:2550
insertTwoZeroBits
__forceinline__ __device__ long long int insertTwoZeroBits(long long int number, int bit1, int bit2)
Definition: QuEST_gpu.cu:106
getQubitBitMask
long long int getQubitBitMask(int *qubits, const int numQubits)
Definition: QuEST_common.c:43
statevec_collapseToOutcomeDistributedSetZero
void statevec_collapseToOutcomeDistributedSetZero(Qureg qureg)
Set all amplitudes in one chunk to 0.
Definition: QuEST_cpu.c:3461
Qureg::chunkId
int chunkId
The position of the chunk of the state vector held by this process in the full state vector.
Definition: QuEST.h:174
densmatr_mixDampingLocal
void densmatr_mixDampingLocal(Qureg qureg, const int targetQubit, qreal damping)
Definition: QuEST_cpu.c:174
densmatr_mixTwoQubitDephasing
void densmatr_mixTwoQubitDephasing(Qureg qureg, const int qubit1, const int qubit2, qreal dephase)
Definition: QuEST_cpu.c:84
ComplexMatrix2::imag
qreal imag[2][2]
Definition: QuEST.h:117
statevec_multiControlledUnitaryLocal
void statevec_multiControlledUnitaryLocal(Qureg qureg, const int targetQubit, long long int ctrlQubitsMask, long long int ctrlFlipMask, ComplexMatrix2 u)
Definition: QuEST_cpu.c:2140
QuEST_cpu_internal.h
statevec_controlledCompactUnitaryLocal
void statevec_controlledCompactUnitaryLocal(Qureg qureg, const int controlQubit, const int targetQubit, Complex alpha, Complex beta)
Definition: QuEST_cpu.c:2069
Qureg::numAmpsPerChunk
long long int numAmpsPerChunk
Number of probability amplitudes held in stateVec by this process In the non-MPI version,...
Definition: QuEST.h:170
densmatr_mixTwoQubitDepolarisingLocal
void densmatr_mixTwoQubitDepolarisingLocal(Qureg qureg, int qubit1, int qubit2, qreal delta, qreal gamma)
Definition: QuEST_cpu.c:385
statevec_initZeroState
void statevec_initZeroState(Qureg qureg)
Definition: QuEST_cpu.c:1396
statevec_initPlusState
void statevec_initPlusState(Qureg qureg)
Definition: QuEST_cpu.c:1406
statevec_createQureg
void statevec_createQureg(Qureg *qureg, int numQubits, QuESTEnv env)
Definition: QuEST_cpu.c:1278
densmatr_mixDepolarisingDistributed
void densmatr_mixDepolarisingDistributed(Qureg qureg, const int targetQubit, qreal depolLevel)
Definition: QuEST_cpu.c:224
densmatr_mixDensityMatrix
void densmatr_mixDensityMatrix(Qureg combineQureg, qreal otherProb, Qureg otherQureg)
Definition: QuEST_cpu.c:889
copyStateToGPU
void copyStateToGPU(Qureg qureg)
In GPU mode, this copies the state-vector (or density matrix) from RAM (qureg.stateVec) to VRAM / GPU...
Definition: QuEST_cpu.c:36
statevec_controlledUnitaryDistributed
void statevec_controlledUnitaryDistributed(Qureg qureg, const 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 ...
Definition: QuEST_cpu.c:2347
alternateNormZeroingSomeAmpBlocks
void alternateNormZeroingSomeAmpBlocks(Qureg qureg, qreal norm, int normFirst, long long int startAmpInd, long long int numAmps, long long int blockSize)
Definition: QuEST_cpu.c:753
QuESTEnv::numRanks
int numRanks
Definition: QuEST.h:202
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 ...
Definition: QuEST_cpu.c:1969
ComplexMatrix4::imag
qreal imag[4][4]
Definition: QuEST.h:128
statevec_compactUnitaryLocal
void statevec_compactUnitaryLocal(Qureg qureg, const int targetQubit, Complex alpha, Complex beta)
Definition: QuEST_cpu.c:1656
densmatr_findProbabilityOfZeroLocal
qreal densmatr_findProbabilityOfZeroLocal(Qureg qureg, const int measureQubit)
Definition: QuEST_cpu.c:3111
normaliseSomeAmps
void normaliseSomeAmps(Qureg qureg, qreal norm, long long int startInd, long long int numAmps)
Definition: QuEST_cpu.c:743
statevec_getEnvironmentString
void statevec_getEnvironmentString(QuESTEnv env, Qureg qureg, char str[200])
Definition: QuEST_cpu.c:1358
ComplexMatrixN::real
qreal ** real
Definition: QuEST.h:139
densmatr_collapseToKnownProbOutcome
void densmatr_collapseToKnownProbOutcome(Qureg qureg, const int measureQubit, int outcome, qreal totalStateProb)
Renorms (/prob) every | * outcome * >< * outcome * | state, setting all others to zero.
Definition: QuEST_cpu.c:784
densmatr_initPureStateLocal
void densmatr_initPureStateLocal(Qureg targetQureg, Qureg copyQureg)
Definition: QuEST_cpu.c:1183
flipBit
__forceinline__ __device__ long long int flipBit(long long int number, int bitInd)
Definition: QuEST_gpu.cu:95
statevec_collapseToKnownProbOutcomeDistributedRenorm
void statevec_collapseToKnownProbOutcomeDistributedRenorm(Qureg qureg, const int measureQubit, const qreal totalProbability)
Renormalise parts of the state vector where measureQubit=0 or 1, based on the total probability of th...
Definition: QuEST_cpu.c:3422
densmatr_mixDepolarisingLocal
void densmatr_mixDepolarisingLocal(Qureg qureg, const int targetQubit, qreal depolLevel)
Definition: QuEST_cpu.c:125
Qureg
Represents a system of qubits.
Definition: QuEST.h:160
densmatr_calcInnerProductLocal
qreal densmatr_calcInnerProductLocal(Qureg a, Qureg b)
computes Tr(conjTrans(a) b) = sum of (a_ij^* b_ij)
Definition: QuEST_cpu.c:957
statevec_pauliXLocal
void statevec_pauliXLocal(Qureg qureg, const int targetQubit)
Definition: QuEST_cpu.c:2464
ComplexMatrixN::imag
qreal ** imag
Definition: QuEST.h:140
statevec_controlledUnitaryLocal
void statevec_controlledUnitaryLocal(Qureg qureg, const int controlQubit, const int targetQubit, ComplexMatrix2 u)
Definition: QuEST_cpu.c:2207
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 b...
Definition: QuEST_cpu.c:2024
statevec_compareStates
int statevec_compareStates(Qureg mq1, Qureg mq2, qreal precision)
Definition: QuEST_cpu.c:1643
QuEST_precision.h
statevec_multiControlledPhaseFlip
void statevec_multiControlledPhaseFlip(Qureg qureg, int *controlQubits, int numControlQubits)
Definition: QuEST_cpu.c:3291
Qureg::stateVec
ComplexArray stateVec
Computational state amplitudes - a subset thereof in the MPI version.
Definition: QuEST.h:179
statevec_unitaryLocal
void statevec_unitaryLocal(Qureg qureg, const int targetQubit, ComplexMatrix2 u)
Definition: QuEST_cpu.c:1900
ComplexMatrix2::real
qreal real[2][2]
Definition: QuEST.h:116
statevec_multiControlledUnitaryDistributed
void statevec_multiControlledUnitaryDistributed(Qureg qureg, const 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,...
Definition: QuEST_cpu.c:2413
Qureg::isDensityMatrix
int isDensityMatrix
Whether this instance is a density-state representation.
Definition: QuEST.h:163
densmatr_mixDephasing
void densmatr_mixDephasing(Qureg qureg, const int targetQubit, qreal dephase)
Definition: QuEST_cpu.c:79
statevec_controlledPauliYDistributed
void statevec_controlledPauliYDistributed(Qureg qureg, const int controlQubit, ComplexArray stateVecIn, ComplexArray stateVecOut, const int conjFac)
Definition: QuEST_cpu.c:2793
ComplexMatrixN::numQubits
int numQubits
Definition: QuEST.h:138
statevec_pauliYLocal
void statevec_pauliYLocal(Qureg qureg, const int targetQubit, const int conjFac)
Definition: QuEST_cpu.c:2647
statevec_multiControlledMultiQubitUnitaryLocal
void statevec_multiControlledMultiQubitUnitaryLocal(Qureg qureg, long long int ctrlMask, int *targs, const int numTargs, ComplexMatrixN u)
Definition: QuEST_cpu.c:1814
statevec_destroyQureg
void statevec_destroyQureg(Qureg qureg, QuESTEnv env)
Definition: QuEST_cpu.c:1316
densmatr_mixTwoQubitDepolarisingLocalPart1
void densmatr_mixTwoQubitDepolarisingLocalPart1(Qureg qureg, int qubit1, int qubit2, qreal delta)
Definition: QuEST_cpu.c:487
statevec_multiRotateZ
void statevec_multiRotateZ(Qureg qureg, long long int mask, qreal angle)
Definition: QuEST_cpu.c:3069
Qureg::numQubitsRepresented
int numQubitsRepresented
The number of qubits represented in either the state-vector or density matrix.
Definition: QuEST.h:165
Qureg::numAmpsTotal
long long int numAmpsTotal
Total number of amplitudes, which are possibly distributed among machines.
Definition: QuEST.h:172
Complex::real
qreal real
Definition: QuEST.h:105
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...
Definition: QuEST_cpu.c:3539
isOddParity
static int isOddParity(long long int number, int qb1, int qb2)
Definition: QuEST_cpu_internal.h:38
statevec_initStateFromSingleFile
int statevec_initStateFromSingleFile(Qureg *qureg, char filename[200], QuESTEnv env)
Definition: QuEST_cpu.c:1593
Complex::imag
qreal imag
Definition: QuEST.h:106
insertZeroBit
__forceinline__ __device__ long long int insertZeroBit(long long int number, int index)
Definition: QuEST_gpu.cu:99
statevec_pauliYDistributed
void statevec_pauliYDistributed(Qureg qureg, ComplexArray stateVecIn, ComplexArray stateVecOut, int updateUpper, const int conjFac)
Rotate a single qubit by +-{{0,-i},{i,0}.
Definition: QuEST_cpu.c:2704
statevec_controlledCompactUnitaryDistributed
void statevec_controlledCompactUnitaryDistributed(Qureg qureg, const 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 ...
Definition: QuEST_cpu.c:2285
Complex
Represents one complex number.
Definition: QuEST.h:103
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.
Definition: QuEST_cpu.c:2894
statevec_controlledPauliYLocal
void statevec_controlledPauliYLocal(Qureg qureg, const int controlQubit, const int targetQubit, const int conjFac)
Definition: QuEST_cpu.c:2740
statevec_reportStateToScreen
void statevec_reportStateToScreen(Qureg qureg, QuESTEnv env, int reportRank)
Definition: QuEST_cpu.c:1334
statevec_setWeightedQureg
void statevec_setWeightedQureg(Complex fac1, Qureg qureg1, Complex fac2, Qureg qureg2, Complex facOut, Qureg out)
Definition: QuEST_cpu.c:3579
densmatr_initPlusState
void densmatr_initPlusState(Qureg qureg)
Definition: QuEST_cpu.c:1153
statevec_initStateOfSingleQubit
void statevec_initStateOfSingleQubit(Qureg *qureg, int qubitId, int outcome)
Initialise the state vector of probability amplitudes such that one qubit is set to 'outcome' and all...
Definition: QuEST_cpu.c:1513
densmatr_mixDampingDistributed
void densmatr_mixDampingDistributed(Qureg qureg, const int targetQubit, qreal damping)
Definition: QuEST_cpu.c:299
statevec_calcInnerProductLocal
Complex statevec_calcInnerProductLocal(Qureg bra, Qureg ket)
Definition: QuEST_cpu.c:1070
ComplexMatrix2
Represents a 2x2 matrix of complex numbers.
Definition: QuEST.h:114
zeroSomeAmps
void zeroSomeAmps(Qureg qureg, long long int startInd, long long int numAmps)
Definition: QuEST_cpu.c:733
densmatr_oneQubitDegradeOffDiagonal
void densmatr_oneQubitDegradeOffDiagonal(Qureg qureg, const int targetQubit, qreal retain)
Definition: QuEST_cpu.c:48
densmatr_calcFidelityLocal
qreal densmatr_calcFidelityLocal(Qureg qureg, Qureg pureState)
computes a few dens-columns-worth of (vec^*T) dens * vec
Definition: QuEST_cpu.c:989