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

Question

我明白在 C 中分配内存以确保二维数组是连续的，但是当我将它传递给 Fortran 时我仍然没有得到预期的结果。下面是我尝试的玩具版本：一个 main.c 文件，它为二维数组分配内存并为每个元素分配一个值，以及一个 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)

Answer 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 给出答案。

Answer 2

foo() 在 Fortran 端接收 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);

并将 print 语句插入到 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 ("Don’t use fake matrices") 以了解后一种方法的完整详细信息（感谢@Olaf 的建议）。

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

Attempting to pass a contiguous dynamic 2d array from C to Fortran

c

fortran

language-interoperability

multidimensional-array