Update sparspak interface

* Base sparspak interface on sparspaklu
* Introduce tests for MultiFloats and ForwardDiff

Author: j-fu

Project.toml

@@ -44,10 +44,12 @@ UnPack = "1"
 julia = "1.6"
 
 [extras]
+MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Pkg", "Random", "SafeTestsets"]
+test = ["Test", "Pkg", "Random", "SafeTestsets", "MultiFloats", "ForwardDiff"]

src/factorization.jl

@@ -503,30 +503,28 @@ end
 
 ## SparspakFactorization is here since it's MIT licensed, not GPL
 
-struct SparspakFactorization <: AbstractFactorization end
+Base.@kwdef struct SparspakFactorization <: AbstractFactorization
+    reuse_symbolic::Bool = true
+end
 
 function init_cacheval(::SparspakFactorization, A, b, u, Pl, Pr, maxiters::Int, abstol,
                        reltol,
                        verbose::Bool, assumptions::OperatorAssumptions)
     A = convert(AbstractMatrix, A)
-    p = Sparspak.Problem.Problem(size(A)...)
-    Sparspak.Problem.insparse!(p, A)
-    Sparspak.Problem.infullrhs!(p, b)
-    s = Sparspak.SparseSolver.SparseSolver(p)
-    return s
+    return sparspaklu(A)
 end
 
 function SciMLBase.solve(cache::LinearCache, alg::SparspakFactorization; kwargs...)
     A = cache.A
     A = convert(AbstractMatrix, A)
     if cache.isfresh
-        p = Sparspak.Problem.Problem(size(A)...)
-        Sparspak.Problem.insparse!(p, A)
-        Sparspak.Problem.infullrhs!(p, cache.b)
-        s = Sparspak.SparseSolver.SparseSolver(p)
-        cache = set_cacheval(cache, s)
+        if cache.cacheval !== nothing && alg.reuse_symbolic
+            fact = sparspaklu!(cache.cacheval, A)
+        else
+            fact = sparspaklu(A)
+        end
+        cache = set_cacheval(cache, fact)
     end
-    Sparspak.SparseSolver.solve!(cache.cacheval)
-    copyto!(cache.u, cache.cacheval.p.x)
-    SciMLBase.build_linear_solution(alg, cache.u, nothing, cache)
+    y = ldiv!(cache.u, cache.cacheval, cache.b)
+    SciMLBase.build_linear_solution(alg, y, nothing, cache)
 end

test/basictests.jl

@@ -1,7 +1,9 @@
-using LinearSolve, LinearAlgebra, SparseArrays
+using LinearSolve, LinearAlgebra, SparseArrays, MultiFloats, ForwardDiff
 using Test
 import Random
 
+const Dual64=ForwardDiff.Dual{Nothing,Float64,1}
+
 n = 8
 A = Matrix(I, n, n)
 b = ones(n)
@@ -124,7 +126,7 @@ end
         @test X * solve(cache) ≈ b1
     end
 
-    @testset "Sparspak Factorization" begin
+    @testset "Sparspak Factorization (Float64)" begin
         A1 = sparse(A / 1)
         b1 = rand(n)
         x1 = zero(b)
@@ -137,6 +139,46 @@ end
         test_interface(SparspakFactorization(), prob1, prob2)
     end
 
+    @testset "Sparspak Factorization (Float64x1)" begin
+        A1 = sparse(A / 1) .|> Float64x1
+        b1 = rand(n)  .|> Float64x1
+        x1 = zero(b)   .|> Float64x1
+        A2 = sparse(A / 2)  .|> Float64x1
+        b2 = rand(n)  .|> Float64x1
+        x2 = zero(b)  .|> Float64x1
+
+        prob1 = LinearProblem(A1, b1; u0 = x1)
+        prob2 = LinearProblem(A2, b2; u0 = x2)
+        test_interface(SparspakFactorization(), prob1, prob2)
+    end
+
+    @testset "Sparspak Factorization (Float64x2)" begin
+        A1 = sparse(A / 1) .|> Float64x2
+        b1 = rand(n)  .|> Float64x2
+        x1 = zero(b)   .|> Float64x2
+        A2 = sparse(A / 2)  .|> Float64x2
+        b2 = rand(n)  .|> Float64x2
+        x2 = zero(b)  .|> Float64x2
+
+        prob1 = LinearProblem(A1, b1; u0 = x1)
+        prob2 = LinearProblem(A2, b2; u0 = x2)
+        test_interface(SparspakFactorization(), prob1, prob2)
+    end
+
+    @testset "Sparspak Factorization (Dual64)" begin
+        A1 = sparse(A / 1) .|> Dual64
+        b1 = rand(n)  .|> Dual64
+        x1 = zero(b)   .|> Dual64
+        A2 = sparse(A / 2)  .|> Dual64
+        b2 = rand(n)  .|> Dual64
+        x2 = zero(b)  .|> Dual64
+
+        prob1 = LinearProblem(A1, b1; u0 = x1)
+        prob2 = LinearProblem(A2, b2; u0 = x2)
+        test_interface(SparspakFactorization(), prob1, prob2)
+    end
+
+    
     @testset "FastLAPACK Factorizations" begin
         A1 = A / 1
         b1 = rand(n)