Skip to content

Conversation

@huanglangwen
Copy link
Contributor

Continuing from #873 and #871 . The goal is to simplify four calc_W! functions into two groups: one for rosenbrock, one for nlnewton. calc_W! should behave the same for OOP and IIP.

@ChrisRackauckas
Copy link
Member

IMO one change we should make while doing this is that the W we want to calculate on should be passed in. See #872 for what happens when doing multithreaded implicit methods. That extra method should just be deleted when this is fixed.

@coveralls
Copy link

coveralls commented Aug 16, 2019

Coverage Status

Coverage decreased (-0.04%) to 78.869% when pulling 58804f1 on huanglangwen:simplify_dutils into 893bc16 on JuliaDiffEq:master.

@codecov
Copy link

codecov bot commented Aug 16, 2019

Codecov Report

Merging #876 into master will decrease coverage by 0.03%.
The diff coverage is 21.05%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #876      +/-   ##
==========================================
- Coverage   80.31%   80.27%   -0.04%     
==========================================
  Files          92       92              
  Lines       30303    30322      +19     
==========================================
+ Hits        24337    24341       +4     
- Misses       5966     5981      +15
Impacted Files Coverage Δ
src/derivative_utils.jl 68.58% <21.05%> (-2.9%) ⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 893bc16...58804f1. Read the comment docs.

@devmotion
Copy link
Member

As soon as the algorithm-specific fields of the non-linear solvers are moved to the cache, there's no need anymore for separate calc_J and calc_J! methods (see #868). It's not possible right now, but that's something we should keep in mind.

is_compos && (integrator.eigen_est = opnorm(J, Inf))
end

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

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(...)

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

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

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

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.

is_compos && (integrator.eigen_est = opnorm(J, Inf))
end

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

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.

J = f.jac(uprev, p, t)
else
cache.uf.t = t
cache.uf.p = p
Copy link
Member

Choose a reason for hiding this comment

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

Great! 👍 Can you remove the updates of uf in the oop version of calc_W!?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I'm thinking of delete oop version of cache_W! directly.

Copy link
Member

Choose a reason for hiding this comment

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

And only use update_W!? Currently I don't like that calc_W! is called calc_W! both for in-place and out-of-place functions but updates W only for the in-place functions and returns a new W for the out-of-place functions. So I thought maybe it could be done similar to calc_J and calc_J!, i.e., we could have an out-of-place variant calc_W and an updating variant calc_W!.

Moreover, I'm not happy with the fact that at the moment there's a lot of reuse strategies hidden in the implementation of calc_W!. Maybe it would be cleaner if calc_W! would just calculate W, without considering any reuse strategies, similar to what we do in calc_J. Ideally then we would have a separate function (such as update_W!) that first figures out if J should be updated (and then calls calc_J!) and then if W should be updated (and then calls calc_W!).

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Agree.

@huanglangwen
Copy link
Contributor Author

I'm a bit confused by the calc_W!. Rosenbrock and extrapolation methods call calc_W!(integrator, cache, dtgamma, repeat_step, W_transform) while methods with nlnewton call calc_W!(nlsolver, integrator, cache, dt, repeat_step, true). So why the first calc_W! has so many logic for newton methods like W_dt = isnewton ? cache.nlsolver.cache.W_dt : dt ?

@devmotion
Copy link
Member

It was one implementation in the beginning and then just copied to a new function when the non-linear solvers were changed. I'm sure many parts can be improved.

@huanglangwen
Copy link
Contributor Author

huanglangwen commented Aug 16, 2019

Why is f.Wfact called every time (if exists) when calc_W! is called?

@kanav99
Copy link
Contributor

kanav99 commented Aug 16, 2019

Check out #878 if that helps us out for this patch

@huanglangwen
Copy link
Contributor Author

I guess I'll wait for #878 . @kanav99 has done most of the work and we just need to append some oop functions to it.

@ChrisRackauckas
Copy link
Member

and this intersects with SciML/DiffEqBase.jl#325 . @devmotion might want to take note of all of these and make sure that works with these 3 planned changes.

@devmotion
Copy link
Member

The non-linear solver changes will break the accesses nlsolver.uf but there's a bunch of things that we have to fix after the DiffEqBase update, so I would go on and merge this PR in it's current state. Tests pass and oop functions can be added in a separate PR, IMO.

@devmotion
Copy link
Member

OK, let's merge this and continue improving derivative_utils.jl in separate PRs.

@devmotion devmotion merged commit e46c08b into SciML:master Aug 22, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants