Search code examples
cfortrandlsymfortran2003fortran-iso-c-binding

How to access (dynamically allocated) Fortran arrays in C


My main question is why arrays do such weird things and whether there is any way at all to do the following in a "clean" way.

I currently have a C program foo.c interfacing a Fortran program bar.f90 via dlopen/dlsym, roughly like in the code below:

foo.c:

#include <dlfcn.h>
#include <stdio.h>

int main()
{
int i, k = 4;
double arr[k];
char * e;

void * bar = dlopen("Code/Test/bar.so", RTLD_NOW | RTLD_LOCAL);

void (*allocArray)(int*);
*(void **)(&allocArray) = dlsym(bar, "__bar_MOD_allocarray");
void (*fillArray)(double*);
*(void **)(&fillArray) = dlsym(bar, "__bar_MOD_fillarray");
void (*printArray)(void);
*(void **)(&printArray) = dlsym(bar, "__bar_MOD_printarray");
double *a = (double*)dlsym(bar, "__bar_MOD_a");

for(i = 0; i < k; i++)
    arr[i] = i * 3.14;

(*allocArray)(&k);
(*fillArray)(arr);
(*printArray)();
for(i = 0; i < 4; i++)
    printf("%f ", a[i]);
printf("\n");

return 0;
}

bar.f90:

module bar

integer, parameter :: pa = selected_real_kind(15, 307)
real(pa), dimension(:), allocatable :: a
integer :: as

contains

subroutine allocArray(asize)
    integer, intent(in) :: asize

    as = asize
    allocate(a(asize))

    return
end subroutine

subroutine fillArray(values)
    real(pa), dimension(as), intent(in) :: values

    a = values
    return
end subroutine

subroutine printArray()
    write(*,*) a
    return
end subroutine

end module

Running main yields

0.0000000000000000        3.1400000000000001        6.2800000000000002        9.4199999999999999     
0.000000 -nan 0.000000 0.000000 

which shows that Fortran allocates the array correctly and even correctly stores the given values but they are not accessible via dlsym anymore (working on that data results in segfaults). I also tried this for fixed-size arrays - the results stay the same.

Does anyone know the reason for this behaviour? I, personally, would have expected things to either work out bidirectional or alternatively not at all - this "Fortran accepting C arrays but not vice versa" makes me wonder whether there's some basic mistake I made on accessing the array from C in this fashion.

The other (and even more important) question is, how to do array accesses like these "the right way". Currently I'm not even sure if sticking to the "Fortran as .so" interface is a good way at all - I think it would also be possible to attempt mixed programming in this case. Nonetheless, the arrays issue remains - I read that this could be solved somehow using the ISO C Binding, but I couldn't figure out how, yet (I haven't worked a lot with Fortran, yet, especially not with said Binding), so help on this issue would be greatly appreciated.

Edit:

Okay, so I read into the ISO C Binding a little more and found a quite useful approach here. Using C_LOC I can get C pointers to my Fortran structures. Unfortunately, pointers to arrays seem to be pointers to pointers and need to be dereferenced in the C code before they can be treates as C arrays - or something like that.

Edit:

Got my program to work now using C binding the way Vladimir F pointed out, at least for the most part. The C file and Fortran files are now linked together, so I can avoid the libdl interface, at least for the Fortran part - I still need to load a dynamic C library, get a function pointer to one of the symbols in there and pass that as a function pointer to Fortran, which later calls that function as part of its calculation. As said function expects double*s [arrays], I couldn't manage to pass my Fortran arrays using C_LOC, strangely - neither C_LOC(array) nor C_LOC(array(1)) passed the correct pointers back to the C function. array(1) did the trick though. Sadly, this isn't the "cleanest" way to do this. If anyone got a hint for me how to do this using the C_LOC function, that would be great. Nonetheless I accept Vladimir F's answer, as I deem it to be the safer solution.


Solution

  • In my opinion it is not good practice to try to access global data in Fortran library. It can be done using COMMON blocks, but they are evil and require statically sized arrays. Generally storage association is a bad thing.

    Never access the module symbols as "__bar_MOD_a" they are compiler specific and not meant to be used directly. Pass poiters using functions and subroutines.

    Pass the array as a subroutine argument. You can also allocate the array in C and pass it to Fortran. What can be also done is getting a pointer to the first element of the array. It will serve es the C pointer to the array.

    My solution, for simplicity without the .so, it is trivial to add it:

    bar.f90

    module bar
     use iso_C_binding
    
    implicit none
    
    integer, parameter :: pa = selected_real_kind(15, 307)
    
    real(pa), dimension(:), allocatable,target :: a
    integer :: as
    
    contains
    
    subroutine allocArray(asize,ptr) bind(C,name="allocArray")
        integer, intent(in) :: asize
        type(c_ptr),intent(out) :: ptr
    
        as = asize
        allocate(a(asize))
    
        ptr = c_loc(a(1))
    end subroutine
    
    subroutine fillArray(values) bind(C,name="fillArray")
        real(pa), dimension(as), intent(in) :: values
    
        a = values
    end subroutine
    
    subroutine printArray()  bind(C,name="printArray")
    
        write(*,*) a
    end subroutine
    
    end module
    

    main.c

    #include <dlfcn.h>
    #include <stdio.h>
    
    int main()
    {
    int i, k = 4;
    double arr[k];
    char * e;
    double *a;
    void allocArray(int*,double**);
    void fillArray(double*);
    void allocArray();
    
    
    for(i = 0; i < k; i++)
        arr[i] = i * 3.14;
    
    allocArray(&k,&a);
    fillArray(arr);
    printArray();
    for(i = 0; i < 4; i++)
        printf("%f ", a[i]);
    printf("\n");
    
    return 0;
    }
    

    compile and run:

    gcc -c -g main.c
    
    gfortran -c -g -fcheck=all bar.f90
    
    gfortran main.o bar.o
    
    ./a.out
    0.0000000000000000        3.1400000000000001        6.2800000000000002        9.4199999999999999     
    0.000000 3.140000 6.280000 9.420000 
    

    Note: There is no reason for the returns in your Fortran subroutines, they only obscure the code.