Clean up sparse handling

- Fix case with sparse b, fixes #248
- Removes Sparspak from the section off code that's only for non-GPL cases. Sparspak.jl isn't GPL-licensed so we can always use it.
- Restructured the default algorithm so that Sparspak is used when GPL is disabled, making it so it can always factorize properly.

Author: ChrisRackauckas

src/common.jl

@@ -111,7 +111,14 @@ function SciMLBase.init(prob::LinearProblem, alg::Union{SciMLLinearSolveAlgorith
     Tc = typeof(cacheval)
 
     A = alias_A ? A : deepcopy(A)
-    b = alias_b ? b : deepcopy(b)
+    b = if b isa SparseArray && !(A isa Diagonal)
+            Array(b) # the solution to a linear solve will always be dense!
+        elseif alias_b 
+            b 
+        else 
+            deepcopy(b)
+        end
+    end
 
     cache = LinearCache{
                         typeof(A),

src/default.jl

@@ -152,6 +152,20 @@ function defaultalg(A, b, ::OperatorAssumptions{true})
     alg
 end
 
+function defaultalg(A::SparseMatrixCSC, b, ::OperatorAssumptions{true})
+    # If GPL libraries are loaded, then use SuiteSparse. Otherwise Sparspak
+    if INCLUDE_SPARSE 
+        if length(b) <= 10000
+            alg = KLUFactorization()
+        else
+            alg = UMFPACKFactorization()
+        end
+    else
+        alg = SparspakFactorization()
+    end
+    alg
+end
+
 function defaultalg(A, b, ::OperatorAssumptions{false})
     QRFactorization()
 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

@@ -5,32 +5,4 @@ function _ldiv!(x::Vector,
                          SuiteSparse.SPQR.QRSparse,
                          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
+end