小编Nic*_*o B的帖子

我需要执行(复杂)矩阵的离散卷积,并在Julia中定义以下函数:

function convolve(M::Array{Complex{Float64},2}, K::Array{Float64,2}, p::Int)
    (n,m) = size(M)
    res = zeros(Complex{Float64},n)
    for k=1:p
        for l=1:n
            res[l] += M[l,k]*K[l,end-p+k]
        end
    end
    return res
end

我这样使用它:

M=complex(rand(2000,2000))
K=rand(2000,2000)
@time convolve(M,K,2000,0)

现在这是相对快速和令人惊讶的,比矢量化版本更快(约3倍),我用它替换内部循环res += M[:,k].*K[:,end-p+k].(我认为这是由于临时阵列的很多记忆分配,我可以忍受).

但是矢量化的MATLAB代码运行速度提高了大约5倍:

function res = convolve(M, K, p)
    n = size(M,1);
    res = zeros(n,1);
    for k=1:p
        res = res + M(:,k).*K(:,end-p+k);
    end
end

我做错了什么？如何让Julia像MATLAB一样快速地执行元素乘法？这是一个索引问题吗？

注意:我已经检查过@code_warntype没有类型犹豫不决的有趣业务(没有Any或Union等),但问题可能更微妙.宏@code_llvm产生了令人惊讶的长输出,但我不是专家所以我很难看到发生了什么.