Calling a fortran function with a derived type containing assumed size arrays

Question

I would like to write functions or subprograms that return a derived type containing arrays with a size that is unknown at compile time. The rank is known. I would like to call these from C. One method on the fortran side is perhaps to use a parameterised derived type. Can I call this from C? Is there another preferred method? Can I make this a function rather than a subroutine? My fortran code so far is given below:

module dt

 use iso_c_binding
 implicit none
 
 type, bind(c) :: results(n)
    integer(c_int), len :: n
    real(c_double), dimension(n) :: myarray
 end type results

contains

 subroutine set_results(res, n)
   use iso_c_binding
   integer(c_int) :: n
   integer(c_int) :: i
   type(results(n)), intent(out) :: res
   do i=1,n
      res%myarray(i) = i
   end do
 end subroutine set_results
 
end module dt


program use_dt

 use dt
 implicit none

 type(results(10)) :: myresults

 call set_results(myresults, 10)
 print *,myresults
 
end program use_dt

Answer 1

Let's assume for the time being that you can allocate storage of the correct size in the calling program (the C++ desktop application). In other words all dimensions are known before entering the Fortran subprogram.

For interoperable types, the easiest option is to use assumed-size arrays, meaning array dimensions are need to be passed explicitly (or are assumed to be known constants):

module interop
use, intrinsic :: iso_c_binding, only: &
   rp => c_double, &
   ip => c_int
implicit none
contains
   !> Set results in an assumed-size array
   !> void set_results_v1(int n, double res[]);
   subroutine set_results_v1(n,res) bind(c)
      integer(ip), value :: n
      real(rp), intent(out) :: res(*)  ! <- contiguous storage
      integer(ip) :: i
      do i = 1, n
         res(i) = i
      end do
   end subroutine
end module

For a 2D array of size n × m, your procedure interface might look like this:

! void fill_random_2d(int n, int m, int lda, double *a)       // in C/C++
! void fill_random_2d(int n, int m, int lda, double a[][lda]) // in C99
   subroutine fill_random_2d(n,m,lda,a) bind(c)
      integer(ip), value :: n, m, lda
      real(rp), intent(inout) :: a(lda,*) !  
      call random_number(a(1:n,1:m))
   end subroutine

The integer lda can be used to account for any array padding. With 2D arrays you need to be careful with differences due to opposite storage order, i.e. Fortran (column-major) vs C/C++ (row-major). Arguably, it's easier to assume column-oriented storage in C and do multi-dimensional indexing manually if necessary.

---

The newer, but potentially more powerful option is to the use the C descriptors introduced in Fortran 2018. Descriptors can be used to create Fortran objects in C/C++ which have no direct equivalent. Descriptors allow you to interoperate also with assumed-shape, allocatable, and pointer arrays.

Using an assumed-shape array, your routine is written as

   !> Set results using assumed-shape array 
   !> (Fortran 2018 C-descriptor needed for interoperation)
   subroutine set_results_v2(res) bind(c)
      real(rp), intent(out) :: res(:) ! <- contiguous or strided storage
      integer :: i
      do i = 1, size(res)
         res(i) = i
      end do
   end subroutine

Here's how the calling the Fortran procedure will look like in C++:

// For GCC compile using
//  g++ -Wall main.cpp interop.o -lgfortran
// where interop.o contains the external routine.
#include <vector>
#include <iostream>

#include <ISO_Fortran_binding.h>

/**
 * Fortran routine prototypes  
 */ 
extern "C" {
void set_results_v2(CFI_cdesc_t *res);
}


int main(int argc, char const *argv[])
{
    const int N = 10;

    std::vector<double> myresult(N);
    
    CFI_CDESC_T(1) res_desc; /* Rank-1 array descriptor */

    // The C descriptor is an opaque struct that stores
    // information about a Fortran objects which have
    // no proper equivalent in C. Such structs can be 
    // used as dummy arguments in inter-language calls. 

    // The descriptor can only be modified by functions provided
    // in the header "ISO_Fortran_binding.h". These functions
    // expect a pointer to the array descriptor.

    auto res = reinterpret_cast<CFI_cdesc_t *>(&res_desc);

    // First we establish a local array of extents
    // In this case we are wrapping a std::vector, which
    // can be seen as a 1-dimensional array.

    CFI_index_t res_extents[1] = { N };

    // Now we establish the C-descriptor.
    // In this case `res` will be a rank-1 array
    // of type double referring to the data
    // in the std::vector `myresult`.

    int err = CFI_establish(
            res,                  /* pointer to the descriptor */
            myresult.data(),      /* pointer to the data */
            CFI_attribute_other,  /* attribute */
            CFI_type_double,      /* type of  */
            sizeof(double),       /* size of element */
            1,                    /* rank */
            res_extents);         /* array extents */

    // In case something went wrong, a non-zero
    // value is returned, which can be used for
    // error-checking

    if (err != CFI_SUCCESS) 
        return 1;

    // Now we have a complete descriptor which
    // can be passed to a Fortran routine.

    set_results_v2(res);

    for(auto v : myresult) {
        std::cout << v << ' ';
    }
    std::cout << '\n';

    return 0;
}

If you'd like to pass a 2-dimensional assumed-shape array you just need to adjust the rank-constant (replace 1 with 2), and adjust the extents, e.g.

    CFI_index_t res_extents[2] = { N, M };

It is also possible to pass array slices, provide custom array bounds, and more. To learn more about this approach I would recommend reading Metcalf, Reid, & Cohen (2018), Modern Fortran Explained, Oxford University Press, pgs. 397-421.

---

Upon reading the OP's comment, here is a take at a solution where the C side is responsible for allocating buffers encapsulated in a struct, and on the Fortran we create operate on a derived type which shadows the C struct:

(I haven't tested this, so apologies if any errors remain)

module work

use, intrinsic :: iso_c_binding, only: &
   c_int, c_double, c_ptr, c_f_pointer

implicit none
private

public :: do_work

!> The C storage struct (memory will be 
!> allocated and deallocated by calls to malloc/free)
! struct cwork {
!  int na, nb;
!  double *a, *b;
! };
type, bind(c) :: cwork
   integer(c_int) :: na, nb
   type(c_ptr) :: a, b
end type

!> A Fortran view of the cwork struct
type :: cwork_view
   integer(c_int) :: na, nb
   real(c_double), pointer :: a(:) => null()
   real(c_double), pointer :: b(:) => null()
end type

interface
   
   !> C helper function to perform memory allocation
   !> Returns 0 if successful, and a negative number 
   !> if something went wrong.

   ! int workspace_alloc(struct cwork *work, int na, int nb);
   integer(c_int) function workspace_alloc(work,na,nb) bind(C)
      import cwork, c_int
      type(cwork), intent(inout) :: work
      integer(c_int), value :: na, nb
   end function

end interface

contains

   !> Helper routine to create Fortran view of C struct
   subroutine make_view(fw,cw,istat) 
      type(cwork_view), intent(out) :: fw
      type(cwork), intent(in) :: cw

      fw%na = cw%na
      call c_f_pointer(cw%a,fw%a,[cw%na])
      
      fw%nb = cw%nb
      call c_f_pointer(cw%b,fw%b,[cw%nb])

   end subroutine

   !> Public routine to be called from C/C++
   !> void do_work(struct cwork *cw, int *ierr);
   subroutine do_work(cw,ierr) bind(c)
      type(cwork), intent(out), target :: cw
      integer(c_int), intent(out) :: ierr

      type(cwork_view) :: fw
      integer :: na, nb, i

      ! maybe na and nb become known at runtime based 
      ! on same other routine (either in C or Fortran)

      ! we call a C routine to allocate the needed memory
      ierr =  workspace_alloc(cw,na,nb)
      if (ierr < 0) then
         ! oops, something went wrong
         return
      end if

      ! if allocation succedded, we call the helper 
      ! routine to associate the Fortran pointers
      call make_view(fw=fw,cw=cw)

      ! -------------------------

      ! finally we can do our work

      do i = 1, fw%na
         fw%a(i) = i
      end do

      do i = fw%nb, 1, -1
         fw%b(i) = i**2
      end od

   end subroutine

end module

In your main C++ app, you would then call the procedure as follows:

struct cwork cw;
int ierr;

// memory allocated within do_work
do_work(&cw,&ierr);

if (ierr)
  printf("do_work failed");

for(int i = 0; i < cw.na; ++i) {
   printf("%d",cw.a[i]);
}

for(int i = 0; i < cw.nb; ++i) {
   printf("%d",cw.b[i]);
}

free(cw.a);
free(cw.b);

You will probably want to use the RAII idiom, and wrap the C struct in a C++ class which will do the deallocation once the struct goes out of scope.

Answer 2

A solution is to declare the arrays as type(c_ptr) in the Fortran derived type, and convert them to Fortran pointers with the c_f_pointer() routine. AFAIK the allocation has to be performed on the C side.

module dt
 use iso_c_binding
 implicit none
 
 type, bind(c) :: results
    integer(c_int) :: n
    type(c_ptr) :: myarray
 end type results

contains

 subroutine set_results(res) bind(c)
   use iso_c_binding
   integer :: i
   type(results), intent(out) :: res
   
   real(c_double), pointer :: myarrayf(:)
   
   call c_f_pointer(res%myarray,myarrayf,[res%n])
   do i=1,res%n
      myarrayf(i) = i
   end do
 end subroutine set_results
 
end module dt

The tedious part is that you have to call c_f_pointer() for each individual array in the derived type.