Filling a preallocated array without allocations

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)

I got

julia> @benchmark loop($kfaux)
BenchmarkTools.Trial: 
  memory estimate:  6.67 KiB
  allocs estimate:  61
  --------------
  minimum time:     3.361 μs (0.00% GC)
  median time:      3.467 μs (0.00% GC)
  mean time:        4.745 μs (15.47% GC)
  maximum time:     860.028 μs (98.64% GC)
  --------------
  samples:          10000
  evals/sample:     8

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)
4 Likes

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)
5 Likes

Can you explain why this is allocation less? Would like to learn :slight_smile:

Kind regards

Sure! :slightly_smiling_face:

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.

2 Likes

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.

Kind regards

“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.

6 Likes

Thanks! It helped a lot! :slight_smile: