From 44909e20c875507cd1d772a8a4c9a15d320a85cd Mon Sep 17 00:00:00 2001
From: Tyson Jones <tyson.jones.input@gmail.com>
Date: Tue, 6 May 2025 00:02:30 +0200
Subject: [PATCH] added calcAmpSum backend

---
 quest/src/api/calculations.cpp    |  7 ++++---
 quest/src/core/accelerator.cpp    | 14 ++++++++++++++
 quest/src/core/accelerator.hpp    |  7 +++++++
 quest/src/core/localiser.cpp      | 21 +++++++++++++++++++++
 quest/src/core/localiser.hpp      |  7 +++++++
 quest/src/cpu/cpu_subroutines.cpp | 28 ++++++++++++++++++++++++++++
 quest/src/cpu/cpu_subroutines.hpp |  7 +++++++
 quest/src/gpu/gpu_subroutines.cpp | 26 ++++++++++++++++++++++++++
 quest/src/gpu/gpu_subroutines.hpp |  7 +++++++
 quest/src/gpu/gpu_thrust.cuh      | 21 +++++++++++++++++++++
 10 files changed, 142 insertions(+), 3 deletions(-)

diff --git a/quest/src/api/calculations.cpp b/quest/src/api/calculations.cpp
index 72c34a14..19927ca3 100644
--- a/quest/src/api/calculations.cpp
+++ b/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);
 }
 
 
diff --git a/quest/src/core/accelerator.cpp b/quest/src/core/accelerator.cpp
index 9f877fba..9426e326 100644
--- a/quest/src/core/accelerator.cpp
+++ b/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 
  */
diff --git a/quest/src/core/accelerator.hpp b/quest/src/core/accelerator.hpp
index 01cf3efc..8309c4fd 100644
--- a/quest/src/core/accelerator.hpp
+++ b/quest/src/core/accelerator.hpp
@@ -147,6 +147,13 @@ using std::vector;
         name = compileval;
 
 
+/*
+ * PR DEMOS
+ */
+
+qcomp accel_statevec_calcAmpSum(Qureg qureg);
+
+
 /*
  * GETTERS 
  */
diff --git a/quest/src/core/localiser.cpp b/quest/src/core/localiser.cpp
index 7cc84ba0..65a8b5f9 100644
--- a/quest/src/core/localiser.cpp
+++ b/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
  */
diff --git a/quest/src/core/localiser.hpp b/quest/src/core/localiser.hpp
index a9b0391a..075c833a 100644
--- a/quest/src/core/localiser.hpp
+++ b/quest/src/core/localiser.hpp
@@ -23,6 +23,13 @@
 using std::vector;
 
 
+/*
+ * PR DEMOS
+ */
+
+qcomp localiser_statevec_calcAmpSum(Qureg qureg);
+
+
 /*
  * GETTERS
  */
diff --git a/quest/src/cpu/cpu_subroutines.cpp b/quest/src/cpu/cpu_subroutines.cpp
index 9c90e08c..132f8989 100644
--- a/quest/src/cpu/cpu_subroutines.cpp
+++ b/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
  */
diff --git a/quest/src/cpu/cpu_subroutines.hpp b/quest/src/cpu/cpu_subroutines.hpp
index 661de150..78d00355 100644
--- a/quest/src/cpu/cpu_subroutines.hpp
+++ b/quest/src/cpu/cpu_subroutines.hpp
@@ -20,6 +20,13 @@
 using std::vector;
 
 
+/*
+ * PR DEMOS
+ */
+
+qcomp cpu_statevec_calcAmpSum(Qureg qureg);
+
+
 /*
  * GETTERS
  */
diff --git a/quest/src/gpu/gpu_subroutines.cpp b/quest/src/gpu/gpu_subroutines.cpp
index 0b0bc46a..1f7315c2 100644
--- a/quest/src/gpu/gpu_subroutines.cpp
+++ b/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
  */
diff --git a/quest/src/gpu/gpu_subroutines.hpp b/quest/src/gpu/gpu_subroutines.hpp
index 04fabe0a..2085ab0f 100644
--- a/quest/src/gpu/gpu_subroutines.hpp
+++ b/quest/src/gpu/gpu_subroutines.hpp
@@ -17,6 +17,13 @@
 using std::vector;
 
 
+/*
+ * PR DEMOS
+ */
+
+qcomp gpu_statevec_calcAmpSum(Qureg qureg);
+
+
 /*
  * GETTERS
  */
diff --git a/quest/src/gpu/gpu_thrust.cuh b/quest/src/gpu/gpu_thrust.cuh
index 7cb45e93..a51e88cf 100644
--- a/quest/src/gpu/gpu_thrust.cuh
+++ b/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
  */