Merge pull request #234 from amontoison/update_krylov

Update the wrapper of Krylov.jl
SciML · Nov 27, 2022 · 8c63893 · 8c63893
2 parents fb4bb2d + 127970d
commit 8c63893
Show file tree

Hide file tree

Showing 6 changed files with 47 additions and 16 deletions.
diff --git a/Project.toml b/Project.toml
@@ -29,7 +29,7 @@ FastLapackInterface = "1"
 GPUArraysCore = "0.1"
 IterativeSolvers = "0.9.2"
 KLU = "0.3.0, 0.4"
-Krylov = "0.7.11, 0.8"
+Krylov = "0.9"
 KrylovKit = "0.5, 0.6"
 RecursiveFactorization = "0.2.8"
 Reexport = "1"

diff --git a/docs/src/solvers/solvers.md b/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:
 

diff --git a/src/LinearSolve.jl b/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
diff --git a/src/default.jl b/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

diff --git a/src/iterative_wrappers.jl b/src/iterative_wrappers.jl
@@ -15,19 +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...)
-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
@@ -89,6 +100,10 @@ function get_KrylovJL_solver(KrylovAlg)
         Krylov.QmrSolver
     elseif (KrylovAlg === Krylov.gmres!)
         Krylov.GmresSolver
+    elseif (KrylovAlg === Krylov.fgmres!)
+        Krylov.FgmresSolver
+    elseif (KrylovAlg === Krylov.gpmr!)
+        Krylov.GpmrSolver
     elseif (KrylovAlg === Krylov.fom!)
         Krylov.FomSolver
     end
@@ -105,6 +120,8 @@ function init_cacheval(alg::KrylovJL, A, b, u, Pl, Pr, maxiters::Int, abstol, re
     solver = if (alg.KrylovAlg === Krylov.dqgmres! ||
                  alg.KrylovAlg === Krylov.diom! ||
                  alg.KrylovAlg === Krylov.gmres! ||
+                 alg.KrylovAlg === Krylov.fgmres! ||
+                 alg.KrylovAlg === Krylov.gpmr! ||
                  alg.KrylovAlg === Krylov.fom!)
         KS(A, b, memory)
     elseif (alg.KrylovAlg === Krylov.minres! ||

diff --git a/test/nonsquare.jl b/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 = Matrix(A) \ b
+@test solve(prob, KrylovJL_CRAIGMR()) ≈ res