T19: general LU ("G"): reference, device kernels in regimes A/B/C, LU solves with A/Aᵀ/Aᴴ, lu/lu!, report

[github-actions]

Closes #19

Task: T19 — opened by the workflow after the implementer run https://github.com/exanauts/SparseDirectSolver.jl/actions/runs/37330630726 ended without a pull request (the App token had expired). The Report block in TASKS.md and the commit message are the implementer's; the body of the PR template is not filled in.

Co-Authored-By: Claude Opus 5.5 noreply@anthropic.com

[claude]

nit (non-blocking): this testset draws rand(T, n) for x0 but has no Random.seed!(666) of its own (every other numeric testset in the file seeds). It is deterministic today only because the preceding testsets seed and run first; add Random.seed!(666) at the top of the body so it also holds when the file is run with a subset of testsets (AGENTS.md test conventions).

[claude]

nit (non-blocking): the signature lines at the top of the forward_sweep!, backward_sweep! and sweep_solve! docstrings (lines 865, 940, 973) still list the old keyword sets; add transpose = false to each so the docstring header matches the method (the prose below already describes it).

[claude]

VERDICT: APPROVE (head 8a49c8c)

Checked

• Scope: only T19 deliverables (reference LU, device LU in regimes A/B/C, LU solves with A/Aᵀ/Aᴴ incl. CSC input, perm_row/perm_col, lu/lu!, uniform batch, ported CUDSS.jl unsymmetric tests). PLAN.md untouched, no Manifest.toml, TASKS.md changed only in the T19 section (status + Report). Public names and parameter/phase strings as in PLAN §1.2–§1.4/§3.1; structure G with view F only.
• Mathematics: A[p, q] = L D U with p = perm[piv], q = perm; the forward sweep gathers by the local row order, the Aᵀ solve uses Uᵀ/Lᵀ with the order restored by the backward sweep; solve_transposed/solve_conjugated give M, Mᵀ, conj(M), conj(Mᵀ) correctly for CSR and CSC input; the refinement residual uses the Mᵀ rows of the G map; logabsdet sign = parity of piv × Π d. Device kernels match the reference step by step (choose / swap over all columns / perturb / rank-1 update of the fully-summed columns and rows / scale / CB update summed in k order): same pivot sequence and summation order, no atomics, no allocation or host sync, @LocalMem sized by Val (two buffers of (LB − reserve) ÷ (2 sizeof(T)), consistent with schedule_elsize). Signed amap offsets, _group_amap run ordering, extend-add of the Uᵀ packed blocks and extract_lu's ipiv[rows[i]] checked.
• Conventions: generic T/INT, host symbolic in Int, dense ops only via the interface, no @allowscalar, errors from src/errors.jl, docstrings on the public API, backend code only in ext/.
• Tests: test_numeric_lu.jl covers reference, device vs reference under three regime mixes, regime-A Int64 maps, perturbation, determinism/allocation budgets, sweeps one at a time against extract_lu, uniform batch, and the public API over INTTYPES; the ported unsymmetric tests mirror CUDSS.jl's (pivot in C/R/N, lu/lu!/ldiv!/backslash, uniform_batch_lu). Loops over BACKENDS/ELTYPES, tol(T), shared helpers; no test weakened, no @test_broken added.
• Report: matches the diff; the regime-C deviation (fused KA kernel instead of vendor getrf + laswp) is the same as T15's accepted deviation and is documented with the T25/#75 follow-up; #84 (found-by-agent) exists and is triaged. CPU CI and the CUDA job are still pending (not a blocker; the pipeline handles them).

Blocking findings: none.

Nits (non-blocking, inline): missing Random.seed!(666) in the perturbation testset; docstring signature lines of forward_sweep!/backward_sweep!/sweep_solve! omit the new transpose keyword.

[claude]

nit (non-blocking): this testset draws x0 = rand(T, n) (line 112) but does not call Random.seed!(666) first, unlike the other testsets in this file. The checks do not depend on x0 so nothing breaks, but AGENTS.md asks every numeric test to seed. (The LDLᵀ perturbation testset has the same omission; fine to leave for a follow-up.)

Suggested change

	n, j = 60, 23
	n, j = 60, 23
	Random.seed!(666)

[claude]

VERDICT: APPROVE (head bd0d8ab)

Reviewed T19 (general LU) against the task text on main, AGENTS.md and PLAN.md §2.3/§2.4/§2.5/§3.3/§3.4.

Scope. Only the T19 deliverables: src/reference/lu.jl, src/numeric/lu.jl, the Uᵀ storage (ufactor/ustack), "G" assembly maps, LU solve kinds in the sweeps (A, Aᵀ, Aᴴ incl. regime-C dense path and CSC input), the Mᵀ rows of the refinement map, perm_row/perm_col/diag/npivots/inertia, lu/lu!/logabsdet, the CUDA lu(::CuSparseMatrixCSR). PLAN.md untouched; TASKS.md changed in a single hunk inside the T19 section; no Manifest.toml; public names, parameter and phase strings as in PLAN §1.2–§1.4/§3.1. Nothing from T20+.

Correctness (read in detail). Reference ref_lu! / _lu_front!: assembly of upper entries via negative offsets, in-block threshold row pivoting (γ over the whole front column, candidates restricted to fully-summed rows), GESP perturbation, per-step update of fully-summed rows/columns and one F₂₂ −= (L₂₁D)U₁₂ with the same summation order as the kernels; extract_lu row placement ipiv[rows[i]] and A[p,q] = LDU, ref_solve_lu! for A and Aᵀ all check out. Device kernels: _lu_get/_lu_set! split of F over the two lower structures, swap over all f columns with distinct storage per column (no races), update index ranges match the reference loops, scale/CB update, regime-A packed layout offsets (_lt_local_cb_offset, _subtree_scatter_lu!) consistent with the LDLᵀ helpers, local memory 2·NE·sizeof(T) + ctl + pv ≤ 48 KiB with the Cholesky reserve (24 control words). Sweeps: _SV_LU (L reordered fwd, U bwd) and _SV_LU_T (Uᵀ fwd, Lᵀ bwd then order restored) are correct with the kernels' LDL/CONJ flags, tr = 'T' on the dense path, diagonal kernel safe (no 2×2 kinds). solve_transposed/solve_conjugated give the right (Mᵀ, conj) combination for CSR and CSC input in all three modes; _grow_refinement keeps n and the 2n+1 row pointers. logabsdet sign = parity of piv × Π sign(d) is right since the symmetric reordering cancels. No atomics in assembly/extend-add, no allocation or host sync in the numeric phase, generic T/INT, errors from src/errors.jl.

Tests. test_numeric_lu.jl covers reference, device-vs-reference on three regime mixes (incl. "C only"), Int64 maps, perturbation (default and regime-C root, pivot_epsilon, scaled), determinism and allocation budgets, sweeps one at a time against extract_lu and every dense impl, uniform batch with ubatch_index, and the public API over INTTYPES (solve_mode 0/1/2 with sub-phases, IR and FGMRES-IR, CSC input, view refusal, lu/lu!). Ported cudss_execution/cudss_generic/cudss_solver/cudss_uniform_batch LU parts match the CUDSS.jl tests. Helpers live in test/utils.jl/test/matrices.jl, loops over BACKENDS/ELTYPES, tolerances from tol(T)/panel_tol(T) (the 1e-6 on the singular perturbation case mirrors the LDLᵀ test). test_api.jl/test_symbolic_schedule.jl adaptations strengthen rather than weaken. No @test_broken added.

Report. Matches the diff; deviations (KA kernel in regime C instead of vendor getrf/laswp, Uᵀ in the L layout instead of U panels, view 'F' only) are stated with rationale and plan-change suggestions; follow-ups have issues (#84 found-by-agent, triaged; #75/T25). CI (Run tests CPU, CUDA, Aqua) is still pending on this head; the Report claims CPU green only and defers CUDA to CI, so no contradiction.

Blocking findings: none.

Non-blocking: one inline nit (missing Random.seed!(666) in the perturbation testset).

[github-actions]

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

[claude]

VERDICT: APPROVE (head 99b4f07)

Reviewed against the T19 text in TASKS.md on main, AGENTS.md and PLAN.md §2.3–§2.5, §3.1–§3.5.

Scope. Only T19 deliverables: reference LU, device kernels (A and B/C), LU solves with A/Aᵀ/Aᴴ, refinement with Mᵀ, lu/lu!, perm_row/perm_col, ported unsymmetric and uniform-batch LU tests. PLAN.md untouched; TASKS.md changed only in the T19 section (marker, Report); no Manifest.toml; CUDSS parameter/phase strings unchanged. Follow-up #84 exists (found-by-agent, triaged).

Correctness (read line by line). The L/Uᵀ split with _lu_get/_lu_set! keeps every access in the lower triangle of a packed or panel front; the pivot step (choose on item 1, full-row swap, perturb, rank-1 update of columns k+1:w × rows k+1:f and rows k+1:w × columns w+1:f, scale) matches _lu_front! of the reference term by term, including operation order. The CB update writes each (ii, jj) once (no atomics, deterministic) and for regime A the packed m×m offset coincides with the trailing columns of the f-triangle. Regime-A local memory: (LB − 256) ÷ (2 sizeof(T)) entries per buffer matches subtree_capacity with elsize = 2 sizeof(T); ctl + pv fit in the 256-byte reserve. Solve: _SV_LU (L forward with the pivot reorder, U backward) and _SV_LU_T (Uᵀ forward, Lᵀ backward with the order restored) are consistent with extract_lu (A[perm[piv], perm] = L D U, below-block L rows at ipiv[rows[i]]); both panels store a unit diagonal, so the non-unit TRSV/trsm are exact. solve_transposed/solve_conjugated give the right (Mᵀ, conj) pair for CSR and CSC inputs and modes 0/1/2; the "G" refinement map rows n+1:2n are the rows of Mᵀ, _grow_refinement keeps n from R. _logabsdet_lu sign = parity of the local row order × Π sign(d). INT conversion of the 2n+1 row pointers goes through the overflow-checked path. No @allowscalar, no allocation or host sync inside the numeric/solve phases, errors from src/errors.jl.

Tests. test_numeric_lu.jl loops over BACKENDS and ELTYPES, seeds 666, uses tol(T)/panel_tol(T) (the 1e-6 on the perturbed singular case has the T15 precedent); helpers added in test/utils.jl and test/matrices.jl. Ported cudss_execution, cudss_generic, cudss_solver, cudss_uniform_batch cover "G". No weakened assertions, no new @test_broken.

Report. Matches the diff (storage, 32 kernel arguments, doubled memory-estimate slots, perm_row of member 1, deviations). CI is still pending on this head; the pipeline handles it.

Non-blocking notes for the owner.

• Deviation to confirm: regime C uses the fused KA kernel instead of vendor getrf/laswp + post-check (same rationale and precedent as T15; follow-up T25/#75). The pivot search is serial on one work item per front.
• test_numeric_lu is not in SPLIT_FILES (≈9.5 min cold standalone); splitting per element type would help CI wall time.
• pivot_sign is silently ignored for "G" (as for Cholesky); fine, but worth one line in the pivot_sign docs.
• The PR body is not from the template (the workflow opened the PR after the App token expired); Closes #19 is present.

[github-actions]

CI-FIX round 2 of 3: https://github.com/exanauts/SparseDirectSolver.jl/actions/runs/37525343332 failed on 99b4f07; Claude is diagnosing it.

[claude]

nit (non-blocking): this follow-up has no issue. getparam(solver, "perm_row") on a uniform batch silently returns member 1's row order (src/solver.jl, the "perm_row" branch of getparam), while every other per-member output ("info", "npivots", "inertia", "pivot_stats", "diag") is returned per member. Either open a found-by-agent issue for a per-member "perm_row" (e.g. an n × nbatch matrix, as "diag" concatenates members) and cite it here, or state in the getparam docstring that it is by design. Not blocking: the docstring table already says "batch member 1".

[claude]

nit (non-blocking): the Report measures this file at ≈ 9.5 min standalone, "mostly kernel compilation", yet it is not in SPLIT_FILES of test/runtests.jl. Every testset here already loops over ELTYPES (none needs RUN_SHARED), so adding "test_numeric_lu" to SPLIT_FILES would run it as four one-eltype parts on separate workers, each compiling only its element type's LU kernels, like test_numeric_ldlt. Worth doing in a follow-up if it shows up on the CI critical path.

[claude]

VERDICT: APPROVE (head 73a0063)

Reviewed T19 (general LU) at head 73a0063 against the T19 text in TASKS.md on main, AGENTS.md and PLAN.md §2.2–§2.5, §3.1–§3.5. This head adds, on top of the previously approved 99b4f07, the owner's growth-aware panel tolerance (growth_tol, factor_growth in test/utils.jl, the AGENTS.md test-convention note) and the updated CPU/CUDA counts in the Report.

Scope. Only T19 deliverables: src/reference/lu.jl, src/numeric/lu.jl, the Uᵀ storage (ufactor/ustack), "G" assembly maps with the signed-offset grouping, schedule_elsize, LU solve kinds in the sweeps (A, Aᵀ, Aᴴ, regime-C dense path, CSC input), the Mᵀ rows of the refinement map for IR and FGMRES-IR, perm_row/perm_col/diag/npivots/inertia/pivot_stats, lu/lu!/logabsdet, the CUDA lu(::CuSparseMatrixCSR), ported CUDSS.jl unsymmetric, generic and uniform-batch LU tests. PLAN.md untouched; TASKS.md changed only in the T19 section; no Manifest.toml; cuDSS parameter and phase strings unchanged; "G" keeps the view-F rule of T05 with #84 (triaged) for the cuDSS difference.

Mathematics. Reference _lu_front!: row interchanges restricted to the fully-summed rows, threshold test against the whole remaining front column, static ±ε perturbation, F₂₂ -= (L₂₁ D) U₁₂ summed in the device order; extract_lu places the rows of L below a diagonal block at ipiv[row] and leaves U columns unpermuted, so A[perm[piv], perm] = L D U and both ref_solve_lu! branches are consistent with it. Device kernels: _lu_choose!/_lu_swap!/_lu_pivot!/_lu_update!/_lu_scale!/_lu_cb_update! reproduce the reference step by step (same expressions, same summation order, pivot order swapped with the rows); every access to the L and Uᵀ structures stays in columns ≤ w; one writer per destination, no atomics, no allocation or host sync in the phase; @localmem sized by Val, the regime-A budget halves NE to match schedule_elsize = 2 sizeof(T). Sweeps: A x = b is L forward with the local row order, diagonal, U = (Uᵀ)ᵀ backward without reordering; Aᵀ x = b is Uᵀ forward, diagonal, Lᵀ backward restoring the order, which is exactly what _sv_piv_apply! (gather by piv forward, scatter to piv backward) and the assembled-numbering gather of the rows below require; Aᴴ goes through solve_conjugated. solve_transposed is right for CSR and CSC input. _logabsdet_lu uses the cycle parity of piv (the symmetric P cancels). _grow_refinement keeps n from W.R, so the 2n-row map is safe.

Conventions. Generic T/INT; host symbolic in Int with the existing overflow-checked cast; errors from src/errors.jl; dense ops through the interface with the fac buffer; no @allowscalar (the piv/d reads in getparam/logabsdet are explicit phase-boundary copies); backend code only in ext/; docstrings on every new public function.

Tests. test/test_numeric_lu.jl covers everything the task lists (reference, device vs reference in three regime mixes with exact pivot sequence, solves A/Aᵀ with nrhs 1 and 5, atomic/deterministic, Int64 maps, perturbation incl. scaled ε, determinism and allocation budgets, sweeps one at a time and every dense impl, uniform batch incl. ubatch_index, public API for both INT with solve_mode 0/1/2, sub-phases, IR/FGMRES-IR, logabsdet, CSC input, lu/lu!); the ported tests match CUDSS.jl modulo the documented generator changes; random_general residuals ≤ tol(T). Loops over BACKENDS/ELTYPES/INTTYPES, Random.seed!(666), tolerances from tol/panel_tol/growth_tol; the pivot sequence is still compared exactly, so growth_tol does not weaken the oracle. test_api.jl and test_symbolic_schedule.jl adaptations follow from the feature; no @test_broken added.

Report. Matches the diff: status [!] with the documented deviations (KA kernel in regime C instead of vendor getrf/laswp, Uᵀ in a second L-layout buffer, pivot_type semantics, view F only), follow-ups #84/#75/#86 cited, counts plausible (CPU full run, test_numeric_lu 1816 CPU / 1788 CUDA on the owner's machine). CI on this head: build green, CPU and CUDA test jobs still pending; the pipeline handles them.

Blocking findings: none.

Non-blocking nits (inline): the uniform-batch perm_row (member 1 only) follow-up has no issue; test_numeric_lu could join SPLIT_FILES.