From 2ef3edd8bde1ba0aad82933838368e0042be2c45 Mon Sep 17 00:00:00 2001 From: Michel Schanen Date: Mon, 20 Jul 2026 11:33:01 -0500 Subject: [PATCH 1/2] [MadNLPGPU] cuDSS: recover inertia under matching; guard matching off ubatch cuDSS (through 0.8) stops reporting a usable inertia once matching is enabled: it returns (0, 0) with info == 0 even though the factorization is correct. InertiaBased/InertiaAuto then bump regularization forever and dive into restoration. Recover the inertia from the sign counts of the factor diagonal ("diag" data parameter) whenever cudss_matching is on, using a buffer cached on the solver. Exact only for 1x1 pivots, which holds for the quasi-definite condensed KKT systems this solver targets. Also skip matching on uniform-batch solvers: setting matching_alg on a ubatch solver fails the analysis phase with CUDSS_STATUS_NOT_SUPPORTED, so the two-stage per-scenario batch solver must never enable it. --- lib/MadNLPGPU/ext/MadNLPGPUCUDAExt/cudss.jl | 29 +++++++++++++++++++-- 1 file changed, 27 insertions(+), 2 deletions(-) diff --git a/lib/MadNLPGPU/ext/MadNLPGPUCUDAExt/cudss.jl b/lib/MadNLPGPU/ext/MadNLPGPUCUDAExt/cudss.jl index 28508a54c..63fe54c04 100644 --- a/lib/MadNLPGPU/ext/MadNLPGPUCUDAExt/cudss.jl +++ b/lib/MadNLPGPU/ext/MadNLPGPUCUDAExt/cudss.jl @@ -47,7 +47,12 @@ function set_cudss_options!(solver::CUDSS.CudssSolver, opt::CudssSolverOptions) CUDSS.cudss_set(solver, "pivot_threshold", opt.cudss_pivot_threshold) end if opt.cudss_matching - if pkgversion(CUDSS) < v"0.8" + if solver.matrix.nbatch > 1 + # cuDSS matching (`matching_alg`) fails the *analysis* phase with + # CUDSS_STATUS_NOT_SUPPORTED on a uniform-batch solver (through cuDSS 0.8), + # so never enable it there (e.g. the two-stage per-scenario batch solver). + Base.@warn "cuDSS matching is not supported on uniform-batch (ubatch) solvers; ignoring `cudss_matching = true` for this batched solver." maxlog = 1 + elseif pkgversion(CUDSS) < v"0.8" CUDSS.cudss_set(solver, "use_matching", 1) if opt.cudss_matching_alg != "default" CUDSS.cudss_set(solver, "matching_alg", opt.cudss_matching_alg) @@ -90,6 +95,7 @@ mutable struct CUDSSSolver{T, V} <: MadNLP.AbstractLinearSolver{T} x_gpu::CUDSS.CudssMatrix{T} b_gpu::CUDSS.CudssMatrix{T} buffer::V + diag::V opt::CudssSolverOptions logger::MadNLP.MadNLPLogger @@ -151,9 +157,13 @@ function CUDSSSolver( # Always allocate it to support dynamic updates to opt.cudss_ir buffer = CuVector{T}(undef, n * nbatch) + # Scratch for the factor diagonal, used to recover the inertia when matching is on + # (cuDSS misreports it as (0, 0)); only the nbatch == 1 path ever queries inertia. + diag = CuVector{T}(undef, n) + return CUDSSSolver( solver, csc, - x_gpu, b_gpu, buffer, + x_gpu, b_gpu, buffer, diag, opt, logger, ) end @@ -202,6 +212,21 @@ function MadNLP.inertia(M::CUDSSSolver) elseif M.opt.cudss_algorithm == MadNLP.LDL # N.B.: cuDSS does not always return the correct inertia. if info == 0 + if M.opt.cudss_matching + # cuDSS (through 0.8) reports inertia (0, 0) whenever matching is enabled, + # even though the factorization is correct — trusting it sends + # InertiaBased/InertiaAuto into an endless regularization bump and then + # restoration. Recover the inertia from the sign counts of the factor + # diagonal instead. Caveat: exact only for 1×1 pivots, which holds for the + # (quasi-definite) condensed KKT family this solver targets. + CUDSS.cudss_set(M.inner, "diag", M.diag) + CUDSS.cudss_get(M.inner, "diag") + d = Array(M.diag) + z = zero(eltype(d)) + npos = count(>(z), d) + nneg = count(<(z), d) + return (npos, n - npos - nneg, nneg) + end (k, l) = CUDSS.cudss_get(M.inner, "inertia") @assert 0 ≤ k + l ≤ n return (k, n - k - l, l) From 90102473cf1bf8384445f779eeccf5bb97e626f4 Mon Sep 17 00:00:00 2001 From: Michel Schanen Date: Thu, 23 Jul 2026 11:11:10 -0500 Subject: [PATCH 2/2] [MadNLPGPU] cuDSS: guard the diag-based inertia recovery behind an explicit LDL check Review feedback on #633 (amontoison): the factor diagonal D only determines the inertia for an LDL' factorization, so hard-check that before dumping it. Extract the recovery into inertia_from_diag so the guard travels with the code, and document why the dumped "diag" (D of the congruence P*S*A*S*P' = L*D*L') carries A's inertia under matching (frapac's question). Co-Authored-By: Claude Fable 5 --- lib/MadNLPGPU/ext/MadNLPGPUCUDAExt/cudss.jl | 31 ++++++++++++++------- 1 file changed, 21 insertions(+), 10 deletions(-) diff --git a/lib/MadNLPGPU/ext/MadNLPGPUCUDAExt/cudss.jl b/lib/MadNLPGPU/ext/MadNLPGPUCUDAExt/cudss.jl index 63fe54c04..3f7eee0e6 100644 --- a/lib/MadNLPGPU/ext/MadNLPGPUCUDAExt/cudss.jl +++ b/lib/MadNLPGPU/ext/MadNLPGPUCUDAExt/cudss.jl @@ -193,6 +193,25 @@ end MadNLP.input_type(::Type{CUDSSSolver}) = :csc MadNLP.default_options(::Type{CUDSSSolver}) = CudssSolverOptions() MadNLP.is_inertia(M::CUDSSSolver) = (M.inner.matrix.nbatch == 1) # Uncomment if MadNLP.LU is supported -- (M.opt.cudss_algorithm ∈ (MadNLP.CHOLESKY, MadNLP.LDL)) + +# Recover the inertia from the sign counts of the factor diagonal D, dumped via the +# cuDSS "diag" data parameter ("Diagonal of the factorized matrix", i.e. D of P·S·A·S·Pᵀ +# = L·D·Lᵀ — a congruence of A, so the sign counts equal A's inertia by Sylvester's law, +# matching scaling/permutation included). Only meaningful for an LDLᵀ factorization, +# hence the hard check. Exact only for 1×1 pivots, which holds for the (quasi-definite) +# condensed KKT family this solver targets. +function inertia_from_diag(M::CUDSSSolver) + @assert M.opt.cudss_algorithm == MadNLP.LDL "the factor diagonal D only determines the inertia for an LDLᵀ factorization" + n = size(M.tril, 1) + CUDSS.cudss_set(M.inner, "diag", M.diag) + CUDSS.cudss_get(M.inner, "diag") + d = Array(M.diag) + z = zero(eltype(d)) + npos = count(>(z), d) + nneg = count(<(z), d) + return (npos, n - npos - nneg, nneg) +end + function MadNLP.inertia(M::CUDSSSolver) @assert M.inner.matrix.nbatch == 1 n = size(M.tril, 1) @@ -216,16 +235,8 @@ function MadNLP.inertia(M::CUDSSSolver) # cuDSS (through 0.8) reports inertia (0, 0) whenever matching is enabled, # even though the factorization is correct — trusting it sends # InertiaBased/InertiaAuto into an endless regularization bump and then - # restoration. Recover the inertia from the sign counts of the factor - # diagonal instead. Caveat: exact only for 1×1 pivots, which holds for the - # (quasi-definite) condensed KKT family this solver targets. - CUDSS.cudss_set(M.inner, "diag", M.diag) - CUDSS.cudss_get(M.inner, "diag") - d = Array(M.diag) - z = zero(eltype(d)) - npos = count(>(z), d) - nneg = count(<(z), d) - return (npos, n - npos - nneg, nneg) + # restoration. Recover the inertia from D instead. + return inertia_from_diag(M) end (k, l) = CUDSS.cudss_get(M.inner, "inertia") @assert 0 ≤ k + l ≤ n