Skip to content

[WIP] Simplify derivative_utils.jl #876

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

Merged
merged 2 commits into from
Aug 22, 2019
Merged

[WIP] Simplify derivative_utils.jl #876

Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
dcff0bd
add calcJ! for oop
huanglangwen Aug 16, 2019
58804f1
simplify calc_J! and add jacobian2W for oop
huanglangwen Aug 16, 2019
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
24 changes: 23 additions & 1 deletion src/derivative_utils.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ jacobian update function, then it will be called for the update. Otherwise,
either ForwardDiff or finite difference will be used depending on the
`jac_config` of the cache.
"""
function calc_J!(integrator, cache::OrdinaryDiffEqMutableCache, is_compos)
function calc_J!(integrator, cache::OrdinaryDiffEqCache, is_compos)
if isdefined(cache, :nlsolver)
calc_J!(cache.nlsolver, integrator, cache, is_compos)
elseif isdefined(cache, :J)
Expand All @@ -106,6 +106,17 @@ function calc_J!(nlsolver::NLSolver, integrator, cache::OrdinaryDiffEqMutableCac
is_compos && (integrator.eigen_est = opnorm(J, Inf))
end

function calc_J!(nlsolver::NLSolver, integrator, cache::OrdinaryDiffEqConstantCache, is_compos)
Copy link
Member

@devmotion devmotion Aug 16, 2019

Choose a reason for hiding this comment

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

There already exists an implementation for out-of-place function, it's called calc_J. So if we need an in-place variant calc_J! (which I'm not sure about) we could just define calc_J! as cache.J = calc_J(...)

@unpack t,dt,uprev,u,f,p = integrator
if DiffEqBase.has_jac(f)
nlsolver.cache.J = f.jac(uprev, p, t)
else
nlsolver.cache.J = jacobian(nlsolver.uf,uprev,integrator)
Copy link
Member

@devmotion devmotion Aug 16, 2019

Choose a reason for hiding this comment

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

We should update uf first, similar to the other implementations of calc_J!. This is currently missing in the out-of-place implementations calc_J.

Copy link
Contributor Author

@huanglangwen huanglangwen Aug 16, 2019

Choose a reason for hiding this comment

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

But why does calc_J produce the right answer without updating uf ?

Copy link
Member

@devmotion devmotion Aug 16, 2019

Choose a reason for hiding this comment

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

Since it's updated in calc_W! in lines such as https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/dcff0bd6809adc640ce910ddbeeba8ed05091431/src/derivative_utils.jl#L525 (but not for p, so that's still not correct). I think we should remove those updates in calc_W! since the update should only happen if we actually need uf, which is only the case in calc_J.

end
integrator.destats.njacs += 1
is_compos && (integrator.eigen_est = opnorm(nlsolver.cache.J, Inf))
end

function calc_J_in_cache!(integrator, cache::OrdinaryDiffEqMutableCache, is_compos)
@unpack t,dt,uprev,u,f,p = integrator
J = cache.J
Expand All @@ -121,6 +132,17 @@ function calc_J_in_cache!(integrator, cache::OrdinaryDiffEqMutableCache, is_comp
is_compos && (integrator.eigen_est = opnorm(J, Inf))
end

function calc_J_in_cache!(integrator, cache::OrdinaryDiffEqConstantCache, is_compos)
Copy link
Member

@devmotion devmotion Aug 16, 2019

Choose a reason for hiding this comment

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

As I said, I'm not sure we want to do this.

@unpack t,dt,uprev,u,f,p = integrator
if DiffEqBase.has_jac(f)
cache.J = f.jac(uprev, p, t)
else
cache.J = jacobian(cache.uf,uprev,integrator)
end
integrator.destats.njacs += 1
is_compos && (integrator.eigen_est = opnorm(cache.J, Inf))
end

"""
WOperator(mass_matrix,gamma,J[;transform=false])

Expand Down