Reverse mode is going to give columns, and I don’t think you need that. Using double forward mode will be the fastest here. You’ll need to mapchildren to remove the tracker information (or use the Flux#zygote branch) and then just forward diff (or use a hyperdual)
If you do want to Forward-over-Reverse for Hess-vec products though, it is implemented in SparseDiffTools.jl
https://github.com/JuliaDiffEq/SparseDiffTools.jl#jacobian-vector-and-hessian-vector-products
but note that our tests don’t show that using Zygote here is the fastest yet 