Unnecessary allocations using Julia update operators

Question

Consider the following function:

function mytest(x, b)
    y = zeros(x[:,:,1])
    for i in 1:length(b)
        y += b[i] * x[:,:,i]
    end
    return y
end

When I run it, I get the following:

x = rand(30,30,100000)
b = rand(100000)
@time mytest(x,b)

elapsed time: 0.571765222 seconds (727837732 bytes allocated, 66.49% gc time)

Why is it allocating so much memory and spending so much time doing garbage collection? The code should be type stable, and I would expect the += operator to perform no re-allocation. However, it seems that it is re-allocating each time it adds two matrices.

Should I consider this to be a bug in Julia? And more important, how can I write this code in a way that does not re-allocate?

EDIT: fixed typo.

Answer 1

@cd98 requested my three-nested-loop solution, which solves the allocation problem but I assume would underperform an equivalent vectorized version. Here it is:

function mytest(x, b)
    d1, d2, d3 = size(x)
    y = zeros(eltype(x), d1, d2)
    for i in 1:d3
        for j in 1:d2
            for k in 1:d1
                y[k,j] += b[i] * x[k,j,i]
            end
        end
    end
    return y
end

x = rand(30,30,100000)
b = rand(100000)
@time mytest(x,b)
@time mytest(x,b)

And the output:

elapsed time: 0.218220119 seconds (767172 bytes allocated)
elapsed time: 0.197181799 seconds (7400 bytes allocated)