计算科学 - 异步元胞自动机的并行 (GPU) 算法 - 吾爱随笔录

异步元胞自动机的并行 (GPU) 算法

计算科学并行计算蒙特卡洛显卡

2021-12-18 02:55:10

我有一组可以被描述为异步元胞自动机的计算模型。这些模型类似于 Ising 模型，但稍微复杂一些。似乎这些模型将受益于在 GPU 而不是 CPU 上运行。不幸的是，并行化这样一个模型并不是很简单，而且我也不清楚如何去做。我知道有关于这个主题的文献，但这似乎都是针对对算法复杂性的细节感兴趣的核心计算机科学家，而不是像我这样只想描述我可以实现的东西的人，并且因此，我觉得它相当难以理解。

为清楚起见，我并不是在寻找一种最佳算法，而是在寻找一种可以在 CUDA 中快速实现的算法，这种算法可能会显着加快我的 CPU 实现速度。在这个项目中，程序员时间比计算机时间更多地是一个限制因素。

我还应该澄清一下，异步元胞自动机与同步元胞自动机是完全不同的东西，用于并行化同步 CA 的技术（例如 Conway 的生命）不容易适应这个问题。不同之处在于，同步 CA 在每个时间步同时更新每个单元，而异步 CA 在每个时间步更新随机选择的本地区域，如下所述。

我希望并行化的模型是在由约 100000 个单元组成的格子（通常是六边形）上实现的（尽管我想使用更多），运行它们的非并行算法如下所示：

随机选择一对相邻的单元格
计算“能量”函数 $\Delta E$ 基于这些细胞周围的当地社区
概率取决于 $e^{-\beta \Delta E}$ （和 $\beta$ 一个参数），要么交换两个单元的状态，要么什么都不做。
无限期地重复上述步骤。

边界条件也有一些复杂性，但我想这些不会对并行化造成太大困难。

值得一提的是，我对这些系统的瞬态动力学感兴趣，而不仅仅是平衡状态，所以我需要与上述具有等效动力学的东西，而不仅仅是接近相同平衡分布的东西。（所以棋盘算法的变体不是我想要的。）

并行化上述算法的主要困难是冲突。因为所有的计算只依赖于晶格的局部区域，所以许多晶格站点可以并行更新，只要它们的邻域不重叠。问题是如何避免这种重叠。我可以想到几种方法，但我不知道哪种方法最好。这些如下：

使用 CPU 生成随机网格站点列表并检查冲突。当网格点的数量等于 GPU 处理器的数量时，或者如果检测到碰撞，则将每组坐标发送到 GPU 单元以更新相应的网格点。这将很容易实现，但可能不会加快速度，因为检查 CPU 上的冲突可能不会比在 CPU 上进行整个更新便宜得多。
将晶格划分为多个区域（每个 GPU 单元一个），并有一个 GPU 单元负责随机选择和更新其区域内的网格单元。但是这个想法有很多我不知道如何解决的问题，最明显的是当一个单位选择一个与其区域边缘重叠的邻域时应该发生什么。
近似系统如下：让时间以离散的步骤进行。把格子分成不同的根据一些预定义的方案在每个时间步上设置一组区域，并让每个 GPU 单元随机选择和更新一对邻域不与区域边界重叠的网格单元。由于边界在每个时间步都发生变化，因此只要区域相对较大，此约束可能不会对动态产生太大影响。这似乎很容易实现并且可能很快，但我不知道它将如何近似动态，或者在每个时间步上选择区域边界的最佳方案是什么。我发现了一些对“块同步元胞自动机”的引用，这可能与这个想法相同，也可能不同。（我不知道，因为该方法的所有描述似乎都是俄语或我无法访问的来源。）

我的具体问题如下：

上述任何算法是否是处理异步 CA 模型的 GPU 并行化的明智方法？
有没有更好的办法？
是否有针对此类问题的现有库代码？
我在哪里可以找到“块同步”方法的清晰英文描述？

进步

我相信我已经想出了一种方法来并行化可能合适的异步 CA。下面概述的算法适用于一次只更新一个单元格的普通异步 CA，而不是像我的那样更新相邻的一对单元格。将其推广到我的具体案例存在一些问题，但我想我知道如何解决它们。但是，由于下面讨论的原因，我不确定它会带来多少速度优势。

这个想法是用行为等效的随机同步 CA (SCA) 替换异步 CA（以下简称 ACA）。为此，我们首先假设 ACA 是一个泊松过程。也就是说，时间连续进行，并且每个单元作为每单位时间执行其更新功能的恒定概率，独立于其他单元。

我们构建了一个 SCA，它的每个单元都存储两个东西：状态 $X_{ij}$ 单元格（即在顺序实现中将存储在每个单元格中的数据）和一个浮点数 $t_{ij}$ 表示下一次更新的（连续）时间。该连续时间与 SCA 的更新步骤不对应。我将后者称为“逻辑时间”。时间值根据指数分布随机初始化： $t_{ij}(0) \sim \operatorname{Exp}(\lambda)$ . （在哪里 $\lambda$ 是一个参数，其值可以任意选择。）

在每个逻辑时间步，SCA 的单元更新如下：

如果，对于任何 $k, l$ 在附近 $i,j$ ，时间 $t_{kl}<t_{ij}$ ，没做什么。
否则，（1）更新状态 $X_{ij}$ 根据各州 $X_{kl}$ 相邻小区，使用与原始 ACA 相同的规则；(2) 生成一个随机值 $\Delta t \sim \operatorname{Exp}(\lambda)$ 并更新 $t_{ij}$ 到 $t_{ij}+\Delta t$ .

我相信这可以保证单元格将按照可以“解码”以对应于原始 ACA 的顺序进行更新，同时避免冲突并允许并行更新某些单元格。但是，由于上面的第一个要点，这意味着大多数 GPU 处理器在 SCA 的每个时间步上大部分都处于空闲状态，这并不理想。

我需要多思考一下这个算法的性能是否可以提高，以及如何扩展这个算法来处理ACA中多个cell同时更新的情况。然而，它看起来很有希望，所以我想我会在这里描述它，以防任何人 (a) 知道文献中的任何类似内容，或者 (b) 可以提供对这些剩余问题的任何见解。

4个回答

我将使用第一个选项，并在之前使用同步 AC 运行（使用 GPU）来检测碰撞，执行一个六边形 AC 的步骤，其规则是中心单元的值 = Sum（邻居），这个 CA 必须有应使用随机选择的单元启动七个状态，并在为每个 gpu 运行更新规则之前验证它们的状态。

示例 1. 共享相邻单元格的值

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

CA 的一个步骤，其规则是六边形中心单元 = 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

示例 2. 将要更新的单元格的值视为另一个单元格的邻居

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

迭代后

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

示例 3. 没有关系

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

迭代后

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

根据您在上述评论中对我的问题的回答，我建议您尝试一种基于锁的方法，其中每个线程在计算实际更新之前尝试锁定它将更新的邻域。

int您可以使用 CUDA 中提供的原子操作以及包含每个单元格的锁的数组来执行此操作，例如lock. 然后每个线程执行以下操作：

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;

    }

请注意，这种方法可能不是最佳的，但它可以提供一个有趣的起点。如果线程之间有很多冲突，即每 32 个线程一个或多个（如每个 warp 一个冲突），那么就会有相当多的分支转移。此外，原子操作可能有点慢，但由于您只进行比较和交换操作，它应该可以扩展。

锁定开销可能看起来很吓人，但实际上只有几个分配和分支，不多。

另请注意，我在i邻居的循环中使用符号快速而松散。

附录：我很傲慢地假设你可以简单地在配对碰撞时退缩。如果不是这种情况，那么您可以将第二行的所有内容都包含在一个while-loop 中，并在最后一个-statementbreak的末尾添加 a。if

然后所有线程都必须等到最后一个线程完成，但如果冲突很少，你应该能够摆脱它。

附录 2：不要试图在这段代码中的任何地方添加调用，__syncthreads()尤其是之前附录中描述的循环版本！在后一种情况下，异步性对于避免重复冲突至关重要。

我是 LibGeoDecomp 的首席开发人员。虽然我同意 vanCompute 的观点，即您可以使用 CA 模拟您的 ACA，但您说得对，这不会非常有效，因为在任何给定步骤中只有少数单元格需要更新。这确实是一个非常有趣的应用程序——而且修改起来很有趣！

我建议您结合 jlopez1967 和 Pedro 提出的解决方案：Pedro 的算法很好地捕获了并行性，但是那些原子锁非常慢。jlopez1967 的解决方案在检测碰撞时很优雅，但是检查所有n单元格时，只有一个较小的子集（我从现在开始假设有一些参数k表示要同时更新的单元格的数量）处于活动状态，显然是禁止的。

__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);
  }
}

在 GPU 上没有良好的全局同步的情况下，您需要为不同的阶段调用多个内核。在 Nvidia 的 Kepler 上，您甚至可以将主循环移至 GPU，但我不认为这会获得太多收益。

这些算法实现了（可配置的）并行度。我想，有趣的问题是当你增加时碰撞是否会影响你的随机分布k。

我建议你看到这个链接http://www.wolfram.com/training/courses/hpc021.html大约 14:15 分钟到视频当然，mathematica 培训他们使用 CUDA 实现细胞自动机, 从那里你可以修改它。

其它你可能感兴趣的问题

上一篇快速找到点集中的粗线下一篇给定一个 SPD 三对角线性系统，我们可以预先计算以使任何三个索引都可以在 O(1) 时间内链接起来吗？