Repository navigation
T16: iterative refinement, solve_mode, interrupt, logging - #79
Conversation
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
claude
Bot
added
the
claude:pr
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
![claude[bot]](https://avatars.githubusercontent.com/in/1236702?s=60&v=4){kind=small}
| r = @index(Group, Linear) | ||
| nr = @localmem RT (WG,) | ||
| nb = @localmem RT (WG,) | ||
| nrhs = size(R, 2) |
There was a problem hiding this comment.
Blocking: wrong nrhs for a row-major strided-vector B. nrhs = size(R, 2) is the capacity of the refinement workspace (max(nrhs, max_rhs(solver.workspace)) at allocation, or the grown size), not the number of right-hand sides of this solve. It is passed to _rhs_get(B, i, r, n, nrhs, Val(BT)), whose vector/Val(true) branch indexes B[r + (i - 1) * nrhs]. That layout is reachable: MatrixDescriptor(T, n, k; transposed = true) + update!(desc, v::AbstractVector) is accepted (update! skips the shape check for 1-D data) and _rhs_data returns (v, true). With ir_tol > 0, after any earlier solve with more right-hand sides (or any _grow_refinement), the kernel reads the wrong entries of B and, under @inbounds, past its end: wrong ‖Bₖ‖ → wrong early exit, or a crash.
residual!, copy_rhs! and add_correction! already receive the true nrhs; do the same here: add an nrhs argument to _residual_norms_kernel! (pass Int(nrhs) from residual_norms!) and drop the size(R, 2) line. A test with a transposed descriptor over a strided vector, ir_tol > 0, after a wider solve, would catch it.
| catch err | ||
| if err isa InterruptedError | ||
| # the panels are partly overwritten: back to "analyzed", a "factorization" must follow | ||
| solver.stage = STAGE_ANALYZED |
There was a problem hiding this comment.
nit (non-blocking): after an interrupted "refactorization" the solver is back at STAGE_ANALYZED but fresh_factorization stays false, so the LinearAlgebra layer's cholesky!/ldlt! (which pick "refactorization" when !fresh_factorization) raise FactorizationError instead of recovering. Setting solver.fresh_factorization = true here makes "clear the flag and factorize again" work from both layers.
| steps = 0 | ||
| for _ in 1:nsteps | ||
| _poll_interrupt(interrupt) | ||
| residual!(W, nzval, X, B; nrhs, transposed, b_transposed, conjugate) | ||
| if tol > 0 | ||
| norms = residual_norms!(W, B; nrhs, transposed = b_transposed) | ||
| rel = _max_relative_residual(norms, nrhs) | ||
| _log(LOG_DEBUG, () -> "refinement: step $steps, relative residual $rel") | ||
| rel <= tol && break | ||
| end | ||
| permute_rhs!(ws.Y, W.R, S.perm; conjugate) | ||
| forward_sweep!(ws, S, N; nrhs, deterministic) | ||
| diagonal_sweep!(ws, S, N; nrhs) | ||
| backward_sweep!(ws, S, N; nrhs) | ||
| add_correction!(X, ws.Y, S.perm; nrhs, transposed, conjugate) | ||
| steps += 1 |
There was a problem hiding this comment.
nit (non-blocking): the Report lists "an interrupted refinement reports ir_n_steps = 0" as an open issue, but it is a two-line fix inside this task's scope: have refine! report the completed steps on interrupt (e.g. catch err; err isa InterruptedError && rethrow(InterruptedError(...)) is overkill; simplest is to pass a Ref/callback or wrap the loop body in try … finally in _refine_phase! so solver.ir_steps is set from a counter the loop updates). Then the data parameter matches "steps performed" in every exit path.
There was a problem hiding this comment.
VERDICT: CHANGES_REQUESTED (head bb32e34)
Checked: scope (only T16 deliverables; PLAN.md untouched; TASKS.md changed only in the T16 section; no Manifest.toml; phase and parameter strings unchanged), the four listed tests in test/test_refinement.jl (all present; the one @test_broken is documented in the test and the Report with a sound argument — symmetric scaling cannot defeat a backward-stable Cholesky — and the 100× criterion is asserted on a perturbed-pivot KKT factorization instead), shared helpers / Random.seed!(666) via runtests / tol(T), the conventions (generic T/INT, host map in Int with an overflow check before the cast, 1-D workgroups, @localmem sized by Val, gather-based residual without atomics, no @allowscalar, errors from src/errors.jl), the mathematics of solve_conjugated for "S"/"H"/"SPD"/"HPD" with CSR and CSC input, the conjugation flow residual → permute → sweeps → add_correction!, the aliasing handling of "solve"/"solve_refinement", the interrupt recovery path, and the Report against the diff (consistent; CI pending: cpu, cuda, review; Aqua build green).
Blocking finding (inline comment):
src/solve/refinement.jl:131—_residual_norms_kernel!usesnrhs = size(R, 2)(workspace capacity) in_rhs_getfor a row-major strided-vectorB. That layout is reachable throughMatrixDescriptor(T, n, k; transposed = true)+update!(desc, vector). Withir_tol > 0and a workspace wider than the currentnrhs(earlier solve with more RHS, or_grow_refinement), the‖Bₖ‖reduction reads the wrong entries and out of bounds under@inbounds: wrong early-exit decision or a crash. Pass the truenrhsas a kernel argument as the other three kernels do.
Non-blocking nits (inline): fresh_factorization is not reset after an interrupted refactorization (LinearAlgebra cholesky!/ldlt! then cannot recover); ir_n_steps after an interrupted refinement reports 0 (listed as an open issue, but fixable in scope).
Accepted deviations: lazy allocation of the refinement workspace at the first refining solve (documented; keeps the MadNLP path allocation-free), ir_tol = 0 in the LinearAlgebra layer (early exit only when the user sets ir_tol), interrupt polled between launch groups.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
Review round 1 addressed in ca7e87e. Blocking
Nits (both fixed)
Rejected: none. Tests (CPU): 🤖 Generated with Claude Code |
|
CI-FIX round 1 of 3: https://github.com/exanauts/SparseDirectSolver.jl/actions/runs/37053617713 failed on bb32e34; Claude is diagnosing it. |
| _log(LOG_DEBUG, () -> "refinement: step $steps, relative residual $rel") | ||
| rel <= tol && break | ||
| end | ||
| permute_rhs!(ws.Y, W.R, S.perm; conjugate) |
There was a problem hiding this comment.
nit (non-blocking): permute_rhs!(ws.Y, W.R, …) derives nrhs from size(W.R, 2), i.e. the workspace capacity, not the nrhs of this solve. When the refinement workspace holds more columns than the current solve (e.g. after a 3-rhs solve, a 1-rhs refinement), the stale columns of W.R are permuted into ws.Y[:, nrhs+1:end] on every step. Correct (the sweeps only touch 1:nrhs, and max_rhs(W) ≤ max_rhs(ws) holds by construction), but it is wasted launches/bandwidth in the IPM loop; permuting view(W.R, :, 1:nrhs) or giving permute_rhs! an explicit nrhs keyword fixes it.
There was a problem hiding this comment.
VERDICT: APPROVE (head ca7e87e)
Checked against TASKS.md T16, AGENTS.md and PLAN.md §1.2–§1.4, §2.5, §3.1 (git diff origin/main...HEAD).
Scope: only the task's deliverables; PLAN.md untouched; TASKS.md changed only in the T16 section (status [!], Report); no Manifest.toml; public names, parameter strings and phase strings unchanged ("solve_refinement", ir_n_steps, ir_tol, solve_mode, user_host_interrupt). Lifting the T13 complex-Hermitian-CSC restriction is in scope (it needed the conjugated solve of solve_mode).
Tests: all four listed tests exist in test/test_refinement.jl; they use the shared helpers (relres, tol(T), thrown, api_matrix, kkt_matrix, random_spd/random_symindef, BACKENDS/ELTYPES/INTTYPES), seeding is done per file in runtests.jl. The single @test_broken (100× in one step on a scaled SPD matrix) is documented in the Report with a sound reason (Cholesky's normwise backward error is invariant under symmetric diagonal scaling, so the unrefined relres is already at rounding level) and the 100× criterion is asserted instead on a KKT factorization with static pivot perturbations, which is the case refinement exists for (issue #71). test_api.jl changes replace the T13–T15 "not yet" assertions with assertions of the results, nothing weakened.
Conventions: host map in Int with an overflow check before casting to INT (_to_index_type); kernels generic in T/INT, 1-D, @localmem sized by Val, no atomics in the residual (gather SpMV) or correction; no @allowscalar; errors from src/errors.jl; docstrings on the new public functions. The host read in residual_norms! only runs when ir_tol > 0, and the refinement storage is allocated by the first refining solve (both documented deviations; the latter mirrors the existing _workspace! growth).
Correctness: solve_conjugated checked for all 12 (structure × input orientation × solve_mode) cases: op(A) is always M or conj(M), solved as conj(M⁻¹ conj(b)) with conjugating permutations; the residual kernel's xor(s < 0, CJ) conjugation is consistent with FullPatternMap.conjflag; add_correction! is conflict-free; X === B in "solve" is handled by copying B before the sweeps on the same queue; interrupted factorization/refactorization resets to "analyzed" with fresh_factorization = true; interrupted refinement reports the completed steps (round-1 findings all addressed).
Report: matches the diff; counts plausible (517 pass / 1 broken for the file). CI: build green, CPU and CUDA test jobs still pending at review time (handled by the pipeline).
Blocking findings: none.
Non-blocking: one inline nit (refine! permutes the capacity columns of W.R rather than nrhs). The PR body's test counts (453) are stale relative to the Report (517). PLAN §3.1 says the LinearAlgebra layer refines "with early exit on ir_tol"; the Report keeps ir_tol = 0 there to avoid a host sync per step — documented, for the owner to confirm.
|
CI-FIX round 2 of 3: https://github.com/exanauts/SparseDirectSolver.jl/actions/runs/37056988218 failed on ca7e87e; Claude is diagnosing it. |
Bitwise comparisons between solves in test_refinement.jl use deterministic_mode = 1: the atomic forward sweep sums in a run-dependent order on GPUs. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
|
CI-FIX round 2 (703a467) Root cause: on the CUDA runner (H200), 5 bitwise Fix: the two testsets that compare whole solves bitwise now set Tests: |
There was a problem hiding this comment.
VERDICT: APPROVE (head 703a467)
Reviewed git diff origin/main...HEAD against TASKS.md T16, AGENTS.md and PLAN.md §1.2–§1.4, §2.5, §3.1. The delta since the approved head ca7e87e is the CI-fix round 2: tests only (deterministic_mode = 1 in the two testsets that compare whole solves bitwise) plus one Report line; no library code changed.
Scope: only T16 deliverables (src/solve/refinement.jl, src/logging.jl, solve_mode conjugation in permute.jl/solver.jl, interrupt polls in factorize.jl/ldlt.jl, LinearAlgebra-layer ir_n_steps = 2, bench/refinement.jl); PLAN.md untouched; TASKS.md changed only in the T16 section (marker [!], Report); no Manifest.toml; phase and parameter strings unchanged. Lifting the T13 complex-Hermitian-CSC restriction is a consequence of the conjugated solve and in scope.
Tests: the four listed tests exist in test/test_refinement.jl and use the shared helpers (relres, tol(T), thrown, api_matrix, kkt_matrix, random_spd/random_symindef, new badly_scaled_spd in test/matrices.jl), BACKENDS/ELTYPES/INTTYPES, seeding via runtests.jl; the file is picked up automatically by TEST_FILES. The single @test_broken (100× in one step on a symmetrically scaled SPD matrix) is documented in the test and the Report with a sound argument, and the 100× criterion is asserted on a perturbed-pivot LDLᵀ factorization instead. The round-2 change does not weaken the "composes exactly" test: the assertion is still ==, bitwise equality between two solves is only defined on the atomic-free path (PLAN §2.5, deterministic_mode = 1), and default-mode compositions remain bitwise-compared on smaller systems in test_api.jl and test_numeric_ldlt.jl. test_api.jl changes replace "not yet" assertions with assertions of the results.
Conventions: host refinement map in Int with an overflow check before casting to INT; kernels generic in T/INT, 1-D, @localmem sized by Val, gather SpMV and conflict-free correction (no atomics), no @allowscalar, errors from src/errors.jl, docstrings on new public functions. The host read of residual_norms! only runs when ir_tol > 0; the refinement storage is allocated by the first refining solve (documented deviation, mirrors _workspace!).
Correctness: solve_conjugated re-derived for all structure × CSR/CSC × solve_mode cases (op(A) is M or conj(M)); residual conjugation xor(s < 0, CJ) consistent with FullPatternMap.conjflag; the conjugated correction conj(M⁻¹ conj(R)) flows through permute_rhs!/sweeps/add_correction! with the right flags; _residual_norms_kernel! now takes the true nrhs (round-1 blocker fixed, regression test present); X === B in "solve" copies B before the sweeps on the same queue; interrupted factorization resets to analyzed with fresh_factorization = true; interrupted refinement reports completed steps via the progress Ref in a finally.
Report: matches the diff; counts (517 pass / 1 broken; 60490 full) are plausible and unchanged by the round-2 edit; the CUDA failure of the previous head is explained and the fix recorded. I could not run Julia in this session; CI on this head is pending (build, cpu, cuda) and is handled by the pipeline.
Blocking findings: none.
Nits (non-blocking): the PR body still quotes the round-0 test counts (453) while the Report says 517; @test_broken r1 <= r0 / 100 could flip to an unexpected pass on a backend with different rounding (measured gain ≤ 42×, so unlikely); permute_rhs!(ws.Y, W.R, …) in refine! still permutes the workspace capacity rather than the solve's nrhs (wasted bandwidth only).
Closes #16
Task: T16 — Iterative refinement, solve sub-phases,
solve_mode, interrupt, loggingWhat was built:
src/solve/refinement.jl: KA gather SpMV over the full-pattern map (symmetric expansion, conjugated mirror for"H"/"HPD"), residual norms reduced on the device,refine!honoringir_n_steps/ir_tol;"solve_refinement"phase,"solve"= sweeps + refinement (the six sub-phases compose bitwise to"solve"); data parameterir_n_steps= steps performed; LinearAlgebra layer defaults toir_n_steps = 2.solve_mode0/1/2 for all symmetric structures via conjugated permutations (op(A)is alwaysMorconj(M)); this also enables complex Hermitian CSC input.user_host_interruptpolled between launch groups, at analysis start and between refinement steps →InterruptedError, solver reset to "analyzed".src/logging.jl:SDS_LOG_LEVEL/SparseDirectSolver.set_log_level!,@info/@debugphase summaries.test/test_refinement.jl,badly_scaled_spdgenerator,bench/refinement.jl(K2 relres table in the Report).Tests:
SDS_TEST_GPU=0 SDS_TEST_ONLY=test_refinement: 453 pass, 0 fail, 1 broken (CPU). FullSDS_TEST_GPU=0suite: 60433 pass, 0 fail, 1 broken. CUDA/AMDGPU: pending CI.Deviations from PLAN.md / TASKS.md: the "≥100× in one step on a badly row-scaled SPD matrix" criterion is unobservable (Cholesky is backward stable under symmetric scaling; measured ≤7×) →
@test_brokenwith explanation; the 100× reduction is asserted on a KKT matrix with static pivot perturbations instead.ir_tolcheck costs one host sync per step, so it only runs whenir_tol > 0; LA layer keepsir_tol = 0. Refinement storage is allocated lazily at the first refining solve. K2 table: case118 reaches 1e-15..1e-18 with 5 steps (better than the cuDSS bar); case1354 stays at 1e-5..1e-7 and needs T21.Follow-up issues opened: none.
🤖 Generated with Claude Code