Skip to content

[MadNLPGPU] cuDSS: recover inertia under matching; guard matching off ubatch #633

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
michel2323 wants to merge 2 commits into
base: master
Choose a base branch
from
Open

[MadNLPGPU] cuDSS: recover inertia under matching; guard matching off ubatch #633

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
2ef3edd
[MadNLPGPU] cuDSS: recover inertia under matching; guard matching off…
michel2323 Jul 20, 2026
9010247
[MadNLPGPU] cuDSS: guard the diag-based inertia recovery behind an ex…
michel2323 Jul 23, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
40 changes: 38 additions & 2 deletions lib/MadNLPGPU/ext/MadNLPGPUCUDAExt/cudss.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

@amontoison amontoison Jul 22, 2026 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

NVIDIA is using an external library for the matching (HSL MC64).
They probably didn't extended it for batch...

@michel2323 michel2323 Jul 23, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That matches what we observe: the ubatch analysis phase itself returns CUDSS_STATUS_NOT_SUPPORTED as soon as matching_alg is set, so the guard just skips matching with a one-time warning instead of failing the whole batched solver.

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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -183,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)
Expand All @@ -202,6 +231,13 @@ 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

@michel2323 michel2323 Jul 21, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@frapac @sshin23 @amontoison Is this what we want to do, or do we want to automatically switch to inertia-free if matching is enabled?

@frapac frapac Jul 21, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think your solution is the most appropriate. Are we sure M.diag is the diagonal factor?

@amontoison amontoison Jul 22, 2026 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes but we should always check that the factorization is an LDL' before dumping the diagonal D.

@michel2323 michel2323 Jul 23, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes — cuDSS's "diag" data parameter returns the "Diagonal of the factorized matrix", i.e. the D of P·S·A·S·Pᵀ = L·D·Lᵀ (a congruence of A, so its sign counts equal A's inertia by Sylvester's law, matching scaling/permutation included).

Verified empirically on an RTX 4080 with cuDSS 0.8 on KKT-like symmetric indefinite matrices ([H Aᵀ; A -δI], n=70 and n=280): the sign counts of the dumped D match the eigenvalue-count inertia exactly, both with matching off (where it also agrees with cuDSS's native "inertia") and with matching on (where native reports (0, 0)). End-to-end, SparseCondensedKKTSystem + InertiaBased + cudss_matching=true now converges in the identical iteration count/objective as with matching off.

@amontoison addressed in 9010247: the recovery is extracted into inertia_from_diag with a hard @assert cudss_algorithm == MadNLP.LDL, so the LDLᵀ check travels with the code (previously it was only implied by the enclosing branch).

# 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 D instead.
return inertia_from_diag(M)
end
(k, l) = CUDSS.cudss_get(M.inner, "inertia")
@assert 0 ≤ k + l ≤ n
return (k, n - k - l, l)
Expand Down
Loading