From 9d9b8cec1792087e175a978b7400a8d4fe6256c4 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Fri, 25 Nov 2022 15:47:54 -0500 Subject: [PATCH 1/3] Update the wrapper of Krylov.jl --- Project.toml | 2 +- src/iterative_wrappers.jl | 7 +++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a8080343a..cc3d4a7f6 100644 --- 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/src/iterative_wrappers.jl b/src/iterative_wrappers.jl index 1b67d17a4..1fabea2e3 100644 --- a/src/iterative_wrappers.jl +++ b/src/iterative_wrappers.jl @@ -17,6 +17,7 @@ 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...) end @@ -89,6 +90,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 +110,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! || From d9ec11219e8415f13a7da0d5529dcc26c6a79f79 Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Fri, 25 Nov 2022 21:09:57 -0500 Subject: [PATCH 2/3] Update the wrapper of Krylov.jl --- docs/src/solvers/solvers.md | 10 ++++++---- src/LinearSolve.jl | 5 +++-- src/default.jl | 7 ++++++- src/iterative_wrappers.jl | 20 +++++++++++++++----- test/nonsquare.jl | 14 ++++++++++---- 5 files changed, 40 insertions(+), 16 deletions(-) diff --git a/docs/src/solvers/solvers.md b/docs/src/solvers/solvers.md index 02cfa42f9..a23db41cb 100644 --- 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 index 10db88a73..1945bde48 100644 --- 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 index bd61f22e4..12c2ada32 100644 --- 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 index 1fabea2e3..f5a56f767 100644 --- a/src/iterative_wrappers.jl +++ b/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 diff --git a/test/nonsquare.jl b/test/nonsquare.jl index 5c9057fd6..10f2a817b 100644 --- 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 = A \ b +@test solve(prob, KrylovJL_CRAIGMR()) ≈ res From 127970d093d326845dd2a265d20dfb41f1546a1f Mon Sep 17 00:00:00 2001 From: Alexis Montoison Date: Sat, 26 Nov 2022 15:13:49 -0500 Subject: [PATCH 3/3] Fix issue in nonsquare.jl --- src/default.jl | 2 +- test/nonsquare.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/default.jl b/src/default.jl index 12c2ada32..f316e7f2a 100644 --- a/src/default.jl +++ b/src/default.jl @@ -71,7 +71,7 @@ end function defaultalg(A::SciMLBase.AbstractDiffEqOperator, b, assumptions::OperatorAssumptions{false}) - m,n = size(A) + m, n = size(A) if m < n KrylovJL_CRAIGMR() else diff --git a/test/nonsquare.jl b/test/nonsquare.jl index 10f2a817b..a555bc54f 100644 --- a/test/nonsquare.jl +++ b/test/nonsquare.jl @@ -22,5 +22,5 @@ res = A \ b A = sprand(n, m, 0.5) b = rand(n) prob = LinearProblem(A, b) -res = A \ b +res = Matrix(A) \ b @test solve(prob, KrylovJL_CRAIGMR()) ≈ res