使用Julia更新运算符进行不必要的分配

Question

考虑以下功能:

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

当我运行它时,我得到以下内容:

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

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

为什么要分配这么多内存并花费这么多时间进行垃圾回收？代码应该是类型稳定的,我希望+=运营商不执行重新分配.但是,似乎每次添加两个矩阵时它都会重新分配.

我应该认为这是朱莉娅的一个错误吗？更重要的是,我如何以不重新分配的方式编写此代码？

编辑:修正错字.

Answer 1

@cd98 请求我的三嵌套循环解决方案，该解决方案解决了分配问题，但我认为其性能低于等效的矢量化版本。这里是：

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)

和输出：

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