Juliaで自作VectorとArray{Float64,1}の速度比較

Juliaの勉強を始めた。
Juliaは数値計算の関数を大量に有しているけど、可変長のArrayをもとにしていて、3x3の行列や3x1ベクトルの単純なものだと自分で定義したほうが高速なのだろうか?という疑問が湧いたので実験してみた。
もしかしたらJuliaは裏でものすごいことをやっていてるかもしれんし。

3x1のベクトルVecを以下のように定義した。

immutable Vec
    x::Float64
    y::Float64
    z::Float64
end
dot(a::Vec, b::Vec) = (a.x*b.x + a.y*b.y + a.z*b.z)
cross(a::Vec, b::Vec) = Vec(a.y*b.z - a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x)

ただ外積やらなんやらを計算する関数を作った。

function cross_loop(N::Int, a::Vec, b::Vec)
    tmp = 0.0
    for i in 1:N
        v = cross(a,b)
        tmp += dot(v,v)
    end 
    tmp/N
end

function cross_loop(N::Int, a::Array{Float64,1}, b::Array{Float64,1})
    tmp = 0.
    for i in 1:N
        v = Base.cross(a,b)
        tmp += Base.dot(v,v) 
    end
    tmp / N
end

Arrayの方はデフォルトで入っているBase.dotやcrossをつかった。

N=10^7くらいにして計算
速度は@timeで計測した。1回目はコンパイルに時間を取られるので2回目の時間を使用。

N = 10^7
a = Vec(1,2,3)
b = Vec(4,5,6)
v = [1.,2.,3.]
u = [4.,5.,6.]

@time cross_loop(N, a, b)
@time cross_loop(N, u, v)

結果はダントツで自作のほうが速い。

elapsed time: 0.013432361 seconds (96 bytes allocated)
elapsed time: 2.079223219 seconds (1280000096 bytes allocated, 43.49% gc time)

ちなみにVecをimmutableでなくtypeで定義すると、それだけでかなり遅くなる。GCが働いているようだ。

elapsed time: 0.419364716 seconds (320000096 bytes allocated, 53.07% gc time)

結論:3次元に特化したベクトル計算では自分で定義してやったほうが速い。