Fortran 类似于 Cython 中的数组切片

Question

我正在寻找一种简洁和高效的方法来获取数组的多维切片，对这些切片执行标量和矩阵运算，然后最终将结果数组保存为另一个数组中的切片。

您可以使用如下语法在 Fortran 中很好地完成此操作：

real*8, dimension(4,4,4,4)               ::  matrix_a
real*8, dimension(4,4)               ::  matrix_b
...
matrix_a(:, 2, :, 4) =  matrix_a(:, 2, :, 4) + (2 * matrix_b(:, :))

我正在尝试寻找在 Cython 中执行此操作的方法。这是我能想到的最好的：

cdef double matrix_a[4][4][4][4]
cdef double matrix_b[4][4]
...
cdef int i0, i1
for i0 in range(4):
  for i1 in range(4):
    matrix_a[i0][1][i1][3] += (2 * matrix_b[i0][i1])

显然，您可以使用 Numpy 数组以非常简洁的方式执行此操作，但这似乎比上面的代码慢几个数量级。有什么方法可以让我两全其美吗？可能我可以通过某种方式调用 Numpy 使其更快？或者可能是某种可以从 Cython 中利用的 C/C++ API。谢谢

Answer 1

好问题。简短的回答是，不。

长答案...这取决于您想要什么以及您愿意为此付出多少努力。没有单一的答案，因为问题非常广泛，而不是因为无法进行操作。

在撰写本文时，没有默认方法可以在 cython 的内存视图上执行 fortran 速度算术。在 Cython 邮件列表上对此进行了多次讨论。这个recent thread总结的很好

现在，有很多方法可以使用 NumPy 数组进行快速算术运算。没有任何一个具有其他所有功能，但有些可能适合您的需要。

第一个：

NumPy 的矢量化写得很好。在类型推断中有一些预先的开销，但是，除了在紧密循环内对小数组的操作之外的任何其他事情，大部分计算时间是运行在 C 中非常优化的循环中。如果你正在做像这样一个简单的就地操作，即使在 Cython 中，最好的方法可能就是做类似 np.add(a, b, out=c) 的事情。如果您有更具体的需求，np.einsum、scipy 的 blas/lapack 包装器和 numpy 的数组缩减方法提供更大的灵活性。

在大多数情况下，如果 a 和 b 是 numpy 数组，例如

a[:,2,:,4] += 2 * b

应该足够了。

即使您在 Cython 中操作内存视图而不是 numpy 数组，您也可以执行以下操作：

import numpy as np
np.add(a, b, out=c)

像这样调用 Python 函数有固定成本，但如果数组很大，则无所谓。

在您的情况下，numpy 仍然较慢的原因有两个。首先，开始循环操作的固定成本很高，因为 numpy 在运行时刻处理大量其类型和形状元数据。对于较大的数组，这种差异会消失，因此，除非您处于紧密循环中，否则这无关紧要。其次，numpy 和 fortran 的默认内存布局不同。这种差异可能只会出现在较大数组的操作中。如果您以 fortran 内存顺序初始化 numpy 数组并且数组很大，numpy 应该与 fortran 相同。

现在，对于那些确实确实处于 numpy 速度不够快或对数组的控制不够精细的人来说，有很多其他选项。以下是一些更有希望的：

在 Cython 的 cdef 函数中用循环编写算术运算。让它接受内存视图作为参数并传入预分配的输出数组。签名可能类似于 cdef inline void add(double[:,:] a, double[:,:] b, double[:,:] out) nogil: 可以在没有 Python 调用开销的情况下对内存视图进行切片，因此您可以将切片传递给您定义的函数，而不会产生 Python 调用的开销。这种方法不是最快的，但老实说，它通常已经足够好了。它还避免了必须手动管理数组的内存布局所带来的大部分痛苦。这也比用原始 C 数组做同样的事情更容易，因为内存视图的内存是动态管理的。
存在多个可识别 numpy 的自动向量化器。考虑尝试 numba, parakeet, or pythran。这些包主要集中在编译特定类型的函数版本，这些函数在 numpy 函数上运行，然后在适用的地方调用这些例程。我对 numba 和 parakeet 的研究结果好坏参半，对 pythran 的研究也不多。据我所知，numba 和 parakeet 目前不进行任何类型的静态表达式分析，尽管这不一定超出它们的范围。 Pythran 看起来很有前途。它的主要卖点是对数组表达式（如 Fortran）进行静态分析。这些库易于使用，因此值得尝试一下，看看它们是否有帮助。他们特别擅长优化循环。
您还可以使用一些 Python 级表达式优化器。参见 numexpr and Theano。这些包，粗略地说，专注于将算术表达式编译成机器代码，而不是整个函数。据我所知，它们没有针对切片进行优化，但您可以将 numpy 数组的切片传递给它们。这些最擅长优化较大数组上的操作，但它们通常可以提供比使用 numpy 数组的直接算术更快的表达式代码。
你可以用fortran实现操作，用Cython封装。在 fortran90 website. I wrote an example of how to do this some time ago in this answer 中描述了如何使用 Fortran 执行此操作。您必须自己管理内存布局，但这并不难，而且会非常快。
您还可以使用各种 C++ 库中的任何一种进行数组操作。不幸的是，它们中的许多（eigen, armadillo, blaze-lib）目前没有对高维数组的内置支持，但您拥有的特定表达式实际上只是跨步二维数组的算术运算，因此它可以正常工作。如果您不介意这个限制，已经有一些初步努力为犰狳制作 Cython 绑定（cy-arma，尽管它们仍未开发。一个支持 n 维数组操作的 c++ 库是 libdynd。它还有一组用 Cython 编写的更好的 python bindings。用 C++ 实现这些操作并通过 Cython 为 Python 包装它们可能是最简单的，但您可能会从包装器中取出 C++ pxd 声明并尝试在 Cython 级别编写您的函数。
有一个名为ceygen的库，它使用eigen作为后端，将内存视图上的算术运算编写为函数调用。我在编译时遇到了麻烦，但从概念上讲，这可能就是您要找的东西。我不确定它是否能很好地处理跨步信息。
您可以自己在 C 中实现这些操作。唯一更好的情况是当您想将 sse 内在函数应用于数组时（或某种比纯循环更好的优化） ).写这样的东西非常耗时，但你可以做到。如果你只是循环，你最好在 Cython 中编写循环并在内存视图上运行。
您可以将 NumPy C API 用于 ufuncs or gufuncs。这允许您进行标量（或数组）操作并生成一个将其应用于更一般数组的函数。目前文档还不是很完整。
这不适用于 Cython 中的基本算术运算，但很快将添加到 scipy 的新功能之一是用于 BLAS 和 LAPACK 的 Cython 级包装器.这将使您能够进行各种线性代数运算，而无需为每个函数调用付出很少的代价。现在，这里是相关的 pull request.

其中一些需要您自己管理内存布局，但这并不难。

如果我处在你的位置，我会尝试选项 1，然后是选项 4。选项 1 会相当快，但如果这还不够，选项 4 可能会更快，尽管现在你将拥有自己担心数组的内存布局。

Fortran 类似于 Cython 中的数组切片

Fortran like array slicing in Cython

python

arrays

fortran

cython

fortran90