From 9a9db911c233174402efde2fa79bbc92bd587b95 Mon Sep 17 00:00:00 2001 From: ChrisRackauckas-Claude Date: Thu, 4 Jun 2026 22:40:44 -0400 Subject: [PATCH 1/3] Add PureUMFPACKFactorization extension wrapping PureUMFPACK.jl Adds an opt-in `PureUMFPACKFactorization <: AbstractSparseFactorization` backed by the pure-Julia UMFPACK port PureUMFPACK.jl, mirroring the existing SuiteSparse `UMFPACKFactorization` and `PureKLUFactorization`. PureUMFPACK is a weak dependency loaded via the new `LinearSolvePureUMFPACKExt` extension; the SuiteSparse `UMFPACKFactorization` and the default polyalgorithm are unchanged. PureUMFPACK has no in-place numeric-refactorization (`lu!`-style) API, so each fresh solve refactorizes via `splu`; `reuse_symbolic`/`check_pattern` are accepted for API parity but cannot skip the symbolic step. Singular factors (zero on U's diagonal under `check=false`) map to `ReturnCode.Infeasible`. PureUMFPACK.jl is not yet registered in General, so the test adds it by URL. Co-Authored-By: Chris Rackauckas --- Project.toml | 6 ++- ext/LinearSolvePureUMFPACKExt.jl | 92 ++++++++++++++++++++++++++++++++ src/LinearSolve.jl | 2 +- src/factorization.jl | 40 ++++++++++++++ test/pureumfpack.jl | 59 ++++++++++++++++++++ test/runtests.jl | 4 ++ 6 files changed, 201 insertions(+), 2 deletions(-) create mode 100644 ext/LinearSolvePureUMFPACKExt.jl create mode 100644 test/pureumfpack.jl diff --git a/Project.toml b/Project.toml index e1a28ee3d..8a5f28d20 100644 --- a/Project.toml +++ b/Project.toml @@ -58,6 +58,7 @@ ParU_jll = "9e0b026c-e8ce-559c-a2c4-6a3d5c955bc9" Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9" PartitionedSolvers = "11b65f7f-80ac-401b-9ef2-3db765482d62" +PureUMFPACK = "b7e1f0a2-3c4d-4e5f-9a0b-1c2d3e4f5a6b" RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4" STRUMPACK_jll = "86fbd0b9-476f-557c-b766-62c724b42d8c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -94,6 +95,7 @@ LinearSolvePETScExt = ["PETSc", "SparseArrays", "SparseMatricesCSR"] LinearSolvePETScMPIExt = ["PETSc", "PartitionedArrays", "SparseArrays", "SparseMatricesCSR"] LinearSolvePartitionedSolversExt = ["PartitionedArrays", "PartitionedSolvers"] LinearSolveParUExt = ["ParU_jll", "SparseArrays"] +LinearSolvePureUMFPACKExt = ["PureUMFPACK", "SparseArrays"] LinearSolvePardisoExt = ["Pardiso", "SparseArrays"] LinearSolveRecursiveFactorizationExt = "RecursiveFactorization" LinearSolveSTRUMPACKExt = ["SparseArrays", "STRUMPACK_jll"] @@ -151,6 +153,7 @@ Pkg = "1.10" PrecompileTools = "1.2" Preferences = "1.4" PureKLU = "1.0.1" +PureUMFPACK = "0.1" Random = "1.10" RecursiveArrayTools = "3.37, 4" RecursiveFactorization = "0.2.26" @@ -189,6 +192,7 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" KrylovPreconditioners = "45d422c2-293f-44ce-8315-2cb988662dec" MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +PureUMFPACK = "b7e1f0a2-3c4d-4e5f-9a0b-1c2d3e4f5a6b" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4" STRUMPACK_jll = "86fbd0b9-476f-557c-b766-62c724b42d8c" @@ -201,4 +205,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["AlgebraicMultigrid", "BandedMatrices", "BlockDiagonals", "CliqueTrees", "ComponentArrays", "FastAlmostBandedMatrices", "FastLapackInterface", "FiniteDiff", "ForwardDiff", "InteractiveUtils", "IterativeSolvers", "KrylovKit", "KrylovPreconditioners", "MultiFloats", "Pkg", "Random", "RecursiveFactorization", "STRUMPACK_jll", "SafeTestsets", "SparseArrays", "Sparspak", "SpecializingFactorizations", "StaticArrays", "Test", "Zygote"] +test = ["AlgebraicMultigrid", "BandedMatrices", "BlockDiagonals", "CliqueTrees", "ComponentArrays", "FastAlmostBandedMatrices", "FastLapackInterface", "FiniteDiff", "ForwardDiff", "InteractiveUtils", "IterativeSolvers", "KrylovKit", "KrylovPreconditioners", "MultiFloats", "Pkg", "PureUMFPACK", "Random", "RecursiveFactorization", "STRUMPACK_jll", "SafeTestsets", "SparseArrays", "Sparspak", "SpecializingFactorizations", "StaticArrays", "Test", "Zygote"] diff --git a/ext/LinearSolvePureUMFPACKExt.jl b/ext/LinearSolvePureUMFPACKExt.jl new file mode 100644 index 000000000..e28bd2c6f --- /dev/null +++ b/ext/LinearSolvePureUMFPACKExt.jl @@ -0,0 +1,92 @@ +module LinearSolvePureUMFPACKExt + +using LinearSolve: LinearSolve, PureUMFPACKFactorization, OperatorAssumptions, + LinearVerbosity +using PureUMFPACK: PureUMFPACK, PureLU, splu +using SparseArrays: SparseArrays, AbstractSparseArray, SparseMatrixCSC, + nonzeros, rowvals, getcolptr +using SciMLOperators: AbstractSciMLOperator, has_concretization +using SciMLLogging: @SciMLMessage +using SciMLBase: SciMLBase, ReturnCode +using LinearAlgebra: diag + +# PureUMFPACK has no preallocatable / in-place-refactorable factorization object +# (no `lu!`-style API): a fresh `splu` recomputes ordering + symbolic + numerics. +# `init_cacheval` returns an empty typed `PureLU` purely to fix the cache field's +# type so `solve!` can store the real factorization; the first solve factorizes. + +function _empty_purelu(::Type{Tv}, ::Type{Ti}) where {Tv, Ti} + Tr = real(Tv) + empty = SparseMatrixCSC{Tv, Ti}(0, 0, Ti[1], Ti[], Tv[]) + return PureLU{Tv, Ti, Tr}(empty, empty, Ti[], Ti[], Tr[], empty) +end + +function LinearSolve.init_cacheval( + alg::PureUMFPACKFactorization, A::AbstractArray, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions + ) + return nothing +end + +function LinearSolve.init_cacheval( + alg::PureUMFPACKFactorization, A::AbstractSparseArray{Tv, Ti}, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions + ) where {Tv, Ti <: Integer} + return _empty_purelu(Tv, Ti) +end + +function LinearSolve.init_cacheval( + alg::PureUMFPACKFactorization, A::AbstractSciMLOperator, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions + ) + if has_concretization(A) + return LinearSolve.init_cacheval( + alg, convert(AbstractMatrix, A), b, u, Pl, Pr, + maxiters, abstol, reltol, verbose, assumptions + ) + else + return nothing + end +end + +function SciMLBase.solve!( + cache::LinearSolve.LinearCache, alg::PureUMFPACKFactorization; kwargs... + ) + A = cache.A + A = convert(AbstractMatrix, A) + if cache.isfresh + Asp = SparseMatrixCSC( + size(A)..., getcolptr(A), rowvals(A), nonzeros(A) + ) + # No symbolic-reuse / in-place-refactor API in PureUMFPACK: a fresh `splu` + # recomputes ordering, symbolic analysis, and numerics together, so each + # fresh solve refactorizes. `reuse_symbolic`/`check_pattern` are accepted + # for API parity with `UMFPACKFactorization` but cannot skip the symbolic + # step here. `check = false` keeps singular matrices from throwing; the + # diagonal check below maps that to `ReturnCode.Infeasible`. + fact = splu(Asp; check = false) + cache.cacheval = fact + cache.isfresh = false + end + + F = LinearSolve.@get_cacheval(cache, :PureUMFPACKFactorization) + # A zero on U's diagonal means a singular (or numerically singular) factor; + # PureUMFPACK with `check = false` produces it instead of throwing. + return if !any(iszero, diag(F.U)) + y = PureUMFPACK.solve(F, cache.b) + copyto!(cache.u, y) + SciMLBase.build_linear_solution( + alg, cache.u, nothing, cache; retcode = ReturnCode.Success + ) + else + @SciMLMessage("Solver failed", cache.verbose, :solver_failure) + SciMLBase.build_linear_solution( + alg, cache.u, nothing, cache; retcode = ReturnCode.Infeasible + ) + end +end + +end # module diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index ae40d1eee..4726e9fd5 100644 --- a/src/LinearSolve.jl +++ b/src/LinearSolve.jl @@ -520,7 +520,7 @@ export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization, GenericLUFactorization, SimpleLUFactorization, RFLUFactorization, ButterflyFactorization, NormalCholeskyFactorization, NormalBunchKaufmanFactorization, UMFPACKFactorization, KLUFactorization, PureKLUFactorization, - SparseColumnPivotedQRFactorization, FastLUFactorization, + PureUMFPACKFactorization, SparseColumnPivotedQRFactorization, FastLUFactorization, FastQRFactorization, SparspakFactorization, DiagonalFactorization, CholeskyFactorization, BunchKaufmanFactorization, CHOLMODFactorization, LDLtFactorization, diff --git a/src/factorization.jl b/src/factorization.jl index b0961658d..70c9e0079 100644 --- a/src/factorization.jl +++ b/src/factorization.jl @@ -1278,6 +1278,46 @@ function init_cacheval( return nothing end +""" +`PureUMFPACKFactorization(; reuse_symbolic = true, check_pattern = true)` + +A pure-Julia port of SuiteSparse's UMFPACK unsymmetric sparse LU solver, provided +by [PureUMFPACK.jl](https://github.com/SciML/PureUMFPACK.jl). It has no SuiteSparse +binary dependency and supports generic element types in addition to +`Float64`/`ComplexF64`. It is the pure-Julia analogue of the SuiteSparse-backed +[`UMFPACKFactorization`](@ref). + +!!! note + + `PureUMFPACKFactorization` is only available once the `PureUMFPACK` package is + loaded (`using PureUMFPACK`). Unlike SuiteSparse UMFPACK, PureUMFPACK has no + in-place numeric-refactorization (`lu!`-style) API: each fresh factorization + recomputes the ordering, symbolic analysis, and numerics together. The + `reuse_symbolic` and `check_pattern` keywords are accepted for API parity with + [`UMFPACKFactorization`](@ref) and control caching of the factorization object + across solves, but no symbolic factorization is shared between numeric refactors. + +## Keyword Arguments + + - `reuse_symbolic`: reuse the cached factorization across solves when the sparsity + pattern is unchanged. Defaults to `true`. + - `check_pattern`: check whether the sparsity pattern changed before reusing the + cached factorization. Defaults to `true`. +""" +Base.@kwdef struct PureUMFPACKFactorization <: AbstractSparseFactorization + reuse_symbolic::Bool = true + check_pattern::Bool = true +end + +function init_cacheval( + alg::PureUMFPACKFactorization, + A, b, u, Pl, Pr, + maxiters::Int, abstol, reltol, + verbose::Union{LinearVerbosity, Bool}, assumptions::OperatorAssumptions + ) + return nothing +end + """ `SparseColumnPivotedQRFactorization(; reuse_symbolic = true, ordering = :default)` diff --git a/test/pureumfpack.jl b/test/pureumfpack.jl new file mode 100644 index 000000000..504cbb609 --- /dev/null +++ b/test/pureumfpack.jl @@ -0,0 +1,59 @@ +using Pkg +# PureUMFPACK.jl is not yet registered in the General registry; add it by URL so +# the extension `LinearSolvePureUMFPACKExt` loads. +if Base.identify_package("PureUMFPACK") === nothing + Pkg.add(url = "https://github.com/SciML/PureUMFPACK.jl.git") +end + +using LinearSolve +import PureUMFPACK +using SparseArrays +using LinearAlgebra +using Test + +@test Base.get_extension(LinearSolve, :LinearSolvePureUMFPACKExt) !== nothing + +n = 40 +A = sprand(n, n, 0.2) + 5I +b = rand(n) + +@testset "PureUMFPACKFactorization: square solve" begin + prob = LinearProblem(A, b) + sol = solve(prob, PureUMFPACKFactorization()) + @test sol.retcode == ReturnCode.Success + @test norm(A * sol.u - b) < 1.0e-9 +end + +@testset "PureUMFPACKFactorization: caching / refactor" begin + prob = LinearProblem(A, b) + cache = init(prob, PureUMFPACKFactorization()) + sol1 = solve!(cache) + @test norm(A * sol1.u - b) < 1.0e-9 + + # New rhs, same factorization (cached) -> reuse + b2 = rand(n) + cache.b = b2 + sol2 = solve!(cache) + @test norm(A * sol2.u - b2) < 1.0e-9 + + # New matrix values, same sparsity pattern -> refactor + A2 = copy(A) + A2.nzval .*= 2.5 + cache.A = A2 + sol3 = solve!(cache) + @test norm(A2 * sol3.u - b2) < 1.0e-9 +end + +@testset "PureUMFPACKFactorization: reuse_symbolic = false" begin + prob = LinearProblem(A, b) + sol = solve(prob, PureUMFPACKFactorization(reuse_symbolic = false)) + @test sol.retcode == ReturnCode.Success + @test norm(A * sol.u - b) < 1.0e-9 +end + +@testset "PureUMFPACKFactorization: singular -> Infeasible" begin + As = sparse([1, 2, 1, 2], [1, 1, 2, 2], [1.0, 2.0, 2.0, 4.0], 2, 2) + prob = LinearProblem(As, [1.0, 2.0]) + sol = solve(prob, PureUMFPACKFactorization()) + @test sol.retcode == ReturnCode.Infeasible +end diff --git a/test/runtests.jl b/test/runtests.jl index 8b91949db..1a4acb6ba 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -143,6 +143,10 @@ else end end + if GROUP == "LinearSolvePureUMFPACK" + @time @safetestset "PureUMFPACK" include("pureumfpack.jl") + end + # ParU_jll requires Julia >= 1.12 (SuiteSparse_jll in older stdlib is incompatible) if GROUP == "LinearSolveParU" && VERSION >= v"1.12.0-" activate_group_env("paru") From c553a208ac73d8057e068d638690a0bf7b868013 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 11 Jun 2026 08:39:43 +0000 Subject: [PATCH 2/3] Update pureumfpack.jl --- test/pureumfpack.jl | 7 ------- 1 file changed, 7 deletions(-) diff --git a/test/pureumfpack.jl b/test/pureumfpack.jl index 504cbb609..ded464aeb 100644 --- a/test/pureumfpack.jl +++ b/test/pureumfpack.jl @@ -1,10 +1,3 @@ -using Pkg -# PureUMFPACK.jl is not yet registered in the General registry; add it by URL so -# the extension `LinearSolvePureUMFPACKExt` loads. -if Base.identify_package("PureUMFPACK") === nothing - Pkg.add(url = "https://github.com/SciML/PureUMFPACK.jl.git") -end - using LinearSolve import PureUMFPACK using SparseArrays From 58a6b3e797e8708d2499b7d8e3eb289b1b75c7e3 Mon Sep 17 00:00:00 2001 From: Christopher Rackauckas Date: Thu, 11 Jun 2026 08:40:30 +0000 Subject: [PATCH 3/3] Update test_groups.toml --- test/test_groups.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/test_groups.toml b/test/test_groups.toml index fd0457df4..66b5de1ab 100644 --- a/test/test_groups.toml +++ b/test/test_groups.toml @@ -88,6 +88,9 @@ versions = ["1", "pre"] runner = "self-hosted" num_threads = 2 +[LinearSolvePureUMFPACK] +versions = ["lts", "1", "pre"] + # Trim ran on GitHub-hosted ubuntu-latest (self-hosted was false for this # group); keep it on the default ubuntu-latest runner. [Trim]