Iso_c_binding: pass an array from c to fortran (Edit: python interop content)

1

Hey everyone,

I’m having trouble figuring out how to pass an array from c to fortran.

Here’s an example of some C and fortran code that illustrates what I’m trying to do, but does not work:

extern void trivial(float*, int);
int main(void)
{
	float f[] = {1,2,3,4};
	float * fptr = &f[0];	

	trivial(fptr, 4);
}

module stuff
        use iso_c_binding
contains
        subroutine trivial(arr, n) bind(c)
                integer(c_int) :: n
                real(c_float)  :: arr(n)

                arr = arr+1
                print *, arr
        end subroutine
end module

It compiles alright, by compiling them separately to .o files and then linking them with gcc -lgfortran, but it segfaults when run. Does anyone know how to get this to work?

Thanks

2

change

to
integer(c_int), value :: n
and it should work

5 Likes
3

:man_facepalming:

Thank you very much. That fixed it

4

@hsnyder and any readers needing to interoperate Fortran codes with a companion C processor:

You may note many of the feature sets in Fortran starting with Fortran 90 thru’ Fortran 2018 can really prove useful in scientific and technical computing contexts, classic examples would be ASSUMED-SHAPE arrays, dummy arguments with ALLOCATABLE attribute, etc. With this in mind, note a “trivial” subprogram in modern Fortran can be as follows where the array descriptors are part of the dummy argument itself thereby making it unnecessary to pass additional arguments:

   ..
   subroutine Fsub( x ) bind(C, name="Fsub")
      ! Argument list
      real(kind=CF), intent(inout) :: x(:)  !<-- assumed-shape
..

Then if your strengths are in programming in C, C++, note with some additional effort, you can layer your C code to interoperate with Fortran’s strengths using the enhanced interoperability facility introduced in Fortran 2018 via “ISO_Fortran_binding.h” header. The standard for Fortran provides details on this, as does the book by Metcalf et al., “Modern Fortran Explained”

Here is a modified example of your code: as you will notice, it is an overkill in simple scenarios, but it can prove handy when it comes to consuming performant Fortran codes with a companion C processor.

#include <stdlib.h>
#include <stdio.h>
#include "ISO_Fortran_binding.h"

extern void Fsub(CFI_cdesc_t *);

int main(void)
{
   enum N { N = 4 };
   float x[N];
   int i;
   CFI_CDESC_T(1) dat; // employ a macro defined in Fortran binding header
   CFI_index_t ext[1]; // employ a type defined in Fortran binding header

   ext[0] = (CFI_index_t)N;
   i = CFI_establish((CFI_cdesc_t *)&dat, x, CFI_attribute_other,
                      CFI_type_float, 0, (CFI_rank_t)1, ext);

   Fsub((CFI_cdesc_t *)&dat );
   printf("In C main: following call to Fsub, x = \n");
   for (i=0; i<N-1; i++) printf("%f, ", x[i]); printf("%f\n", x[N-1]);

   return EXIT_SUCCESS;
}

module m
   use, intrinsic :: iso_c_binding, only : CF => c_float
contains
   subroutine Fsub( x ) bind(C, name="Fsub")
      ! Argument list
      real(kind=CF), intent(inout) :: x(:)
      integer :: i
      print *, "In Fsub: shape(x) = ", shape(x)
      x = [( real(i,kind=CF), i=lbound(x,dim=1),ubound(x,dim=1) )]
   end subroutine 
end module

C:\Temp>ifort /c /standard-semantics /warn:all /stand:f18 m.f90
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.1 Build 20201112_000000
Copyright (C) 1985-2020 Intel Corporation. All rights reserved.

C:\Temp>cl /c /EHsc /W3 c.c
Microsoft (R) C/C++ Optimizing Compiler Version 19.26.28806 for x64
Copyright (C) Microsoft Corporation. All rights reserved.

c.c

C:\Temp>link c.obj m.obj /subsystem:console /out:p.exe
Microsoft (R) Incremental Linker Version 14.26.28806.0
Copyright (C) Microsoft Corporation. All rights reserved.

C:\Temp>p.exe
In Fsub: shape(x) = 4
In C main: following call to Fsub, x =
1.000000, 2.000000, 3.000000, 4.000000

C:\Temp>

1 Like
5

I was wondering a few days ago on what were the motives to introduce the ISO_Fortran_binding.h which opens up the array API? Was this a feature requested by Fortran application developers, or was it ultimately pushed by compiler developers?

Do you know of any publically available efforts to build a higher level API into Fortran objects?

What I have in mind is something like:

// fortran_array.cpp

#include "ISO_Fortran_binding.h"
#include <iostream>

template<typename T>
struct FortranArray
{
    const CFI_cdesc_t *obj; // A const reference to rank 2 Fortran array (can be discontiguous with regular strides).

    FortranArray(const CFI_cdesc_t *obj_) : obj(obj_) {}

    // Returns the i'th component of the j'th column in the array:
    inline T operator()(const size_t i, const size_t j) const
    {
        char *elt = (char *) obj->base_addr;
        elt += (i*obj->dim[0].sm + j*obj->dim[1].sm);
        return *(T *) elt;
    }

    size_t size(const size_t d) {
        return obj->dim[d].extent;
    }

};

extern "C" {

void print_array(CFI_cdesc_t *array) {

    // View into a Fortran array
    FortranArray<double> a(array);

    for (size_t i = 0; i < a.size(0); i++) {
        for (size_t j = 0; j < a.size(1); j++) {
            std::cout << a(i,j) << " ";
        }
        std::cout << "\n";
    }
}

}

The driver program in Fortran:

! test.f90

module print_mod
  use, intrinsic :: iso_c_binding, only: c_double
  implicit none

  interface
    subroutine print_array(A) bind(c,name="print_array")
      import c_double
      real(c_double), intent(in) :: A(:,:)
    end subroutine
  end interface

end module

program test

  use print_mod
  implicit none

  real(c_double) :: A(3,2)

  A = reshape([1,2,3,4,5,6],shape(A))
  
  call print_array(A)

end program

Compiling, linking, and output:

$ g++ fortran_array.cpp -c
$ gfortran fortran_array.o test.f90 -o test -lstdc++
$ ./test
1 4 
2 5 
3 6

Edit: I added a template parameter to the FortranArray struct.

1 Like
6

My question above was addressed to @FortranFan.

7

Re: your first question about enhanced interoperability with C in Fortran 2018, the WG5 site can provide you with a document trail on this work post Fortran 2003, especially with what a focused development with a technical specification (TS) project 29113. You will notice there was motivation from the users, particularly the MPI community, to make it easier to access from C the Fortran APIs with many of the advantageous features of Fortran such as dummy arguments with assumed-shape arrays and other attributes such as ALLOCATABLE.

As to your other question, I’m not aware of any publicly available efforts to encapsulate Fortran standard types (and methods) in ISO_Fortran_binding.h.

8

@FortranFan’s explanation inspired me to play around a bit. I came up with the following, as an interesting way to convert Python arrays to Fortran arrays, using pybind11 (which is a header only library that makes it easy to call c++ from python).

I’m aware of the existence of f2py, which would eliminate the need for a c++ wrapper, but at work we use pybind11 already for wrapping C kernels to be called from python, and I’ve had mixed experiences with f2py, f90wrap, and such tools. Not saying they’re bad, but for me this worked better.

trivial.f90

module trivial_fortran
        use iso_c_binding
contains
        subroutine trivial(arr) bind(c)
                real(c_float), intent(inout)      :: arr(:)

                arr = arr + 1
        end subroutine
end module

pybind_fortran.h

#pragma once

#include<pybind11/pybind11.h>
#include<pybind11/numpy.h>

#include<stdexcept>
#include<string>

extern "C" {
#include<ISO_Fortran_binding.h>
}

#define CONCAT(a,b) a ## b
#define FARRAY_PTR     CFI_cdesc_t*
#define FARRAY         CFI_CDESC_T(CFI_MAX_RANK)
#define FARRAY_TYPE_T  CFI_type_t
#define FARRAY_TYPE(x) CONCAT(CFI_type_,x)


void py_array_to_fortran(pybind11::array & pyarr, FARRAY_PTR farr, FARRAY_TYPE_T type)
{
	int rc = CFI_establish((CFI_cdesc_t*)farr, 
                           pyarr.request().ptr, 
                           CFI_attribute_other, 
                           type, 
                           0, 
                           (CFI_rank_t)pyarr.request().shape.size(), 
                           (const CFI_index_t *)pyarr.request().shape.data());
	// error check omitted for brevity
}

triv_wrap.cpp

#include<pybind11/pybind11.h>
#include<pybind11/numpy.h>

#include "pybind_fortran.h"

namespace py = pybind11;


extern "C" void trivial(FARRAY_PTR);
void trivial_wrap(py::array_t<float, py::array::c_style> arr)
{

	FARRAY dat;
	py_array_to_fortran(arr, (FARRAY_PTR) &dat, FARRAY_TYPE(float));

	trivial((FARRAY_PTR) &dat);
}

PYBIND11_MODULE(triv_wrap, m) {
	m.def("trivial", &trivial_wrap);
}

test.py

import triv_wrap as tw
import numpy as n

a = n.arange(9, dtype=n.float32)
print(a)
tw.trivial(a)
print(a)

Makefile

triv_wrap.so:
	g++ -fpic -c $(shell python3 -m pybind11 --includes) triv_wrap.cpp 
	gfortran -fpic -c trivial.f90
	g++ -shared triv_wrap.o trivial.o -o triv_wrap.so -lgfortran

clean:
	rm *.o *.mod *.so -f

test:
	python3 test.py

.PHONY: clean test

Running make && make test:

python3 test.py
[0. 1. 2. 3. 4. 5. 6. 7. 8.]
[1. 2. 3. 4. 5. 6. 7. 8. 9.]
2 Likes
9

I’ve been trying to learn Fortran-C interoperation and I’ve come across a rather simple question/realization about arrays.

  1. Is it correct that the result of a bind(c) function cannot be an array? For example, the following is not possible?
function foo() bind(C)
  import :: c_int
  integer(c_int), dimension(10) :: foo
end function
  1. Accordingly, a function with the following prototype cannot be called from Fortran, correct?
char *concatenate_substrings(const char *substring1, const char *substring2);
  1. The only option is then to convert to a void function and make a subroutine interface?

Thanks.

10

I am not sure if function result can somehow be a pointer or not, but I know that if one passes an array pointer from c, it can then be accessed from fortran side. That is how I did it when communicating between Python-Fortran. (I had to use c_ptr type and intrinsic subroutine c_f_pointer on Fortran side)

Edit. Sorry, I realize now that I didn’t read your post well and I see that I didn’t answer to your question repeated what you already know.

11

@HugoMVale ,

No, that’s incorrect.

You may consume a function in C which returns a pointer to some memory that you may see as an array. Note you can refer to any number of books and blogs and other material to know the treatment of “arrays” in C is different than what you may envision with Fortran. Rather, you will notice there are quite a few similarities in how arrays are applied in C code with the semantics of pointers.

You can use such info to interoperate said C function with Fortran, but note that:

  1. you are interoperating effectively with an opaque address i.e., void * pointer. Thus presuppositions as to the underlying type of data in a given address comes into play that is in essence casting blindly. This can be unsafe.
  2. Also, on the C side, such a function has to either allocate memory dynamically and in which case the caller has to free such memory. Or, the function has to work with a predefined buffer. Again, this leads to further presuppositions.
  3. At a bare minimum for reasonably safe coding practice, the API toward the interoperation shall communicate a measure of the array size.

Here is a naive example:

Click to see code
  • C library function
#include <stdlib.h>
#include <stdio.h>

char* foo( size_t );

char* foo( size_t len ) {
   char *r = malloc(len);
   int irc = snprintf(r, len, "Hello World!");
   return r;
}
  • Fortran caller that presumes the need to “free” the memory
   use, intrinsic :: iso_c_binding, only : c_char, c_size_t, c_ptr, c_f_pointer
   interface
      function foo( len ) result(r) bind(C, name="foo")
         import :: c_char, c_size_t, c_ptr
         ! Argument list
         integer(c_size_t), value, intent(in) :: len
         ! Function result
         type(c_ptr) :: r
      end function 
      subroutine free( ptr ) bind(C, name="free")
      ! C stdlib function
      ! void free(void* ptr);
         import :: c_ptr
         ! Argument list
         type(c_ptr), intent(in), value :: ptr
      end subroutine
   end interface
   type(c_ptr) :: cs
   integer(c_size_t) :: len
   len = 13
   cs = foo( len )
   block
      character(kind=c_char, len=len), pointer :: s
      call c_f_pointer( cptr=cs, fptr=s )
      print *, s
      s => null()
   end block
   call free( cs )
end
  • Program response using a processor
C:\temp>gfortran -c c.c

C:\temp>gfortran -c -ffree-form p.f

C:\temp>gfortran p.o c.o -o p.exe

C:\temp>p.exe
 Hello World!

C:\temp>
1 Like
12

Re Q1, isn’t it illegal to return an array (not a pointer) with bind-C function…? (according to the result of compiler explorer and this stackoverflow page) Maybe possible to return an array if it is included into a derived type (“Option 4” in a comment of the SO page)?

module test_m
use iso_c_binding
contains
function foo() result(res) bind(C)
integer(c_int) :: res(10)
end function
end module

> ifx, ifort
error #8523: A result variable of a BIND(C) procedure cannot be an array.
> gfortran
Error: Return type of BIND(C) function 'res' at (1) cannot be an array
13

W.r.t Q1, you are correct. It is not possible to return an array. The reason is it would violate item (2a) on what procedures are interpoerable,

A Fortran procedure interface is interoperable with a C function prototype if
(1) the interface has the BIND attribute,
(2) either
(a) the interface describes a function whose result is a scalar variable that is interoperable with the result of the prototype or
(b) the interface describes a subroutine and the prototype has a result type of void

The complete rules (for F2018) are found under section 18.3.6 of the following document: https://j3-fortran.org/doc/year/18/18-007r1.pdf (PDF, 8.86 MB)

Including the array in a scalar derived type is possible, if the derived type is also interoperable (say if the length of the array is fixed). If the members of the derived type aren’t interoperable, then an “opaque” handle-based solution must be used.

W.r.t Q2, that is incorrect. A char * result is a scalar variable interoperable with type(c_ptr). As @FortranFan has noted, there are some unspecified assumptions here about who is responsible for the memory of the char * result.

1 Like
14

When it comes to interoperability with a C function, a point to keep in mind always is the function prototype in C itself!

In the context of arrays in C, considering my points re: pointers and arrays in C upthread, ask yourself when has anyone been able to return anything other than a pointer itself as a function result when the intent might have been an array?!!!

And if you can’t, why expect the same when it comes to an INTERFACE for an interoperable function in Fortran?!

15

Hi FortranFan, though I can’t understand your message very well (partly because I can’t understand subtle nuances(?) because I am not a native speaker), I guess your meaning of “array” here is probably different from my meaning of “array”. In your case, I guess yours is referring to the underlying array data, while mine is referring to the “array” as defined in the C language. (For example, even C functions cannot return an “array”, it can effectively return the underlying data in various way, as in “Options” mentioned in the following SO page. I guess your comments are probably related to how such Options work (from the Fortran side)).

By the way, I am interested in why such a “restriction” existed (exists) historically in the C language… Even if the return variable is treated with stack, if a type containing a fixed-size array is OK to return, it seems that the latter can cause a similar problem for small stack size (in old machines). So the “restriction” comes more about the type system of C…? (Indeed, I don’t know how to write a C function that returns an array of two elements of int ). Though I’ve never tried, F77 also had a similar limitation before…?

It is also interesting that newer languages (than C) can return fixed-size arrays directly, like…

– modern Fortran [compiler explorer]
– D [compiler explorer]

16

Yes, an f77 function could not return an array. That only became possible with f90 and later.

I’m not fluent with C, but my understanding has always been that “arrays” in C are just an alternative notation for pointer arithmetic. For example, if you have an int array a and an int index j, then the C programmer can write a[j] or j[a] and they both evaluate to the same expression *(a+j). Right?

1 Like
17

Thank you, @FortranFan , @ivanpribec.

Regarding que context of Q1: I had skimmed over chapter 19 of Metcalf et al. and had (incorrectly) not picked up any obvious restrictions against returning arrays via functions. For sure, all the examples in the book were with subroutines, but I just naively assumed that if this is legal:

subroutine foo(result) bind(C)
  import :: c_int
  integer(c_int), dimension(10), intent(out) :: result
end subroutine

so must this:

function foo() bind(C) ! ILLEGAL
 import :: c_int
 integer(c_int), dimension(10) :: foo
end function

So, when the linter yelled at me, I even thought the linter had a bug! 🫨

The “trick” explained by @FortranFan works well, and I’ve implemented it here, Case 2.

I have two additional questions. In a interface such as:

subroutine add_vectors(length, vec1, vec2, result) bind(C)
    import :: c_int
    implicit none
    integer(c_int), value :: length
    integer(c_int), dimension(*), intent(in) :: vec1, vec2
    integer(c_int), intent(out) :: result(*)
end subroutine
  1. Is there any difference between these declarations?
...
    integer(c_int), dimension(*), intent(in) :: vec1, vec2
    integer(c_int), dimension(:), intent(in) :: vec1, vec2
    integer(c_int), dimension(length), intent(in) :: vec1, vec2
...
  1. Do the intent() attributes really do something? I feel good about adding them, but I wonder if they are actually doing something. Metcalf et al. only uses intent() on a few particular cases, e.g. when passing type(c_ptr) as input…
18

Note the second listing with an assumed-shape array (.. dimension(:), .. :: vec1 ..) is not directly interoperable with a C function and data types inherent to C. You can look into extended interoperability features in Fortran starting with Fortran 2018 for actions on the C side to achieve interoperability with assumed-shape arrays.

The third option with explicit-shape array (dimension(length)..), though interoperable, can give a false sense of security regarding the dimensions of the array and the memory. In actuality, it works pretty much the same as your first listing.

On the Fortran side, they serve as programmer aids with the all too important compile-time checks if the semantic restrictions are violated e.g., when an argument intended as output must be defined in an interoperable procedure but the input doesn’t matter, marking it as INTENT(OUT) can help the programmer.even though such a thing is irrelevant to C.

On the C side, look into const function parameters, you can find some parallels and do some pairing with INTENT(IN) on the Fortran side with interoperable procedures.

Note there is no direct and common enforcement, you have two different systems doing their own things independently, the so-called Fortran processor and the companion C processor.

19

On the fortran side, this assumed shape declaration requires an explicit interface. That is because a different argument association is required, one that allows the type, kind, and ranks of the actual and dummy arguments to match and for array extents to be passed to the dummy arguments. None of that can occur with the usual conventions with a C calling program.

20

On the other hand, assuming the procedure is implemented in Fortran, debug builds instrumented with runtime checks such as -fcheck=bounds (gfortran), -check bounds (ifort, ifx), will only work with the explicit-shape array. This doesn’t prevent misuse at the call site, but is still better than nothing IMO.

The second aspect I find nice is the relation between dummy arguments is captured by the declarations, and not just a comment (or “docstring” as Python calls them). An IDE can perhaps do something with this information.