Route reducible sparse-J WOperators to the sparse LHL solver (cost-model-gated default) - #1247
Closed
ChrisRackauckas wants to merge 3 commits into
Closed
Route reducible sparse-J WOperators to the sparse LHL solver (cost-model-gated default)#1247ChrisRackauckas wants to merge 3 commits into
ChrisRackauckas wants to merge 3 commits into
Conversation
…olver LHLFactorization.jl's `LHLFactorization` algorithm now handles a sparse `J` through the block-triangular sparse LHL of that package's SparseArrays + PureKLU extension (both are LinearSolve dependencies, so the extension is always available here). For a sparse `J`, `init_cacheval` builds the factorization with `lhl(J)` — analyze + reduce — and holds it in the same `LHLCache`; the solve, the shift move (`update_gamma!`) and the re-reduction drive it through the same verbs (`lhl_ldiv!`, `lhl_shift!`, `lhl!`, `lhl_refine!`), made polymorphic over the dense `LHLWorkspace` and the sparse factorization by `lhl_isreduced`, `_lhl_size1` and `_lhl_do_reduce!`. `defaultalg` selects `LHLFactorization` for a sparse-`J` `WOperator` only where the sparse solver is expected to win — a *reducible* pattern, decided cheaply by `lhl_prefers_sparse` (its symmetric block triangular form has more than one block, and its largest irreducible block is within the dense-kernel range). One big irreducible block is KLU's regime and is left to the sparse-LU default; the per-block `kernel = :auto` cost model then makes the fine choice at factorization time. A complex `gamma` on a real sparse `J` keeps the reduction real (the Radau split), as it already does for a dense `J`. Needs LHLFactorization ≥ 2.1 (the release adding the sparse extension and the `lhl_isreduced` / `lhl_prefers_sparse` hooks). The adjoint-reuse hook is left as `_NoAdjointFactorizationReuse` for now: LHLFactorization's `lhl_ldivH!` currently fails its own adjoint tests, so wiring it would return wrong answers — tracked separately. Tests in test/Core/lhl.jl: a reducible sparse WOperator solves and defaults to LHL, `update_gamma!` reuses the reduction, one big irreducible block is not defaulted to LHL (but opt-in still solves), and a complex gamma on a real sparse J. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018WbHDQHK9XRA9xUqTmW2N2
…ction Now that LHLFactorization's `lhl_ldivH!` is correct on every ISA (SciML/LHLFactorization.jl#16), the reverse pass reuses the reduction instead of refactorizing `adjoint(A)`. `LHLFactorization` moves from `_NoAdjointFactorizationReuse` to `_CustomAdjointFactorizationReuse`: - `_custom_adjoint_factorization_solve` syncs the cached reduction+shift to `A` (`_lhl_sync!`) and solves `Aᴴ x = b = Z⁻ᴴ(σI+τH)⁻ᴴZᴴ b` with `lhl_ldivH!` — the same three O(n²) phases as the forward solve, no refactorization — then mirrors the forward iterative refinement with `lhl_refineH!`. Handles a vector or a batched (matrix) right-hand side, real or complex shift. - `_custom_can_reuse_adjoint_factorization` is true only for the dense workspace; the block-triangular sparse LHL has no adjoint kernel, so its adjoint is solved against the assembled system `σI + τJ` with a fresh sparse LU (no reuse). Two subtleties this works around: - A `WOperator`'s `adjoint` does not track `update_gamma!` — its forward `mul!` reads γ live, but `A'` keeps the γ it was built with. Refinement therefore forms its residual against the assembled `σI + τJ`, not `A'`. - The generic reverse fallback (`solve(LinearProblem(adjoint(A), b), alg)`) cannot factorize an `adjoint(WOperator)`, which is why the sparse branch assembles the matrix here rather than declining to the fallback. Tests in test/Core/lhl.jl: reuse over several γ, across `update_gamma!`, complex shift on a real J, a fully complex J, a bare matrix, a batched RHS, and the sparse decline-to-assembled path — each checked against `adjoint(σI+τJ) \ b`. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018WbHDQHK9XRA9xUqTmW2N2
The adjoint reuse path relies on a correct `lhl_ldivH!` on every ISA. The 128-bit/NEON adjoint kernels were fixed across 2.1.0 (#16) and 2.1.1 (#17-#19, which generalized the fold count and fixed the 2- and 8-lane sweeps), so pin the minimum to 2.1.1 to guarantee correctness wherever the reverse pass runs. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_018WbHDQHK9XRA9xUqTmW2N2
ChrisRackauckas
force-pushed
the
feat/lhl-sparse
branch
from
August 23, 2026 11:18
8477324 to
ab439b4
Compare
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What and why
LHLFactorization(theWOperatorre-shift algorithm, #1215) currently accepts only a denseJ, and its reverse pass refactorizesadjoint(A)from scratch. This teaches it — and the default solver — a sparseJ, and lets the reverse pass reuse the reduction too.Sparse J
Via the block-triangular sparse LHL added to LHLFactorization.jl in SciML/LHLFactorization.jl#15 (a
SparseArrays+PureKLUextension; both are already LinearSolve dependencies). For a sparseJ,init_cachevalbuilds the factorization withlhl(J)and holds it in the sameLHLCache; the solve, the shift move (update_gamma!) and re-reduction drive it through the same verbs, made polymorphic over the denseLHLWorkspaceand the sparse factorization by three small hooks (lhl_isreduced,_lhl_size1,_lhl_do_reduce!). A complexgammaon a real sparseJkeeps the reduction real.defaultalgselectsLHLFactorizationfor a sparse-JWOperatoronly where it is expected to win — a reducible pattern, decided cheaply bylhl_prefers_sparse(J)(its symmetric block triangular form has more than one block and its largest irreducible block is within the dense-kernel range). One big irreducible block is KLU's regime and is left to the sparse-LU default; the per-blockkernel = :autocost model makes the fine choice at factorization time.Adjoint reuse (now included)
LHLFactorizationmoves from_NoAdjointFactorizationReuseto_CustomAdjointFactorizationReuse. The reverse pass now solvesAᴴ x = b = Z⁻ᴴ(σI+τH)⁻ᴴZᴴ bwithlhl_ldivH!— the same three O(n²) phases as the forward solve, no refactorization — with the forward iterative refinement mirrored bylhl_refineH!. Vector or batched (matrix) RHS, real or complex shift.This was previously omitted because LHLFactorization's
lhl_ldivH!gave wrongFloat64answers on 128-bit-SIMD targets; SciML/LHLFactorization.jl#16 fixes that, so it is safe to wire now. The dense workspace reuses the reduction; the block-triangular sparse factorization has no adjoint kernel, so its adjoint is solved against the assembledσI + τJwith a fresh sparse LU (correct, no reuse).Two subtleties handled: a
WOperator'sadjointdoes not trackupdate_gamma!(forwardmul!reads γ live,A'keeps the build-time γ), so refinement forms its residual against the assembledσI + τJ, notA'; and the generic reverse fallback cannot factorize anadjoint(WOperator), which is why the sparse branch assembles the matrix rather than declining to it.Tests
test/Core/lhl.jl: a reducible sparseWOperatorsolves and defaults to LHL;update_gamma!reuses the reduction; one big irreducible block is not defaulted to LHL (opt-in still solves); a complexgammaon a real sparseJ; and adjoint reuse across several γ, acrossupdate_gamma!, complex shift on a real J, a fully complex J, a bare matrix, a batched RHS, and the sparse decline-to-assembled path — each checked againstadjoint(σI+τJ) \ b. The full existing dense suite is unchanged.Dependency
Project.tomlcompat:LHLFactorization = "2.1.1"— the release adding the sparse extension (#15), thelhl_isreduced/lhl_prefers_sparsehooks, and the complete adjoint kernels (#16 plus the aarch64/NEON follow-ups #17–#19). It is registered; CI passes once the package-server mirrors propagate it.🤖 Generated with Claude Code
https://claude.ai/code/session_018WbHDQHK9XRA9xUqTmW2N2