我注意到我的几个项目中有一个小模式，其中我对一些代码进行了矢量化和并行化。最初的算法采用“推”的方式，而优化后的算法采用“拉”的方式。在本文中，我将描述我的意思，虽然这主要是为了展示一些漂亮的视频、图片和演示。
一个很好的起点是阿贝尔沙堆模型，就像我之前的许多模型一样，它在一段时间内完全吸引了我的注意。它是一个细胞自动机，每个细胞都是一堆沙粒-一个沙堆。在每一步，任何超过四个颗粒的沙堆都会将一个颗粒泄漏到它的四个4个相连的邻居中，而不考虑这些相邻单元格中的颗粒数量。边缘的细胞将它们的颗粒溢出到被遗忘的地方，这些颗粒就不复存在了。
随着多余的沙子落到边缘，模型最终达到不稳定状态，所有的桩都有三个或更少的颗粒。然而，直到灯光达到稳定，各种有趣的图案都在元胞自动机中荡漾。在某些情况下，最终的图案本身是美丽和有趣的。
Numberphile有一个很棒的视频，描述了如何在重复配置上形成一个组(也是)。简而言之，对于任何给定的网格大小，都有一个稳定的身份配置，当“添加”到组中的任何其他元素时，它将稳定回到该元素。身份结构本身就是一种分形，其本身一直是研究的热点。
计算身份配置实际上就是从某些开始配置开始运行模拟直到完成几次。以下是计算64x64身份配置过程的动画：
作为分形，网格越大，观察到的自相似图案就越多。网上有很多样本，我能找到的最大的就是维基百科上的3000x3000。但是我想看一个更大的，该死的！因此，跳到最后，我最终计算了这个10000x10000身份配置：
我使用OpenMP跨内核并行化，使用SIMD在线程内并行化。每个线程一次操作32个沙堆，为了计算身份沙堆，每个沙堆只需要3位状态，因此在相同的硬件上，这可能会增加到一次85个沙堆。输出格式是我以前的支柱Netpbm，包括视频输出。
那么，我说的推拉是什么意思？模拟沙堆的天真方法如下所示：
对于沙堆中的每个i，{if input[i]<；4{output[i]=input[i]}Else{output[i]=input[i]-4对于邻居中的每个j{output[j]=output[j]+1}。
当算法检查每个单元格时，它会将结果推入八个单元格。如果我们使用并发，这意味着执行的多线程可能会改变同一个单元，这需要同步锁、原子等。大量的同步是性能的丧钟。线程将花费所有时间争夺相同的资源，即使这只是虚假的共享。
对于沙堆中的每个i，{if input[i]<；4{output[i]=input[i]}对于邻居中的每个j，其他{output[i]=input[i]-4}{if input[j]>；=4{output[i]=output[i]+1}。
每个线程只修改一个单元-它负责更新的单元-所以不需要同步。它是着色器友好的，如果你看过我对康威的生活游戏的WebGL实现，应该听起来很熟悉。它本质上是相同的算法。如果你在网上搜索各种阿贝尔沙堆的参考资料，你最终会看到Cameron Fish在2017年发表的一篇关于在GPU上运行沙堆模拟的论文。他引用了我的WebGL生活游戏的文章，把一切都画上了圆圈。我们当时通过电子邮件进行了交谈，他与我分享了他的互动模拟。
向量化该算法很简单：一次加载多个桩，每个SIMD通道一个桩，并使用掩码实现分支。在我的代码中，我还展开了循环。为了避免SIMD代码中的边界检查，我用零填充状态数据结构，以便edgecells具有静态邻居，不再是特殊的。
回到过去，其中一个很酷的图形技巧是火焰动画。它很容易在有限的硬件上实现。事实上，计算它的最明显的方法是直接在帧缓冲区中，例如在VGA缓冲区中，没有外部状态。
屏幕底部有一个热源，算法自下而上运行，随机向上传播热量：
函数rand(min，max)//对于从底部开始的每个x，y，以[min，max]为单位的随机整数{buf[y-1][x+rand(-1，1)]=buf[y][x]-rand(0，1)}
作为推送算法，它在单线程上运行良好，但不能很好地转换到现代视频硬件上。因此，将其转换为Pull算法！
对于每个x，y{sx=x+rand(-1，1)sy=y+rand(1，2)output[y][x]=input[sy][sx]-rand(0，1)}。
细胞将火从底部向上拉起。尽管这一次有了acatch：这个算法将会有微妙的不同结果。
在原始版本中，只有一个状态缓冲器，因此火焰可以在一次传递中向上传播多次。我在这里的补偿是让火焰立刻传播得更远。
在原始版本中，火焰只传播到另一个细胞。在这个版本中，两个细胞可能从同一火焰中取出，克隆它。
最后很难区分，所以这就解决了问题。
rand()函数中仍然存在潜在的争用，但这可以通过接受x和y作为输入的散列函数来解决。
对这篇文章有什么评论吗？通过发送电子邮件至~Skeeto/[email protected][邮件列表礼仪]开始我的公共收件箱中的讨论，或查看现有讨论。