[使用带有多维数组的cffi从python调用fortran dll

Question

我使用一个包含微分方程求解器的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)

由于具有一维数组的函数可以工作，所以我认为多维数组是个问题。

非常感谢任何帮助，锡洛

Answer 1

无法访问您所引用的dll会使给出明确的答案变得复杂，但是，dll的文档和所提供的Python脚本可能足以诊断问题。您的示例中至少有两个问题：

• C头接口：

  您的文档链接清楚地说明了函数的C头接口应为什么样。我不太熟悉C，Python的cffi或cdef，但是函数接口中a的参数声明似乎不对。如文档中所述，功能接口中的double** a（指向双精度指针的指针）最有可能是double a[]或double* a（双精度指针的指针）。

• 用Fortran顺序定义2d Numpy数组：

  1. 请注意，您的Numpy数组ax1..5是一维数组，因为这些数组只有一维order='F'，并且order='C'就内存布局和访问而言是等效的。因此，在此处指定order='F'可能不会达到预期的效果（对于多维数组，Fortran使用列主排序的方式）。
  2. [ax变量是Numpy数组的元组，而不是2d Numpy数组，因此在内存中的表示形式（将数据传递到Fortran dll时最重要）与2d数组非常不同。 >

寻求解决方案

我的第一步是更正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]]