Le dimanche 25 mai 2014 à 07:39 -0700, Christoph Ortner a écrit : > > > > I just started experimenting with Julia, as a platform for rapid > prototyping for some more exotic molecular modelling applications > where no black-box software exists. So far I like very much what I > find. To test its performance, I implemented a simple Lennard-Jones > energy and force assembly 5 times: > > 1. Julia1 : a "pretty" code, how I would like to write i > 2. Julia2 : an ugly code with some performance optimisation > 3. C : an even uglies C code > 4. Matlab1 : pretty matlab code > 5. Matlab2 : ugly matlab code > > > Typical runtimes for 30 particles (s): > Julia1: 1.37 > Julia2: 0.74 > C: 0.0267 > Matlab1: 12.51 > Matlab2: 1.91 > > > Good news 1: Julia 1 kill Matlab1 and clearly beats Matlab2 as well. > Good news 2: Julia2 is only about twice as fast as Julia1 > Bad news: I cannot get better than a factor of 3 of C, which is not > bad, but not the 70-99% that the Julia webpage promises :) > (and this C code is even called from Julia, otherwise I would get > another little bit of speed-up) > > > Any advise on further performance improvements would be great > appreciated? (JULIA2 is copied below, I can provide the other codes > as well if helpful) > (grouping the multiplications more cleverly makes virtually no > difference) > > > Many thanks, > Christoph > > > > > JULIA2 > function energy_inline(x) > E = 0.0 > dE = zeros(size(x)) > N = size(x,2) > d = size(x,1) > for n = 1:(N-1) > for k = (n+1):N > r = x[:,k] - x[:,n]
Currently, AFAIK, this command is going to make a copy of the two slices, and a allocate a new array to store the result. See [1] for a probable future fix. So I think you should be able to speed this up by allocating r in advance to avoid allocating in each iteration. Then you would fill each element of r using a for loop (as I guess you're doing in C). > s = 0.; for i=1:d; 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)) * r > dE[:, k] += dJ > dE[:, n] -= dJ Same for all these element-wise operations. Have a look at NumericExtensions.jl for other ways of avoiding allocation of temporary arrays. Regards > end > end > return E, dE > end > 1: https://github.com/JuliaLang/julia/pull/5556
