Allow the use of a combined function + Jacobian

[astrozot]

I need to solve a relatively complicated function (which depends on changing parameters) several thousands or millions of times, and I am trying to obtain a code that is as fast as possible. I am able to provide the Jacobian, but I would like to be able to compute the gradient and the function itself within the same function to avoid repeating, costly computations.

However, I do not see how one can use this package (and the sub package SimpleNonlinearSolve) with a setting where f and df share common computations. Ideally, I would like to be able to provide three separate functions: f, df, and fdf, with the last one used to compute at the same time the function and its Jacobian.

See also

https://discourse.julialang.org/t/specifying-analytical-jacobian-in-nonlinearproblem-a-la-optim-only-fg/115991

for a similar discussion.

P.S. For performance reasons, since my problem is in low dimensions I am using non-allocating functions and StaticArrays, so memoization is not an option.

[Nidish96]

Bump.

I would definitely like to see this for my applications!

[Nidish96]

I wanted to add how I'm currently using this in my applications (since I don't usually write different functions for residue and jacobian):

R = zeros(mdl.Ndofs*Nhc);
dRdU = zeros(mdl.Ndofs*Nhc, mdl.Ndofs*Nhc);
dRdw = zeros(mdl.Ndofs*Nhc);
HBRESFUN!(Uw0, mdl, Fl, h, N; R=R, dRdU=dRdU, dRdw=dRdw)

fun = NonlinearFunction((r,u,p)->HBRESFUN!([u;p], mdl, Fl, h, N; R=r),
    jac=(J,u,p)->HBRESFUN!([u;p], mdl, Fl, h, N; dRdU=J),
    paramjac=(Jp,u,p)->HBRESFUN!([u;p], mdl, Fl, h, N; dRdw=Jp));

prob = NonlinearProblem(fun, Uw0[1:end-1], Uw0[end]);
sol = solve(prob, show_trace=Val(true));

Here, the function "HBRESFUN!" does the appropriate computations based on which of R,dRdU,dRdw are not nothing. So I'm basically giving the same as inputs to NonlinearFunction.
It would be really great if an inbuilt way of doing this in a neater fashion is available.