Take the 2-minute tour ×

I'm doing a Monte-Carlo simulation of some particles particles. There are several bottle necks in my code but the main one is that in some of the tries I make, I need to update all the particles properties. The code is written in c++ and currently I have several loops to achieve that:
1. a loop to store the old properties of all the particles and update the new properties.
2. a 2D loop of interactions.
3. another 2D loop of interactions (I can't combine it with the first one).
4. a loop to store accept the step/a loop to reject the step.

I am hoping to remove step 4 using swap but I can't find a way to do so. All the particles are a class which has several elements named properties and nextProperties or oldProperties. How would you approach that?

share|improve this question
This is very hard to answer without seeing some code. Can produce a simple example that demonstrates your problem? –  Björn Pollex Oct 20 '11 at 7:44
You could try timing specific sections and seeing where, exactly, your bottleneck is. Once you've found that, posting the slow code will help the good people at SO help you find a solution. –  Carl Winder Oct 20 '11 at 8:00

3 Answers 3

It sounds like you could use double buffering. Basically, you'd maintain pointers to two arrays of particle objects — call them, say, accepted and trial. At the beginning of a trial, you copy the properties of the particles on the accepted array to those on the trial array, and make any modifications you want. If the trial is successful, you then just swap the pointers, so that what used to be the trial array becomes accepted and vice versa.

Also, you say that only some of your trials involve costly updates. If so, you might be interested in techniques like fast variable dragging or ensemble updating.

share|improve this answer
I don't how your c++ skills are, do you know if I can hold a pointer to iterator? –  Yotam Oct 20 '11 at 17:28
If you have vector<Foo> v1 and vector<Foo> v2, then v1.swap(v2) should be cheap. –  msandiford Oct 20 '11 at 22:15

Steps #2 and #3 are going to be O(N^2), whereas #1 and #4 are O(N), so you should focus on those.

If you absolutely need to do a calculation between every pair of particles, there isn't much you can do. But most likely, particles that are more than a certain distance apart will have negligible affect on each other, or perhaps you only need to deal with the nearest k particles (for a fixed k). In this case, an octtree (or kd-tree, aabb tree, or something similar) is your best bet for reducing the number of pair-wise calculations.

In particular, you might want to look into the Barnes-Hut method, which is used to reduce the complexity of gravitational calculations.

share|improve this answer
Step #3 is electrostatic which is similar to gravity in many ways. In step 4 I'm doing a lot of copying which, I think, are time consuming, I'll look on your method though. –  Yotam Oct 20 '11 at 16:46

I think the main bottleneck is your 2D interaction. If you are doing N^2 interactions, you need to probably think if you can work with a mean-field approximation, basically splitting the computational domain in cells, computing a mean-field for each cell for each iteration, and then only compute the interactions with particles in the same cell, plus a mean field correction for surrounding cells.

Here you can read more details regarding this optimization and when it is relevant

share|improve this answer
The first N^2 step is handling this cell wise. The second not (and I'm pretty sure that this is the real bottleneck). However, I need to check how acceptable is this for my purposes. –  Yotam Oct 20 '11 at 17:08
it depends on the reason why you are doing the second N^2 pass and what the pass is doing, then we can help you. Can you explain this step with more detail? –  lurscher Oct 20 '11 at 17:15
Step #3 is summing over all the charged particles in my system. When there are 1000 of those, the code is really slow. I'm also using a method called Lekner summation to deal with periodic property of the system. This method force me to use a map. –  Yotam Oct 20 '11 at 17:55
is there a reason you are not using cells to compute step #3? if you are trying to get an error estimator, you could increase the threshold minimal cells around each particle where the mean-field approximation kicks in, and use the difference to get your estimation –  lurscher Oct 20 '11 at 18:01
‎The main reason I'm not using the cells method is that I need to ask my adviser if that's O.K. This will only happen at Sunday though. –  Yotam Oct 20 '11 at 18:44

Your Answer


