When FORTRAN met OpenCL – A Love Story? (of Integration)

Both of the projects Lukas and I are working on, here at NBI Copenhagen, are based on FORTRAN. Either Photon-Plasma simulations or CESM for climate modeling (as most climate modeling software around). Certainly, a vast majority of scientific code is also written in FORTRAN.

One of the roles in our projects is to add GPGPU capability through OpenCL to accelerate various routines. For those who are familiar, OpenCL is a C based API. Therefore, to integrate OpenCL into FORTRAN we have to bridge with C somehow. Anyway, the two of us are not fluent native FORTRAN speakers from elementary school and feel more comfortable in C for advanced programming and due to other technical issues with libraries being used we have to stick with C/C++.

The rest of the post describes the process of our development and early integration. Those who will be patient reading through without skipping, will enjoy and hopefully benefit, to the end.

As an old Hebrew saying: “If you worked hard and achieved, believe it” (is possible).
Sometimes good things come by chance.

FORTRAN and C

Are pretty much different languages. Although similar in mind and programming constructs, sometimes the difference can make them very distinct from each other. FORTRAN suits scientific computing very well, but when it comes to working with hardware issues, C is better and closer to the “metal”. One good example is C support for pointers, while in FORTRAN they have different meaning and implementation.

Using OpenCL from FORTRAN

As mentioned earlier, OpenCL is C API, and can be accessible directly from C, but not from FORTRAN. OpenCL in turn, is usually used in a kernel execution model, where some specific accelerated routine is implemented on the GPU and then interfaced from host code (C/C++ etc.).

Therefore, to use OpenCL from FORTRAN, we have several approaches:

Option 1
Write algorithms in C, and supply FORTRAN with a simple function to call and perform the computation, something like calling: “multiply_matrices“.

The advantage is by hiding all the burden from FORTRAN required to initialize the device, pointers, GPU buffers and other management. The code is simpler (it’s C and native to OpenCL) and can be easily maintained. However, doing so incorporates C source into our codebase, so we have to manage two different languages.

How is it possible to call a C function from FORTRAN?
Assume the following FORTRAN function call:

! Define variables.
real :: a(12, 12)
real :: b(12, 10)
real :: c(12, 10)

! Call compute function.
call multiply_matrix(a, b, c)

Where a, b and c denote a FORTRAN matrix (2D array) and the operation is simply: c = a * b.

Now, a matching C implementation of the function above would be:

void multiply_matrix_(float *a, float *b, float *c)
{ // Compute part.
}

Those with a sharp eye, could notice the additional underscore (_) in the C function name. It’s not a mistake, but a supported way to write a function in C that can be recognized later by the FORTRAN linker. It works with the compilers we tested so far (GNU, Intel and others).

Users still have to be warned, as FORTRAN stores arrays in column based ordering, while C/C++ in row based. This can lead to some serious problems if not handled correctly – while matrices occupy the same memory size, elements are not indexed in the same order!

Name Mangling is the process of converting a function name from our source code, into native representation. We will get later to the mangling point – it will play an important role soon.

Short Note: the above is impossible in C++ (compared to C). C function names are saved the way they are written, but C++ does more sophisticated things to allow function overriding (the function above could be saved as: _Z__multiply_matrix_F_F_F, not that human readable and also compiler dependent).

Such an approach was already taken 3-4 years ago, when porting cloud micro-physics code in WRF climate model from FORTRAN to CUDA (by NVIDIA and NCAR). Most of the code was written in C and then interfaced from FORTRAN.

Option 2
We can create a wrapper for frequently used OpenCL API functions and expose them to FORTRAN.

For example, a function that queries which devices are available in the system: clGetDeviceIDs(…), can be wrapped in a similar fashion:

cl_int clGetDeviceIDs_(cl_platform_id *platform, 
    cl_device_type *dev_type, 
    cl_uint *num_entries, 
    cl_device_id *devices, 
    cl_uint *num_devices)
{ // Call the real clGetDeviceIDs function and pass parameters.
}

Again, for those with a sharp eye and familiar with the OpenCL API, might wonder why the function accepts all its parameters as pointers, even for scalar values? That’s another difference between FORTRAN and C – FORTRAN passes all parameter as pointers (by reference).

Doing so for all the OpenCL functions can be a tedious task, and the same pitfalls from Option 1 still remain here, having to maintain C code together with FORTRAN. However, we have greater flexibility in our FORTRAN code.

Pointer Handling
Something that was left behind, intentionally.

When wrapping OpenCL API to FORTRAN we will have to have them, where option 1 didn’t require explicit use of pointers at all. So, how to have C pointers in FORTRAN? FORTRAN has no support for a type that changes its size according to 32/64 bit processor architecture (like size_t or void*). To do that, we have to introduce a new type to represent a pointer:

!
! Define a pointer type based on compiler target architecture.
! The following definitions are based on the GNU C Preprocessor and can be
! seen on: http://gcc.gnu.org/onlinedocs/cpp/Common-Predefined-Macros.html
!
! Any definition of more complex types can be based on that (for example
! representing handles to a platform, device, context, etc.
type NATIVE_PTR
    !sequence
#if __SIZEOF_POINTER__ == 4
    integer*4 :: ptr
#elif  __SIZEOF_POINTER__ == 8
    integer*8 :: ptr
#endif
end type NATIVE_PTR

Another option exists for using something called cray-pointers, developed in LLNL, but less elegant for our purposes. This feature is also supported by all major compilers.

I’ll save some surprise for later.

What to do?
Well, the approach we chose to take is a combination of the above, clearly we need FORTRAN management capabilities for GPU resources, but the algorithmic flow can be written in C for now.

But hold on! The story doesn’t end here.

A Good Software Engineer

Well, up to this point we have a set of OpenCL interface functions wrapped for FORTRAN in C, and some algorithmic code.

To make things more organized and robust, the next step was to create a FORTRAN module for OpenCL interface definitions. Similar to a C header file or a library, it allows code reuse to consume these definitions later:

subroutine my_function
use clfortran
! Now we can call all the OpenCL functions :)
end subroutine my_function

A Disaster
Was waiting around the corner. Using modules, function names are mangled differently than we prepared for, so two functions with the same name can coexist in different modules. Good for them, for us it was a headache.

Now, our C interface should have looked like this:

cl_int clfortran_MOD__clGetDeviceIDs(cl_platform_id *platform, 
    cl_device_type *dev_type, 
    cl_uint *num_entries, 
    cl_device_id *devices, 
    cl_uint *num_devices)
{ // Call the real clGetDeviceIDs function and pass parameters.
}

You would say, OK, only a prefix and can be handled. But that’s not the only problem, as the prefix given, differs from compiler to compiler! (Intel ifort gave something else).

Full of good spirit we went on to study how to figure the prefix the compiler would give in compile time, so we can handle it dynamically for all compilers. It’s something that can be done, although in perspective, might be better to avoid.

Last Notes (Happy End)

Upon searching for a solution, we came through two results:

1. A project called FORTRANCL, which provides a C based wrapper to most OpenCL functions and for FORTRAN. Similar to the approach we took in the first place.
2. The ISO_C_BINDING extension to FORTRAN, which aims to provide better FORTRAN <-> C interface support. It was introduced in FORTRAN 2000/2003.

The FORTRAN C extension, was in the back of our minds, but really shadowed, as FORTRANER’s are always afraid not to break any FORTRAN 90/95 conformance. This extension allows using the same types in FORTRAN as we define them in C, and the more important, it provides control over the name mangling of function names!!!

Our interface definition code to call OpenCL API looks something like this:

integer(c_int) function clGetDeviceIDs(platform, device_type, num_entries, devices, num_devices) BIND(C, NAME='clGetDeviceIDs')

    USE ISO_C_BINDING

    ! Define parameters.
    type(c_ptr), VALUE :: platform
    integer(c_int64_t), VALUE :: device_type
    integer(c_int), VALUE :: num_entries
    type(c_ptr), VALUE :: devices
    integer(c_int), intent(out) :: num_devices
end function

Very much like a DllImport in C#.

Advantages:

1. OpenCL Standard Conformance – All parameters are defined according to the original standard and can be used from FORTRAN.
2. No C code is involved – it’s all pure FORTRAN definitions (the compiler takes care of handling any issues behind the scenes, passing the correct parameters etc).

The decision to conform to OpenCL C API was to ease the life of programmers, so they can refer to the same (and only) manual and figure parameter meanings easily.

So, we are now in an ongoing process to extend the CLFORTRAN API and expose more OpenCL functions. Hopefully it would allow other scientists take advantage of the GPGPU power from their existing codes. Being a standardized and supported way, it works on GNU, Intel and IBM XL compilers (Bluegene, AIX and Linux) to name a few, and also “one step for mankind into the year 2000“.

If you’d notice our module naming, luckily since the beginning there was no conflict with the FORTRAN CL project :).

A short “Hello, World!” example of how to use OpenCL from FORTRAN (queries devices and prints their names):

program test_cl
    use clfortran
    use ISO_C_BINDING
    implicit none

    integer(c_int) :: err
    integer(c_int) :: num_platforms
    integer(c_int) :: i
    integer(c_int64_t) :: device_type
    integer(c_int) :: num_devices
    type(c_ptr), allocatable, target :: platform_ids(:)

    ! Get platform count to allocate an array.
    print *, clGetPlatformIDs(0, C_NULL_PTR, num_platforms)
    print *, 'Num Platforms: ', num_platforms

    ! Allocate an array to hold platforms.
    allocate(platform_ids(num_platforms))

    ! Get platforms IDs.
    err = clGetPlatformIDs(num_platforms, C_LOC(platform_ids), num_platforms)

    ! Loop over all platforms and query devices.
    do i = 1, num_platforms
        ! Iterate over platforms and get number of devices.
        print *, 'Platform: ', i

        ! Get device count.
        device_type = -1        ! ALL, a constant will be added later
        print *, clGetDeviceIDs(platform_ids(i), device_type, 0, C_NULL_PTR, num_devices)
        print *, 'Num Devices: ', num_devices
    end do
end program test_cl

In future posts we will discuss more about usage of OpenCL through our projects and the way it was integrated into the original source code.