Allocating memory in C for a Fortran allocatable

Question

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?

Answer 1

The ISO_C_BINDING "equivalent" code:

c++ code:

extern "C" int size;
extern "C" int* c_a;
extern "C" void hello();
int main(int args, char** argv)
{
  size = 10; 
  c_a = new int[size];
  for(int i=0; i<size; ++i) c_a[i] = i; 
  hello(); 
  return 0;
}

fortran code:

  MODULE MYMOD
    USE, INTRINSIC :: ISO_C_BINDING
    IMPLICIT NONE
    INTEGER, BIND(C) :: SIZE
    TYPE (C_PTR), BIND(C) :: C_A 
    INTEGER(C_INT), POINTER :: A(:)
    SAVE
  END MODULE

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

  SUBROUTINE HELLO() BIND(C)
    USE, INTRINSIC :: ISO_C_BINDING
    USE MYMOD
    IMPLICIT NONE
    CALL C_F_POINTER(C_A,A,(/SIZE/))
    PRINT *,"A(1): ",A(1)
    PRINT *,"A(2): ",A(2)
    CALL TEST(A)
  END SUBROUTINE

Output:

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

Answer 2

The internal representation of fortran arrays is very different than the one used in C/C++.

Fortran uses descriptors that start with a pointer to the array data, and followed by element type size, number of dimensions, some padding bytes, an internal 32/64 bit byte sequence indicating various flags such as pointer, target, allocatable, can be deallocated, etc. Most of these flags are not documented (at least in ifort that I have worked with), and at the end is a sequence of records, each describing the number of elements in the corresponding dimension, distance between elements, etc.

To 'see' an externally created array from fortran, you'd need to create such descriptors in C/C++, but, it does not end there because fortran also makes copies of them in the startup code of each subroutine before it gets to the first one of your statements, depending on indicators like 'in', 'out, 'inout', and other indicators used in the fortran array declaration.

Arrays within a type declared with specific sizes map well (again in ifort) to corresponding C struct members of the same type and number of elements, but pointer and allocatable type members are really descriptors in the type that need to be initialized to the correct values in all their fields so fortran can 'see' the allocatable value. This is at best tricky and dangerous, since the fortran compiler may generate copy code for arrays in undocumented ways for optimization purposes, but it needs to 'see' all the involved fortran code to do so. Anything coming outise of the fortran domain, is not known and can result in unexpected behavior.

Your best bet is to see if gfortran supports something like iso_c_binding and define such interfaces for your fortran code, and then use iso_c_binding intrinsics to map the C_PTR pointers to fortran pointers to types, arrays, etc.

You can also pass a pointer to a one-dimensional array of char, and its size, and this works for strings mostly as long as the size is passed by value as last argument (again, compiler and compiler-flag dependent).

Hope this helps.

EDIT: changed 'ifort's iso_c_binding' to 'iso_c_binding after Vladimir's comment - thanks!

Answer 3

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.