T15: device LDLᵀ/LDLᴴ in regimes A/B/C, diagonal solve, pivot statistics

[michel2323]

Closes #15

Task: T15 — GPU LDLᵀ/LDLᴴ in regimes A/B/C, diagonal solve, statistics

Recovered by the owner from the implementer run 37032539796: the run took 81 minutes, the GitHub App token had expired when it pushed, so Claude printed the commit as a format-patch into the execution log; it was reassembled (sha256 47a160d8…647ae6 verified), applied on main with git am --3way and pushed as task/T15-gpu-ldlt. Commit c1906f1 of the runner = 9f738bd here. The Report block in TASKS.md is the implementer's.

What was built: src/numeric/ldlt.jl: front_ldlt_kernel! (regimes B/C, fused assembly) and subtree_ldlt_kernel! (regime A, local memory) with the T14 in-front Bunch–Kaufman pivoting, threshold acceptance, perturbation and pivot_sign, same pivot sequence as ref_ldlt!; reduce_stats!, factorize_ldlt!; factorize! dispatches on the structure. Solve: local pivot orders in the sweeps, Lᵀ for complex symmetric "S", diagonal_sweep! with 1×1/2×2 blocks. API: structures "S"/"H", "solve_diag", getparam npivots/inertia/pivot_stats, diag = D, pivot_sign checks, ldlt/ldlt!, logabsdet for LDLᵀ, CUDA ldlt(::CuSparseMatrixCSR) and wrappers. Tests: test/test_numeric_ldlt.jl, ported "Symmetric -- Hermitian" parts of cudss_execution and cudss_generic, cudss_solver for "S"/"H".

Tests: SDS_TEST_GPU=0 julia --project=. -e 'using Pkg; Pkg.test()' on the runner (Julia 1.13.1): 59877 pass / 0 fail / 0 broken (11 min 54 s); test_numeric_ldlt 1822, ported 8104; with CI's flags test_numeric_ldlt + test_api 3768 / 0. CUDA: pending CI on the PR.

Deviations from PLAN.md / TASKS.md (Report, status [!]):

• No vendor sytrf for regime C: sytrf picks its own pivot order on F₁₁ alone, so D/pivot_kind would no longer equal the reference as the task's first test requires. Regime-C fronts run the KA kernel for "S" and "H". Owner decision pending on Device LDLᵀ: serial pivot search, one workgroup per regime-C front, no vendor sytrf (T25 follow-up from T15) #75.
• Performance gaps filed as Device LDLᵀ: serial pivot search, one workgroup per regime-C front, no vendor sytrf (T25 follow-up from T15) #75 for T25: serial pivot search (13.5 s vs 12.9 s reference on kkt_matrix(3000,1000,1e-8), CPU backend), one workgroup per regime-C front, regime B keeps F₁₁ in global memory.
• Reference change: Bunch–Kaufman's α computed in the working precision so the device matches it and Metal (no Float64) can run it.
• Regime-A kernel cut to 32 arguments (at 34 the KA CPU launch is not specialized and allocates); Numeric has three new fields; older tests that expected "S"/"H", "solve_diag", npivots, inertia, pivot_stats to be refused now assert their values.

Follow-up issues opened: #75.

🤖 Generated with Claude Code

https://claude.ai/code/session_014NfiFN5Ukdcnbh87uff2Dg

---

Generated by Claude Code

[claude]

nit (non-blocking): opts.pivot_sign is parsed to a host Vector{Int} (_parse_int_vector), so this is a copyto! from Vector{Int} into a device Int8 vector. On the CPU backend that is a plain Base copy; on CUDA/AMDGPU it relies on GPUArrays' element-type-converting host→device path (and on the host conversion allocating per factorization). Converting once, e.g. copyto!(N.psign, convert(Vector{Int8}, ps)) or storing Vector{Int8} in Options, avoids depending on that fallback. CI on the GPU runners will show whether it matters.

[claude]

nit (non-blocking): for Hermitian fronts the device forces the diagonal of the contribution block to be real here, while ref_ldlt! updates F₂₂ with mul! and keeps the (rounding-level) imaginary parts. Both drop the imaginary part at pivot time (T(real(x))), but _fabs(F, c, c) in the pivot tests sees abs(x) of the full complex value in the reference and of the real-only value here. That is one more O(eps) source of a near-tie pivot difference between the device and the oracle on top of the FMA one the Report lists; worth recording in the Report's open issues so a future bitwise mismatch is diagnosed quickly.

[claude]

nit (non-blocking): this testset and the ones after it ('D'/'N', determinism, solve sweeps, MadNLP loop, public API, complex symmetric) draw with rand without a Random.seed!(666) at the start of the testset. runtests.jl seeds once per file, so the draws are deterministic in a full run but depend on the testset order and on ELTYPES/BACKENDS; test_reference_ldlt.jl reseeds per testset for that reason. Adding the seed at the top of each testset keeps the per-case draws fixed if a testset is reordered or a backend is added.

[claude]

nit (non-blocking): for "SPD"/"HPD" this reduces the per-front stats whenever "inertia"/"npivots"/"pivot_stats" is read, including after a failed factorization (info ≠ 0), where the counts of the fronts after the failure are stale/garbage. _need_stage(STAGE_FACTORIZED) does not check solver.info. cuDSS behaves similarly, so this is not blocking, but documenting "valid only when info == 0" in the getparam table (or returning the counts of the successful prefix only) would avoid a surprise in MadNLP's inertia logic.

[claude]

VERDICT: APPROVE (head 9f738bd)

Checked

• Scope: only T15 deliverables; PLAN.md untouched; TASKS.md changed only in the T15 section (status [!], Report); no Manifest.toml; phase and parameter strings ("solve_diag", "npivots", "inertia", "pivot_stats", "diag", "pivot_sign", structures "S"/"H") as in PLAN 1.2-1.4; ldlt / ldlt! mirror CUDSS.jl ("S" real, "H" complex).
• Tests: test/test_numeric_ldlt.jl covers every listed test (D, pivot_kind, piv and stats equal to ref_ldlt! on the T14 matrices under default, A+B+C and C-only analyses; residuals within tol(T); exact inertia; perturbation and pivot_sign including a zero pivot in a regime-C root front, asserted with snparent == 0 and takes_c_path; bitwise refactorization determinism and deterministic solves on every backend; the MadNLP loop through the handle layer on the three T14 orderings). Ported "Symmetric -- Hermitian" parts of cudss_execution and cudss_generic present; loops over BACKENDS, ELTYPES, INTTYPES; tol(T); shared helpers in test/utils.jl and test/matrices.jl. No @test_broken. The adapted older assertions ("S"/"H", "solve_diag", npivots/inertia/pivot_stats now asserted instead of refused) are required by the new features; nothing loosened.
• Conventions: generic T and INT; no @allowscalar; 1-D workgroups, @LocalMem sized by constants/Val; no atomics, no device allocation and no host synchronization in factorize_ldlt! and the sweeps (totals copied only in getparam); errors from src/errors.jl; CUDA code only in ext/; docstrings on the new public functions.
• Correctness, traced against src/reference/ldlt.jl: Bunch-Kaufman choice, threshold tests, _best_1x1 fallback, perturbation sign (psign[perm[piv[g]]] = the reference signof(lp[k])), the lower-triangle symmetric interchange (LAPACK-style, conj on the mirrored entries, race-free per row), 1x1/2x2 updates and scaling in the same operation order as the reference, W = L21 D in the contribution-block update, the regime-A in-place trailing-triangle offset (matches subtree.jl), inertia/n2x2/nzero counting, stats reduction. Solve: the local pivot order is applied after the pulls of the children and before the TRSV, restored after the backward TRSV; the scratch rows 1:w of the U buffer of a front are indeed unused by the sweeps; the 2x2 diagonal solve is the correct inverse; transposed (not conjugated) solves for complex symmetric on both the kernel and the dense path; logabsdet from D.
• Report: matches the diff (no vendor sytrf, panel in global memory, BUNCH_KAUFMAN_ALPHA in working precision, three new Numeric fields, 32-argument kernel, 64 B allocation). Follow-up #75 exists and is triaged. Counts are plausible (+7151 over T14) but could not be rerun here (running the suite needs an approval this unattended session cannot grant); the CPU-only and CUDA CI jobs are still pending and the pipeline handles them. The red "fix" check is a cancelled, superseded "Claude review" run, not a test job.

Deviation for the owner (documented, not blocking): the task text asks for regime-C "S" fronts via vendor sytrf with a post-check and a KA redo; the PR runs the KA in-front kernel for every regime-C front because the first listed test (D and pivot_kind equal to the T14 reference) cannot hold with the pivot sequence of sytrf. The suggested plan change and #75 are recorded in the Report.

Blocking findings: none.

Non-blocking inline notes: the copyto! of pivot_sign from a host Vector{Int} into a device Int8 vector relies on an element-type-converting copy; forcing a real diagonal in the Hermitian contribution-block update is a second O(eps) source of near-tie pivot differences against the oracle; several testsets do not reseed; inertia/npivots for Cholesky after a failed factorization return stale counts.

[github-actions]

CI-FIX round 1 of 3: https://github.com/exanauts/SparseDirectSolver.jl/actions/runs/37042438119 failed on 9f738bd; Claude is diagnosing it.

[claude]

CI fix round 1 (run 37042438119, commit 103d790)

Root cause: only test-gpu (cuda) failed: 12 of 66438 checks, all in perturbation and pivot_sign (CUDA, Float32) (test/test_numeric_ldlt.jl:96 and :122). Those two assertions compared the device D bitwise with the host reference ref_ldlt! (Nh.d == Nr.d). On CUDA, Float32 Schur-complement updates get contracted to FMA, so the regular (non-perturbed) pivots of the SPD part differ from the host by rounding. piv and pivot_kind are unaffected. The CPU job and every other CUDA eltype passed. The library code is not at fault: every other device-vs-reference check in the file already compares D with d_error(Nh, Nr) <= panel_tol(T).

Fix (test only, library unchanged): these two places now check piv and pivot_kind with ==, D within panel_tol(T) (100·eps) like the rest of the file, and the perturbed D entries still bitwise against the reference. The perturbed entries are what this test is about, and the existing exact ±ε check on row j is still there. No test was removed or marked broken. The T15 Report notes the change.

Verification: SDS_TEST_GPU=0 Pkg.test(): 59973 pass / 0 fail / 0 broken. I could not run CUDA here, so the CUDA result will come from CI on this push.

[claude]

VERDICT: APPROVE (head 103d790)

Checked (full re-read of the diff against origin/main, the T15 task text, AGENTS.md and PLAN §1.4/§1.7/§2.4/§3.3)

• Scope: only T15 deliverables; PLAN.md untouched; TASKS.md changed only in the T15 section (status [!], Report, CI-fix note); no Manifest.toml; phase and parameter strings ("solve_diag", "npivots", "inertia", "pivot_stats", "diag", "pivot_sign", structures "S"/"H") as in PLAN §1.2–§1.4; ldlt/ldlt! mirror CUDSS.jl ("S" real, "H" complex); CUDA wrappers in ext/ only.
• Correctness: src/numeric/ldlt.jl reproduces ref_ldlt! step for step on the stored lower triangle (_lt_choose = _choose_pivot incl. the BK α/σ tests, threshold acceptance against rows to f, _best_1x1 fallback; LAPACK-style symmetric interchange with conjugated mirrors, disjoint writes per work item; 1×1/2×2 D storage d[g], d[n+g]; W = L₂₁D with the 2×2 blocks and conj/plain mirrors for "H" vs complex "S"; identity block of L for 2×2 pivots; perturbation sign from psign[perm[piv[g]]] = the reference signof(lp[k]); stats incl. negative-determinant 2×2 counting one of each). Packed-front offsets of regime A match _subtree_write!/_packed; the stack block offsets of regimes B/C match _front_syrk!/_extend_add_child!. Solve: local pivot order applied before the forward TRSV and undone after the backward TRSV through rows 1:w of the fronts Uslice (unused by the update buffers and the parent pull);Lᵀfor complex symmetric;_sv_diag_kernel!identical to the reference 2×2 solve, one writer per block. No atomics in assembly/extend-add/update, fixed accumulation order: deterministic. No@allowscalar, no vendor calls outside ext/, @LocalMemsized byVal, kernels ≤ 32 arguments, no hard-coded Float64/Intin device code; index arithmetic inIntfromINT/Int32` sources.
• Tests: every listed test is present and loops over BACKENDS/ELTYPES (and INTTYPES where relevant): D/piv/pivot_kind/stats vs the reference under default, A+B+C and C-only analyses, residuals ≤ tol(T), exact inertia vs eigenvalues, perturbation and pivot_sign with the zero pivot forced into a regime-C root front (snparent == 0, takes_c_path), bitwise refactorization determinism and deterministic solves on every backend, the MadNLP loop through the handle layer on the three T14 orderings, the ported "Symmetric -- Hermitian" execution/generic tests and cudss_solver for "S"/"H". No @test_broken. Helpers in test/utils.jl/test/matrices.jl; tol(T)/panel_tol(T); runtests.jl seeds per file.
• CI-fix round 1 (9f738bd → 103d790): test-only. The two bitwise Nh.d == Nr.d checks in "perturbation and pivot_sign" became d_error ≤ panel_tol(T) plus bitwise equality of the perturbed entries, consistent with the rest of the file; the exact ±ε check at row j, piv/pivot_kind equality and the inertia checks stay. Bitwise equality of the whole of D was never the contract (the device accumulates in a different order than the references BLAS mul!` even on the CPU), so this is not a weakened test. The Report records it.
• Report: claims match the diff; the previous heads test-github-cpuonlyjob passed and the CUDA job failed exactly the 12 checks the Report describes. Deviation (no vendorsytrffor regime C, both B and C throughfront_ldlt_kernel!) is documented with its reason (the tasks first test requires the reference pivot sequence, which sytrf cannot give) and tracked in #75 (found-by-agent, triaged), with suggested plan changes. Status [!] is appropriate; the owner should decide whether PLAN §2.4/§3.3 adopt the KA regime-C path or whether a later task reintroduces sytrf with a documented different pivot sequence.

Blocking findings: none.

Non-blocking (nits, from the previous round, still open): copyto!(N.psign, ps) converts/copies the host pivot_sign on every factorization (a host→device copy before the kernels; caching it or storing Int8 in Options avoids it); the device forces the Hermitian CB diagonal real while the reference does not (one more O(eps) near-tie source worth listing under open issues); per-testset Random.seed!(666) in test_numeric_ldlt.jl; "inertia"/"npivots" for Cholesky after a failed factorization return stale counts (document info == 0). New: the Report line "Julia 1.10 not run locally (CI)" is stale, the package requires Julia 1.13.

CI on this head (Run tests CPU and CUDA jobs) is still pending; the pipeline gates the merge on it.