Update the wrapper of Krylov.jl

Author: amontoison

docs/src/solvers/solvers.md

@@ -178,10 +178,12 @@ IterativeSolversJL(args...;
 
 ### Krylov.jl
 
-- `KrylovJL_CG(args...;kwargs...)`: A generic CG implementation
-- `KrylovJL_GMRES(args...;kwargs...)`: A generic GMRES implementation
-- `KrylovJL_BICGSTAB(args...;kwargs...)`: A generic BICGSTAB implementation
-- `KrylovJL_MINRES(args...;kwargs...)`: A generic MINRES implementation
+- `KrylovJL_CG(args...;kwargs...)`: A generic CG implementation for Hermitian and positive definite linear systems
+- `KrylovJL_MINRES(args...;kwargs...)`: A generic MINRES implementation for Hermitian linear systems
+- `KrylovJL_GMRES(args...;kwargs...)`: A generic GMRES implementation for square non-Hermitian linear systems
+- `KrylovJL_BICGSTAB(args...;kwargs...)`: A generic BICGSTAB implementation for square non-Hermitian linear systems
+- `KrylovJL_LSMR(args...;kwargs...)`: A generic LSMR implementation for least-squares problems
+- `KrylovJL_CRAIGMR(args...;kwargs...)`: A generic CRAIGMR implementation for least-norm problems
 
 The general algorithm is:
 

src/LinearSolve.jl

@@ -74,9 +74,10 @@ export LUFactorization, SVDFactorization, QRFactorization, GenericFactorization,
 
 export LinearSolveFunction
 
-export KrylovJL, KrylovJL_CG, KrylovJL_GMRES, KrylovJL_BICGSTAB, KrylovJL_MINRES,
+export KrylovJL, KrylovJL_CG, KrylovJL_MINRES, KrylovJL_GMRES,
+       KrylovJL_BICGSTAB, KrylovJL_LSMR, KrylovJL_CRAIGMR,
        IterativeSolversJL, IterativeSolversJL_CG, IterativeSolversJL_GMRES,
        IterativeSolversJL_BICGSTAB, IterativeSolversJL_MINRES,
-       KrylovKitJL, KrylovKitJL_CG, KrylovKitJL_GMRES, KrylovJL_LSMR
+       KrylovKitJL, KrylovKitJL_CG, KrylovKitJL_GMRES
 
 end

src/default.jl

@@ -71,7 +71,12 @@ end
 
 function defaultalg(A::SciMLBase.AbstractDiffEqOperator, b,
                     assumptions::OperatorAssumptions{false})
-    KrylovJL_LSMR()
+    m,n = size(A)
+    if m < n
+        KrylovJL_CRAIGMR()
+    else
+        KrylovJL_LSMR()
+    end
 end
 
 # Handle ambiguity

src/iterative_wrappers.jl

@@ -15,20 +15,30 @@ function KrylovJL(args...; KrylovAlg = Krylov.gmres!,
                     args, kwargs)
 end
 
-KrylovJL_CG(args...; kwargs...) = KrylovJL(args...; KrylovAlg = Krylov.cg!, kwargs...)
-KrylovJL_GMRES(args...; kwargs...) = KrylovJL(args...; KrylovAlg = Krylov.gmres!, kwargs...)
-KrylovJL_FGMRES(args...; kwargs...) = KrylovJL(args...; KrylovAlg = Krylov.fgmres!, kwargs...)
-function KrylovJL_BICGSTAB(args...; kwargs...)
-    KrylovJL(args...; KrylovAlg = Krylov.bicgstab!, kwargs...)
+function KrylovJL_CG(args...; kwargs...)
+    KrylovJL(args...; KrylovAlg = Krylov.cg!, kwargs...)
 end
+
 function KrylovJL_MINRES(args...; kwargs...)
     KrylovJL(args...; KrylovAlg = Krylov.minres!, kwargs...)
 end
 
+function KrylovJL_GMRES(args...; kwargs...)
+    KrylovJL(args...; KrylovAlg = Krylov.gmres!, kwargs...)
+end
+
+function KrylovJL_BICGSTAB(args...; kwargs...)
+    KrylovJL(args...; KrylovAlg = Krylov.bicgstab!, kwargs...)
+end
+
 function KrylovJL_LSMR(args...; kwargs...)
     KrylovJL(args...; KrylovAlg = Krylov.lsmr!, kwargs...)
 end
 
+function KrylovJL_CRAIGMR(args...; kwargs...)
+    KrylovJL(args...; KrylovAlg = Krylov.craigmr!, kwargs...)
+end
+
 function get_KrylovJL_solver(KrylovAlg)
     KS = if (KrylovAlg === Krylov.lsmr!)
         Krylov.LsmrSolver

test/nonsquare.jl

@@ -3,18 +3,24 @@ using SparseArrays, LinearAlgebra
 
 m, n = 13, 3
 
-A = rand(m, n);
-b = rand(m);
+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);
+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
+
+A = sprand(n, m, 0.5)
+b = rand(n)
+prob = LinearProblem(A, b)
+res = A \ b
+@test solve(prob, KrylovJL_CRAIGMR()) ≈ res