I would rather say that the explicit-shape convention (n) is the programmer telling the compiler that the length is compatible with whatever is the actual argument. Once the programmer does that, the compiler can then enforce bounds checking with that array shape and whole array operations are allowed, including i/o with the whole array. Those things are not allowed with dummy arguments with assumed size declarations (*).
Again, thank you @FortranFan, @RonShepard, @ivanpribec, for the helpful explanations.
The context of my question – call C from Fortran – probably got lost along the posts. The interface I mentioned is an “import” interface, i.e. one that is meant to bind to an existing C function. Therefore, IMO dimension(:) and dimension(length) are meaningless. What surprises me is that the compiler does not complain, giving one a sense of looseness. Maybe this is so, because a interface itself does not know what it is meant for (import or export), but I would still argue the compiler could figure that out… Anyway, I believe I now know the proper way to do it, even if the compiler does not enforce it.
I am trying to deal with another case, where I call a C function apply_function_to_vector() that takes another function square as argument. This is what my code looks like - I suppose I am missing a step to properly create a pointer to square. Can you please help?
program fprogram
use, intrinsic :: iso_fortran_env, only : real64
use, intrinsic :: iso_c_binding, only: c_double, c_int, c_ptr, c_funptr, c_funloc
implicit none
interface
! this is the C function to be called from Fortran
subroutine apply_function_to_vector(f, input_vector, output_vector, length) bind(C)
import :: c_double, c_int, c_funptr
implicit none
type(c_funptr), value :: f
real(c_double), dimension(*), intent(in) :: input_vector
real(c_double), dimension(*), intent(out) :: output_vector
integer(c_int), value :: length
end subroutine
end interface
integer, parameter :: length = 3
real(real64), dimension(3) :: input_vector, output_vector
input_vector = [1d0, 2d0, 3d0]
call apply_function_to_vector(c_funloc(square), input_vector, output_vector, length)
print *, "Output Vector:"
print *, output_vector
contains
! this is the Fortran function to be passed as argument to `apply_function_to_vector`
real(real64) function square(x)
real(real64), intent(in) :: x
square = x**2
end function
end program fprogram
Calling convention mismatch between C and Fortran also applies to function pointers passed around the different processors.
In your cmodule.h the function pointer has interface:
typedef double (*func_ptr)(double);
while your Fortran function square is the equivalent to
double (*func_ptr)(double *);
Try changing the Fortran function to accept real(real64), value :: x. This should fix the issue.
1 Like
In addition, the Fortran callback should be marked as bind(c)
1 Like
Good point. I had added it myself, but forgot to mention it. Thanks for pointing it out.
1 Like
One downside of the interface that accepts a scalar function is that due to the cross-language calls, you are unlikely to benefit from inlining and/or vectorization. Maybe inter-procedural optimization can help, but we’d need to verify this to to be sure.
The other way is to modify the callback, and make it the duty of the callback to evaluate all elements:
// apply_vfun.c
#include <stdio.h>
// evaluate y_i = f(x_i), for i = 0, n-1
typedef void (*vfun)(int,
const double *restrict,
double *restrict);
void apply_vfun(vfun f, int n, const double *x, double *y)
{
printf("applying vector fun for n = %d\n", n);
f(n, x, y);
}
Here is IMO a good reason why the dimension(n) syntax is not meaningless – it allows the use of array expressions in the Fortran callback, which are amenable to auto-vectorization. If array syntax is not enough, then either an elemental function or an explicit loop with !$omp simd will do the trick.
! demo.f90
program demo
use, intrinsic :: iso_c_binding
implicit none
abstract interface
subroutine vfun(n,x,y) bind(c)
import c_int, c_double
implicit none
integer(c_int), value :: n
real(c_double), intent(in) :: x(n)
real(c_double), intent(out) :: y(n)
end subroutine
end interface
interface
subroutine apply_vfun(f,n,x,y) bind(c)
import c_int, c_double, vfun
implicit none
procedure(vfun) :: f
integer(c_int), value :: n
real(c_double), intent(in) :: x(n)
real(c_double), intent(out) :: y(n)
end subroutine
end interface
integer, parameter :: n = 32
real(c_double) :: x(n), y(n)
integer :: i
call random_number(x)
! y := f(x)
call apply_vfun(my_vfun,n,x,y)
write(*,'(A)') "Results:"
write(*,'(*(G0,2X,G0,/))') [(x(i),y(i),i=1,n)]
contains
! Callback passed to C
subroutine my_vfun(n,x,y) bind(c)
integer(c_int), value :: n
real(c_double), intent(in) :: x(n) ! Note the x(n) syntax, ...
real(c_double), intent(out) :: y(n)
! ...that allows us to use array expression.
y = sin(x) * x**2
! if x(*) were used, then we'd need to use the
! explicit x(1:n) syntax instead
end subroutine
end program
If you were to look at this in Compiler Explorer, and compiler options like -O3 -march=skylake, my_vfun contains the hot loop,
.L9:
vmovupd ymm0, YMMWORD PTR [r14]
add r14, 32
add r15, 32
vmulpd ymm2, ymm0, ymm0
vmovapd YMMWORD PTR [rbp-80], ymm2
call _ZGVdN4v_sin
vmulpd ymm0, ymm0, YMMWORD PTR [rbp-80]
vmovupd YMMWORD PTR [r15-32], ymm0
cmp r14, rbx
jne .L9
that uses the 256-bit AVX registers (we evaluate 4 elements at the same time), including a vector implementation of the sine function.
This is one the reasons I like Fortran.
2 Likes