[WIP] Operators and Integration with cudensitymat

CMakeLists.txt

@@ -151,6 +151,9 @@ endif()
 if(NOT CUTENSORNET_ROOT)
   SET(CUTENSORNET_ROOT "$ENV{CUQUANTUM_INSTALL_PREFIX}")
 endif()
+if(NOT CUDENSITYMAT_ROOT)
+  SET(CUDENSITYMAT_ROOT "$ENV{CUQUANTUM_INSTALL_PREFIX}")
+endif()
 if(NOT CUTENSOR_ROOT)
   SET(CUTENSOR_ROOT "$ENV{CUTENSOR_INSTALL_PREFIX}")
 endif()

docs/sphinx/api/languages/cpp_api.rst

@@ -223,6 +223,53 @@ Utilities
 .. doxygentypedef:: cudaq::real 
 
 .. doxygenfunction:: cudaq::range(std::size_t)
+
+Dynamics 
+=========
+
+.. .. doxygenclass:: cudaq::EvolveResult
+    :members:
+
+.. .. doxygenclass:: cudaq::AsyncEvolveResult
+    :members:
+
+.. doxygenclass:: cudaq::operator_sum
+    :members:
+
+.. doxygenclass:: cudaq::product_operator
+    :members:
+
+.. doxygenclass:: cudaq::scalar_operator
+    :members:
+
+.. doxygenclass:: cudaq::matrix_operator
+    :members:
+
+.. doxygenclass:: cudaq::OperatorArithmetics
+    :members:
+
+.. doxygenclass:: cudaq::MatrixArithmetics
+    :members:
+
+.. doxygenclass:: cudaq::Schedule
+    :members:
+
+.. doxygenclass:: cudaq::operators
+    :members:
+
+.. doxygenclass:: cudaq::pauli
+    :members:
+
+.. .. doxygenfunction:: cudaq::evolve(Operator hamiltonian, std::map<int,int> dimensions, Schedule schedule, bool store_intermediate_states)
+
+.. .. doxygenfunction:: cudaq::evolve(Operator hamiltonian, std::map<int,int> dimensions, Schedule schedule, std::vector<Operator> collapse_operators, std::vector<Operator> observables, bool store_intermediate_states)
+
+.. .. doxygenfunction:: cudaq::evolve(Operator hamiltonian, std::map<int,int> dimensions, Schedule schedule, state initial_state, std::vector<Operator> collapse_operators, std::vector<Operator> observables, bool store_intermediate_states)
+
+.. .. doxygenfunction:: cudaq::evolve(Operator hamiltonian, std::map<int,int> dimensions, Schedule schedule, std::vector<state> initial_states, std::vector<Operator> collapse_operators, std::vector<Operator> observables, bool store_intermediate_states)
+
+.. .. doxygenfunction:: cudaq::evolve_async
+
 
 Namespaces 
 ===========

runtime/common/EvolveResult.h

@@ -108,8 +108,8 @@ class evolve_result {
     return final_expectation_values;
   }
 
-  std::optional<std::vector<std::vector<observe_result>>>
-  get_expectation_values() {
+  const std::optional<std::vector<std::vector<observe_result>>> &
+  get_expectation_values() const {
     return expectation_values;
   }
 

runtime/cudaq/CMakeLists.txt

@@ -40,7 +40,7 @@ if (CUDA_FOUND)
     PRIVATE .)
 
   target_link_libraries(${LIBRARY_NAME}
-    PUBLIC dl cudaq-spin cudaq-common cudaq-nlopt cudaq-ensmallen
+    PUBLIC dl cudaq-spin cudaq-operators cudaq-common cudaq-nlopt cudaq-ensmallen
     PRIVATE nvqir fmt::fmt-header-only CUDA::cudart_static)
 
   target_compile_definitions(${LIBRARY_NAME} PRIVATE CUDAQ_HAS_CUDA)
@@ -52,7 +52,7 @@ else()
     PRIVATE .)
 
   target_link_libraries(${LIBRARY_NAME}
-    PUBLIC dl cudaq-spin cudaq-common cudaq-nlopt cudaq-ensmallen
+    PUBLIC dl cudaq-spin cudaq-operators cudaq-common cudaq-nlopt cudaq-ensmallen
     PRIVATE nvqir fmt::fmt-header-only)
 endif()
 
@@ -61,6 +61,7 @@ add_subdirectory(algorithms)
 add_subdirectory(platform)
 add_subdirectory(builder)
 add_subdirectory(domains)
+add_subdirectory(dynamics)
 
 install(TARGETS ${LIBRARY_NAME} EXPORT cudaq-targets DESTINATION lib)
 

runtime/cudaq/base_integrator.h

@@ -0,0 +1,81 @@
+/****************************************************************-*- C++ -*-****
+ * Copyright (c) 2022 - 2025 NVIDIA Corporation & Affiliates.                  *
+ * All rights reserved.                                                        *
+ *                                                                             *
+ * This source code and the accompanying materials are made available under    *
+ * the terms of the Apache License 2.0 which accompanies this distribution.    *
+ ******************************************************************************/
+
+#pragma once
+
+#include "base_time_stepper.h"
+#include "operators.h"
+#include "schedule.h"
+#include <map>
+#include <memory>
+#include <vector>
+
+namespace cudaq {
+template <typename TState>
+class BaseIntegrator {
+protected:
+  std::map<std::string, double> integrator_options;
+  TState state;
+  double t;
+  // std::map<int, int> dimensions;
+  // std::shared_ptr<Schedule> schedule;
+  // std::shared_ptr<operator_sum<HandlerTy>> hamiltonian;
+  std::shared_ptr<BaseTimeStepper<TState>> stepper;
+  // std::vector<std::shared_ptr<operator_sum<HandlerTy>>> collapse_operators;
+
+  virtual void post_init() = 0;
+
+public:
+  /// @brief Default constructor
+  BaseIntegrator() = default;
+
+  /// @brief Constructor to initialize the integrator with a state and time
+  /// stepper.
+  /// @param initial_state Initial quantum state.
+  /// @param t0 Initial time.
+  /// @param stepper Time stepper instance.
+  BaseIntegrator(TState &&initial_state, double t0,
+                 std::shared_ptr<BaseTimeStepper<TState>> stepper)
+      : state(std::move(initial_state)), t(t0), stepper(std::move(stepper)) {
+    if (!this->stepper) {
+      throw std::runtime_error("Time stepper is not initialized.");
+    }
+  }
+
+  virtual ~BaseIntegrator() = default;
+
+  /// @brief Set the initial state and time
+  void set_state(const TState &initial_state, double t0 = 0.0) {
+    state = initial_state;
+    t = t0;
+  }
+
+  /// @brief Set an option for the integrator
+  void set_option(const std::string &key, double value) {
+    integrator_options[key] = value;
+  }
+
+  /// @brief Set the system parameters (dimensions, schedule, and operators)
+  // void set_system(const std::map<int, int> &dimensions,
+  //                 std::shared_ptr<Schedule> schedule,
+  //                 std::shared_ptr<operator_sum<HandlerTy>> hamiltonian,
+  //                 std::vector<std::shared_ptr<operator_sum<HandlerTy>>>
+  //                     collapse_operators = {}) {
+  //   this->dimensions = dimensions;
+  //   this->schedule = schedule;
+  //   this->hamiltonian = hamiltonian;
+  //   this->collapse_operators = collapse_operators;
+  // }
+
+  /// @brief Perform integration to the target time.
+  virtual void integrate(double target_time) = 0;
+
+  /// @brief Get the current time and state.
+  std::pair<double, const TState &> get_state() const { return {t, state}; }
+};
+} // namespace cudaq

runtime/cudaq/base_time_stepper.h

@@ -0,0 +1,24 @@
+/****************************************************************-*- C++ -*-****
+ * Copyright (c) 2022 - 2025 NVIDIA Corporation & Affiliates.                  *
+ * All rights reserved.                                                        *
+ *                                                                             *
+ * This source code and the accompanying materials are made available under    *
+ * the terms of the Apache License 2.0 which accompanies this distribution.    *
+ ******************************************************************************/
+
+#pragma once
+
+namespace cudaq {
+template <typename TState>
+class BaseTimeStepper {
+public:
+  virtual ~BaseTimeStepper() = default;
+
+  /// @brief Compute the next time step for the given quantum state.
+  /// @param state The quantum state to evolve.
+  /// @param t Current time.
+  /// @param step_size Time step size.
+  /// @return The updated quantum state after stepping.
+  virtual TState compute(TState &state, double t, double step_size) = 0;
+};
+} // namespace cudaq

runtime/cudaq/cudm_error_handling.h

@@ -0,0 +1,30 @@
+/****************************************************************-*- C++ -*-****
+ * Copyright (c) 2022 - 2025 NVIDIA Corporation & Affiliates.                  *
+ * All rights reserved.                                                        *
+ *                                                                             *
+ * This source code and the accompanying materials are made available under    *
+ * the terms of the Apache License 2.0 which accompanies this distribution.    *
+ ******************************************************************************/
+
+#pragma once
+#include <cudensitymat.h>
+#include <fmt/core.h>
+#include <stdexcept>
+
+#define HANDLE_CUDM_ERROR(x)                                                   \
+  {                                                                            \
+    const auto err = x;                                                        \
+    if (err != CUDENSITYMAT_STATUS_SUCCESS) {                                  \
+      throw std::runtime_error(fmt::format("[cudaq] %{} in {} (line {})", err, \
+                                           __FUNCTION__, __LINE__));           \
+    }                                                                          \
+  }
+
+#define HANDLE_CUDA_ERROR(x)                                                   \
+  {                                                                            \
+    const auto err = x;                                                        \
+    if (err != cudaSuccess) {                                                  \
+      throw std::runtime_error(fmt::format("[cuda] %{} in {} (line {})", err,  \
+                                           __FUNCTION__, __LINE__));           \
+    }                                                                          \
+  }

runtime/cudaq/cudm_expectation.h

@@ -0,0 +1,35 @@
+/****************************************************************-*- C++ -*-****
+ * Copyright (c) 2022 - 2025 NVIDIA Corporation & Affiliates.                  *
+ * All rights reserved.                                                        *
+ *                                                                             *
+ * This source code and the accompanying materials are made available under    *
+ * the terms of the Apache License 2.0 which accompanies this distribution.    *
+ ******************************************************************************/
+
+#pragma once
+#include <complex>
+#include <cudensitymat.h>
+namespace cudaq {
+
+class cudm_expectation {
+  cudensitymatHandle_t m_handle{nullptr};
+  cudensitymatOperator_t m_hamOp{nullptr};
+  cudensitymatExpectation_t m_expectation{nullptr};
+  cudensitymatWorkspaceDescriptor_t m_workspace{nullptr};
+
+public:
+  cudm_expectation(cudensitymatHandle_t handle, cudensitymatOperator_t op);
+  cudm_expectation(const cudm_expectation &) = delete;
+  cudm_expectation &operator=(const cudm_expectation &) = delete;
+  cudm_expectation(cudm_expectation &&src) {
+    std::swap(m_handle, src.m_handle);
+    std::swap(m_hamOp, src.m_hamOp);
+    std::swap(m_expectation, src.m_expectation);
+    std::swap(m_workspace, src.m_workspace);
+  }
+  ~cudm_expectation();
+  void prepare(cudensitymatState_t state);
+  std::complex<double> compute(cudensitymatState_t state, double time);
+};
+
+} // namespace cudaq

runtime/cudaq/cudm_helpers.h

@@ -0,0 +1,98 @@
+/****************************************************************-*- C++ -*-****
+ * Copyright (c) 2022 - 2025 NVIDIA Corporation & Affiliates.                  *
+ * All rights reserved.                                                        *
+ *                                                                             *
+ * This source code and the accompanying materials are made available under    *
+ * the terms of the Apache License 2.0 which accompanies this distribution.    *
+ ******************************************************************************/
+
+#pragma once
+
+#include "cudaq/operators.h"
+#include "cudaq/utils/tensor.h"
+#include <cudensitymat.h>
+#include <iostream>
+#include <map>
+#include <stdexcept>
+#include <vector>
+
+namespace cudaq {
+class cudm_helper {
+public:
+  explicit cudm_helper(cudensitymatHandle_t handle);
+  ~cudm_helper();
+
+  // Matrix flattening
+  static std::vector<std::complex<double>>
+  flatten_matrix(const matrix_2 &matrix);
+
+  // State Operations
+  void scale_state(cudensitymatState_t state, double scale_factor,
+                   cudaStream_t stream);
+
+  // Compute Lindblad Operator
+  cudensitymatOperator_t
+  compute_lindblad_operator(const std::vector<matrix_2> &c_ops,
+                            const std::vector<int64_t> &mode_extents);
+
+  // Convert operator sum to cudensitymat operator
+  template <typename HandlerTy>
+  cudensitymatOperator_t convert_to_cudensitymat_operator(
+      const std::map<std::string, std::complex<double>> &parameters,
+      const operator_sum<HandlerTy> &op,
+      const std::vector<int64_t> &mode_extents);
+
+  // Construct Liouvillian
+  cudensitymatOperator_t construct_liouvillian(
+      const operator_sum<cudaq::matrix_operator> &op,
+      const std::vector<operator_sum<cudaq::matrix_operator> *>
+          &collapse_operators,
+      const std::vector<int64_t> &mode_extents,
+      const std::map<std::string, std::complex<double>> &parameters,
+      bool is_master_equation);
+
+  // Construct Liouvillian
+  cudensitymatOperator_t construct_liouvillian(
+      const cudensitymatOperator_t &hamiltonian,
+      const std::vector<cudensitymatOperator_t> &collapse_operators,
+      double gamma);
+
+  // Helper Functions
+  std::map<int, int>
+  convert_dimensions(const std::vector<int64_t> &mode_extents);
+  std::vector<int64_t>
+  get_subspace_extents(const std::vector<int64_t> &mode_extents,
+                       const std::vector<int> &degrees);
+
+  // Callback Wrappers
+  static cudensitymatWrappedScalarCallback_t
+  _wrap_callback(const scalar_operator &scalar_op);
+  static cudensitymatWrappedTensorCallback_t
+  _wrap_tensor_callback(const matrix_operator &op);
+
+  // Elementary Operator Functions
+  void append_scalar_to_term(cudensitymatOperatorTerm_t term,
+                             const scalar_operator &scalar_op);
+  cudensitymatElementaryOperator_t create_elementary_operator(
+      const cudaq::matrix_operator &elem_op,
+      const std::map<std::string, std::complex<double>> &parameters,
+      const std::vector<int64_t> &mode_extents);
+  void append_elementary_operator_to_term(
+      cudensitymatOperatorTerm_t term,
+      const std::vector<cudensitymatElementaryOperator_t> &elem_ops,
+      const std::vector<std::vector<int>> &degrees, bool is_dagger);
+
+  // GPU memory management
+  static void *
+  create_array_gpu(const std::vector<std::complex<double>> &cpu_array);
+  static void destroy_array_gpu(void *gpu_array);
+
+private:
+  cudensitymatHandle_t handle;
+};
+
+extern template cudensitymatOperator_t
+cudm_helper::convert_to_cudensitymat_operator<cudaq::matrix_operator>(
+    const std::map<std::string, std::complex<double>> &,
+    const operator_sum<cudaq::matrix_operator> &, const std::vector<int64_t> &);
+} // namespace cudaq