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

Question

I understand that extra care needs to be taken when allocating memory in C to ensure a 2d array is contiguous, but I still don't get expected outcomes when I pass it to Fortran. Below is a toy version of my attempt: A main.c file which allocates memory for a 2d array and assigns a value to each element, and a foo.f90 file that prints out elements of the 2d array.

#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

As output, I expect a list of array elements. Instead, I get the following output

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

Answer 1

foo() on the Fortran side receives x as an explicit-shape array, so we need to pass the address of the first element of x. So, we modify the prototype as

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

and pass the address of the first element as

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

Also, we probably need to use integer(c_int) rather than integer(c_long) to match int on the C side, such that

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

(On my computer, use of integer(c_long) gives segmentation fault because rows and cols are not passed correctly.)

Further, to make it easy to check the correspondence of array elements, I modified the test values as

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

and inserted print statements into the Fortran code as

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

Then, GCC-6 on my computer (OSX 10.9) gives

 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

As a side note, the following also seems to work (rather than creating x manually):

_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 );

but I'm not very sure whether this use of Array is standard-conforming... (in C).

---

[ EDIT ]

Using modern C, it seems possible to declare x directly as a rectangular array with dynamic size with memory allocated on the heap, such that:

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

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

(Please note that the variable names "cols" and "rows" refer to those on the Fortran side, so they appear counter-intuitive on the C side.)

With this x, we can proceed the same way as above to call foo_(), i.e.,

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

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

Please see the Jens Gustedts' blog ("Don’t use fake matrices") for full details of the latter approach (and thanks to @Olaf for suggestions).

Answer 2

You said so yourself: contiguous!

You do not allocate a contiguous array. To allocate one you must write:

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

Unfortunately, you must now write the following in C to index 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;

This assumes the matrix is in row-major order (row after row after row laid out contiguous in memory).

Note: in C99 there are Variable Length Arrays where the compiler manages x[i][j] properly. I do not have C99, but maybe another user can give an answer with VLAs.