[Numpy-discussion] Rookie problems - Why is C-code much faster?
Mads Ipsen
mpi at osc.kiku.dk
Wed Feb 22 03:56:04 CST 2006
On Tue, 21 Feb 2006, Tim Hochberg wrote:
> Mads Ipsen wrote:
>
> >On Tue, 21 Feb 2006, Tim Hochberg wrote:
> >
> >
> >
> >>This all makes perfect sense, but what happended to box? In your
> >>original code there was a step where you did some mumbo jumbo and box
> >>and rint. Namely:
> >>
> >>
> >
> >It's a minor detail, but the reason for this is the following
> >
> >Suppose you have a line with length of box = 10 with periodic boundary
> >conditions (basically this is a circle). Now consider two points x0 = 1
> >and x1 = 9 on this line. The shortest distance dx between the points x0
> >and x1 is dx = -2 and not 8. The calculation
> >
> > dx = x1 - x0 ( = +8)
> > dx -= box*rint(dx/box) ( = -2)
> >
> >will give you the desired result, namely dx = -2. Hope this makes better
> >sense. Note that fmod() won't work since
> >
> > fmod(dx,box) = 8
> >
> >
> I think you could use some variation like "fmod(dx+box/2, box) - box/2"
> but rint seems better.
>
> >Part of my original post was concerned with the fact, that I initially was
> >using around() from numpy for this step. This was terribly slow, so I made
> >some custom changes and added rint() from the C-math library to the numpy
> >module, giving a speedup factor about 4 for this particular line in the
> >code.
> >
> >Best regards // Mads
> >
> >
> >
> OK, that all makes sense. You might want to try the following, which
> factors out all the divisions and half the multiplies by box and
> produces several fewer temporaries. Note I replaced x**2 with x*x,
> which for the moment is much faster (I don't know if you've been
> following the endless yacking about optimizing x**n, but x**2 will get
> fast eventually). Depending on what you're doing with r2, you may be
> able to avoid the last multiple by box as well.
>
>
> # Loop over all particles
> xbox = x/box
> ybox = y/box
> for i in range(n-1):
> dx = xbox[i+1:] - xbox[i]
> dy = ybox[i+1:] - ybox[i]
> dx -= rint(dx)
> dy -= rint(dy)
> r2 = (dx*dx + dy*dy)
> r2 *= box
>
>
>
> Regards,
>
> -tim
>
Thanks Tim,
I am only a factor 2.5 slower than the C loop now, thanks to your
suggestions.
// Mads
More information about the Numpy-discussion
mailing list