Add PureUMFPACKFactorization extension (PureUMFPACK.jl)

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"]

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

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,

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)`
 

test/pureumfpack.jl

@@ -0,0 +1,52 @@
+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

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")

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]