Profiling report: PDLP on the 10 largest MIPLIB2017 benchmark instances (L40S, CuSparseMatrixCSR, Float32)

[gdalle-bot]

Opened by a coding agent at @gdalle's request — please review before acting on it. Everything below is measured, not inferred; the harness, raw result.json files, CUPTI traces, Julia profiles and nvidia-smi samples are on master01-hpc under ~/pdlp-profiling/ (see README.md there for the run order).

What was measured

Hardware / software. One NVIDIA L40S (sm_89, 44.99 GiB, 142 SMs, 864 GB/s theoretical peak, 350 W cap, SM clock pinned at 2520 MHz throughout), driver 590.48.01 / CUDA 13.1, Julia 1.12.7, CUDA.jl CUDACore 6.3.0, cuSPARSE 12.7.3.

Configuration. PDLP(Float32, Int32, cuSPARSE.CuSparseMatrixCSR; backend = CUDABackend(), time_limit = 120.0), all other parameters at their defaults (ruiz_iter = 10, check_every = 100, termination_reltol = 1e-4, record_error_history = true). One fresh Julia process per instance, exclusive access to the GPU.

Instances. The 10 largest by nnz(A) in the MIPLIB2017 benchmark set, found by reading all 240 instances through MathOptBenchmarkInstances.jl. (8 of the 240 could not be read at all — read_miplib2017_instance lowercases the name, so mixed-case files like eilA101-2.mps.gz are not found; all 8 are far too small to belong in the top 10. Worth an issue on that repo.)

Three instruments, as requested.

1. CUDA.@profile (CUPTI): full host-API + device-activity trace over a ~6 s window, with NVTX.@range markers around steps / termination_check! / restart_check! / restart!. Device operations are attributed to a phase through their correlation ID → launching API call → enclosing NVTX range. Note CUPTI inflates absolute host times, so only shares are quoted from it.
2. Julia's sampling profiler over the same window, restricted to samples whose backtrace contains the solver loop.
3. nvidia-smi sampled at 250 ms for the whole 120 s of each clean (unprofiled) run.

Plus per-call microbenchmarks with CUDA.@elapsed (median of 20), which are the numbers quoted in µs.

nsys/ncu are not installed on the cluster, so the CUPTI-based integrated profiler is the timeline source.

Headline

One number dominates everything else in this report: the step size is wrong, and it is wrong because spectral_norm does not return a spectral norm. Kernel-level work is in reasonable shape — cuSPARSE SpMV runs at 63–110 % of the L40S's DRAM bandwidth (above 100 % where the matrix fits in the 48 MB L2) and is 47–97 % of all GPU time — so the throughput wins available from micro-optimisation are worth maybe 10–35 %, whereas the step size is worth up to 45× on the same instance. The rest of the findings are real and worth fixing, but they should be read in that order of magnitude.

Also note these are the ten largest instances, i.e. the most favourable case for the current design. The fixed per-iteration overheads below get proportionally worse on the other 230.

---

1. spectral_norm returns a provably wrong value, and the step size follows it — refines #95

fixed_stepsize sets η = 0.9 / spectral_norm(A, At) and never changes it again (reset_stepsize! only zeroes η_sum). So this single scalar fixes the step size for the whole solve.

spectral_norm calls powm!(KᵀK, x0) with no tol/maxiter, and x0 is a raw randn! vector (norm ≈ √n) rather than the normalised guess powm! documents. IterativeSolvers' default tol = eps(T)·size(B,2)^3 is astronomically loose in Float32 (≈ 1.0e8 for n = 10⁵), so the iteration stops before the first normalisation ever takes effect. powm! then returns θ = x₀ᵀKᵀKx₀ = ‖Kx₀‖², and spectral_norm hands back √θ = ‖Kx₀‖ — a quantity that scales with ‖x₀‖ ≈ √n, rather than the Rayleigh quotient ‖Kx₀‖/‖x₀‖.

That claim is not a guess. ‖A‖₂ ≤ √(‖A‖₁·‖A‖∞) is a rigorous upper bound, and on 3 of the 4 audited instances the returned value exceeds it:

instance	spectral_norm returns	rigorous upper bound √(‖A‖₁‖A‖∞)	300-iteration normalised power method	η used	η from converged ‖A‖
square47	29.659	5.826 ❌	1.0000	0.0303	0.900
roi5alpha10n8	23.170	6.292 ❌	1.0000	0.0388	0.900
supportcase7	29.709	13.901 ❌	0.9999	0.0303	0.900
neos-4647030-tutaki	0.9436	1.916 ✅	0.9968	0.9538	0.903

Across all ten instances the ratio (converged ‖A‖) / (returned value) ranges from 0.0027 (neos-5114902-kasavu, i.e. a 370× overestimate) to 1.53 (neos-3402454-bohle, an underestimate that puts η‖A‖ ≈ 1.38 > 1, violating the PDHG stability condition). The error is unpredictable in both magnitude and sign — see Table 8.

Because Ruiz + Chambolle–Pock equilibration leaves ‖A‖₂ ≈ 1, the practical effect is usually η about 25–30× too small.

What that costs — A/B at identical settings, 120 s budget

instance	η default	result	η corrected	result
roi5alpha10n8	0.0388	OPTIMAL, 382,100 iters, 33.4 s	0.900	OPTIMAL, 9,500 iters, 0.75 s
square47	0.0303	TIME_LIMIT, rel. err 9.19e-01	0.900	TIME_LIMIT, rel. err 1.38e-01
neos-4647030-tutaki	0.9538	OPTIMAL, 21,800 iters, 3.0 s	0.903	OPTIMAL, 32,600 iters, 4.2 s
supportcase7	0.0303	TIME_LIMIT, rel. err 4.68e-03	0.900	TIME_LIMIT, rel. err 1.11e+04 💥

A 45× speedup on roi5alpha10n8 and a 6.7× better residual on square47, from one scalar.

But the fix is not just "compute ‖A‖ correctly"

supportcase7 diverges with the mathematically-correct η = 0.9/‖A‖, while the accidental 30×-too-small η converges nicely. In other words the broken estimator is currently acting as an erratic safety factor that happens to paper over the absence of the adaptive step size that PDLP actually specifies. Fixing spectral_norm on its own would be a regression on at least one of these ten instances.

Suggested fix, in order:

1. Fix spectral_norm: normalise x0, pass an explicit tol and maxiter, check powm!'s convergence flag, guard sqrt against a negative λ, and cross-check against √(‖A‖₁‖A‖∞) (both are cheap SpMVs against abs.(nzval) and give a free correctness assertion).
2. Then add the adaptive/linesearch step size from Applegate et al. §3.1 (what cuPDLP/cuPDLPx use), or at minimum a divergence guard that halves η and restarts when the KKT residual grows by orders of magnitude. Step 1 without step 2 is not safe.
3. Consider making the safety factor invnorm_scaling do real work in Float32 (0.9 leaves little margin once ‖A‖ itself is a Float32 estimate).

---

2. Five KKT evaluations per check where two would do — extends #100

Per outer check (check_every = 100), kkt_errors! runs five times:

caller	evaluates	needed?
termination_check!	sol → stats.err	yes
restart_check! → best_error!!	sol → err_current	no — exact duplicate of the line above (#100)
restart_check! → best_error!!	sol_avg → err_avg	yes
restart_check! → best_error!!	sol_last → err_last	no — see below
restart_check! → best_error!!	sol_avg_last → err_avg_last	no — see below

Verified on the GPU (redundancy_check.jl): termination_check! and restart_check! return bit-identical errors for the same iterate, e.g. on neos-4647030-tutaki both give (19.947240829467773, 0.0008696960867382586, 53675.42578125).

The err_last / err_avg_last pair exists only to produce abs_candidate_last for the no_local_progress = candidate > candidate_last test. But restart_stats.abs_candidate from the previous call already holds a candidate error — caching it would remove two full evaluations. It does change the semantics from "one iteration ago" to "one check ago", which is what the reference implementations actually compare, and is arguably more correct given the criterion is only evaluated every 100 iterations anyway.

Each kkt_errors! is 2 SpMVs + ~9 elementwise passes + 7 reductions. Removing three of the five is worth, per check:

	square47	square41	neos-2075418-temuka	neos-4647030-tutaki	roi5alpha10n8
3 × kkt_errors!	3725 µs	2688 µs	1886 µs	1011 µs	899 µs
share of a 100-iteration check	5.6 %	8.0 %	8.1 %	9.4 %	12.0 %

(Full table in the appendix; the check machinery as a whole is 8.5 %–19.7 % of wall clock, mean 13.0 %.)

2b. primal_scale and dual_scale are problem constants, recomputed every time

rescaled_combined_bounds = @. scratch.y = inv(D1.diag) * combine(lc, uc)
err.primal_scale = colnorm!!(err.primal_scale, rescaled_combined_bounds)
...
rescaled_obj = @. scratch.x = inv(D2.diag) * c
err.dual_scale = colnorm!!(err.dual_scale, rescaled_obj)

These depend only on D1, lc, uc, D2, c — none of which move during a solve. Confirmed bit-identical over 1000 iterations on two instances:

after  200 iterations: primal_scale = 32085.00195   dual_scale = 38.41656876
after 1000 iterations: primal_scale = 32085.00195   dual_scale = 38.41656876

They cost 2 of the ~9 elementwise passes and 2 of the 7 reductions (with their 2 blocking syncs) in every one of the 5 calls — 10 wasted reductions, 10 wasted vector passes and 10 wasted syncs per check. Hoist them into initialize and store them next to the preconditioner.

---

3. 14.8 GPU operations and ~290 CUDA API calls per PDHG iteration

The CUPTI trace gives exactly the same op count per iteration on all ten instances:

per PDHG iteration	count	source
cusparse::csrmv_v3_kernel	2.10	the two SpMVs (the 0.10 is the checks' share)
cusparse::csr_partition_kernel	2.10	cuSPARSE recomputing its row partitioning on every call
[set device memory]	2.11	cuSPARSE zeroing its workspace on every call
gpu_broadcast_kernel_linear	5.46	clamp ×2, 2x-xₖ, colaxpby! ×2
[copy device to device]	2.01	copy!(sol_avg_last, sol_avg) (#111)
nrm2 / mapreduce	0.71	colnorm!! / colsum!!
[copy device to pageable]	0.36	scalar readbacks (see §4)
total	14.85

For comparison, a fused PDHG step needs ~4: SpMV, fused primal+extrapolation+average, SpMV, fused dual+average.

Where that shows up in GPU time (share of GPU-busy in the profiled window):

	square47 (27.3M nnz)	neos-5114902-kasavu (4.2M)	roi5alpha10n8 (2.4M)	supportcase7 (2.8M)
cuSPARSE SpMV	96.7 %	70.2 %	58.5 %	47.2 %
cuSPARSE row partitioning	1.0 %	5.6 %	12.0 %	15.0 %
elementwise broadcasts	1.3 %	14.1 %	16.1 %	21.0 %
D2D copies (sol_avg_last)	0.4 %	5.5 %	4.4 %	5.6 %
cuSPARSE workspace memset	0.3 %	2.9 %	4.2 %	5.2 %
reductions	0.3 %	1.6 %	4.0 %	5.1 %

And GPU-busy / wall for the whole profiled region drops from 93 % (square47) to 32 % (supportcase7) — the smaller instances are launch-bound, not compute-bound.

3a. cuSPARSE mv! rebuilds everything on every call

cuSPARSE.jl's mv! constructs a fresh CuSparseMatrixDescriptor and two dense-vector descriptors, queries cusparseSpMV_bufferSize, and calls with_workspace — which allocates and frees a CuVector{UInt8} per call (cuMemAllocFromPoolAsync 2.6/iter, cuMemFreeAsync 2.1/iter in the trace). It never calls cusparseSpMV_preprocess, even though mm! already calls cusparseSpMM_preprocess.

Since A and Aᵀ are fixed for the whole solve, CoolPDLP could build the descriptors once, size and keep the workspace once, call cusparseSpMV_preprocess once per matrix, and then issue only cusparseSpMV per iteration. That removes 2.10 partition kernels + 2.11 memsets + ~2 device alloc/free pairs per iteration — worth ~20 % of GPU time on supportcase7/roi5alpha10n8 and ~1 % on square47.

3b. copy!(sol_avg_last, sol_avg) every iteration — confirms #111

Exactly 2.01 device-to-device copies per iteration, 0.4 %–5.6 % of GPU time, for a value read once every 100 iterations. Since update_average! already writes sol_avg in place, the cheapest fix is to write the new average into sol_avg_last and swap the two references (as sol/sol_last already do) — that also turns a copy + in-place broadcast into a single fused broadcast. If §2's caching lands, sol_avg_last may disappear entirely.

update_average! alone is 3 %–29 % of step! (Table 2), rising as instances get smaller.

3c. Fuse the elementwise passes

At_y = mul!(scratch.x, At, y)
@. sol.x = clamp(x - τ * (c - At_y), lv, uv)     # pass 1 over n
xdiff = @. scratch.x = 2sol.x - x                 # pass 2 over n — reads what pass 1 just wrote

sol.x and xdiff can be produced by one kernel writing two outputs, and the primal half of colaxpby! folds into the same pass; likewise for the dual side. That takes 5.46 broadcasts/iteration down to 2.

---

4. Every colnorm!! / colsum!! drains the pipeline — confirms #89

In the non-batched path colnorm!!(::Number, v) = norm(v) and colsum!!(::Number, v) = sum(v) return CPU scalars, so each one is a blocking device→host transfer. The trace shows 0.36 [copy device to pageable memory] events per iteration — matching 7 per kkt_errors! call × 5 calls per 100 iterations, plus the two in primal_weight_update!! at each restart. That is ~35 pipeline drains per check.

Julia's profiler makes this vivid. On square47, of the 830 samples inside the solver loop:

frame	share of solver samples
kkt_errors!	77 %
↳ colnorm!! → cuBLAS.nrm2 → norm	74 %
↳ CUDACore.nonblocking_synchronize	71 %
__futex_abstimed_wait_common (self)	71 %

The CPU is blocked rather than busy, so on a GPU-bound instance the first drain simply absorbs the queued backlog. But the cost is real when the GPU is not saturated: an isolated norm on a device vector measures 26–30 µs against a 7.2 µs kernel-launch floor and ~1.5 µs for the reduction kernel itself. This is part of why the smallest instances sit at 32–38 % GPU-busy.

Contrast with the tiny reference instance I used to validate the harness (neos5, 63×63): there 83 % of solver samples are inside step! issuing kernels, and GPU-busy/wall is 14 %.

Suggested fix (this is #89): keep the reductions in 0-dimensional device arrays, do the absolute!!/relative!! arithmetic on-device, and read back exactly one scalar per check — the termination decision. That collapses 35 syncs per check into 1. Better still, the 7 reductions inside kkt_errors! are over just two vectors (x-shaped and y-shaped) and could be fused into two multi-output reduction kernels.

Related, the host-side chatter is significant on small instances: cuCtxGetCurrent is called 128.7 times per iteration and cuStreamIsCapturing 87.6 times (9.0 % and 6.5 % of host CUDA-API time on supportcase7). That's CUDA.jl/GPUArrays per-operation bookkeeping, not CoolPDLP's code, but it scales with the op count from §3 — halving GPU ops halves this too.

---

5. GPUSparseMatrixCSR's own SpMV kernel is one-thread-per-row, which is a lottery

The package's spmv_csr! assigns one thread per row. Measured against cuSPARSE on the same preconditioned matrices:

instance	rows	nnz/row	cuSPARSE A·x	GPUSparseMatrixCSR A·x	ratio
neos-3402454-bohle	2,897,380	3.1	130 µs	32 µs	0.24× (4× faster)
neos-5114902-kasavu	961,170	4.4	67 µs	164 µs	2.5×
square47	61,591	444	389 µs	18,412 µs	47×
roi5alpha10n8	4,665	508	24 µs	4,284 µs	182×
supportcase7	6,532	436	25 µs	6,120 µs	249×

The kernel is excellent when rows are short and numerous (a coalesced, perfectly balanced case) and catastrophic when rows are long — one thread serially walking 400+ non-zeros with a strided gather. Since GPUSparseMatrixCSR is the portable path (AMD/Intel/Apple, where there is no cuSPARSE), this is worth fixing: a warp-per-row (CSR-Vector) variant, or a row-block/adaptive partitioning chosen from nnz/row at construction time, would recover most of it. Selecting between scalar and vector kernels on nnz/nrows alone would already cover both ends of this table.

---

6. Two defaults that quietly cap a 120 s run

• max_kkt_passes = 10^5. On square47 this fires at 72 s of a 120 s budget (ITERATION_LIMIT), and the faster instances reach it in single-digit seconds — supportcase7 does 1.29 M iterations in 120 s, i.e. 13× the default cap. Every measurement in this report was taken with max_kkt_passes = 10^9 so that time_limit is the only stopping rule. Consider scaling the default with the time limit, or documenting that the two limits interact. (kkt_passes also doesn't count the 5 evaluations per check — that's kkt_passes undercounts the actual number of KKT evaluations performed #99.)
• Preprocessing counts against time_limit. solve sets starting_time = time() before preprocess, and the MOI wrapper goes through solve. Ruiz + Chambolle–Pock runs entirely on the CPU and costs 4.3 s–9.5 s here, i.e. 3.6 %–7.9 % of a 120 s budget gone before the first PDHG step. Combined with GPU-friendly preconditioning #54 (GPU-friendly preconditioning) that is worth reclaiming. (initialize measures 12.8–13.3 s, but that is almost entirely first-call Julia compilation — the spectral-norm computation it wraps takes 26–780 ms — so it is a one-off per session, not per solve.)

---

7. GPU memory is not a constraint, which is an opportunity

nvidia-smi sampled at 250 ms over each full 120 s run:

	smallest (roi5alpha10n8)	largest (square47)
peak memory.used	587 MiB	1003 MiB
of 46,068 MiB available	1.3 %	2.2 %
mean GPU utilisation	68 %	94 %
mean power	171 W	255 W (350 W cap)
SM clock	2520 MHz (pinned, no throttling anywhere)	2520 MHz

Roughly 250–300 MiB of that is the CUDA context; A and Aᵀ in CSR/Float32/Int32 account for 2·(8·nnz + 4·(m+1)) ≈ 437 MiB on square47, and the solver state (5 PrimalDualSolutions + scratch ≈ 12 vectors) is only ~7.5 MiB.

So even the largest instance in MIPLIB2017 uses 2 % of an L40S. Every "trade memory for speed" fix in this report — cached cuSPARSE descriptors and workspaces, precomputed constant scales, a second matrix layout for the portable kernel — is essentially free. Power never approaches the cap and clocks never throttle, which also confirms the runs are latency/bandwidth-bound rather than thermally limited.

---

Suggested order of work

#	change	measured payoff	existing issue
1	Fix spectral_norm (normalise x0, explicit tol/maxiter, convergence check, √(‖A‖₁‖A‖∞) sanity assert) together with an adaptive step size or divergence guard	up to 45×; but see the supportcase7 divergence — do not ship the first half alone	#95
2	Reuse stats.err in restart_check!; cache the previous candidate error instead of recomputing err_last/err_avg_last	5.6 %–12 % of wall clock	#100
3	Hoist primal_scale / dual_scale into initialize	10 reductions + 10 vector passes + 10 blocking syncs per check	—
4	Keep reductions on-device; one readback per check	dominates the CPU profile; 35 → 1 syncs per check	#89
5	Cache cuSPARSE descriptors + workspace, call cusparseSpMV_preprocess once per matrix	~20 % of GPU time on the smaller instances, ~1 % on the largest	—
6	Pointer-swap instead of copy!(sol_avg_last, sol_avg); fuse the average update into the PDHG kernels	0.4 %–5.6 % of GPU time, plus 3 fewer kernels/iteration	#111
7	Warp-per-row (or adaptive) spmv_csr! for GPUSparseMatrixCSR	up to 249× on long-row matrices; keep the scalar kernel for short-row ones	—
8	Reconsider max_kkt_passes default; move preconditioning to the GPU / start the clock after it	3.6 %–7.9 % of a 120 s budget	#54, #99

Reproducing

~/pdlp-profiling/          # on master01-hpc
├── profile_instance.jl    # the measurement
├── run_array.sh           # sbatch --array=1-10%1, gpu-lvmt
├── specnorm_audit.jl      # §1, including the A/B
├── redundancy_check.jl    # §2
├── aggregate.py           # the tables below
└── results/<instance>/    # result.json, cuda_profile_summary.txt,
                           # julia_profile_*.txt, nvidia_smi.csv, log.txt

---

Appendix A — full measurement tables (all 10 instances)

Table 1 — instances and end-to-end behaviour (120 s budget)

instance	rows m	cols n	nnz	read s	precond s	iters	µs/iter	restarts	status	final rel. KKT err
square47	61,591	95,030	27,329,856	27.0	9.5	167,400	717	18	TIME_LIMIT	9.09e-01
square41	40,160	62,234	13,566,426	15.4	6.8	305,900	392	19	TIME_LIMIT	4.44e-01
neos-3402454-bohle	2,897,380	2,904	8,953,800	16.8	6.7	5,200	485	12	OPTIMAL	9.17e-05
neos-2075418-temuka	349,602	122,304	7,610,261	8.3	5.9	455,900	263	20	TIME_LIMIT	1.00e+00
neos-5052403-cygnet	38,268	32,868	4,898,304	6.2	5.1	669,900	176	32	OPTIMAL	9.91e-05
supportcase19	10,713	1,429,098	4,287,094	9.3	5.2	658,200	182	26	TIME_LIMIT	1.00e+00
neos-5114902-kasavu	961,170	710,164	4,240,376	7.6	5.5	648,000	185	28	TIME_LIMIT	2.98e-02
neos-4647030-tutaki	8,382	12,600	3,953,388	5.8	4.4	21,800	178	21	OPTIMAL	3.54e-05
supportcase7	6,532	138,844	2,845,545	4.8	4.4	1,326,700	90	47	TIME_LIMIT	1.87e-02
roi5alpha10n8	4,665	106,150	2,370,224	5.2	4.3	270,300	90	28	OPTIMAL	9.98e-05

Table 2 — where the wall clock goes, per outer check of 100 iterations (microbenchmarks)

instance	step! µs	100 steps µs	termination_check! µs	restart_check! µs	check overhead	kkt_errors! µs	update_average! µs	update_average! share of step!
square47	683	61258	918	4766	8.5%	1242	20	3%
square41	369	30011	580	3041	10.8%	896	20	6%
neos-3402454-bohle	298	28217	597	2430	9.7%	605	28	9%
neos-2075418-temuka	216	20275	560	2541	13.3%	629	24	11%
neos-5052403-cygnet	116	11075	350	1383	13.5%	352	20	18%
supportcase19	153	14756	387	1541	11.6%	390	23	15%
neos-5114902-kasavu	155	14984	391	1538	11.4%	392	27	17%
neos-4647030-tutaki	95	9042	336	1327	15.5%	337	20	22%
supportcase7	91	8670	334	1312	16.0%	336	19	21%
roi5alpha10n8	70	6014	301	1177	19.7%	300	20	29%

Table 3 — GPU utilisation during the profiled window

instance	GPU busy / wall	device ops / iter	host API calls / iter	nvidia-smi util mean	mem used MiB	power W (cap 350)	SM clock MHz
square47	93%	14.8	300	94%	1003	255	2520
square41	87%	14.8	301	91%	779	244	2520
neos-3402454-bohle	88%	14.8	300	51%	843	164	2520
neos-2075418-temuka	86%	14.8	295	85%	715	233	2520
neos-5052403-cygnet	80%	14.8	287	85%	651	221	2520
supportcase19	81%	14.8	288	86%	715	251	2520
neos-5114902-kasavu	81%	14.8	288	86%	715	245	2520
neos-4647030-tutaki	77%	14.8	286	63%	619	160	2520
supportcase7	32%	14.8	279	58%	587	157	2520
roi5alpha10n8	38%	14.8	278	68%	587	171	2520

Table 4 — SpMV efficiency

instance	A·x µs	Aᵀ·y µs	achieved BW A·x	achieved BW Aᵀ·y	% of 864 GB/s peak	GPUSparseMatrixCSR A·x µs	KA / cuSPARSE
square47	389	352	564 GB/s	624 GB/s	65%	18412	47.32x
square41	177	195	616 GB/s	561 GB/s	71%	11192	63.18x
neos-3402454-bohle	130	139	729 GB/s	598 GB/s	84%	32	0.24x
neos-2075418-temuka	102	94	627 GB/s	671 GB/s	73%	549	5.36x
neos-5052403-cygnet	70	29	569 GB/s	1381 GB/s	66%	185	2.66x
supportcase19	63	67	632 GB/s	688 GB/s	73%	2592	40.82x
neos-5114902-kasavu	67	65	668 GB/s	674 GB/s	77%	164	2.46x
neos-4647030-tutaki	58	24	544 GB/s	1349 GB/s	63%	982	16.82x
supportcase7	25	23	951 GB/s	1061 GB/s	110%	6120	249.04x
roi5alpha10n8	24	22	825 GB/s	922 GB/s	95%	4284	181.92x

Table 5 — GPU time by kernel category (share of total GPU-busy time in the profiled window)

instance	SpMV (cusparse csrmv)	SpMV setup (csr_partition)	elementwise broadcast	D2D copies (copy!/copyto!)	memset (cuSPARSE buffer)	reductions (norm/sum)	scalar readback D2H	other
square47	96.7%	1.0%	1.3%	0.4%	0.3%	0.3%	0.1%	0.0%
square41	93.8%	1.9%	2.3%	0.6%	0.6%	0.6%	0.1%	0.0%
neos-3402454-bohle	79.8%	3.6%	10.2%	3.2%	2.1%	0.8%	0.1%	0.0%
neos-2075418-temuka	86.8%	3.4%	5.3%	2.1%	1.2%	1.0%	0.2%	0.0%
neos-5052403-cygnet	85.0%	4.6%	5.6%	1.5%	1.6%	1.3%	0.3%	0.0%
supportcase19	71.7%	4.9%	14.6%	4.2%	2.8%	1.6%	0.3%	0.0%
neos-5114902-kasavu	70.2%	5.6%	14.1%	5.5%	2.9%	1.6%	0.3%	0.0%
neos-4647030-tutaki	82.6%	4.9%	6.7%	1.8%	1.8%	1.8%	0.4%	0.0%
supportcase7	47.2%	15.0%	21.0%	5.6%	5.2%	5.1%	1.0%	0.0%
roi5alpha10n8	58.5%	12.0%	16.1%	4.4%	4.2%	4.0%	0.8%	0.0%
mean	77.2%	5.7%	9.7%	2.9%	2.3%	1.8%	0.4%	0.0%

Table 6 — calls per PDHG iteration, by kernel category

instance	SpMV (cusparse csrmv)	SpMV setup (csr_partition)	elementwise broadcast	D2D copies (copy!/copyto!)	memset (cuSPARSE buffer)	reductions (norm/sum)	scalar readback D2H	other
square47	2.10	2.10	5.46	2.01	2.11	0.71	0.36	0.00
square41	2.10	2.10	5.46	2.01	2.11	0.71	0.36	0.00
neos-3402454-bohle	2.10	2.10	5.46	2.01	2.11	0.71	0.36	0.00
neos-2075418-temuka	2.10	2.10	5.46	2.01	2.11	0.71	0.36	0.00
neos-5052403-cygnet	2.10	2.10	5.46	2.01	2.11	0.71	0.36	0.00
supportcase19	2.10	2.10	5.46	2.01	2.11	0.71	0.36	0.00
neos-5114902-kasavu	2.10	2.10	5.46	2.01	2.11	0.71	0.36	0.00
neos-4647030-tutaki	2.10	2.10	5.46	2.01	2.11	0.71	0.36	0.00
supportcase7	2.10	2.10	5.46	2.01	2.11	0.71	0.36	0.00
roi5alpha10n8	2.10	2.10	5.46	2.01	2.11	0.71	0.36	0.00

Table 7 — NVTX phase breakdown (share of GPU-busy time / of host time)

instance	steps (dev / host)	termination_check! (dev / host)	restart_check! (dev / host)	restart! (dev / host)
square47	94.8% / 19.4%	1.0% / 72.6%	4.1% / 7.9%	0.0% / 0.1%
square41	94.5% / 27.5%	1.1% / 60.7%	4.4% / 11.6%	0.0% / 0.2%
neos-3402454-bohle	94.0% / 34.0%	1.2% / 55.6%	4.7% / 10.2%	0.1% / 0.2%
neos-2075418-temuka	94.0% / 45.3%	1.2% / 44.5%	4.7% / 9.9%	0.1% / 0.3%
neos-5052403-cygnet	93.5% / 62.5%	1.3% / 24.2%	5.1% / 13.0%	0.1% / 0.3%
supportcase19	93.0% / 60.6%	1.4% / 26.7%	5.5% / 12.4%	0.1% / 0.3%
neos-5114902-kasavu	93.1% / 62.4%	1.3% / 24.8%	5.4% / 12.6%	0.1% / 0.3%
neos-4647030-tutaki	93.3% / 72.0%	1.1% / 14.8%	5.6% / 12.9%	0.1% / 0.3%
supportcase7	88.9% / 81.5%	2.2% / 4.5%	8.7% / 13.7%	0.2% / 0.4%
roi5alpha10n8	90.3% / 81.7%	1.9% / 4.4%	7.7% / 13.6%	0.1% / 0.4%

Table 8 — setup cost and spectral-norm estimate

instance	preprocess (Ruiz, CPU) s	initialize s	spectral_norm default	spectral_norm maxiter=500	ratio tight/default	problem data on GPU MiB	state on GPU MiB
square47	9.5	12.77	29.659	1	0.0337	1282	54
square41	6.8	12.99	26.1568	1	0.0382	642	54
neos-3402454-bohle	6.7	12.86	0.596563	0.915313	1.5343	450	54
neos-2075418-temuka	5.9	12.96	89.2866	0.956015	0.0107	386	54
neos-5052403-cygnet	5.1	12.98	15.1585	1	0.0660	258	54
supportcase19	5.2	12.98	59.5996	0.7964	0.0134	226	54
neos-5114902-kasavu	5.5	13.06	369.33	0.979904	0.0027	226	54
neos-4647030-tutaki	4.4	13.34	0.943598	0.999991	1.0598	194	54
supportcase7	4.4	12.98	29.7089	1	0.0337	162	54
roi5alpha10n8	4.3	12.79	23.1701	1	0.0432	130	54

Appendix B — kernel and CUDA-API breakdowns

Top device operations (3 representative instances)

square47 (2000 iterations profiled, 1339.9 ms GPU-busy)

kernel	calls	calls/iter	total ms	% GPU	mean µs
cuSPARSE SpMV (csrmv_v3_kernel)	4200	2.10	1295.27	96.7%	308.40
cuSPARSE SpMV row partitioning (recomputed each call)	4200	2.10	13.34	1.0%	3.18
broadcast: colaxpby! — running-average update	4000	2.00	5.71	0.4%	1.43
device-to-device copy — copy!(sol_avg_last, sol_avg)	4012	2.01	4.70	0.4%	1.17
memset — cuSPARSE SpMV workspace	4212	2.11	4.42	0.3%	1.05
broadcast: clamp — primal/dual projection	2000	1.00	3.55	0.3%	1.77
broadcast: clamp — primal/dual projection	2000	1.00	3.43	0.3%	1.72
broadcast: elementwise (PDHG update)	2000	1.00	2.97	0.2%	1.49
cuBLAS nrm2 — colnorm!!	824	0.41	2.97	0.2%	3.60
device-to-host scalar readback (forces a sync)	712	0.36	0.87	0.1%	1.22
mapreduce — colsum!!/sum	300	0.15	0.68	0.1%	2.27
mapreduce — colsum!!/sum	300	0.15	0.53	0.0%	1.77

neos-5114902-kasavu (2000 iterations profiled, 302.1 ms GPU-busy)

kernel	calls	calls/iter	total ms	% GPU	mean µs
cuSPARSE SpMV (csrmv_v3_kernel)	4200	2.10	211.95	70.2%	50.46
cuSPARSE SpMV row partitioning (recomputed each call)	4200	2.10	16.95	5.6%	4.04
device-to-device copy — copy!(sol_avg_last, sol_avg)	4012	2.01	16.48	5.5%	4.11
broadcast: colaxpby! — running-average update	4000	2.00	13.59	4.5%	3.40
broadcast: clamp — primal/dual projection	2000	1.00	10.46	3.5%	5.23
memset — cuSPARSE SpMV workspace	4212	2.11	8.68	2.9%	2.06
broadcast: clamp — primal/dual projection	2000	1.00	8.58	2.8%	4.29
broadcast: elementwise (PDHG update)	2000	1.00	6.26	2.1%	3.13
cuBLAS nrm2 — colnorm!!	824	0.41	3.05	1.0%	3.70
mapreduce — colsum!!/sum	300	0.15	1.10	0.4%	3.67
device-to-host scalar readback (forces a sync)	712	0.36	0.86	0.3%	1.20
broadcast: safeprod_left — dual objective (kkt_errors!)	200	0.10	0.81	0.3%	4.06

supportcase7 (2000 iterations profiled, 85.4 ms GPU-busy)

kernel	calls	calls/iter	total ms	% GPU	mean µs
cuSPARSE SpMV (csrmv_v3_kernel)	4200	2.10	40.31	47.2%	9.60
cuSPARSE SpMV row partitioning (recomputed each call)	4200	2.10	12.82	15.0%	3.05
broadcast: colaxpby! — running-average update	4000	2.00	5.83	6.8%	1.46
device-to-device copy — copy!(sol_avg_last, sol_avg)	4012	2.01	4.75	5.6%	1.18
memset — cuSPARSE SpMV workspace	4212	2.11	4.43	5.2%	1.05
broadcast: clamp — primal/dual projection	2000	1.00	4.10	4.8%	2.05
broadcast: clamp — primal/dual projection	2000	1.00	3.29	3.8%	1.64
broadcast: elementwise (PDHG update)	2000	1.00	3.20	3.7%	1.60
cuBLAS nrm2 — colnorm!!	824	0.41	3.09	3.6%	3.75
device-to-host scalar readback (forces a sync)	712	0.36	0.86	1.0%	1.21
mapreduce — colsum!!/sum	300	0.15	0.72	0.8%	2.41
mapreduce — colsum!!/sum	300	0.15	0.54	0.6%	1.81

Top host-side CUDA API calls (largest and smallest)

square47

API call	calls	calls/iter	total ms	% host
cuStreamSynchronize	1424	0.7	1027.88	80.7%
cudaLaunchKernel_v7000	9224	4.6	72.11	5.7%
cuMemcpyDtoDAsync_v2	4012	2.0	42.61	3.3%
cuLaunchKernelEx	11524	5.8	39.56	3.1%
cuCtxGetCurrent	279102	139.6	16.66	1.3%
cudaMemsetAsync_v3020	4200	2.1	16.10	1.3%
cuStreamIsCapturing	175240	87.6	10.75	0.8%
cuMemcpyDtoHAsync_v2	712	0.4	10.22	0.8%
cuMemAllocFromPoolAsync	5212	2.6	10.07	0.8%
cuStreamQuery	31305	15.7	7.58	0.6%

supportcase7

API call	calls	calls/iter	total ms	% host
cuLaunchKernelEx	11524	5.8	36.97	22.9%
cudaLaunchKernel_v7000	9224	4.6	29.98	18.5%
cuMemcpyDtoDAsync_v2	4012	2.0	17.53	10.8%
cudaMemsetAsync_v3020	4200	2.1	14.85	9.2%
cuCtxGetCurrent	257334	128.7	14.58	9.0%
cuStreamIsCapturing	175240	87.6	10.56	6.5%
cuMemAllocFromPoolAsync	5212	2.6	7.80	4.8%
cuMemcpyDtoHAsync_v2	712	0.4	7.03	4.3%
cuMemFreeAsync	4200	2.1	5.76	3.6%
cuOccupancyMaxPotentialBlockSize	11824	5.9	5.64	3.5%