Default to PureKLU for generic-eltype sparse LU (e.g. BigFloat)#1037
Merged
ChrisRackauckas merged 2 commits intoJun 13, 2026
Merged
Conversation
PureUMFPACKFactorization (added in SciML#1004) is an AbstractSparseFactorization but was missing from resolve.jl's sparse-conversion lists, so the Re-solve testset exercised it with a dense matrix and errored in do_factorization. Add it alongside the other sparse LU factorizations. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
For non-BLAS sparse element types (e.g. BigFloat), the default sparse LU was SparspakFactorization, which requires `using Sparspak` and whose symbolic-factorization reuse produces invalid factorizations when the sparsity pattern changes across solves (the per-solve dropzeros of the nonstructural_zeros path). That manifested as stalled BVP Newton iterations for BigFloat problems (SciML/BoundaryValueDiffEq.jl#509). PureKLU is a pure-Julia, hard dependency that factors any `Number` element type and correctly re-analyzes on pattern changes, so make it the default sparse LU for generic eltypes. Benchmarks (BigFloat, n=10..200) show PureKLU is competitive-to-faster than Sparspak (e.g. n=200: 0.219s vs 0.307s single solve) with comparable accuracy, and it needs no extra `using`. - Generic-eltype `init_cacheval` for PureKLU and SparseColumnPivotedQR so the default polyalgorithm's KLU and SPQR slots are concretely typed for non-BLAS eltypes (the latter is required for the singular-LU QR fallback to `setfield!` a real factorization). - `defaultalg(::AbstractSparseMatrixCSC)` for square systems now selects the KLU (PureKLU) slot. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
4 tasks
ChrisRackauckas-Claude
pushed a commit
to ChrisRackauckas-Claude/LinearSolve.jl
that referenced
this pull request
Jun 14, 2026
The original PR widened the `{T, Int64}` / `{T, Int32}` PureKLU `init_cacheval`
methods from `Union{Float64, ComplexF64}` to `Union{Real, Complex}` to let
ForwardDiff duals dispatch on the direct dual solve path. After rebasing onto
main, that is unnecessary: SciML#1037 ("Default to PureKLU for generic-eltype sparse
LU") added a catch-all `where {T <: Number, Ti <: Integer}` method that already
builds the correct empty cache for any number eltype (duals included). The
widened specializations only duplicated it, so revert that file to main.
Verified: test/core/forwarddiff_overloads.jl (incl. the PureKLU sparse-dual
cases) passes with this file unchanged from main.
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
This was referenced Jun 14, 2026
ChrisRackauckas
added a commit
that referenced
this pull request
Jun 14, 2026
…arspak) since #1037 (#1043) * Update DefaultsLoading test: sparse BigFloat defaults to PureKLU, not Sparspak #1037 changed the default sparse LU for generic (non-BLAS) element types from SparspakFactorization to PureKLU (a hard dependency), removing the "SparspakFactorization required" error branch from `defaultalg`. The DefaultsLoading test still asserted that error was thrown when solving a sparse BigFloat problem before `using Sparspak`, so it failed on master with "No exception thrown" (defaults_loading.jl:30). Update the test to the new intended behavior: a generic-eltype sparse problem now routes to `KLUFactorization` and solves without Sparspak loaded. Closes the `@test_throws` into a positive assertion (defaultalg === KLUFactorization, solve returns Vector{BigFloat}). Fixes #1042. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> * Apply suggestion from @ChrisRackauckas --------- Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
ChrisRackauckas
added a commit
that referenced
this pull request
Jun 14, 2026
* Opt PureKLU and RFLU out of split dual AD path * Update forwarddiff_overloads.jl * Apply suggestion from @ChrisRackauckas * Fix direct dual solve path for PureKLU and non-dual b (#1024) Two fixes for the direct dual opt-out path: 1. PureKLUFactorization's sparse init_cacheval methods were restricted to Float64/ComplexF64 eltypes, so Dual-eltype sparse matrices fell through to the nothing fallback and solve! crashed on cacheval.nzval. PureKLU's kernels are pure Julia and generic over Union{Real, Complex}, so widen the specializations accordingly. 2. The direct dual path built its inner problem from dual_A/dual_b as stored, but when only A carries duals b is a plain array, and __init takes the solution eltype from b, so the dual solution could not be stored (MethodError: no method matching Float64(::Dual)). This also affected the pre-existing GenericLUFactorization opt-out. Promote b to the cache's dual type in _solve_direct_dual!, and only take the direct path when A itself carries duals: with duals just in b, the split path (one primal factorization + partials back-solves) is strictly cheaper and works for all algorithms. Adds test coverage for the duals-only-in-A and duals-only-in-b cases for GenericLUFactorization, RFLUFactorization, and PureKLUFactorization. Co-authored-by: Chris Rackauckas (Claude) <accounts@chrisrackauckas.com> Co-authored-by: Claude Fable 5 <noreply@anthropic.com> * Drop redundant PureKLU init_cacheval widening (subsumed by #1037) The original PR widened the `{T, Int64}` / `{T, Int32}` PureKLU `init_cacheval` methods from `Union{Float64, ComplexF64}` to `Union{Real, Complex}` to let ForwardDiff duals dispatch on the direct dual solve path. After rebasing onto main, that is unnecessary: #1037 ("Default to PureKLU for generic-eltype sparse LU") added a catch-all `where {T <: Number, Ti <: Integer}` method that already builds the correct empty cache for any number eltype (duals included). The widened specializations only duplicated it, so revert that file to main. Verified: test/core/forwarddiff_overloads.jl (incl. the PureKLU sparse-dual cases) passes with this file unchanged from main. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> --------- Co-authored-by: Christopher Rackauckas <accounts@chrisrackauckas.com> Co-authored-by: Claude Fable 5 <noreply@anthropic.com>
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
2 participants
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.
Summary
For non-BLAS sparse element types (e.g.
BigFloat), the default sparse LU wasSparspakFactorization. That has two problems:using Sparspak(otherwise the default errors).nonstructural_zerosper-solvedropzerospath does. This manifested as stalled BVP Newton iterations for BigFloat problems (Misc test group red on master: BigFloat BVP solves stall — regression from LinearSolve v3.85 nonstructural-zeros reduction BoundaryValueDiffEq.jl#509).This PR makes PureKLU (a pure-Julia, hard dependency that factors any
Numberelement type and correctly re-analyzes on pattern changes) the default sparse LU for generic eltypes.Changes
init_cachevalforPureKLUFactorizationandSparseColumnPivotedQRFactorization, so the default polyalgorithm's KLU and SPQR slots are concretely typed for non-BLAS eltypes (the SPQR one is required so the singular-LU → column-pivoted-QR fallback cansetfield!a real factorization into its slot).defaultalg(::AbstractSparseMatrixCSC)for square systems now selects the KLU (PureKLU) slot instead of Sparspak. Nousing Sparspakrequired.test/core/default_algs.jl): BigFloat sparse default chooses KLU, solves accurately, re-solves correctly across numeric updates, and falls back to QR on a singular system.PureUMFPACKFactorizationtest fix (pre-existing master failure; shared with Fix Sparspak symbolic reuse invalidated by per-solve dropzeros #1036).Benchmark (BigFloat, PureKLU vs Sparspak)
Single solve (factorize+solve) and 5× same-pattern resolve, seconds:
PureKLU is competitive-to-faster (notably faster at larger sizes) with comparable accuracy, and needs no extra
using.Test plan
test/core/default_algs.jlpasses (incl. new BigFloat tests)test/core/resolve.jlpassessolve, caching re-solve, and singular→QR fallback verifiedMIRKN4BigFloatSecondOrderBVProblemfrom BoundaryValueDiffEq.jl#509 now returnsSuccess(was Stalled) with this changeRelationship to #1036
#1036 hardens
SparspakFactorizationitself against the same pattern-change bug (defense-in-depth for users who explicitly request Sparspak). This PR addresses the root cause for the default path by not selecting Sparspak for generic eltypes at all. They are complementary; this PR is the more fundamental fix for the BVP issue.