是否可以加速去矢量化的 Julia 代码答案

【问题标题】：Is it possible to speed up devectorized Julia codes是否可以加速去矢量化的 Julia 代码
【发布时间】：2015-08-12 15:09:00
【问题描述】：

作为我脚本的一部分。我有一些代码如下（devectorized julia - as possible）

for kk=1:n # Main loop
    for j=1:m
        rhs[j]=2*u0[j]-alf*dt*u1[j]-2*mu*u2[j];    
    end
    c=lhs\rhs'; #c: coefficients to be obtained

    u2=c'*h;
    u1=c'*p.-c'*f;
    u0=c'*Q-c'*f*x;

    for j=1:m
        for i=1:m
            lhs[j,i]=2*(Q[i,j]-x[j]*f[i])+alf*dt*(p[i,j]-f[i])+eps*dt*(Q[i,j]-x[j]*f[i])*u1[j]+eps*u0[j]*dt*(p[i,j]-f[i])-2*mu*h[i,j];
        end
    end

end

其中 h、p、Q、lhs 是mxm 矩阵； u0、u1、u2、rhs 和 x 是 1xm 数组，alf、dt、mu、eps 是标量常数，f、c 是 mx1 数组。我在脚本开始时预先分配了矩阵和数组。上述代码的向量化形式如下

for kk=1:n # Main loop

    rhs=2*u0-alf*dt*u1-2*mu*u2;    

    c=lhs\rhs'; #c coefficients to be obtained

    u2=c'*h;
    u1=c'*p.-c'*f;
    u0=c'*Q-c'*f*x;

    lhs=2*(Q-f*x)+alf*dt*(p.-f)+eps*dt*(Q-f*x).*u1+eps*dt*u0.*(p.-f)-2*mu*h;
    lhs=lhs';

end

例如，对于 n=100 和 m=64，经过的时间如下：

去矢量化的 Julia：1.8 秒

矢量化 Julia：0.2 秒

矢量化 numpy：0.04 秒

矢量化 julia 代码比去矢量化 julia 代码快大约 9 倍，矢量化 python 代码比矢量化 julia 代码快大约 5 倍。

对于 n=500 和 m=256

去矢量化的 Julia：85.589233013 秒

矢量化 Julia：8.232898003 秒

矢量化 numpy：1.62000012398 秒

我的问题：在这种情况下是否可以提高 Julia 的性能？

【问题讨论】：

您应该查看文档的performance tips 部分。具体来说，查看accessing arrays according to memory order 的部分并在函数内部而不是在顶层运行所有代码。
你似乎在全局范围内运行，这就是为什么矢量化版本更快，因为它们更快地脱离全局范围
@spencerlyon2：将所有代码放入函数中可以很好地提高性能。谢谢。

标签： performance julia

【解决方案1】：

我认为也可以像这样对u0,u1,u2 的计算进行去矢量化：

function vectorized()
    m = [1.0 2.0 3.0; 1.0 2.0 3.0; 1.0 2.0 3.0]
    c = [1.0, 2.0, 3.0]

    for i in 1:100000
        x1 = c'*m
        x2 = c'*m
        x3 = c'*m
    end

    return
end

function vectime(N)
    timings = Array(Float64, N)

    # Force compilation
    vectorized()

    for itr in 1:N
        timings[itr] = @elapsed vectorized()
    end

    return timings
end

println("vectorized=",mean(vectime(20)))

function devectorized()
    m = [1.0 2.0 3.0; 1.0 2.0 3.0; 1.0 2.0 3.0]
    c = [1.0, 2.0, 3.0]
    x1 = [0.0, 0.0, 0.0]
    x2 = [0.0, 0.0, 0.0]
    x3 = [0.0, 0.0, 0.0]
    mx = 3
    for i in 1:100000
        for k in 1:mx
            for kk in 1:mx
              x1[k]=x1[k]+c[k]*m[k,kk];
              x2[k]=x2[k]+c[k]*m[k,kk];
              x3[k]=x3[k]+c[k]*m[k,kk];
            end
        end
    end
    return
end

function dvectime(N)
    timings = Array(Float64, N)

    # Force compilation
    devectorized()

    for itr in 1:N
        timings[itr] = @elapsed devectorized()
    end

    return timings
end

println("devectorized=",mean(dvectime(20)))

以上代码结果：

vectorized=0.17680755404999998
devectorized=0.00441064295

【讨论】：