考虑 Julia 中的以下 4 个函数:它们都选择/计算矩阵 A
的随机列,并将常数乘以该列添加到向量 z
。
slow1
和 fast1
之间的区别在于 z
的更新方式,对于 slow2
和 也是如此fast2
.
1
函数和 2
函数之间的区别在于矩阵 A
是传递给函数还是动态计算。
奇怪的是,对于 1
函数,fast1
更快(正如我在使用 BLAS 而不是 +=
时所期望的那样) ,但对于 2
函数,slow1
更快。
在这台计算机上,我得到以下计时(每个函数的第二次运行):
@time slow1(A, z, 10000);
0.172560 seconds (110.01 k allocations: 940.102 MB, 12.98% gc time)
@time fast1(A, z, 10000);
0.142748 seconds (50.07 k allocations: 313.577 MB, 4.56% gc time)
@time slow2(complex(float(x)), complex(float(y)), z, 10000);
2.265950 seconds (120.01 k allocations: 1.529 GB, 1.20% gc time)
@time fast2(complex(float(x)), complex(float(y)), z, 10000);
4.351953 seconds (60.01 k allocations: 939.410 MB, 0.43% gc time)
对此行为有解释吗?还有一种方法可以让 BLAS 比 +=
更快?
M = 2^10
x = [-M:M-1;]
N = 2^9
y = [-N:N-1;]
A = cis( -2*pi*x*y' )
z = rand(2*M) + rand(2*M)*im
function slow1(A::Matrix{Complex{Float64}}, z::Vector{Complex{Float64}}, maxiter::Int)
S = [1:size(A,2);]
for iter = 1:maxiter
idx = rand(S)
col = A[:,idx]
a = rand()
z += a*col
end
end
function fast1(A::Matrix{Complex{Float64}}, z::Vector{Complex{Float64}}, maxiter::Int)
S = [1:size(A,2);]
for iter = 1:maxiter
idx = rand(S)
col = A[:,idx]
a = rand()
BLAS.axpy!(a, col, z)
end
end
function slow2(x::Vector{Complex{Float64}}, y::Vector{Complex{Float64}}, z::Vector{Complex{Float64}}, maxiter::Int)
S = [1:length(y);]
for iter = 1:maxiter
idx = rand(S)
col = cis( -2*pi*x*y[idx] )
a = rand()
z += a*col
end
end
function fast2(x::Vector{Complex{Float64}}, y::Vector{Complex{Float64}}, z::Vector{Complex{Float64}}, maxiter::Int)
S = [1:length(y);]
for iter = 1:maxiter
idx = rand(S)
col = cis( -2*pi*x*y[idx] )
a = rand()
BLAS.axpy!(a, col, z)
end
end
更新:
分析slow2
:
2260 task.jl; anonymous; line: 92
2260 REPL.jl; eval_user_input; line: 63
2260 profile.jl; anonymous; line: 16
2175 /tmp/axpy.jl; slow2; line: 37
10 arraymath.jl; .*; line: 118
33 arraymath.jl; .*; line: 120
5 arraymath.jl; .*; line: 125
46 arraymath.jl; .*; line: 127
3 complex.jl; cis; line: 286
3 complex.jl; cis; line: 287
2066 operators.jl; cis; line: 374
72 complex.jl; cis; line: 286
1914 complex.jl; cis; line: 287
1 /tmp/axpy.jl; slow2; line: 38
84 /tmp/axpy.jl; slow2; line: 39
5 arraymath.jl; +; line: 96
39 arraymath.jl; +; line: 98
6 arraymath.jl; .*; line: 118
34 arraymath.jl; .*; line: 120
分析fast2
:
4288 task.jl; anonymous; line: 92
4288 REPL.jl; eval_user_input; line: 63
4288 profile.jl; anonymous; line: 16
1 /tmp/axpy.jl; fast2; line: 47
1 random.jl; rand; line: 214
3537 /tmp/axpy.jl; fast2; line: 48
26 arraymath.jl; .*; line: 118
44 arraymath.jl; .*; line: 120
1 arraymath.jl; .*; line: 122
4 arraymath.jl; .*; line: 125
53 arraymath.jl; .*; line: 127
7 complex.jl; cis; line: 286
3399 operators.jl; cis; line: 374
116 complex.jl; cis; line: 286
3108 complex.jl; cis; line: 287
2 /tmp/axpy.jl; fast2; line: 49
748 /tmp/axpy.jl; fast2; line: 50
748 linalg/blas.jl; axpy!; line: 231
奇怪的是,尽管到目前为止函数是相同的,但 col
的计算时间却有所不同。
但 +=
仍然比 axpy!
更快。
最佳答案
现在 Julia 0.6 已经发布了,还有更多信息。要将向量乘以标量,至少有四种选择。根据 Tim 的建议,我使用了 BenchmarkTool 的 @btime 宏。事实证明,循环融合(最朱利安的编写方式)与调用 BLAS 相当。这是 Julia 开发人员值得自豪的事情!
using BenchmarkTools
function bmark(N)
a = zeros(N);
@btime $a *= -1.;
@btime $a .*= -1.;
@btime LinAlg.BLAS.scal!($N, -1.0, $a, 1);
@btime scale!($a, -1.);
end
以及 10^5 个数字的结果。
julia> bmark(10^5);
78.195 μs (2 allocations: 781.33 KiB)
35.102 μs (0 allocations: 0 bytes)
34.659 μs (0 allocations: 0 bytes)
34.664 μs (0 allocations: 0 bytes)
分析回溯显示,scale!
只是在后台调用 blas
,因此它们应该给出相同的最佳时间。
关于julia - BLAS.axpy!比 Julia 中的 += 慢,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/32610062/