Support measurement sampling seed for cutensornet backends

runtime/nvqir/cutensornet/CMakeLists.txt

@@ -60,8 +60,8 @@ set(CUTENSORNET_PATCH ${CMAKE_MATCH_1})
 
 set(CUTENSORNET_VERSION ${CUTENSORNET_MAJOR}.${CUTENSORNET_MINOR}.${CUTENSORNET_PATCH})
 message(STATUS "Found cutensornet version: ${CUTENSORNET_VERSION}")
-# We need cutensornet v2.3+
-if (${CUTENSORNET_VERSION} VERSION_GREATER_EQUAL "2.3")
+# We need cutensornet v2.5.0+
+if (${CUTENSORNET_VERSION} VERSION_GREATER_EQUAL "2.5")
   set (BASE_TENSOR_BACKEND_SRS  
         simulator_cutensornet.cpp 
         tensornet_spin_op.cpp

runtime/nvqir/cutensornet/mps_simulation_state.cpp

@@ -26,9 +26,11 @@ std::size_t MPSSimulationState::getNumQubits() const {
 MPSSimulationState::MPSSimulationState(std::unique_ptr<TensorNetState> inState,
                                        const std::vector<MPSTensor> &mpsTensors,
                                        ScratchDeviceMem &inScratchPad,
-                                       cutensornetHandle_t cutnHandle)
+                                       cutensornetHandle_t cutnHandle,
+                                       std::mt19937 &randomEngine)
     : m_cutnHandle(cutnHandle), state(std::move(inState)),
-      m_mpsTensors(mpsTensors), scratchPad(inScratchPad) {}
+      m_mpsTensors(mpsTensors), scratchPad(inScratchPad),
+      m_randomEngine(randomEngine) {}
 
 MPSSimulationState::~MPSSimulationState() { deallocate(); }
 
@@ -266,8 +268,8 @@ MPSSimulationState::getAmplitude(const std::vector<int> &basisState) {
   }
 
   if (getNumQubits() > 1) {
-    TensorNetState basisTensorNetState(basisState, scratchPad,
-                                       state->getInternalContext());
+    TensorNetState basisTensorNetState(
+        basisState, scratchPad, state->getInternalContext(), m_randomEngine);
     // Note: this is a basis state, hence bond dim == 1
     std::vector<MPSTensor> basisStateTensors = basisTensorNetState.factorizeMPS(
         1, std::numeric_limits<double>::min(),
@@ -369,7 +371,8 @@ static Eigen::MatrixXcd reshapeStateVec(const Eigen::VectorXcd &stateVec) {
 
 MPSSimulationState::MpsStateData MPSSimulationState::createFromStateVec(
     cutensornetHandle_t cutnHandle, ScratchDeviceMem &inScratchPad,
-    std::size_t size, std::complex<double> *ptr, int bondDim) {
+    std::size_t size, std::complex<double> *ptr, int bondDim,
+    std::mt19937 &randomEngine) {
   const std::size_t numQubits = std::log2(size);
   // Reverse the qubit order to match cutensornet convention
   auto newStateVec = TensorNetState::reverseQubitOrder(
@@ -386,8 +389,8 @@ MPSSimulationState::MpsStateData MPSSimulationState::createFromStateVec(
     MPSTensor stateTensor;
     stateTensor.deviceData = d_tensor;
     stateTensor.extents = std::vector<int64_t>{2};
-    auto state = TensorNetState::createFromMpsTensors({stateTensor},
-                                                      inScratchPad, cutnHandle);
+    auto state = TensorNetState::createFromMpsTensors(
+        {stateTensor}, inScratchPad, cutnHandle, randomEngine);
     return {std::move(state), std::vector<MPSTensor>{stateTensor}};
   }
 
@@ -460,7 +463,7 @@ MPSSimulationState::MpsStateData MPSSimulationState::createFromStateVec(
   mpsTensors.emplace_back(stateTensor);
   assert(mpsTensors.size() == numQubits);
   auto state = TensorNetState::createFromMpsTensors(mpsTensors, inScratchPad,
-                                                    cutnHandle);
+                                                    cutnHandle, randomEngine);
   return {std::move(state), mpsTensors};
 }
 
@@ -486,16 +489,17 @@ MPSSimulationState::createFromSizeAndPtr(std::size_t size, void *ptr,
       MPSTensor stateTensor{d_tensor, mpsExtents};
       mpsTensors.emplace_back(stateTensor);
     }
-    auto state = TensorNetState::createFromMpsTensors(mpsTensors, scratchPad,
-                                                      m_cutnHandle);
-    return std::make_unique<MPSSimulationState>(std::move(state), mpsTensors,
-                                                scratchPad, m_cutnHandle);
+    auto state = TensorNetState::createFromMpsTensors(
+        mpsTensors, scratchPad, m_cutnHandle, m_randomEngine);
+    return std::make_unique<MPSSimulationState>(
+        std::move(state), mpsTensors, scratchPad, m_cutnHandle, m_randomEngine);
   }
-  auto [state, mpsTensors] = createFromStateVec(
-      m_cutnHandle, scratchPad, size,
-      reinterpret_cast<std::complex<double> *>(ptr), MPSSettings().maxBond);
-  return std::make_unique<MPSSimulationState>(std::move(state), mpsTensors,
-                                              scratchPad, m_cutnHandle);
+  auto [state, mpsTensors] =
+      createFromStateVec(m_cutnHandle, scratchPad, size,
+                         reinterpret_cast<std::complex<double> *>(ptr),
+                         MPSSettings().maxBond, m_randomEngine);
+  return std::make_unique<MPSSimulationState>(
+      std::move(state), mpsTensors, scratchPad, m_cutnHandle, m_randomEngine);
 }
 
 MPSSettings::MPSSettings() {

runtime/nvqir/cutensornet/mps_simulation_state.h

@@ -35,7 +35,8 @@ class MPSSimulationState : public cudaq::SimulationState {
   MPSSimulationState(std::unique_ptr<TensorNetState> inState,
                      const std::vector<MPSTensor> &mpsTensors,
                      ScratchDeviceMem &inScratchPad,
-                     cutensornetHandle_t cutnHandle);
+                     cutensornetHandle_t cutnHandle,
+                     std::mt19937 &randomEngine);
 
   MPSSimulationState(const MPSSimulationState &) = delete;
   MPSSimulationState &operator=(const MPSSimulationState &) = delete;
@@ -84,7 +85,7 @@ class MPSSimulationState : public cudaq::SimulationState {
                                          ScratchDeviceMem &inScratchPad,
                                          std::size_t size,
                                          std::complex<double> *data,
-                                         int bondDim);
+                                         int bondDim, std::mt19937 &randomEngine);
 
   /// Retrieve the MPS tensors
   std::vector<MPSTensor> getMpsTensors() const { return m_mpsTensors; }
@@ -105,6 +106,7 @@ class MPSSimulationState : public cudaq::SimulationState {
   // This speeds up sequential state amplitude accessors for small states.
   static constexpr std::size_t g_maxQubitsForStateContraction = 30;
   std::vector<std::complex<double>> m_contractedStateVec;
+  std::mt19937 &m_randomEngine;
 };
 
 } // namespace nvqir

runtime/nvqir/cutensornet/simulator_cutensornet.cpp

@@ -289,8 +289,8 @@ void SimulatorTensorNetBase::setToZeroState() {
   const auto numQubits = m_state->getNumQubits();
   m_state.reset();
   // Re-create a zero state of the same size
-  m_state =
-      std::make_unique<TensorNetState>(numQubits, scratchPad, m_cutnHandle);
+  m_state = std::make_unique<TensorNetState>(numQubits, scratchPad,
+                                             m_cutnHandle, m_randomEngine);
 }
 
 void SimulatorTensorNetBase::swap(const std::vector<std::size_t> &ctrlBits,

runtime/nvqir/cutensornet/simulator_cutensornet.h

@@ -93,10 +93,8 @@ class SimulatorTensorNetBase : public nvqir::CircuitSimulatorBase<double> {
   std::unique_ptr<TensorNetState> m_state;
   std::unordered_map<std::string, void *> m_gateDeviceMemCache;
   ScratchDeviceMem scratchPad;
-  // Note: cutensornet sample API uses an internal random engine that doesn't
-  // support random seed. This engine only affects the mid-circuit measurements
-  // whereby this simulator generates a random probability value.
-  // See also: https://github.com/NVIDIA/cuda-quantum/issues/895
+  // Random number generator for generating 32-bit numbers with a state size of
+  // 19937 bits for measurements.
   std::mt19937 m_randomEngine;
 };
 

runtime/nvqir/cutensornet/simulator_mps_register.cpp

@@ -47,8 +47,8 @@ class SimulatorMPS : public SimulatorTensorNetBase {
       throw std::invalid_argument(
           "[SimulatorMPS simulator] Incompatible state input");
     if (!m_state) {
-      m_state = TensorNetState::createFromMpsTensors(casted->getMpsTensors(),
-                                                     scratchPad, m_cutnHandle);
+      m_state = TensorNetState::createFromMpsTensors(
+          casted->getMpsTensors(), scratchPad, m_cutnHandle, m_randomEngine);
     } else {
       // Expand an existing state: Append MPS tensors
       // Factor the existing state
@@ -76,8 +76,8 @@ class SimulatorMPS : public SimulatorTensorNetBase {
                                      tensorSizeBytes, cudaMemcpyDefault));
         tensors.emplace_back(MPSTensor(mpsTensor, extents));
       }
-      m_state = TensorNetState::createFromMpsTensors(tensors, scratchPad,
-                                                     m_cutnHandle);
+      m_state = TensorNetState::createFromMpsTensors(
+          tensors, scratchPad, m_cutnHandle, m_randomEngine);
     }
   }
 
@@ -164,13 +164,13 @@ class SimulatorMPS : public SimulatorTensorNetBase {
     LOG_API_TIME();
     if (!m_state) {
       if (!ptr) {
-        m_state = std::make_unique<TensorNetState>(numQubits, scratchPad,
-                                                   m_cutnHandle);
+        m_state = std::make_unique<TensorNetState>(
+            numQubits, scratchPad, m_cutnHandle, m_randomEngine);
       } else {
         auto [state, mpsTensors] = MPSSimulationState::createFromStateVec(
             m_cutnHandle, scratchPad, 1ULL << numQubits,
             reinterpret_cast<std::complex<double> *>(const_cast<void *>(ptr)),
-            m_settings.maxBond);
+            m_settings.maxBond, m_randomEngine);
         m_state = std::move(state);
       }
     } else {
@@ -194,14 +194,14 @@ class SimulatorMPS : public SimulatorTensorNetBase {
                                        cudaMemcpyHostToDevice));
           tensors.emplace_back(MPSTensor(mpsTensor, extents));
         }
-        m_state = TensorNetState::createFromMpsTensors(tensors, scratchPad,
-                                                       m_cutnHandle);
+        m_state = TensorNetState::createFromMpsTensors(
+            tensors, scratchPad, m_cutnHandle, m_randomEngine);
       } else {
         // Non-zero state needs to be factorized and appended.
         auto [state, mpsTensors] = MPSSimulationState::createFromStateVec(
             m_cutnHandle, scratchPad, 1ULL << numQubits,
             reinterpret_cast<std::complex<double> *>(const_cast<void *>(ptr)),
-            m_settings.maxBond);
+            m_settings.maxBond, m_randomEngine);
         auto tensors =
             m_state->factorizeMPS(m_settings.maxBond, m_settings.absCutoff,
                                   m_settings.relCutoff, m_settings.svdAlgo);
@@ -214,8 +214,8 @@ class SimulatorMPS : public SimulatorTensorNetBase {
         mpsTensors.front().extents = extents;
         // Combine the list
         tensors.insert(tensors.end(), mpsTensors.begin(), mpsTensors.end());
-        m_state = TensorNetState::createFromMpsTensors(tensors, scratchPad,
-                                                       m_cutnHandle);
+        m_state = TensorNetState::createFromMpsTensors(
+            tensors, scratchPad, m_cutnHandle, m_randomEngine);
       }
     }
   }
@@ -224,16 +224,17 @@ class SimulatorMPS : public SimulatorTensorNetBase {
     LOG_API_TIME();
 
     if (!m_state || m_state->getNumQubits() == 0)
-      return std::make_unique<MPSSimulationState>(std::move(m_state),
-                                                  std::vector<MPSTensor>{},
-                                                  scratchPad, m_cutnHandle);
+      return std::make_unique<MPSSimulationState>(
+          std::move(m_state), std::vector<MPSTensor>{}, scratchPad,
+          m_cutnHandle, m_randomEngine);
 
     if (m_state->getNumQubits() > 1) {
       std::vector<MPSTensor> tensors =
           m_state->factorizeMPS(m_settings.maxBond, m_settings.absCutoff,
                                 m_settings.relCutoff, m_settings.svdAlgo);
       return std::make_unique<MPSSimulationState>(std::move(m_state), tensors,
-                                                  scratchPad, m_cutnHandle);
+                                                  scratchPad, m_cutnHandle,
+                                                  m_randomEngine);
     }
 
     auto [d_tensor, numElements] = m_state->contractStateVectorInternal({});
@@ -244,7 +245,7 @@ class SimulatorMPS : public SimulatorTensorNetBase {
 
     return std::make_unique<MPSSimulationState>(
         std::move(m_state), std::vector<MPSTensor>{stateTensor}, scratchPad,
-        m_cutnHandle);
+        m_cutnHandle, m_randomEngine);
   }
 
   virtual ~SimulatorMPS() noexcept {

runtime/nvqir/cutensornet/simulator_tensornet_register.cpp

@@ -47,22 +47,22 @@ class SimulatorTensorNet : public SimulatorTensorNetBase {
 
   std::unique_ptr<cudaq::SimulationState> getSimulationState() override {
     LOG_API_TIME();
-    return std::make_unique<TensorNetSimulationState>(std::move(m_state),
-                                                      scratchPad, m_cutnHandle);
+    return std::make_unique<TensorNetSimulationState>(
+        std::move(m_state), scratchPad, m_cutnHandle, m_randomEngine);
   }
 
   void addQubitsToState(std::size_t numQubits, const void *ptr) override {
     LOG_API_TIME();
     if (!m_state) {
       if (!ptr) {
-        m_state = std::make_unique<TensorNetState>(numQubits, scratchPad,
-                                                   m_cutnHandle);
+        m_state = std::make_unique<TensorNetState>(
+            numQubits, scratchPad, m_cutnHandle, m_randomEngine);
       } else {
         auto *casted =
             reinterpret_cast<std::complex<double> *>(const_cast<void *>(ptr));
         std::span<std::complex<double>> stateVec(casted, 1ULL << numQubits);
-        m_state = TensorNetState::createFromStateVector(stateVec, scratchPad,
-                                                        m_cutnHandle);
+        m_state = TensorNetState::createFromStateVector(
+            stateVec, scratchPad, m_cutnHandle, m_randomEngine);
       }
     } else {
       if (!ptr) {
@@ -85,9 +85,9 @@ class SimulatorTensorNet : public SimulatorTensorNetBase {
       throw std::invalid_argument(
           "[Tensornet simulator] Incompatible state input");
     if (!m_state) {
-      m_state = TensorNetState::createFromOpTensors(in_state.getNumQubits(),
-                                                    casted->getAppliedTensors(),
-                                                    scratchPad, m_cutnHandle);
+      m_state = TensorNetState::createFromOpTensors(
+          in_state.getNumQubits(), casted->getAppliedTensors(), scratchPad,
+          m_cutnHandle, m_randomEngine);
     } else {
       // Expand an existing state:
       //  (1) Create a blank tensor network with combined number of qubits

runtime/nvqir/cutensornet/tensornet_state.cpp

@@ -15,8 +15,9 @@ namespace nvqir {
 
 TensorNetState::TensorNetState(std::size_t numQubits,
                                ScratchDeviceMem &inScratchPad,
-                               cutensornetHandle_t handle)
-    : m_numQubits(numQubits), m_cutnHandle(handle), scratchPad(inScratchPad) {
+                               cutensornetHandle_t handle, std::mt19937 &randomEngine)
+    : m_numQubits(numQubits), m_cutnHandle(handle), scratchPad(inScratchPad),
+      m_randomEngine(randomEngine) {
   const std::vector<int64_t> qubitDims(m_numQubits, 2);
   HANDLE_CUTN_ERROR(cutensornetCreateState(
       m_cutnHandle, CUTENSORNET_STATE_PURITY_PURE, m_numQubits,
@@ -25,8 +26,8 @@ TensorNetState::TensorNetState(std::size_t numQubits,
 
 TensorNetState::TensorNetState(const std::vector<int> &basisState,
                                ScratchDeviceMem &inScratchPad,
-                               cutensornetHandle_t handle)
-    : TensorNetState(basisState.size(), inScratchPad, handle) {
+                               cutensornetHandle_t handle, std::mt19937 &randomEngine)
+    : TensorNetState(basisState.size(), inScratchPad, handle, randomEngine) {
   constexpr std::complex<double> h_xGate[4] = {0.0, 1.0, 1.0, 0.0};
   constexpr auto sizeBytes = 4 * sizeof(std::complex<double>);
   void *d_gate{nullptr};
@@ -43,8 +44,8 @@ TensorNetState::TensorNetState(const std::vector<int> &basisState,
 }
 
 std::unique_ptr<TensorNetState> TensorNetState::clone() const {
-  return createFromOpTensors(m_numQubits, m_tensorOps, scratchPad,
-                             m_cutnHandle);
+  return createFromOpTensors(m_numQubits, m_tensorOps, scratchPad, m_cutnHandle,
+                             m_randomEngine);
 }
 
 void TensorNetState::applyGate(const std::vector<int32_t> &controlQubits,
@@ -166,6 +167,16 @@ TensorNetState::sample(const std::vector<int32_t> &measuredBitIds,
     HANDLE_CUTN_ERROR(cutensornetSamplerConfigure(
         m_cutnHandle, sampler, CUTENSORNET_SAMPLER_OPT_NUM_HYPER_SAMPLES,
         &numHyperSamples, sizeof(numHyperSamples)));
+
+    // Generate a random seed from the backend simulator's random engine.
+    // Note: Even after a random seed setting at the user's level,
+    // consecutive `cudaq::sample` calls will still return different results
+    // (yet deterministic), i.e., the seed that we send to cutensornet should
+    // not be the user's seed.
+    const int32_t rndSeed = m_randomEngine();
+    HANDLE_CUTN_ERROR(cutensornetSamplerConfigure(
+        m_cutnHandle, sampler, CUTENSORNET_SAMPLER_CONFIG_DETERMINISTIC,
+        &rndSeed, sizeof(rndSeed)));
   }
 
   // Prepare the quantum circuit sampler
@@ -591,12 +602,12 @@ std::complex<double> TensorNetState::computeExpVal(
 
 std::unique_ptr<TensorNetState> TensorNetState::createFromMpsTensors(
     const std::vector<MPSTensor> &in_mpsTensors, ScratchDeviceMem &inScratchPad,
-    cutensornetHandle_t handle) {
+    cutensornetHandle_t handle, std::mt19937 &randomEngine) {
   LOG_API_TIME();
   if (in_mpsTensors.empty())
     throw std::invalid_argument("Empty MPS tensor list");
   auto state = std::make_unique<TensorNetState>(in_mpsTensors.size(),
-                                                inScratchPad, handle);
+                                                inScratchPad, handle, randomEngine);
   std::vector<const int64_t *> extents;
   std::vector<void *> tensorData;
   for (const auto &tensor : in_mpsTensors) {
@@ -613,10 +624,11 @@ std::unique_ptr<TensorNetState> TensorNetState::createFromMpsTensors(
 /// operators.
 std::unique_ptr<TensorNetState> TensorNetState::createFromOpTensors(
     std::size_t numQubits, const std::vector<AppliedTensorOp> &opTensors,
-    ScratchDeviceMem &inScratchPad, cutensornetHandle_t handle) {
+    ScratchDeviceMem &inScratchPad, cutensornetHandle_t handle,
+    std::mt19937 &randomEngine) {
   LOG_API_TIME();
   auto state =
-      std::make_unique<TensorNetState>(numQubits, inScratchPad, handle);
+      std::make_unique<TensorNetState>(numQubits, inScratchPad, handle, randomEngine);
   for (const auto &op : opTensors)
     if (op.isUnitary)
       state->applyGate(op.controlQubitIds, op.targetQubitIds, op.deviceData,
@@ -641,14 +653,13 @@ TensorNetState::reverseQubitOrder(std::span<std::complex<double>> stateVec) {
   return ket;
 }
 
-std::unique_ptr<TensorNetState>
-TensorNetState::createFromStateVector(std::span<std::complex<double>> stateVec,
-                                      ScratchDeviceMem &inScratchPad,
-                                      cutensornetHandle_t handle) {
+std::unique_ptr<TensorNetState> TensorNetState::createFromStateVector(
+    std::span<std::complex<double>> stateVec, ScratchDeviceMem &inScratchPad,
+    cutensornetHandle_t handle, std::mt19937 &randomEngine) {
   LOG_API_TIME();
   const std::size_t numQubits = std::log2(stateVec.size());
   auto state =
-      std::make_unique<TensorNetState>(numQubits, inScratchPad, handle);
+      std::make_unique<TensorNetState>(numQubits, inScratchPad, handle, randomEngine);
 
   // Support initializing the tensor network in a specific state vector state.
   // Note: this is not intended for large state vector but for relatively small