Expose PureKLU pivot tolerance option

Project.toml

@@ -170,7 +170,7 @@ PartitionedSolvers = "0.3"
 Pkg = "1.10"
 PrecompileTools = "1.2"
 Preferences = "1.4"
-PureKLU = "1.1"
+PureKLU = "1.4"
 PureUMFPACK = "0.1"
 Random = "1.10"
 RecursiveArrayTools = "3.37, 4"

ext/LinearSolveSparseArraysExt.jl

@@ -653,10 +653,11 @@ function SciMLBase.solve!(
                         size(A)..., getcolptr(A), rowvals(A),
                         nonzeros(A)
                     ),
-                    check = false, use_fma = alg.use_fma,
+                    check = false, use_fma = alg.use_fma, tol = alg.tol,
                     fully_preallocated = alg.fully_preallocated
                 )
             else
+                # Refactorization reuses pivot ordering; `tol` does not apply
                 fact = PureKLU.klu!(cacheval, nonzeros(A), check = false)
             end
         else
@@ -668,7 +669,7 @@ function SciMLBase.solve!(
                     size(A)..., getcolptr(A), rowvals(A),
                     nonzeros(A)
                 ),
-                check = false, use_fma = alg.use_fma,
+                check = false, use_fma = alg.use_fma, tol = alg.tol,
                 fully_preallocated = alg.fully_preallocated
             )
         end

src/factorization.jl

@@ -1411,7 +1411,7 @@ function init_cacheval(
 end
 
 """
-`PureKLUFactorization(; reuse_symbolic = true, check_pattern = true, use_fma = true, fully_preallocated = nothing)`
+`PureKLUFactorization(; reuse_symbolic = true, check_pattern = true, use_fma = true, fully_preallocated = nothing, tol = 0.001)`
 
 A pure-Julia port of SuiteSparse's KLU sparse LU solver, provided by
 [PureKLU.jl](https://github.com/SciML/PureKLU.jl). It has no SuiteSparse binary
@@ -1438,12 +1438,17 @@ in the default polyalgorithm.
     SuiteSparse `KLUFactorization`. Defaults to `true`.
   - `fully_preallocated`: PureKLU's `fully_preallocated` option. `nothing` (default) lets
     PureKLU choose automatically based on the maximum block size.
+  - `tol`: Pivot on a column's diagonal instead of largest entry if it is at least `tol` times
+    larger in magnitude. Set `tol = 1.0` for partial pivoting, and `tol = 0.0` to always use the
+    diagonal. Only applies to the initial factorization; refactorizations reuse the existing
+    pivot ordering. Defaults to `0.001`.
 """
 Base.@kwdef struct PureKLUFactorization <: AbstractSparseFactorization
     reuse_symbolic::Bool = true
     check_pattern::Bool = true
     use_fma::Bool = true
     fully_preallocated::Union{Bool, Nothing} = nothing
+    tol::Float64 = 0.001
 end
 
 function init_cacheval(

test/Core/basictests.jl

@@ -210,6 +210,22 @@ end
         X = sprand(n, n, 0.8)
         cache.A = X
         @test X * solve!(cache) ≈ b1
+
+        # Partial pivoting tolerance option
+        @test PureKLUFactorization().tol == 0.001 # default value
+        function rowperm(tol)
+            A = sparse([1.0e-3 1.0; 1.0 1.0])
+            b = [1.0, 2.0]
+            alg = PureKLUFactorization(; tol)
+            cache = SciMLBase.init(LinearProblem(A, b), alg)
+            solve!(cache)
+            return cache.cacheval.p # row permutation vector
+        end
+        @test rowperm(0.0) == [1, 2] # below threshold: pivot on diagonal entry
+        @test rowperm(1.0e-4) == [1, 2] # below threshold: pivot on diagonal entry
+        @test rowperm(1.0e-3) == [1, 2] # at threshold: pivot on diagonal entry
+        @test rowperm(1.0e-2) == [2, 1] # above threshold: pivot on largest entry
+        @test rowperm(1.0) == [2, 1] # above threshold: pivot on largest entry
     end
 
     @testset "RFLUFactorization multi-RHS" begin