BLAS.axpy! plus lent que + = dans Julia

Question

Considérons les 4 fonctions suivantes dans Julia: Ils choisissent/calculent tous une colonne aléatoire d'une matrice A et ajoute une fois cette colonne à un vecteur.

La différence entre slow1 et fast1 comment z est mis à jour et même pour slow2 et fast2. La différence entre les fonctions 1 et 2 est de savoir si la matrice A est passée aux fonctions ou calculée à la volée.

La chose étrange est que, pour les 1 fonctions, fast1 est plus rapide (comme je l'attendre lors de l'utilisation BLAS au lieu de +=), mais pour les fonctions 2slow1 est plus rapide. Sur cet ordinateur, je reçois les horaires suivants (pour la deuxième manche de chaque fonction):

@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) 

Y at-il une explication à ce comportement? Et un moyen de rendre BLAS plus rapide que +=?

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 

Mise à jour: profilage 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 

profilage 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 

Bizarrement, le temps de calcul de col diffère même si les fonctions sont identiques à ce point . Mais += est encore relativement plus rapide que axpy!.

Answer 1

Plus d'infos maintenant que julia 0.6 est sorti. Pour multiplier un vecteur par un scalaire en place, il y a au moins quatre options. Suivant les suggestions de Tim, j'ai utilisé la macro @btime. Il s'avère que la fusion de boucles, la façon la plus julienne de l'écrire, est à égalité avec l'appel de BLAS. C'est quelque chose dont les développeurs de julia peuvent être fiers!

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 

Et les résultats pour 10^5 nombres.

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) 

Le profilage backtrace montre que scale! appelle simplement blas en arrière-plan, donc ils devraient donner le même meilleur temps.