[michel2323]

Refs #194

Task: not a TASKS.md task; follow-up of #193 on top of #192. Stacked on PR #192 (base maint/kkt-solver-0.4.0). The branch will be rebased once #192 is updated or merged: #192 gets one more commit with review fixes that is not pushed yet.

What was built:

• Static pivot replacement in the fused L D Lᵀ (FusedSolver, structure "S"): pivot_epsilon and pivot_sign are applied at all three fused pivot sites (regime-A subtrees, narrow panels, in-kernel wide panels), with the same rule as the stock no-pivot kernel. Replaced pivots are counted in "npivots" and take the requested sign in "inertia". The default pivot_epsilon is now the stock one; it was 0 before, and setting pivot_epsilon = 0 restores the old behaviour. A NaN or Inf pivot still fails the factorization.
• Matching and scaling on fused "S" (matching_alg, including "algo5" and "algo6"): the factorization runs on scaled values, the 2×2 pivot pairs from the matching shape the ordering, and the solve undoes the scaling inside the permutation kernels, so no extra launch. "perm_matching", "scale_row" and "scale_col" are forwarded.
• KKTSolver stays on the fused pipeline when matching or pivot replacement is requested. 0.4.0: KKTSolver front end, MadNLP binding moves to MadNLP #192 sent these cases to the stock pipeline. Now only pivoting = :bk and backends without forward progress use the stock pipeline.
• ir_n_steps on the fused solve: plain refinement against the caller's unscaled values, with no host synchronization, and captured by the solve graph. It is also a new KKTSolver keyword that works on both pipelines. On fused, ir_tol > 0 and ir_mode = "fgmres" raise NotSupportedError.
• Docs (lib/fused.md, lib/kkt.md, man/madnlp.md, man/expectations.md), docstrings and CHANGELOG are updated.

Evidence (measured with CUDA on the owner's machine; timings are single wall-clock runs, not Slurm benchmarks)

60 recorded KKT matrices from a CO2-capture flowsheet solve. Settings: static pivoting + algo5, ε = 1e-10, pivot_sign. Every row had 0 wrong inertias and 0 failures.

refinement steps	fused median / max residual	stock median / max residual	fused fact+solve	stock fact+solve
0	8.6e-07 / 5.9e-06	1.1e-06 / 1.2e-05	1.23 ms	5.82 ms
1	2.1e-12 / 2.4e-10	2.1e-12 / 1.3e-09	1.51 ms	6.19 ms
2	6.2e-13 / 2.1e-11	5.7e-13 / 2.7e-11	1.79 ms	6.45 ms

Live MadNLP design solve of the flowsheet. Preset: ShiftedPreset, with the options patched in memory (MadNLP's binding does not pass ir_n_steps yet).

configuration	pipeline	iterations	termination	factorizations	linear solver time
cuDSS default (matching + its pivoting)	—	88	optimal 1e-8	139	0.58 s
SDS static + algo6, 2 refinement steps	fused	104	optimal 1e-8	208	0.86 s
SDS static + algo6, 1 refinement step	fused	104	optimal 1e-8	208	1.42 s
SDS static + algo6, 2 refinement steps	stock	104	optimal 1e-8	208	2.6 s
SDS static + algo6, 1 refinement step	stock	274	optimal 1e-8	452	8.4 s
SDS Bunch–Kaufman + algo6 (preset of MadNLP PR #728)	stock	108	optimal 1e-8	153	3.9 s
SDS static + algo6, no refinement	fused	1178	search direction too small (failed)	2525	115.8 s
SDS static + algo6, no refinement	stock	3000	iteration limit (failed)	3727	293.8 s
SDS static + algo6, no refinement, analysis+matching redone every factorization	stock	90	optimal 1e-8	166	16.7 s
cuDSS, matching, pivot_type 'N', no refinement	—	3000	iteration limit (failed)	3865	127.8 s

What the live solve shows:

• Static pivoting without refinement fails in both SDS pipelines and in cuDSS.
• The cause is the matching going stale: recomputing it at every factorization gives 90 iterations.
• Two refinement steps make up for the stale matching. With them, the fused pipeline converges in 104 iterations and spends 0.86 s in the linear solver, against 2.6 s on the stock pipeline and 3.9 s for the Bunch–Kaufman preset.

Open points:

• One refinement step gives 104 iterations on fused but 274 on stock. This is not explained yet.
• A gap to cuDSS remains in the factorization count (208 vs 139).
• MadNLP's binding and preset change (an sds_ir option, static + algo6 + 2 steps) is a follow-up on MadNLP PR #728, once this PR is green.

Tests: new testsets in test/test_fused.jl and test/test_kkt.jl cover:

• a replaced zero pivot with both signs and two ε values;
• pivot replacement at every kernel site, checked against the stock factor (pivot kinds exact, D within growth_tol);
• matching algo5/algo6 with refactorization, bitwise reproducibility, graph replay and allocation budget;
• ir_n_steps with matching and pivot replacement, in place and out of place, against the stock pipeline, plus the graph replay of the corrections;
• KKTSolver with matching, pivot_epsilon, pivot_sign, deferred analysis and ir_n_steps on every pipeline the backend runs.

Local status: only the new and changed testsets were smoke-run, on CPU and CUDA for Float32 and Float64. The full suite (all four element types, CPU, CUDA, AMDGPU, Aurora) is left to CI.

Two tests assert that the unrefined residual on kkt_slack_matrix(T, 200, 100, 0) under algo5 is above tol(T), which shows the test is not vacuous: without refinement the solve really is inaccurate. Margin measured over 5 seeds, on the CPU and CUDA backends and on both pipelines (always 50 replaced pivots):

• Float32: r0 is 10–30× tol(T);
• Float64: r0 is 3e4–1e5× tol(T);
• after one step: about 1e-6 (Float32) and 1e-15 (Float64).

The inaccuracy comes from the ε perturbation of the 50 replaced pivots, not from rounding, so I do not expect AMDGPU or oneAPI to close the 10× Float32 margin. I have not run it there.

Deviations from PLAN.md / TASKS.md:

• No PLAN.md deviation. Behaviour changes of the released FusedSolver:
  • The fused "S" default pivot_epsilon changed from 0 to the stock default (listed under "Changed" in CHANGELOG).
  • wides_potrf = true (an internal, non-default plan option) refuses a nonzero LDL pivot_epsilon with NotSupportedError.

Follow-up issues opened: none (MadNLP binding: follow-up on MadNLP PR #728).

🤖 Generated with Claude Code

https://claude.ai/code/session_01MsLzHkzHJjM2g5P6ZyZKhr

[claude]

VERDICT: APPROVE (head c11f803)

Not a TASKS.md task (follow-up of #193/#194 on top of the merged #192), so the task-section rules do not apply; reviewed against the issue, AGENTS.md and PLAN.md §1.2–§1.4/§3.1.

Checked

• Scope: src/fused.jl, src/kkt.jl, tests, docs, CHANGELOG only. PLAN.md, TASKS.md and Manifest.toml untouched. Parameter strings (pivot_epsilon, pivot_sign, matching_alg, ir_n_steps, perm_matching, scale_row, scale_col) and the phase strings are the cuDSS ones.
• Static pivot: _ldl_static_pivot is the stock _lt_select! rule with pivot_type = N (requested sign, else the sign of d, + for 0), differing only on NaN/Inf, which the docstring, lib/kkt.md and the CHANGELOG state. pivot_sign is read in factor-column order through S.perm, as the stock kernel does (psign[perm[g]]); the symmetric matching keeps cperm = S.perm, so no second permutation is missing. The pure function keeps the factorization bitwise deterministic.
• All three fused pivot sites write d, pivot_kind and the replaced pivot where the later reads expect it: the subtree keeps D in d (unit diagonal in buf), the narrow panel writes the diagonal slot back in _front_ldl_scale! after the barrier (the loop only touches i, k > j, as the subtree loop does), the in-kernel wide panel writes it under li == 1 before the barrier. npivots/inertia are device counts at the phase boundary, no @allowscalar.
• Matching: _factor_values rescales into the same st.values buffer at every (re)factorization, so the captured factorization graph stays valid (F.gnz === F.nzval); the solve passes the solve_io permutations and scalings into the existing permutation kernels, no extra launch. _fused_pairs! applies the KKTSolver rule (pairs with jobs 5/6 only) at every analysis.
• Refinement: residual from vec(solver.A.nzval) (unscaled, the user values), corrections through fused = fs, ir_tol = 0 and ir_mode = "ir" enforced at analysis and at every refining solve (no host sync), workspace allocated at analysis or by setparam! outside the solve phase, aliasing X === B handled by the W.Bc copy, ir_steps semantics as on DirectSolver.
• KKTSolver: only :bk and backends without forward progress go to the stock pipeline; pivot_sign, pivot_epsilon, ir_n_steps set on either inner solver; negative ir_n_steps rejected by the option parser.
• Conventions: generic T/INT (eps cast to real(T) at launch, Int8 signs), @Const on psg, no atomics, no allocation inside the numeric/solve phases, errors from src/errors.jl, no backend code outside ext/.
• Tests: new testsets loop over BACKENDS and eltypes_among((Float32, Float64)) (the fused pipeline is real-only), Random.seed!(666), tol/growth_tol/weak_pivot_tol from the helpers; the removed NotSupportedError assertions are replaced by positive tests of the same options; the site-coverage test checks all four kernel sites; graph replay and allocation budgets covered. One nit inline (a hard-coded 1e-6 tolerance).
• PR body vs diff: every behaviour claim matches the code; the body states honestly that only the new and changed testsets were run locally (CPU and CUDA, Float32/Float64) and leaves the full suite to CI.
• CI: docs built; the CPU, CUDA, AMDGPU and Aurora legs are still pending at review time, so nothing contradicts the body. Note for the owner: two tests assert the unrefined residual is above tol(T) (kkt_slack_matrix(T, 200, 100, 0) under algo5); the body measures a 10–30× margin for Float32 on CPU/CUDA, unmeasured on AMDGPU/oneAPI.

Blocking findings: none.