Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

reinit! calls precs three times #527

Closed
j-fu opened this issue Aug 19, 2024 · 5 comments · Fixed by #547
Closed

reinit! calls precs three times #527

j-fu opened this issue Aug 19, 2024 · 5 comments · Fixed by #547
Assignees
Labels
bug Something isn't working

Comments

@j-fu
Copy link
Contributor

j-fu commented Aug 19, 2024

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.

@j-fu j-fu added the bug Something isn't working label Aug 19, 2024
@oscardssmith
Copy link
Contributor

good catch. I added the reinit before the set property so I forgot about the interaction. fix incoming

oscardssmith added a commit to oscardssmith/LinearSolve.jl that referenced this issue Aug 19, 2024
@oscardssmith
Copy link
Contributor

a reuse_precs kwarg seems reasonable for this function.

@ChrisRackauckas
Copy link
Member

Why would it actually make the new prec instead of just setting a flag to make a prec and then make at the next solve step?

@oscardssmith
Copy link
Contributor

oh, we could possibly used cache.isfresh for this... that's interesting

@ChrisRackauckas
Copy link
Member

Yes because you only have to update the precs when A changes, which is when isfresh

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
3 participants