Hello,

I have a class of `ODEFunction`

s for which I know how to exploit a lot of the internal structure for efficient calculations of Jacobians. I am a bit confused about how I may teach the `SciML`

-universe about those. Generally, I think I have the following three use-cases:

- Linear solve within stiff solvers (nonsingular mass matrix)
- parameter estimation with SciMLSensitivity
- initial state estimation with SciMLSensitivity

According to the `ODEFunction`

docstring there are 4 relevant parameters

`jac`

(and`jac_prototype`

) for the matrix jacobian`paramjac`

(but no`paramjac_prototype`

?) matrix jacobian with respect to parameters`jvp`

and`vjp`

The `SciMLSensitivity`

docs suggest to supply `vjp`

and then use `ZygoteVJP`

to make sure the solver uses those. However this part of the docs seem out of date as signature for the custom `vjp`

does not match the one from the `ODEFunction`

docstring, so I am not sure if this is up to date.

Also, shouldn’t there be something like `paramvjp`

for the parameter estimation case? Will the `vjp`

definition also be used inside the stiff solvers or does this only work for SciMLSensitivity and the linear solve will look for the `jac`

? Should I define all of those options above?

An alternative option might be create a “lazy” `SciMLOperator`

for my `jac`

and `paramjac`

without actually building the full jacobian.

Another completely different option that comes to mind is to define custom ChainRules rules for my function. But will those be actually used inside the package?

It seems like everything I need is there, but I am a bit overwhelmed by the different options. So please help me understand how those different approaches relate to each other. Thanks!