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
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.