fix sparspak solves in multithreaded loop (#760)

Author: j-fu

ext/LinearSolveSparspakExt.jl

@@ -41,7 +41,7 @@ function SciMLBase.solve!(
         cache::LinearSolve.LinearCache, alg::SparspakFactorization; kwargs...)
     A = cache.A
     if cache.isfresh
-        if cache.cacheval !== nothing && alg.reuse_symbolic
+        if !(cache.cacheval === PREALLOCATED_SPARSEPAK) && alg.reuse_symbolic
             fact = sparspaklu!(LinearSolve.@get_cacheval(cache, :SparspakFactorization),
                 SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),
                     nonzeros(A)))

test/basictests.jl

@@ -679,3 +679,34 @@ end
     u = solve!(cache)
     @test norm(u - u0, Inf) < 1.0e-8
 end
+
+@testset "ParallelSolves" begin
+    n=1000
+    @info "ParallelSolves: Threads.nthreads()=$(Threads.nthreads())"
+    A_sparse = 10I - sprand(n, n, 0.01)
+    B = [rand(n), rand(n)]
+    U = [A_sparse \ B[i] for i in 1:2]
+    sol = similar(U)
+
+    Threads.@threads for i in 1:2
+        sol[i] = solve(LinearProblem(A_sparse, B[i]), UMFPACKFactorization())
+    end
+
+    for i in 1:2
+        @test sol[i] ≈ U[i]
+    end
+    
+    Threads.@threads for i in 1:2
+        sol[i] = solve(LinearProblem(A_sparse, B[i]), KLUFactorization())
+    end
+    for i in 1:2
+        @test sol[i] ≈ U[i]
+    end
+
+    Threads.@threads for i in 1:2
+        sol[i] = solve(LinearProblem(A_sparse, B[i]), SparspakFactorization())
+    end
+    for i in 1:2
+        @test sol[i] ≈ U[i]
+    end
+end