Skip to content

update precs when cache.isfresh #533

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.

Sign up for GitHub

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

Closed
wants to merge 1 commit into from
Closed

update precs when cache.isfresh #533

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions ext/LinearSolvePardisoExt.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,11 @@ function SciMLBase.solve!(cache::LinearSolve.LinearCache, alg::PardisoJL; kwargs
@unpack A, b, u = cache
A = convert(AbstractMatrix, A)
if cache.isfresh
if hasproperty(alg, :precs) && !isnothing(alg.precs)
Pl, Pr = cache.alg.precs(x, cache.p)
cache.Pl = Pl
cache.Pr = Pr
end
Comment on lines +135 to +139
Copy link
Contributor

@j-fu j-fu Aug 22, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this needed here ? Pl, Pr are not used by Pardiso.

Copy link
Member Author

@oscardssmith oscardssmith Aug 22, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

oh, is krylov the only one that uses them? I thought there was another...

Copy link
Member

@ChrisRackauckas ChrisRackauckas Aug 22, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

IterativeSolvers, KrylovKit.

phase = alg.cache_analysis ? Pardiso.NUM_FACT : Pardiso.ANALYSIS_NUM_FACT
Pardiso.set_phase!(cache.cacheval, phase)
Pardiso.pardiso(cache.cacheval, A, eltype(A)[])
Expand Down
30 changes: 2 additions & 28 deletions src/common.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,19 +84,8 @@ mutable struct LinearCache{TA, Tb, Tu, Tp, Talg, Tc, Tl, Tr, Ttol, issq, S}
end

function Base.setproperty!(cache::LinearCache, name::Symbol, x)
if name === :A
if hasproperty(cache.alg, :precs) && !isnothing(cache.alg.precs)
Pl, Pr = cache.alg.precs(x, cache.p)
setfield!(cache, :Pl, Pl)
setfield!(cache, :Pr, Pr)
end
if name === :A || name === :p
setfield!(cache, :isfresh, true)
elseif name === :p
if hasproperty(cache.alg, :precs) && !isnothing(cache.alg.precs)
Pl, Pr = cache.alg.precs(cache.A, x)
setfield!(cache, :Pl, Pl)
setfield!(cache, :Pr, Pr)
end
elseif name === :b
# In case there is something that needs to be done when b is updated
update_cacheval!(cache, :b, x)
Expand Down Expand Up @@ -224,20 +213,7 @@ function SciMLBase.reinit!(cache::LinearCache;
u = cache.u,
p = nothing,
reinit_cache = false,)
(; alg, cacheval, abstol, reltol, maxiters, verbose, assumptions, sensealg) = cache

precs = (hasproperty(alg, :precs) && !isnothing(alg.precs)) ? alg.precs : DEFAULT_PRECS
Pl, Pr = if isnothing(A) || isnothing(p)
if isnothing(A)
A = cache.A
end
if isnothing(p)
p = cache.p
end
precs(A, p)
else
(cache.Pl, cache.Pr)
end
(; alg, cacheval, abstol, reltol, maxiters, verbose, assumptions, sensealg, Pl, Pr) = cache
isfresh = true

if reinit_cache
Expand All @@ -250,8 +226,6 @@ function SciMLBase.reinit!(cache::LinearCache;
cache.b = b
cache.u = u
cache.p = p
cache.Pl = Pl
cache.Pr = Pr
cache.isfresh = true
end
end
Expand Down
5 changes: 5 additions & 0 deletions src/iterative_wrappers.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -226,6 +226,11 @@ end

function SciMLBase.solve!(cache::LinearCache, alg::KrylovJL; kwargs...)
if cache.isfresh
if hasproperty(alg, :precs) && !isnothing(alg.precs)
Pl, Pr = cache.alg.precs(x, cache.p)
Copy link
Contributor

@j-fu j-fu Aug 22, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here we may need

Pl, Pr = cache.alg.precs(cache.A, cache.p)

cache.Pl = Pl
cache.Pr = Pr
end
Copy link
Contributor

@j-fu j-fu Aug 22, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And generally, could we have it like this ?

function SciMLBase.solve!(cache::LinearCache, alg::KrylovJL; reuse_precs=false, kwargs...)
    if cache.isfresh
        if hasproperty(alg, :precs) && !isnothing(alg.precs)
            if !reuse_precs || (cache.Pl == nothing && cache.Pr == nothing)
                Pl, Pr = cache.alg.precs(cache.A, cache.p)
                cache.Pl = Pl
                cache.Pr = Pr
            end
        end
 ...

EDIT: I meant reuse_precs=false as default.

Copy link
Contributor

@j-fu j-fu Aug 22, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also, in the moment this leads to precs being called twice: once in init and once here.

Copy link
Member

@ChrisRackauckas ChrisRackauckas Aug 22, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

init should set fresh to false.

Copy link
Member Author

@oscardssmith oscardssmith Aug 22, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ah that's tricky. I think I will need to add an extra flag to fix this.

Copy link
Member Author

@oscardssmith oscardssmith Aug 22, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChrisRackauckas really? I think the code relies on it being true to set up some of the C libraries.

Copy link
Member Author

@oscardssmith oscardssmith Aug 24, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the problem with the current way is that it makes the fallback slow. if we didn't have the "fast path", we could set cache.isfresh=false which would recover the speed of the fast path for all types, as long as the user inits on a matrix they actually want to solve.

Copy link
Member

@ChrisRackauckas ChrisRackauckas Aug 24, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You could just set it to false on the slow path.

Copy link
Member Author

@oscardssmith oscardssmith Aug 24, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the current approach is better by 1 factorization for users that create a LinearCache with garbage data, but those users will almost always be doing lots of factorizations, minimizing the advantage. OTOH, the current approach has overhead for users who just call solve (since the fake factorization isn't free), or users of the cache form who init on data they want to solve eventually.

Copy link
Member Author

@oscardssmith oscardssmith Aug 24, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

we could set it to false on the slow path, but doing so would be pretty hard since init_cacheval doesn't get the cache.

Copy link
Member Author

@oscardssmith oscardssmith Aug 24, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

and actually users that are initing with fake matrices can init with a 1x1 themselves, which will completely remove the penalty of my suggestion

solver = init_cacheval(alg, cache.A, cache.b, cache.u, cache.Pl, cache.Pr,
cache.maxiters, cache.abstol, cache.reltol, cache.verbose,
cache.assumptions, zeroinit = false)
Expand Down