производительность — Наивный перевод C ++ на Юлию, есть ли место для оптимизации?

Всякий раз, когда я сравниваю производительность Джулии с чудовищным C ++, я всегда поражаюсь. По своему опыту, я часто нахожу Джулию наравне с С ++ или даже быстрее. Мой последний эксперимент сравнивал скорость Джулии с Библиотека C ++ Armadillo используя свои микро тесты. В большинстве случаев я обнаружил, что Джулия быстрее, чем С ++, в таком числовом коде (осторожно: это мои собственные заявления, конечно же, YMMV).

Теперь я сделал наивный перевод этот код C ++ пересечения лучей и треугольников в Джулию и ожидал получить такие же результаты, как мои предыдущие эксперименты. Запустив оба кода, я получил следующий результат для C ++ (скомпилированный с помощью: gcc -std = c ++ 11 -lm -O3 -ffast-math):

Total intersection tests:  100,000,000
Hits:                      4,930,610 ( 4.93%)
Misses:                   95,069,390 (95.07%)

Total time:                       1.51 seconds
Millions of tests per second:    66.18

а для Юлии (версия 0.6.0 на Windows 10) я получил:

Total intersection tests:  100000000
Hits:             5221395  ( 5.221395%)
Misses:           94778605 (94.778605%)

Total time:         2.144723297 seconds
Millions of tests per second:   46.6260613384851

Ниже приведен построчный перевод на Юлию без каких-либо попыток оптимизации. Что не так с моим кодом? или Юлия struct менее эффективны? Где это 1.4X родом из? и как я могу оптимизировать этот код? Я знаю, что использую двойную точность, но это не проблема для Юлии.

struct vec3
x::Float64
y::Float64
z::Float64
end

mutable struct ray
orig::vec3
dir ::vec3
end

sub3(a::vec3, b::vec3) = vec3(a.x - b.x, a.y - b.y, a.z - b.z)
dot3(a::vec3, b::vec3) = a.x*b.x +  a.y*b.y + a.z*b.z
len(v::vec3) = sqrt( v.x * v.x + v.y * v.y + v.z * v.z  )
normalize(v) = begin l = length(v); vec3( v[1]/l, v[2]/l, v[3]/l ) end
cross3(a::vec3, b::vec3) = vec3(a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x)

# Ray-triangle intersection routine
function rayTriangleIntersect(r::ray, v0::vec3, v1::vec3, v2::vec3)::Float64
v0v1 = sub3(v1, v0);
v0v2 = sub3(v2, v0);
pvec = cross3(r.dir, v0v2);
detc = dot3(v0v1, pvec);

(detc < 0.000001) && return -Inf

invDet = 1.0 / detc;
tvec = sub3(r.orig, v0);
u = dot3(tvec, pvec) * invDet;

(u < 0 || u > 1) && return -Inf

qvec = cross3(tvec, v0v1);
v = dot3(r.dir, qvec) * invDet;

(v < 0 || u + v > 1) && return -Inf

return dot3( v0v2, qvec ) * invDet;
end

function randomSphere()::vec3
r1  = rand();
r2  = rand();
lat = acos(2r1 - 1) - pi/2;
lon = 2pi * r2;
return vec3( cos(lat)cos(lon), cos(lat)sin(lon), sin(lat) )
end

function generateRandomTriangles(numTriangles::Int64)#::Array{vec3}(3numTriangles)
vertices = Array{vec3}(3numTriangles)
for i = 1:3numTriangles
vertices[i] = randomSphere();
end
return vertices
end

function main_ray()
const NUM_RAYS = 1000;
const NUM_TRIANGLES = 100 * 1000;

vertices = generateRandomTriangles(NUM_TRIANGLES);
const numVertices = NUM_TRIANGLES * 3;

numHit  = 0;
numMiss = 0;

o  = vec3(0.0, 0.0, 0.0)
d  = vec3(0.0, 0.0, 0.0)
p1 = vec3(0.0, 0.0, 0.0)
r  = ray(o,d)

tTotal = @elapsed for i = 1:NUM_RAYS
r.orig = randomSphere();
p1     = randomSphere();
r.dir  = sub3(p1,r.orig);

for j = 0:div(numVertices,3) - 1
t = rayTriangleIntersect(r, vertices[j*3 + 1],
vertices[j*3 + 2],
vertices[j*3 + 3]);
t >= 0? numHit+=1 : numMiss+=1
end
end

numTests = NUM_RAYS * NUM_TRIANGLES;
hitPerc  = numHit / numTests * 100.0
missPerc = numMiss / numTests * 100.0
mTestsPerSecond = numTests / tTotal / 10^6

println("Total intersection tests:  $numTests");
println("  Hits:\t\t\t    $numHit  ( $hitPerc%)");
println("  Misses:\t\t    $numMiss ($missPerc%)\n");
println("Total time:\t\t\t$tTotal seconds");
println("Millions of tests per second:\t$mTestsPerSecond\n");
end

main_ray();