Route duals-in-b PureKLU solves to a native LinearCache (uses PureKLU 1.1 mixed ldiv!) by ChrisRackauckas-Claude · Pull Request #1069 · SciML/LinearSolve.jl

ChrisRackauckas-Claude · 2026-06-29T11:39:45Z

Adds native handling of duals-only-in-b PureKLU solves: a DualBLinearProblem (A primal, b dual) is routed to a plain LinearCache and solved with PureKLU's primal-factor / Dual-RHS ldiv!, instead of the split DualLinearCache. Came out of the #1064 / #1065 discussion (cc @hersle).

Depends on PureKLU 1.1.0 (SciML/PureKLU.jl#67), now merged and registered — the PureKLU = "1.1" compat resolves directly.

What's here

init(::DualBLinearProblem, ::PureKLUFactorization) = __init(prob, alg) — routes duals-in-b to a native LinearCache. Type-stable: dispatch is purely on the problem subtype (A primal / b dual) + alg type, so init always returns a LinearCache here (verified with @inferred). This sidesteps the Union{LinearCache, DualLinearCache} instability that caused the value-based opt-out to be reverted (commit 25f159c), because the subtype already encodes "A primal, b dual." Scope is narrow: only PureKLU; A-dual and mixed-A/b problems still take the split path.

Plus a PureKLU [compat] bump to 1.1.

Scope note (this does NOT fix #1064 by itself)

#1064 is an A-dual problem (the Rosenbrock W carries duals). This PR only covers the b-dual case, so it is additive, not a replacement for the A-dual fix (#1066). The two cover different cells:

dual location	handled by
b only (A primal)	this PR + PureKLU 1.1
A / A-and-b (the #1064 case)	#1066 (route A-dual PureKLU → split)

Why it's worth it

Simpler + a latent fix: the native cache reuses the factorization across b-only reinit!s for free. The split path's reinit! unconditionally sets linear_cache.isfresh = true, so it refactorizes even when only b changed. For a fixed operator with a varying dual RHS (cache reused), the native path is 5–15× faster on higher-fill matrices for that reason. (When A changes every step, both refactorize and it's a tie — a simplification + fixed-A win, not a claim of being fundamentally faster than split.)
Non-allocating: the solve allocates nothing after warmup — the n×(N+1) real RHS buffer lives on the PureKLU factorization (PureKLU#67, @allocated == 0).

Verification (Julia 1.12, Linux)

test/Core/forwarddiff_overloads.jl passes end-to-end against registered PureKLU 1.1.0; integration testset 12/12 (correctness through the public solve, routing returns LinearCache, @inferred type-stability). Runic --check clean.

Please ignore until reviewed by @ChrisRackauckas.

🤖 Generated with Claude Code

A factorization is a linear operator, so `K \ Dual(v, p₁…p_N) = Dual(K\v, K\p₁ … K\p_N)`. PureKLU stores its L/U in the factor element type, so a `KLUFactorization{<:AbstractFloat}` can backsolve a Dual right-hand side by solving the primal value and each partial column in one multi-RHS `klu_solve!`, then repacking — without ever promoting `A` to a Dual. This keeps the factorization in fast Float64 arithmetic and lets the duals ride only through the back-substitution, which is the single-solve analogue of the split path's per-partial back-solves. It makes `solve(LinearProblem(A_float64, b_dual), PureKLU)` work natively (previously a MethodError on `ldiv!(::KLUFactorization{Float64}, ::Vector{<:Dual})`). The factor eltype is constrained to `AbstractFloat` so it never collides with PureKLU's own `ldiv!(::KLUFactorization{Tv}, ::VecOrMat{Tv})` on the dual-factor path — a `Dual` is `<: Real` but not `<: AbstractFloat`, so the dual-factor case stays with PureKLU's method. Adds a correctness test (value + every partial column matches the dense Dual solve, factor stays Float64) for chunk sizes 1/2/3. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

With the mixed-type `ldiv!` (primal KLU factor against a Dual RHS) in place, a `DualBLinearProblem` (duals only in b, A primal) no longer needs the split DualLinearCache machinery: it can be solved on a plain `LinearCache` whose factorization stays in Float64, with the duals carried through the back-solve. Add an `init(::DualBLinearProblem, ::PureKLUFactorization)` opt-out that returns `__init(prob, alg)` directly. This is type-stable — dispatch is purely on the problem subtype (b-dual / A-plain) and the alg type, so `init` always returns a `LinearCache` for this method (verified with `@inferred`), unlike the reverted value-based opt-out in SciML#1041-era code. It is also simpler than the split path and gets correct factorization reuse across b-only `reinit!`s for free, which the split path's `reinit!` does not (it unconditionally marks the inner cache fresh, forcing a refactorization even when only b changed). Scope: only PureKLU is routed, since it is the only sparse factorization with the mixed `ldiv!`. A-dual and mixed A/b problems are untouched and still take the split path. Adds routing + type-stability + correctness tests. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

The mixed-type `ldiv!` (primal KLU factor backsolving a Dual RHS) belongs in PureKLU, which owns `KLUFactorization` — it is the analogue of PureKLU's existing real-factor / Complex-RHS `ldiv!`, avoids type piracy here, and there it is non-allocating (the real `n × (N+1)` RHS reuses a buffer on the factorization). See SciML/PureKLU.jl#67. This drops the stop-gap `ldiv!` (and the `PureKLU` import) from this extension and bumps the PureKLU [compat] to 1.1, which provides it. The duals-in-b routing (`init(::DualBLinearProblem, ::PureKLUFactorization)` -> native `LinearCache`) is unchanged and now relies on PureKLU's `ldiv!`. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

…solves The duals-in-b routing needs the primal-factor/Dual-RHS `ldiv!` from PureKLU 1.1.0 (SciML/PureKLU.jl#67), which is not yet registered. Add a temporary `[sources]` entry pointing at the PR branch so `PureKLU = "1.1"` resolves in CI (Julia >= 1.11; ignored on older Julia). To be removed once PureKLU 1.1.0 is in the General registry. Verified locally: a fresh resolve with this `[sources]` pulls PureKLU v1.1.0. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

PureKLU 1.1.0 (with the primal-factor/Dual-RHS `ldiv!`, SciML/PureKLU.jl#67) is in the General registry, so the `PureKLU = "1.1"` compat resolves directly. Remove the temporary `[sources]` override that pointed at the PR branch. Verified: forwarddiff_overloads.jl passes end-to-end against the registered PureKLU 1.1.0. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>

DualBLinearProblem's A slot is `Union{Number, AbstractArray}`, which a Dual-eltype A also matches, so `DualLinearProblem <: DualBLinearProblem`. The b-only PureKLU opt-out introduced in #1069 is more specific in the alg argument than the general dual-path `init`, so both-Dual PureKLU problems were routed to a plain `LinearCache` and `solve!` then failed with a FieldError on `dual_linear_cache` (test/Core/forwarddiff_overloads.jl "PureKLU direct dual path reuses inner LinearCache"). Add a more-specific `DualLinearProblem` method restoring the direct dual path; the mixed-type ldiv! only covers primal A. Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com> Co-authored-by: Claude Fable 5 <noreply@anthropic.com>

The repo-wide Runic Format Check fails on main because this file landed unformatted in aabf007 (SciML#1069); format it here alongside the dispatch fix for the same commit's fallout. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

…0) (#1072) * Support batched (matrix) right-hand sides in factorization algorithms `solve(LinearProblem(A, B))` with `B::AbstractMatrix` now computes the equivalent of `A \ B`: `u0` is initialized as a `size(A, 2) × size(B, 2)` matrix and the factorization-based algorithms solve all columns against a single factorization of `A`. This is a breaking change (v4.0.0): previously a matrix `b` initialized a flattened vector `u` and errored downstream. - `__init_u0_from_Ab` gains matrix-`b` methods (including SMatrix disambiguation). - Iterative/Krylov algorithms (`KrylovJL`, `SimpleGMRES`, and the other `AbstractKrylovSubspaceMethod`s, plus the default algorithm when it selects a Krylov method for an operator) throw an informative `ArgumentError` at `init` time for matrix `b`; `SimpleLUFactorization` (vector-only workspace) does the same. - The default polyalgorithm's unused Krylov cacheval slots initialize to `nothing` for matrix `b` so structured-matrix defaults (Diagonal, SymTridiagonal, ...) work. - `_check_residual_safety` sizes its residual buffer shape-generically. - Default-algorithm size heuristics use `size(b, 1)` instead of `length(b)` so batched right-hand sides don't inflate the apparent problem size (same for the sparse `use_klulike_sparse_structure` heuristic). - MKL/OpenBLAS/AppleAccelerate LU copy solution columns correctly for matrix `b` in the rectangular (`m > n`) path. - SparseArrays extension: CHOLMOD pre-1.12 `_ldiv!` fallback widened to matrices, SPQR gains a matrix `_ldiv!` via `\` (no in-place matrix ldiv! on any Julia version), and SparseColumnPivotedQR solves matrix `b` column-by-column against the one factorization. Part of #552 Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh * Add batched right-hand side tests Covers dense LU/GenericLU/QR (NoPivot + ColumnNorm)/SVD/Cholesky/ BunchKaufman/NormalCholesky vs A \ B for Float64 and ComplexF64, the default algorithm (dense, sparse, structured, non-square with k != n), sparse UMFPACK/KLU/CHOLMOD and the non-square sparse QR default, cache reuse via `cache.b = B2` and `cache.A = A2`, failure retcodes on singular A (including the default's QR fallback), residual-safety with matrix B, informative errors for Krylov methods and SimpleLUFactorization, static arrays, BigFloat, and single-column-matrix vs vector consistency. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh * Document batched right-hand sides Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01EYp371jx6LurezUDhKcYRh * Widen in-repo LinearSolve compat for v4 and source docs env locally All sublibrary, test, qa, GPU, Trim, and docs environments pinned LinearSolve = "3[.x.y]", so every CI job died at resolution against the 4.0.0 bump ("empty intersection between LinearSolve@4.0.0 and project compatibility"). Widen each to include 4. The docs environment also gains [sources] entries for LinearSolve and LinearSolveAutotune (the registered Autotune restricts LinearSolve to 3, so the docs build must use the local copies, matching the pattern the test environments already use). Verified locally: docs, lib/LinearSolveAutotune, and test/qa environments all instantiate against 4.0.0. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> * Runic-format forwarddiff_overloads.jl The repo-wide Runic Format Check fails on main because this file landed unformatted in aabf007 (#1069); format it here alongside the dispatch fix for the same commit's fallout. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> * Allow AlgebraicMultigrid 0.5 so the test target resolves against v4 Every registered AlgebraicMultigrid >= 0.6 caps LinearSolve at 3, so the root test target (which all CI test groups resolve) was unsatisfiable against the 4.0.0 bump. AMG 0.5.x predates its LinearSolve dependency and exposes the identical API surface the extension uses (SmoothedAggregationAMG, RugeStubenAMG, aspreconditioner, ruge_stuben); the AlgebraicMultigridJL testset passes on 0.5.1 (verified locally, including the tight-tolerance case). The resolver will pick AMG 2.x again automatically once AlgebraicMultigrid widens its own compat; the compat comment marks the entry for removal then. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> * Retrigger CI (self-hosted runner disconnect on Downgrade Core) Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> * Drop LinearSolve v3 from in-repo lower bounds All in-repo environments and both registered sublibraries now require LinearSolve 4: the environments always develop the local copy, and v4-era releases of LinearSolveAutotune/LinearSolvePyAMG should require the batched-RHS LinearSolve rather than straddle the major boundary. Minor-bump both sublibraries (Autotune 1.12.0, PyAMG 1.2.0) for the dependency floor raise. All environment classes re-verified to resolve. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> --------- Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com> Co-authored-by: Claude Fable 5 <noreply@anthropic.com>

ChrisRackauckas added 3 commits June 29, 2026 07:24

ChrisRackauckas-Claude changed the title ~~Mixed-type ldiv! (primal KLU factor / Dual RHS) + native duals-in-b routing~~ Route duals-in-b PureKLU solves to a native LinearCache (uses PureKLU 1.1 mixed ldiv!) Jun 29, 2026

ChrisRackauckas added 2 commits June 29, 2026 12:21

ChrisRackauckas reviewed Jun 30, 2026

View reviewed changes

Comment thread ext/LinearSolveForwardDiffExt.jl Outdated

Apply suggestion from @ChrisRackauckas

ec308de

ChrisRackauckas marked this pull request as ready for review June 30, 2026 09:58

ChrisRackauckas merged commit aabf007 into SciML:main Jun 30, 2026
6 of 13 checks passed

This was referenced Jul 3, 2026

Batched (matrix) right-hand sides: LinearProblem(A, B) ≡ A \ B (v4.0.0) #1072

Merged

Fix PureKLU Dual-A dispatch stolen by DualBLinearProblem opt-out #1073

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Uh oh!

Route duals-in-b PureKLU solves to a native LinearCache (uses PureKLU 1.1 mixed ldiv!)#1069

Route duals-in-b PureKLU solves to a native LinearCache (uses PureKLU 1.1 mixed ldiv!)#1069
ChrisRackauckas merged 6 commits into
SciML:mainfrom
ChrisRackauckas-Claude:mixed-ldiv-dual-rhs

ChrisRackauckas-Claude commented Jun 29, 2026 •

edited

Loading

Uh oh!

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Uh oh!

Uh oh!

Conversation

ChrisRackauckas-Claude commented Jun 29, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

What's here

Scope note (this does NOT fix #1064 by itself)

Why it's worth it

Verification (Julia 1.12, Linux)

Uh oh!

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

ChrisRackauckas-Claude commented Jun 29, 2026 •

edited

Loading