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

TASKS.md

@@ -1983,7 +1983,7 @@ whenever `atomic_add` is `false`. Close #36 in this task's PR.
 
 ---
 
-## T15 — GPU LDLᵀ/LDLᴴ in regimes A/B/C, diagonal solve, statistics   `[ ]`
+## T15 — GPU LDLᵀ/LDLᴴ in regimes A/B/C, diagonal solve, statistics   `[!]`
 
 **Reads**: PLAN §2.4, §3.3, §2.6 (sytrf), §1.4.
 
@@ -2010,6 +2010,135 @@ whenever `atomic_add` is `false`. Close #36 in this task's PR.
 * MadNLP-style loop from T14 on the GPU; the ported CUDSS.jl symmetric tests
   pass for `T ∈ ELTYPES`.
 
+### Report
+
+- Status: [!] (done on the CPU backend; regime C runs the KA kernel for `"S"` too, no vendor `sytrf`, see deviations;
+  CUDA from CI)
+- What was built:
+  - `src/numeric/ldlt.jl` (new): the T14 in-front factorization as KA kernels, pivot for pivot.
+    `front_ldlt_kernel!` (regimes B and C): one workgroup (`LDLT_WORKGROUP = 128`) per front of a launch group,
+    fused zero/scatter/owner-pull extend-add (the T09 device functions), then the pivoted factorization of the
+    `w` fully-summed columns in place in the panel, then `F₂₂ ← F₂₂ − (L₂₁D)L₂₁ᴴ` into the front's packed block on
+    the update stack (`W = L₂₁D` formed on the fly from `d`/`pivot_kind`, 2×2 blocks included), unit diagonal.
+    `subtree_ldlt_kernel!` (regime A): the T11 subtree walk with packed fronts and blocks in `@localmem`, the same
+    pivot steps on the packed front, the block update in place in the trailing triangle. Per pivot step (a loop of
+    `w` iterations with a uniform trip count; a 2×2 pivot leaves the last ones idle): work item 1 chooses the pivot
+    (`_lt_choose` = the reference's `_choose_pivot`, `_accept_1x1`, `_accept_2x2`, `_best_1x1` on the lower
+    triangle), the workgroup swaps rows/columns (LAPACK-style symmetric interchange of the lower triangle), work
+    item 1 stores D (perturbation `±ε` with the `pivot_sign` of the original row, `psign[perm[piv[g]]]`), the pivot
+    kinds and the counts, the workgroup applies the rank-1/rank-2 update to the remaining fully-summed columns, then
+    scales the pivot columns into L. Pivot choice, control words and per-front counts live in `@localmem` and go
+    to `numeric.stats`. No atomics, deterministic. `_LDLTDevice{R, H}` carries `pivot_type`, `pivot_threshold`,
+    `pivot_epsilon` and the scaled flag (`H`: Hermitian/real vs complex symmetric as a type parameter).
+    `factorize_ldlt!(numeric, symbolic, nzval; opts)` (copies `pivot_sign` to `numeric.psign`, `abs_max!` for
+    `pivot_epsilon_alg = "algo1"`, one launch per regime-A class and per B/C launch group, `reduce_stats!`; returns 0,
+    no host synchronization), `reduce_stats!` (one workgroup, `@localmem` tree reduction of the per-front statistics
+    into `numeric.totals`), `pivot_totals` (reads them), `abs_max!`. `factorize!(N, S, nzval; impl, opts)`
+    dispatches to it for `"S"`/`"H"`.
+  - `Numeric` gains `totals` (6 `Int64`), `psign` (`n` `Int8`) and `aux` (one `T`: `max|aᵢⱼ|`); `allocate_numeric`,
+    `host_numeric`, `memory_estimates` slot 11 updated. `ref_ldlt!` fills `totals` too.
+  - `src/solve/sweeps.jl`: the forward kernel reorders a front's slice of `Y` by its local pivot order before its
+    TRSV, the backward kernel restores it after its TRSV (`Val(LDL)`; scratch: rows `1:w` of the front's
+    gather-list buffer in `U`, unused before), as `ref_solve_ldlt!` does; `Val(CONJ) = false` solves with `Lᵀ` for
+    complex symmetric `"S"`; the regime-C path gets a pivot-order kernel each way and `'T'` dense calls for complex
+    `"S"`. The unit diagonal is stored, so the TRSVs divide by 1 (exact) and need no change. `diagonal_sweep!` is
+    one launch (1×1 and 2×2 blocks of D); identity for Cholesky.
+  - API (`src/solver.jl`, `src/generic.jl`, CUDA extension): structures `"S"` (complex: complex symmetric) and
+    `"H"` accepted (`"G"` refused, T19); `execute!("solve_diag", …)`; `getparam` for `"npivots"` (`INT`), `"inertia"`
+    (`Tuple{INT, INT}`, as CUDSS.jl), `"pivot_stats"` (NamedTuple) from the device totals (Cholesky: reduced on
+    demand, `(npos, 0)`), `"diag"` = `d[1:n]` for LDLᵀ; `setparam!(solver, "pivot_sign", v)` checks the length
+    (host or device vector), the next factorization copies it to the device; `ldlt(A::CSR; view, check)`,
+    `ldlt!(solver, A)`, `logabsdet` for LDLᵀ/LDLᴴ (`det D`, 2×2 blocks, sign ±1 or complex); CUDA extension:
+    `ldlt(::CuSparseMatrixCSR)` and its `Symmetric`/`Hermitian` wrappers.
+  - Tests: `test/test_numeric_ldlt.jl` (new); `test/ported/cudss_execution.jl` and `cudss_generic.jl` gain the
+    "Symmetric -- Hermitian" parts; `cudss_solver.jl` runs `"S"`/`"H"` too and checks `"npivots"`/`"inertia"`;
+    helpers `ldlt_setup`, `d_error`, `ldlt_alloc_budget`, `madnlp_inertia_loop_device` in `test/utils.jl`,
+    `solve_alloc_budget(…; ldlt)`; `MADNLP_ORDERINGS` moved from `test_reference_ldlt.jl` to `test/matrices.jl`.
+- Tests: `SDS_TEST_GPU=0 julia --project=. -e 'using Pkg; Pkg.test()'` (Julia 1.13.1, ubuntu-latest, CPU):
+  59877 pass / 0 fail / 0 broken (11 min 54 s); test_numeric_ldlt 1822 pass (31 s in the suite, ~6 min
+  standalone, mostly compilation), test_ported 8104 pass. Under CI's flags (`coverage = true`, `--check-bounds=yes`)
+  test_numeric_ldlt + test_api: 3768 pass / 0 fail. Julia 1.10 not run locally (CI). Covered, per backend and `T ∈ ELTYPES`
+  (`"S"` real, `"H"` complex):
+  - Device factor vs `ref_ldlt!` on the T14 matrices (`random_symindef(400,0.01)`, `kkt(300,100,1e-2)`,
+    `kkt(300,100,1e-8)` default and interleaved ordering, `random_symindef(300,0.02)`, indefinite-H KKT δ = 0,
+    the 2×2 case `hessian_scale = 1e-3` interleaved, complex symmetric for complex `T`), each under three analyses
+    (default, regimes A+B+C with small budgets and C thresholds, C only): `piv == `, `pivot_kind ==`, `stats ==`,
+    totals `==`, D and panels within `panel_tol(T)` (100 eps) of the reference; inertia equals the eigenvalue
+    signs on the `n ≤ 300` cases; `relres ≤ tol(T)` for `nrhs ∈ {1, 5}`, atomic and deterministic sweeps
+    (`kkt(300,100,1e-8)` ComplexF32 after one refinement step, as in T14 / #66).
+  - Regime-A groups with Int64 maps; views `'U'`/`'F'` give `==` factors.
+  - `singular_block_matrix` (structurally zero and stored zero, both `pivot_sign` signs) under the default analysis
+    and with the zero pivot in a regime-C root front (`"algo2"`; asserted `snparent == 0`, `takes_c_path`):
+    `piv`, kinds `==` the reference, D within `panel_tol(T)`, perturbed D entries `==` the reference, perturbed entry `±ε` exactly at row `j`, inertia, consistent-RHS relres ≤ 1e-6;
+    `pivot_epsilon`, scaled `"algo1"` (device `max|a|`), `'N'`, `'D'`; wrong `pivot_sign` length.
+  - `'D'`/`'N'` on quasi-definite KKT (no 2×2, no perturbation, inertia `(nh, nj)`, `'N'` no interchange).
+  - Determinism: refactorization with new values and back gives `==` factor, D, `piv`; deterministic solves `==`;
+    no allocation beyond the budgets in `factorize!` (64 B measured) and the solve (0 B) on the CPU backend.
+  - Sweeps one at a time (permute, forward, diagonal, backward) `==` the solve, against `extract_ldlt`'s `L` and
+    `D`; every dense implementation on the regime-C path.
+  - MadNLP-style inertia loop through the handle layer (`DirectSolver`, `nzval` rewritten, `"refactorization"`,
+    `getparam("inertia")`, `"npivots"`) for the three T14 orderings: same iterations and `δw` as the reference loop,
+    final solve `relres ≤ tol(T)`.
+  - Public API for `INT ∈ INTTYPES`: phases, `npivots`/`inertia`/`pivot_stats` (types, values against the
+    reference, phase-order errors), `diag`, `getparam!`, the five sub-phases `==` `"solve"`, `logabsdet` against the
+    dense one (also complex symmetric), `pivot_sign` set through the solver (length error, device vector, perturbed
+    sign and inertia), `ldlt(::CSR)`.
+  - Ported: `cudss_execution` "Symmetric -- Hermitian" (views F/L/U × pivoting C/R/N, refactorization after a
+    diagonal shift), `cudss_generic` `ldlt`/`ldlt!`/`ldiv!`/`\`, `cudss_solver` for `"S"`/`"H"`; `T ∈ ELTYPES`,
+    `INT ∈ INTTYPES`.
+  - CUDA/AMDGPU: pending CI on the PR. CI fix round 1: the first CUDA run failed 12 checks in "perturbation and
+    pivot_sign (CUDA, Float32)", which compared the whole of D bitwise with the host reference; nvcc contracts the
+    Float32 Schur updates to FMA, so the non-perturbed pivots differ by rounding. These checks now compare D within
+    `panel_tol(T)` like every other device-vs-reference test, and the perturbed entries (the subject of the test)
+    still bitwise. CPU after the fix: 59973 pass / 0 fail / 0 broken.
+- Measurements (ubuntu-latest KA CPU backend, Float64, best of 3; device = `factorize!` on the CPU backend, so only
+  launch counts and equality are meaningful, the GPU numbers are owed):
+
+  ```text
+  matrix                         fronts A/B/C  launches  ref_ldlt!   device    ref_solve!  sweep_solve!  same piv
+  random_symindef(2000,0.002)     194/33/23       34      159 ms     292 ms     0.72 ms     1.45 ms      yes
+  kkt(3000,1000,1e-8)             151/80/35       31    12.9 s      13.5 s      6.0 ms      6.8 ms       yes
+  kkt(3000,1000,1e-8), C only       0/0/302       37     2.06 s      7.7 s      5.4 ms      5.8 ms       yes
+  laplacian2d(100,100) - 3I       3305/27/1       20     13.3 ms    31.8 ms     2.1 ms      5.3 ms       yes
+  ```
+
+  The 13 s of the KKT case are the reference's `_best_1x1` fallback scan (`O(w²f)` per front) on fronts with
+  unpaired dual rows, inherited by the device kernel (#75). `@allocated factorize!` (LDLᵀ, second call): 64 B.
+- Deviations from PLAN.md / this task:
+  - **No vendor `sytrf` for regime C.** Regime-B and regime-C groups both run the fused KA `front_ldlt_kernel!`
+    (one workgroup per front), for `"S"` as for `"H"`. LAPACK/cuSOLVER `sytrf` factors `F₁₁` alone with its own
+    Bunch–Kaufman sequence (no `pivot_threshold` against the rows of `F₂₁`, its own tie-breaking, blocked updates),
+    so `D`/`pivot_kind` would not equal the reference, which this task's first test requires, and the "post-check
+    and KA redo" needs a copy of the assembled block (`sytrf` overwrites it). The large-front performance (vendor
+    BLAS-3 for `L₂₁` and the contribution block, or `sytrf` with a documented different sequence) is follow-up #75
+    (T25).
+  - Regime B keeps the panel in global memory (pivot decisions, control words and counts in `@localmem`), the
+    pivot search is serial on work item 1 (the reference's logic verbatim); both are T25 items (#75).
+  - Bunch–Kaufman `α` is now `real(T)(α)` in the reference too (was a Float64 constant, so Float32 comparisons ran
+    in Float64): the device kernels need the same arithmetic, and Metal has no Float64.
+  - `getparam(solver, "inertia")` for SPD/HPD returns `(positive pivots, 0)` and `"npivots"` 0 (not specified).
+  - `Numeric` gains three more fields (`totals`, `psign`, `aux`; PLAN §3.2).
+  - The first regime-A kernel had 34 arguments; its KA CPU launch was not specialized (dynamic dispatch, ~400 B per
+    workgroup). The Hermitian flag is now a type parameter of the pivot parameters and `n` comes from
+    `length(perm)`, so the kernel has 32 arguments and allocates nothing (32 is what was verified, not the exact
+    limit).
+  - Test adaptations required by the new features: `test_api.jl` no longer expects `"S"`/`"H"` (now `"G"`),
+    `"solve_diag"`, `"npivots"`, `"inertia"`, `"pivot_stats"` to raise `NotSupportedError`, and asserts their values
+    instead; `cudss_solver.jl` likewise. No assertion was loosened.
+- Open issues / follow-ups:
+  - #75 (found-by-agent): serial pivot search with the `O(w²f)` fallback scan, one workgroup per regime-C front, no
+    vendor `sytrf`, `F₁₁` not in local memory (T25).
+  - On a GPU, the pivot sequence can differ from the CPU reference if FMA contraction changes the rounding of a
+    near-tie in the pivot tests; the equality tests use well-separated generators, CI will tell.
+  - T16: `"solve_diag"` exists now; the complex `solve_mode` and refinement remain. `ir_n_steps` still warns.
+- Suggested plan changes:
+  - PLAN §2.4 / §3.3: regime C of `"S"` uses the KA in-front kernel (same pivot sequence as the oracle); vendor
+    calls only for the BLAS-3 parts (`L₂₁`, contribution block), or `sytrf` when equality with the reference is
+    not required.
+  - PLAN §3.2: `Numeric` also holds `totals`, `psign`, `aux`.
+  - PLAN §2.7: KA CPU launches with many arguments (34 here) are not specialized and allocate; keep kernels at
+    ≤ 32 arguments.
+
 ---
 
 ## T16 — Iterative refinement, solve sub-phases, `solve_mode`, interrupt, logging   `[ ]`

ext/SparseDirectSolverCUDAExt.jl

@@ -109,6 +109,19 @@ LinearAlgebra.cholesky(A::Symmetric{T, <:CuSparseMatrixCSR{T}}, p::NoPivot = NoP
 LinearAlgebra.cholesky(A::Hermitian{T, <:CuSparseMatrixCSR{T}}, p::NoPivot = NoPivot(); check::Bool = false) where {T} =
     cholesky(CSR(A.data), p; view = A.uplo, check)
 
+"""
+    ldlt(A::CuSparseMatrixCSR; view = 'F', check = false) -> DirectSolver
+    ldlt(Symmetric(A::CuSparseMatrixCSR)) / ldlt(Hermitian(A::CuSparseMatrixCSR))
+
+LDLᵀ/LDLᴴ factorization of `A` on the GPU (≅ CUDSS.jl's `ldlt`); the wrappers
+pass their `uplo` as the view. See `ldlt(::CSR)`.
+"""
+LinearAlgebra.ldlt(A::CuSparseMatrixCSR; view::Char = 'F', check::Bool = false) = ldlt(CSR(A); view, check)
+LinearAlgebra.ldlt(A::Symmetric{T, <:CuSparseMatrixCSR{T}}; check::Bool = false) where {T <: Real} =
+    ldlt(CSR(A.data); view = A.uplo, check)
+LinearAlgebra.ldlt(A::Hermitian{T, <:CuSparseMatrixCSR{T}}; check::Bool = false) where {T} =
+    ldlt(CSR(A.data); view = A.uplo, check)
+
 # ---------------------------------------------------------------------------
 # vendor dense bindings (see `src/dense/vendor.jl` for the contracts)
 

src/SparseDirectSolver.jl

@@ -60,6 +60,7 @@ include("numeric/assembly.jl")
 include("numeric/front.jl")
 include("numeric/subtree.jl")
 include("numeric/factorize.jl")
+include("numeric/ldlt.jl")
 include("numeric/extract.jl")
 
 # solve phase on the device (PLAN §2.5)

src/generic.jl

@@ -1,5 +1,5 @@
 # Julia-native layer (PLAN §1.5, §3.1) on top of the handle layer: `cholesky`,
-# `cholesky!`, `ldiv!`, `\`, `logabsdet`, `diag`, `nnz`. Methods take the
+# `cholesky!`, `ldlt`, `ldlt!`, `ldiv!`, `\`, `logabsdet`, `diag`, `nnz`. Methods take the
 # in-package `CSR`; the CUDA extension adds `CuSparseMatrixCSR` and its
 # `Symmetric`/`Hermitian` wrappers, as CUDSS.jl does. A `SparseMatrixCSC` is not
 # accepted by `cholesky` (that method belongs to CHOLMOD): wrap it with `CSR(A)`.
@@ -33,6 +33,40 @@ function _check_info(solver::DirectSolver, check::Bool)
     return nothing
 end
 
+"""
+    ldlt(A::CSR; view = 'F', check = false) -> DirectSolver
+
+LDLᵀ (real `T`, structure `"S"`) or LDLᴴ (complex `T`, `"H"`) factorization of
+the sparse symmetric/Hermitian matrix `A` on its backend (≅ CUDSS.jl's `ldlt`):
+a [`DirectSolver`](@ref) after `"analysis"` and `"factorization"`
+(synchronized), with in-front Bunch–Kaufman pivoting and static perturbation of
+tiny pivots (read `getparam(solver, "npivots")` and `"inertia"`). `view` as in
+[`cholesky`](@ref). The factorization does not fail (`"info"` stays 0), so
+`check` has no effect; it is accepted for symmetry with `cholesky`.
+"""
+function LinearAlgebra.ldlt(A::CSR{T}; view::Char = 'F', check::Bool = false) where {T}
+    solver = DirectSolver(A, T <: Real ? "S" : "H", view)
+    execute!("analysis", solver, nothing, nothing)
+    execute!("factorization", solver, nothing, nothing; asynchronous = false)
+    _check_info(solver, check)
+    return solver
+end
+
+"""
+    ldlt!(solver::DirectSolver, A; check = false) -> solver
+
+Factorize the new values `A` reusing the analysis of `solver`, as
+[`cholesky!`](@ref) (`"factorization"` the first time after an analysis,
+`"refactorization"` afterwards, synchronized).
+"""
+function LinearAlgebra.ldlt!(solver::DirectSolver, A; check::Bool = false)
+    update!(solver, A)
+    phase = solver.fresh_factorization ? "factorization" : "refactorization"
+    execute!(phase, solver, nothing, nothing; asynchronous = false)
+    _check_info(solver, check)
+    return solver
+end
+
 """
     cholesky!(solver::DirectSolver, A; check = false) -> solver
 
@@ -86,8 +120,8 @@ Base.:\(solver::DirectSolver{T}, B::Union{Array{Complex{T}, 1}, Array{Complex{T}
 """
     diag(solver::DirectSolver)
 
-Diagonal of the factor `L` in factor order (the `"diag"` data parameter), a
-vector on the solver's backend.
+Diagonal of the factor `L` (Cholesky) or of `D` (LDLᵀ/LDLᴴ) in factor order
+(the `"diag"` data parameter), a vector on the solver's backend.
 """
 LinearAlgebra.diag(solver::DirectSolver) = getparam(solver, "diag")
 
@@ -101,17 +135,45 @@ SparseArrays.nnz(solver::DirectSolver) = Int(getparam(solver, "lu_nnz"))
 """
     logabsdet(solver::DirectSolver) -> (log|det A|, sign)
 
-From the diagonal of the Cholesky factor: `log|det A| = 2 Σ log Lₖₖ`, `sign =
-one(T)` (`A` is positive definite). Raises [`FactorizationError`](@ref) when the
+Cholesky: from the diagonal of the factor, `log|det A| = 2 Σ log Lₖₖ`, `sign =
+one(T)` (`A` is positive definite). LDLᵀ/LDLᴴ: `det A = det D` (the
+permutations cancel, `det L = 1`), the product of the 1×1 pivots and of the
+determinants of the 2×2 blocks; `sign` is `±1` (real or Hermitian) or `det D /
+|det D|` (complex symmetric). With perturbed pivots (`"npivots" > 0`) this is
+the determinant of `A + E`. Raises [`FactorizationError`](@ref) when the
 factorization failed (`"info" ≠ 0`).
 """
 function LinearAlgebra.logabsdet(solver::DirectSolver{T}) where {T}
     solver.stage >= STAGE_FACTORIZED && solver.info == 0 ||
         throw(FactorizationError(solver.info, "logabsdet needs a successful factorization"))
+    _is_ldlt_structure(solver.structure) && return _logabsdet_ldlt(solver)
     d = Array(getparam(solver, "diag"))
     return 2 * sum(x -> log(abs(x)), d; init = zero(real(T))), one(T)
 end
 
+function _logabsdet_ldlt(solver::DirectSolver{T}) where {T}
+    N = solver.numeric
+    n = size(solver, 1)
+    d, kind = Array(N.d), Array(N.pivot_kind)
+    herm = !(T <: Complex) || solver.structure == STRUCTURE_HERMITIAN
+    la, sgn = zero(real(T)), one(T)
+    k = 1
+    while k <= n
+        if kind[k] == PIVOT_KIND_2X2_FIRST
+            b = d[n + k]
+            x = d[k] * d[k + 1] - (herm ? conj(b) : b) * b
+            herm && (x = T(real(x)))
+            k += 2
+        else
+            x = d[k]
+            k += 1
+        end
+        la += log(abs(x))
+        sgn *= herm ? T(sign(real(x))) : x / abs(x)
+    end
+    return la, sgn
+end
+
 """
     logdet(solver::DirectSolver)
 

src/numeric/extract.jl

@@ -4,7 +4,8 @@
 """
     host_numeric(numeric) -> Numeric{T, Vector{T}}
 
-Host copy of `numeric` (panels, D, update stack, workspace, statistics, status, pivot order and kinds; the
+Host copy of `numeric` (panels, D, update stack, workspace, statistics and totals, status, pivot order, kinds and sign
+requests, `aux`; the
 plan is shared), usable by [`ref_solve!`](@ref) and [`extract_L`](@ref).
 A host `numeric` is returned as is.
 """
@@ -13,7 +14,8 @@ function host_numeric(N::Numeric{T}) where {T}
     factor, d, stack, work = Array(N.factor), Array(N.d), Array(N.stack), Array(N.work)
     stats, info, piv, pivot_kind = Array(N.stats), Array(N.info), Array(N.piv), Array(N.pivot_kind)
     return Numeric{T, Vector{T}, Vector{Int64}, Vector{Int32}, Vector{Int8}}(factor, d, stack, work, stats, info, piv,
-                                                                             pivot_kind, N.plan)
+                                                                             pivot_kind, Array(N.totals),
+                                                                             Array(N.psign), Array(N.aux), N.plan)
 end
 
 """