SciML / LinearSolve.jl

LinearSolve.jl: High-Performance Unified Interface for Linear Solvers in Julia. Easily switch between factorization and Krylov methods, add preconditioners, and all in one interface.
https://docs.sciml.ai/LinearSolve/stable/
Other
244 stars 52 forks source link

`reinit!` calls `precs` three times #527

Open j-fu opened 3 weeks ago

j-fu commented 3 weeks ago

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.

oscardssmith commented 3 weeks ago

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

oscardssmith commented 3 weeks ago

a reuse_precs kwarg seems reasonable for this function.

ChrisRackauckas commented 3 weeks ago

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 commented 3 weeks ago

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

ChrisRackauckas commented 3 weeks ago

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