Fix non-square sparse matrix preallocations for defaults

ext/LinearSolveSparseArraysExt.jl

@@ -27,13 +27,6 @@ function LinearSolve.init_cacheval(alg::RFLUFactorization,
     nothing, nothing
 end
 
-function LinearSolve.init_cacheval(
-        alg::QRFactorization, A::Symmetric{<:Number, <:SparseMatrixCSC}, b, u, Pl, Pr,
-        maxiters::Int, abstol, reltol, verbose::Bool,
-        assumptions::OperatorAssumptions)
-    return nothing
-end
-
 function LinearSolve.handle_sparsematrixcsc_lu(A::AbstractSparseMatrixCSC)
     lu(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A), nonzeros(A)),
         check = false)
@@ -71,22 +64,43 @@ const PREALLOCATED_UMFPACK = SparseArrays.UMFPACK.UmfpackLU(SparseMatrixCSC(0, 0
     Int[], Float64[]))
 
 function LinearSolve.init_cacheval(
-        alg::UMFPACKFactorization, A::SparseMatrixCSC{Float64, Int}, b, u,
+        alg::Union{LUFactorization, GenericLUFactorization}, A::AbstractSparseArray{<:Number, <:Integer}, b, u,
+        Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Bool, assumptions::OperatorAssumptions)
+    nothing
+end
+
+function LinearSolve.init_cacheval(
+        alg::UMFPACKFactorization, A::AbstractArray, b, u,
+        Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Bool, assumptions::OperatorAssumptions)
+    nothing
+end
+
+function LinearSolve.init_cacheval(
+        alg::Union{LUFactorization, GenericLUFactorization}, A::AbstractSparseArray{Float64, Int64}, b, u,
         Pl, Pr,
         maxiters::Int, abstol, reltol,
         verbose::Bool, assumptions::OperatorAssumptions)
     PREALLOCATED_UMFPACK
 end
 
 function LinearSolve.init_cacheval(
-        alg::UMFPACKFactorization, A::AbstractSparseArray{Float64}, b, u, Pl, Pr,
+        alg::Union{LUFactorization, GenericLUFactorization}, A::AbstractSparseArray{Float64, Int32}, b, u,
+        Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Bool, assumptions::OperatorAssumptions)
+    SparseArrays.UMFPACK.UmfpackLU(SparseMatrixCSC{Float64, Int32}(0, 0, [Int32(1)], Int32[], Float64[]))
+end
+
+function LinearSolve.init_cacheval(
+        alg::UMFPACKFactorization, A::AbstractSparseArray{Float64, Int}, b, u, Pl, Pr,
         maxiters::Int, abstol,
         reltol,
         verbose::Bool, assumptions::OperatorAssumptions)
-    A = convert(AbstractMatrix, A)
-    zerobased = SparseArrays.getcolptr(A)[1] == 0
-    return SparseArrays.UMFPACK.UmfpackLU(SparseMatrixCSC(size(A)..., getcolptr(A),
-        rowvals(A), nonzeros(A)))
+    PREALLOCATED_UMFPACK
 end
 
 function SciMLBase.solve!(
@@ -129,21 +143,27 @@ const PREALLOCATED_KLU = KLU.KLUFactorization(SparseMatrixCSC(0, 0, [1], Int[],
     Float64[]))
 
 function LinearSolve.init_cacheval(
-        alg::KLUFactorization, A::SparseMatrixCSC{Float64, Int}, b, u, Pl,
+        alg::KLUFactorization, A::AbstractArray, b, u, Pl,
         Pr,
         maxiters::Int, abstol, reltol,
         verbose::Bool, assumptions::OperatorAssumptions)
+    nothing
+end
+
+function LinearSolve.init_cacheval(
+        alg::KLUFactorization, A::AbstractSparseArray{Float64, Int64}, b, u, Pl, Pr,
+        maxiters::Int, abstol,
+        reltol,
+        verbose::Bool, assumptions::OperatorAssumptions)
     PREALLOCATED_KLU
 end
 
 function LinearSolve.init_cacheval(
-        alg::KLUFactorization, A::AbstractSparseArray{Float64}, b, u, Pl, Pr,
+        alg::KLUFactorization, A::AbstractSparseArray{Float64, Int32}, b, u, Pl, Pr,
         maxiters::Int, abstol,
         reltol,
         verbose::Bool, assumptions::OperatorAssumptions)
-    A = convert(AbstractMatrix, A)
-    return KLU.KLUFactorization(SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),
-        nonzeros(A)))
+    KLU.KLUFactorization(SparseMatrixCSC{Float64, Int32}(0, 0, [Int32(1)], Int32[], Float64[]))
 end
 
 # TODO: guard this against errors
@@ -249,6 +269,34 @@ function LinearSolve.defaultalg(
     end
 end
 
+# SPQR Handling
+function LinearSolve.init_cacheval(
+        alg::QRFactorization, A::AbstractSparseArray{<:Number, <:Integer}, b, u,
+        Pl, Pr,
+        maxiters::Int, abstol, reltol,
+        verbose::Bool, assumptions::OperatorAssumptions)
+    nothing
+end
+
+function LinearSolve.init_cacheval(alg::QRFactorization, A::SparseMatrixCSC{Float64, Int}, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol, verbose::Bool,
+        assumptions::OperatorAssumptions)
+    LinearSolve.ArrayInterface.qr_instance(convert(AbstractMatrix, A), alg.pivot)
+end
+
+function LinearSolve.init_cacheval(alg::QRFactorization, A::SparseMatrixCSC{Float64, Int32}, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol, verbose::Bool,
+        assumptions::OperatorAssumptions)
+    LinearSolve.ArrayInterface.qr_instance(convert(AbstractMatrix, A), alg.pivot)
+end
+
+function LinearSolve.init_cacheval(
+        alg::QRFactorization, A::Symmetric{<:Number, <:SparseMatrixCSC}, b, u, Pl, Pr,
+        maxiters::Int, abstol, reltol, verbose::Bool,
+        assumptions::OperatorAssumptions)
+    return nothing
+end
+
 LinearSolve.PrecompileTools.@compile_workload begin
     A = sprand(4, 4, 0.3) + I
     b = rand(4)

ext/LinearSolveSparspakExt.jl

@@ -17,18 +17,23 @@ function LinearSolve.init_cacheval(
 end
 
 function LinearSolve.init_cacheval(
-        ::SparspakFactorization, A::AbstractSparseMatrixCSC, b, u, Pl, Pr, maxiters::Int, abstol,
+        ::SparspakFactorization, A::AbstractSparseMatrixCSC{Tv, Ti}, b, u, Pl, Pr, maxiters::Int, abstol,
         reltol,
-        verbose::Bool, assumptions::OperatorAssumptions)
-    A = convert(AbstractMatrix, A)
-    if A isa SparseArrays.AbstractSparseArray
-        return sparspaklu(
-            SparseMatrixCSC(size(A)..., getcolptr(A), rowvals(A),
-                nonzeros(A)),
-            factorize = false)
+        verbose::Bool, assumptions::OperatorAssumptions) where {Tv, Ti}
+
+    if size(A,1) == size(A,2)
+        A = convert(AbstractMatrix, A)
+        if A isa SparseArrays.AbstractSparseArray
+            return sparspaklu(
+                SparseMatrixCSC{Tv, Ti}(size(A)..., getcolptr(A), rowvals(A),
+                    nonzeros(A)),
+                factorize = false)
+        else
+            return sparspaklu(SparseMatrixCSC(0, 0, [one(Ti)], Ti[], eltype(A)[]),
+                factorize = false)
+        end
     else
-        return sparspaklu(SparseMatrixCSC(0, 0, [1], Int[], eltype(A)[]),
-            factorize = false)
+        PREALLOCATED_SPARSEPAK
     end
 end
 

test/default_algs.jl

@@ -144,3 +144,31 @@ cache.A = [2.0 1.0
 sol = solve!(cache)
 
 @test !SciMLBase.successful_retcode(sol.retcode)
+
+## Non-square Sparse Defaults 
+# https://github.com/SciML/NonlinearSolve.jl/issues/599
+A = SparseMatrixCSC{Float64, Int64}([
+    1.0 0.0
+    1.0 1.0
+])
+b = ones(2)
+A2 = hcat(A,A)
+prob = LinearProblem(A, b)
+@test SciMLBase.successful_retcode(solve(prob))
+
+prob2 = LinearProblem(A2, b)
+@test SciMLBase.successful_retcode(solve(prob2))
+
+A = SparseMatrixCSC{Float64, Int32}([
+    1.0 0.0
+    1.0 1.0
+])
+b = ones(2)
+A2 = hcat(A,A)
+prob = LinearProblem(A, b)
+@test SciMLBase.successful_retcode(solve(prob))
+
+@info "This test"
+
+prob2 = LinearProblem(A2, b)
+@test SciMLBase.successful_retcode(solve(prob2))