GPU-Initiated Discrete Simulated Bifurcation: Low-Latency Requests and Streaming Dense Couplings

GPU-Initiated Discrete Simulated Bifurcation: Low-Latency Requests and Streaming Dense Couplings YAOCHENG CHEN1 1

Instituto de Astronomia, Geofísica e Ciências Atmosféricas (IAG-USP), Universidade de São Paulo, Rua do Matão 1226, São Paulo, SP 05508-090, Brazil (e-mail: yaochengchen@usp.br)

Corresponding author: Yaocheng Chen (e-mail: yaochengchen@usp.br).

arXiv:2610.06622v1 [cs.DC] 5 Oct 2026

Yaocheng Chen is supported by a FAPESP Postdoctoral Fellowship (Grant No. 25/25586-5).

ABSTRACT GPU-based optimization faces two communication bottlenecks: coordinating frequent requests and delivering dense models that exceed device memory. We present a discrete simulated bifurcation (dSB) architecture that addresses both through NVIDIA DOCA GPUNetIO. For resident models, a persistent service receives field updates, executes each solve within one GPU thread block, and returns the result. Exact integer coupling sums, GPU work queues, and batched transmission keep the receive–solve–reply path on the device without a dedicated CPU data-path core. In comparisons with socket-based servers using the same solver, the largest latency gains occur under concurrent load. As the offered load increases from 400 to 800 thousand requests per second, median round-trip latency rises by only 6%. At the highest tested load, median and 99th-percentile latencies are 189 and 218 µs, compared with 288 and 609 µs for the tuned persistent CPU proxy across repeated runs. For models larger than device memory, a streaming solver retains dynamical state on the GPU and reuses incoming coupling tiles across replicas. It evaluates ten-million-variable dense binary matrices at approximately 307 Gb/s, consuming a 12.5-TB logical matrix through a 64-MiB packet buffer. Ground-state recovery on planted instances and agreement with reference executions verify the computation. Together, the two modes scale dSB to concurrent requests and dense models beyond GPU memory. INDEX TERMS DOCA GPUNetIO, GPUDirect RDMA, GPU-initiated networking, Ising machine, low latency, persistent kernel, portfolio optimization, simulated bifurcation.

I. INTRODUCTION

Many combinatorial decisions admit Ising or quadratic unconstrained binary optimization (QUBO) formulations [1], [2]. Simulated bifurcation (SB) solves these formulations through the dynamics of coupled nonlinear oscillators [3], [4]. SB is also used for repeated decisions driven by incoming observations. Examples include market-graph search for pairs trading, detection-to-track assignment in automotive perception, time– frequency resource allocation in wireless networks, and online beam selection [5]–[8]. Reported systems include a 33-µs market-feed-to-order trading loop, vehicle-mountable tracking above 20 frames/s, and wireless resource allocation below 0.5 ms. The trading, tracking, and online beam-control systems use ballistic SB (bSB), while the wireless resource-allocation study reports an SB machine. These applications motivate efficient repeated execution; here we study the communication path for discrete SB (dSB). In a dynamic decision process, instance ℓ has energy Eℓ (s) = − 12 s⊤ Jℓ s − h⊤ ℓ s. Manuscript for review

(1)

The instance index ℓ denotes a new decision, separately from the integration steps within one solve. Some applications reuse Jℓ and update only the field hℓ ; others update the quadratic coefficients as observations change. In online blackbox optimization, feedback trains a quadratic surrogate, and an embedded Ising machine selects the next action from the updated model [8]. The representation of Jℓ determines the cost of constructing, storing, and delivering each instance, so model updates and communication consume part of the decision budget. Graphics processing units (GPUs) offer a programmable environment for combining model construction, numerical optimization, and application processing. A conventional GPU service receives a request on the CPU, transfers its inputs, dispatches GPU work, and returns the result through the host network stack. For a short solve, this communication and dispatch can occupy a substantial fraction of the round trip. GPU-initiated networking provides a direct device-side path: GPUnet, GPUrdma, GPU-Ether, and related systems 1

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

established the communication mechanisms [9]–[12], and DOCA GPUNetIO exposes receive and transmit operations to CUDA kernels [13]. An optimization service must additionally coordinate model reuse, persistent scheduling, and solver execution without host dispatch between integration steps. Higher network bandwidth also increases the demand for coordination, while gains in general-purpose single-thread performance have slowed [14]. Moving packet handling and solver dispatch onto the device reduces the host coordination required as arrival rates and queue counts grow. NVIDIA reports this coordination cost even for CPU-driven DPDK pipelines that already place packets directly in GPU memory [15]. These workloads expose two complementary scaling limits. At small problem sizes, scaling is primarily a request-rate problem: the model fits in GPU memory, but repeated host coordination can limit the number of concurrent decisions. At large problem sizes, scaling becomes a coefficient-delivery problem: the dense model requires O(n2 ) storage, while the O(Bn) state of B solver replicas can remain resident. The first regime calls for efficient receive–solve–reply coordination; the second calls for sustained delivery of coefficients to the solver. Our approach is to bring incoming coefficients directly to GPU-resident dSB state and perform frequent scheduling on the device. A resident-model service accepts small field updates and executes each complete trajectory within one persistent worker block. We compare it with socket-based launchper-request and CPU-proxied persistent servers using the same solver, protocol, and client, measuring latency under load and CPU data-path requirements. A streaming solver consumes nonresident coupling tiles while retaining positions, momenta, spins, and accumulators on the GPU. Each coupling evaluation requires a matrix pass. The request service schedules complete solves on the GPU; the streaming implementation schedules tiles on the GPU with host coordination at pass boundaries. The contributions are: • A dSB solver that assigns each request to one block, retains

dynamical state in registers and models in shared memory, and computes exact coupling sums with XOR/popcount or DP4A (Section IV). • A resident-model protocol and GPU-resident receive–solve– reply pipeline with bounded work queues and batched transmission (Sections III and V). • An evaluation against two socket-based CPU-mediated implementations across 19 request configurations, examining dedicated CPU data-path cores, host-state dependence, and latency under load (Sections VI and VII). • A streaming dSB implementation that evaluates dense n = 107 binary coupling matrices at 306–307 Gb/s for one or eight replicas, with planted-instance and referenceexecution validation (Section VIII). II. BACKGROUND AND RELATED WORK Manuscript for review

A. DISCRETE SIMULATED BIFURCATION

SB minimizes the Ising energy X X E(s) = − 12 Jij si sj − hi si , i̸=j

s ∈ {−1, +1}n , (2)

i

by integrating a classical system of positions xi and momenta yi through a bifurcation [3], [4]. In the discrete variant the coupling force depends only on sign(x). With detuning δ = 1, pump schedule pk = k/(K − 1) for K > 1 (and p0 = 0 for K = 1), coupling gain ξ and time step ∆t, one step is P ai = j Jij sign(xj ) + hi , (3)
yi ← yi + −(δ − pk )xi + ξai ∆t, xi ← xi + ∆t δ yi .

(4)

followed by the inelastic wall: if |xi | > 1 then xi ← sign(xi ) and yi ← 0. The answer is sign(x) after K steps; B independent replicas are integrated from different initial states and the lowest-energy one is returned. We use sign(0) = −1 and √ ξ = 0.5 n − 1/∥J ∥F by default, as in [16]. B. GPU AND FPGA SB IMPLEMENTATIONS

Field-programmable gate array (FPGA) implementations of SB exploit specialized arithmetic pipelines and tightly integrated I/O to achieve low latency, with multi-chip designs extending the coupling computation [4], [17], [18]. Trading systems demonstrate the value of this specialization in a complete decision loop [19]–[21]. GPU implementations offer a complementary advantage: solver variants, coefficient formats, and surrounding numerical or learning stages can be changed within a common software environment. The companion GPU implementation [16] provides whole-device and bit-packed execution paths. Here, an entire dSB request is confined to one block, allowing independent requests to execute concurrently while persistent device-side agents handle communication. C. GPU-CENTRIC NETWORKING

GPUnet and GPUrdma support communication initiated by GPU threads [9], [10], and GPUDirect Async investigates device-driven communication on InfiniBand [22]. GPU-Ether and Janus provide GPU-native packet I/O and virtualized packet-processing support on Ethernet [11], [12]. DOCA GPUNetIO supplies device APIs for receiving into GPU memory and posting sends [13], [23]; persistent kernels keep these device-side agents available between requests [24]. NCCL’s GPU-initiated networking also includes a GPUNetIObacked transport [25], and a recent survey reviews GPU-centric communication methods [26]. Kernel bypass and GPU initiation address different parts of the communication path. DPDK bypasses the operatingsystem network stack [27]; GPUdev can combine it with GPUDirect RDMA and persistent GPU processing while the CPU manages network queues and notifies the GPU [28]. GPUNetIO additionally permits CUDA kernels to control receive and transmit operations directly [15]. Placing packets directly in GPU memory therefore does not by itself take the 2

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

TABLE 1. Message types. All fields after the Ethernet header are little-endian and naturally aligned.

IV. REQUEST-PER-BLOCK SOLVER A. MAPPING

Message

Payload

LOAD

one fragment ( ≤1408 B) of a packed coupling model; header: model id, n, format, offset, fragment index and count, total size, generation, default ξ none; sent when all fragments of a generation have arrived optional field h (n signed 16-bit integers); header: model id, n, K , B, seed, ξ, ∆t best spin configuration, ⌈n/32⌉ words; header: status, 2E, best replica, worker id, server-side timestamps (queue, solve, transmit)

LOAD_ACK SOLVE RESULT

CPU out of orchestration. Host-driven designs also support high-throughput bulk ingest: Romein reaches high Ethernet rates for a GPU workload with a CPU-driven receive path [29]. Our work integrates these communication mechanisms with coefficient delivery, resident solver state, and persistent scheduling for dSB. The evaluation covers request-rate scaling against socket-based servers using the same solver and modelsize scaling through streamed coupling evaluation. III. SERVICE PROTOCOL

Every message is one raw Ethernet frame with a dedicated EtherType, a 64-byte header including the 14-byte Ethernet header, and a payload (Table 1). Raw frames were chosen over UDP so that the network interface controller (NIC) can select traffic by EtherType alone and the GPU does not have to parse IP headers; the design carries over to UDP by matching on a port instead. A SOLVE frame for n = 128 is 320 bytes and the largest frame in the request-service tests, a model fragment, is 1472 bytes, so the standard MTU suffices. A model is uploaded once and remains resident in a GPU model store of 64 slots of 512 KiB until it is replaced by a new generation. Three packed formats are supported (with W = ⌈n/32⌉): PM1, one bit plane of nW words for dense ±1 couplings; TERNARY, a plane for Jij = +1, a plane for Jij = −1 and a row constant |Pi | − |Mi |; and INT8, dense int8 rows. The field travels with each request. This matches problems in which the quadratic part changes slowly and the linear part changes quickly. For mean–variance portfolio selection with a cardinality penalty [30], min q x ⊤ Σx − µ⊤ x + λ(1⊤ x − k)2 ,

x∈{0,1}n

(5)

the substitution x = (1+
Jij = −Qij /2 for i ̸=⊤j and P s)/2 gives Jii = 0 and hi = − 12 Q + c i with Q = qΣ + λ11 and j ij c = −µ−2λk1. With fixed covariance and penalty parameters, J is reused, and each update of the expected returns µ gives a new h. The client quantizes J to int8 and h to int16 with a common scale. All solver and reference evaluations use the resulting quantized Ising objective, with cardinality enforced through the penalty term. Manuscript for review

One request is solved by one thread block of T = 256 threads. With np = ⌈n/32⌉ · 32, the Bnp pairs (replica b, row r) are distributed cyclically, pair p = bnp + r to thread p mod T , and each thread keeps the x and y of its at most P = 16 pairs in registers for the whole trajectory; hence Bnp ≤ TP = 4096. Because np and T are multiples of 32, the 32 lanes of a warp always hold 32 consecutive rows of the same replica. After the update the new signs of a warp are therefore exactly one 32-bit word of that replica’s sign vector: __ballot_sync forms the word and lane 0 stores it, without atomics or shuffles. Sign vectors are double buffered in shared memory, so a step needs a single __syncthreads(): all warps read buffer k mod 2 and write buffer (k + 1) mod 2. Shared memory holds the coupling rows, padded to a stride Wp = 4 · odd ≥ W words so that the 16-byte loads of eight consecutive lanes fall in distinct banks, the two sign buffers of B rows each, the field and the per-replica energies. The footprint is S = 4np Wp c + 8BWp + 4np + 8B + S0 ,

(6)

with c = 1 for PM1, c = 2 (plus 4np for the row constants) for TERNARY, and S0 = 1024 B for the request header and the reply staging area; for INT8 the coupling term is np Ls with a byte stride Ls = 16 · odd ≥ np and the sign buffers are 2BLs bytes. The server fixes one S for all workers, which sets how many blocks are resident per streaming multiprocessor (SM) and so the number of requests in service at once; with S = 48 KiB two blocks fit on each of the 188 SMs of the server GPU, 376 in all. B. EXACT COUPLING SUMS

Proposition 1. Let J be integer with zero diagonal, s ∈ {±1}n and vj = 1[sj = +1] packed into words. (i) If Jij = 2uij − 1 for j ̸= iP (bit uij = 1 for Jij = +1) with uii = 0 and zero padding, then j Jij sj = n − 2 popc(ui ⊕ v) + si . (ii) If Jij ∈ {−1, 0, 1} with Pij = 1[Jij = 1] and Mij = 1[Jij = −1], then X
Jij sj = 2 popc(Pi ∧ v) − popc(Mi ∧ v) j

− |Pi | − |Mi | . (iii) If |Jij | ≤ 127, |hi | < 215 and n ≤ 2048, the int32 DP4A sum ai of (3) is exact and |ai | < 224 , so ai is exact in FP32. Consequently, if (4) is evaluated with the same IEEE roundto-nearest operations in the same order, the GPU and a CPU reference produce identical trajectories, spins and energies. Proof. (i) A column j < n, j ̸= i contributes +1 when uij = vj and −1 otherwise; the 32W − n padding columns have u = v = 0 and contribute +1 each; the diagonal, stored as uii = 0, contributes +1 if vi = 0 and −1 otherwise, that is −si . Summing over all 32W columns gives 32W − 2 popc(ui ⊕ v), and P removing the padding and diagonal terms gives the identity. (ii) j Pij sj = 2 popc(Pi ∧ v) − |Pi | and likewise for M . (iii) |ai | ≤ 127 · 2047 + 215 < 224 . The final statement 3

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

follows because ai is the same integer on both sides and every remaining operation is a correctly rounded IEEE operation on identical operands.

their replies sent. Up to St requests are in flight. Models are reassembled on the CPU and copied to the GPU model store. B. CPU-PROXIED PERSISTENT KERNEL

The GPU update uses explicitly rounded multiplication, addition, and subtraction (__fmul_rn, __fadd_rn, and __fsub_rn). These operations prevent fused multiply– add contraction. The host reference is compiled with -ffp-contract=off. Initial states are drawn from a counter-based hash of (seed, replica, row), so the reference can regenerate them. The self-test and the client use this property for the energy and reference-output checks specified in Section VI. The equivalence statement uses identical initial states, parameters, tie breaking, rounding, and treatment of subnormal values. C. COST OF ONE REQUEST

A PM1 step issues Bn · 4⌈W /4⌉ population counts per block (twice as many AND/POPC pairs for TERNARY), and an INT8 step Bn np /4 DP4A instructions. With an instruction throughput r per SM per cycle and clock frequency f , the solve time is approximately  ops  tsolve ≈ tload + K + tbar , (7) r f /β where ops is the instruction count per step and β is the number of blocks sharing the SM, tbar the cost of a block barrier and tload the one-time copy of the model rows from L2 into shared memory. Table 3 gives the measured values: the time per step is 0.5–0.6 µs for the smallest blocks, where the barrier and the instruction latency rather than the throughput of the SM set the pace, and grows to 4.7 µs for n = 128, B = 32, INT8. A 100-step request on 128 variables with up to eight replicas therefore costs 60–150 µs, so the 10–19 µs of host-side work that Section V removes is a noticeable fraction of a short request. V. SERVER ARCHITECTURES

Figure 1 shows the three servers compared in this paper. The two persistent designs share a worker loop and the ring organization shown in Fig. 1(c): a work ring of C = 4096 items, each pointing at a frame and carrying a sequence number and a receive timestamp, and a reply ring of 4096 slots of 512 B. A worker takes a ticket t with one atomic increment, waits for item t mod C to carry sequence number t + 1, copies the request header and field into shared memory, solves, allocates a reply slot, writes the reply and publishes a ready word that contains the slot’s sequence number and the reply length. What differs is who fills the work ring and who drains the reply ring. A. LAUNCH PER REQUEST

The conventional design uses a nonblocking AF_PACKET socket to receive each frame, copies it into pinned memory, and on one of St = 4 CUDA streams issues a host-to-device copy, a one-block launch of the solver, a device-to-host copy of the reply and an event; completed events are polled and Manuscript for review

The second design removes launches and explicit copies but keeps the CPU on the data path. The CPU writes each frame into a pool of host-mapped pinned buffers and publishes a work item in a host-mapped ring; 376 persistent worker blocks, running the same worker loop as the GPU-initiated server, read the request over PCIe and write the reply into a host-mapped reply ring, which the CPU polls and sends. The same mapped-buffer organization would also work with kernel-bypass networking; here it uses the kernel socket path. C. GPU-INITIATED SERVER

The third design keeps the request data path on the GPU after startup. The host creates a DOCA Flow root pipe that steers the service EtherType to one receive queue, a cyclic GPUNetIO receive queue whose 16,384-slot buffer (64 MiB) lives in GPU memory, and a send queue of depth 4096 whose completions are also handled on the GPU; the reply ring is allocated in GPU memory and registered with the NIC through a DMA-BUF. It then launches two kernels on two streams, an I/O kernel of two 32-thread blocks first, so that it is resident before the workers fill the GPU, and the worker kernel with Nw = 372 blocks second (four of the 376 resident slots are left free so that the device is never completely booked by persistent blocks). The kernels communicate only through the rings and need no cooperative launch. • The receive warp (block 0 of the I/O kernel) calls the warpscope receive of GPUNetIO in a loop, taking up to 64 frames per call with a timeout of 500 ns: the call returns when 64 frames have arrived or when the timeout has expired. The timeout cannot be zero, because on this platform a zero timeout blocks until a frame arrives and the warp would never see the shutdown flag on an idle link; our first runs used 5 µs. It reads the 64-byte header of each frame with L2coherent loads, validates it, and for a SOLVE frame copies the frame out of the NIC’s cyclic buffer into a frame buffer owned by the work slot (a 320-byte warp copy for n = 128), stamps the receive time from the GPU global timer and publishes the item. For a LOAD fragment the warp copies the payload into the model store and, when all fragments of a generation have arrived, publishes the model and queues an acknowledgement. • Worker blocks run the common worker loop; the worker that solves a request stamps its claim time and the end of the solve into the reply. • The transmit thread (thread 0 of block 1 of the I/O kernel) scans the reply ring in order and stamps the time each reply has waited since the end of its solve. It then writes one send descriptor per ready reply, up to 32, requests a completion only for the last, rings the doorbell once, and waits for that completion before it releases the slots. The NVIDIA gpunetio_simple_send sample offers this sequence and a simpler one in which every reply is posted with its own 4

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

(c) GPU-initiated: persistent I/O and worker kernels, no CPU on the path DMA RX cyclic buffer ConnectX-7 RX / TX queues (GPU memory)

RX warp (I/O kernel)

work ring tickets

worker blocks one request each

reply ring (NIC-registered)

LOAD TX thread (I/O kernel)

model store

send queue: doorbell, NIC reads reply by DMA

socket recv

pinned copy H2D

launch 1 block

socket recv

host-mapped work ring

persistent workers

D2H event

sendto

host-mapped reply ring, poll

(a) launch per request

sendto

(b) CPU-proxied persistent kernel

FIGURE 1. The three servers. All share the wire format, the client and the per-request solver; they differ only in how a frame reaches a worker block and how the reply leaves. Blue: CPU; green: GPU; orange: NIC.

send call and completion notification; our first version used the simpler one, and Section VII-F reports both, because the choice sets the server’s capacity. Three rules make the rings safe without locks. First, the NIC recycles its receive buffer at line rate regardless of how far the workers are behind, so the receiver copies each request before publishing it. Second, a worker acknowledges its ticket in the slot after its reply is published, and the receiver reuses a slot only when the previous ticket of that slot has been acknowledged; otherwise the request is dropped and counted, which bounds the backlog at C = 4096 requests. Third, reply slots are allocated after the solve rather than at arrival, so a long request does not delay the replies of shorter requests that finish after it. Producers write an item and then its sequence number behind a memory fence, and consumers read requests and models with L2-coherent loads (__ldcg), because the data were written by another SM. Model fragments are tracked by index in a per-generation bitmap, so a retransmitted fragment is ignored and a retransmission after completion re-sends the acknowledgement. Host monitoring and GPU occupancy affect latency. The host must be able to observe the server (counters, shutdown) without touching the data path: the request-service counters live in managed memory and are read once per second, and the shutdown flag lives in host-mapped memory. In a first version every spinning block read that flag on every iteration; 372 blocks issuing PCIe reads in a tight loop starved the NIC’s own DMA and doorbells and raised the round-trip time to milliseconds. Each block now polls the flag no more than once every 0.5 ms, reducing PCIe traffic. Four block slots are left free for driver work.

Manuscript for review

D. WHERE THE TIME GOES

The round-trip time seen by a client decomposes as T = twire + tin + tq + tsolve + tout .

(8)

For the launch path tin contains the NIC’s delivery to the operating-system kernel, the socket receive, a pinned copy, the host-to-device copy and the launch, and tout the deviceto-host copy, the completion poll and the socket send. For the proxy path they contain the operating-system receive path, the socket calls and the PCIe round trips of the mapped rings. For the GPU-initiated path they reduce to a direct memory access (DMA) write, the receive poll, a ring hand-off through L2 and a doorbell. Every server stamps its residence time into the reply. The GPU-initiated server stamps tq (receive to claim), tsolve and the transmit wait with the GPU global timer; the two baselines stamp the interval from the return of recv() to the call of sendto() with the CPU clock and report it as tsolve plus a host term. The remainder, T minus recorded server residence, includes wire transit and the client’s network processing. Its server-side contribution follows the timestamp boundaries: NIC-to-GPU DMA and the receive-call wait for GPUNetIO, or kernel processing before recv() returns and after sendto() is called for the baselines. Section VII-D uses this decomposition to interpret the end-to-end measurements. VI. EXPERIMENTAL METHODOLOGY A. PLATFORM

The server is an x86-64 workstation with eight NVIDIA RTX PRO 6000 Blackwell Workstation Edition GPUs, of which one is used, and a ConnectX-7 NIC; the client is a dual-socket GH200 system, whose GPUs are idle in the request-service campaign and one of which is the generator of the streaming mode (Section VIII); the two are connected back to back by 400 GbE (Table 2). The server GPU is attached with a 256 MiB 5

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

TABLE 2. Configuration of the measurement campaign.

Item

Value

Server GPU

RTX PRO 6000 Blackwell Workstation Edition (sm_120, 188 SMs), driver 595.58.03 (reports CUDA 13.2 support) ConnectX-7, firmware 28.48.1000, MLNX_OFED 26.04-0.8.6, DOCA 3.4.0112 Linux 5.14.0-570.62.1.el9 (RHEL 9.6); baseline network thread pinned to one core GH200 (2× 72 Neoverse-V2 cores), Linux 7.0.0-1019-nvidia-64k, ConnectX-7 firmware 28.43.8002; TX and RX threads on two dedicated cores, AF_PACKET 400 GbE direct attach, MTU 9000 (frames ≤1472 B in the request service, 4096 B frames with 4032 B of payload in the streaming mode) generator on GH200 GPU 1: 4 send queues (1–8 in one sweep), 96 synthesis blocks by default; receiver: one receive queue of 16,384 slots of 4 KiB (64 MiB), one receive warp, 372 worker blocks of 512 threads; B ∈ {1, 8} 256 threads, 48 KiB shared memory, 2 per SM: 372 (GPUNetIO), 376 (proxy) 4 CUDA streams default; tuned (no processor idle states, performance governor, no interrupt moderation, interrupts pinned) 100,000 requests after 1,000 warm-up; 291 runs 500 requests per single-launch point; 200,000 (GPU rings) and 100,000 (CPU rings) per pipeline point

Server NIC Server OS Client

Link

Streaming mode

Worker blocks Launch server Host states Network runs GPU-only runs

BAR1 window (resizable BAR disabled), which is why the receive ring is 64 MiB. The client and the baseline servers pin their network threads to dedicated cores. The CPU-mediated servers are measured in two host states. In the default state interrupt moderation, processor idle states and the interrupt affinity of the NICs are left at the distribution’s defaults. In the tuned state a script of the repository disables every processor idle state deeper than C0 on all cores, selects the performance governor, turns interrupt moderation, GRO and LRO off and pins the NIC’s interrupts next to the polling core, on the server and on the client alike. The GPU-initiated server is reported with the tuned client. In a preliminary comparison using perreply sends, its reference-request round trip was 168.5 µs in the default state and 168.6 µs in the tuned state, so we report it in one host state. The server uses PCIe 5.0, and its topology places the GPU in NUMA node 1. GPU-to-NIC paths are reported as SYS for mlx5_0–mlx5_3 and NODE for mlx5_4–mlx5_7: the former traverse the inter-socket interconnect, and the latter traverse host bridges within the same NUMA node. All reported paths therefore include host-side interconnects. NVIDIA’s GPUDirect RDMA guidance identifies a path confined to PCIe switches as the preferred topology [23]. GPU/NIC placement Manuscript for review

and PCIe switching are consequently important configuration choices for scaling the streaming receiver. B. WORKLOADS

Three instance families cover the three formats: dense ±1 Sherrington–Kirkpatrick couplings (PM1), sparse ternary couplings at 5% density (TERNARY), and the int8 portfolio model of Section III with a four-factor covariance, a cardinality target of n/10 and a pool of 256 field vectors generated from AR(1) expected returns and reused cyclically (INT8). The GPU-only sweep covers n ∈ {64, 128, 192, 256, 512}, B ∈ {1, 8, 32} and K ∈ {50, 100, 200, 400} where Eq. (6) admits them, 120 configurations in all. The network sweep uses a scaling grid of 19 configurations: INT8 with n ∈ {64, 128, 192} and K ∈ {50, 100, 200, 400} at B = 8; PM1 with n ∈ {64, 128, 256, 512} at B = 8, K = 100; TERNARY at n = 256; and INT8 n = 128 with B ∈ {1, 32}. The reference configuration for the load and distribution experiments is INT8, n = 128, B = 8, K = 100. Every request has its own seed. The supplied request sweep uses the client defaults ∆t = 1.0 and base seed 1, with model seed 1000 + m for model index m and request seed 7919 + i for request index i. The model’s default coupling gain is used. With these settings the archived sweep is deterministic and can be reproduced. C. MEASUREMENTS

GPU only. dsb_bench measures, without a NIC, (i) the inkernel solve time of one request from the GPU global timer, (ii) the launch-to-completion time from CUDA events and (iii) the host round trip H2D + launch + D2H from the CPU clock; it also saturates the persistent pipeline with a GPU generator warp and a GPU consumer warp that use the same rings as the GPU-initiated server, and the proxied pipeline with a CPU producer and consumer that, like the proxy server, drain all ready replies in each iteration and keep up to 4096 requests outstanding. Network. dsb_client measures the interval from the client timestamp immediately before request-frame construction and sendto to the timestamp after recv in closed loop with 1, 4, 16 and 64 outstanding requests and in open loop with Poisson arrivals at 10 to 800 kreq/s. The client sends from one thread; at 700 and 800 kreq/s offered with the tuned client its own send timestamps show that it emitted 596–681 kreq/s, which is therefore the highest load of the campaign. The first 1000 requests of every run are discarded. We report the 50th, 99th and 99.9th percentiles; the servers’ own timestamps split each round trip into server residence and the remainder of Section V. By construction, the launch and proxy servers each occupy one busy-polling CPU core on the data path and the GPU-initiated server none. Checks are sampled by request index: energy is recomputed from returned spins for up to the first 3000 requests of each run (including warm-up), and a CPU integration checks the final result for the first 20.

6

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

E. REPORTING POLICY

The archive contains 291 network runs covering 246 distinct settings. For each server, workload, arrival rate, and window, the analysis uses the latest run; the unloaded comparison uses 19 settings per server configuration. For GPUNetIO and the tuned proxy at 400–800 kreq/s offered, we report medians across repeated runs with the current code. Table 11 gives the run counts and ranges. Latency percentiles are computed within each run before aggregation. The configured Poisson rate is termed ‘‘offered load,’’ and the measured completion rate is termed ‘‘achieved rate.’’ At the two highest target rates, the client injects 596–681 kreq/s. Latency percentiles are computed over completed requests, and loss is reported separately. VII. RESULTS

The repository includes the benchmark scripts, analysis scripts, and per-request records used to produce these results. A. CORRECTNESS

The self-test compares the GPU with the CPU reference on eleven configurations across the three formats, through both the single-launch and persistent paths. All eleven configurations and 44 replies through the persistent ring passed before the measurement campaign. Across the 291 network runs, the client recomputed 850,874 reply energies and compared 5,820 final results with CPU-reference integrations, with no mismatch on any server.

H2D + launch + D2H − solve (μs)

D. BASELINE SCOPE

All three servers share the solver, workloads, protocol, and client. Both baselines use the kernel socket path, with tuning of processor idle states, interrupts, and offloads. The comparison measures the combined effects of the network stack and CPU– GPU coordination (Section IX-E).

pm1 50 er ov

a he

d=

s

ternary

e olv

int8

20 er ov

10

a he

d=

% 10

s of

e olv

5 100

1000

in-kernel solve time (μs) FIGURE 2. Launch-path overhead (H2D + launch + D2H minus the in-kernel solve) against the solve time, one point per configuration, no network. The overhead is 10–19 µs for every configuration.

sustains 1.9–2.0 Mreq/s, the rate at which one generator lane can publish items; 62 of the 120 configurations reach this ceiling. For longer solves the workers become the limit: the throughput reaches 73–96% of Nw /tsolve , the ideal rate for the GPU-only benchmark’s Nw = 375 workers, the shortfall being the slowdown of two blocks sharing an SM. The CPUproxied pipeline sustains 165–1025 kreq/s. For solves below 150 µs it spends 1.0–1.6 µs of CPU time per request, and the producer and consumer loops are the limit; for the reference configuration the ratio is 1983 to 658 kreq/s, 3.0×. For the longest solves both pipelines approach the worker bound. The GPU-only benchmark reserves one of the 376 resident block slots for its generator/consumer control work, leaving 375 solver workers; the GPUNetIO server reserves four slots for communication and uses 372. D. END-TO-END LATENCY

B. COST OF ONE REQUEST ON THE GPU

Table 3 and Fig. 2 separate the in-kernel solve time from the cost of the launch path without the network. Over the 120 configurations the solve takes 31 µs to 1.9 ms and grows linearly with K . Launching the block costs a further 3.7–7.2 µs (median 6.0) and the full host round trip of host-to-device (H2D) copy, launch, and device-to-host (D2H) copy costs 10.0–19.3 µs above the solve (median 13.0), across the tested requests. Relative to the solve, this overhead adds 60% to the smallest request (31 µs), 8% to the reference portfolio request (152 µs), and 1% to a 1-ms request. These measurements use an idle host with a polling thread and exclude network I/O; Section VII-D reports the end-to-end latency.

Table 5 reports unloaded round-trip time (RTT) for the five server configurations. For the reference request, the GPUNetIO median is 166.4 µs, compared with 175.6 µs for tuned proxy and 178.5 µs for tuned launch. The improvement over the faster tuned baseline is approximately 5%. The recorded solve times are nearly identical (149.9–152.9 µs), so the difference comes from the request path around the solver. Table 6 reports component percentiles; Fig. 3 uses component means for an additive decomposition. On GPUNetIO, the median work-ring and transmit waits are 2.0 and 0.8 µs. The baselines’ median host terms are 11.3–16.1 µs, which is roughly what the GPU-initiated path saves in server residence. 1) Host-state dependence

C. PIPELINE CAPACITY WITHOUT THE NETWORK

Table 4 compares the saturated throughput of the GPU-resident rings, fed by a GPU generator warp, with that of the hostmapped rings fed by a CPU thread. For every configuration whose solve is short enough, the GPU-resident pipeline Manuscript for review

Table 7 separates recorded residence from the residual RTT. The default-host CPU-mediated runs exhibit fast and slow residual modes, with the latter contributing roughly 100 µs. Tuning idle states, interrupt moderation, affinity, and offloads brings the residual close to 11 µs throughout the grid. The 7

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

TABLE 3. One request on the server GPU, no network (median over 500 requests). All times are in µs. The last two columns are the launch-to-completion time and the H2D + launch + D2H round trip, each minus the in-kernel solve, as the median over K (range in parentheses).

in-kernel solve, p50 (µs)

launch

H2D+launch+D2H

format

n

B

K =50

100

200

400

−solve

−solve (µs)

pm1 pm1 pm1 pm1 pm1 pm1 pm1 pm1 pm1 pm1 pm1 pm1 ternary ternary ternary ternary ternary ternary ternary ternary ternary ternary int8 int8 int8 int8 int8 int8 int8 int8

64 64 64 128 128 128 192 192 256 256 512 512 64 64 64 128 128 128 192 192 256 256 64 64 64 128 128 128 192 192

1 8 32 1 8 32 1 8 1 8 1 8 1 8 32 1 8 32 1 8 1 8 1 8 32 1 8 32 1 8

31.3 41.7 89.9 32.2 57.9 166.6 34.9 81.2 35.4 100.8 50.9 224.4 35.5 47.3 103.0 36.4 67.0 197.9 39.8 99.6 41.1 121.9 34.5 48.5 111.3 38.1 82.5 254.3 44.9 131.1

58.6 77.7 168.7 59.8 108.1 313.3 64.7 152.9 65.2 190.3 93.4 428.7 62.8 84.0 188.3 64.6 120.3 370.4 70.6 183.8 72.6 227.0 60.9 87.1 206.8 66.5 151.7 487.1 78.8 246.2

112.8 149.6 326.5 114.9 208.7 606.1 124.3 296.2 124.9 369.1 178.1 837.2 117.7 157.5 358.9 120.9 227.1 713.8 132.0 352.2 135.5 437.1 113.5 164.4 398.0 123.4 290.1 953.2 146.7 476.4

221.1 293.2 641.5 225.1 409.7 1192.5 243.4 582.6 244.4 726.7 347.8 1654.7 227.2 304.4 699.5 233.5 440.3 1408.1 254.1 689.2 261.5 857.4 218.8 318.9 780.3 237.2 566.9 1885.1 282.3 936.9

6.0 5.8 6.1 6.1 5.7 6.3 6.0 6.1 6.5 6.6 5.6 6.0 5.7 5.7 5.9 6.0 5.7 6.4 6.0 6.0 6.5 6.0 6.2 6.0 6.2 6.2 5.8 6.2 6.0 5.8

13.7 (11.4–18.8) 12.4 (12.2–15.8) 14.2 (11.8–17.6) 13.5 (12.5–17.9) 14.1 (12.2–16.0) 15.5 (11.8–19.3) 13.1 (11.1–14.7) 12.0 (11.3–17.2) 12.3 (11.4–14.5) 12.3 (11.7–17.7) 15.5 (12.7–18.6) 13.3 (12.2–16.5) 15.2 (11.9–18.6) 14.5 (14.3–15.2) 12.5 (11.7–13.3) 13.0 (11.3–14.4) 14.7 (11.9–17.8) 13.1 (11.1–17.4) 16.0 (11.4–17.1) 13.3 (11.4–16.5) 12.2 (10.0–14.3) 11.6 (10.0–16.7) 13.4 (11.5–15.6) 12.3 (11.3–13.0) 13.4 (11.7–18.5) 11.8 (11.2–12.2) 14.0 (12.0–17.1) 17.2 (11.4–19.2) 12.9 (11.5–16.6) 14.1 (11.9–16.5)

TABLE 4. Saturated persistent pipelines without the network, B = 8, default host state. The GPU-only benchmark uses Nw = 375 workers (the networked GPUNetIO server uses 372). Nw /tsolve is the bound at the single-request solve time. The CPU proxy drains all ready replies in each iteration and keeps up to 4096 requests outstanding.

GPU-resident rings (GPU producer/consumer)

Manuscript for review

host-mapped rings (CPU proxy)

format

n

B

K

kreq/s

Nw /tsolve (kreq/s)

p99 (µs)

kreq/s

p99 (µs)

pm1 pm1 pm1 pm1 pm1 pm1 pm1 pm1 pm1 pm1 ternary ternary ternary ternary ternary ternary ternary ternary int8 int8 int8 int8 int8 int8

64 64 128 128 192 192 256 256 512 512 64 64 128 128 192 192 256 256 64 64 128 128 192 192

8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8

100 400 100 400 100 400 100 400 100 400 100 400 100 400 100 400 100 400 100 400 100 400 100 400

1990 1224 1981 868 1974 599 1841 483 823 211 2030 1154 1974 808 1891 497 1502 394 2000 1007 1983 536 1139 291

4826 1279 3468 915 2453 644 1971 516 875 227 4464 1232 3117 852 2040 544 1652 437 4303 1176 2471 661 1523 400

543.0 3347.7 577.3 4735.8 621.0 6889.4 2230.7 8547.3 5014.1 19706.2 585.9 3557.6 593.8 5087.7 2165.2 8342.2 2754.8 10471.1 569.1 4082.3 667.9 7696.7 3685.3 14531.4

769 595 721 558 658 472 623 418 459 211 1025 560 674 546 624 418 578 353 997 583 658 446 517 289

5489.9 7378.3 6127.4 8124.5 6710.1 9394.0 7078.5 10966.7 9553.2 19639.4 4283.5 7814.1 6548.8 8717.9 7059.2 11325.1 7621.1 13794.6 4309.6 7512.1 6732.1 10090.5 8517.0 14467.6

8

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

TABLE 5. Unloaded round-trip time in µs (closed loop, one outstanding request, 100,000 requests per run), measured on the client. ‘‘Default’’ and ‘‘tuned’’ are the two host states of Section VI; the GPU-initiated server uses the batched send path. Bold: the lowest median of the row.

proxy, default

launch, tuned

proxy, tuned

GPUNetIO

format

n

B

K

launch, default p50

p99

p50

p99

p50

p99

p50

p99

p50

p99

int8 int8 int8 int8 int8 int8 int8 int8 int8 int8 int8 int8 int8 int8 pm1 pm1 pm1 pm1 ternary

64 64 64 64 128 128 128 128 192 192 192 192 128 128 64 128 256 512 256

8 8 8 8 8 8 8 8 8 8 8 8 1 32 8 8 8 8 8

50 100 200 400 50 100 200 400 50 100 200 400 100 100 100 100 100 100 100

71.7 110.4 283.1 422.0 109.9 183.2 407.0 680.5 245.0 357.6 591.0 1051.3 96.9 601.7 108.3 217.1 306.0 541.1 345.6

78.3 121.3 295.0 442.1 115.3 188.4 418.9 692.2 259.4 371.5 603.0 1067.6 102.2 615.4 174.2 234.3 320.3 555.8 358.6

73.9 113.5 278.4 436.8 108.6 264.8 406.5 685.6 245.2 364.7 595.8 1063.3 94.1 596.2 103.2 220.9 304.6 527.8 346.5

80.8 208.2 292.8 449.9 114.9 280.6 420.9 700.7 260.4 377.0 610.1 1077.1 100.2 626.4 109.4 236.8 320.6 550.2 363.1

67.7 106.6 182.9 338.3 104.5 178.5 313.8 586.5 154.2 266.0 498.4 958.1 92.6 510.5 104.4 129.4 215.7 450.9 252.7

70.6 108.9 192.1 341.1 107.4 180.8 316.4 590.3 157.6 270.0 502.2 967.9 95.8 513.4 107.2 132.5 219.0 453.9 256.2

70.9 110.0 188.5 345.5 104.7 175.6 314.4 595.6 155.2 271.5 505.0 971.1 90.2 507.0 100.0 130.6 215.3 437.0 256.3

75.6 114.5 192.9 350.0 109.1 180.3 319.3 600.5 160.1 276.4 509.8 976.3 94.6 529.3 104.3 135.0 220.1 451.5 263.3

62.1 101.6 181.4 340.8 95.9 166.4 307.1 588.0 146.4 262.8 500.2 959.8 82.9 500.4 92.1 124.6 206.8 431.5 249.6

65.1 103.5 183.5 343.8 98.1 168.2 310.4 590.6 149.3 265.7 504.2 964.4 86.8 525.9 94.7 130.1 210.5 447.4 255.5

TABLE 6. Decomposition of the unloaded round trip for the reference configuration (int8, n = 128, B = 8, K = 100), in µs, from the servers’ timestamps. ‘‘Server I/O’’ is the CPU time between recv() and sendto() minus the solve for the baselines, and the transmit wait for the GPU-initiated server. ‘‘Residual’’ is RTT minus the recorded server residence; its timing boundaries differ by server and it is not a direct wire-latency measurement.

server I/O

residual

server

p50

queue p99

p50

solve p99

p50

p99

p50

p99

p50

p99

launch, default proxy, default launch, tuned proxy, tuned GPUNetIO

0.0 0.0 0.0 0.0 2.0

0.0 0.0 0.0 0.0 2.1

151.7 152.6 151.7 152.9 149.9

152.0 155.0 152.0 155.4 151.1

16.1 11.3 16.1 12.0 0.8

17.0 14.8 17.0 15.5 2.0

15.5 101.1 10.7 10.5 13.5

20.5 115.9 12.6 12.5 15.7

183.2 264.8 178.5 175.6 166.4

188.4 280.6 180.8 180.3 168.2

182 μs

launch, default

260 μs

proxy, default launch, tuned

176 μs

proxy, tuned

176 μs

barely changed between the default and tuned host states (Section VI). The client is still CPU-based, so the load experiments below keep its receive path tuned. 2) Comparison with tuned baselines

166 μs

GPUNetIO 0

50

100

150

200

250

300

mean time per request (μs) solve (same kernel)

server I/O (CPU or GPU)

GPU queue

unmeasured residual

FIGURE 3. Mean per-request timing components for the reference request. Means are used because they add to the mean RTT; component medians do not generally sum to the median RTT. The residual includes uninstrumented communication and endpoint work.

large default-host gap therefore reflects the host configuration, and we use the tuned servers for the comparison. The GPU-initiated receive–solve–reply path bypasses the host network stack and its scheduling, and its unloaded RTT Manuscript for review

RTT

The GPU-initiated service has the lowest recorded median in 15 of 19 settings, improving on the faster tuned baseline by 1.5– 9.2 µs; the remaining four settings differ by 1.5–2.5 µs in favor of a baseline. Its recorded server residence is lower throughout the grid. For the reference request, the GPUNetIO residual is 13.5 µs, approximately 3 µs above the tuned baselines; the residuals use different timestamp boundaries (Section V). 3) Interpreting the timing decomposition

Relative to GPUNetIO, the tuned launch and proxy servers have median residence differences of 9.6 and 10.7 µs across the 19 settings. These differences combine dispatch, copies or mapped-memory access, and the different timestamp boundaries, and they are close to the launch-path overhead measured without the network (Section VII-B).

9

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

launch, default proxy, default

launch, tuned proxy, tuned

1.0

fraction of requests

fraction of requests

1.0

GPUNetIO

0.8 0.6 0.4 0.2 0.0

0.8 0.6 0.4 0.2 0.0

160

180

200

220

240

260

280

300

320

10

100

round-trip time (μs)

RTT − server residence (μs)

FIGURE 4. Distribution of the unloaded round-trip time (left) and of the round trip minus the server’s residence (right) for the reference configuration, 100,000 requests per server. On the default host the proxy’s residual lies in its higher-latency mode for 95% of the requests of this run; on the tuned host both CPU-mediated servers have a single mode near 10.5 µs. The GPU-initiated distribution is shown for comparison.

TABLE 7. Server residence (as stamped by the server) and the rest of the round trip (RTT minus residence), medians in µs, closed loop with one outstanding request.

Manuscript for review

launch, default

proxy, default

launch, tuned

proxy, tuned

GPUNetIO

format

n

B

K

server

rest

server

rest

server

rest

server

rest

server

rest

int8 int8 int8 int8 int8 int8 int8 int8 int8 int8 int8 int8 int8 int8 pm1 pm1 pm1 pm1 ternary

64 64 64 64 128 128 128 128 192 192 192 192 128 128 64 128 256 512 256

8 8 8 8 8 8 8 8 8 8 8 8 1 32 8 8 8 8 8

50 100 200 400 50 100 200 400 50 100 200 400 100 100 100 100 100 100 100

57.6 96.4 180.3 328.1 94.3 167.8 303.2 576.0 143.5 256.1 487.6 953.8 82.0 499.8 94.3 118.9 204.9 440.0 241.8

14.1 14.0 103.2 93.5 15.5 15.5 103.8 104.1 101.5 101.5 103.5 99.5 14.9 101.9 14.0 98.1 101.1 101.7 103.8

59.7 99.0 177.4 334.4 94.0 163.9 304.6 584.4 144.1 260.5 493.7 960.0 79.3 495.1 89.1 120.2 204.3 426.3 245.2

14.1 14.0 100.8 102.8 14.4 101.1 102.3 101.4 101.2 104.5 102.0 103.5 14.6 102.7 14.0 101.2 100.8 102.2 101.7

57.1 96.0 172.3 327.6 94.0 167.8 303.1 575.8 143.3 254.4 487.4 946.6 81.9 499.7 93.9 118.8 204.7 439.9 241.6

10.5 10.5 10.4 10.6 10.5 10.7 10.6 10.7 10.9 11.0 11.0 11.0 10.7 10.7 10.5 10.6 10.9 11.0 11.1

60.2 99.4 177.9 334.9 94.0 164.9 303.4 584.8 144.1 260.5 493.9 959.9 79.5 496.2 89.4 119.8 204.2 426.0 245.2

10.6 10.6 10.5 10.5 10.7 10.5 10.8 10.7 10.9 10.9 11.0 11.1 10.6 10.7 10.5 10.7 11.0 10.8 10.9

48.7 88.3 167.7 326.9 82.9 152.7 293.5 574.7 132.8 249.4 483.7 946.0 68.8 485.8 78.7 109.1 193.0 417.9 235.2

13.3 13.3 13.6 13.7 13.0 13.5 13.5 13.3 13.5 13.4 16.4 13.7 14.0 13.7 13.3 15.3 13.7 13.4 13.9

10

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

E. SCALING WITH SOLVE COST

Fig. 5 and Table 8 follow median RTT as the solve grows, through K at fixed n and through n at fixed K . With the tuned hosts, the CPU-mediated-to-GPUNetIO ratios remain near unity, spanning approximately 0.99–1.14. Dispatch savings matter most for short requests; as the common solver occupies more of the round trip, the ratios approach one. The larger separation under default host settings follows the slow receive mode discussed above. Without load, the main gains are therefore the freed server CPU core and lower host sensitivity; the larger latency difference appears under concurrent load. F. LATENCY UNDER LOAD

Tables 9 and 10 and Fig. 6 report the reference configuration under closed-loop concurrency and open-loop arrivals. With the default client, the residual RTT increases under load across all servers while their residence times remain comparatively stable. We therefore use the tuned client below. At 400–800 kreq/s offered, Table 10 and Fig. 6 report medians over repeated runs; Table 11 lists the counts and ranges. At 800 kreq/s offered, GPUNetIO completes 618– 628 kreq/s and the tuned proxy 620–668 kreq/s without loss; the client limits injection at this load. GPUNetIO’s median RTT is then 189 µs against 288 µs for the tuned proxy, and its 99th-percentile RTT 218 µs against 609 µs, reductions by factors of 1.5 and 2.8. At 600 kreq/s the corresponding values are 188 against 240 µs and 216 against 327 µs; at 400 kreq/s the two are still close (178 against 193 µs and 207 against 229 µs). Between 400 and 800 kreq/s offered, GPUNetIO’s median RTT rises by 6% and the proxy’s by 50%. At the two highest rates, where the client is saturated, individual runs of both servers occasionally show millisecondscale 99th percentiles; at 700 kreq/s this places the median 99th percentile for GPUNetIO (423.5 µs) slightly above that of the proxy (381 µs). The launch server reaches its four-stream limit much earlier, near 24 kreq/s. The server timestamps help explain this behavior. From 400 to 800 kreq/s offered, GPUNetIO’s median server residence stays at 156–159 µs and its 99th percentile at 180–186 µs, whereas the proxy’s rise from 173 to 235 µs and from 208 to 317 µs. In closed loop, the same pipeline also supports a higher completion rate than the proxy at the largest tested window. Batched transmission keeps reply submission from limiting the request service. At the highest offered load, per-reply sends complete 334 kreq/s with a median RTT of 24 ms and 41% request loss (Table 10). Batched descriptor posting avoids this bottleneck at the tested rates. VIII. STREAMING DENSE COUPLINGS A. WORKLOAD AND TRANSFER VOLUME

The request service sends small field updates to a resident model. When the coupling matrix exceeds GPU memory, the streaming solver receives the matrix tile by tile. Positions, momenta, signs, and row accumulators remain on the receiving GPU, where workers consume arriving tiles to evaluate the dSB coupling force. A remote matrix generator provides a Manuscript for review

reproducible source at large scale. We measure coefficient reception together with coupling evaluation, including reuse of each tile across replicas. A dense one-bit representation has logical size n2 /8 bytes. Frames contain a 64-byte header, including Ethernet, and Pf = 4032 payload bytes, representing Cb = 32,256 couplings from one row. The final column block is padded, so the transmitted payload is V (n) = n ⌈n/Cb ⌉ Pf . (9) For n = 107 , the logical matrix is 12.5 TB, whereas a pass carries 12.53952 TB in 3.11 × 109 frames. For n = 104 , padding expands 12.5 MB to 40.32 MB. All streaming rates below count transmitted payload, including padding. For payload rate R in bits/s, the transfer estimate is tpass ≈ 8V /R. A 4096-byte frame occupies 4120 bytes including FCS, preamble, and interpacket gap, giving an ideal payload ceiling of 400(4032/4120) = 391.46 Gb/s. On the current PCIe 5.0 server, the observed 307 Gb/s reaches approximately 78% of this payload ceiling while performing the coupling accumulation; the topology and bandwidth scaling directions are discussed in Section IX-C. B. GENERATOR AND RECEIVER

√ The streaming runs use seed 42, ∆t = 1.1, and ξ = 0.5/ n. The n = 107 , B = 1 run uses two updates; its B = 8 runs use ten. Reference comparisons at n = 104 and 105 use 20 updates. The generator runs on one GH200 GPU. Figure 7 summarizes the data and control paths. Synthesis blocks form frames in column-block-major order, and sender blocks post batches over one to eight queues. A pass ends with one PASS_END per queue; the next pass starts after an acknowledgement. Random couplings are generated symmetrically from a hash of the seed and the unordered pair of indices. Mattis couplings are Jij = ηi ηj for i ̸= j, where η is a hashed sign vector. The transmitted diagonal is +1 and is removed during evaluation. Mattis instances have a known ground state; random instances exercise generic coupling accumulation. The receiver has one cyclic queue containing 16,384 slots of 4096 bytes. A receive warp publishes bursts to 372 persistent worker blocks, each with 512 threads. Workers read frames in place, load the corresponding sign tile of each replica into shared memory, and accumulate XOR/popcount results into int32 row accumulators. At n = 107 and B = 8, the accumulators require 320 MB. The coupling sum is (b)

ai

(b)

= n − 2ci

(b)

− si ,

(10)

(b)

where ci is the mismatch count over the full padded row. Padding bits in both operands are zero. The diagonal correction has the opposite sign from Proposition 1(i), because that proposition stores a zero diagonal bit whereas the stream stores a one. For the measured range n ≤ 107 , the count fits int32 and |ai | ≤ n − 1 < 224 , so conversion of the coupling sum to FP32 remains exact. The energy reduction uses 64-bit integers. The host coordinates pass boundaries: it observes end markers, waits for published bursts to finish, and requests 11

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

1000

launch, tuned proxy, tuned

median round-trip time (μs)

median round-trip time (μs)

launch, default proxy, default

500 300 200

100 64

128

GPUNetIO

1000

500 300 200

100

256

64

128

steps K (int8 portfolio, n = 128, B = 8)

256

512

problem size n (pm1, B = 8, K = 100)

FIGURE 5. Median unloaded round-trip time versus integration steps K (left: INT8, n = 128, B = 8) and problem size n (right: PM1, B = 8, K = 100). Dashed lines and open markers denote default-host CPU-mediated servers; solid lines and filled markers denote tuned CPU-mediated servers and GPUNetIO. Both axes use logarithmic scales.

TABLE 8. Median round-trip time of the CPU-mediated servers relative to the GPU-initiated server, on the default and on the tuned host, with the GPU-initiated server’s residence as the measure of solve cost.

GPUNetIO

Manuscript for review

default host

tuned host

format

n

B

K

server p50 (µs)

launch

proxy

launch

proxy

int8 int8 int8 int8 int8 int8 int8 int8 int8 int8 int8 int8 int8 int8 pm1 pm1 pm1 pm1 ternary

64 64 64 64 128 128 128 128 192 192 192 192 128 128 64 128 256 512 256

8 8 8 8 8 8 8 8 8 8 8 8 1 32 8 8 8 8 8

50 100 200 400 50 100 200 400 50 100 200 400 100 100 100 100 100 100 100

48.7 88.3 167.7 326.9 82.9 152.7 293.5 574.7 132.8 249.4 483.7 946.0 68.8 485.8 78.7 109.1 193.0 417.9 235.2

1.16 1.09 1.56 1.24 1.15 1.10 1.33 1.16 1.67 1.36 1.18 1.10 1.17 1.20 1.18 1.74 1.48 1.25 1.38

1.19 1.12 1.53 1.28 1.13 1.59 1.32 1.17 1.68 1.39 1.19 1.11 1.14 1.19 1.12 1.77 1.47 1.22 1.39

1.09 1.05 1.01 0.99 1.09 1.07 1.02 1.00 1.05 1.01 1.00 1.00 1.12 1.02 1.13 1.04 1.04 1.05 1.01

1.14 1.08 1.04 1.01 1.09 1.06 1.02 1.01 1.06 1.03 1.01 1.01 1.09 1.01 1.09 1.05 1.04 1.01 1.03

12

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

TABLE 9. Closed loop with 1–64 outstanding requests, reference configuration.

outstanding

kreq/s

server residence (µs) p50

p99

p50

p99

p99.9

launch, default

1 4 16 64

5.5 21.4 24.0 24.1

167.8 161.3 162.6 162.4

168.8 169.8 166.5 166.4

183.2 178.6 667.0 2658.7

188.4 278.7 765.0 2760.1

197.6 324.1 772.6 2772.7

proxy, default

1 4 16 64

3.8 20.4 66.6 160.0

163.9 163.4 167.4 175.4

168.2 172.0 198.2 223.6

264.8 178.9 237.1 416.1

280.6 293.5 349.5 548.1

300.2 322.3 440.6 619.3

launch, tuned

1 4 16 64

5.7 22.7 24.4 24.4

167.8 162.0 162.0 162.0

168.7 170.2 165.6 165.7

178.5 175.9 654.3 2618.3

180.8 184.8 661.9 2629.5

184.4 190.2 665.7 2635.2

proxy, tuned

1 4 16 64

5.7 22.9 91.3 339.1

164.9 163.2 162.5 168.9

169.3 168.2 182.8 198.1

175.6 174.0 174.2 185.7

180.3 180.0 194.6 215.7

186.5 194.0 204.3 223.9

GPUNetIO

1 4 16 64

6.0 23.9 95.4 363.1

152.7 153.4 153.1 155.2

154.0 156.0 174.4 178.9

166.4 167.2 166.9 174.5

168.2 170.2 188.5 198.8

177.0 178.1 196.2 205.4

GPUNetIO, per-reply send

1 4 16 64

5.9 23.6 94.2 326.5

152.9 153.3 152.6 165.8

154.3 154.7 171.7 195.7

168.6 169.2 169.1 191.4

172.4 173.7 188.6 221.5

178.2 187.3 199.0 268.6

99th-percentile round-trip time (μs)

median round-trip time (μs)

server

10000 launch, default proxy, default launch, tuned proxy, tuned GPUNetIO

1000

GPUNetIO, per-reply send

10

100

configured Poisson rate (kreq/s)

client RTT (µs)

10000

1000

10

100

configured Poisson rate (kreq/s)

FIGURE 6. Median (left) and 99th-percentile (right) round-trip time against the offered Poisson load, reference configuration; crosses mark the first run with lost requests. The four-stream launch server completes about 24 kreq/s under overload. At the highest offered load, GPUNetIO with per-reply sends completes 334 kreq/s with 41% loss; batched GPUNetIO and the tuned proxy have no observed loss in the retained runs. Points of GPUNetIO and the tuned proxy at 400–800 kreq/s are medians over repeated runs. At the highest offered load, client injection falls below the configured rate.

Manuscript for review

13

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

TABLE 10. Open loop with Poisson arrivals, reference configuration, 100,000 requests per run. ‘‘Lost’’ are requests still unanswered when the client completes its final receive drain, which waits up to 2 s after sending the last request. ‘‘GPUNetIO, per-reply send’’ is the same server with the per-reply send path of the vendor sample (and the 5-µs receive timeout of the first runs). For GPUNetIO and the tuned proxy at 400–800 kreq/s, each entry is the median over repeated runs (Table 11). At 700 and 800 kreq/s offered the client itself emitted 596–681 kreq/s.

server

offered

achieved

lost

kreq/s

kreq/s

%

server residence (µs) p50

p99

p50

client RTT (µs) p99

p99.9

launch, default

10 50 100 200 400 800

10.1 24.1 23.9 23.6 23.0 22.5

0.0 51.8 76.0 88.0 94.0 96.1

161.1 162.0 162.1 162.2 162.3 162.8

174.4 166.2 166.6 166.5 166.7 167.4

249.8 11844.4 11984.0 11918.4 11887.1 11898.8

425.8 12397.1 12352.3 12185.4 12099.1 12124.2

549.6 12507.5 12433.8 12246.8 12131.7 12152.1

proxy, default

10 50 100 200 400 800

10.1 50.2 100.5 200.9 401.5 499.2

0.0 0.0 0.0 0.0 0.0 27.5

163.2 175.1 178.0 185.4 208.5 334.8

188.9 220.3 224.8 237.6 290.8 494.1

180.4 592.5 587.2 474.5 481.4 959.1

387.8 1479.1 1072.8 742.8 636.9 1169.2

501.7 1631.9 1181.6 795.0 691.2 1205.0

launch, tuned

10 50 100 200 400 800

10.1 24.3 24.2 23.9 23.4 22.8

0.0 51.3 75.7 87.8 93.9 96.1

160.9 161.6 161.8 161.7 161.8 161.8

170.0 165.0 165.2 165.7 165.6 165.7

172.3 11467.3 11504.1 11544.8 11585.3 11598.5

264.1 11589.9 11608.7 11617.6 11629.4 11646.0

315.7 11600.4 11618.8 11621.7 11635.8 11650.1

proxy, tuned

10 50 100 200 400 600 700 800

10.1 50.3 100.5 201.0 401.8 602.4 627.9 641.6

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

162.9 161.6 162.2 164.5 173.4 207.0 222.1 235.2

176.7 185.9 188.3 193.3 207.6 264.4 287.4 316.8

173.6 172.9 174.2 178.4 192.6 239.8 263.2 288.3

188.1 197.3 200.5 208.1 229.0 327.0 381.2 608.5

201.4 204.5 207.2 215.2 242.3 363.0 614.0 699.1

GPUNetIO

10 50 100 200 400 600 700 800

10.1 50.3 100.5 201.0 401.9 602.5 628.2 622.5

0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

153.4 153.1 153.3 154.4 156.2 159.0 159.0 158.7

156.8 166.6 178.2 178.7 183.6 186.1 185.7 180.1

166.7 166.5 167.7 170.9 177.8 188.1 188.8 188.6

171.2 180.8 193.0 195.7 207.0 215.8 423.5 217.9

178.3 195.6 200.1 202.4 213.7 235.3 508.3 222.9

GPUNetIO, per-reply send

10 50 100 200 400 800

10.1 50.3 100.5 201.0 401.5 333.7

0.0 0.0 0.0 0.0 0.0 40.9

153.4 152.9 153.4 157.3 249.9 23323.4

160.6 177.7 180.5 186.0 484.9 24793.2

169.1 169.5 171.4 178.2 311.5 24307.5

176.9 194.6 199.0 209.1 537.0 25902.0

195.3 203.7 206.7 218.8 564.6 25946.6

TABLE 11. Repeated open-loop runs of the reference configuration with the current code (GPUNetIO and tuned proxy). Client RTT range over the runs, median in parentheses. One GPUNetIO run at 700 kreq/s, in which the client dropped 1.1% of the replies although the server transmitted all of them, is excluded.

server

offered

runs

kreq/s

Manuscript for review

achieved

client RTT (µs), range over runs (median)

kreq/s

p50

p99

p99.9

GPUNetIO

400 600 700 800

4 3 2 4

402 602 623–633 618–628

177–179 (178) 187–188 (188) 188–190 (189) 187–189 (189)

204–209 (207) 216–321 (216) 206–641 (423) 195–3127 (218)

213–214 (214) 220–531 (235) 216–800 (508) 199–3326 (223)

proxy, tuned

400 600 700 800

4 3 3 4

402 602 617–646 620–668

191–194 (193) 236–246 (240) 262–286 (263) 239–314 (288)

226–232 (229) 313–347 (327) 354–1427 (381) 408–876 (609)

236–247 (242) 353–491 (363) 387–1598 (614) 443–1032 (699)

14

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

C. VALIDATION

Validation combines stream accounting, planted-solution checks, and reference execution. The receiver compares parsed stream-frame totals with the end markers and records dropped bursts. Mattis instances provide the exact ground-state value 2E = −n(n − 1) for checking the accumulated energy and solver integration. At n = 105 , the Mattis and random energy/checksum records match a resident GPU reference byte for byte; at n = 104 , the random records also match a CPU reference. These records contain the best-replica energy and a 64-bit spin-bitmap checksum at each recorded step.

GH200 generator GPU Synthesize tiles; batch sends

Acknowledgement / next pass

400 GbE link 4032 payload bytes / frame

RTX PRO 6000 GPU RX buffer and receive warp Workers: XOR / popcount coupling sums

D. MEASURED RESULTS Host pass control

Row accumulators Evaluate energy; update state* Refresh resident signs

pass complete

GPU sends PASS_ACK Release next matrix pass

*Final pass: energy evaluation only; no state update.

FIGURE 7. Streaming implementation. Solid arrows show coefficient delivery and GPU computation. Dashed arrows show pass-boundary coordination and the GPU-sent acknowledgement that releases the generator for the next pass. Host control follows reception and draining of outstanding work. Each of the K updates requires a matrix pass, followed by one final-energy pass.

an in-kernel update through mapped control memory. Workers evaluate the energy of the incoming signs, advance the dynamical state when an update is requested, and compute a checksum of the resulting sign bitmap. Thus an updating record contains pre-update energy and post-update checksum. A run of K updates requires K + 1 matrix passes, with the last pass evaluating the final energy without advancing the state. Packet reception and coupling accumulation execute on the GPU, with host coordination at pass boundaries. The reported receive-pass interval ends when the host observes the end markers, which it polls every 2 ms. Full iteration time also includes draining outstanding work, updating the state, and acknowledging the pass; these control operations use host polling. Work-ring acknowledgements protect work-item reuse. The cyclic packet buffer is managed separately, so workers must finish consuming each burst before the NIC wraps around. At the measured rate, the buffer holds approximately 1.7 ms of arriving frames. The evaluated implementation uses mapped host memory for its receive/control counters.

Manuscript for review

Table 12 reports the streaming measurements. The n = 107 , B = 1 Mattis run passes all recorded counter checks and reaches 2E = −99,999,990,000,000 after its first update. Its later passes take 326.733 s at approximately 307.03 Gb/s. The 12.5-TB matrix is never stored on the receiver. At n = 106 , the reported runs with B = 1 and B = 8 and one to eight sender queues have similar rates near 307 Gb/s. This supports amortizing matrix transmission across the tested replicas. The runs at n = 105 and 104 serve as reference checks. At n = 104 , a pass carries 40 MB, equivalent to about 1 ms at 307 Gb/s. The host polls for end markers every 2 ms, so the recorded 2.1-ms interval includes a substantial polling delay. At n = 107 with B = 8, the Mattis run achieves 306.94–307.01 Gb/s and all eight replicas reach the known ground-state energy. The random instance achieves 306.20– 306.41 Gb/s; its best final energy is 2E = −41,168,653,412 after ten updates. Each update requires a 327-s pass. The random-instance run tests coupling evaluation over ten updates; it does not establish convergence. Every pass contains the expected 3.11 × 109 parsed stream frames, with zero reported missing stream frames and zero dropped bursts. At B = 8, each received coefficient contributes to eight coupling evaluations. The similar payload rates for one and eight replicas show that this reuse amortizes matrix transfer across the tested replica batch. IX. DISCUSSION A. REQUEST-RATE AND MODEL-SIZE SCALING

The resident-model and streaming modes place different demands on the same communication path. Small field updates require low scheduling overhead across many concurrent solves. Dense coupling tiles require sustained bandwidth while the solver state stays resident. In both cases, GPU workers consume coefficients directly from the receive path. For resident models, the main gain appears under load. The unloaded latency is close to the tuned proxy’s, but GPUNetIO’s median RTT rises by 6% from 400 to 800 kreq/s offered, compared with 50% for the proxy. Each CPU-mediated baseline dedicates a polling core to receiving requests and returning results; GPUNetIO performs that work on the device after setup. Keeping packet handling on the GPU also reduces dependence on host scheduling and interrupt settings. 15

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

TABLE 12. Streaming observations from the preserved records. Rate and receive-pass time are ranges after excluding pass 0, whose interval includes the start-up of the generator, and, for n ≤ 105 , pass 1, which contains a one-time delay of about 10 ms. Payload includes padding. All listed passes have matching parsed stream counts and zero dropped bursts. Reference matches concern energy/checksum records.

Instance

n

B

Nq

Payload (Gb/s)

Receive pass (s)

Evidence

Mattis Mattis random Mattis Mattis random Mattis random random

107 107 107 106 106 106 105 105 104

1 8 8 8 1 8 8 8 8

4 4 4 1–8 4 4 4 4 4

307.0 306.9–307.0 306.2–306.4 307.0–307.1 307.0 307.0–307.1 308.8–309.4 308.8–311.5 153.6–153.7

326.73 326.75–326.83 327.39–327.61 3.361–3.362 3.362 3.361–3.362 0.0417–0.0418 0.0414–0.0418 0.0021

2E = −n(n − 1) at step 1 2E = −n(n − 1) at step 1 2E = −4.12 × 1010 at step 10 2E = −n(n − 1) 2E = −n(n − 1) energy reported GPU records match GPU records match CPU records match

The request-per-block mapping allows independent solves to run concurrently, and batched descriptor posting drains their completed replies. Both are needed to keep communication and solver execution at comparable rates. For model-size scaling, streaming separates quadratic coefficient storage from the resident solver state. A dense one-bit matrix requires O(n2 ) storage, whereas the dynamical state, sign bitmaps, and accumulators require O(Bn). At n = 107 , the 12.5-TB logical matrix is consumed through a 64-MiB packet buffer, with 320 MB of row accumulators for eight replicas. Device memory holds the solver state and receive buffers; the time per update is then governed by coefficient-delivery bandwidth. Each incoming tile is reused across replicas, but every update requires another matrix pass. B. EXTERNAL COEFFICIENTS AND CHANGING MODELS

Streaming separates coefficient production from solver execution. Potential sources include remote memory, storage arrays, and simulation or learning processes that construct successive models. An online application can update hℓ while reusing Jℓ , replace a resident model between requests, or supply a new sequence of dense tiles when the full model cannot be retained. When the coefficients admit a low-rank or sparse representation, preserving that structure offers an additional reduction in model traffic. The measured 306–307 Gb/s corresponds to approximately 38.3 GB/s of coefficient payload. For scale, a Samsung 9100 PRO PCIe 5.0 NVMe SSD specifies sequential reads up to 14.8 GB/s [31]; feeding the measured stream from such storage would call for parallel drives. The coefficient source must sustain this bandwidth to feed the receiver continuously. C. TOPOLOGY AND BANDWIDTH SCALING

The current PCIe 5.0 server sustains 307 Gb/s. Its reported SYS/NODE GPU-to-NIC paths traverse host interconnects, making GPU/NIC placement one target for improving this rate. A direct PCIe-switch path, sufficient upstream bandwidth, and balanced queue placement are candidates for reducing the gap to the 400-GbE link’s 391.46-Gb/s payload ceiling. The similar rates for one and eight replicas suggest room to increase delivery bandwidth while retaining replica reuse. Manuscript for review

For payload volume V and rate R, streaming incurs a transfer term 8V /R per pass. A trajectory of K updates uses K + 1 passes, including final-energy evaluation. The present n = 107 configuration takes approximately 327 s per receive pass. At the ideal payload ceiling of the same 400-GbE link, the transfer term is approximately 256 s; an 800-GbE link with the same frame efficiency reduces this ideal term to approximately 128 s. These estimates describe bandwidth scaling at fixed model size. A distributed extension could partition matrix rows across G receivers, each sustaining payload rate R. With sufficient aggregate source bandwidth, the ideal transfer term becomes 8V /(GR), plus synchronization and sign exchange; each replica’s sign bitmap contains n/8 bytes. PCIe 6.0 doubles the nominal bandwidth of PCIe 5.0 at the same lane count, increasing the per-lane data rate from 32 to 64 GT/s [32]. Combining this wider host interconnect with 800-Gb/s networking expands the available coefficientdelivery bandwidth. Newer ConnectX platforms provide this combination, and the BlueField-4 data processing unit (DPU) integrates PCIe 6.0, an 800-Gb/s network interface, programmable processing, and infrastructure acceleration [33], [34]. The coefficient source, PCIe fabric, and GPU consumer must sustain the higher network rate. A DPU could handle communication control alongside GPU coupling evaluation. D. SCALING BEYOND PER-PACKET CPU CONTROL

The slowdown in single-thread processor performance growth strengthens the case for reducing frequent host coordination as network and PCIe bandwidth increase [14]. This concern also applies after kernel bypass: in NVIDIA’s Aerial example, CPU coordination became a bottleneck as the number of radio units increased in a DPDK GPUdev pipeline, motivating GPUcontrolled reception with GPUNetIO [15]. Queue service, completion handling, and CPU–GPU notification all contribute to this coordination cost. At the 4120-byte wire occupancy used here, 400 and 800 Gb/s correspond to approximately 12.1 and 24.3 million packets/s, respectively. A serial packetprocessing stage would therefore have only 82.4 and 41.2 ns per packet to keep pace. Batching, multiple cores, and offloads help host-mediated designs meet this demand. GPU initiation moves frequent queue service and packet consumption to device-side agents, where each received tile is reused across 16

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

parallel solver replicas. This reduces CPU–GPU notifications, while throughput still depends on topology, batching, queue placement, and GPU occupancy. E. IMPLEMENTATION SCOPE

The request evaluation uses one platform and two socketbased baselines. A matched DPDK/GPUdev baseline using the same solver, packet format, and batching policy is a useful next comparison to separate GPU-initiation gains from kernel-bypass gains. Deployment extensions include synchronized replacement of in-use models, retransmission and per-frame sequence validation, and distributing streaming work across multiple receivers. The current request bounds remain one frame per solve, one 512-KiB model slot, and the register/shared-memory limits of Section IV. X. CONCLUSION

We implemented dSB as a GPU-initiated request service and a solver for streamed dense couplings. The request service keeps receive, solve, and reply operations on the GPU, using one worker block per request and batched transmission. It requires no dedicated server CPU data-path core. From 400 to 800 kreq/s offered, median RTT rises by 6%; at the highest tested load, median and 99th-percentile latencies are lower than those of the tuned persistent CPU proxy by factors of 1.5 and 2.8 across repeated runs. The streaming solver retains O(Bn) state and consumes the O(n2 ) coefficients tile by tile. It evaluates ten-millionvariable binary coupling matrices at 306–307 Gb/s for one or eight replicas. Known Mattis ground states and smallerscale reference comparisons validate the execution path. The two modes extend dSB along different axes: concurrent decisions with a resident model, and larger models supplied as coefficient streams. Faster interconnects and multiple receivers could reduce the matrix-pass time, provided the coefficient source and GPU consumers sustain the aggregate bandwidth. CODE AND DATA AVAILABILITY

The implementation and experimental data are available in the author’s GitHub repository at https://github.com/ yaochengchen/dsb-gpunetio. The companion GPU solver is available at https://github.com/yaochengchen/dSB_GPUs_ 3Million_N [16]. ACKNOWLEDGMENT

The author acknowledges support from FAPESP (grants 25/25586-5 and 2021/02120-0). Generative AI tools (OpenAI ChatGPT, Anthropic Claude and Google Gemini) were used to assist with software development, and ChatGPT for language editing; the author reviewed and verified all AI-assisted output and is responsible for the manuscript and the software. CONFLICT OF INTEREST

The author declares no conflicts of interest. Manuscript for review

REFERENCES [1] A. Lucas, ‘‘Ising formulations of many NP problems,’’ Frontiers in Physics, vol. 2, art. 5, 2014. https://doi.org/10.3389/fphy.2014.00005. [2] F. Glover, G. Kochenberger, and Y. Du, ‘‘Quantum Bridge Analytics I: A tutorial on formulating and using QUBO models,’’ 4OR, vol. 17, no. 4, pp. 335–371, 2019. https://doi.org/10.1007/s10288-019-00424-y. [3] H. Goto, K. Tatsumura, and A. R. Dixon, ‘‘Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems,’’ Science Advances, vol. 5, no. 4, eaav2372, 2019. https://doi.org/10.1126/sciadv.aav2372. [4] H. Goto, K. Endo, M. Suzuki, Y. Sakai, T. Kanao, Y. Hamakawa, R. Hidaka, M. Yamasaki, and K. Tatsumura, ‘‘High-performance combinatorial optimization based on classical mechanics,’’ Science Advances, vol. 7, no. 6, eabe7953, 2021. https://doi.org/10.1126/sciadv.abe7953. [5] K. Tatsumura, R. Hidaka, J. Nakayama, T. Kashimata, and M. Yamasaki, ‘‘Pairs-trading system using quantum-inspired combinatorial optimization accelerator for optimal path search in market graphs,’’ IEEE Access, vol. 11, pp. 104406–104416, 2023. https://doi.org/10.1109/ACCESS.2023.3316727. [6] K. Tatsumura, Y. Hamakawa, M. Yamasaki, K. Oya, and H. Fujimoto, ‘‘Enhancing vehicle-mountable multiple object tracking systems with embeddable Ising machines,’’ Nature Communications, vol. 17, art. 584, 2026. https://doi.org/10.1038/s41467-025-67282-7. [7] Toshiba, ‘‘Toshiba brings simulated bifurcation machine to resource control technology for 5G base stations,’’ Apr. 2024. [Online]. Available: https://www.global.toshiba/ww/ technology/corporate/rdc/rd/topics/24/2404-02.html. [8] T. Kashimata, Y. Hamakawa, M. Yamasaki, and K. Tatsumura, ‘‘Real-time black-box optimization for dynamic discrete environments using embedded Ising machines,’’ Nature Communications, vol. 17, art. 8085, 2026. https://doi.org/10.1038/s41467-026-76069-3. [9] M. Silberstein, S. Kim, S. Huh, X. Zhang, Y. Hu, A. Wated, and E. Witchel, ‘‘GPUnet: Networking abstractions for GPU programs,’’ ACM Transactions on Computer Systems, vol. 34, pp. 1–31, 2016. https://doi.org/10.1145/2963098. [10] F. Daoud, A. Watad, and M. Silberstein, ‘‘GPUrdma: GPU-side library for high performance networking from GPU kernels,’’ in Proceedings of the 6th International Workshop on Runtime and Operating Systems for Supercomputers (ROSS), pp. 1–8, 2016. https://doi.org/10.1145/2931088.2931091. [11] C. Jung, S. Kim, I. Yeom, H. Woo, and Y. Kim, ‘‘GPU-Ether: GPU-native packet I/O for GPU applications on commodity Ethernet,’’ in IEEE Conference on Computer Communications (INFOCOM), pp. 1–10, 2021. https://doi.org/10.1109/INFOCOM42981.2021.9488699. [12] C. Jung, S. Kim, Y. Kim, and I. Yeom, ‘‘Virtualizing GPU direct packet I/O on commodity Ethernet to accelerate GPU-NFV,’’ Journal of Network and Computer Applications, vol. 206, art. 103480, 2022. https://doi.org/10.1016/j.jnca.2022.103480. [13] NVIDIA, DOCA GPUNetIO, DOCA SDK documentation. [Online]. Available: https://docs.nvidia.com/doca/sdk/. Accessed Sep. 29, 2026. [14] J. L. Hennessy and D. A. Patterson, ‘‘A new golden age for computer architecture,’’ Communications of the ACM, vol. 62, no. 2, pp. 48–60, 2019. https://doi.org/10.1145/3282307. [15] E. Agostini, ‘‘Inline GPU packet processing with NVIDIA DOCA GPUNetIO,’’ NVIDIA Technical Blog, Dec. 14, 2022, updated May 8, 2023. [Online]. Available: https://developer.nvidia.com/blog/inline-gpu-packetprocessing-with-nvidia-doca-gpunetio/. Accessed Oct. 4, 2026. 17

Y. Chen: GPU-Initiated Discrete Simulated Bifurcation

[16] Y. Chen, ‘‘dSB_GPUs_3Million_N,’’ GitHub repository. [Online]. Available: https://github.com/yaochengchen/dSB_GPUs_3Million_N. Accessed Oct. 4, 2026. [17] K. Tatsumura, M. Yamasaki, and H. Goto, ‘‘Scaling out Ising machines using a multi-chip architecture for simulated bifurcation,’’ Nature Electronics, vol. 4, pp. 208–217, 2021. https://doi.org/10.1038/s41928-021-00546-4. [18] T. Kashimata, M. Yamasaki, R. Hidaka, and K. Tatsumura, ‘‘Efficient and scalable architecture for multiple-chip implementation of simulated bifurcation machines,’’ IEEE Access, vol. 12, pp. 36606–36621, 2024. https://doi.org/10.1109/ACCESS.2024.3374089. [19] K. Tatsumura, R. Hidaka, J. Nakayama, T. Kashimata, and M. Yamasaki, ‘‘Real-time trading system based on selections of potentially profitable, uncorrelated, and balanced stocks by NP-hard combinatorial optimization,’’ IEEE Access, vol. 11, pp. 120023–120033, 2023. https://doi.org/10.1109/ACCESS.2023.3326816. [20] M. Yamasaki, R. Hidaka, Y. Sakai, J. Yamaguchi, Y. Nakamura, H. Goto, and K. Tatsumura, ‘‘Live demonstration: Capturing short-lived currency arbitrage opportunities with a simulated bifurcation algorithm-based trading system,’’ in IEEE International Symposium on Circuits and Systems (ISCAS), 2020. https://doi.org/10.1109/ISCAS45731.2020.9180679. [21] K. Tatsumura, ‘‘Large-scale combinatorial optimization in real-time systems by FPGA-based accelerators for simulated bifurcation,’’ in Proceedings of the 11th International Symposium on Highly Efficient Accelerators and Reconfigurable Technologies (HEART), pp. 1–6, 2021. https://doi.org/10.1145/3468044.3468045. [22] E. Agostini, D. Rossetti, and S. Potluri, ‘‘GPUDirect Async: Exploring GPU synchronous communication techniques for InfiniBand clusters,’’ Journal of Parallel and Distributed Computing, vol. 114, pp. 28–45, 2018. https://doi.org/10.1016/j.jpdc.2017.12.007. [23] NVIDIA, GPUDirect RDMA, CUDA documentation. [Online]. Available: https://docs.nvidia.com/cuda/gpudirect-rdma/. Accessed Sep. 29, 2026. [24] K. Gupta, J. A. Stuart, and J. D. Owens, ‘‘A study of persistent threads style GPU programming for GPGPU workloads,’’ in Innovative Parallel Computing (InPar), 2012. https://doi.org/10.1109/InPar.2012.6339596. [25] K. Hamidouche, J. Bachan, P. Markthub et al., ‘‘GPU-initiated networking for NCCL,’’ arXiv:2511.15076, 2025. https://arxiv.org/abs/2511.15076. [26] D. Unat, I. Turimbetov, M. K. T. Issa, D. Sağbili, F. Vella, D. De Sensi, and I. Ismayilov, ‘‘The landscape of GPU-centric communication,’’ ACM Computing Surveys, vol. 58, no. 12, pp. 1–36, 2026. https://doi.org/10.1145/3813799. [27] DPDK Project, Data Plane Development Kit. [Online]. Available: https://www.dpdk.org/. Accessed Sep. 29, 2026. [28] E. Agostini, ‘‘Boosting inline packet processing using DPDK and GPUdev with GPUs,’’ NVIDIA Technical Blog, Apr. 28, 2022. [Online]. Available: https://developer.nvidia.com/blog/optimizing-inline-packetprocessing-using-dpdk-and-gpudev-with-gpus/. Accessed Oct. 4, 2026. [29] J. W. Romein, ‘‘Breaking the I/O barrier: 1.2 Tb/s Ethernet packet processing on a GPU,’’ in Euro-Par 2025: Parallel Processing, Lecture Notes in Computer Science, vol. 15902, pp. 225–238, 2025. https://doi.org/10.1007/978-3-031-99872-0_16. [30] H. Markowitz, ‘‘Portfolio selection,’’ The Journal of Finance, vol. 7, no. 1, pp. 77–91, 1952.

Manuscript for review

[31] Samsung, ‘‘9100 PRO: Gen5 PCIe NVMe SSD,’’ product specifications. [Online]. Available: https://images.samsung. com/is/content/samsung/p6pim/us/mz-vap2t0bam/specsheet/us-specsheet-mz-vap2t0b-am-548821014.pdf. Accessed Oct. 3, 2026. [32] PCI-SIG, ‘‘PCI Express 6.0 specification.’’ [Online]. Available: https://pcisig.com/pci-express-6.0-specification. Accessed Oct. 3, 2026. [33] NVIDIA, ‘‘Ethernet SuperNICs,’’ product documentation. [Online]. Available: https://www.nvidia.com/enus/networking/products/ethernet/supernic/. Accessed Oct. 3, 2026. [34] NVIDIA, ‘‘Scaling agentic AI factories through extreme co-design with NVIDIA BlueField,’’ NVIDIA Technical Blog. [Online]. Available: https://developer.nvidia.com/blog/scaling-agentic-ai-factoriesthrough-extreme-co-design-with-nvidia-bluefield/. Accessed Oct. 3, 2026. YAOCHENG CHEN received the B.S. degree in physics and electronics from Central South University in 2015, and the M.S. and Ph.D. degrees in physics from National Taiwan University in 2017 and 2025, respectively. He is a FAPESP postdoctoral fellow at IAG-USP, Universidade de São Paulo, Brazil. Previously, he was an algorithm scientist at Jingdong, China; a project researcher and visiting scholar at the Quantum Science Center of GuangdongHong Kong-Macao Greater Bay Area; and a KASI postdoctoral fellow at the Korea Astronomy and Space Science Institute (KASI), Republic of Korea. His interests include high-energy astroparticle physics, extragalactic astrophysics, GPU computing and FPGA systems. He is a member of the ASTRI Mini-Array and TAROGE collaborations, the CTAO Consortium and the Space Array team. His previous collaborations include ANITA/PUEO, ARA and BURSTT; he was also a member of the ALMA Total Power GPU Spectrometer (TPGS) team. ORCID: https://orcid.org/0000-0002-8967-4911.

18

Related documents

• GPU-Initiated Discrete Simulated Bifurcation: Low-Latency Requests and Streaming Dense Couplings#1215718
• Empirical Analysis of GPU Frequency Behavior Under ML Workloads#353030
• EMA: Elastic and Performance Transparent Memory Across GPUs#1063251
• Echo: Merging Host Device Buffers to Avoid Redundant Data Movement on Unified Memory SoCs#1335821
• Echo: Merging Host Device Buffers to Avoid Redundant Data Movement on Unified Memory SoCs#1169660