ChSonnabend · Oct 7, 2026 · Oct 7, 2026 · Oct 7, 2026 · Oct 8, 2026 · Oct 8, 2026
@@ -19,9 +19,11 @@ ROOT_STANDARD_LIBRARY_PACKAGE(ROOTTMVASofie
   NO_GLOB_HEADERS
   HEADERS
    TMVA/RModel.hxx
+   TMVA/RGPUModel.hxx
    TMVA/SOFIE_common.hxx
   SOURCES
     src/RModel.cxx
+    src/RGPUModel.cxx
     src/SOFIE_common.cxx
     src/SOFIE_common_helpers.cxx
   DEPENDENCIES
@@ -32,3 +34,5 @@ ROOT_STANDARD_LIBRARY_PACKAGE(ROOTTMVASofie
 if(test_tmva_sofie)
   ROOT_ADD_TEST_SUBDIRECTORY(test)
 endif()
+
+target_link_libraries(ROOTTMVASofie PRIVATE ${CMAKE_DL_LIBS})
@@ -191,3 +191,59 @@ parser.CheckModel("example_model.ONNX");
     - [TMVA_SOFIE_RDataFrame_JIT](https://github.com/root-project/root/blob/master/tutorials/machine_learning/TMVA_SOFIE_RDataFrame_JIT.C)
     - [TMVA_SOFIE_RSofieReader](https://github.com/root-project/root/blob/master/tutorials/machine_learning/TMVA_SOFIE_RSofieReader.C)
 
+
+## Experimental CUDA/HIP dense inference
+
+GPU parsing uses SOFIE's existing operator registry and builds an `RModel`.
+`RModel::MakeGPUModel()` runs the existing operator initialization and shape
+inference, then dispatches `ROperator::LowerGPU()`. The existing Gemm and Relu
+operators implement the lowering; other operators reject it explicitly.
+The GPU execution plan currently remains a sequential dense-layer plan.
+FP16 constants are promoted exactly to FP32 for the shared operator representation;
+the execution plan retains FP16 input, output and intermediate tensor precision.
+Direct conversion of an existing FP32 `RModel` is also supported.
+
+
+`RModelParser_ONNX::ParseGPU()` parses sequential, rank-two `Gemm`/`Relu`
+networks with embedded FP32 or FP16 weights. `Gemm` supports constant matrix
+weights, optional vector bias, alpha/beta, transB=0/1 and transA=0. Other graphs
+are rejected rather than falling back to CPU. FP16 uses FP32 accumulation and
+rounds every layer output to FP16. This is an initial scalar-kernel implementation,
+not a tuned BLAS/tensor-core backend.
+
+```cpp
+#include <TMVA/RModelParser_ONNX.hxx>
+
+auto model = TMVA::Experimental::SOFIE::RModelParser_ONNX{}.ParseGPU("network.onnx");
+using Model = TMVA::Experimental::SOFIE::RGPUModel;
+model.Compile({Model::Backend::HIP, "/path/to/hipcc", "gfx90a"});
+auto session = model.CreateSession(nativeStream, deviceId, maxBatch);
+// The framework allocates an aligned device arena of session->WorkspaceSize() bytes.
+session->SetWorkspace(deviceArena, arenaBytes);
+session->Infer(deviceInput, deviceOutput, batch);
+```
+
+Use `Backend::CUDA`, `nvcc`, and the device architecture (e.g. `sm_80`) for CUDA.
+`Compile()` runs the compiler once using an argument vector, loads a private
+shared library, and removes its temporary files. The runtime compiler must be
+available on the execution host; ROOT itself needs no GPU toolkit to build this
+host API. Compilation currently requires a POSIX host. Compiler failures include
+the compiler output. No compilation takes place in `Infer()`.
+
+Streams are passed as opaque native CUDA/HIP stream handles. The current device
+must match the supplied device ID. The arena is 256-byte aligned and contains
+weights plus two intermediate buffers sized for `maxBatch`. Input, output and
+arena must be disjoint. Their precision must match the model. No CUDA/HIP device
+allocation is emitted. `SetWorkspace()` uploads weights on the saved stream;
+`Infer()` enqueues work without synchronizing. GPU faults are observed by the
+caller's normal stream/event synchronization. The framework must finish queued
+work before destroying a session or reclaiming its storage. Separate concurrent
+streams need separate sessions and workspaces; they may share a compiled model.
+
+Tests (not run as part of implementation): with ROOT testing enabled,
+`TestSofieGPU` checks parsing, rejection and capacity handling without a GPU.
+Set `SOFIE_GPU_TEST_BACKEND=CUDA` or `HIP` to additionally build
+`TestSofieGPUDevice`, which needs the corresponding toolkit and compiler. It
+compiles models when the test runs and checks six-layer FP32/FP16 inference,
+external nonblocking streams, several batch sizes, and invalid workspaces.
+The device test skips when no GPU is available.
@@ -0,0 +1,80 @@
+// Author: Christian Sonnabend, 2026
+// For the licensing terms see $ROOTSYS/LICENSE.
+#ifndef TMVA_SOFIE_RGPUMODEL
+#define TMVA_SOFIE_RGPUMODEL
+
+#include <cstddef>
+#include <memory>
+#include <string>
+#include <vector>
+
+namespace TMVA::Experimental::SOFIE {
+
+/// Experimental sequential Gemm/Relu backend. Generated device code has no ROOT dependency.
+class RGPUModel {
+public:
+   enum class Backend {
+      CUDA,
+      HIP
+   };
+   enum class Precision {
+      Float32,
+      Float16
+   };
+   struct CompileOptions {
+      Backend backend = Backend::CUDA;
+      std::string compiler;     ///< Executable path; defaults to nvcc or hipcc.
+      std::string architecture; ///< Required: e.g. sm_80 or gfx90a.
+   };
+   struct Layer {
+      std::size_t inputs = 0, outputs = 0;
+      std::vector<float> weights, bias; ///< Output-major weights, FP32 accumulation for both precisions.
+      float alpha = 1.f, beta = 1.f;
+      bool relu = false;
+   };
+   class Session;
+
+   RGPUModel(Precision precision, std::vector<Layer> layers, std::size_t fixedBatch = 0);
+   std::size_t InputSize() const;
+   std::size_t OutputSize() const;
+   Precision GetPrecision() const { return fPrecision; }
+   std::size_t WorkspaceSize(std::size_t maxBatch) const;
+   static constexpr std::size_t WorkspaceAlignment() { return 256; }
+   std::string GenerateSource(Backend backend) const;
+   void Compile(const CompileOptions &options);
+   std::unique_ptr<Session> CreateSession(void *stream, int device, std::size_t maxBatch) const;
+
+private:
+   struct Program;
+   Precision fPrecision;
+   std::vector<Layer> fLayers;
+   std::vector<float> fWeights;
+   std::size_t fFixedBatch, fWidth = 0;
+   std::shared_ptr<Program> fProgram;
+};
+
+/// Stream, input, output and workspace are borrowed. The caller must complete queued work
+/// before reusing their storage or destroying the session. A session is not thread-safe.
+class RGPUModel::Session {
+public:
+   Session(const Session &) = delete;
+   Session &operator=(const Session &) = delete;
+   ~Session();
+   /// Rebind after the framework allocates/recycles its arena; uploads weights on the saved stream.
+   void SetWorkspace(void *workspace, std::size_t bytes);
+   void Infer(const void *input, void *output, std::size_t batch) const;
+   std::size_t WorkspaceSize() const { return fBytes; }
+
+private:
+   friend class RGPUModel;
+   Session(std::shared_ptr<Program>, std::vector<float>, void *, int, std::size_t, std::size_t, std::size_t,
+           std::size_t, std::size_t);
+   std::shared_ptr<Program> fProgram;
+   std::vector<float> fWeights;
+   void *fStream, *fWorkspace = nullptr;
+   int fDevice;
+   std::size_t fMaxBatch, fFixedBatch, fBytes, fInputBytesPerRow, fOutputBytesPerRow;
+};
+
+} // namespace TMVA::Experimental::SOFIE
+#endif
@@ -1,6 +1,7 @@
 #ifndef TMVA_SOFIE_RMODEL
 #define TMVA_SOFIE_RMODEL
 
+#include "TMVA/RGPUModel.hxx"
 #include "TMVA/SOFIE_common.hxx"
 
 #include "Rtypes.h" // for ClassDefNV
@@ -236,6 +237,7 @@ public:
    template<class T>
    std::vector<T> GetTensorData(const std::string & name);
 
+   RGPUModel MakeGPUModel(RGPUModel::Precision precision, std::size_t fixedBatch = 0);
    void Initialize(int batchSize = -1, bool verbose = false);
    void Initialize(const std::map<std::string,size_t> & inputParams, bool verbose = false);
 

@@ -2,6 +2,8 @@
 #define TMVA_SOFIE_ROPERATOR
 
 #include "TMVA/SOFIE_common.hxx"
+#include "TMVA/RGPUModel.hxx"
+#include <stdexcept>
 
 #include <ROOT/RSpan.hxx>
 
@@ -19,6 +21,10 @@ public:
    virtual std::vector<std::string> GetBlasRoutines() { return {}; }
    virtual std::vector<std::string> GetStdLibs() { return {}; }
    virtual void Initialize(RModel&) = 0;
+   virtual void LowerGPU(RModel &, std::vector<RGPUModel::Layer> &) const
+   {
+      throw std::runtime_error("SOFIE GPU: operator has no GPU lowering");
+   }
    virtual std::string Generate(std::string OpName) = 0;  //expect unique opName for each operator within the same RModel
    // generate initialization code for session constructor
    virtual std::string GenerateInitCode() { return "";}

@@ -147,7 +147,39 @@ namespace SOFIE{
          return DoShapeInference<Dim>(input);
       }
 
-
+      void LowerGPU(RModel &model, std::vector<RGPUModel::Layer> &layers) const override
+      {
+         if (model.GetTensorType(fNB) != ETensorType::FLOAT ||
+             (!fNC.empty() && model.GetTensorType(fNC) != ETensorType::FLOAT))
+            throw std::runtime_error("SOFIE GPU: lowering requires FP32 weight storage");
+         const auto a = model.GetTensorShape(fNA);
+         const auto b = model.GetTensorShape(fNB);
+         if (a.size() != 2 || b.size() != 2 || fAttrTransA || !model.IsInitializedTensor(fNB))
+            throw std::runtime_error("SOFIE GPU: Gemm requires rank-two input and constant weights");
+         RGPUModel::Layer layer;
+         layer.inputs = a[1];
+         layer.outputs = b[fAttrTransB ? 0 : 1];
+         if (b[fAttrTransB ? 1 : 0] != layer.inputs)
+            throw std::runtime_error("SOFIE GPU: inconsistent Gemm dimensions");
+         layer.alpha = fAttrAlpha;
+         layer.beta = fAttrBeta;
+         layer.relu = fActivation == EActivationType::RELU;
+         const auto *weights = static_cast<const float *>(model.GetInitializedTensorData(fNB).get());
+         layer.weights.resize(layer.inputs * layer.outputs);
+         for (size_t n = 0; n < layer.outputs; ++n)
+            for (size_t k = 0; k < layer.inputs; ++k)
+               layer.weights[n * layer.inputs + k] =
+                  weights[fAttrTransB ? n * layer.inputs + k : k * layer.outputs + n];
+         layer.bias.assign(layer.outputs, 0.f);
+         if (!fNC.empty()) {
+            const auto c = model.GetTensorShape(fNC);
+            if (!model.IsInitializedTensor(fNC) || c != std::vector<size_t>{layer.outputs})
+               throw std::runtime_error("SOFIE GPU: Gemm requires constant vector bias");
+            const auto *bias = static_cast<const float *>(model.GetInitializedTensorData(fNC).get());
+            layer.bias.assign(bias, bias + layer.outputs);
+         }
+         layers.push_back(std::move(layer));
+      }
 
       void Initialize(RModel& model) override {
          //TODO: propagate A or B as specified by ONNX standard

@@ -29,6 +29,13 @@ public:
          fOutputTensorNames = { fNY };
       }
 
+      void LowerGPU(RModel &, std::vector<RGPUModel::Layer> &layers) const override
+      {
+         if (layers.empty() || layers.back().relu)
+            throw std::runtime_error("SOFIE GPU: Relu must follow Gemm");
+         layers.back().relu = true;
+      }
+
    void Initialize(RModel& model) override {
       if (model.CheckIfTensorAlreadyExist(fNX) == false){   //input must be a graph input, or already initialized intermediate tensor
          throw std::runtime_error("TMVA SOFIE Relu Op Input Tensor " + fNX + " is not found in model");