I have a collection of computational models that could be described as asynchronous cellular automata. These models resemble the Ising model, but are slightly more complicated. It seems as if such models would benefit from being run on a GPU rather than a CPU. Unfortunately it isn't quite straightforward to parallelise such a model, and it isn't at all clear to me how to go about it. I'm aware that there is literature on the subject, but it all seems to be aimed at hardcore computer scientists who are interested in the details of algorithmic complexity, rather than someone like me who just wants a description of something I can implement, and consequently I find it rather inpenetrable.
For clarity, I'm not looking for an optimal algorithm so much as something I can rapidly implement in CUDA that's likely to give a significant speedup over my CPU implementation. Programmer time is much more of a limiting factor than computer time in this project.
I should also clarify that an asynchronous cellular automaton is a rather different thing from a synchronous one, and techniques for parallelising synchronous CAs (such as Conway's life) cannot easily be adapted to this problem. The difference is that a synchronous CA updates every cell simultaneously at every time step, whereas an asynchronous one updates a randomly chosen local region at every time step as outlined below.
The models I wish to parallelise are implemented on a lattice (usually a hexagonal one) consisting of ~100000 cells (though I'd like to use more), and the non-parallelised algorithm for running them looks like this:
Choose a neighbouring pair of cells at random
Calculate an "energy" function $\Delta E$ based on a local neighbourhood surrounding these cells
With a probability that depends on $e^{-\beta \Delta E}$ (with $\beta$ a parameter), either swap the states of the two cells or do nothing.
Repeat the above steps indefinitely.
There are also some complications to do with the boundary conditions, but I imagine these won't pose much difficulty for parallelisation.
It's worth mentioning that I'm interested in the transient dynamics of these systems rather than just the equilibrium state, so I need something that has equivalent dynamics to the above, rather than just something that will approach the same equilibrium distribution. (So variations of the chequerboard algorithm are not what I'm looking for.)
The main difficulty in parallelising the above algorithm is collisions. Because all the calculations depend only on a local region of the lattice, it's possible for many lattice sites to be updated in parallel, as long as their neighbourhoods aren't overlapping. The question is how to avoid such overlaps. I can think of several ways, but I don't know which if any is the best one to implement. These are as follows:
Use the CPU to generate a list of random grid sites and check for collisions. When the number of grid sites equals the number of GPU processors, or if a collision is detected, send each set of coordinates to a GPU unit to update the corresponding grid site. This would be easy to implement but probably wouldn't give much of a speed up, since checking for collisions on the CPU would probably not be all that much cheaper than doing the whole update on the CPU.
Divide the lattice up into regions (one per GPU unit), and have one GPU unit responsible for randomly selecting and updating grid cells within its region. But there are many issues with this idea that I don't know how to solve, the most obvious being what exactly should happen when a unit chooses a neighbourhood overlapping the edge of its region.
Approximate the system as follows: let time proceed in discrete steps. Divide the lattice up into a different set of regions on every time step according to some pre-defined scheme, and have each GPU unit randomly select and update a pair of grid cells whose neighbourhood does not overlap the region's boundary. Since the boundaries change every time step this constraint might not affect the dynamics too much, as long as the regions are relatively large. This seems easy to implement and likely to be fast, but I don't know how well it will approximate the dynamics, or what is the best scheme for choosing the region boundaries on each time step. I found some references to "block-synchronous cellular automata", which may or may not be the same as this idea. (I don't know because it seems that all the descriptions of the method are either in Russian or are in sources to which I don't have access.)
My specific questions are as follows:
Are any of the above algorithms a sensible way to approach GPU parallelisation of an asynchronous CA model?
Is there a better way?
Is there existing library code for this type of problem?
Where can I find a clear English-language description of the "block-synchronous" method?
Progress
I believe I have come up with a way to parallelise an asynchronous CA that might be suitable. The algorithm outlined below is for a normal asynchronous CA that updates only one cell at a time, rather than a neighbouring pair of cells as mine does. There are some issues with generalising it to my specific case, but I think I have an idea how to resolve them. However, I'm not sure how much of a speed benefit it will give, for reasons discussed below.
The idea is to replace the asynchronous CA (henceforth ACA) with a stochastic synchronous CA (SCA) that behaves equivalently. To do this we first imagine that the ACA is a Poisson process. That is, time proceeds continuously, and each cell as a constant probability per unit time of performing its update function, independently of the other cells.
We construct an SCA whose cells each store two things: the state $X_{ij}$ of the cell (i.e. the data that would be stored in each cell in the sequential implementation), and a floating point number $t_{ij}$ representing the (continuous) time at which it will next update. This continuous time does not correspond to the update steps of the SCA. I will refer to the latter as "logical time". The time values are initialised randomly according to an exponential distribution: $t_{ij}(0) \sim \operatorname{Exp}(\lambda)$. (Where $\lambda$ is a parameter whose value can be chosen arbitrarily.)
At each logical time step, the cells of the SCA are updated as follows:
If, for any $k, l$ in the neighbourhood of $i,j$, the time $t_{kl}<t_{ij}$, do nothing.
Otherwise, (1) update the state $X_{ij}$ according to the states $X_{kl}$ of the neighbouring cells, using the same rule as the original ACA; and (2) generate a random value $\Delta t \sim \operatorname{Exp}(\lambda)$ and update $t_{ij}$ to $t_{ij}+\Delta t$.
I believe this guarantees that the cells will be updated in an order that can be "decoded" to correspond to the original ACA, while avoiding collisions and allowing some cells to be updated in parallel. However, because of the first bullet point above, it means that most of the GPU processors will be mostly idle on each time step of the SCA, which is less than ideal.
I need to give some more thought to whether the performance of this algorithm can be improved, and how to extend this algorithm to deal with the case where multiple cells are updated simultaneously in the ACA. However, it looks promising so I thought I would describe it here in case anyone (a) knows of anything similar in the literature, or (b) can offer any insight into these remaining issues.
exp()
) so I wouldn't have thought it makes much sense to spread it over multiple threads. I think it's better (and easier for me) to try and update multiple pairs in parallel, with one pair per thread. $\endgroup$