Closed
Description
Describe the bug 🐞
In reinit!
. preconditioner construction via precs is called several times: once in the method itself, and then by cache.A=A
and cache.p=p
(triggered by setproperty
).
Expected behavior
I would expect that the preconditioner is updated only once.
Minimal Reproducible Example 👇
using ILUZero, LinearSolve, SparseArrays, LinearAlgebra
struct ILUZero_ILU0 end
function (::ILUZero_ILU0)(A, p=nothing)
println("new prec")
(ilu0(SparseMatrixCSC(A)),I)
end
n=100
u0=ones(n)
println("Original problem:")
A=spdiagm(1 => fill(-1.0, n - 1), 0 => fill(100.0,n), -1 => fill(-1.0, n - 1))
pr=LinearProblem(A,A*u0)
solver=KrylovJL_CG(precs=ILUZero_ILU0(), ldiv=true)
cache=init(pr,solver)
u=solve!(cache)
@show norm(u-u0)
println("Updated problem:")
for i=1:size(A,2)
A[i,i]+=100
end
reinit!(cache;A,b=A*u0)
solve!(cache)
@show norm(u-u0)
Output
Original problem:
new prec
norm(u - u0) = 0.0
Updated problem:
new prec
new prec
new prec
norm(u - u0) = 0.0
julia> pkgversion(LinearSolve)
v"2.33.0"
julia> VERSION
v"1.11.0-rc2"
Discussion
I think that the setproperty
overload for cache.A
and cache.p
is not
necessary. In any case I find it a bit dangerous.
As a feature request, I would like to see a kwarg keep_precs
, allowing
to reinit!
without updating Pl, Pr in a transparent manner.
I could try a corresponding PR once this is fixed.