From b1ea31521ac399e4887fa1381fc0f291e94baea6 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Mon, 29 Jun 2026 07:24:44 -0400 Subject: [PATCH 1/6] Add mixed-type ldiv! for primal KLU factor against a Dual RHS MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- ext/LinearSolveForwardDiffExt.jl | 37 ++++++++++++++++++++++++++++++ test/Core/forwarddiff_overloads.jl | 29 +++++++++++++++++++++++ 2 files changed, 66 insertions(+) diff --git a/ext/LinearSolveForwardDiffExt.jl b/ext/LinearSolveForwardDiffExt.jl index 4acc7e9ad..dce66bb77 100644 --- a/ext/LinearSolveForwardDiffExt.jl +++ b/ext/LinearSolveForwardDiffExt.jl @@ -7,6 +7,7 @@ using LinearAlgebra using SparseArrays using ForwardDiff using ForwardDiff: Dual, Partials +using PureKLU: PureKLU using SciMLBase using RecursiveArrayTools using SciMLLogging @@ -704,4 +705,40 @@ function update_partials_list!(partial_matrix::SparseMatrixCSC, list_cache) return list_cache end +# Backsolve a real (primal) KLU factorization against a Dual right-hand side. +# +# 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{<:Real}` can +# solve the primal value and every partial column without ever promoting `A` to a Dual — +# i.e. the factorization stays in fast `Float64` arithmetic and the duals only ride +# through the back-substitution. This is the single-solve analogue of the split path's +# per-partial back-solves: it lets a duals-in-`b` problem be solved natively on a primal +# factorization (see #1064 discussion). PureKLU's matrix solve handles the `N+1` columns +# in one `klu_solve!` call. +# +# Note: this is a deliberate, narrowly-typed addition over `PureKLU.AbstractKLUFactorization` +# (a hard dependency) and `ForwardDiff.Dual`; it does not promote `A` and is only reachable +# when a primal KLU factor meets a Dual RHS. The factor eltype is constrained to +# `AbstractFloat` (which a `Dual` is *not*, despite `Dual <: Real`) so this never collides +# with PureKLU's own `ldiv!(::KLUFactorization{Tv}, ::VecOrMat{Tv})` on the dual-factor path. +function LinearAlgebra.ldiv!( + K::PureKLU.AbstractKLUFactorization{Tv}, + b::AbstractVector{Dual{T, V, N}} + ) where {Tv <: AbstractFloat, T, V, N} + n = length(b) + rhs = Matrix{Tv}(undef, n, N + 1) + @inbounds for i in 1:n + rhs[i, 1] = ForwardDiff.value(b[i]) + p = ForwardDiff.partials(b[i]) + for j in 1:N + rhs[i, j + 1] = p[j] + end + end + PureKLU.solve!(K, rhs) # in-place multi-RHS solve; reuses the primal factorization + @inbounds for i in 1:n + b[i] = Dual{T, V, N}(rhs[i, 1], Partials{N, V}(ntuple(j -> V(rhs[i, j + 1]), Val(N)))) + end + return b +end + end diff --git a/test/Core/forwarddiff_overloads.jl b/test/Core/forwarddiff_overloads.jl index 978a4faa3..6ac154836 100644 --- a/test/Core/forwarddiff_overloads.jl +++ b/test/Core/forwarddiff_overloads.jl @@ -280,6 +280,35 @@ plain_A = ForwardDiff.value.(A) prob = LinearProblem(sparse(plain_A), b) @test ≈(solve(prob, PureKLUFactorization()), plain_A \ b, rtol = 1.0e-9) +# Mixed-type ldiv!: a primal (Float64) KLU factorization backsolving a Dual RHS +# without promoting A. The factor stays Float64; duals ride through the +# back-substitution (value + each partial column solved in one multi-RHS solve). +@testset "PureKLU primal factor \\ Dual RHS (mixed ldiv!)" begin + Asp = sparse(2.0I, 5, 5) + sparse(plain_A[1, 1] * 0.0I, 5, 5) + for i in 1:4 + Asp[i, i + 1] = 0.3 + Asp[i + 1, i] = 0.2 + end + for nchunk in (1, 2, 3) + bd = [ + ForwardDiff.Dual{Nothing, Float64, nchunk}( + Float64(i), ForwardDiff.Partials(ntuple(k -> sin(i + k), nchunk)) + ) for i in 1:5 + ] + cache = LinearSolve.__init(LinearProblem(Asp, bd), PureKLUFactorization()) + @test eltype(cache.A) == Float64 # A not promoted + u = solve!(cache).u + uref = Matrix{eltype(bd)}(Asp) \ bd + @test isapprox(ForwardDiff.value.(u), ForwardDiff.value.(uref); rtol = 1.0e-10) + @test all( + isapprox( + ForwardDiff.partials(u[i], j), ForwardDiff.partials(uref[i], j); + rtol = 1.0e-8, atol = 1.0e-12 + ) for i in 1:5, j in 1:nchunk + ) + end +end + A, b = h([ForwardDiff.Dual(5.0, 1.0, 0.0), ForwardDiff.Dual(5.0, 0.0, 1.0)]) prob = LinearProblem(sparse(A), sparse(b)) From a063f08db1ad07c654e45eae36d1b07f3a0c181b Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Mon, 29 Jun 2026 07:39:15 -0400 Subject: [PATCH 2/6] Route duals-in-b PureKLU solves to a native LinearCache MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 #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 --- ext/LinearSolveForwardDiffExt.jl | 11 +++++++++++ test/Core/forwarddiff_overloads.jl | 9 +++++++++ 2 files changed, 20 insertions(+) diff --git a/ext/LinearSolveForwardDiffExt.jl b/ext/LinearSolveForwardDiffExt.jl index dce66bb77..d97c80342 100644 --- a/ext/LinearSolveForwardDiffExt.jl +++ b/ext/LinearSolveForwardDiffExt.jl @@ -285,6 +285,17 @@ function SciMLBase.init(prob::DualAbstractLinearProblem, alg::SparspakFactorizat return __init(prob, alg, args...; kwargs...) end +# Duals only in b (A is primal): route PureKLU to a plain LinearCache and solve natively. +# The mixed-type `ldiv!` (primal KLU factor against a Dual RHS) keeps the factorization in +# Float64 and pushes the duals through the back-substitution, so there is no reason to build +# the split DualLinearCache here. 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. It also gets correct factorization reuse across b-only `reinit!`s for free, +# which the split path's `reinit!` does not (it always marks the inner cache fresh). +function SciMLBase.init(prob::DualBLinearProblem, alg::PureKLUFactorization, args...; kwargs...) + return __init(prob, alg, args...; kwargs...) +end + # NOTE: Removed the runtime conditional for DefaultLinearSolver that checked for # GenericLUFactorization. Now always use __dual_init for type stability. function SciMLBase.init(prob::DualAbstractLinearProblem, alg::DefaultLinearSolver, args...; kwargs...) diff --git a/test/Core/forwarddiff_overloads.jl b/test/Core/forwarddiff_overloads.jl index 6ac154836..24dd6f5aa 100644 --- a/test/Core/forwarddiff_overloads.jl +++ b/test/Core/forwarddiff_overloads.jl @@ -307,6 +307,15 @@ prob = LinearProblem(sparse(plain_A), b) ) for i in 1:5, j in 1:nchunk ) end + + # Duals-only-in-b is routed to a plain LinearCache (native solve), not the split + # DualLinearCache, and that routing is type-stable. + Adual, bdual = h([ForwardDiff.Dual(5.0, 1.0, 0.0), ForwardDiff.Dual(5.0, 0.0, 1.0)]) + plain_Asp = sparse(ForwardDiff.value.(Adual)) + bprob = LinearProblem(plain_Asp, bdual) + @test init(bprob, PureKLUFactorization()) isa LinearSolve.LinearCache + @test (@inferred init(bprob, PureKLUFactorization())) isa LinearSolve.LinearCache + @test ≈(solve(bprob, PureKLUFactorization()), Adual \ bdual, rtol = 1.0e-9) end A, b = h([ForwardDiff.Dual(5.0, 1.0, 0.0), ForwardDiff.Dual(5.0, 0.0, 1.0)]) From 90baeedb66e00e8d5e1f660809de593fa8bddedf Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Mon, 29 Jun 2026 12:09:27 -0400 Subject: [PATCH 3/6] Move the primal-factor/Dual-RHS ldiv! to PureKLU; depend on PureKLU 1.1 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- Project.toml | 2 +- ext/LinearSolveForwardDiffExt.jl | 46 +++++--------------------------- 2 files changed, 8 insertions(+), 40 deletions(-) diff --git a/Project.toml b/Project.toml index f3f8715ee..611681111 100644 --- a/Project.toml +++ b/Project.toml @@ -158,7 +158,7 @@ PartitionedSolvers = "0.3" Pkg = "1.10" PrecompileTools = "1.2" Preferences = "1.4" -PureKLU = "1.0.1" +PureKLU = "1.1" PureUMFPACK = "0.1" Random = "1.10" RecursiveArrayTools = "3.37, 4" diff --git a/ext/LinearSolveForwardDiffExt.jl b/ext/LinearSolveForwardDiffExt.jl index d97c80342..54c2244d5 100644 --- a/ext/LinearSolveForwardDiffExt.jl +++ b/ext/LinearSolveForwardDiffExt.jl @@ -7,7 +7,6 @@ using LinearAlgebra using SparseArrays using ForwardDiff using ForwardDiff: Dual, Partials -using PureKLU: PureKLU using SciMLBase using RecursiveArrayTools using SciMLLogging @@ -286,9 +285,10 @@ function SciMLBase.init(prob::DualAbstractLinearProblem, alg::SparspakFactorizat end # Duals only in b (A is primal): route PureKLU to a plain LinearCache and solve natively. -# The mixed-type `ldiv!` (primal KLU factor against a Dual RHS) keeps the factorization in -# Float64 and pushes the duals through the back-substitution, so there is no reason to build -# the split DualLinearCache here. This is type-stable: dispatch is purely on the problem +# PureKLU's mixed-type `ldiv!` (primal KLU factor against a Dual RHS, from PureKLUForwardDiffExt) +# keeps the factorization in Float64 and pushes the duals through the back-substitution, so +# there is no reason to build the split DualLinearCache here. 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. It also gets correct factorization reuse across b-only `reinit!`s for free, # which the split path's `reinit!` does not (it always marks the inner cache fresh). @@ -716,40 +716,8 @@ function update_partials_list!(partial_matrix::SparseMatrixCSC, list_cache) return list_cache end -# Backsolve a real (primal) KLU factorization against a Dual right-hand side. -# -# 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{<:Real}` can -# solve the primal value and every partial column without ever promoting `A` to a Dual — -# i.e. the factorization stays in fast `Float64` arithmetic and the duals only ride -# through the back-substitution. This is the single-solve analogue of the split path's -# per-partial back-solves: it lets a duals-in-`b` problem be solved natively on a primal -# factorization (see #1064 discussion). PureKLU's matrix solve handles the `N+1` columns -# in one `klu_solve!` call. -# -# Note: this is a deliberate, narrowly-typed addition over `PureKLU.AbstractKLUFactorization` -# (a hard dependency) and `ForwardDiff.Dual`; it does not promote `A` and is only reachable -# when a primal KLU factor meets a Dual RHS. The factor eltype is constrained to -# `AbstractFloat` (which a `Dual` is *not*, despite `Dual <: Real`) so this never collides -# with PureKLU's own `ldiv!(::KLUFactorization{Tv}, ::VecOrMat{Tv})` on the dual-factor path. -function LinearAlgebra.ldiv!( - K::PureKLU.AbstractKLUFactorization{Tv}, - b::AbstractVector{Dual{T, V, N}} - ) where {Tv <: AbstractFloat, T, V, N} - n = length(b) - rhs = Matrix{Tv}(undef, n, N + 1) - @inbounds for i in 1:n - rhs[i, 1] = ForwardDiff.value(b[i]) - p = ForwardDiff.partials(b[i]) - for j in 1:N - rhs[i, j + 1] = p[j] - end - end - PureKLU.solve!(K, rhs) # in-place multi-RHS solve; reuses the primal factorization - @inbounds for i in 1:n - b[i] = Dual{T, V, N}(rhs[i, 1], Partials{N, V}(ntuple(j -> V(rhs[i, j + 1]), Val(N)))) - end - return b -end +# NOTE: the primal-KLU-factor / Dual-RHS `ldiv!` that the duals-in-b routing above relies on +# lives in PureKLU itself (`PureKLUForwardDiffExt`, non-allocating via a buffer on the +# factorization), since PureKLU owns the factorization type. See the PureKLU compat bound. end From ec373739c96f19d7866099ac778ee9b069312747 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Mon, 29 Jun 2026 12:21:11 -0400 Subject: [PATCH 4/6] Pull PureKLU 1.1.0 from PR #67 via [sources] so #1069 CI resolves 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 --- Project.toml | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/Project.toml b/Project.toml index 611681111..cc07afcee 100644 --- a/Project.toml +++ b/Project.toml @@ -3,6 +3,13 @@ uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" version = "3.87.0" authors = ["SciML"] +# TEMPORARY: the duals-in-b PureKLU routing needs the primal-factor/Dual-RHS `ldiv!` added +# in PureKLU 1.1.0 (SciML/PureKLU.jl#67). Until that is registered, pull it from the PR so CI +# can resolve `PureKLU = "1.1"` (Julia >= 1.11; on older Julia `[sources]` is ignored). Remove +# this block once PureKLU 1.1.0 is in the General registry. +[sources] +PureKLU = {url = "https://github.com/ChrisRackauckas-Claude/PureKLU.jl", rev = "dual-rhs-ldiv"} + [deps] AMD = "14f7f29c-3bd6-536c-9a0b-7339e30b5a3e" ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" From 992c1b80e5c806bdffc2a471b46f9e137933ec5d Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Mon, 29 Jun 2026 12:49:03 -0400 Subject: [PATCH 5/6] Drop temporary [sources] now that PureKLU 1.1.0 is registered 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 --- Project.toml | 7 ------- 1 file changed, 7 deletions(-) diff --git a/Project.toml b/Project.toml index cc07afcee..611681111 100644 --- a/Project.toml +++ b/Project.toml @@ -3,13 +3,6 @@ uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" version = "3.87.0" authors = ["SciML"] -# TEMPORARY: the duals-in-b PureKLU routing needs the primal-factor/Dual-RHS `ldiv!` added -# in PureKLU 1.1.0 (SciML/PureKLU.jl#67). Until that is registered, pull it from the PR so CI -# can resolve `PureKLU = "1.1"` (Julia >= 1.11; on older Julia `[sources]` is ignored). Remove -# this block once PureKLU 1.1.0 is in the General registry. -[sources] -PureKLU = {url = "https://github.com/ChrisRackauckas-Claude/PureKLU.jl", rev = "dual-rhs-ldiv"} - [deps] AMD = "14f7f29c-3bd6-536c-9a0b-7339e30b5a3e" ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" From ec308debc0ace18b619dbe3a0870e102298eddf8 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Tue, 30 Jun 2026 09:58:27 +0000 Subject: [PATCH 6/6] Apply suggestion from @ChrisRackauckas --- ext/LinearSolveForwardDiffExt.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/ext/LinearSolveForwardDiffExt.jl b/ext/LinearSolveForwardDiffExt.jl index 54c2244d5..3d902d4fc 100644 --- a/ext/LinearSolveForwardDiffExt.jl +++ b/ext/LinearSolveForwardDiffExt.jl @@ -716,8 +716,4 @@ function update_partials_list!(partial_matrix::SparseMatrixCSC, list_cache) return list_cache end -# NOTE: the primal-KLU-factor / Dual-RHS `ldiv!` that the duals-in-b routing above relies on -# lives in PureKLU itself (`PureKLUForwardDiffExt`, non-allocating via a buffer on the -# factorization), since PureKLU owns the factorization type. See the PureKLU compat bound. - end