Can a scalar C function be called from Fortran with array arguments?

c, fortran, fortran-iso-c-binding, fortran2003, interop

Solution

After some searching, I found that this is both possible and quite straightforward, using Fortran 2003 and the `iso_c_binding` intrinsic module; I thought I should document this here. An interface with the `bind` attribute can still be `pure`, so it can be referenced within a Fortran procedure which is itself `elemental`.

The following is a short example with the C99 macro `isnan`, which returns a non-zero integer value if its argument is NaN. Of course, this can also be done for any user-defined C/C++ function that is free of side effects.

elemental function isnan(val)
  use, intrinsic :: iso_c_binding
  implicit none

  real(c_double), intent(in) :: val
  logical(c_bool) :: isnan

  interface
    pure function isnan_C(val) bind(C, name = 'isnan')
      import
      ! Pass the parameter by value:
      real(c_double), value, intent(in) :: val
      integer(c_int) :: isnan_C
    end function
  end interface

  isnan = isnan_C(val) /= 0_c_int
end function

Example output using an array with some NaN values:

program main
  use, intrinsic :: iso_c_binding
  implicit none

  real(c_double) :: nan(2,2) = 0.
  nan(:,2) = nan(1,1)/0.

  write(*,'(4F6.2 / 4L6)') nan, isnan(nan) ! Output: 0.00 0.00  NaN  NaN
                                           !            F    F    T    T
end program

Problem

I recently came across a situation where I wanted to call a C function from Fortran, because of a useful snippet of C code. For convenience in array operations, I wanted to be able to call this function with array arguments as well, but it only accepted scalar arguments. In Fortran, one can of course simply declare a procedure to be `elemental` to achieve this, and one can declare an interface to a C procedure with `bind(C)`. However, since C does not have the concept of elemental procedures, the Fortran (2008) standard rules out this combination: C1246 An elemental procedure shall not have the BIND attribute. So, can this functionality be achieved in Fortran?

Original source