使用带多维数组的 cffi 从 python 调用 fortran dll
calling a fortran dll from python using cffi with multidimensional arrays
我使用的 dll 包含微分方程求解器以及其他有用的数学工具。不幸的是,这个 dll 是用 Fortran 语言编写的。我的程序是用 python 3.7 编写的,我使用 spyder 作为 IDE。
我成功地从 dll 中调用了简单的函数。但是,我似乎无法让需要多维数组的函数工作。
这是我尝试调用的函数的在线文档:
https://www.nag.co.uk/numeric/fl/nagdoc_fl26/html/f01/f01adf.html
如果我执行以下代码,内核会在没有错误消息的情况下死机:
import numpy as np
import cffi as cf
ffi=cf.FFI()
lib=ffi.dlopen("C:\Windows\SysWOW64\DLL20DDS")
ffi.cdef("""void F01ADF (const int *n, double** a, const int *lda, int *ifail);""")
#Integer
nx = 4
n = ffi.new('const int*', nx)
lda = nx + 1
lda = ffi.new('const int*', lda)
ifail = 0
ifail = ffi.new('int*', ifail)
#matrix to be inversed
ax1 = np.array([5,7,6,5],dtype = float, order = 'F')
ax2 = np.array([7,10,8,7],dtype = float, order = 'F')
ax3 = np.array([6,8,10,9],dtype = float, order = 'F')
ax4 = np.array([5,7,9,10], dtype = float, order = 'F')
ax5 = np.array([0,0,0,0], dtype = float, order = 'F')
ax = (ax1,ax2,ax3,ax4,ax5)
#Array
zx = np.zeros(nx, dtype = float, order = 'F')
a = ffi.cast("double** ", zx.__array_interface__['data'][0])
for i in range(lda[0]):
a[i] = ffi.cast("double* ", ax[i].__array_interface__['data'][0])
lib.F01ADF(n, a, lda, ifail)
由于函数与一维数组一起工作,我假设多维数组是问题所在。
非常感谢任何形式的帮助,
蒂洛
无法访问您引用的 dll 会使给出明确的答案变得复杂,但是,dll 的文档和提供的 Python 脚本可能足以诊断问题。您的示例中至少有两个问题:
C头接口:
您的文档 link 清楚地说明了函数的 C 头文件接口应该是什么样子。我不太精通C,Python的cffi或cdef,但是你的函数接口中a
的参数声明似乎是错误的。函数接口中的 double** a
(指向 double 的指针)最有可能是 double a[]
或 double* a
(指向 double 的指针),如文档中所述。
用 Fortran 顺序定义二维 Numpy 数组:
- 请注意,您的 Numpy 数组
ax1..5
是一维数组,因为数组只有一维 order='F'
和 order='C'
在内存布局和访问方面是等效的。因此,在此处指定 order='F'
可能不会产生预期的效果(Fortran 对多维数组使用列优先排序)。
- 变量
ax
是 Numpy 数组的元组,而不是 2d Numpy 数组,因此在内存中的表示形式非常不同(这在将数据传递给 Fortran dll 时至关重要)二维数组。
寻求解决方案
我的第一步是更正 C 头文件接口。接下来,我将使用 Fortran 排序将 ax
声明为具有二维的适当 Numpy 数组,然后将其转换为适当的数据类型,如本例所示:
#file: test.py
import numpy as np
import cffi as cf
ffi=cf.FFI()
lib=ffi.dlopen("./f01adf.dll")
ffi.cdef("""void f01adf_ (const int *n, double a[], const int *lda, int *ifail);""")
# integers
nx = 4
n = ffi.new('const int*', nx)
lda = nx + 1
lda = ffi.new('const int*', lda)
ifail = 0
ifail = ffi.new('int*', ifail)
# matrix to be inversed
ax = np.array([[5, 7, 6, 5],
[7, 10, 8, 7],
[6, 8, 10, 9],
[5, 7, 9, 10],
[0, 0, 0, 0]], dtype=float, order='F')
# operation on matrix using dll
print("BEFORE:")
print(ax.astype(int))
a = ffi.cast("double* ", ax.__array_interface__['data'][0])
lib.f01adf_(n, a, lda, ifail)
print("\nAFTER:")
print(ax.astype(int))
出于测试目的,请考虑使用以下具有 the same interface as your actual dll 的 Fortran 子例程来替代您的 dll。它只会将 10**(i-1) 添加到输入数组 a
的第 i 列。这将允许检查 Python 和 Fortran 之间的接口是否按预期工作,以及数组 a
的预期元素是否在以下操作:
!file: f01adf.f90
Subroutine f01adf(n, a, lda, ifail)
Integer, Intent (In) :: n, lda
Integer, Intent (Inout) :: ifail
Real(Kind(1.d0)), Intent (Inout) :: a(lda,*)
Integer :: i
print *, "Fortran DLL says: Hello world!"
If ((n < 1) .or. (lda < n+1)) Then
! Input variables not conforming to requirements
ifail = 2
Else
! Input variables acceptable
ifail = 0
! add 10**(i-1) to the i'th column of 2d array 'a'
Do i = 1, n
a(:, i) = a(:, i) + 10**(i-1)
End Do
End If
End Subroutine
编译 Fortran 代码,然后 运行 建议的 Python 脚本,得到以下输出:
> gfortran -O3 -shared -fPIC -fcheck=all -Wall -Wextra -std=f2008 -o f01adf.dll f01adf.f90
> python test.py
BEFORE:
[[ 5 7 6 5]
[ 7 10 8 7]
[ 6 8 10 9]
[ 5 7 9 10]
[ 0 0 0 0]]
Fortran DLL says: Hello world!
AFTER:
[[ 6 17 106 1005]
[ 8 20 108 1007]
[ 7 18 110 1009]
[ 6 17 109 1010]
[ 1 10 100 1000]]
我使用的 dll 包含微分方程求解器以及其他有用的数学工具。不幸的是,这个 dll 是用 Fortran 语言编写的。我的程序是用 python 3.7 编写的,我使用 spyder 作为 IDE。
我成功地从 dll 中调用了简单的函数。但是,我似乎无法让需要多维数组的函数工作。
这是我尝试调用的函数的在线文档: https://www.nag.co.uk/numeric/fl/nagdoc_fl26/html/f01/f01adf.html
如果我执行以下代码,内核会在没有错误消息的情况下死机:
import numpy as np
import cffi as cf
ffi=cf.FFI()
lib=ffi.dlopen("C:\Windows\SysWOW64\DLL20DDS")
ffi.cdef("""void F01ADF (const int *n, double** a, const int *lda, int *ifail);""")
#Integer
nx = 4
n = ffi.new('const int*', nx)
lda = nx + 1
lda = ffi.new('const int*', lda)
ifail = 0
ifail = ffi.new('int*', ifail)
#matrix to be inversed
ax1 = np.array([5,7,6,5],dtype = float, order = 'F')
ax2 = np.array([7,10,8,7],dtype = float, order = 'F')
ax3 = np.array([6,8,10,9],dtype = float, order = 'F')
ax4 = np.array([5,7,9,10], dtype = float, order = 'F')
ax5 = np.array([0,0,0,0], dtype = float, order = 'F')
ax = (ax1,ax2,ax3,ax4,ax5)
#Array
zx = np.zeros(nx, dtype = float, order = 'F')
a = ffi.cast("double** ", zx.__array_interface__['data'][0])
for i in range(lda[0]):
a[i] = ffi.cast("double* ", ax[i].__array_interface__['data'][0])
lib.F01ADF(n, a, lda, ifail)
由于函数与一维数组一起工作,我假设多维数组是问题所在。
非常感谢任何形式的帮助, 蒂洛
无法访问您引用的 dll 会使给出明确的答案变得复杂,但是,dll 的文档和提供的 Python 脚本可能足以诊断问题。您的示例中至少有两个问题:
C头接口:
您的文档 link 清楚地说明了函数的 C 头文件接口应该是什么样子。我不太精通C,Python的cffi或cdef,但是你的函数接口中
a
的参数声明似乎是错误的。函数接口中的double** a
(指向 double 的指针)最有可能是double a[]
或double* a
(指向 double 的指针),如文档中所述。用 Fortran 顺序定义二维 Numpy 数组:
- 请注意,您的 Numpy 数组
ax1..5
是一维数组,因为数组只有一维order='F'
和order='C'
在内存布局和访问方面是等效的。因此,在此处指定order='F'
可能不会产生预期的效果(Fortran 对多维数组使用列优先排序)。 - 变量
ax
是 Numpy 数组的元组,而不是 2d Numpy 数组,因此在内存中的表示形式非常不同(这在将数据传递给 Fortran dll 时至关重要)二维数组。
- 请注意,您的 Numpy 数组
寻求解决方案
我的第一步是更正 C 头文件接口。接下来,我将使用 Fortran 排序将 ax
声明为具有二维的适当 Numpy 数组,然后将其转换为适当的数据类型,如本例所示:
#file: test.py
import numpy as np
import cffi as cf
ffi=cf.FFI()
lib=ffi.dlopen("./f01adf.dll")
ffi.cdef("""void f01adf_ (const int *n, double a[], const int *lda, int *ifail);""")
# integers
nx = 4
n = ffi.new('const int*', nx)
lda = nx + 1
lda = ffi.new('const int*', lda)
ifail = 0
ifail = ffi.new('int*', ifail)
# matrix to be inversed
ax = np.array([[5, 7, 6, 5],
[7, 10, 8, 7],
[6, 8, 10, 9],
[5, 7, 9, 10],
[0, 0, 0, 0]], dtype=float, order='F')
# operation on matrix using dll
print("BEFORE:")
print(ax.astype(int))
a = ffi.cast("double* ", ax.__array_interface__['data'][0])
lib.f01adf_(n, a, lda, ifail)
print("\nAFTER:")
print(ax.astype(int))
出于测试目的,请考虑使用以下具有 the same interface as your actual dll 的 Fortran 子例程来替代您的 dll。它只会将 10**(i-1) 添加到输入数组 a
的第 i 列。这将允许检查 Python 和 Fortran 之间的接口是否按预期工作,以及数组 a
的预期元素是否在以下操作:
!file: f01adf.f90
Subroutine f01adf(n, a, lda, ifail)
Integer, Intent (In) :: n, lda
Integer, Intent (Inout) :: ifail
Real(Kind(1.d0)), Intent (Inout) :: a(lda,*)
Integer :: i
print *, "Fortran DLL says: Hello world!"
If ((n < 1) .or. (lda < n+1)) Then
! Input variables not conforming to requirements
ifail = 2
Else
! Input variables acceptable
ifail = 0
! add 10**(i-1) to the i'th column of 2d array 'a'
Do i = 1, n
a(:, i) = a(:, i) + 10**(i-1)
End Do
End If
End Subroutine
编译 Fortran 代码,然后 运行 建议的 Python 脚本,得到以下输出:
> gfortran -O3 -shared -fPIC -fcheck=all -Wall -Wextra -std=f2008 -o f01adf.dll f01adf.f90
> python test.py
BEFORE:
[[ 5 7 6 5]
[ 7 10 8 7]
[ 6 8 10 9]
[ 5 7 9 10]
[ 0 0 0 0]]
Fortran DLL says: Hello world!
AFTER:
[[ 6 17 106 1005]
[ 8 20 108 1007]
[ 7 18 110 1009]
[ 6 17 109 1010]
[ 1 10 100 1000]]