Hard to test without some way of generating data, but I'd guess the problem is
dE[:, k] += dJ dE[:, n] -= dJ which (unfortunately) actually expands to dE[:, k] = dE[:, k] + dJ dE[:, n] = dE[:, n] - dJ so there is a new array being created on the right-hand-side - it isn't an inplace operation. If you broke that out into a for loop, you should see a speed up. After that, run the profiler. 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. Finally - Julia 0.2 or 0.3? On Sunday, May 25, 2014 10:39:51 AM UTC-4, Christoph Ortner wrote: > > > 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] > > 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 > > end > > end > > return E, dE > > end > > >
