Search code examples
c++fortranmixed-codeallocatable-array

Allocating memory in C for a Fortran allocatable


We are trying to take over the memory allocation of a legacy Fortran code (+100,000 lines of code) in C++, because we are using a C library for partitioning and allocating distributed memory on a cluster. The allocatable variables are defined in modules. When we call subroutines that use these modules the index seems to be wrong (shifted by one). However, if we pass the same argument to another subroutine we get what we expect. The following simple example illustrates the issue:

hello.f95:

 MODULE MYMOD
    IMPLICIT NONE
    INTEGER, ALLOCATABLE, DIMENSION(:) :: A
    SAVE
  END MODULE

  SUBROUTINE TEST(A)
    IMPLICIT NONE
    INTEGER A(*)
    PRINT *,"A(1): ",A(1)
    PRINT *,"A(2): ",A(2)
  END

  SUBROUTINE HELLO()
    USE MYMOD
    IMPLICIT NONE
    PRINT *,"A(1): ",A(1)
    PRINT *,"A(2): ",A(2)
    CALL TEST(A)
  end SUBROUTINE HELLO

main.cpp

extern "C" int* __mymod_MOD_a; // Name depends on compiler
extern "C" void hello_();      // Name depends on compiler

int main(int args, char** argv)
{
  __mymod_MOD_a = new int[10];
  for(int i=0; i<10; ++i) __mymod_MOD_a[i] = i;
  hello_();
  return 0;
}

We are compiling with:

gfortran -c hello.f95; c++ -c main.cpp; c++ main.o hello.o -o main -lgfortran;

Output from running ./main is

 A(1):            1
 A(2):            2
 A(1):            0
 A(2):            1

As you can see the output of A is different, though both subroutines printed A(1) and A(2). Thus, it seems that HELLO starts from A(0) and not A(1). This is probably due to that ALLOCATE has never been called directly in Fortran so that it is not aware of the bounds of A. Any work arounds?


Solution

  • Fortran array dummy arguments always start at the lower bound defined in the subroutine. Their lower bound is not retained during the call. Therefore the argument A in TEST() will always start at one. If you wish it to start from 42, you must do:

    INTEGER A(42:*)
    

    Regarding the allocation, you are playing with fire. It is much better to use Fortran pointers for this.

    integer, pointer :: A(:)
    

    You can then set the array to point to a C buffer by

    use iso_c_binding
    
    call c_f_pointer(c_ptr, a, [the dimensions of the array])
    

    where c_ptr is of type(c_ptr), interoperable with void *, which also comes from iso_c_binding.

    ---Edit--- Once I see that @Max la Cour Christensen implemented what I sketched above, I see I misunderstood the output of your code. The descriptor was indeed wrong, though I didn't write anything plain wrong. The solution above still applies.