Fortran and C Plus Plus Interoperability with YAKL - mrnorman/YAKL GitHub Wiki

We provide the gator_mod Fortran module to interface with YAKL's internal device allocator. To use it, you'll first need to make all "automatic" fortran arrays into allocatable arrays and explicitly allocate them:

real var1( lbnd_x:ubnd_x , lbnd_y:ubnd_y , lbnd_z:ubnd_z )
real var2(nx,ny,nz)

will become:

real, allocatable :: var1(:,:,:)
real, allocatable :: var2(:,:,:)
...
allocate(var1( lbnd_x:ubnd_x , lbnd_y:ubnd_y , lbnd_z:ubnd_z ))
allocate(var2(nx,ny,nz))
...
deallocate(var1)
deallocate(var2)

Next, to interface with gator_mod, you'll need to transform all allocatables into pointer, contiguous. The contiguous attribute is recommended because some compilers will perform more optimizations on Fortran pointers when it is present. The resulting code is:

use gator_mod, only: gator_allocate, gator_deallocate

real, pointer, contiguous :: var1(:,:,:)
real, pointer, contiguous :: var2(:,:,:)
...
call gator_allocate( var1 , (/ ubnd_x-lbnd_x+1 , ubnd_y-lbnd_y+1 , ubnd_z-lbnd_z+1 /) , &
                            (/ lbnd_x , lbnd_y , lbnd_z /) )
call gator_allocate( var2 , (/nx,ny,nz/) )
...
call gator_deallocate(var1)
call gator_deallocate(var2)

Note that YAKL does support non-1 lower bounds in Fortran in gator_allocate(). The first set of indices specify the total extents of each dimension. The second set of indices specify the lower bounds of each dimension in case they aren't the default 1 lower bounds of Fortran. No, it's not the most convenient looking syntax for the user, but it is easier to implement this way :).

If a Fortran routine uses module data, when porting to C++, it is often easiest to mirror the Fortran practice. In this case, it is best to pass to allocate the Fortran data with gator_allocate and then immediately pass it to a C++ wrapping function, which will wrap it in unmanaged (non-owned) YAKL Array classes. With the above data, you will have the following code to do that:

module cpp_interface_mod
  use iso_c_binding
  interface
    subroutine wrap_arrays(var1,var2) bind(C,name="wrap_arrays")
      implicit none
      real(c_float), dimension(*) :: var1, var2
    end subroutine wrap_arrays
  end interface
contains
  ! All scalars can be directly bound to C
  real(c_float) , bind(C) :: scalar1, scalar2, scalar3
  
  ! Parameters cannot use bind(C), but rather must be replicated in
  ! a C++ header file with the constexpr keyword
  integer(c_int), parameter :: nx=DOM_NX, ny=DOM_NY, nz=DOM_NZ
  integer(c_int), parameter :: lbnd_x=-1, ubnd_x=nx+2
  integer(c_int), parameter :: lbnd_y=-1, ubnd_y=ny+2
  integer(c_int), parameter :: lbnd_z=-1, ubnd_z=nz+2
  
  ! Pass all arrays to a wrapper function to wrap them in Fortran-style Array objects in C++
  ! It's best to do this even for automatic arrays as well
  real(c_float) :: var1( lbnd_x:ubnd_x , lbnd_y:ubnd_y , lbnd_z:ubnd_z )
  real(c_float) :: var2( nx , ny , nz )
end module cpp_interface_mod

Fortran assumed size arrays (i.e., dimension(*)) are very convenient because an array of any dimension can be legally passed to the subroutine. Since Fortran passes by reference by default, the location of the data is passed to this array.

All variables in Fortran modules need to be moved to iso_c_binding types such as: real(c_float), real(c_double), integer(c_int), and logical(c_bool). All module-level scalars can directly be bound to C with bind(C). Arrays, however, need to be passed to a wrapping routine.

In C++, you would have a header file, fortran_data.h with the following:

#pragma once

int constexpr nx=DOM_NZ;
int constexpr ny=DOM_NZ;
int constexpr nz=DOM_NZ;
int constexpr lbnd_x=-1; 
int constexpr lbnd_y=-1; 
int constexpr lbnd_z=-1; 
int constexpr ubnd_x=nx+2; 
int constexpr ubnd_y=ny+2; 
int constexpr ubnd_z=nz+2; 

typedef yakl::Array<float,3,yakl::memDevice,yakl::styleFortran> real3d;

extern "C" void wrap_arrays(float *var1, float *var2);

// Declare Array wrappers defined in fortran_data.cpp
extern real3d var1, var2;

// Declare external scalars defined in Fortran code
extern int scalar1, scalar2, scalar3;  

#endif

And you would have a source file, fortran_data.cpp, with:

#include "fortran_data.h"

real3d var1, var2;

extern "C" void wrap_arrays(float *var1_p, float *var2_p) {
  // These create un-owned YAKL Arrays using existing allocations from Fortran
  var1 = real3d( "var1" , var1_p , {lbnd_x:ubnd_x} , {lbnd_y:ubnd_y} , {lbnd_x:ubnd_x} );
  var2 = real3d( "var2" , var2_p , nx, ny, nz   );
}

Notice that because of the YAKL Fortran-style Array class, you do not need to change the way you index these arrays in the C++ code at all. It's column-major ordering like Fortran, it defaults to lower bounds of 1, and it supports non-1 lower bounds as well. In the end, you often have C++ code that looks nearly identical to your previous Fortran code, with the exception of any Fortran intrinsics not supported by YAKL.

Any time Fortran data is passed by parameter, you can use the un-owned Array constructors to wrap them just as seen above in the wrap_arrays functions.

We've already seen how to pass Fortran arrays to C++. For scalars, though, consider the following Fortran function:

module mymod
contains
  function blah(n,x,which,y) result(z)
    integer(c_int) , intent(in   ) :: n
    real(c_float)  , intent(in   ) :: x
    logical(c_bool), intent(in   ) :: which
    real(c_double) , intent(  out) :: y
    real(c_float)                  :: z
    ! code
    ! y = ...
  end function
end module

When you port this to C++, you'll change the Fortran code to the following:

module mymod
  interface 
    function blah(n,x,which,y) result(z)  bind(C, name="blah")
      integer(c_int) , intent(in   ), value :: n
      real(c_float)  , intent(in   ), value :: x
      logical(c_bool), intent(in   ), value :: which
      real(c_double) , intent(  out)        :: y
      real(c_float)                         :: z
    end function
  end interface
contains
  ! You remove the code from here
  ! And you only put the function *header* above (no code)
end module

And you'll have a C++ function:

extern "C" float blah( int n , real x , bool which , double &y ) {
  // code
  // real y = ...
  return y;
}

Notice that in the Fortran interface, any scalar with intent(in) must be passed by value and not by reference (i.e., the value fortran keyword must be added). However, since the scalar y has intent(out), we must pass it by reference, meaning we do not use the value keyword. On the C++ side of things, we accept n, x, and which as a simple int, real, and bool. But we accept y as &y, meaning when we change its value, the change lives outside the function like we want it to.

Regarding return values for functions, do not return a reference like you do with dummy arguments. If you added the & to the function return, you would essentially be returning a reference to a locally scoped variable, which won't have any meaning outside the function. The iso_c_binding handles things appropriately for you, so just return the appropriate type by value in C++, and map it to a simple return of the same type in the Fortran interface block.

The following table can help you convert Fortran parameter types to Fortran interfaces, and C++ dummy arguments:

Fortran datatype	Fortran interface	C++ dummy argument
integer(c_int), intent(in)	integer(c_int), value	int
real(c_double), intent(out)	real(c_double)	double &
real(c_float), dimension(...)	real(c_float), dimension(*)	float *
logical(c_bool), dimension(...)	logical(c_bool), dimension(*)	bool *