Some clarifications:
Should we want to parallelize this rootfinding business with the same approach as for the ODEs
When you write “parallelize”, you mean create a parallelized implementation of the root-solver, comparable to the one we’re trying to build with the ODE integrator? I’m certainly interested! This is why we (Bill and I) considered using the KINSOL solver instead of Eigen’s. That said, the design should be agnostic to which solver we use.
we would need to make the default implementation use AD to get the Jacobian
I’m not sure what you mean. For the root-solver, we are not solely using AD. If we did that, we’d have to template the entire dogleg function (which is what I did for the current matrix_exp() function). We do use jacobian() to get Jx and Jy, where I’m using Michael’s notation. Which Jacobian would the user be passing? ∂h(y)/∂y ?