diff --git a/tmva/sofie/CMakeLists.txt b/tmva/sofie/CMakeLists.txt
index 84ddfdc7092d3..3232a9a90e976 100644
--- a/tmva/sofie/CMakeLists.txt
+++ b/tmva/sofie/CMakeLists.txt
@@ -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})
diff --git a/tmva/sofie/README.md b/tmva/sofie/README.md
index 144944e8f58e7..021b833a02c68 100644
--- a/tmva/sofie/README.md
+++ b/tmva/sofie/README.md
@@ -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.
diff --git a/tmva/sofie/inc/TMVA/RGPUModel.hxx b/tmva/sofie/inc/TMVA/RGPUModel.hxx
new file mode 100644
index 0000000000000..d8e4d8c055f77
--- /dev/null
+++ b/tmva/sofie/inc/TMVA/RGPUModel.hxx
@@ -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
diff --git a/tmva/sofie/inc/TMVA/RModel.hxx b/tmva/sofie/inc/TMVA/RModel.hxx
index dc1b604dc76f3..46cd82eddd0c4 100644
--- a/tmva/sofie/inc/TMVA/RModel.hxx
+++ b/tmva/sofie/inc/TMVA/RModel.hxx
@@ -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);
 
diff --git a/tmva/sofie/inc/TMVA/ROperator.hxx b/tmva/sofie/inc/TMVA/ROperator.hxx
index 13c501b0905d6..08d5df49c9fa4 100644
--- a/tmva/sofie/inc/TMVA/ROperator.hxx
+++ b/tmva/sofie/inc/TMVA/ROperator.hxx
@@ -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 "";}
diff --git a/tmva/sofie/inc/TMVA/ROperator_Gemm.hxx b/tmva/sofie/inc/TMVA/ROperator_Gemm.hxx
index e177ef97605ef..49b1ba3332f4f 100644
--- a/tmva/sofie/inc/TMVA/ROperator_Gemm.hxx
+++ b/tmva/sofie/inc/TMVA/ROperator_Gemm.hxx
@@ -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
diff --git a/tmva/sofie/inc/TMVA/ROperator_Relu.hxx b/tmva/sofie/inc/TMVA/ROperator_Relu.hxx
index 3bd959de6f726..c9fea36273134 100644
--- a/tmva/sofie/inc/TMVA/ROperator_Relu.hxx
+++ b/tmva/sofie/inc/TMVA/ROperator_Relu.hxx
@@ -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");
diff --git a/tmva/sofie/src/RGPUModel.cxx b/tmva/sofie/src/RGPUModel.cxx
new file mode 100644
index 0000000000000..3c033444fd307
--- /dev/null
+++ b/tmva/sofie/src/RGPUModel.cxx
@@ -0,0 +1,304 @@
+// Author: Christian Sonnabend, 2026
+// For the licensing terms see $ROOTSYS/LICENSE.
+#include "TMVA/RGPUModel.hxx"
+
+#include <algorithm>
+#include <cerrno>
+#include <cmath>
+#include <cstdint>
+#include <cstdlib>
+#include <filesystem>
+#include <fstream>
+#include <iomanip>
+#include <limits>
+#include <locale>
+#include <sstream>
+#include <stdexcept>
+#include <utility>
+#ifndef _WIN32
+#include <dlfcn.h>
+#include <fcntl.h>
+#include <spawn.h>
+#include <sys/wait.h>
+#include <unistd.h>
+extern char **environ;
+#endif
+
+namespace TMVA::Experimental::SOFIE {
+namespace {
+std::size_t Multiply(std::size_t a, std::size_t b)
+{
+   if (b && a > std::numeric_limits<std::size_t>::max() / b)
+      throw std::overflow_error("SOFIE GPU: tensor size overflow");
+   return a * b;
+}
+std::size_t Add(std::size_t a, std::size_t b)
+{
+   if (a > std::numeric_limits<std::size_t>::max() - b)
+      throw std::overflow_error("SOFIE GPU: workspace size overflow");
+   return a + b;
+}
+std::size_t Align(std::size_t bytes)
+{
+   return Add(bytes, 255) & ~std::size_t(255);
+}
+} // namespace
+
+struct RGPUModel::Program {
+   using Upload = const char *(*)(void *, int, void *, const float *, std::size_t);
+   using Invoke = const char *(*)(void *, int, void *, const void *, void *, std::size_t, std::size_t);
+   void *library = nullptr;
+   Upload upload = nullptr;
+   Invoke invoke = nullptr;
+   ~Program()
+   {
+#ifndef _WIN32
+      if (library)
+         dlclose(library);
+#endif
+   }
+};
+
+RGPUModel::RGPUModel(Precision precision, std::vector<Layer> layers, std::size_t fixedBatch)
+   : fPrecision(precision), fLayers(std::move(layers)), fFixedBatch(fixedBatch)
+{
+   if (fLayers.empty())
+      throw std::invalid_argument("SOFIE GPU: empty network");
+   if (precision != Precision::Float32 && precision != Precision::Float16)
+      throw std::invalid_argument("SOFIE GPU: unsupported precision");
+   std::size_t previous = fLayers.front().inputs;
+   for (const auto &layer : fLayers) {
+      if (!layer.inputs || !layer.outputs || previous != layer.inputs ||
+          layer.weights.size() != Multiply(layer.inputs, layer.outputs) || layer.bias.size() != layer.outputs ||
+          !std::isfinite(layer.alpha) || !std::isfinite(layer.beta))
+         throw std::invalid_argument("SOFIE GPU: invalid dense layer");
+      fWeights.insert(fWeights.end(), layer.weights.begin(), layer.weights.end());
+      fWeights.insert(fWeights.end(), layer.bias.begin(), layer.bias.end());
+      fWidth = std::max(fWidth, layer.outputs);
+      previous = layer.outputs;
+   }
+}
+std::size_t RGPUModel::InputSize() const
+{
+   return fLayers.front().inputs;
+}
+std::size_t RGPUModel::OutputSize() const
+{
+   return fLayers.back().outputs;
+}
+std::size_t RGPUModel::WorkspaceSize(std::size_t maxBatch) const
+{
+   if (!maxBatch || (fFixedBatch && maxBatch != fFixedBatch))
+      throw std::invalid_argument("SOFIE GPU: invalid maximum batch size");
+   const auto tensor = Align(Multiply(Multiply(maxBatch, fWidth), fPrecision == Precision::Float16 ? 2 : 4));
+   return Add(Align(Multiply(fWeights.size(), sizeof(float))), Multiply(2, tensor));
+}
+
+std::string RGPUModel::GenerateSource(Backend backend) const
+{
+   if (backend != Backend::CUDA && backend != Backend::HIP)
+      throw std::invalid_argument("SOFIE GPU: unsupported backend");
+   const std::string api = backend == Backend::HIP ? "hip" : "cuda";
+   std::ostringstream out;
+   out.imbue(std::locale::classic());
+   out << std::scientific << std::setprecision(std::numeric_limits<float>::max_digits10);
+   out << "#include <" << (backend == Backend::HIP ? "hip/hip_runtime.h" : "cuda_runtime.h") << ">\n"
+       << "#include <" << (backend == Backend::HIP ? "hip/hip_fp16.h" : "cuda_fp16.h") << ">\n"
+       << "#include <cstddef>\n#include <cstdint>\n#include <cstdlib>\n"
+       << "using Scalar = " << (fPrecision == Precision::Float16 ? "__half" : "float") << ";\n"
+       << "using Stream = " << api << "Stream_t;\n";
+   out << R"cpp(
+__device__ float read(float x) { return x; }
+__device__ float read(__half x) { return __half2float(x); }
+__device__ void write(float &x, float y) { x = y; }
+__device__ void write(__half &x, float y) { x = __float2half_rn(y); }
+template <size_t K, size_t N, bool Relu>
+__global__ void dense(const Scalar *input, const float *weights, const float *bias,
+                      Scalar *output, size_t batch, float alpha, float beta)
+{
+   for (size_t i = size_t(blockIdx.x) * blockDim.x + threadIdx.x; i < batch * N;
+        i += size_t(blockDim.x) * gridDim.x) {
+      const size_t row = i / N, col = i % N;
+      float sum = 0.f;
+      for (size_t k = 0; k < K; ++k) sum += read(input[row * K + k]) * weights[col * K + k];
+      Scalar value;
+      write(value, alpha * sum + beta * bias[col]);
+      if (Relu && read(value) < 0.f) write(value, 0.f);
+      output[i] = value;
+   }
+}
+)cpp";
+   out << "static const char *checkDevice(int expected) { int actual = -1; auto e = " << api
+       << "GetDevice(&actual); if (e != " << api << "Success) return " << api
+       << "GetErrorString(e); return actual == expected ? nullptr : \"SOFIE GPU: wrong active device\"; }\n";
+   out << "extern \"C\" const char *sofie_upload(void *s, int device, void *arena, const float *weights, size_t bytes) "
+          "{\n"
+       << "if (auto e = checkDevice(device)) return e;\n"
+       << "auto e = " << api << "MemcpyAsync(arena, weights, bytes, " << api
+       << "MemcpyHostToDevice, static_cast<Stream>(s));\n"
+       << "return e == " << api << "Success ? nullptr : " << api << "GetErrorString(e); }\n";
+   out << "extern \"C\" const char *sofie_infer(void *s, int device, void *arena, const void *input, void *output, "
+          "size_t batch, size_t maxBatch) {\n"
+       << "if (auto e = checkDevice(device)) return e;\n"
+       << "auto stream = static_cast<Stream>(s); auto weights = static_cast<const float *>(arena);\n"
+       << "size_t stride = (maxBatch * " << fWidth << " * sizeof(Scalar) + 255) & ~size_t(255);\n"
+       << "auto a = reinterpret_cast<Scalar *>(static_cast<char *>(arena) + " << Align(fWeights.size() * sizeof(float))
+       << ");\n"
+       << "auto b = reinterpret_cast<Scalar *>(reinterpret_cast<char *>(a) + stride);\n";
+   std::size_t offset = 0;
+   for (std::size_t i = 0; i < fLayers.size(); ++i) {
+      const auto &l = fLayers[i];
+      const auto count = l.inputs * l.outputs;
+      out << "dense<" << l.inputs << ", " << l.outputs << ", " << (l.relu ? "true" : "false")
+          << "><<<unsigned((batch * " << l.outputs << " + 255) / 256 > 65535 ? 65535 : (batch * " << l.outputs
+          << " + 255) / 256), 256, 0, stream>>>("
+          << (i == 0 ? "static_cast<const Scalar *>(input)" : (i % 2 ? "a" : "b")) << ", weights + " << offset
+          << ", weights + " << offset + count << ", "
+          << (i + 1 == fLayers.size() ? "static_cast<Scalar *>(output)" : (i % 2 ? "b" : "a")) << ", batch, " << l.alpha
+          << "f, " << l.beta << "f);\n"
+          << "if (auto e = " << api << "GetLastError(); e != " << api << "Success) return " << api
+          << "GetErrorString(e);\n";
+      offset += count + l.outputs;
+   }
+   out << "return nullptr;\n}\n";
+   return out.str();
+}
+
+void RGPUModel::Compile(const CompileOptions &options)
+{
+   if (fProgram)
+      throw std::logic_error("SOFIE GPU: model already compiled");
+   const auto source = GenerateSource(options.backend);
+   if (options.architecture.empty() ||
+       options.architecture.find_first_not_of("abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ0123456789_:+-") !=
+          std::string::npos)
+      throw std::invalid_argument("SOFIE GPU: provide a CUDA/HIP architecture");
+#ifndef _WIN32
+   std::string pattern = (std::filesystem::temp_directory_path() / "sofie-gpu-XXXXXX").string();
+   if (!mkdtemp(pattern.data()))
+      throw std::runtime_error("SOFIE GPU: cannot create compilation directory");
+   struct TemporaryDirectory {
+      std::filesystem::path path;
+      ~TemporaryDirectory()
+      {
+         std::error_code e;
+         std::filesystem::remove_all(path, e);
+      }
+   } directory{pattern};
+   const auto input = directory.path / "model.cu", library = directory.path / "model.so",
+              log = directory.path / "compiler.log";
+   {
+      std::ofstream file(input);
+      file << source;
+      if (!file)
+         throw std::runtime_error("SOFIE GPU: cannot write source");
+   }
+   const bool hip = options.backend == Backend::HIP;
+   std::vector<std::string> args{options.compiler.empty() ? (hip ? "hipcc" : "nvcc") : options.compiler, "-std=c++17",
+                                 "-O2", "-shared"};
+   if (hip) {
+      args.emplace_back("-fPIC");
+      args.emplace_back("--offload-arch=" + options.architecture);
+   } else {
+      args.emplace_back("-Xcompiler=-fPIC");
+      args.emplace_back("-arch=" + options.architecture);
+   }
+   args.insert(args.end(), {input.string(), "-o", library.string()});
+   std::vector<char *> argv;
+   for (auto &arg : args)
+      argv.push_back(arg.data());
+   argv.push_back(nullptr);
+   posix_spawn_file_actions_t actions;
+   if (posix_spawn_file_actions_init(&actions))
+      throw std::runtime_error("SOFIE GPU: cannot initialize compiler process");
+   int error =
+      posix_spawn_file_actions_addopen(&actions, STDOUT_FILENO, log.c_str(), O_WRONLY | O_CREAT | O_TRUNC, 0600);
+   if (!error)
+      error = posix_spawn_file_actions_adddup2(&actions, STDOUT_FILENO, STDERR_FILENO);
+   pid_t pid = 0;
+   if (!error)
+      error = posix_spawnp(&pid, argv[0], &actions, nullptr, argv.data(), environ);
+   posix_spawn_file_actions_destroy(&actions);
+   if (error)
+      throw std::runtime_error("SOFIE GPU: cannot start compiler " + args[0] + " (error " + std::to_string(error) +
+                               ")");
+   int status = 0;
+   pid_t waited;
+   do {
+      waited = waitpid(pid, &status, 0);
+   } while (waited < 0 && errno == EINTR);
+   if (waited < 0 || !WIFEXITED(status) || WEXITSTATUS(status)) {
+      std::ifstream file(log);
+      std::ostringstream diagnostic;
+      diagnostic << file.rdbuf();
+      throw std::runtime_error("SOFIE GPU: compilation failed\n" + diagnostic.str());
+   }
+   auto program = std::make_shared<Program>();
+   program->library = dlopen(library.c_str(), RTLD_NOW | RTLD_LOCAL);
+   if (!program->library)
+      throw std::runtime_error(std::string("SOFIE GPU: cannot load model: ") + dlerror());
+   program->upload = reinterpret_cast<Program::Upload>(dlsym(program->library, "sofie_upload"));
+   program->invoke = reinterpret_cast<Program::Invoke>(dlsym(program->library, "sofie_infer"));
+   if (!program->upload || !program->invoke)
+      throw std::runtime_error("SOFIE GPU: missing generated entry point");
+   fProgram = std::move(program);
+#else
+   throw std::runtime_error("SOFIE GPU: runtime compilation currently requires a POSIX host");
+#endif
+}
+
+RGPUModel::Session::Session(std::shared_ptr<Program> program, std::vector<float> weights, void *stream, int device,
+                            std::size_t maxBatch, std::size_t fixedBatch, std::size_t bytes, std::size_t inputBytes,
+                            std::size_t outputBytes)
+   : fProgram(std::move(program)),
+     fWeights(std::move(weights)),
+     fStream(stream),
+     fDevice(device),
+     fMaxBatch(maxBatch),
+     fFixedBatch(fixedBatch),
+     fBytes(bytes),
+     fInputBytesPerRow(inputBytes),
+     fOutputBytesPerRow(outputBytes)
+{
+   if (device < 0)
+      throw std::invalid_argument("SOFIE GPU: invalid device index");
+}
+RGPUModel::Session::~Session() = default;
+std::unique_ptr<RGPUModel::Session> RGPUModel::CreateSession(void *stream, int device, std::size_t maxBatch) const
+{
+   if (!fProgram)
+      throw std::logic_error("SOFIE GPU: compile the model before creating a session");
+   const auto elementSize = fPrecision == Precision::Float16 ? 2u : 4u;
+   return std::unique_ptr<Session>(new Session(fProgram, fWeights, stream, device, maxBatch, fFixedBatch,
+                                               WorkspaceSize(maxBatch), Multiply(InputSize(), elementSize),
+                                               Multiply(OutputSize(), elementSize)));
+}
+void RGPUModel::Session::SetWorkspace(void *workspace, std::size_t bytes)
+{
+   if (!workspace || reinterpret_cast<std::uintptr_t>(workspace) % WorkspaceAlignment() || bytes < fBytes)
+      throw std::invalid_argument("SOFIE GPU: missing, undersized or unaligned workspace");
+   fWorkspace = nullptr;
+   if (auto error = fProgram->upload(fStream, fDevice, workspace, fWeights.data(), fWeights.size() * sizeof(float)))
+      throw std::runtime_error(error);
+   fWorkspace = workspace;
+}
+void RGPUModel::Session::Infer(const void *input, void *output, std::size_t batch) const
+{
+   if (batch > fMaxBatch || (batch && fFixedBatch && batch != fFixedBatch))
+      throw std::invalid_argument("SOFIE GPU: batch outside the prepared range");
+   if (!batch)
+      return;
+   if (!fWorkspace || !input || !output || input == output)
+      throw std::invalid_argument("SOFIE GPU: provide workspace and distinct input/output buffers");
+   auto overlaps = [](const void *a, std::size_t na, const void *b, std::size_t nb) {
+      const auto x = reinterpret_cast<std::uintptr_t>(a), y = reinterpret_cast<std::uintptr_t>(b);
+      return x <= y ? y - x < na : x - y < nb;
+   };
+   const auto inputBytes = Multiply(batch, fInputBytesPerRow), outputBytes = Multiply(batch, fOutputBytesPerRow);
+   if (overlaps(input, inputBytes, output, outputBytes) || overlaps(input, inputBytes, fWorkspace, fBytes) ||
+       overlaps(output, outputBytes, fWorkspace, fBytes))
+      throw std::invalid_argument("SOFIE GPU: input, output and workspace must not overlap");
+   if (auto error = fProgram->invoke(fStream, fDevice, fWorkspace, input, output, batch, fMaxBatch))
+      throw std::runtime_error(error);
+}
+} // namespace TMVA::Experimental::SOFIE
diff --git a/tmva/sofie/src/RModel.cxx b/tmva/sofie/src/RModel.cxx
index 22bd4318044ef..0b6761cc8ab10 100644
--- a/tmva/sofie/src/RModel.cxx
+++ b/tmva/sofie/src/RModel.cxx
@@ -659,6 +659,26 @@ void RModel::CheckAndFlushIntermediateMemory(std::span<const std::string_view> o
    }
 }
 
+RGPUModel RModel::MakeGPUModel(RGPUModel::Precision precision, std::size_t fixedBatch)
+{
+   if (fInputTensorNames.size() != 1)
+      throw std::runtime_error("SOFIE GPU: expected a single input");
+   Initialize();
+   std::vector<RGPUModel::Layer> layers;
+   std::string current = fInputTensorNames.front();
+   for (const auto &op : fOperators) {
+      const auto inputs = op->GetOpInputTensors();
+      const auto outputs = op->GetOpOutputTensors();
+      if (inputs.empty() || inputs.front() != current || outputs.size() != 1)
+         throw std::runtime_error("SOFIE GPU: unsupported graph connectivity");
+      op->LowerGPU(*this, layers);
+      current = outputs.front();
+   }
+   if (fOutputTensorNames.size() != 1 || fOutputTensorNames.front() != current)
+      throw std::runtime_error("SOFIE GPU: expected the sequential graph output");
+   return RGPUModel(precision, std::move(layers), fixedBatch);
+}
+
 void RModel::Initialize(int batchSize, bool verbose) {
    std::map<std::string, size_t> inputParams;
    if (batchSize > 0) {
diff --git a/tmva/sofie/test/CMakeLists.txt b/tmva/sofie/test/CMakeLists.txt
index 2826b6b823c25..e675ce0472245 100644
--- a/tmva/sofie/test/CMakeLists.txt
+++ b/tmva/sofie/test/CMakeLists.txt
@@ -169,3 +169,20 @@ if (ROOT_TORCH_FOUND AND ROOT_ONNX_FOUND AND NOT broken_onnx)
   endif()
 endif()
 
+
+ROOT_ADD_GTEST(TestSofieGPU TestSofieGPU.cxx LIBRARIES ROOTTMVASofieParser)
+target_include_directories(TestSofieGPU PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/../inc)
+set(SOFIE_GPU_TEST_BACKEND "" CACHE STRING "Run SOFIE device tests with CUDA or HIP")
+set_property(CACHE SOFIE_GPU_TEST_BACKEND PROPERTY STRINGS "" CUDA HIP)
+if(SOFIE_GPU_TEST_BACKEND STREQUAL "CUDA")
+  find_package(CUDAToolkit REQUIRED)
+  ROOT_ADD_GTEST(TestSofieGPUDevice TestSofieGPUDevice.cxx LIBRARIES ROOTTMVASofie CUDA::cudart)
+  target_compile_definitions(TestSofieGPUDevice PRIVATE SOFIE_TEST_COMPILER="${CUDAToolkit_NVCC_EXECUTABLE}")
+elseif(SOFIE_GPU_TEST_BACKEND STREQUAL "HIP")
+  find_package(hip CONFIG REQUIRED)
+  find_program(SOFIE_HIPCC hipcc HINTS "${hip_ROOT}/bin" REQUIRED)
+  ROOT_ADD_GTEST(TestSofieGPUDevice TestSofieGPUDevice.cxx LIBRARIES ROOTTMVASofie hip::host)
+  target_compile_definitions(TestSofieGPUDevice PRIVATE SOFIE_TEST_HIP SOFIE_TEST_COMPILER="${SOFIE_HIPCC}")
+elseif(SOFIE_GPU_TEST_BACKEND)
+  message(FATAL_ERROR "SOFIE_GPU_TEST_BACKEND must be empty, CUDA, or HIP")
+endif()
diff --git a/tmva/sofie/test/TestSofieGPU.cxx b/tmva/sofie/test/TestSofieGPU.cxx
new file mode 100644
index 0000000000000..14f66f61abf09
--- /dev/null
+++ b/tmva/sofie/test/TestSofieGPU.cxx
@@ -0,0 +1,112 @@
+#include "TMVA/ROperator.hxx"
+#include "TMVA/ROperator_Gemm.hxx"
+// For the licensing terms see $ROOTSYS/LICENSE.
+#include "TMVA/RGPUModel.hxx"
+#include "TMVA/RModelParser_ONNX.hxx"
+#include <gtest/gtest.h>
+#include <cstdint>
+#include <cstring>
+#include <limits>
+#include <sstream>
+
+using namespace TMVA::Experimental::SOFIE;
+namespace {
+std::string Varint(std::uint64_t v)
+{
+   std::string s;
+   do {
+      s += char((v & 127) | (v >= 128 ? 128 : 0));
+      v >>= 7;
+   } while (v);
+   return s;
+}
+std::string Number(int field, std::uint64_t n)
+{
+   return Varint(field * 8) + Varint(n);
+}
+std::string Bytes(int field, const std::string &s)
+{
+   return Varint(field * 8 + 2) + Varint(s.size()) + s;
+}
+std::string TensorInfo(const char *name, int type, int width)
+{
+   auto shape = Bytes(1, Bytes(2, "batch")) + Bytes(1, Number(1, width));
+   return Bytes(1, name) + Bytes(2, Bytes(1, Number(1, type) + Bytes(2, shape)));
+}
+std::string TinyModel(bool half, const std::string &domain = "", const std::string &op = "Gemm")
+{
+   const int type = half ? 10 : 1;
+   std::string raw = half ? std::string("\0\x3c\0\x40", 4) : std::string("\0\0\x80\x3f\0\0\0\x40", 8);
+   auto weight = Number(1, 2) + Number(1, 1) + Number(2, type) + Bytes(8, "w") + Bytes(9, raw);
+   auto node = Bytes(1, "x") + Bytes(1, "w") + Bytes(2, "y") + Bytes(4, op) + Bytes(7, domain);
+   auto graph =
+      Bytes(1, node) + Bytes(5, weight) + Bytes(11, TensorInfo("x", type, 2)) + Bytes(12, TensorInfo("y", type, 1));
+   return Number(1, 8) + Bytes(7, graph);
+}
+} // namespace
+TEST(SofieGPU, ParseFloatAndHalf)
+{
+   RModelParser_ONNX parser;
+   for (bool half : {false, true}) {
+      std::istringstream input(TinyModel(half));
+      auto model = parser.ParseGPU(input);
+      EXPECT_EQ(model.InputSize(), 2u);
+      EXPECT_EQ(model.OutputSize(), 1u);
+      EXPECT_EQ(model.GetPrecision(), half ? RGPUModel::Precision::Float16 : RGPUModel::Precision::Float32);
+      EXPECT_EQ(model.WorkspaceSize(3) % model.WorkspaceAlignment(), 0u);
+      const auto cuda = model.GenerateSource(RGPUModel::Backend::CUDA);
+      const auto hip = model.GenerateSource(RGPUModel::Backend::HIP);
+      EXPECT_NE(cuda.find("cuda_runtime.h"), std::string::npos);
+      EXPECT_NE(hip.find("hip/hip_runtime.h"), std::string::npos);
+      EXPECT_EQ(cuda.find("TMVA/"), std::string::npos);
+      EXPECT_EQ(hip.find("Malloc"), std::string::npos);
+      EXPECT_THROW(model.CreateSession(nullptr, 0, 3), std::logic_error);
+   }
+}
+TEST(SofieGPU, RejectInvalidGraphs)
+{
+   RModelParser_ONNX parser;
+   for (const auto &bytes : {TinyModel(false, "custom"), TinyModel(false, "", "Conv"), std::string("invalid")}) {
+      std::istringstream input(bytes);
+      EXPECT_THROW(parser.ParseGPU(input), std::runtime_error);
+   }
+}
+TEST(SofieGPU, CapacityAndShapes)
+{
+   RGPUModel::Layer layer{2, 1, {1.f, 2.f}, {0.f}};
+   RGPUModel model(RGPUModel::Precision::Float32, {layer});
+   EXPECT_THROW(model.WorkspaceSize(0), std::invalid_argument);
+   EXPECT_THROW(model.WorkspaceSize(std::numeric_limits<std::size_t>::max()), std::overflow_error);
+   RGPUModel fixed(RGPUModel::Precision::Float32, {layer}, 4);
+   EXPECT_THROW(fixed.WorkspaceSize(3), std::invalid_argument);
+   EXPECT_GT(fixed.WorkspaceSize(4), 0u);
+   layer.weights.pop_back();
+   EXPECT_THROW((RGPUModel(RGPUModel::Precision::Float32, {layer})), std::invalid_argument);
+}
+
+TEST(SofieGPU, UsesSharedOperatorRegistry)
+{
+   RModelParser_ONNX parser;
+   bool called = false;
+   parser.RegisterOperator("Gemm", [&](RModelParser_ONNX &, const onnx::NodeProto &) -> std::unique_ptr<ROperator> {
+      called = true;
+      throw std::runtime_error("custom registry entry");
+   });
+   std::istringstream input(TinyModel(false));
+   EXPECT_THROW(parser.ParseGPU(input), std::runtime_error);
+   EXPECT_TRUE(called);
+}
+
+TEST(SofieGPU, LowerExistingRModel)
+{
+   RModel model;
+   model.AddInputTensorInfo("x", ETensorType::FLOAT, std::vector<size_t>{1, 2});
+   model.AddInputTensorName("x");
+   float weights[] = {1.f, 2.f};
+   model.AddInitializedTensor("w", ETensorType::FLOAT, {2, 1}, weights);
+   model.AddOperator(std::make_unique<ROperator_Gemm<float>>(1.f, 1.f, 0, 0, "x", "w", "y"));
+   model.AddOutputTensorNameList({"y"});
+   auto gpu = model.MakeGPUModel(RGPUModel::Precision::Float32);
+   EXPECT_EQ(gpu.InputSize(), 2u);
+   EXPECT_EQ(gpu.OutputSize(), 1u);
+}
diff --git a/tmva/sofie/test/TestSofieGPUDevice.cxx b/tmva/sofie/test/TestSofieGPUDevice.cxx
new file mode 100644
index 0000000000000..371be074f70e1
--- /dev/null
+++ b/tmva/sofie/test/TestSofieGPUDevice.cxx
@@ -0,0 +1,119 @@
+// For the licensing terms see $ROOTSYS/LICENSE.
+#include "TMVA/RGPUModel.hxx"
+#include <gtest/gtest.h>
+#include <algorithm>
+#include <cmath>
+#include <stdexcept>
+#include <vector>
+#ifdef SOFIE_TEST_HIP
+#include <hip/hip_runtime.h>
+#include <hip/hip_fp16.h>
+#define GPU(name) hip##name
+#else
+#include <cuda_runtime.h>
+#include <cuda_fp16.h>
+#define GPU(name) cuda##name
+#endif
+
+using namespace TMVA::Experimental::SOFIE;
+namespace {
+struct Stream {
+   GPU(Stream_t) value {};
+   Stream()
+   {
+      if (GPU(StreamCreateWithFlags)(&value, GPU(StreamNonBlocking)))
+         throw std::runtime_error("stream creation failed");
+   }
+   ~Stream()
+   {
+      GPU(StreamSynchronize)(value);
+      GPU(StreamDestroy)(value);
+   }
+};
+struct Buffer {
+   void *data = nullptr;
+   explicit Buffer(std::size_t bytes)
+   {
+      if (GPU(Malloc)(&data, bytes))
+         throw std::bad_alloc();
+   }
+   ~Buffer() { GPU(Free)(data); }
+};
+float Value(float x)
+{
+   return x;
+}
+float Value(__half x)
+{
+   return __half2float(x);
+}
+} // namespace
+
+TEST(SofieGPUDevice, ExternalStreamAndWorkspace)
+{
+   int count = 0;
+   if (GPU(GetDeviceCount)(&count) != GPU(Success) || !count)
+      GTEST_SKIP() << "No GPU available";
+   int device = -1;
+   ASSERT_EQ(GPU(GetDevice)(&device), GPU(Success));
+#ifdef SOFIE_TEST_HIP
+   hipDeviceProp_t properties;
+#else
+   cudaDeviceProp properties;
+#endif
+   ASSERT_EQ(GPU(GetDeviceProperties)(&properties, device), GPU(Success));
+#ifdef SOFIE_TEST_HIP
+   const RGPUModel::CompileOptions options{RGPUModel::Backend::HIP, SOFIE_TEST_COMPILER, properties.gcnArchName};
+#else
+   const RGPUModel::CompileOptions options{RGPUModel::Backend::CUDA, SOFIE_TEST_COMPILER,
+                                           "sm_" + std::to_string(properties.major) + std::to_string(properties.minor)};
+#endif
+   for (bool half : {false, true}) {
+      RGPUModel::Layer first{2, 2, {1.f, -2.f, .5f, 1.f}, {.25f, -.5f}, 1.f, 1.f, true};
+      RGPUModel::Layer last{2, 1, {2.f, -1.f}, {.125f}};
+      std::vector<RGPUModel::Layer> layers{first};
+      for (int i = 0; i < 4; ++i)
+         layers.push_back({2, 2, {1.f, 0.f, 0.f, 1.f}, {0.f, 0.f}, 1.f, 1.f, true});
+      layers.push_back(last);
+      RGPUModel model(half ? RGPUModel::Precision::Float16 : RGPUModel::Precision::Float32, std::move(layers));
+      model.Compile(options);
+      EXPECT_THROW(model.Compile(options), std::logic_error);
+      Stream stream;
+      Buffer arena(model.WorkspaceSize(7)), input(14 * sizeof(float)), output(7 * sizeof(float));
+      auto session = model.CreateSession(stream.value, device, 7);
+      EXPECT_THROW(session->SetWorkspace(arena.data, 1), std::invalid_argument);
+      EXPECT_THROW(session->SetWorkspace(static_cast<char *>(arena.data) + 1, model.WorkspaceSize(7)),
+                   std::invalid_argument);
+      session->SetWorkspace(arena.data, model.WorkspaceSize(7));
+      EXPECT_THROW(session->Infer(input.data, output.data, 8), std::invalid_argument);
+      EXPECT_THROW(session->Infer(arena.data, output.data, 1), std::invalid_argument);
+      EXPECT_THROW(session->Infer(input.data, input.data, 1), std::invalid_argument);
+      auto run = [&](auto scalar) {
+         using T = decltype(scalar);
+         std::vector<T> host(14), result(7);
+         for (std::size_t i = 0; i < host.size(); ++i)
+            host[i] = T(float(i) * .125f - .625f);
+         for (std::size_t batch : {1u, 3u, 7u}) {
+            ASSERT_EQ(
+               GPU(MemcpyAsync)(input.data, host.data(), batch * 2 * sizeof(T), GPU(MemcpyHostToDevice), stream.value),
+               GPU(Success));
+            session->Infer(input.data, output.data, batch);
+            ASSERT_EQ(
+               GPU(MemcpyAsync)(result.data(), output.data, batch * sizeof(T), GPU(MemcpyDeviceToHost), stream.value),
+               GPU(Success));
+            ASSERT_EQ(GPU(StreamSynchronize)(stream.value), GPU(Success));
+            for (std::size_t row = 0; row < batch; ++row) {
+               const float x = Value(host[2 * row]), y = Value(host[2 * row + 1]);
+               const float a = std::max(0.f, Value(T(x - 2.f * y + .25f)));
+               const float b = std::max(0.f, Value(T(.5f * x + y - .5f)));
+               EXPECT_NEAR(Value(result[row]), Value(T(2.f * a - b + .125f)), half ? .002f : 1.e-6f);
+            }
+         }
+      };
+      if (half)
+         run(__half{});
+      else
+         run(float{});
+      ASSERT_EQ(GPU(StreamSynchronize)(stream.value), GPU(Success));
+   }
+}
diff --git a/tmva/sofie_parsers/CMakeLists.txt b/tmva/sofie_parsers/CMakeLists.txt
index 12c0ccba85612..8bd392fc2cacf 100644
--- a/tmva/sofie_parsers/CMakeLists.txt
+++ b/tmva/sofie_parsers/CMakeLists.txt
@@ -27,6 +27,7 @@ ROOT_STANDARD_LIBRARY_PACKAGE(ROOTTMVASofieParser
     src/ParseConv.cxx
     src/ParseConvTranspose.cxx
     src/ParseGemm.cxx
+    src/ParseGPU.cxx
     src/ParseGRU.cxx
     src/ParseIdentity.cxx
     src/ParseLeakyRelu.cxx
diff --git a/tmva/sofie_parsers/inc/TMVA/RModelParser_ONNX.hxx b/tmva/sofie_parsers/inc/TMVA/RModelParser_ONNX.hxx
index 11c128c80511b..3e7bfc2c103ed 100644
--- a/tmva/sofie_parsers/inc/TMVA/RModelParser_ONNX.hxx
+++ b/tmva/sofie_parsers/inc/TMVA/RModelParser_ONNX.hxx
@@ -2,6 +2,7 @@
 #define TMVA_SOFIE_RMODELPARSER_ONNX
 
 #include "TMVA/RModel.hxx"
+#include "TMVA/RGPUModel.hxx"
 
 #include <memory>
 #include <functional>
@@ -57,6 +58,10 @@ private:
 
 
 public:
+   /// Parse sequential Gemm/Relu models for CUDA/HIP code generation.
+   RGPUModel ParseGPU(const std::string &filename);
+   RGPUModel ParseGPU(std::istream &input);
+
    // Register an ONNX operator
    void RegisterOperator(const std::string &name, ParserFuncSignature func);
 
diff --git a/tmva/sofie_parsers/src/ParseGPU.cxx b/tmva/sofie_parsers/src/ParseGPU.cxx
new file mode 100644
index 0000000000000..74395b89ae5eb
--- /dev/null
+++ b/tmva/sofie_parsers/src/ParseGPU.cxx
@@ -0,0 +1,148 @@
+// Author: Christian Sonnabend, 2026
+// For the licensing terms see $ROOTSYS/LICENSE.
+#include "TMVA/RModelParser_ONNX.hxx"
+#include "TMVA/RGPUModel.hxx"
+#include "TMVA/ROperator.hxx"
+#include "onnx.hxx"
+
+#include <cmath>
+#include <cstdint>
+#include <cstring>
+#include <limits>
+#include <stdexcept>
+#include <unordered_map>
+#include <unordered_set>
+
+namespace TMVA::Experimental::SOFIE {
+namespace {
+float HalfToFloat(std::uint16_t bits)
+{
+   const int exponent = (bits >> 10) & 31, fraction = bits & 1023;
+   float value = exponent == 31
+                    ? (fraction ? std::numeric_limits<float>::quiet_NaN() : std::numeric_limits<float>::infinity())
+                    : std::ldexp(float(exponent ? fraction + 1024 : fraction), exponent ? exponent - 25 : -24);
+   return bits & 0x8000 ? -value : value;
+}
+std::vector<float> ReadWeights(const onnx::TensorProto &tensor, int type)
+{
+   if (tensor.data_type() != type || tensor.data_location() != onnx::TensorProto::DEFAULT ||
+       !tensor.external_data().empty())
+      throw std::runtime_error("SOFIE GPU: weights must be embedded and match the input precision");
+   std::size_t count = 1;
+   for (int i = 0; i < tensor.dims_size(); ++i) {
+      const auto dim = tensor.dims(i);
+      if (dim <= 0 || std::size_t(dim) > std::numeric_limits<std::size_t>::max() / count)
+         throw std::runtime_error("SOFIE GPU: invalid initializer shape");
+      count *= dim;
+   }
+   const auto &raw = tensor.raw_data();
+   const std::size_t width = type == onnx::TensorProto::FLOAT ? 4 : 2;
+   if (raw.empty()) {
+      if (type == onnx::TensorProto::FLOAT && tensor.float_data().size() == count)
+         return tensor.float_data();
+      if (type == onnx::TensorProto::FLOAT16 && tensor.int32_data().size() == count) {
+         std::vector<float> values;
+         for (auto bits : tensor.int32_data())
+            values.push_back(HalfToFloat(std::uint16_t(bits)));
+         return values;
+      }
+   }
+   if (raw.size() % width || raw.size() / width != count)
+      throw std::runtime_error("SOFIE GPU: invalid initializer data length");
+   std::vector<float> values(count);
+   for (std::size_t i = 0; i < count; ++i) {
+      std::uint32_t bits = 0;
+      for (std::size_t b = 0; b < width; ++b)
+         bits |= std::uint32_t(static_cast<unsigned char>(raw[i * width + b])) << (8 * b);
+      if (width == 2)
+         values[i] = HalfToFloat(std::uint16_t(bits));
+      else
+         std::memcpy(&values[i], &bits, 4);
+   }
+   return values;
+}
+} // namespace
+RGPUModel RModelParser_ONNX::ParseGPU(const std::string &filename)
+{
+   std::ifstream file(filename, std::ios::binary);
+   if (!file)
+      throw std::runtime_error("SOFIE GPU: cannot open " + filename);
+   return ParseGPU(file);
+}
+RGPUModel RModelParser_ONNX::ParseGPU(std::istream &input)
+{
+   auto proto = LoadModel(input);
+   if (!proto)
+      throw std::runtime_error("SOFIE GPU: malformed ONNX input");
+   const auto &graph = proto->graph();
+   if (graph.input_size() != 1 || graph.output_size() != 1 || graph.node_size() == 0)
+      throw std::runtime_error("SOFIE GPU: expected one input, one output and a nonempty sequential graph");
+   const auto &in = graph.input(0).type().tensor_type(), &out = graph.output(0).type().tensor_type();
+   const int type = in.elem_type();
+   if ((type != onnx::TensorProto::FLOAT && type != onnx::TensorProto::FLOAT16) || out.elem_type() != type ||
+       in.shape().dim_size() != 2 || out.shape().dim_size() != 2 || in.shape().dim(1).dim_value() <= 0)
+      throw std::runtime_error("SOFIE GPU: expected matching FP32/FP16 rank-two tensors");
+   const auto &batch = in.shape().dim(0), &outputBatch = out.shape().dim(0);
+   if (batch.value_case() != outputBatch.value_case() || batch.dim_value() != outputBatch.dim_value() ||
+       batch.dim_param() != outputBatch.dim_param() || (batch.dim_param().empty() && batch.dim_value() <= 0))
+      throw std::runtime_error("SOFIE GPU: invalid or inconsistent batch dimension");
+   fTensorTypeMap.clear();
+   fFusedOperators.clear();
+   RModel model;
+   // GPU tensors retain their requested precision; operator shape inference uses promoted FP32 values.
+   model.AddInputTensorInfo(graph.input(0).name(), ETensorType::FLOAT,
+                            std::vector<size_t>{1, size_t(in.shape().dim(1).dim_value())});
+   model.AddInputTensorName(graph.input(0).name());
+   RegisterTensorType(graph.input(0).name(), ETensorType::FLOAT);
+   std::unordered_set<std::string> names{graph.input(0).name()};
+   for (int i = 0; i < graph.initializer_size(); ++i) {
+      const auto &tensor = graph.initializer(i);
+      if (!names.insert(tensor.name()).second)
+         throw std::runtime_error("SOFIE GPU: duplicate initializer name");
+      auto values = ReadWeights(tensor, type);
+      std::vector<size_t> shape;
+      for (int d = 0; d < tensor.dims_size(); ++d)
+         shape.push_back(size_t(tensor.dims(d)));
+      model.AddInitializedTensor(tensor.name(), ETensorType::FLOAT, shape, values.data());
+      RegisterTensorType(tensor.name(), ETensorType::FLOAT);
+   }
+   std::string current = graph.input(0).name();
+   std::vector<size_t> nodes;
+   for (int i = 0; i < graph.node_size(); ++i)
+      nodes.push_back(size_t(i));
+   for (size_t i = 0; i < nodes.size(); ++i) {
+      const auto &node = graph.node(int(i));
+      if ((!node.domain().empty() && node.domain() != "ai.onnx") || node.output_size() != 1 || node.input_size() < 1 ||
+          node.input(0) != current || !names.insert(node.output(0)).second)
+         throw std::runtime_error("SOFIE GPU: unsupported graph connectivity/domain at " + node.name());
+      if (node.op_type() == "Gemm") {
+         if (node.input_size() < 2 || node.input_size() > 3)
+            throw std::runtime_error("SOFIE GPU: invalid Gemm input count");
+         std::unordered_set<std::string> attrs;
+         for (int a = 0; a < node.attribute_size(); ++a) {
+            const auto &attr = node.attribute(a);
+            const bool scalar = (attr.name() == "alpha" || attr.name() == "beta") &&
+                                attr.type() == onnx::AttributeProto::FLOAT;
+            const bool transpose = (attr.name() == "transA" || attr.name() == "transB") &&
+                                   attr.type() == onnx::AttributeProto::INT;
+            if (!attrs.insert(attr.name()).second || (!scalar && !transpose))
+               throw std::runtime_error("SOFIE GPU: unsupported or duplicate Gemm attribute " + attr.name());
+         }
+      } else if (node.op_type() == "Relu" && (node.input_size() != 1 || node.attribute_size())) {
+         throw std::runtime_error("SOFIE GPU: invalid Relu inputs or attributes");
+      }
+      auto op = ParseOperator(i, graph, nodes, {});
+      if (!op)
+         throw std::runtime_error("SOFIE GPU: unexpected operator fusion");
+      model.AddOperator(std::move(op));
+      current = node.output(0);
+   }
+   model.AddOutputTensorNameList({graph.output(0).name()});
+   auto gpu = model.MakeGPUModel(type == onnx::TensorProto::FLOAT ? RGPUModel::Precision::Float32
+                                                                  : RGPUModel::Precision::Float16,
+                                 batch.dim_param().empty() ? batch.dim_value() : 0);
+   if (current != graph.output(0).name() || out.shape().dim(1).dim_value() != std::int64_t(gpu.OutputSize()))
+      throw std::runtime_error("SOFIE GPU: inconsistent output shape");
+   return gpu;
+}
+} // namespace TMVA::Experimental::SOFIE
diff --git a/tmva/sofie_parsers/src/onnx.hxx b/tmva/sofie_parsers/src/onnx.hxx
index ba7bc079ced38..3dcb4cf086a38 100644
--- a/tmva/sofie_parsers/src/onnx.hxx
+++ b/tmva/sofie_parsers/src/onnx.hxx
@@ -504,6 +504,7 @@ private:
 class NodeProto {
 public:
    const std::string &op_type() const { return fOpType; }
+   const std::string &domain() const { return fDomain; }
    const std::string &name() const { return fName; }
    int input_size() const { return int(fInput.size()); }
    const std::string &input(int i) const { return fInput[i]; }
@@ -523,6 +524,7 @@ public:
          case 2: fOutput.push_back(detail::Str(r.ReadLen())); break;
          case 3: fName = detail::Str(r.ReadLen()); break;
          case 4: fOpType = detail::Str(r.ReadLen()); break;
+         case 7: fDomain = detail::Str(r.ReadLen()); break;
          case 5: {
             auto s = r.ReadLen();
             fAttribute.emplace_back();
@@ -536,7 +538,7 @@ public:
 
 private:
    std::vector<std::string> fInput, fOutput;
-   std::string fName, fOpType;
+   std::string fName, fOpType, fDomain;
    std::vector<AttributeProto> fAttribute;
 };