Improve applyCompMatr summation accuracy

quest/src/cpu/cpu_subroutines.cpp

@@ -607,7 +607,8 @@ void cpu_statevec_anyCtrlAnyTargDenseMatr_sub(Qureg qureg, ConstList64 ctrls, Co
 
                 // i = nth local index where ctrls are active and targs form value k
                 qindex i = setBits(i0, targs.data(), numTargBits, k); // loop may be unrolled
-                amps[i] = getCpuQcomp(0, 0);
+                cpu_qcomp sum = getCpuQcomp(0, 0);
+                cpu_qcomp correction = getCpuQcomp(0, 0);
 
                 // loop may be unrolled
                 for (qindex j=0; j<numTargAmps; j++) {
@@ -624,18 +625,19 @@ void cpu_statevec_anyCtrlAnyTargDenseMatr_sub(Qureg qureg, ConstList64 ctrls, Co
                     if constexpr (ApplyConj)
                         elem = conj(elem);
 
-                    amps[i] += elem * cache[j];
-
-                    /// @todo
-                    /// qureg.cpuAmps[i] is being serially updated by only this thread,
-                    /// so is a candidate for Kahan summation for improved numerical
-                    /// stability. Explore whether this is time-free and worthwhile!
-                    ///
-                    /// BEWARE that Kahan summation may be incompatible with
-                    /// the commutator tricks used in base_qcomp's (ancestor
-                    /// of cpu_qcomp) arithmetic operator overloads. Check
-                    /// base_qcomp.hpp before implementing compensation.
+                    cpu_qcomp term = elem * cache[j];
+                    qreal re = term.re - correction.re;
+                    qreal nextRe = sum.re + re;
+                    correction.re = (nextRe - sum.re) - re;
+                    sum.re = nextRe;
+
+                    qreal im = term.im - correction.im;
+                    qreal nextIm = sum.im + im;
+                    correction.im = (nextIm - sum.im) - im;
+                    sum.im = nextIm;
                 }
+
+                amps[i] = sum;
             }
         }
     }

tests/unit/operations.cpp

@@ -40,6 +40,8 @@
 #include "tests/utils/macros.hpp"
 #include "tests/utils/random.hpp"
 
+#include <cmath>
+#include <limits>
 #include <tuple>
 
 using std::tuple;
@@ -1315,6 +1317,42 @@ TEST_ALL_CTRL_OPERATIONS( PauliGadget, any, pauligad, nullptr );
 TEST_ALL_CTRL_OPERATIONS( CompMatr1, one, compmatr, nullptr );
 TEST_ALL_CTRL_OPERATIONS( CompMatr2, two, compmatr, nullptr );
 TEST_ALL_CTRL_OPERATIONS( CompMatr,  any, compmatr, nullptr );
+
+TEST_CASE( "applyCompMatr uses compensated summation", TEST_CATEGORY_OPS ) {
+
+    int numTargs = 4;
+    qindex numAmps = getPow2(numTargs);
+    int targets[] = {0,1,2,3};
+
+    Qureg qureg = createQureg(numTargs);
+    CompMatr matrix = createCompMatr(numTargs);
+
+    qvector amps(numAmps, qcomp(1, 0));
+    setQuregAmps(qureg, 0, amps.data(), amps.size());
+
+    for (qindex i=0; i<matrix.numRows; i++)
+        for (qindex j=0; j<matrix.numRows; j++)
+            matrix.cpuElems[i][j] = 0;
+
+    qreal large = std::ldexp(qreal(1), std::numeric_limits<qreal>::digits);
+    matrix.cpuElems[0][0] = qcomp( large,  large);
+    for (qindex j=1; j<matrix.numRows-1; j++)
+        matrix.cpuElems[0][j] = qcomp(1, -1);
+    matrix.cpuElems[0][matrix.numRows-1] = qcomp(-large, -large);
+    syncCompMatr(matrix);
+
+    setQuESTValidationEpsilon(0);
+    applyCompMatr(qureg, targets, numTargs, matrix);
+    setQuESTValidationEpsilonToDefault();
+
+    qcomp amp = getQuregAmp(qureg, 0);
+    REQUIRE( std::real(amp) == qreal(numAmps - 2) );
+    REQUIRE( std::imag(amp) == -qreal(numAmps - 2) );
+
+    destroyCompMatr(matrix);
+    destroyQureg(qureg);
+}
+
 TEST_ALL_CTRL_OPERATIONS( DiagMatr1, one, diagmatr, nullptr );
 TEST_ALL_CTRL_OPERATIONS( DiagMatr2, two, diagmatr, nullptr );
 TEST_ALL_CTRL_OPERATIONS( DiagMatr,  any, diagmatr, nullptr );