NumPy用户指南 >使用NumPy C-API > 使用Python作为胶水
许多人喜欢说Python是一种很棒的粘合语言。希望本章将使您确信这是事实。最初将Python运用于科学的人通常是使用Python将超级计算机上运行的大型应用程序代码粘合在一起的人。用Python编写的代码不仅比Shell脚本或Perl更好,此外,易于扩展Python的能力使创建专门用于解决问题的新类和类型相对容易。从这些早期贡献者的互动中,Numeric成为了一个类似于数组的对象,可用于在这些应用程序之间传递数据。
随着Numeric的成熟和发展成为NumPy,人们已经能够直接在NumPy中编写更多代码。通常,此代码可以快速用于生产环境,但是有时仍然需要访问已编译的代码。要么从算法中获得最后的效率,要么使访问以C / C ++或Fortran编写的广泛使用的代码更容易。
本章将回顾许多可用于访问以其他编译语言编写的代码的工具。有许多资源可用于学习从Python调用其他已编译的库,本章的目的不是使您成为专家。主要目标是使您了解其中的一些可能性,以便您了解“ Google”的知识,以了解更多信息。
虽然Python是一种伟大的语言和高兴的,它的动态性结果代码的开销,可能会导致一些代码(即 for循环中的原始计算)比用静态编译语言编写的等效代码慢10-100倍。此外,由于在计算过程中创建和销毁了临时数组,这可能导致内存使用量超出必要数量。对于许多类型的计算需求,通常无法避免额外的速度降低和内存消耗(至少对于代码中时间或内存至关重要的部分)。因此,最常见的需求之一是从Python代码调出快速的机器代码例程(例如,使用C / C ++或Fortran进行编译)。这相对容易做到这一事实,是Python成为科学和工程编程如此出色的高级语言的重要原因。
它们是调用已编译代码的两种基本方法:编写扩展模块,然后使用import命令将其导入到Python,或使用ctypes 模块直接从Python调用共享库子例程。编写扩展模块是最常见的方法。
警告
如果不小心,从Python调用C代码可能会导致Python崩溃。本章中的方法都不是免疫的。您必须了解有关NumPy和所使用的第三方库处理数据的方式的一些知识。
在编写扩展模块中讨论了扩展模块。与已编译代码交互的最基本方法是编写扩展模块并构造一个调用已编译代码的模块方法。为了提高可读性,您的方法应利用该PyArg_ParseTuple
调用在Python对象和C数据类型之间进行转换。对于标准C数据类型,可能已经有一个内置转换器。对于其他人,您可能需要编写自己的转换器并使用"O&"
格式字符串,该字符串允许您指定一个函数,该函数将用于执行从Python对象到所需C结构的转换。
一旦执行了对适当的C结构和C数据类型的转换,包装器的下一步就是调用基础函数。如果基础函数是C或C ++,则这很简单。但是,为了调用Fortran代码,您必须熟悉如何使用编译器和平台从C / C ++调用Fortran子例程。这可能会因平台和编译器的不同而有所差异(这是f2py使接口Fortran代码的工作变得更加简单的另一个原因),但通常涉及名称的下划线修饰和所有变量都是通过引用传递(即,所有参数都是指针)的事实。
手工生成的包装器的优点是,您可以完全控制C库的使用和调用方式,从而可以实现精简且紧密的接口,并且开销最少。缺点是您必须编写,调试和维护C代码,尽管大多数代码都可以使用来自其他扩展模块的“剪切,粘贴和修改”的悠久技术进行修改。因为调用附加C代码的过程相当合理,所以已经开发了代码生成过程来简化此过程。这些代码生成技术之一随NumPy一起分发,并允许与Fortran和(简单)C代码轻松集成。下一部分将简要介绍此软件包f2py。
F2py允许您自动构造一个扩展模块,该模块与Fortran 77/90/95代码中的例程相接口。它具有解析Fortran 77/90/95代码并自动为其遇到的子例程生成Python签名的能力,或者您可以通过构造接口定义文件(或修改由f2py生成的一个)来指导子例程与Python的接口。 )。
引入f2py的最简单方法可能是提供一个简单示例。这是名为的文件中包含的子例程之一
add.f
:
C
SUBROUTINE ZADD(A,B,C,N)
C
DOUBLE COMPLEX A(*)
DOUBLE COMPLEX B(*)
DOUBLE COMPLEX C(*)
INTEGER N
DO 20 J = 1, N
C(J) = A(J)+B(J)
20 CONTINUE
END
该例程仅将元素添加到两个连续的数组中,然后将结果放在第三个数组中。所有三个阵列的内存必须由调用例程提供。f2py可以自动生成此例程的一个非常基本的接口:
f2py -m add add.f
假设您的搜索路径设置正确,您应该能够运行此命令。此命令将在当前目录中生成一个名为addmodule.c的扩展模块。现在可以像其他任何扩展模块一样从Python编译和使用此扩展模块。
您还可以让f2py编译add.f并编译其产生的扩展模块,仅留下可以从Python导入的共享库扩展文件:
f2py -c -m add add.f
此命令在当前目录中保留一个名为add。{ext}的文件(其中{ext}是平台上python扩展模块的适当扩展名-pyd 等)。然后可以从Python导入此模块。它将为add中的每个子例程包含一个方法(zadd,cadd,dadd,sadd)。每个方法的文档字符串包含有关如何调用模块方法的信息:
>>> import add
>>> print add.zadd.__doc__
zadd - Function signature:
zadd(a,b,c,n)
Required arguments:
a : input rank-1 array('D') with bounds (*)
b : input rank-1 array('D') with bounds (*)
c : input rank-1 array('D') with bounds (*)
n : input int
默认接口是将fortran代码非常直观地翻译成Python。现在,Fortran数组参数必须是NumPy数组,整数参数应该是整数。接口将尝试将所有参数转换为所需的类型(和形状),如果不成功,则发出错误。但是,由于它对参数的语义一无所知(例如C是输出,n应该真正匹配数组的大小),因此有可能以可能导致Python崩溃的方式滥用此函数。例如:
>>> add.zadd([1,2,3], [1,2], [3,4], 1000)
在大多数系统上都会导致程序崩溃。在幕后,这些列表将被转换为适当的数组,但随后告知底层的add循环以循环的方式超出已分配内存的边界。
为了改善界面,应提供指令。这是通过构造接口定义文件来完成的。通常最好从f2py可以生成的接口文件开始(从中获取其默认行为)。要使f2py生成接口文件,请使用-h选项:
f2py -h add.pyf -m add add.f
此命令将文件add.pyf保留在当前目录中。该文件对应于zadd的部分是:
subroutine zadd(a,b,c,n) ! in :add:add.f
double complex dimension(*) :: a
double complex dimension(*) :: b
double complex dimension(*) :: c
integer :: n
end subroutine zadd
通过放置intent指令并检查代码,可以相当多地清理接口,直到Python模块方法更易于使用和更强大为止。
subroutine zadd(a,b,c,n) ! in :add:add.f
double complex dimension(n) :: a
double complex dimension(n) :: b
double complex intent(out),dimension(n) :: c
integer intent(hide),depend(a) :: n=len(a)
end subroutine zadd
Intent指令intent(out)用于告诉f2py这c
是一个输出变量,应该在传递给基础代码之前由接口创建。intent(hide)指令告诉f2py不允许用户指定变量n
,而是从的大小获取它a
。将取决于(a
)指令是必要告诉f2py n的值取决于输入(所以它不会尝试创建变量n,直至创建一个变量)。
修改后add.pyf
,可以通过add.f95
和编译生成新的python模块文件add.pyf
:
f2py -c add.pyf add.f95
新接口具有文档字符串:
>>> import add
>>> print add.zadd.__doc__
zadd - Function signature:
c = zadd(a,b)
Required arguments:
a : input rank-1 array('D') with bounds (n)
b : input rank-1 array('D') with bounds (n)
Return objects:
c : rank-1 array('D') with bounds (n)
现在,可以以更强大的方式调用该函数:
>>> add.zadd([1,2,3],[4,5,6])
array([ 5.+0.j, 7.+0.j, 9.+0.j])
注意自动转换为正确的格式。
还可以通过将变量指令作为特殊注释放入原始fortran代码中来自动生成nice接口。因此,如果我将源代码修改为包含:
C
SUBROUTINE ZADD(A,B,C,N)
C
CF2PY INTENT(OUT) :: C
CF2PY INTENT(HIDE) :: N
CF2PY DOUBLE COMPLEX :: A(N)
CF2PY DOUBLE COMPLEX :: B(N)
CF2PY DOUBLE COMPLEX :: C(N)
DOUBLE COMPLEX A(*)
DOUBLE COMPLEX B(*)
DOUBLE COMPLEX C(*)
INTEGER N
DO 20 J = 1, N
C(J) = A(J) + B(J)
20 CONTINUE
END
然后,我可以使用以下命令编译扩展模块:
f2py -c -m add add.f
函数add.zadd的结果签名与之前创建的签名完全相同。如果原始源代码包含和A(N)
而不是and ,则只需将注释行放在源代码中,就可以获得(几乎)相同的接口
。唯一的区别是这是一个可选输入,默认为。A(*)
B
C
INTENT(OUT) :: C
N
A
为了与要讨论的其他方法进行比较。这是使用固定平均滤波器对双精度浮点数的二维数组进行滤波的函数的另一个示例。从此示例中应该清楚知道使用Fortran索引到多维数组的优点。
SUBROUTINE DFILTER2D(A,B,M,N)
C
DOUBLE PRECISION A(M,N)
DOUBLE PRECISION B(M,N)
INTEGER N, M
CF2PY INTENT(OUT) :: B
CF2PY INTENT(HIDE) :: N
CF2PY INTENT(HIDE) :: M
DO 20 I = 2,M-1
DO 40 J=2,N-1
B(I,J) = A(I,J) +
$ (A(I-1,J)+A(I+1,J) +
$ A(I,J-1)+A(I,J+1) )*0.5D0 +
$ (A(I-1,J-1) + A(I-1,J+1) +
$ A(I+1,J-1) + A(I+1,J+1))*0.25D0
40 CONTINUE
20 CONTINUE
END
可以使用以下代码编译此代码并将其链接到名为filter的扩展模块中:
f2py -c -m filter filter.f
这将在当前目录中生成一个名为filter.so的扩展模块,该模块具有一个名为dfilter2d的方法,该方法返回输入的过滤版本。
f2py程序是用Python编写的,可以在您的代码内部运行以在运行时编译Fortran代码,如下所示:
from numpy import f2py
with open("add.f") as sourcefile:
sourcecode = sourcefile.read()
f2py.compile(sourcecode, modulename='add')
import add
源字符串可以是任何有效的Fortran代码。如果要保存扩展模块源代码,则可以通过source_fn
关键字为编译功能提供合适的文件名。
如果要分发f2py扩展模块,则只需包括.pyf文件和Fortran代码。NumPy中的distutils扩展允许您完全根据此接口文件定义扩展模块。setup.py
允许分发add.f
模块的有效文件(作为软件包的一部分,
f2py_examples
以便作为加载f2py_examples.add
)是:
def configuration(parent_package='', top_path=None)
from numpy.distutils.misc_util import Configuration
config = Configuration('f2py_examples',parent_package, top_path)
config.add_extension('add', sources=['add.pyf','add.f'])
return config
if __name__ == '__main__':
from numpy.distutils.core import setup
setup(**configuration(top_path='').todict())
使用以下命令即可轻松安装新软件包:
pip install .
假设您具有使用您所使用的Python版本写入主sitepackage目录的适当权限。为了使生成的包起作用,您需要创建一个名为的文件__init__.py
(与相同的目录中add.pyf
)。请注意,扩展模块完全根据add.pyf
和add.f
文件定义。.pyf文件到.c文件的转换由numpy.disutils处理。
接口定义文件(.pyf)是您如何微调Python和Fortran之间的接口。在numpy / f2py / docs目录中找到关于f2py的体面文档,无论您在系统上安装了NumPy的位置如何(通常在站点软件包下)。在https://scipy-cookbook.readthedocs.io的“与其他语言交互”标题下,还提供了有关使用f2py的更多信息(包括如何使用它包装C代码)。
链接已编译代码的f2py方法是目前最复杂和集成的方法。它允许将Python与已编译的代码完全分开,同时仍允许扩展模块的单独分发。唯一的缺点是它需要存在Fortran编译器才能使用户安装代码。但是,由于存在免费编译器g77,gfortran和g95以及高质量的商用编译器,因此这种限制并不是特别繁重。我认为,Fortran仍然是编写快速清晰的代码以进行科学计算的最简单方法。它以最直接的方式处理复数和多维索引。但是请注意,某些Fortran编译器将无法优化代码以及良好的手写C代码。
Cython是Python方言的编译器,它添加(可选)静态类型以提高速度,并允许将C或C ++代码混合到模块中。它产生可以用Python代码编译和导入的C或C ++扩展。
如果您正在编写一个扩展模块,该模块还将包含您自己的很多算法代码,那么Cython是一个很好的选择。它的功能之一就是能够轻松快速地处理多维数组。
请注意,Cython仅是扩展模块生成器。与f2py不同,它不包含用于编译和链接扩展模块的自动工具(必须以通常的方式完成)。它确实提供了一个改进的distutils类,该类build_ext
使您可以从.pyx
源代码构建扩展模块。因此,您可以写一个setup.py
文件:
from Cython.Distutils import build_ext
from distutils.extension import Extension
from distutils.core import setup
import numpy
setup(name='mine', description='Nothing',
ext_modules=[Extension('filter', ['filter.pyx'],
include_dirs=[numpy.get_include()])],
cmdclass = {'build_ext':build_ext})
当然,仅当您在扩展模块中使用NumPy数组时才添加NumPy include目录(这是我们假定您使用Cython的条件)。NumPy中的distutils扩展还支持自动生成扩展模块并从.pyx
文件链接它。它的工作原理是,如果用户未安装Cython,则它将查找具有相同文件名但.c
扩展名的文件,然后使用该扩展名而不是.c
再次尝试生成该文件。
如果仅使用Cython来编译标准的Python模块,那么您将获得一个C扩展模块,该模块的运行速度通常比等效的Python模块要快一些。通过使用cdef
关键字静态定义C变量,可以进一步提高速度。
让我们看一下我们之前看到的两个示例,看看如何使用Cython来实现它们。这些示例使用Cython 0.21.1编译到扩展模块中。
这是一个名为Cython模块的一部分,该模块add.pyx
实现了我们之前使用f2py实现的复杂加法功能:
cimport cython
cimport numpy as np
import numpy as np
# We need to initialize NumPy.
np.import_array()
#@cython.boundscheck(False)
def zadd(in1, in2):
cdef double complex[:] a = in1.ravel()
cdef double complex[:] b = in2.ravel()
out = np.empty(a.shape[0], np.complex64)
cdef double complex[:] c = out.ravel()
for i in range(c.shape[0]):
c[i].real = a[i].real + b[i].real
c[i].imag = a[i].imag + b[i].imag
return out
此模块显示了如何使用该cimport
语句从numpy.pxd
Cython附带的标头中加载定义。看起来NumPy被导入了两次;cimport
仅使NumPy C-API可用,而常规规则则import
在运行时导致Python样式的导入,并可以调用熟悉的NumPy Python API。
该示例还演示了Cython的“类型化内存视图”,类似于C级别的NumPy数组,从某种意义上来说,它们是形状和跨步的数组,它们知道自己的范围(不同于通过裸指针寻址的C数组)。语法表示具有任意跨度的双精度的一维数组(向量)。一个连续的int数组将是,而一个浮点数矩阵将是。double complex[:]
int[::1]
float[:, :]
显示为cython.boundscheck
带有注释的装饰器,该装饰器根据功能打开或关闭内存视图访问的边界检查。我们可以使用它来进一步加快代码速度,但要牺牲安全性(或在进入循环之前进行手动检查)。
除了视图语法外,该函数可立即由Python程序员读取。变量的静态类型i
是隐式的。除了视图语法,我们还可以使用Cython的特殊NumPy数组语法,但最好使用视图语法。
我们使用Fortran创建的二维示例在Cython中也很容易编写:
cimport numpy as np
import numpy as np
np.import_array()
def filter(img):
cdef double[:, :] a = np.asarray(img, dtype=np.double)
out = np.zeros(img.shape, dtype=np.double)
cdef double[:, ::1] b = out
cdef np.npy_intp i, j
for i in range(1, a.shape[0] - 1):
for j in range(1, a.shape[1] - 1):
b[i, j] = (a[i, j]
+ .5 * ( a[i-1, j] + a[i+1, j]
+ a[i, j-1] + a[i, j+1])
+ .25 * ( a[i-1, j-1] + a[i-1, j+1]
+ a[i+1, j-1] + a[i+1, j+1]))
return out
该二维平均滤波器运行很快,因为循环位于C中,并且指针计算仅在需要时执行。如果上面的代码被编译为module image
,则img
可以使用以下代码非常快速地过滤出2-d图像:
import image
out = image.filter(img)
关于代码,有两点需要注意:首先,不可能将内存视图返回给Python。而是out
先创建一个NumPy数组,然后使用b
该数组的视图进行计算。其次,b
键入视图。这意味着二维数组具有连续的行,即C矩阵顺序。明确指定顺序可以加快某些算法的速度,因为它们可以跳过跨步计算。double[:, ::1]
Cython是几种科学的Python库(包括Scipy,Pandas,SAGE,scikit-image和scikit-learn)以及XML处理库LXML的首选扩展机制。语言和编译器维护良好。
使用Cython有几个缺点:
在编写自定义算法时,有时在包装现有的C库时,需要对C有一定的了解。特别是,在使用C内存管理(malloc
和朋友)时,很容易引入内存泄漏。但是,只需编译一个重命名为的Python模块就.pyx
可以加快速度,而添加一些类型声明可以在某些代码中显着提高速度。
很容易在Python和C之间失去清晰的分隔,这使得将C代码重用于其他与Python不相关的项目变得更加困难。
Cython生成的C代码很难读取和修改(通常会在烦人但无害的警告下进行编译)。
Cython生成的扩展模块的一大优势是它们易于分发。总而言之,Cython是一种非常强大的工具,可用于粘合C代码或快速生成扩展模块,因此不应忽视。这对于不能或不会编写C或Fortran代码的人特别有用。
C类型 是stdlib中包含的Python扩展模块,它允许您直接从Python调用共享库中的任意函数。这种方法使您可以直接从Python与C代码进行交互。这为Python提供了许多库供使用。但是,缺点是编码错误很容易导致丑陋的程序崩溃(就像C语言中那样),因为对参数进行的类型或界限检查很少。当数组数据作为指向原始内存位置的指针传入时,尤其如此。这样,您就有责任使子例程不会访问实际阵列区域之外的内存。但,
因为ctypes方法将原始接口公开给已编译的代码,所以它并不总是能够容忍用户的错误。强大地使用ctypes模块通常会涉及额外的Python代码层,以便检查传递给基础子例程的对象的数据类型和数组范围。这一额外的检查层(更不用说从ctypes对象到ctypes本身执行的C数据类型的转换)将使该接口比手写扩展模块接口慢。但是,如果被调用的C例程正在做大量工作,则此开销可以忽略不计。如果您是一名优秀的Python程序员,但C语言技能较弱,那么ctypes是将有用的接口写入(共享的)已编译代码库的简便方法。
要使用ctypes,您必须
有一个共享库。
加载共享库。
将python对象转换为ctypes可理解的参数。
使用ctypes参数从库中调用函数。
Python int / long,字符串和unicode对象会根据需要自动转换为等效的ctypes参数。None对象也会自动转换为NULL指针。所有其他Python对象必须转换为特定于ctypes的类型。围绕此限制有两种方法可以使ctype与其他对象集成。
不要设置函数对象的argtypes属性,也不要_as_parameter_
为要传递的对象定义
方法。该
_as_parameter_
方法必须返回一个Python int,它将直接传递给函数。
将argtypes属性设置为一个列表,该列表的条目包含带有名为from_param的类方法的对象,该方法知道如何将您的对象转换为ctypes可以理解的对象(int / long,字符串,unicode或具有该_as_parameter_
属性的对象)。
NumPy会优先使用第二种方法和第二种方法,因为它可以更安全。ndarray的ctypes属性返回一个对象,该对象的_as_parameter_
属性返回一个整数,该整数表示与之关联的ndarray的地址。结果,可以将此ctypes属性对象直接传递给期望指向ndarray中数据指针的函数。调用者必须确保ndarray对象的类型,形状正确,并且设置了正确的标志,否则,如果传入不适当数组的数据指针,则可能导致崩溃。
为了实现第二种方法,NumPy ndpointer
在numpy.ctypeslib
模块中提供了class-factory函数。这个classfactory函数产生一个合适的类,该类可以放在ctypes函数的argtypes属性条目中。该类将包含from_param方法,ctypes将使用该方法将传递给函数的任何ndarray转换为ctypes识别的对象。在此过程中,转换将检查用户在调用中指定的ndarray的任何属性ndpointer
。可以检查的ndarray的方面包括所传递的任何数组上的标志的数据类型,维数,形状和/或状态。from_param方法的返回值是数组的ctypes属性(因为它包含_as_parameter_
指向数组数据区域的属性)可以直接由ctypes使用。
ndarray的ctypes属性还具有其他属性,当将有关数组的其他信息传递给ctypes函数时,这些属性可能很方便。属性data,
shape和stride可以提供与数组的数据区域,形状和步幅相对应的ctype兼容类型。data属性返回一个c_void_p
代表指向数据区域的指针的。shape和stride属性每个都返回一个ctypes整数数组(如果为0-d数组,则返回None表示NULL指针)。数组的基本ctype是与平台上的指针大小相同的ctype整数。还有一些方法
data_as({ctype})
,和shape_as(<base ctype>)
strides_as(<base
ctype>)
。这些返回数据作为您选择的ctype对象,并使用您选择的基础基本类型返回shape / strides数组。为了方便起见,该ctypeslib
模块还包含c_intp
一个ctypes整数数据类型,其大小c_void_p
与平台上的大小相同
(如果未安装ctypes,则其值为None)。
该函数可以作为已加载共享库的属性或项目来访问。因此,如果./mylib.so
有一个名为的函数
cool_function1
,我可以通过以下方式访问此函数:
lib = numpy.ctypeslib.load_library('mylib','.')
func1 = lib.cool_function1 # or equivalently
func1 = lib['cool_function1']
在ctypes中,默认情况下,函数的返回值设置为'int'。可以通过设置函数的restype属性来更改此行为。如果该函数没有返回值(“ void”),请对None类型使用None:
func1.restype = None
如前所述,您还可以设置函数的argtypes属性,以使ctypes在调用函数时检查输入参数的类型。使用ndpointer
工厂函数来生成现成的类,以检查新函数的数据类型,形状和标志。该ndpointer
功能具有签名
ndpointer
(dtype = None,ndim = None,shape = None,flags = None )¶None
不检查具有该值的关键字参数。指定关键字可在转换为ctypes兼容对象时强制检查ndarray的这一方面。dtype关键字可以是任何理解为数据类型对象的对象。ndim关键字应为整数,shape关键字应为整数或整数序列。flags关键字指定传入的任何数组上所需的最小标志。可以将其指定为以逗号分隔的要求字符串,指示要求位或的“整数”或从变量的flags属性返回的flags对象。具有必要要求的阵列。
在argtypes方法中使用ndpointer类可以使使用ctypes和ndarray的数据区域调用C函数的安全性大大提高。您可能仍然希望将函数包装在其他Python包装器中,以使其易于使用(隐藏一些明显的参数并使某些参数输出参数)。在此过程中,requires
NumPy中的函数对于从给定输入返回正确类型的数组可能很有用。
在此示例中,我将展示如何使用ctypes来实现以前使用其他方法实现的加法函数和过滤器函数。第一,它实现了算法的C代码所包含的功能zadd
,dadd
,sadd
,cadd
,和dfilter2d
。该zadd
函数是:
/* Add arrays of contiguous data */
typedef struct {double real; double imag;} cdouble;
typedef struct {float real; float imag;} cfloat;
void zadd(cdouble *a, cdouble *b, cdouble *c, long n)
{
while (n--) {
c->real = a->real + b->real;
c->imag = a->imag + b->imag;
a++; b++; c++;
}
}
具有类似的代码cadd
,dadd
以及sadd
用于处理复杂的浮点,双精度和浮点数据类型,分别为:
void cadd(cfloat *a, cfloat *b, cfloat *c, long n)
{
while (n--) {
c->real = a->real + b->real;
c->imag = a->imag + b->imag;
a++; b++; c++;
}
}
void dadd(double *a, double *b, double *c, long n)
{
while (n--) {
*c++ = *a++ + *b++;
}
}
void sadd(float *a, float *b, float *c, long n)
{
while (n--) {
*c++ = *a++ + *b++;
}
}
该code.c
文件还包含以下功能dfilter2d
:
/*
* Assumes b is contiguous and has strides that are multiples of
* sizeof(double)
*/
void
dfilter2d(double *a, double *b, ssize_t *astrides, ssize_t *dims)
{
ssize_t i, j, M, N, S0, S1;
ssize_t r, c, rm1, rp1, cp1, cm1;
M = dims[0]; N = dims[1];
S0 = astrides[0]/sizeof(double);
S1 = astrides[1]/sizeof(double);
for (i = 1; i < M - 1; i++) {
r = i*S0;
rp1 = r + S0;
rm1 = r - S0;
for (j = 1; j < N - 1; j++) {
c = j*S1;
cp1 = j + S1;
cm1 = j - S1;
b[i*N + j] = a[r + c] +
(a[rp1 + c] + a[rm1 + c] +
a[r + cp1] + a[r + cm1])*0.5 +
(a[rp1 + cp1] + a[rp1 + cm1] +
a[rm1 + cp1] + a[rm1 + cp1])*0.25;
}
}
}
该代码相对于Fortran等效代码可能具有的优势是,它可以任意跨步(即非连续数组),并且根据编译器的优化能力,其运行速度也可能更快。但是,它显然比中的简单代码复杂filter.f
。此代码必须编译到共享库中。在我的Linux系统上,使用以下命令完成此操作:
gcc -o code.so -shared code.c
它将在当前目录中创建一个名为code.so的shared_library。在Windows上,不要忘记__declspec(dllexport)
在每个函数定义之前的行的void前面添加,或编写一个