Wrapping Fortran Derived Types : Part 1

May 23, 2019 7 min read

Abstract

Fortran C/C++/Python interoperability is an important consideration while designing scientific HPC codes. Fortran to a large extent still dominates in terms of lines of code in scientific software due to its legacy, reliable mathematical libraries and many years of research into the design of its compilers. Therefore it becomes mandatory to support an array of codes and libraries developed in mixed languages say Fortran, C and C++. It is counter productive to rewrite the codes and revalidate just for the sake of using a single monolithic language.

Some shortcomings of `iso_c_binding`

The Fortran 2003 standard introduced the iso_c_binding module to formalise known C interoperability hacks into Fortran language. It still leaves out a dark area, namely wrapping user defined types when the type has allocatable or pointer attributes. A simple example is shown below:

1type :: my_type
2  real(c_double), pointer :: my_data(:)
3end type my_type

Note that one could also include an allocatable array inside the user defined type such as:

1type :: my_type
2  real(c_double), allocatable :: my_data(:)
3end type my_type

The user is advised against using this approach since this is not flexible enough to expose the data outside of Fortran because a pointer variable cannot be made to point to this target array my_data if it was declared allocatable. For the very reason the following piece of code is invalid:

1type :: my_type
2  real(c_double), allocatable, target :: my_data(:)
3end type my_type

When compiled in gfortran one gets following error:

     real(c_double), allocatable, target :: my_data(:)
                                   1
Error: Attribute at (1) is not allowed in a TYPE definition

Another approach one would be tempted to take is to use the bind(C) technique from iso_c_binding as shown below:

1type, bind(C) :: my_type
2  real(c_double), allocatable, target :: my_data(:)
3end type my_type

But this results in a similar compiler error as the target approach:

   type, bind(C) :: my_type
                          2
     real(c_double), pointer :: my_data(:)
                                         1
Error: Component 'my_data' at (1) cannot have the POINTER attribute because it is a member of the BIND(C) derived type 'my_type' at (2)

Technique 1: Opaque object handle hack

A famous trick employed by programmers to wrap user defined Fortran types into C is the use of scratch space handles or opaque objects. Simply put the C side allocates scratch space in memory (adequate enough to store the Fortran type) and pass that as a handle to Fortran. This handle can then be safely stored in the C side and passed around to Fortran functions to do computation. This has been give in detail in the work of Pletzer et. al., “Exposing Fortran derived types to C and other languages”, comp. sci. eng. 2008. The object storing the scratch memory is called the opaque object to the Fortran type.

The opaque object technique applied to our simple problem is shown below:

 1! @file test.f90
 2! Fortran side code
 3module typedef
 4  use iso_c_binding
 5  type :: my_type
 6    real(c_double), pointer :: my_data(:)
 7  end type my_type
 8end module typedef
 9
10subroutine c_opaque_alloc(c_obj, n)
11  use typedef
12  implicit none
13  integer(c_int), intent(in) :: n
14  type(my_type) :: c_obj
15  allocate(c_obj%my_data(n))
16end subroutine c_opaque_alloc
17
18subroutine c_opaque_free(c_obj)
19  use typedef
20  implicit none
21  type(my_type) :: c_obj
22  deallocate(c_obj%my_data)
23end subroutine c_opaque_free

In the C side the following code is used to allocate the user defined type:

 1/** @file test.c
 2  * C side code
 3  */
 4#define FC_GLOBAL_(name,NAME) name##_
 5#define OPAQUE_ALLOC FC_GLOBAL_(c_opaque_alloc, C_OPAQUE_ALLOC)
 6#define OPAQUE_FREE FC_GLOBAL_(c_opaque_free, C_OPAQUE_FREE)
 7
 8#define OPAQUE_STORAGE_SIZE 64
 9
10void OPAQUE_ALLOC(char *c_obj, int *n);
11void OPAQUE_FREE(char *c_obj);
12
13int main(int nargs, char *args[] ) {
14  char c_obj[OPAQUE_STORAGE_SIZE];
15  int n = 100;
16  OPAQUE_ALLOC(c_obj, &n);
17  OPAQUE_FREE(c_obj);
18  return 0;
19}

The value of the OPAQUE_STORAGE_SIZE is an ambiguous one and varies from one compiler to another. Pletzer et. al. gives a long list of pointer array sizes used by various compilers based on their supported processor architecture. But the author does not give a reliable method to determine this size using existing Fortran features. This is one of the main contributions of this work where we provide a reliable method to obtain this size. Fortran 2008 standard defines a new inquiry function called STORAGE_SIZE. As per the standards STORAGE_SIZE(A, [KIND]) returns a scalar integer with the kind type parameter specified by KIND (or default integer type if KIND is missing). The result value is the size expressed in bits for an element of an array that has the dynamic type and type parameters of A. Note that a dummy type dummy_type object was created in order to obtain the internal storage size in bits of this data structure. The routine would throw an error if one passed simple the my_type. Currently there does not seem to be a better way to avoid this dummy type creation.

1subroutine c_storage_size( my_size )
2  use typedef
3  implicit none
4  type(my_type) :: dummy_obj
5  integer(c_int) :: my_size
6  my_size = STORAGE_SIZE(dummy_obj, c_int)
7  my_size = my_size / 8 ! return size in bytes
8end subroutine c_storage_size

 1#define FC_GLOBAL_(name,NAME) name##_
 2#define OPAQUE_ALLOC FC_GLOBAL_(c_opaque_alloc, C_OPAQUE_ALLOC)
 3#define OPAQUE_FREE FC_GLOBAL_(c_opaque_free, C_OPAQUE_FREE)
 4#define OPAQUE_SIZE FC_GLOBAL_(c_storage_size, C_STORAGE_SIZE)
 5
 6void OPAQUE_ALLOC(char *c_obj, int *n);
 7void OPAQUE_FREE(char *c_obj);
 8void OPAQUE_SIZE(int *size);
 9
10int main(int nargs, char *args[] ) {
11  char *c_obj;
12  int my_size;
13  OPAQUE_SIZE(&my_size);
14  c_obj = (char *)malloc(my_size);
15  int n = 100;
16  OPAQUE_ALLOC(c_obj, &n);
17  OPAQUE_FREE(c_obj);
18  free(c_obj);
19  return 0;
20}

In addition, the contents of the types can be indirectly exposed to C using Fortran call back functions (Getter/Setter type). While such an approach can be good for scalar values it is inefficient for array types because two copies of the arrays must exist in memory (one in the C side and the other in the Fortran side). In addition, the data copy/move from one to the other is redundant. An example of this approach is provided nevertheless:

 1subroutine c_opaque_array_copy(c_obj, n, vec)
 2  use typedef
 3  implicit none
 4  type(my_type) :: c_obj
 5  integer(c_int), intent(in) :: n
 6  real(c_double), intent(in) :: vec(n)
 7  allocate(c_obj%my_data(n))
 8  c_obj%my_data(1:n) = vec(1:n)
 9  write(*,*) c_obj%my_data
10end subroutine c_opaque_array_copy

 1#include <stdlib.h>
 2
 3#define FC_GLOBAL_(name,NAME) name##_
 4#define OPAQUE_ARRAY_COPY FC_GLOBAL_(c_opaque_array_copy, C_OPAQUE_ARRAY_COPY)
 5#define OPAQUE_FREE FC_GLOBAL_(c_opaque_free, C_OPAQUE_FREE)
 6#define OPAQUE_SIZE FC_GLOBAL_(c_storage_size, C_STORAGE_SIZE)
 7#define SIZE 100
 8
 9void OPAQUE_FREE(char *c_obj);
10void OPAQUE_SIZE(int *size);
11void OPAQUE_ARRAY_COPY(char *c_obj, int *n, double *vec);
12
13int main(int nargs, char *args[] ) {
14  char *c_obj;
15  int my_size, i, n=SIZE;
16  double vec[SIZE];
17
18  for(i=0; i<SIZE; ++i) vec[i] = -1.0;
19  OPAQUE_SIZE(&my_size);
20  c_obj = (char *)malloc(my_size);
21  OPAQUE_ARRAY_COPY(c_obj, &n, vec);
22  OPAQUE_FREE(c_obj);
23  free(c_obj);
24  return 0;
25}

One can improve upon this recipe and also return handles to the arrays/scalars inside the opaque handle to avoid redundant storage. C can then use these array handles to access the array data without any redundant copy/move.

1subroutine c_opaque_array(c_obj, c_array)
2  use typedef
3  implicit none
4  type(my_type) :: c_obj
5  type(c_ptr), intent(out) :: c_array
6  c_array = c_loc(c_obj%my_data(1))
7end subroutine c_opaque_array

 1#define FC_GLOBAL_(name,NAME) name##_
 2#define OPAQUE_ALLOC FC_GLOBAL_(c_opaque_alloc, C_OPAQUE_ALLOC)
 3#define OPAQUE_ARRAY FC_GLOBAL_(c_opaque_array, C_OPAQUE_ARRAY)
 4#define OPAQUE_FREE FC_GLOBAL_(c_opaque_free, C_OPAQUE_FREE)
 5#define OPAQUE_SIZE FC_GLOBAL_(c_storage_size, C_STORAGE_SIZE)
 6#define SIZE 100
 7
 8void OPAQUE_ALLOC(char *c_obj, int *n);
 9void OPAQUE_FREE(char *c_obj);
10void OPAQUE_SIZE(int *size);
11void OPAQUE_ARRAY(char *c_obj, double **vec);
12
13int main(int nargs, char *args[] ) {
14  char *c_obj;
15  int my_size, i, n=SIZE;
16  double *vec;
17
18  OPAQUE_SIZE(&my_size);
19  c_obj = (char *)malloc(my_size);
20  OPAQUE_ALLOC(c_obj, &n);
21  OPAQUE_ARRAY(c_obj, &vec);
22  for(i=0; i<SIZE; ++i)
23    vec[i] = -1.0;
24  OPAQUE_FREE(c_obj);
25  free(c_obj);
26  return 0;
27}

While the STORAGE_SIZE primitive gave us a portable and safe solution to wrap user defined types, some older Fortran compiler might not support the Fortran 2008 standards yet. In such situations one should resort to an alternative approach using TRANSFER which will be presented in the second part of this work.

Access to the source code with Makefiles to be uploaded soon to GitHub after testing against Intel, GCC, PGI, and LLVM compilers.

Wrapping Fortran Derived Types : Part 1

Abstract

Some shortcomings of iso_c_binding

Technique 1: Opaque object handle hack

Some shortcomings of `iso_c_binding`