I have a function
f! that operates in-place in a
D dimensional Vector, used in DifferentialEquations. For example
@inline @inbounds function f!(du, u, p, t) σ = p; ρ = p; β = p du = σ*(u-u) du = u*(ρ-u) - u du = u*u - β*u return nothing end
However, I am referring to the general case where I have a general
f! that operates in-place on a
Now I want to modify this function, so that it operates on
k such vectors “in parallel”. These Vectors are currently the columns of a matrix
S, but the storage format does not matter and can be changed.
What I Do is
function create_parallel_eom(f!, S) k = size(S) veom! = (du, u, p, t) -> begin for j in 1:k f!(view(du, :, j), view(u, :, j), p, t) end return end return veom! end
and then use this
veom! to create a new
ODEProblem. If I try to do this for 2 parallel vectors, I get a slowdown of 150x instead of the “theoretical” 2x, simply because the vector field function is absurdly expensive compared to calculating it twice without views.
Any tips for making this more performant? At the moment I am looking at