Merge pull request #188 from SciML/nonsquare_operators

Default non-square operators to LSMR

Author: ChrisRackauckas

src/LinearSolve.jl

@@ -58,6 +58,6 @@ export LinearSolveFunction
 export KrylovJL, KrylovJL_CG, KrylovJL_GMRES, KrylovJL_BICGSTAB, KrylovJL_MINRES,
        IterativeSolversJL, IterativeSolversJL_CG, IterativeSolversJL_GMRES,
        IterativeSolversJL_BICGSTAB, IterativeSolversJL_MINRES,
-       KrylovKitJL, KrylovKitJL_CG, KrylovKitJL_GMRES
+       KrylovKitJL, KrylovKitJL_CG, KrylovKitJL_GMRES, KrylovJL_LSMR
 
 end

src/default.jl

@@ -69,6 +69,11 @@ function defaultalg(A::SciMLBase.AbstractDiffEqOperator, b,
     KrylovJL_GMRES()
 end
 
+function defaultalg(A::SciMLBase.AbstractDiffEqOperator, b,
+                    assumptions::OperatorAssumptions{false})
+    KrylovJL_LSMR()
+end
+
 # Handle ambiguity
 function defaultalg(A::GPUArraysCore.AbstractGPUArray, b::GPUArraysCore.AbstractGPUArray,
                     ::OperatorAssumptions{true})

src/factorization.jl

@@ -5,6 +5,14 @@ function _ldiv!(x::Vector, A::Factorization, b::Vector)
     ldiv!(A, x)
 end
 
+# Specialize QR for the non-square case
+# Missing ldiv! definitions: https://github.com/JuliaSparse/SparseArrays.jl/issues/242
+function _ldiv!(x::Vector,
+                A::Union{SparseArrays.QR, LinearAlgebra.QRCompactWY,
+                         SuiteSparse.SPQR.QRSparse}, b::Vector)
+    x .= A \ b
+end
+
 function SciMLBase.solve(cache::LinearCache, alg::AbstractFactorization; kwargs...)
     if cache.isfresh
         fact = do_factorization(alg, cache.A, cache.b, cache.u)
@@ -90,7 +98,7 @@ end
 
 function do_factorization(alg::QRFactorization, A, b, u)
     A = convert(AbstractMatrix, A)
-    if alg.inplace
+    if alg.inplace && (!(A isa SparseMatrixCSC) || VERSION >= v"1.8-")
         fact = qr!(A, alg.pivot)
     else
         fact = qr(A) # CUDA.jl does not allow other args!

src/iterative_wrappers.jl

@@ -24,6 +24,10 @@ function KrylovJL_MINRES(args...; kwargs...)
     KrylovJL(args...; KrylovAlg = Krylov.minres!, kwargs...)
 end
 
+function KrylovJL_LSMR(args...; kwargs...)
+    KrylovJL(args...; KrylovAlg = Krylov.lsmr!, kwargs...)
+end
+
 function get_KrylovJL_solver(KrylovAlg)
     KS = if (KrylovAlg === Krylov.lsmr!)
         Krylov.LsmrSolver

test/nonsquare.jl

@@ -0,0 +1,20 @@
+using LinearSolve, Test
+using SparseArrays, LinearAlgebra
+
+m, n = 13, 3
+
+A = rand(m, n);
+b = rand(m);
+prob = LinearProblem(A, b)
+res = A \ b
+@test solve(prob).u ≈ res
+@test solve(prob, QRFactorization()) ≈ res
+@test solve(prob, KrylovJL_LSMR()) ≈ res
+
+A = sprand(m, n, 0.5);
+b = rand(m);
+prob = LinearProblem(A, b)
+res = A \ b
+@test solve(prob).u ≈ res
+@test solve(prob, QRFactorization()) ≈ res
+@test solve(prob, KrylovJL_LSMR()) ≈ res

test/runtests.jl

@@ -21,6 +21,7 @@ end
 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
 end
 
 if GROUP == "LinearSolveCUDA"