Parallel (GPU) algorithms for asynchronous cellular automata

Question

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:

1. Choose a neighbouring pair of cells at random

2. Calculate an "energy" function $\Delta E$ based on a local neighbourhood surrounding these cells

3. 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.

4. 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.

Answer 1

I would use the first option and would use a synchronous AC run before (using the GPU), to detect collisions, execute a step of a hexagonal AC whose rule is the value of the center cell = Sum (neighbors), This CA must have seven states should be initiated with randomly selected cell, and their status verified before running the update rule for each gpu.

Sample 1. The value of a neighboring cell is shared

0 0 0 0 0 0 0

  0 0 1 0 0 0

0 0 0 0 0 0 0

  0 0 0 1 0 0

0 0 0 0 0 0 0

a step of a CA whose rule is hexagonal central cell = Sum (neighbors)

0 0 1 1 0 0 0

  0 1 1 1 0 0

0 0 1 2 1 0 0

  0 0 1 1 1 0

0 0 0 1 1 0 0

Sample 2. The value of a cell to update is taken into account as a neighbor on the other

0 0 0 0 0 0 0

  0 0 1 0 0 0

0 0 0 1 0 0 0

  0 0 0 0 0 0

0 0 0 0 0 0 0

After the iteration

0 0 1 1 0 0 0

  0 1 2 2 0 0

0 0 2 2 1 0 0

  0 0 1 1 0 0

0 0 0 0 0 0 0

Sample 3. There is no relationship

  0 0 0 0 0 0

0 0 1 0 0 0 0

  0 0 0 0 0 0

0 0 0 0 0 0 0

  0 0 0 1 0 0

0 0 0 0 0 0 0

After the iteration

  0 1 1 0 0 0

0 1 1 1 0 0 0

  0 1 1 0 0 0

0 0 0 1 1 0 0

  0 0 1 1 1 0

0 0 0 1 1 0 0

Answer 2

Following your answers to my questions in the comments above, I would suggest you try a lock-based approach in which each thread tries to lock-down the neighbourhood it will be updating before computing the actual update.

You can do this using the atomic operations provided for in CUDA, and an array of int containing the locks for each cell, e.g. lock. Each thread then does the following:

ci, cj = choose a pair at random.

int locked = 0;

/* Try to lock the cell ci. */
if ( atomicCAS( &lock[ci] , 0 , 1 ) == 0 ) {

    /* Try to lock the cell cj. */
    if ( atomicCAS( &lock[cj] , 0 , 1 ) == 0 ) {

        /* Now try to lock all the neigbourhood cells. */
        for ( cn = indices of all neighbours )
            if ( atomicCAS( &lock[cn] , 0 , 1 ) != 0 )
                break;

        /* If we hit a break above, we have to unroll all the locks. */
        if ( cn < number of neighbours ) {
            lock[ci] = 0;
            lock[cj] = 0;
            for ( int i = 0 ; i < cn ; i++ )
                lock[i] = 0;
            }

        /* Otherwise, we've successfully locked-down the neighbourhood. */
        else
            locked = 1;

        }

    /* Otherwise, back off. */
    else
        lock[ci] = 0;
    }

/* If we got everything locked-down... */
if ( locked ) {

    do whatever needs to be done...

    /* Release all the locks. */
    lock[ci] = 0;
    lock[cj] = 0;
    for ( int i = 0 ; i < cn ; i++ )
        lock[i] = 0;

    }

Note that this approach is probably not the most optimal, but it could provide an interesting starting point. If there are a lot of collisions between threads, i.e. one or more per 32 threads (as in one collision per warp), then there will be a fair bit of branch diversion. Also, atomic operations can be a bit slow, but since you are only doing compare-and-swap operations, it should scale ok.

The locking overhead may look intimidating, but it's really only a few assignments and branches, not much more.

Note also that I'm being fast and loose with notation in the loops of i over the neighbours.

Addendum: I was cavalier enough to assume that you could simply back off when to pairs collide. If this is not the case, then you can wrap everything as of the second line in a while-loop and add a break at the end of the final if-statement.

All the threads will then have to wait until the last one is done, but if collisions are rare, you should be able to get away with it.

Addendum 2: Do not be tempted to add calls to __syncthreads() anywhere in this code, especially it the looping version described in the previous addendum! Asynchronicity is essential in avoiding repeated collisions in the latter case.

Answer 3

I'm the lead developer of LibGeoDecomp. While I agree with vanCompute that you could emulate your ACA with a CA, you're right that this would not be very efficient, as only few cells in any given step are meant to be updated. This is indeed a very interesting application -- and fun to tinker with!

I'd suggest you to combine the solutions proposed by jlopez1967 and Pedro: Pedro's algorithm captures the parallelism well, but those atomic locks are awfully slow. The solution of jlopez1967 is elegant when it comes to detecting collisions, but checking all n cells, when only a smaller subset (I'll from now on assume that there is some parameter k which denotes the number of cells to be updated simultaneously) are active, is clearly prohibitive.

__global__ void markPoints(Cell *grid, int gridWidth, int *posX, int *posY)
{
    int id = blockIdx.x * blockDim.x + threadIdx.x;
    int x, y;
    generateRandomCoord(&x, &y);
    posX[id] = x;
    posY[id] = y;
    grid[y * gridWidth + x].flag = 1;
}

__global__ void checkPoints(Cell *grid, int gridWidth, int *posX, int *posY, bool *active)
{
    int id = blockIdx.x * blockDim.x + threadIdx.x;
    int x = posX[id];
    int y = posY[id];
    int markedNeighbors = 
        grid[(y - 1) * gridWidth + x + 0].flag +
        grid[(y - 1) * gridWidth + x + 1].flag +
        grid[(y + 0) * gridWidth + x - 1].flag +
        grid[(y + 0) * gridWidth + x + 1].flag +
        grid[(y + 1) * gridWidth + x + 0].flag +
        grid[(y + 1) * gridWidth + x + 1].flag;
    active[id] = (markedNeighbors > 0);
}


__global__ void update(Cell *grid, int gridWidth, int *posX, int *posY, bool *active)
{
    int id = blockIdx.x * blockDim.x + threadIdx.x;
    int x = posX[id];
    int y = posY[id];
    grid[y * gridWidth + x].flag = 0;
    if (active[id]) {
        // do your fancy stuff here
    }
}

int main() 
{
  // alloc grid here, update up to k cells simultaneously
  int n = 1024 * 1024;
  int k = 1234;
  for (;;) {
      markPoints<<<gridDim,blockDim>>>(grid, gridWidth, posX, posY);
      checkPoints<<<gridDim,blockDim>>>(grid, gridWidth, posX, posY, active);
      update<<<gridDim,blockDim>>>(grid, gridWidth, posX, posY, active);
  }
}

In the absence of good global synchronization on the GPU, you need to invoke multiple kernels for the different phases. On Nvidia's Kepler you can move even the main loop to the GPU, but I don't expect that to gain much.

The algorithms achieves a (configurable) degree of parallelism. I guess, the interesting question is whether the collisions will affect your random distribution when you increase k.

Answer 4

I suggest to you that you see this link http://www.wolfram.com/training/courses/hpc021.html approximately 14:15 minutes into the video of course, mathematica training where they make an implementation of a cellular automata using CUDA, from there and you can modify it.