> Thank you both for the suggestions. I've re-written the code >> JULIA3, with
> amazing results. JULIA2 was the previous "optimised" code.
>
>
Test 1, J2: 0.751340928, J3: 0.008927998, C: 0.007420171; max-error = 0.0
Test 2, J2: 0.71344446, J3: 0.009042345, C: 0.007583811; max-error = 0.0
Test 3, J2: 0.719000046, J3: 0.008970343, C: 0.007626061; max-error = 0.0
Test 4, J2: 0.707967183, J3: 0.008979525, C: 0.007572056; max-error = 0.0
Test 5, J2: 0.730254977, J3: 0.009012892, C: 0.007649305; max-error = 0.0
. . .. (repeating the test gives consistent results; C is gcc with -O3)
The new code and the test-code are copied below. Of course this means I
have to write C-style codes in Julia to get this sort of performance. Why
does Julia not optimise
dE[:, k] += dJ
dE[:, n] -= dJ
to
for i = 1:d
dE[i, k] += dJ * r[i]
dE[i, n] -= dJ * r[i]
end
?
> I think you should also replace your (s*s*s*s*s) with s^5 - it'll
> automatically do the "right thing", and I'd be surprised if that is slower.
If I revert to s^5, etc, then I lose 2 orders of magnitude.
I will look into NumericExtensions and the profiler next.
Thank you again for the help.
Christoph
function energy_julia3(x)
N = size(x,2); d = size(x,1)
E = 0.0; dE = zeros(d, N)
r = zeros(d);
dJ = 0.; s = 0.
for n = 1:(N-1)
for k = (n+1):N
s = 0.
for i = 1:d
r[i] = x[i,k]-x[i,n]
s += r[i]*r[i]
end
E += 1./(s*s*s*s*s*s) - 2. / (s*s*s)
dJ = -12. * (1./(s*s*s*s*s*s*s) - 1./(s*s*s*s))
for i = 1:d
dE[i, k] += dJ * r[i]
dE[i, n] -= dJ * r[i]
end
end
end
return E, dE
end
function meshgrid{T}(vx::AbstractVector{T}, vy::AbstractVector{T})
m, n = length(vy), length(vx)
vx = reshape(vx, 1, n)
vy = reshape(vy, m, 1)
(repmat(vx, m, 1), repmat(vy, 1, n))
end
function lj_test_juliaopt(N)
x = linspace(0, N, N+1)
x, y = meshgrid(x, x)
x = [x[:] y[:]]'
for n = 1:10
tic(); Ej2, dEj2 = energy_julia2(x); t2 = toq();
tic(); Ej3, dEj3 = energy_julia3(x); t3 = toq();
tic();
dEc = zeros(size(x))
Ec = ccall( (:energy, "./libljtest_c"), Cdouble, (Ptr{Cdouble},
Ptr{Cdouble}, Cint, Cint), x, dEc, size(x,2), size(x,1))
tc = toq();
error = max( abs(Ej2-Ej3), abs(Ej2-Ec), norm(dEj2[:]-dEj3[:], Inf),
norm(dEj2[:]-dEc[:], Inf) )
println("Test ", n, ", J2: ", t2, ", J3: ", t3, ", C: ", tc, ";
max-error = ", error)
end
end
