Merge pull request #249 from SciML/sparse

Clean up sparse handling

Author: ChrisRackauckas

src/common.jl

@@ -98,6 +98,15 @@ function SciMLBase.init(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorith
                         kwargs...)
     @unpack A, b, u0, p = prob
 
+    A = alias_A ? A : deepcopy(A)
+    b = if b isa SparseArrays.AbstractSparseArray && !(A isa Diagonal)
+        Array(b) # the solution to a linear solve will always be dense!
+    elseif alias_b
+        b
+    else
+        deepcopy(b)
+    end
+
     u0 = if u0 !== nothing
         u0
     else
@@ -110,9 +119,6 @@ function SciMLBase.init(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorith
     isfresh = true
     Tc = typeof(cacheval)
 
-    A = alias_A ? A : deepcopy(A)
-    b = alias_b ? b : deepcopy(b)
-
     cache = LinearCache{
                         typeof(A),
                         typeof(b),

src/default.jl

@@ -53,7 +53,7 @@ end
     end
 else
     function defaultalg(A::SparseMatrixCSC, b, ::OperatorAssumptions{true})
-        KrylovJL_GMRES()
+        SparspakFactorization()
     end
 end
 

src/factorization.jl

@@ -500,3 +500,33 @@ function SciMLBase.solve(cache::LinearCache, alg::FastQRFactorization{P};
     y = ldiv!(cache.u, cache.cacheval.factors, cache.b)
     SciMLBase.build_linear_solution(alg, y, nothing, cache)
 end
+
+## SparspakFactorization is here since it's MIT licensed, not GPL
+
+struct SparspakFactorization <: AbstractFactorization 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
+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)
+    end
+    Sparspak.SparseSolver.solve!(cache.cacheval)
+    copyto!(cache.u, cache.cacheval.p.x)
+    SciMLBase.build_linear_solution(alg, cache.u, nothing, cache)
+end

src/factorization_sparse.jl

@@ -6,31 +6,3 @@ function _ldiv!(x::Vector,
                          SuiteSparse.CHOLMOD.Factor}, b::Vector)
     x .= A \ b
 end
-
-struct SparspakFactorization <: AbstractFactorization 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
-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)
-    end
-    Sparspak.SparseSolver.solve!(cache.cacheval)
-    copyto!(cache.u, cache.cacheval.p.x)
-    SciMLBase.build_linear_solution(alg, cache.u, nothing, cache)
-end

test/runtests.jl

@@ -22,6 +22,7 @@ if GROUP == "All" || GROUP == "Core"
     @time @safetestset "Basic Tests" begin include("basictests.jl") end
     @time @safetestset "Zero Initialization Tests" begin include("zeroinittests.jl") end
     @time @safetestset "Non-Square Tests" begin include("nonsquare.jl") end
+    @time @safetestset "SparseVector b Tests" begin include("sparse_vector.jl") end
 end
 
 if GROUP == "LinearSolveCUDA"

test/sparse_vector.jl

@@ -0,0 +1,43 @@
+using SparseArrays
+using LinearSolve
+using LinearAlgebra
+
+# Constructing sparse array 
+function hess_sparse(x::Vector{T}) where {T}
+    return [
+        -sin(x[1] + x[2]) + 1,
+        -sin(x[1] + x[2]),
+        -sin(x[1] + x[2]),
+        -sin(x[1] + x[2]) + 1.0,
+        1.0,
+        1.0,
+        12.0 * x[5]^2 + 1.0,
+        1.0,
+    ]
+end
+rowval = [1, 1, 2, 2, 3, 4, 5, 6]
+colval = [1, 2, 1, 2, 3, 4, 5, 6]
+
+# Constructing sparse vec
+function grad_sparse(x::Vector{T}) where {T <: Number}
+    return [cos(x[1] + x[2]), cos(x[1] + x[2]), 2 * x[3], 1 / 2, 4 * x[5]^3, 1 / 2]
+end
+gradinds = [1, 2, 3, 4, 5, 6]
+
+# Forming the matrix and vector
+x0 = [
+    0.7853981648713337,
+    0.7853981693418342,
+    1.023999999999997e-7,
+    -1.0,
+    0.33141395338218227,
+    -1.0,
+]
+n = length(x0)
+hess_mat = sparse(rowval, colval, hess_sparse(x0), n, n)
+grad_vec = sparsevec(gradinds, grad_sparse(x0), n)
+
+# # Converting grad_vec to dense succeds in solving
+prob = LinearProblem(hess_mat, grad_vec)
+linsolve = init(prob)
+@test solve(linsolve).u ≈ hess_mat \ Array(grad_vec)