added calcAmpSum backend

quest/src/api/calculations.cpp

@@ -134,9 +134,10 @@ qreal calcRealAmpSum(Qureg qureg) {
     validate_quregFields(qureg, __func__);
     validate_quregHasEvenNumQubits(qureg, __func__);
 
-    /// @todo
-    /// implement and call new backend function
-    return -1;
+    // logic is identical for both statevectors and density matrices
+    qcomp value = localiser_statevec_calcAmpSum(qureg);
+
+    return std::real(value);
 }
 
 

quest/src/core/accelerator.cpp

@@ -182,6 +182,20 @@ using std::min;
 
 
 
+/*
+ * PR DEMOS
+ */
+
+
+qcomp accel_statevec_calcAmpSum(Qureg qureg) {
+
+    return (qureg.isGpuAccelerated)?
+        gpu_statevec_calcAmpSum(qureg):
+        cpu_statevec_calcAmpSum(qureg);
+}
+
+
+
 /*
  * GETTERS 
  */

quest/src/core/accelerator.hpp

@@ -147,6 +147,13 @@ using std::vector;
         name = compileval;
 
 
+/*
+ * PR DEMOS
+ */
+
+qcomp accel_statevec_calcAmpSum(Qureg qureg);
+
+
 /*
  * GETTERS 
  */

quest/src/core/localiser.cpp

@@ -461,6 +461,27 @@ void exchangeAmpsToBuffersWhereQubitsAreInStates(Qureg qureg, int pairRank, vect
 
 
 
+/*
+ * PR DEMOS
+ */
+
+
+qcomp localiser_statevec_calcAmpSum(Qureg qureg) {
+
+    // each node computes the sum of its local amps
+    // in an embarrassingly parallel fashion
+    qcomp out =  accel_statevec_calcAmpSum(qureg);
+
+    // node partial sums are combined, and every
+    // node receives a copy of the full sum (consensus)
+    if (qureg.isDistributed)
+        comm_reduceAmp(&out);
+
+    return out;
+}
+
+
+
 /*
  * GETTERS
  */

quest/src/core/localiser.hpp

@@ -23,6 +23,13 @@
 using std::vector;
 
 
+/*
+ * PR DEMOS
+ */
+
+qcomp localiser_statevec_calcAmpSum(Qureg qureg);
+
+
 /*
  * GETTERS
  */

quest/src/cpu/cpu_subroutines.cpp

@@ -40,6 +40,34 @@ using std::vector;
 
 
 
+/*
+ * PR DEMOS
+ */
+
+
+qcomp cpu_statevec_calcAmpSum(Qureg qureg) {
+
+    // we brazenly perform this all-state reduction without
+    // a single thought to finite precision effects - arrest me! 
+
+    // separately reduce real and imag components to make MSVC happy
+    qreal outRe = 0;
+    qreal outIm = 0;
+
+    // every local amplitude contributes to the sum
+    qindex numIts = qureg.numAmpsPerNode;
+
+    #pragma omp parallel for reduction(+:outRe,outIm) if(qureg.isMultithreaded)
+    for (qindex n=0; n<numIts; n++) {
+        outRe += std::real(qureg.cpuAmps[n]);
+        outIm += std::imag(qureg.cpuAmps[n]);
+    }
+
+    return qcomp(outRe, outIm);
+}
+
+
+
 /*
  * GETTERS
  */

quest/src/cpu/cpu_subroutines.hpp

@@ -20,6 +20,13 @@
 using std::vector;
 
 
+/*
+ * PR DEMOS
+ */
+
+qcomp cpu_statevec_calcAmpSum(Qureg qureg);
+
+
 /*
  * GETTERS
  */

quest/src/gpu/gpu_subroutines.cpp

@@ -65,6 +65,32 @@ using std::vector;
 
 
 
+/*
+ * PR DEMOS
+ */
+
+
+qcomp gpu_statevec_calcAmpSum(Qureg qureg) {
+
+#if COMPILE_CUDA || COMPILE_CUQUANTUM
+
+    // there is no cuQuantum function for this facility,
+    // nor is there a need to implement a custom kernel;
+    // we will leverage an exting Thrust facility
+    cu_qcomp out =  thrust_statevec_calcAmpSum(qureg);
+
+    // the GPU backend uses a CUDA-equivalent of 'qcomp'
+    // which we must cast back and forth between (ew!)
+    return toQcomp(out);
+
+#else
+    error_gpuSimButGpuNotCompiled();
+    return -1;
+#endif
+}
+
+
+
 /*
  * GETTERS
  */

quest/src/gpu/gpu_subroutines.hpp

@@ -17,6 +17,13 @@
 using std::vector;
 
 
+/*
+ * PR DEMOS
+ */
+
+qcomp gpu_statevec_calcAmpSum(Qureg qureg);
+
+
 /*
  * GETTERS
  */

quest/src/gpu/gpu_thrust.cuh

@@ -635,6 +635,27 @@ struct functor_setRandomStateVecAmp : public thrust::unary_function<qindex,cu_qc
 
 
 
+/*
+ * PR DEMOS
+ */
+
+
+cu_qcomp thrust_statevec_calcAmpSum(Qureg qureg) {
+
+    // beware that we must explicitly instantiate a 
+    // qcomp type here, rather than pass reduce() a
+    // literal, which would cause a silent bug (grr)
+    cu_qcomp init = getCuQcomp(0, 0);
+
+    cu_qcomp out = thrust::reduce(
+        getStartPtr(qureg), getEndPtr(qureg), 
+        init, thrust::plus<cu_qcomp>());
+
+    return out;
+}
+
+
+
 /*
  * MATRIX INITIALISATION
  */