Skip to content

Commit aad2263

Browse files
Fix default Krylov preconditioner dispatch
Dispatch on the selected inner Krylov workspace so DefaultLinearSolver forwards supplied left and right preconditioners. Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com> Co-Authored-By: Claude <noreply@anthropic.com> Claude-Session: https://chatgpt.com/codex/tasks/01a03a17-ad6f-7131-82fc-d0fd57ea6512
1 parent 5dcf04d commit aad2263

2 files changed

Lines changed: 38 additions & 12 deletions

File tree

src/iterative_wrappers.jl

Lines changed: 12 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -680,41 +680,41 @@ function SciMLBase.solve!(cache::LinearCache, alg::KrylovJL; kwargs...)
680680
ldiv = true, history = true, filtered_kwargs...,
681681
)
682682

683-
if cache.cacheval isa Krylov.CgWorkspace
683+
if cacheval isa Krylov.CgWorkspace
684684
N !== I &&
685685
@SciMLMessage(
686686
"$(alg.KrylovAlg) doesn't support right preconditioning.",
687687
verbose, :no_right_preconditioning
688688
)
689689
Krylov.krylov_solve!(args...; M, kwargs...)
690-
elseif cache.cacheval isa Krylov.GmresWorkspace
690+
elseif cacheval isa Krylov.GmresWorkspace
691691
Krylov.krylov_solve!(args...; M, N, restart = alg.gmres_restart > 0, kwargs...)
692-
elseif cache.cacheval isa Krylov.FgmresWorkspace
692+
elseif cacheval isa Krylov.FgmresWorkspace
693693
Krylov.krylov_solve!(args...; M, N, kwargs...)
694-
elseif cache.cacheval isa Krylov.BicgstabWorkspace
694+
elseif cacheval isa Krylov.BicgstabWorkspace
695695
Krylov.krylov_solve!(args...; M, N, kwargs...)
696-
elseif cache.cacheval isa Krylov.MinresWorkspace
696+
elseif cacheval isa Krylov.MinresWorkspace
697697
N !== I &&
698698
@SciMLMessage(
699699
"$(alg.KrylovAlg) doesn't support right preconditioning.",
700700
verbose, :no_right_preconditioning
701701
)
702702
Krylov.krylov_solve!(args...; M, kwargs...)
703-
elseif cache.cacheval isa Krylov.BlockGmresWorkspace
703+
elseif cacheval isa Krylov.BlockGmresWorkspace
704704
Krylov.krylov_solve!(args...; M, N, restart = alg.gmres_restart > 0, kwargs...)
705-
elseif cache.cacheval isa Krylov.BlockMinresWorkspace
705+
elseif cacheval isa Krylov.BlockMinresWorkspace
706706
N !== I &&
707707
@SciMLMessage(
708708
"$(alg.KrylovAlg) doesn't support right preconditioning.",
709709
verbose, :no_right_preconditioning
710710
)
711711
Krylov.krylov_solve!(args...; M, kwargs...)
712-
elseif cache.cacheval isa Krylov.LsmrWorkspace ||
713-
cache.cacheval isa Krylov.LsqrWorkspace ||
714-
cache.cacheval isa Krylov.LslqWorkspace
712+
elseif cacheval isa Krylov.LsmrWorkspace ||
713+
cacheval isa Krylov.LsqrWorkspace ||
714+
cacheval isa Krylov.LslqWorkspace
715715
Krylov.krylov_solve!(args...; M, N, kwargs...)
716-
elseif cache.cacheval isa Krylov.CglsWorkspace ||
717-
cache.cacheval isa Krylov.CrlsWorkspace
716+
elseif cacheval isa Krylov.CglsWorkspace ||
717+
cacheval isa Krylov.CrlsWorkspace
718718
N !== I &&
719719
@SciMLMessage(
720720
"$(alg.KrylovAlg) doesn't support right preconditioning.",

test/Core/default_algs.jl

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,15 @@
11
using LinearSolve, RecursiveFactorization, LinearAlgebra, SparseArrays, Test
22
using SciMLOperators: FunctionOperator, MatrixOperator, WOperator, has_concretization
33

4+
struct CountingIdentityPreconditioner
5+
calls::Base.RefValue{Int}
6+
end
7+
8+
function LinearAlgebra.ldiv!(y, P::CountingIdentityPreconditioner, x)
9+
P.calls[] += 1
10+
return copyto!(y, x)
11+
end
12+
413
@test LinearSolve.defaultalg(nothing, zeros(3)).alg === LinearSolve.DefaultAlgorithmChoice.GenericLUFactorization
514
prob = LinearProblem(rand(3, 3), rand(3))
615
solve(prob)
@@ -710,3 +719,20 @@ end
710719
LinearSolve.DefaultAlgorithmChoice.LHLFactorization
711720
@test solve(LinearProblem(Wdense, bw)).u ref rtol = 1.0e-8
712721
end
722+
723+
@testset "default matrix-free solver uses supplied preconditioners" begin
724+
n = 4
725+
Jm = rand(n, n) + n * I
726+
b = rand(n)
727+
fj(v, u, p, t) = Jm * v
728+
fj(w, v, u, p, t) = mul!(w, Jm, v)
729+
Jfree = FunctionOperator(fj, zeros(n), zeros(n); islinear = true)
730+
Wfree = WOperator{true}(I, 0.1, Jfree, zeros(n))
731+
Pl = CountingIdentityPreconditioner(Ref(0))
732+
Pr = CountingIdentityPreconditioner(Ref(0))
733+
734+
solve(LinearProblem(Wfree, b), nothing; Pl, Pr)
735+
736+
@test Pl.calls[] > 0
737+
@test Pr.calls[] > 0
738+
end

0 commit comments

Comments
 (0)