如何编写与维度无关的代码?

计算科学 算法 软件 正则
2021-11-26 21:37:21

我经常发现自己为给定操作/算法的一维、二维和三维版本编写了非常相似的代码。维护所有这些版本可能会变得乏味。简单的代码生成工作得相当好,但似乎认为必须有更好的方法。

有没有一种相对简单的方法来编写一次操作并将其推广到更高或更低的维度?

一个具体的例子是:假设我需要计算谱空间中速度场的梯度。在三个维度上,Fortran 循环看起来像:

do k = 1, n
  do j = 1, n
    do i = 1, n
      phi(i,j,k) = ddx(i)*u(i,j,k) + ddx(j)*v(i,j,k) + ddx(k)*w(i,j,k)
    end do
  end do
end do

ddx适当定义数组的位置(也可以用矩阵乘法来做到这一点。)二维流的代码几乎完全相同,除了:第三维从循环、索引和组件数量中删除。有没有更好的表达方式?

另一个例子是:假设我在三维网格上逐点定义流体速度。要将速度内插到任意位置(即不对应于网格点),可以在所有三个维度上连续使用一维Neville算法(即降维)。给定简单算法的一维实现,是否有一种简单的方法来进行降维?

4个回答

该问题强调大多数“普通”编程语言(至少是 C、Fortran)不允许您干净地执行此操作。一个额外的限制是您需要符号方便良好的性能。

因此,与其编写特定于维度的代码,不如考虑编写一个生成特定于维度的代码的代码。这个生成器是维度无关的,即使计算代码不是。换句话说,您在符号和表达计算的代码之间添加了一层推理。C++ 模板相当于同一件事:好的一面,它们是内置在语言中的。缺点是写起来有点麻烦。这减少了如何实际实现代码生成器的问题。

OpenCL 允许您在运行时相当干净地进行代码生成。它还在“外部控制程序”和“内部循环/内核”之间进行了非常清晰的划分。外部的生成程序对性能的限制要小得多,因此最好用一种舒适的语言编写,比如 Python。这就是我对PyOpenCL将如何被使用的希望——对于更新的无耻插件感到抱歉。

你看看 deal.II ( http://www.dealii.org/ ) 是如何做到的——在那里,维度独立性是库的核心,并且被建模为大多数数据类型的模板参数。例如,参见第 4 步教程程序中的与维度无关的拉普拉斯求解器:

http://www.dealii.org/developer/doxygen/deal.II/step_4.html

也可以看看

https://github.com/dealii/dealii/wiki/Frequently-Asked-Questions#why-use-templates-for-the-space-dimension

这可以用任何语言通过以下粗略的心理原型来完成:

  1. 创建每个维度的范围列表(我认为类似于 MATLAB 中的 shape())
  2. 在每个维度中创建您当前位置的列表。
  3. 在每个维度上编写一个循环,其中包含一个循环,该循环根据外部循环更改谁的大小。

从那里开始,这是一个与您的某种语言的语法作斗争以保持您的代码符合 nd 的问题。

在编写了一个n 维流体动力学求解器后,我发现拥有一种支持将类似对象的列表解包为函数参数的语言会很有帮助。即 a = (1,2,3) f(a*) -> f(1,2,3)。此外,高级迭代器(例如numpy 中的ndenumerate)使代码更简洁。

你可以写一个算法n1×n2×n3网格,然后通过设置一些专门来降低维度nj=1.