Skip to content

Commit 53fcf8f

Browse files
authored
Stop cg_lanczos! resolving to the shifted workspace (#1251)
`get_KrylovJL_solver` tested `KrylovAlg === Krylov.cg_lanczos!` twice. The first branch returned `CgLanczosShiftWorkspace`, so it won, and the second branch returning `CgLanczosWorkspace` was unreachable. `cg_lanczos!` is listed in the `KrylovJL` docstring as a supported `KrylovAlg`, so this is reachable from documented usage: solve(LinearProblem(A, b), KrylovJL(KrylovAlg = Krylov.cg_lanczos!)) MethodError: no method matching CgLanczosShiftWorkspace(::Matrix{Float64}, ::Vector{Float64}) Dropping the first branch lets it fall through to the right workspace. The solve then succeeds, to a relative error of 5.9e-9 on a symmetric positive definite system. `cg_lanczos_shift!` and `cgls_lanczos_shift!` are deliberately left unmapped. Both exist in Krylov 0.10, but they take a `shifts` argument that this wrapper never passes, and neither appears in the docstring's list, so they keep hitting "Invalid Krylov method detected" rather than being advertised and broken. Spotted in #554, which is otherwise superseded: block GMRES and MINRES are already interfaced through `BlockGmresWorkspace` and `BlockMinresWorkspace` for a matrix right-hand side.
1 parent 7d64eae commit 53fcf8f

2 files changed

Lines changed: 30 additions & 2 deletions

File tree

src/iterative_wrappers.jl

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -373,8 +373,6 @@ function get_KrylovJL_solver(KrylovAlg)
373373
Krylov.CrmrWorkspace
374374
elseif (KrylovAlg === Krylov.cg!)
375375
Krylov.CgWorkspace
376-
elseif (KrylovAlg === Krylov.cg_lanczos!)
377-
Krylov.CgLanczosShiftWorkspace
378376
elseif (KrylovAlg === Krylov.cgls!)
379377
Krylov.CglsWorkspace
380378
elseif (KrylovAlg === Krylov.cg_lanczos!)

test/Core/nonsquare.jl

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -345,3 +345,33 @@ end
345345
end
346346
end
347347
end
348+
349+
# `get_KrylovJL_solver` matched `cg_lanczos!` twice. The first branch handed back the
350+
# *shifted* workspace, so the documented `KrylovAlg = Krylov.cg_lanczos!` threw a
351+
# `MethodError` out of the workspace constructor, and the second branch, the correct
352+
# one, was unreachable. See SciML/LinearSolve.jl#554.
353+
@testset "documented KrylovAlg values map to their own workspace" begin
354+
for (f, W) in (
355+
(Krylov.cg!, Krylov.CgWorkspace),
356+
(Krylov.cg_lanczos!, Krylov.CgLanczosWorkspace),
357+
(Krylov.cgls!, Krylov.CglsWorkspace),
358+
(Krylov.minres!, Krylov.MinresWorkspace),
359+
(Krylov.gmres!, Krylov.GmresWorkspace),
360+
)
361+
@test LinearSolve.get_KrylovJL_solver(f) === W
362+
end
363+
364+
# The shifted variants need a `shifts` argument that this wrapper never passes, and
365+
# the docstring does not list them, so they stay unmapped rather than advertised.
366+
@test_throws ErrorException LinearSolve.get_KrylovJL_solver(Krylov.cg_lanczos_shift!)
367+
368+
# `cg_lanczos!` wants a symmetric positive definite system.
369+
Random.seed!(554)
370+
nl = 24
371+
Xl = rand(nl, nl)
372+
Al = Xl' * Xl + nl * I
373+
bl = rand(nl)
374+
sol = solve(LinearProblem(Al, bl), KrylovJL(KrylovAlg = Krylov.cg_lanczos!))
375+
@test SciMLBase.successful_retcode(sol)
376+
@test sol.u Al \ bl rtol = 1.0e-6
377+
end

0 commit comments

Comments
 (0)