Merge pull request #212 from QuEST-Kit/develop

Bug fixes

QuEST/CMakeLists.txt

@@ -109,6 +109,10 @@ endif()
 
 if (GPUACCELERATED)
     find_package(CUDA REQUIRED)
+    # Stop nvcc sending c compile flags through using -Xcompiler and breaking
+    # on compilation of a cpp file receiving -std=c99. In long term should figure 
+    # out why CMAKE_C_FLAGS and not CMAKE_CXX_FLAGS are being sent through to a cpp file
+    set(CUDA_PROPAGATE_HOST_FLAGS FALSE)
 endif()
 
 

QuEST/include/QuEST.h

@@ -269,12 +269,12 @@ void destroyQureg(Qureg qureg, QuESTEnv env);
  * return from functions.
  * 
  * One can instead use getStaticComplexMatrixN() to create a ComplexMatrixN struct 
- * in the stack (which doesn't need to be explicitly freed).
+ * in the stack (which doesn't need to be later destroyed).
  *
  * @ingroup type
  * @param[in] numQubits the number of qubits of which the returned ComplexMatrixN will correspond
  * @returns a dynamic ComplexMatrixN struct, that is one where the .real and .imag
- *  fields are arrays kept in the heap and must be freed.
+ *  fields are arrays kept in the heap and must be later destroyed.
  * @author Tyson Jones
  */
 ComplexMatrixN createComplexMatrixN(int numQubits);
@@ -2017,9 +2017,9 @@ void mixDepolarising(Qureg qureg, const int targetQubit, qreal prob);
  *      if \p qureg is not a density matrix,
  *      or if \p targetQubit is outside [0, \p qureg.numQubitsRepresented),
  *      or if \p prob is not in [0, 1]
- * @author Nicolas Vogt of HQS
- * @author Ania Brown (patched)
- * @author Tyson Jones (doc)
+ * @author Nicolas Vogt of HQS (local CPU)
+ * @author Ania Brown (GPU, patched local CPU)
+ * @author Tyson Jones (distributed, doc)
  */
 void mixDamping(Qureg qureg, const int targetQubit, qreal prob);
 

QuEST/src/CPU/QuEST_cpu.c

@@ -81,11 +81,6 @@ void densmatr_mixDephasing(Qureg qureg, const int targetQubit, qreal dephase) {
     densmatr_oneQubitDegradeOffDiagonal(qureg, targetQubit, retain);
 }
 
-void densmatr_mixDampingDephase(Qureg qureg, const int targetQubit, qreal dephase) {
-    qreal retain=sqrt(1-dephase);
-    densmatr_oneQubitDegradeOffDiagonal(qureg, targetQubit, retain);
-}
-
 void densmatr_mixTwoQubitDephasing(Qureg qureg, const int qubit1, const int qubit2, qreal dephase) {
     qreal retain=1-dephase;
 
@@ -304,11 +299,11 @@ void densmatr_mixDepolarisingDistributed(Qureg qureg, const int targetQubit, qre
 void densmatr_mixDampingDistributed(Qureg qureg, const int targetQubit, qreal damping) {
     qreal retain=1-damping;
     qreal dephase=sqrt(1-damping);
-    // first do dephase part. 
-    // TODO -- this might be more efficient to do at the same time as the depolarise if we move to
-    // iterating over all elements in the state vector for the purpose of vectorisation
-    // TODO -- if we keep this split, move this function to densmatr_mixDepolarising()
-    densmatr_mixDampingDephase(qureg, targetQubit, dephase);
+
+    // multiply the off-diagonal (|0><1| and |1><0|) terms by sqrt(1-damping)
+    densmatr_oneQubitDegradeOffDiagonal(qureg, targetQubit, dephase);
+
+    // below, we modify the diagonals terms which require |1><1| to |0><0| communication
 
     long long int sizeInnerBlock, sizeInnerHalfBlock;
     long long int sizeOuterColumn, sizeOuterHalfColumn;
@@ -387,7 +382,6 @@ void densmatr_mixDampingDistributed(Qureg qureg, const int targetQubit, qreal da
     }    
 }
 
-// @TODO
 void densmatr_mixTwoQubitDepolarisingLocal(Qureg qureg, int qubit1, int qubit2, qreal delta, qreal gamma) {
     const long long int numTasks = qureg.numAmpsPerChunk;
     long long int innerMaskQubit1 = 1LL << qubit1;

QuEST/src/CPU/QuEST_cpu_distributed.c

@@ -1382,25 +1382,43 @@ void statevec_swapQubitAmps(Qureg qureg, int qb1, int qb2) {
 void statevec_multiControlledTwoQubitUnitary(Qureg qureg, long long int ctrlMask, const int q1, const int q2, ComplexMatrix4 u) {
     int q1FitsInNode = halfMatrixBlockFitsInChunk(qureg.numAmpsPerChunk, q1);
     int q2FitsInNode = halfMatrixBlockFitsInChunk(qureg.numAmpsPerChunk, q2);
-
-    if (q1FitsInNode && q2FitsInNode) {
+        
+    if (q1FitsInNode && q2FitsInNode) {    
         statevec_multiControlledTwoQubitUnitaryLocal(qureg, ctrlMask, q1, q2, u);
 
     } else if (q1FitsInNode) {
         int qSwap = (q1 > 0)? q1-1 : q1+1;
+
+        // ensure ctrl == qSwap, ensure ctrlMask updates under the swap
+        if (maskContainsBit(ctrlMask, qSwap))
+            ctrlMask = flipBit(flipBit(ctrlMask, q2), qSwap);
+
         statevec_swapQubitAmps(qureg, q2, qSwap);
         statevec_multiControlledTwoQubitUnitaryLocal(qureg, ctrlMask, q1, qSwap, u);
         statevec_swapQubitAmps(qureg, q2, qSwap);
 
     } else if (q2FitsInNode) {
         int qSwap = (q2 > 0)? q2-1 : q2+1;
+
+        // ensure ctrl == qSwap, ensure ctrlMask updates under the swap
+        if (maskContainsBit(ctrlMask, qSwap))
+            ctrlMask = flipBit(flipBit(ctrlMask, q1), qSwap);
+
         statevec_swapQubitAmps(qureg, q1, qSwap);
         statevec_multiControlledTwoQubitUnitaryLocal(qureg, ctrlMask, qSwap, q2, u);
         statevec_swapQubitAmps(qureg, q1, qSwap);
 
     } else {
+        // we know with certainty that both q1 and q2 >= 2
         int swap1 = 0;
         int swap2 = 1;
+
+        // if ctrl == swap1 or swap2, ensure ctrlMask updates under the swap
+        if (maskContainsBit(ctrlMask, swap1))
+            ctrlMask = flipBit(flipBit(ctrlMask, swap1), q1);
+        if (maskContainsBit(ctrlMask, swap2))
+            ctrlMask = flipBit(flipBit(ctrlMask, swap2), q2);
+
         statevec_swapQubitAmps(qureg, q1, swap1);
         statevec_swapQubitAmps(qureg, q2, swap2);
         statevec_multiControlledTwoQubitUnitaryLocal(qureg, ctrlMask, swap1, swap2, u);

QuEST/src/GPU/QuEST_gpu.cu

@@ -351,7 +351,7 @@ int GPUExists(void){
 }
 
 QuESTEnv createQuESTEnv(void) {
-    // init MPI environment
+
     if (!GPUExists()){
         printf("Trying to run GPU code with no GPU available\n");
         exit(EXIT_FAILURE);

QuEST/src/QuEST_common.c

@@ -446,9 +446,9 @@ void statevec_multiRotatePauli(
     }
 }
 
-void applyPauliProd(Qureg workspace, int* targetQubits, enum pauliOpType* pauliCodes, int numTargets) {
+/* produces both pauli|qureg> or pauli * qureg (as a density matrix) */
+void statevec_applyPauliProd(Qureg workspace, int* targetQubits, enum pauliOpType* pauliCodes, int numTargets) {
 
-    // produces both pauli|qureg> or pauli * qureg (as a density matrix)
     for (int i=0; i < numTargets; i++) {
         // (pauliCodes[i] == PAULI_I) applies no operation
         if (pauliCodes[i] == PAULI_X)
@@ -464,7 +464,7 @@ void applyPauliProd(Qureg workspace, int* targetQubits, enum pauliOpType* pauliC
 qreal statevec_calcExpecPauliProd(Qureg qureg, int* targetQubits, enum pauliOpType* pauliCodes, int numTargets, Qureg workspace) {
 
     statevec_cloneQureg(workspace, qureg);
-    applyPauliProd(workspace, targetQubits, pauliCodes, numTargets);
+    statevec_applyPauliProd(workspace, targetQubits, pauliCodes, numTargets);
 
     // compute the expected value
     qreal value;
@@ -505,11 +505,11 @@ void statevec_applyPauliSum(Qureg inQureg, enum pauliOpType* allCodes, qreal* te
         Complex zero = (Complex) {.real=0, .imag=0};
 
         // outQureg += coef paulis(inQureg)
-        applyPauliProd(inQureg, targs, &allCodes[t*numQb], numQb);
-        setWeightedQureg(coef, inQureg, iden, outQureg, zero, outQureg); 
+        statevec_applyPauliProd(inQureg, targs, &allCodes[t*numQb], numQb);
+        statevec_setWeightedQureg(coef, inQureg, iden, outQureg, zero, outQureg); 
 
         // undero paulis(inQureg), exploiting XX=YY=ZZ=I
-        applyPauliProd(inQureg, targs, &allCodes[t*numQb], numQb);
+        statevec_applyPauliProd(inQureg, targs, &allCodes[t*numQb], numQb);
     }
 }