Julia でループ展開しても、あまり速くない？

疎行列と密ベクトルの積を計算するような以下のコードを書いてみた。

==== Ax.jl =====

function Ax(n, density, seed)
    srand(seed)
    println("Generating A")
    time_A = @elapsed @time A = sprand(n,n,density)
    x = rand(n,1)
    println("Computing b = A*x")
    time_b = @elapsed @time b = A*x

    println("Computing b2 = A*x")
    time_b2 = @elapsed @time begin
        b2 = zeros(n,1)

        for j=1:n
            const xj = x[j]
            for index_i = A.colptr[j]:(A.colptr[j+1]-1)
                b2[A.rowval[index_i]] += A.nzval[index_i]*xj
            end
        end
    end

    println("norm(b-b2) = ", norm(b-b2))

    println("time (A) = ", time_A, " time (b) = ", time_b, " time(b2) = ", time_b2)

    return A, x, b, b2

end

==============

実行結果は以下の通り。

julia> include("Ax.jl"); (A, x, b,b2) = Ax(2000,0.8,1024); (A,x,b,b2) = Ax(2000,0.8,1024);
Generating A
elapsed time: 0.173287485 seconds (76853144 bytes allocated, 24.33% gc time)
Computing b = A*x
elapsed time: 0.011240257 seconds (16064 bytes allocated)
Computing b2 = A*x
elapsed time: 2.677362552 seconds (498951072 bytes allocated, 9.45% gc time)
norm(b-b2) = 0.0
time (A) = 0.173453157 time (b) = 0.011356873 time(b2) = 2.677506154
Generating A
elapsed time: 0.141352727 seconds (76853144 bytes allocated, 12.73% gc time)
Computing b = A*x
elapsed time: 0.010914271 seconds (16064 bytes allocated)
Computing b2 = A*x
elapsed time: 2.708701738 seconds (498951072 bytes allocated, 9.03% gc time)
norm(b-b2) = 0.0
time (A) = 0.141499582 time (b) = 0.011023018 time(b2) = 2.708858507

b の計算はもともとあるもので、b2 のほうがループを展開して計算している。

ちなみに、２回実行しているのは、Julia のコンパイルで時間がかかると遅くなるかもしれない、ということで２回実行してみた。

うーむ、やはり numpy のように、できるだけ行列やベクトルの演算に持ち込んだ方が

速いのだろうか？

そうそう、@elapsed, @time マクロは begin - end で挟むと、その間の時間を測定してくれる。

@elapsed マクロのほうは返り値で経過時間を返してくれるのでそれを代入しているが、@time マクロのほうが細かい情報が表示される。