试图将一个连续的动态二维数组从 C 传递到 Fortran

Posted 2023-03-10

【中文标题】试图将一个连续的动态二维数组从 C 传递到 Fortran

【英文标题】：Attempting to pass a contiguous dynamic 2d array from C to Fortran

【发布时间】：2016-12-09 11:58:32

【问题描述】：

我知道extra care needs to be taken 在 C 中分配内存以确保二维数组是连续的时，但是当我将它传递给 Fortran 时仍然没有得到预期的结果。下面是我尝试的玩具版本：为 2d 数组分配内存并为每个元素分配值的 main.c 文件，以及打印出 2d 数组元素的 foo.f90 文件。

#include <stdio.h>
#include <stdlib.h>

void foo_(double **,int *,int *);

int main() 
    int i, j, cols, rows;
    double *_x, **x;

    cols = 3;
    rows = 2;

//  Allocate memory
    _x = malloc(rows*cols*sizeof(double));
    x = malloc(cols*sizeof(double *));
    for(i=0;i<cols;i++)
        x[i] = &(_x[rows*i]);

//  Generate elements for the 2d array 
    for(i=0;i<cols;i++)
    for(j=0;j<rows;j++)
        x[i][j] = i*j;

//  Call Fortran subroutine foo
    foo_(x,&rows,&cols);

    return EXIT_SUCCESS;

foo.h

subroutine foo(x,rows,cols)
    use iso_c_binding
    implicit none

    integer(c_long), intent(in)                :: rows,cols
    real(c_double), intent(in), dimension(rows,cols) :: x

    integer :: i,j

    do i = 1,cols
    do j = 1,rows
        print *, j,i,x(j,i)
    end do
    end do

end subroutine

作为输出，我需要一个数组元素列表。相反，我得到以下输出

           1           1   1.1654415706619996E-316
           2           1   1.1654423611670330E-316
Segmentation fault (core dumped)

【问题讨论】：

您的 C 代码中有 no 二维数组，没有任何东西可以表示或指向一个！指针不是数组。如果您需要二维数组，请使用一个！

双指针可以表示二维数组，实际上有@Olaf，我认为问题出在 foo_(double **,int *,int *);不应该只是 foo_(double **,int,int); ?

@Olaf 你的评论要么是错误的，要么是愚蠢的，以至于在教学上毫无用处。在 C 编程语言中，存储在堆中的多维数组在运行时保留其内存，指针用作管理地址的簿记工具（cs.swarthmore.edu/~newhall/unixhelp/C_arrays.html）。

@francescalus 好点。我已编辑问题以包含输出。

@ThomasKelly：1) C 标准不强制使用特定的内存管理模型进行动态内存分配。 2）你的断言是完全错误的。锯齿状数组不是使用二维数组的正确方法。这是一个完全不同的数据结构！您可以在 C 中轻松分配 2D 数组。您可能想要更新到 17 岁的 C99，或者 - 更好 - 到 C11（也已经 5 岁）的 C 标准。告诉你使用 double ** 的人正在教古老的 C。3）这可能是你的问题的原因之一。

【参考方案1】：

你自己也是这么说的：连续！

您没有分配一个连续的数组。要分配一个，你必须写：

//  Allocate memory
double *x = malloc(rows*cols*sizeof(double));

不幸的是，您现在必须用 C 编写以下代码来索引 x：

//  Generate elements for the 2d array 
for(i=0;i<cols;i++)
    for(j=0;j<rows;j++)
        *(x+j*cols+i) = i*j;

这假设矩阵是按行优先的顺序（在内存中连续排列的一行接一行）。

注意：在 C99 中有可变长度数组，编译器在其中正确管理 x[i][j]。我没有 C99，但也许其他用户可以用 VLA 给出答案。

【讨论】：

这是类型错误。 x 是 double ** 并且您的索引是错误的。为什么不使用二维数组呢？这完全取决于 Fortran ABI 的预期。

问题标记为 C，即标准 C，即 C11。虽然 C11 使 VLA 成为可选项，但没有任何现代 C11 兼容编译器不支持它们（使现有语言功能成为可选项是委员会维护向下兼容性实践的重大突破，可能让过时编译器的维护者更好地出售它们）。跨度>

【参考方案2】：

Fortran 端的foo() 接收x 作为显式形状数组，因此我们需要传递x 的第一个元素的地址。所以，我们将原型修改为

void foo_( double *, int *, int * );

并将第一个元素的地址传递为

foo_( _x, &rows, &cols );
// or
// foo_( &( _x[0] ), &rows, &cols );
// or
// foo_( &( x[0][0] ), &rows, &cols );

另外，我们可能需要使用integer(c_int) 而不是integer(c_long) 来匹配C 端的int，这样

integer(c_int), intent(in) :: rows,cols

（在我的计算机上，使用 integer(c_long) 会导致分段错误，因为 rows 和 cols 未正确传递。）

另外，为了方便检查数组元素的对应关系，我将测试值修改为

for(i=0;i<cols;i++)
for(j=0;j<rows;j++)
    x[i][j] = (i+1) + 1000*(j+1);

并将打印语句插入到 Fortran 代码中

print *, "rows = ", rows
print *, "cols = ", cols

然后，我的计算机上的 GCC-6 (OSX 10.9) 给出了

 rows =            2
 cols =            3
           1           1   1001.0000000000000     
           2           1   2001.0000000000000     
           1           2   1002.0000000000000     
           2           2   2002.0000000000000     
           1           3   1003.0000000000000     
           2           3   2003.0000000000000

作为旁注，以下似乎也有效（而不是手动创建x）：

_x = malloc( rows * cols * sizeof(double) );

typedef double (*Array)[ rows ];
Array x = (Array) _x;

// set x[i][j] the same way

foo_( _x, &rows, &cols );
// or
// foo_( &( x[0][0] ), &rows, &cols );

但我不太确定 Array 的这种用法是否符合标准...（在 C 中）。

---

[编辑]

使用现代 C，似乎可以将 x 直接声明为具有动态大小的矩形数组，并在堆上分配内存，这样：

double (* x)[rows] = malloc( cols * rows * sizeof(double) );

// or
// double (* x)[rows] = malloc( sizeof( double[cols][rows] ) );

（请注意，变量名称“cols”和“rows”指的是 Fortran 端的变量名，因此它们在 C 端显得违反直觉。）

有了这个x，我们可以像上面一样继续调用foo_()，即，

// set x[i][j] via loops

foo_( &( x[0][0] ), &rows, &cols );

请参阅Jens Gustedts' blog（“不要使用假矩阵”）了解后一种方法的完整详细信息（感谢@Olaf 的建议）。

【讨论】：

这可能是真的 Fortran ABI 使用一维数组来连接 C 代码，但应该首先验证这一点。毕竟，这仍然与 OP 使用的锯齿状数组非常不同，也不是 2D 数组。由于它具有相同的（考虑正确的行/列顺序）布局，它应该可以工作，但它仍然具有不同的类型并且可能导致细微的错误，resp。并发症，需要其他不必要的演员..

OP 使用指向指针的指针，resp。一个指针数组。通常称为“锯齿状数组”，无论二级数组的大小是否相同。他在块中分配子数组的事实并没有改变这一点。仅 C 语言的正确方法是使用 2D 数组，但正如我所写，我不确定 Fortran 期望什么。

@Olaf 第一级数组 (_x) 是线性且连续的，对吗？然后，我的理解是x（=二级数组）提供了对_x指向的内存的矩形访问。传递给 Fortran 的是 malloc() 分配的地址，它是线性的（我假设）。另外，我的理解是“锯齿状数组”是每一行或每一列都是独立分配的，所以内存是不连续的。

但是因为 OP 使用变量名“rows”和“cols”有点令人困惑（为了匹配 Fortran 主要列并避免转置），所以乍一看代码对我来说不是很可读（对不起，只是我的意见......）

是的。 Jens 是这里的成员之一，通常可以信任 C 标准问题。也不知道他写博客。虽然 C99 现在已经 17 岁了，但这篇文章很有道理。