Skip to content

reinit! calls precs three times #527

Closed
@j-fu

Description

@j-fu

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.

Metadata

Metadata

Assignees

Labels

bugSomething isn't working

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions