I have been switching over a 2D flux computation method to use CartesianIndices so it will work for arrays of any dimension, ala https://julialang.org/blog/2016/02/iteration/
Unfortunately, I’ve been running into a bunch of memory allocation “gotchas” and I don’t understand what I need to do to avoid them. Here’s an example
function test(a) σ = 0. R = CartesianIndex(size(a)[1:end-1]); o = oneunit(R) @time for I ∈ 2o:R σ += a[I,1]-a[I-o,2]^2 # actually, a more complex function... end return σ end
But when I run this I get slow code with a giant memory issue
julia> test(rand(256,256,2)); 0.046946 seconds (585.23 k allocations: 13.891 MiB)
Please, could someone let me know:
- What is going on here?
- How am I supposed to do this elegantly?