OPT2026: 18th Annual Workshop on Optimization for Machine Learning
Distributed Linear Programming on GPU Clusters at Extreme Scale Arnaud Deza* Santanu Dey Pascal Van Hentenryck
ADEZA 3@ GATECH . EDU SANTANU . DEY @ ISYE . GATECH . EDU PVH @ GATECH . EDU
arXiv:2609.09108v1 [math.OC] 8 Sep 2026
H. Milton Stewart School of Industrial and Systems Engineering Georgia Institute of Technology, Atlanta, Georgia, USA
Abstract Large linear programs can exceed the memory of a single compute node. Although first-order methods replace sparse factorizations with GPU-suited matrix–vector products, other solver phases can reintroduce a single-node memory limit. We present S HARD LP, a distributed GPU LP solver that keeps the matrix and primal–dual state partitioned from sharded input through solution output. On the Google PDLP benchmark, S HARD LP reaches the published criterion on nine of eleven instances, compared with eight in the published CPU PDLP study. On the largest benchmark, eight H200 GPUs solve a 1.185-billion-variable, 6.338-billion-nonzero LP in 9.9 minutes; the published CPU experiment reports 21.06 hours on different hardware. Beyond this benchmark, separately checked multi-node solves reach up to 13.604 billion variables and 40.807 billion nonzeros, while validated executions span up to 76 GPUs across 29 compute nodes. For column-partitioned solves, support-aware communication skips GPUs that store no coefficients for a row; on an LP with 2.76 billion nonzeros, it cuts modeled communication by 92.97% and improves solver time by 1.27×–1.52×.
1. Introduction Linear programming (LP) is a core model in large-scale optimization. Modern simplex and interiorpoint methods are highly effective, but at very large scale their factorization-based linear algebra can become memory intensive and difficult to parallelize [5]. PDHG-based first-order methods have a different computational profile: their dominant operations are sparse matrix–vector products, projections, and vector operations, which map naturally to GPUs and distributed memory. Primal–dual hybrid gradient (PDHG), also known as the Chambolle–Pock method, is the firstorder method underlying Primal-Dual Linear Programming (PDLP) [3, 12]. PDLP combines the basic iteration with diagonal scaling, presolve, adaptive step sizes, restart, and feasibility polishing [4, 5]. We use the eleven released instances from the large-scale Google PDLP study [5] as our principal benchmark. Its largest instance contains 1.185 billion variables and 6.338 billion nonzeros. Recent GPU implementations show that PDHG maps well to accelerators, but scaling a general sparse LP solver across multiple compute nodes remains much less developed. Among the closest systems, D-PDLP [19] distributes the two PDHG matrix products across eight H100 GPUs within a single node, while MPAX [21] leaves efficient distributed sparse-data sharding as future work. Our question is therefore not only how to distribute the PDHG iteration, but how to turn it into a complete * Corresponding author: [email protected]
D ISTRIBUTED LP ON GPU C LUSTERS
solver when neither the matrix nor the primal–dual state fits on one node. A broader comparison with prior GPU and distributed optimization systems is given in Appendix E. We make three contributions. First, S HARD LP keeps the matrix and primal–dual state distributed during every solver stage—persistent ownership—so no phase needs a full copy on one node. Second, it reaches the Google PDLP study’s published criterion on nine of eleven benchmark instances and separately validates multi-node solves up to 13.604 billion variables and 40.807 billion nonzeros. Third, we show that communication can erase multi-GPU speedups and introduce support-aware communication, which skips GPUs with no coefficients for a row. On Design Match, this cuts modeled communication by 92.97% and improves solver time by 1.27×–1.52×.
2. From a distributed iteration to a distributed solver 2.1. PDHG and the D-PDLP matrix decomposition Consider a linear program in interval form, with c ∈ Rn , min c⊤ x x∈X
Ax ∈ S,
s.t.
X = [ℓv , uv ],
A ∈ Rm×n .
S = [ℓc , uc ],
(1)
Here, X gives the variable bounds and S the lower and upper bounds on each row activity. At a high level, PDHG alternates a primal update that uses A⊤ y and a dual update that uses Ax. Suppressing diagonal preconditioning, following [3, 12], define the base map Tτ,σ (w) = (b x, yb) for w = (x, y) by x b = projX x − τ (c − A⊤ y) , (2) yb = ΦS,σ (y, A(2b x − x)) .
(3)
Here τ and σ are the primal and dual step sizes, projX clips to the variable bounds, and ΦS,σ denotes the coordinatewise dual update for the row bounds. S HARD LP does not introduce a new PDHG variant. Its optimization logic follows cuPDLPx [22]; within each restart epoch, it applies the reflected-Halpern update wk+1 =
k+1 1 (1 + γ)Tτ,σ (wk ) − γwk + w0 , k+2 k+2
γ ∈ [0, 1],
(4)
where γ controls the reflection and a restart resets the anchor w0 and epoch counter. Scaling, step selection, primal-weight adaptation, and restart follow Lu et al. [22]; their control logic uses scalar reductions. To distribute the two matrix products, S HARD LP adopts D-PDLP’s two-dimensional matrix decomposition [19], with one Message Passing Interface (MPI) process, or rank, per GPU. We arrange the p = RC ranks in an R × C process grid, with R constraint (matrix-row) partitions and C variable (matrix-column) partitions. For each experiment, we choose the R × C grid before sharding so that the estimated matrix block and replicated vector slices fit in each GPU’s memory. Row set Ir and column set Jc define the block Arc = AIr ,Jc . The corresponding vector blocks are yr = yIr and xc = xJc . The two products are
⊤
A y
Jc
=
R X
A⊤ rc yr ,
(Ax)Ir =
r=1
C X c=1
2
Arc xc .
(5)
D ISTRIBUTED LP ON GPU C LUSTERS
(a) Rooted lifecycle one compute node full A, x, y input · scale · recover scatter gather PDHG on GPUs 0 3 1 2
per-node memory ceiling
(b) Persistent ownership 3 compute nodes × 2 GPUs/node; one MPI rank per GPU node 0 node 1 node 2 xJ0 xJ1 xJ2 y I0 y I1
input
O(nnz(A) + m + n)
A00
A01
A02
(Ax)I0
A10
A11
A12
(Ax)I1
(AT y)J0
(AT y)J1
(AT y)J2
PDHG
scale
polish
row reductions assemble Ax
column reductions assemble AT y output recover
same ownership throughout
Figure 1: Why persistent ownership matters. (a) Centralized solver phases impose an O(nnz(A) + m + n) per-node memory ceiling even when PDHG is distributed. (b) In a 2 × 3 process grid, rank (r, c) stores Arc ; row/column reductions form Ax and A⊤ y while ownership remains sharded from input through output. Rank (r, c) stores Arc . Ranks sharing c form a process column: they replicate xc and combine partial A⊤ y products. Ranks sharing r form a process row: they replicate yr and combine partial Ax products. Storage per rank is O(nnz(Arc ) + n/C + m/R). Up to this point, the decomposition is D-PDLP’s. S HARD LP’s contribution is to preserve this ownership beyond the two matrix products and throughout the rest of the solve, as illustrated in Figure 1. 2.2. Persistent ownership beyond the iteration Distributing only Ax and A⊤ y is not enough: if scaling, polishing, recovery, or another solver stage gathers the full problem on one node, that node still sets the memory limit. S HARD LP therefore keeps the same partition throughout the run. Each rank reads a prepartitioned prepared shard; constructing these shards from a monolithic model is an offline step excluded from timing. Scaling, restart, termination, polishing, recovery, and output remain distributed. General presolve is disabled. A distributed pass fixes variables forced to zero by globally singleton zero-equality rows; in both arms of the RQ3 Design Match communication experiments, this removes 1.123 billion nonzeros (40.69%), leaving 1.637 billion nonzeros. The RQ1 experiments do not use this reduction. Appendix C.1 gives the full lifecycle and recovery details.
3. Computational Results We organize the experiments around three questions: RQ1: Can S HARD LP solve the Google PDLP benchmark? RQ2: How far can it scale beyond that benchmark? RQ3: Does communication limit multi-GPU scaling, and when does support-aware communication help? All experiments ran on Georgia Tech’s PACE Phoenix Cluster, with one MPI rank per GPU. A compute node is one physical server. The validated runs use 1–76 GPUs and up to 29 nodes; Appendix A gives the hardware, interconnect, software, and memory details. Validation and timing. For RQ1, we use the Google study’s published scaled-LP criterion [5]. For every other reported solve, a separate checker evaluates the exported primal–dual solution against the 3
D ISTRIBUTED LP ON GPU C LUSTERS
original unscaled LP. Acceptance requires all nine quantities in Equation (6), covering normalized primal feasibility, stationarity, dual-sign admissibility, and relative primal–dual gap, to be finite and at most 10−6 . Solver time covers optimization and feasibility polishing when enabled; end-to-end time additionally includes sharded input, initialization, recovery, and output. Offline shard construction and the final independent check are excluded. Appendix C gives the full protocol. 3.1. RQ1: Google PDLP benchmark S HARD LP – SECONDS
Problem size Instance
m
n
Published CPU – HOURS
nnz(A) GPUs/nodes solver end-to-end Google PDLP
TSP-Gaia-100M 162.935M 1.185B 6.338B Design Match 22.000M 40.000M 2.760B QAP-THO-150 6.705M 249.784M 1.006B World Shipping 15.304M 228.868M 688.659M Mediterranean 7.491M 208.479M 628.927M Production Inventory 4.651M 18.271M 500.050M TSP-Gaia-10M 17.017M 60.602M 475.702M Supply Chain 2.210M 201.000M 403.000M QAP-WIL-100 1.980M 49.015M 198.020M Heat Source Easy 15.625M 31.628M 125.000M Heat Source Hard 15.625M 31.628M 125.000M
8/1 2/1 4/2 2/1 2/1 1/1 1/1 1/1 2/1
593 s 232 s 133 s 1,921 s 5,839 s 2,305 s 243 s 162 s 15.9 s – –
933 s 373 s 160 s 2,000 s 5,890 s 2,358 s 335 s 224 s 44.6 s
21.06 h 9.34 h – 53.81 h 32.62 h – 2.99 h 18.90 h 0.28 h 59.97 h –
Gurobi – 32.4 h – – – 5.8 h 28.1 h 2.5 h – – –
Table 1: Google PDLP benchmark results [5]. S HARD LP runtimes are reported in seconds; published CPU runtimes are reported in hours. A dash marks a target not reached. Instances are ordered by decreasing nnz(A). S HARD LP reaches the published PDLP criterion on nine of the eleven benchmark instances, versus eight in the published CPU PDLP study (Table 1). On TSP-Gaia-100M, eight H200 GPUs require 592.6 seconds of solver time and 933.4 seconds end to end; the published 32-core CPU experiment reports 21.06 hours. 3.2. RQ2: Scaling beyond the Google benchmark Table 2 pushes the same distributed solve beyond the Google benchmark. It includes KDD12 [13, 20], QAPLIB [1, 11] and MS1 matching [18] cases, plus a deterministic multicommodity-flow family used for the largest scale tests. MCF5B is validated in all three runs on 40 GPUs across ten compute nodes. Larger instances extend the results to 25.000 billion nonzeros on 56 GPUs and 40.807 billion nonzeros on 76 GPUs across 29 nodes. The latter has 13.604 billion variables and solves in 1,905 seconds. The KDD12, QAPLIB and MS1 cases show that the distributed path is not limited to the synthetic flow family. Within the comparison scope in Appendix E, MCF13.60B is, to our knowledge, the largest reported GPU LP solve with a separately validated complete primal–dual solution.
4
D ISTRIBUTED LP ON GPU C LUSTERS
S HARD LP – SECONDS
Problem size Instance
source
MCF13.60B MCF8.35B MCF8B MCF5B tai256c AJ KDD12 L1-SVM MS1 matching
synthetic synthetic synthetic synthetic QAPLIB public data public data
m
n
nnz(A)
GPUs/nodes
solver
end-to-end
2.040M 24.951M 101.059M 63.162M 33.424M 149.639M 43.144M
13.604B 8.353B 8.000B 5.000B 2.131B 259.012M 1.309B
40.807B 25.000B 16.042B 10.027B 8.557B 3.442B 2.617B
76 mixed / 29 56 mixed / 25 12 H200 / 3 40 RTX / 10 8 BW / 1 8 H200 / 1 32 V100 / 16
1,905 s 2,933 s 14,905 s 560 s 717 s 37,294 s 3,721 s
2,057 s 3,045 s 15,510 s 634 s 920 s 37,587 s 3,760 s
Table 2: Validated scale experiments beyond the Google benchmark. All rows are separately validated solves and are ordered by decreasing nnz(A). BW and RTX denote RTX PRO 6000 Blackwell and Quadro RTX 6000, respectively. Mixed GPU allocations are detailed in Appendix A. 3.3. RQ3: Does communication limit multi-GPU scaling? Adding GPUs reduces local matrix work but adds communication. On Mediterranean Shipping, 4,000 iterations take 107.2 seconds on one RTX PRO 6000 Blackwell GPU and 112.7 seconds on four GPUs: the solver loop is slightly slower, even though end-to-end time improves from 281.5 to 168.1 seconds (Figure 2a). Thus, smaller per-GPU matrix blocks do not automatically make the solver loop faster. With a column partition, each GPU stores only some columns of A. A dense reduction still involves every GPU when forming Ax, even when a GPU has no nonzeros in a row. Support-aware communication skips those GPUs for that row. Because their contribution is zero, the exact-arithmetic PDHG update is unchanged; Appendix D gives the details.
Speedup vs. 1 GPU
(a) Fixed-work scaling Mediterranean Shipping end-to-end 1.67×
(b) Benefit of support-aware communication 1.0× = no improvement
Design Match
1.52×
1.6×
Quadro RTX 6000 · 8 GPUs / 4 nodes
1.4×
Design Match V100 · 8 GPUs / 4 nodes
1.2×
solver loop 0.95×
1.0× 1
2
93.0% modeled reduction
1.27× 93.0% modeled reduction
QAP-WIL
1.02×
V100 · 2 GPUs / 1 node
38.4% modeled reduction
4
1.0× 1.2× 1.4× 1.6×
GPUs on one compute node Runtime speedup: dense / support-aware
Figure 2: Communication effects. (a) Speedup for a fixed 4,000 Mediterranean Shipping iterations relative to one GPU. (b) Geometric-mean speedup of support-aware over dense communication; 1× means no change. Figure 2b shows that the benefit depends on how many GPUs actually contain each row. On Design Match, modeled communication falls by 92.97%, and support-aware communication is 5
D ISTRIBUTED LP ON GPU C LUSTERS
1.269× faster on V100s and 1.518× faster on Quadro RTX 6000s. On QAP-WIL, the modeled reduction is only 38.35% and the speedup is 1.024×. Every paired run passes the separate originalspace 10−6 check. Pairing and timing details are in Appendix A.
4. Limitations and Future Work Limitations. S HARD LP assumes that the LP has already been partitioned into shards; converting a monolithic model to this format is currently an offline preprocessing step. In addition, support-aware communication is implemented only for column partitions with fixed matrix support. These choices are sufficient for the experiments in this paper, but they leave room for a more general distributed solver pipeline. Future work. Future work will focus on scaling to still larger LPs and reducing both communication and memory use. This includes better partitioning, broader distributed presolve and postsolve, and extending support-aware communication to more general process grids. We are also interested in applying the same ideas to other large-scale optimization problems, including quadratic and semidefinite programs, and in using S HARD LP as a subroutine inside algorithms for very large discrete optimization problems, such as decomposition methods or large-scale heuristics.
5. Conclusion Previous work distributes the main PDHG computations across GPUs. S HARD LP goes further by keeping the problem and solution distributed throughout the solve. Persistent ownership enables solutions meeting the published criterion on nine of eleven Google benchmarks and separately checked LPs with 13.604 billion variables and 40.807 billion nonzeros, with validated runs spanning up to 29 compute nodes. Smaller per-GPU blocks alone do not ensure speed: communication must also follow matrix sparsity. On Design Match, support-aware communication cuts modeled data movement by 92.97% and yields 1.27×–1.52× solver-time speedups. Together, these results isolate two distinct scale bottlenecks: memory ownership and communication.
6
D ISTRIBUTED LP ON GPU C LUSTERS
References [1] Warren P. Adams and Terri A. Johnson. Improved linear programming-based lower bounds for the quadratic assignment problem. In Panos M. Pardalos and Henry Wolkowicz, editors, Quadratic Assignment and Related Problems, volume 16 of DIMACS Series on Discrete Mathematics and Theoretical Computer Science, pages 43–75. American Mathematical Society, Providence, RI, 1994. [2] David Applegate, Robert Bixby, Vašek Chvátal, and William Cook. Concorde TSP solver, 2020. URL https://www.math.uwaterloo.ca/tsp/concorde.html. [3] David Applegate, Mateo Dı́az, Oliver Hinder, Haihao Lu, Miles Lubin, Brendan O’Donoghue, and Warren Schudy. Practical large-scale linear programming using primal-dual hybrid gradient. Advances in Neural Information Processing Systems, 34:20243–20257, 2021. [4] David Applegate, Oliver Hinder, Haihao Lu, and Miles Lubin. Faster first-order primal-dual methods for linear programming using restarts and sharpness. Mathematical Programming, 201(1):133–184, 2023. doi: 10.1007/s10107-022-01901-9. [5] David Applegate, Mateo Dı́az, Oliver Hinder, Haihao Lu, Miles Lubin, Brendan O’Donoghue, and Warren Schudy. PDLP: A practical first-order method for large-scale linear programming. Mathematical Programming Computation, 2026. doi: 10.1007/s12532-026-00309-2. Published online. [6] Kinjal Basu, Amol Ghoting, Rahul Mazumder, and Yao Pan. ECLIPSE: An extreme-scale linear program solver for web-applications. In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 704–714, 2020. [7] Aharon Ben-Tal, Alexander Goryashko, Elana Guslitzer, and Arkadi Nemirovski. Adjustable robust solutions of uncertain linear programs. Mathematical Programming, 99(2):351–376, 2004. doi: 10.1007/s10107-003-0454-y. [8] Amanda Bienz, William Gropp, and Luke Olson. Node-aware sparse matrix–vector multiplication. Journal of Parallel and Distributed Computing, 130:166–178, 2019. doi: 10.1016/j.jpdc.2019.03.016. [9] Berit D. Brouer, J. Fernando Alvarez, Christian E. M. Plum, David Pisinger, and Mikkel M. Sigurd. A base integer programming model and benchmark suite for liner-shipping network design. Transportation Science, 48(2):281–312, 2014. doi: 10.1287/trsc.2013.0471. [10] A. G. A. Brown, A. Vallenari, T. Prusti, et al. Gaia data release 2: Summary of the contents and survey properties. Astronomy & Astrophysics, 616:A1, 2018. doi: 10.1051/0004-6361/ 201833051. [11] Rainer E. Burkard, Stefan E. Karisch, and Franz Rendl. QAPLIB—a quadratic assignment problem library. Journal of Global Optimization, 10:391–403, 1997. doi: 10.1023/A: 1008293323270.
7
D ISTRIBUTED LP ON GPU C LUSTERS
[12] Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision, 40:120–145, 2011. doi: 10.1007/s10851-010-0251-1. [13] Chih-Chung Chang and Chih-Jen Lin. LIBSVM: A library for support vector machines. ACM Transactions on Intelligent Systems and Technology, 2(3):27:1–27:27, 2011. doi: 10.1145/ 1961189.1961199. [14] William Cook. 100,000,000 stars, 2024. URL https://www.math.uwaterloo.ca/ tsp/star/star100m.html. [15] William Cook. 10,000,000 stars, 2024. URL https://www.math.uwaterloo.ca/ tsp/star/star10m.html. [16] Torsten Hoefler and Jesper Larsson Träff. Sparse collective operations for MPI. In 2009 IEEE International Symposium on Parallel and Distributed Processing, pages 1–8, 2009. doi: 10.1109/IPDPS.2009.5160935. [17] Lannie Dalton Hough, Emir Gencer, Hoffmann Muki, and Abhinav Bhatele. Adaptive spaceefficient collectives for dynamic and unstructured sparsity on GPU platforms. arXiv:2607.04676, 2026. [18] Mohsen Koohi Esfahani, Sebastiano Vigna, Paolo Boldi, Hans Vandierendonck, and Peter Kilpatrick. MS-BioGraphs: Trillion-scale sequence similarity graph datasets. IEEE DataPort, 2024. URL https://doi.org/10.21227/gmd9-1534. [19] Hongpei Li, Yicheng Huang, Huikang Liu, Dongdong Ge, and Yinyu Ye. D-PDLP: Scaling PDLP to distributed multi-GPU systems. arXiv:2601.07628, 2026. [20] LIBSVM Project. LIBSVM data: Classification (binary class), 2026. URL https: //www.csie.ntu.edu.tw/˜cjlin/libsvmtools/datasets/binary.html. KDD2012 data description; accessed 2026-09-01. [21] Haihao Lu, Zedong Peng, and Jinwen Yang. MPAX: Mathematical programming in JAX. arXiv:2412.09734, 2024. [22] Haihao Lu, Zedong Peng, and Jinwen Yang. cuPDLPx: A further enhanced GPU-based first-order solver for linear programming. arXiv:2507.14051, 2025. [23] Aida Rahmattalabi, Gregory Dexter, Sanjana Garg, Qinquan Song, Shenyinying Tu, Yuan Gao, Zhipeng Wang, and Rahul Mazumder. Large-scale regularized matching on GPU clusters. arXiv:2606.07777, 2026. [24] Bora Uçar and Cevdet Aykanat. Revisiting hypergraph models for sparse matrix partitioning. SIAM Review, 49(4):595–603, 2007. doi: 10.1137/060662459. [25] Huasha Zhao and John Canny. Sparse AllReduce: Efficient scalable communication for power-law data. arXiv:1312.3020, 2013.
8
D ISTRIBUTED LP ON GPU C LUSTERS
[26] José R. Zubizarreta. Using mixed integer programming for matching in an observational study of kidney failure after surgery. Journal of the American Statistical Association, 107(500): 1360–1371, 2012. doi: 10.1080/01621459.2012.703874.
Appendix A. Computational environment and hardware All experiments were executed on Georgia Tech’s PACE Phoenix cluster under Slurm, with one MPI process per GPU. Table 3 summarizes the node classes used, using retained GPU inventories and the corresponding Slurm node configurations. All solved rows in Table 1 use H200 GPUs; Table 2 identifies the GPU families for the scale experiments. The GPU/node counts in these tables describe GPUs used by the solver, which may be fewer than the GPUs installed in an allocated node. Table 3: PACE Phoenix node classes used in the reported experiments. GPU memory is the capacity reported by nvidia-smi; host RAM is Slurm’s configured node memory, rounded to the nearest GiB. CPU counts are configured cores per node. These are node capacities, not per-job reservations. GPU
memory/GPU installed (GiB) GPUs/node
NVIDIA H200
140.4
RTX PRO 6000 Blackwell Server Edition Quadro RTX 6000 Tesla V100 PCIe
95.6
Host CPU (cores/node)
8
Intel Xeon Platinum 8562Y+ (64) 8 Intel Xeon, Granite Rapids (64)
24.0 16.0 / 32.0
4 Intel Xeon Gold 6226 (24) 2 Intel Xeon Gold 6226 (24)
host RAM (GiB) 2,015 2,015 376 / 754 376 / 754
Software and communication. The distributed C/CUDA runs use CUDA 12.9.1 and CUDAaware OpenMPI 4.1.8. GPU communication uses NVIDIA Collective Communications Library (NCCL) 2.26.5, with builds targeting each GPU architecture. MPI manages the process grid and control reductions; NCCL exchanges GPU buffers. Retained nvidia-smi topo -m output shows NVLink connectivity between the H200 GPUs (NV18), and PCIe paths on the Blackwell, Quadro RTX 6000, and V100 nodes used here. Inter-node transport is run-specific. MCF8B, MCF5B, MS1, and the Design Match communication comparisons disable NCCL’s InfiniBand transport (NCCL IB DISABLE=1) and use its socket path. The two-node QAP-THO-150 benchmark run enables the InfiniBand transport. These settings do not establish a common physical-link bandwidth or a cluster-wide RDMA performance claim. Allocation and timing conditions. MCF8.35B combines 42 Quadro RTX 6000 and 14 V10032GB GPUs across 25 nodes. MCF13.60B combines 60 Quadro RTX 6000, 14 V100-32GB, and two L40S GPUs across 29 nodes. Both use column partitions, one MPI rank per GPU, and socket transport. Jobs reserve GPUs, CPU cores, and host memory through Slurm; whole-node exclusivity is not assumed. For example, MCF8B uses four of the eight H200s on each of three nodes, and the Quadro RTX 6000 Design Match comparisons use two of four GPUs per node. The Mediterranean 1/2/4-GPU fixed-work measurements use exclusive Blackwell nodes, with all three GPU counts tested within each allocation. Dense/support-aware pairs run sequentially on the same allocation and GPU 9
D ISTRIBUTED LP ON GPU C LUSTERS
set, with order reversed across pairs. Figure 2b reports geometric means over two order-balanced pairs per Design Match platform and eight QAP-WIL pairs across four allocations. For Design Match, the timing metric is solver time through the first validated checkpoint; for QAP-WIL, it is elapsed time through validation. MCF4.15B uses two shared nodes; its timing is descriptive rather than a scaling measurement. Table 4: Physical footprint of selected checked executions. Sizes are GiB (230 bytes); input and primal–dual output are totals across all shards. Peak GPU memory is the largest sampled device-memory use on any participating GPU. BW denotes RTX PRO 6000 Blackwell Server Edition. Instance
GPU
TSP-Gaia-100M H200 tai256c AJ BW MCF8B H200
process grid
GPUs/nodes
1×8 1×8 1 × 12
input
8/1 121.59 8/1 148.23 12/3 380.69
peak GPU output memory/GPU 18.87 32.00 119.96
80.95 69.39 79.78
Table 4 reports physical storage and device-memory use. Peak memory is the maximum of periodic nvidia-smi memory.used samples, matched to the participating GPU UUIDs across the allocation. It is a sampled device-level quantity, not an exact allocator high-water mark.
Appendix B. Instance origins Here m, n, and nnz(A) denote constraints, variables, and explicitly stored matrix coefficients. Table 5 summarizes the released dimensions and construction of the Google PDLP benchmark [5].1 We use synthetic for generated, application-inspired models and identify LPs constructed from a named public data source separately. MCF8.35B uses 8,303 commodities, 1,000 factories per commodity, 1,000 warehouses, and five stores. MCF13.60B uses 34 commodities, 20,000 factories per commodity, 20,000 warehouses, and five stores. The latter has more nonzeros but fewer rows, reducing row-vector communication in the column-partitioned solve.
1. The scaled instance files are available from Oliver Hinder’s benchmark download page; generators and detailed construction notes for the subset with released construction code are available in the companion GitHub repository.
10
D ISTRIBUTED LP ON GPU C LUSTERS
Table 5: Meaning and construction of the eleven Google benchmark instances. Instance Design Match
m
n
nnz(A) origin and model
22,000,135
40,000,000 2,760,000,000 Synthetic covariate-balancing statistical-matching LP from the application class studied by Zubizarreta [26]. Heat Source Easy 15,625,000 31,628,008 125,000,000 Synthetic inverse heat-source problem with temperature observations and linear PDE constraints. Heat Source Hard 15,625,000 31,628,008 125,000,000 The same synthetic inverse problem with fewer measurements and more true and candidate source locations. Mediterranean Shipping 7,490,593 208,479,461 628,927,462 LP relaxation of a liner-shipping mixed-integer program (MIP) built from the LINERLIB Mediterranean benchmark [9]. Production Inventory 4,650,850 18,270,600 500,049,700 Synthetic robust production–inventory model based on Ben-Tal et al. [7], with randomized data. QAP–THO–150 6,705,300 249,783,750 1,005,795,000 Adams–Johnson LP relaxation [1] of the QAPLIB tho150 benchmark [11]. QAP–WIL–100 1,980,200 49,015,000 198,020,000 Adams–Johnson LP relaxation [1] of the QAPLIB wil100 benchmark [11]. Supply Chain 2,210,100 201,000,100 403,000,100 Synthetic multicommodity-flow model for large-retailer supply-chain planning. TSP–Gaia–100M 162,934,799 1,184,557,727 6,337,834,450 Public-data-derived TSP lower-bound LP on the 100 million nearest stars in Gaia DR2 [10, 14]; degree constraints plus 62,934,799 cuts collected with Concorde [2]. TSP–Gaia–10M 17,016,681 60,601,996 475,701,996 Public-data-derived counterpart on the ten million nearest stars in Gaia DR2 [10, 15]; degree constraints plus 7,016,681 cuts collected with Concorde [2]. World Shipping 15,304,282 228,867,510 688,658,522 LP relaxation of a liner-shipping MIP built from the LINERLIB World benchmark [9].
11
D ISTRIBUTED LP ON GPU C LUSTERS
Table 6: Origins and validation status of the beyond-benchmark experiments. Instance KDD12 L1-SVM
m
n
149,639,105
259,012,009
nnz(A) origin, model, and status
3,441,699,415 Public-data-derived LP from the LIBSVM kdd12.xz click-log data, using a no-intercept, C = 1, split-variable ℓ1 -regularized hinge-loss model [13, 20]; one validated run. MS1 matching 43,144,218 1,308,742,322 2,617,484,644 Maximum-weight fractional matching LP from the MS-BioGraphs MS1 protein sequence-similarity graph [18]; one variable per undirected non-loop edge, unit vertex capacities, and 0 ≤ xe ≤ 1. Self-loops and reverse duplicates are removed; one original-LP checked run. tai256c AJ 33,423,872 2,130,804,736 8,556,511,232 Benchmark-derived, symmetry-reduced Adams–Johnson/substitution LP generated from QAPLIB tai256c [1, 11]; one separately checked execution. MCF4.15B 52,423,770 4,149,986,400 8,321,930,400 Deterministic synthetic multicommodity-flow model; one separately checked execution on eight RTX PRO 6000 Blackwell GPUs across two shared compute nodes, used only to demonstrate feasibility, not performance or scaling. MCF5B 63,161,790 5,000,032,800 10,026,520,800 Deterministic synthetic multicommodity-flow model with commodity-local factory–warehouse–store blocks; three separately checked executions. MCF8B 101,059,055 8,000,067,600 16,042,463,600 Larger member of the same deterministic synthetic family; two separately checked executions. MCF8.35B 24,950,515 8,352,818,000 25,000,333,000 Same deterministic MCF family; one checked execution on 56 GPUs across 25 nodes. MCF13.60B 2,040,170 13,604,080,000 40,807,480,000 Same family with larger commodity blocks; one checked execution on 76 GPUs across 29 nodes. MCF8.5B capacity 107,374,470 8,500,010,400 17,044,994,400 Capacity-only run from the same family on 70 GPUs; complete solution vectors were not exported, so the run is not counted as a solve.
12
D ISTRIBUTED LP ON GPU C LUSTERS
Appendix C. Distributed solver stages and solution validation C.1. Data layout across solver stages The solve begins from prepared matrix shards. Conversion from a monolithic input (e.g., an MPS file) to these shards is performed offline and is not included in reported times. In an R × C grid, each process holds Arc , xc , and yr , with vector replicas only within the corresponding grid column or row. Table 7 specifies how this ownership is retained. Table 7: Implemented distributed lifecycle. “Global reduction” means a scalar or partition-sized collective, never a complete matrix or primal–dual gather. Phase
Ownership and communication
Prepared-shard input
Each process reads its matrix block and matching row/column data directly. Offline construction and conversion are excluded; no claim is made for arbitrary monolithic-input conversion. Scaling Local nonzeros form row and column statistics; grid reductions form the diagonal scales, which remain partitioned with the vectors. Optimization Local products plus the two reductions in Equation (5); Halpern, projection, and vector updates act on local blocks. Restart/termination Processes reduce norms and scalar decision statistics; latest, average, anchor, and reflected states retain the same block layout. Singleton-zero One exact elimination pass identifies zero-equality singleton rows, fixes their columns presolve at zero, removes those entries locally, and stores a distributed pivot map. This is not arbitrary general presolve. Feasibility Primal and dual phases reuse the distributed operator; candidate selection and restart polishing decisions use collective statistics while phase states remain local. Recovery, output, The pivot map reconstructs eliminated solution coordinates on their replicas; each process and check writes owned primal, dual, and reduced-cost slices. The separate checker consumes those shards without a root-vector gather.
For a singleton pivot (j, s, asj ), recovery forms the original-unit reduced cost r = c − A⊤ y and applies xj = 0, ys ← ys + rj /asj , and rj = 0. This recovery covers primal–dual solutions; infeasibility and unboundedness rays are not implemented. Prepared blocks use contiguous row and column intervals in a fixed coordinate order; the partition does not change during a solve. C.2. Validation criteria and timing For Table 1, general presolve and the singleton-zero reduction are disabled, and the scaling from the Google PDLP study is used. Both arms of the later Design Match communication comparisons use the same singleton-zero reduction. The published rule requires absolute ℓ∞ primal and stationarity residuals at most 10−8 on the scaled LP, primal/dual bound and sign membership, and a relative primal–dual gap at most 10−2 . Our finite-precision membership tests use a 10−8 tolerance [5]. For every other run counted as solved, we use a uniform acceptance tolerance of 10−6 . The separate checker reads the postsolved, sharded solution in the coordinates of the original unscaled LP. Its inputs are the primal vector x, the row-dual vector y, and the exported reduced-cost vector r. It tests four requirements: primal feasibility, stationarity, admissible signs for y and r, and agreement between the primal and dual objectives.
13
D ISTRIBUTED LP ON GPU C LUSTERS
For an interval I = [ℓ, u], define
Interval quantities.
s(I) := max {0} ∪ {|b| : b ∈ {ℓ, u}, b finite} ,
vI (t) := dist(t, I),
D(I) := {q ∈ R : q > 0 ⇒ ℓ > −∞,
q < 0 ⇒ u < +∞},
ψI (q) := inf qt. t∈I
Here, vI (t) is the violation of the interval, and s(I) is a scale computed from its finite endpoints. The set D(I) gives the multiplier signs allowed by the interval: a positive multiplier requires a finite lower bound, and a negative multiplier requires a finite upper bound. The function ψI (q) is the contribution of that interval to the dual objective. Stationarity and objective values.
We form sign-admissible copies of the reported dual quantities,
ȳi = projD(Si ) (yi ),
r̄j = projD(Xj ) (rj ).
These projections are used only when evaluating stationarity and the dual objective. They do not hide invalid reported signs: the distances between (y, r) and (ȳ, r̄) are checked separately below. Using the projected quantities, define the stationarity error e = c − A⊤ ȳ − r̄, the primal objective p = c0 + c⊤ x, and the dual objective d = c0 +
X
ψSi (ȳi ) +
i
X
ψXj (r̄j ),
j
where c0 = 0 if the model has no objective constant. For each variable and row, define the primal violations and their data scales: δjx = vXj (xj ),
bxj = s(Xj ),
δic = vSi ((Ax)i ),
bci = s(Si ).
Thus, δ x measures violations of the variable bounds x ∈ X, while δ c measures violations of the row bounds Ax ∈ S. Nine acceptance criteria.
The checker evaluates
∥δ x ∥∞ , 1 + ∥bx ∥∞
g2 = max
g1 =
δic
j
δjx , 1 + bxj
g3 =
∥δ c ∥2 , 1 + ∥bc ∥2
|ej | , 1 + |cj | |p − d| g7 = max dist(yi , D(Si )), g8 = max dist(rj , D(Xj )), g9 = . i j 1 + |p| + |d|
g4 = max i
, 1 + bci
g5 =
∥e∥2 , 1 + ∥c∥2
g6 = max
(6)
j
The pairs (g1 , g2 ), (g3 , g4 ), and (g5 , g6 ) measure variable-bound feasibility, row feasibility, and stationarity, respectively. In each pair, the first quantity is a normwise residual and the second is the worst coordinatewise normalized residual. Criteria g7 and g8 measure violations of the required dual and reduced-cost signs, and g9 is the relative primal–dual objective gap. The added 1 in each denominator keeps the normalization well defined when the corresponding data or objective is zero. 14
D ISTRIBUTED LP ON GPU C LUSTERS
A candidate is accepted if and only if every computed quantity is finite and max gq ≤ 10−6 .
1≤q≤9
Table 8 displays g3 , g5 , and g9 in the “rel. primal,” “rel. stationarity,” and “rel. gap” columns. Its “max. criterion” column is the maximum over all nine tests, including the six not printed separately. The separate checker evaluates the tests in original coordinates using long-double accumulation. Table 8: Iterations, runtimes, and separate original-space checks for the scale experiments. Convergence is evaluated every 200 iterations, giving iteration counts in multiples of 200. Times are solver/end-to-end. The maximum criterion is the largest of the nine quantities in Equation (6). Instance
iter.
solver / end-to-end (s)
rel. primal
rel. stationarity
rel. gap
max. criterion
tai256c AJ 3,200 717 / 920 3.62 × 10−10 8.87 × 10−17 5.82 × 10−10 1.26 × 10−9 KDD12 L1-SVM 74,600 37,294 / 37,587 1.82 × 10−10 1.23 × 10−10 4.57 × 10−10 6.72 × 10−7 MS1 matching 31,800 3,721 / 3,760 4.10 × 10−11 1.74 × 10−16 5.27 × 10−15 9.10 × 10−8 MCF4.15B 1,200 490 / 812 1.01 × 10−9 3.43 × 10−16 7.25 × 10−9 2.86 × 10−8 MCF5B 1,200 560 / 634 1.01 × 10−9 8.21 × 10−15 7.26 × 10−9 2.86 × 10−8 MCF8B 1,200 14,905 / 15,510 1.01 × 10−9 2.92 × 10−16 7.27 × 10−9 2.86 × 10−8 MCF8.35B 4,800 2,933 / 3,045 3.15 × 10−10 1.25 × 10−13 2.19 × 10−9 2.19 × 10−9 MCF13.60B 4,600 1,905 / 2,057 1.04 × 10−10 1.32 × 10−13 6.74 × 10−8 6.74 × 10−8
Solver time is measured inside the iterative solver and includes residual evaluations, restart, and feasibility polishing when enabled. End-to-end time adds rank-local input, initialization, solution recovery, and output, but excludes offline conversion and the separate final check. In Table 1, TSP-Gaia-10M, Production Inventory, QAP-WIL-100, and Supply Chain are medians of three runs; the remaining solved rows use one checked run. Table 2 uses one checked run each for KDD12 L1-SVM, tai256c AJ, and MS1 matching, the componentwise median of three MCF5B runs, and the first of two checked MCF8B runs, whose end-to-end times were 15,509.9 and 15,519.7 seconds. Table 8 additionally includes one checked MCF4.15B run; it used two shared compute nodes, so its timing is descriptive rather than a scaling measurement. MCF8.35B and MCF13.60B each report one cold-start execution without presolve or feasibility polishing. Internal tolerances are 10−8 and 10−7 , respectively; both pass the same independent 10−6 check. The Google study used a six-day limit [5]; accordingly, Table 1 compares coverage and reports raw runtimes rather than a controlled hardware comparison. The published CPU experiments used either a 16-core AMD EPYC 7302 compute node with 256 GB of memory or a 32-core Intel Xeon Platinum 8352Y compute node with 1 TB. PDLP used 16 or 32 threads, while Gurobi 11.0.2 used its default thread count and termination rule. Table 1 reports the fastest published Gurobi time among barrier, primal simplex, and dual simplex [5]. MS1 uses a 1 × 32 layout over sixteen nodes with two V100 GPUs each (16- and 32-GiB variants), socket transport, a cold start, and no presolve or feasibility polishing. The exported solution passes both the original-LP checker and an independent original-source audit at 10−6 . Its end-to-end time is the solver-process wall time, excluding the separate final check. MCF8B uses a 1 × 12 layout over three compute nodes, without presolve or feasibility polishing. Hardware, software versions, interconnect settings, and allocation details for all experiments are reported in Appendix A.
15
D ISTRIBUTED LP ON GPU C LUSTERS
Appendix D. Support-aware update equivalence and communication model For a 1 × p column partition, define the participant set of row i by Pi := j ∈ {0, . . . , p − 1} : nnz(Ai,Jj ) > 0 ,
ki := |Pi |.
(7)
Thus, Pi contains exactly the GPUs whose column shards store a nonzero coefficient in row i. (a) Row incidence
(b) Dense update
column shard j=0 j=1 j=2 j=3
ai
Ai,J0
0
Ai,J2
0
(0)
0
ai
(2)
0
ai
Pi = {0, 2}
(0)
(c) Participant update (2)
0
ai
0
(0)
ai
P (j) si = 3j=0 ai + each rank computes yi = Φi (yi , si ) yi+
yi+
yi+
yi+
communicate and store on all four ranks
(2)
–
–
ai (0)
(2)
owner: si = ai + ai yi+ = Φi (yi , si ) yi+
–
yi+
–
communicate only on Pi
Figure 3: Dense and support-aware updates for a row supported on column shards 0 and 2. The dense path reduces the row activity and stores the updated dual coordinate on every rank. The participant path reduces at an owner and returns the update only to ranks whose local A⊤ y product uses that coordinate. The reported participant runs use a 1 × p grid and apply the construction to the Ax reduction and dual update; A⊤ y follows its unchanged path. The communication plan is built once after all supportchanging preprocessing, from scaled local compressed sparse row (CSR) support, and checked for agreement across the row communicator. Rows with the same participant set are packed together. Full-participation buckets use NCCL Reduce/Broadcast, whereas partial buckets use grouped NCCL point-to-point transfers through the lowest participating rank; rank zero owns structurally empty rows. The sparsity pattern remains fixed during a solve, so a structural change or repartitioning would require rebuilding the plan. Extending this path to a general R × C grid remains future work. Exact-arithmetic equivalence. Assume exact arithmetic and that both executions use the same arithmetic, parameters, and restart, evaluation, checkpoint, and polishing rules. At the start of each participant epoch, the dense and participant executions must have the same fully materialized solver state—including current, average, anchor, and reflected primal–dual components. The support of A must then remain fixed, every Pi in Equation (7) must be the exact row support over ranks, and each nonempty set must have an owner oi ∈ Pi holding the authoritative yi . For j ∈ / Pi , Ai,Jj = 0. Let x e denote the vector supplied to the Ax product (for Equation (3), x e = 2b x − x), and define P P (j) (j) (j) (j) ai = Ai,Jj x eJj ; then ai = 0 whenever j ∈ / Pi . Thus j ai = j∈Pi ai , which gives the same scalar row activity to the ith coordinate of ΦS,σ . The owner computes the same coordinate update, and dissemination reaches every rank whose transpose product contains that coordinate. The local Halpern and vector updates therefore match. Before evaluation, adaptive restart, checkpoint selection, polishing, or any other operation requiring replicated dual state, the participant execution reconstructs every required coordinate; global statistics count each authoritative coordinate once. Induction over these epochs then proves whole-solver equality at synchronization points. If Pi = ∅, activity is identically zero and a designated owner evaluates the ith dual update with zero row activity and no operator communication. 16
D ISTRIBUTED LP ON GPU C LUSTERS
Floating-point reductions can use different summation orders. Across the shared evaluation checkpoints of the V100 Design Match paired runs, the largest absolute differences between the reported dense and participant fields are 1.87 × 10−12 for relative primal residual, 1.98 × 10−11 for relative stationarity, and 1.30 × 10−15 for relative gap. For p > 1, count one scalar transferred across one rank-to-rank hop as one logical scalar-hop. Let S be the ms rows assigned to the participant path, h the number of dual updates, b the number of full-state synchronization boundaries, and q the number of dual vectors synchronized at each boundary. Under the ring model, the corresponding counts are Hdense = 2h(p−1)ms ,
Hpart = 2h
X (ki −1)+ +qb(p−1)ms , i∈S
ρtot = 1−
Hpart , (8) Hdense
where (t)+ = max{t, 0}. The factor of two counts the reduction and dissemination legs of each dual update; the boundary term counts a one-way full-state synchronization. Figure 2b reports reductions in this logical scalar-hop count, not measured bytes on physical links. Packing, imbalance, topology, and latency are excluded, so the model need not predict elapsed time. For Design Match, P p = 8, ms = 22,000,135, h = 27,081, b = 138, q = 1, and i (ki − 1)+ = 10,426,445, giving the reported 92.9748%. Hypergraph partitioning, configured sparse collectives, and node-aware sparse matrix–vector multiplication provide the related sparse-communication work [8, 16, 17, 24, 25]. We do not claim a new generic collective: the contribution is the solver-specific integration of a fixed producer/consumer plan with the nonlinear ΦS,σ update, empty-row ownership, reflected and restart state, polishing, checking, and output boundaries.
Appendix E. Related work and comparison scope For this comparison, we focus on GPU LP solvers that use a general sharded sparse-matrix representation and solver kernels, rather than application-specific operators or decompositions. We count a run only if it exports a complete primal–dual solution that passes the independent checker. D-PDLP evaluates its two-dimensional decomposition on a single compute node containing eight H100 GPUs connected by NVLink. Its largest-variable instance has 127.5 million variables and 257.5 million nonzeros; another reaches 690 million nonzeros [19]. MPAX reports dense multi-GPU PDHG at approximately 900 million represented nonzeros on up to four H100 GPUs, with efficient distributed sparse-data sharding left as future work [21]. NVIDIA’s cuOpt 26.08 release notes report 2.5×–8.8× PDLP speedups on eight NVLink-connected B200 GPUs; they do not identify a cross-node execution.2 The cuOpt FAQ gives single-H100 capacity examples up to two billion nonzeros.3 A specialized regularized-matching solver spans sixteen H100 GPUs across two compute nodes and reports models up to one billion nonzeros [23]. It uses a source-decomposable formulation and a dual-gradient method specialized to matching. ECLIPSE reports a 1012 -variable structured web LP by optimizing the smooth dual of a ridge-perturbed model with distributed matrix–vector products [6]. Among the general sparse GPU LP solvers discussed above, MCF13.60B is, to our knowledge, the largest reported validated solve, with 40.807 billion nonzeros. This is a scale claim, not a claim 2. NVIDIA cuOpt 26.08 release notes. 3. NVIDIA cuOpt FAQ.
17
D ISTRIBUTED LP ON GPU C LUSTERS
of being the first multi-node GPU optimization method or of a matched speedup over the systems above. Additionally, the 70-GPU capacity experiment fits a deterministic flow model with 8.500 billion variables, 107.374 million constraints, and 17.045 billion nonzeros across twenty compute nodes. It ran 280 iterations in 217.46 seconds with 339 MiB free on the fullest GPU, but did not export the complete vectors needed for a separate primal–dual check. The time characterizes this capacity execution, not a checked solve or a time-to-solution result.
18