Hello, does anybody know if it is possible to fill a preallocated array without allocations?
In other words, is there a way of performing this function loop without allocations
mutable struct KF
F::Array{Float64, 3}
function KF()
new(rand(2, 2, 100))
end
end
function loop(kfaux::KF)
@inbounds for t in 40:100
kfaux.F[:, :, t] = kfaux.F[:, :, t-1]
end
return
end
kfaux = KF()
using BenchmarkTools
@benchmark loop($kfaux)
Sure–one way to do it is just to write out the loop:
julia> function loop(kfaux::KF)
for t in 40:100
for j in axes(kfaux.F, 2)
for i in axes(kfaux.F, 1)
kfaux.F[i, j, j] = kfaux.F[i, j, t - 1]
end
end
end
end
loop (generic function with 1 method)
julia> @btime loop($kfaux)
496.222 ns (0 allocations: 0 bytes)
In addition, if the dimensions of your example is representative of what you are really trying to do, maybe it makes sense to swap Array{T, 3} to Vector{SMatrix} using StaticArrays
using StaticArrays
mutable struct KFS
F::Vector{SMatrix{2,2,Float64,4}}
function KFS()
new([@SMatrix rand(2,2) for _ in 1:100])
end
end
function loop(kfaux::KFS)
F = kfaux.F
@inbounds for t in 40:100
F[t] = F[t-1]
end
return
end
Benchmarks:
julia> using BenchmarkTools
julia> @btime loop($kfaux) # Original version
2.503 ÎĽs (61 allocations: 6.67 KiB)
julia> @btime loop′($kfaux) # Robin's version
459.447 ns (0 allocations: 0 bytes)
julia> @btime loop($kfauxs)
18.602 ns (0 allocations: 0 bytes)
When you write x[:, :, t], you are making a copy of that slice of the array, so this line:
x[:, :, t] = x[:, :, t - 1]
makes a copy of x[:, :, t - 1] and then assigns that to x[:, :, t].
You can use a view instead:
x[:, :, t] .= @view x[:, :, t - 1]
but it does not make a huge difference in your particular case because the sizes of the first two dimensions of your array are so small.
A loop doesn’t require creating a copy or a view at all, so it’s a nice solution. I would also second @fredrikekre’s suggestion: if your data structure is actually a collection of N 2x2 matrices, then it’s much nicer to work with it in that way, and using StaticArrays makes that very efficient.
Thanks, so loops are actually pretty smart, nice to know. Used to that vectorizing would be faster, so that is why I sometimes have a hard time grasping these kind of things, but I really liked your explanation.
“Vectorizing” in that sense is only necessary in languages where loops are slow (matlab, python, etc.). After all, a vectorized operation in numpy or matlab is actually implemented via a loop in a fast language (like C). Julia is a fast language, so you can choose a vectorized operation or a loop based on whichever one makes your particular task easier to write, read, or understand.