Assume we need to write fortran and wish to call C code that has already been written and cannot be changed (e.g. a precompiled library). There are different levels of complexity possible depending on what is needed. Eg:
I have picked a function that gives us \(\pi\) on integration. Out of all such integrals, one with definite limits is required. Also, absence of singularities in the integrand is desirable. One such integral is: $$ \int_{0}^{1} \frac{4}{1+x^2} dx=\pi$$ The integrand is well behaved in the region of interest and will do the job for now.
There are, of course, more superior and efficient ways to calculate \(\pi\) numerically/computationally. For instance, the Ramanujan-type formulas or the "spigot" algorithms like the Bailey–Borwein–Plouffe formula or the Gauss-Legendre algorithm (different from the Gauss-Legendre quadrature), some of which have been used for getting a world record number of digits of \(\pi\).
Obviously, we need to know how the function is to be called. Thus, we require the name of the function, its return type, the list of arguments and the corresponding data types. The function prototype for the example here is given below. I have taken this C code from the wonderful www.mymathlib.com website.
/*a=lower limit
b=upper limit
f=pointer to function to be integrated
*/
long double
Gauss_Legendre_Integration_10pts( long double a, long double b, long double (*f)(long double) )
{...
}
The whole fortran program is given below with discussion following.
program Integrate
use, intrinsic :: iso_c_binding,cl=>c_long_double
implicit none
interface
function gl10(a, b, f) &
bind(c,name='Gauss_Legendre_Integration_10pts') !the function name case-sensitive because C is case-sensitive
import::cl,c_funptr
real(cl), intent(in), value :: a, b
type(c_funptr), intent(in), value :: f
real(cl) gl10
endfunction
function f1(x) bind(c)
import::cl
real(cl),intent(in), value ::x
real(cl) f1
endfunction
endinterface
type(c_funptr) :: f1_cptr
real(cl) :: a,b
a=0; b=1
f1_cptr = c_funloc (f1) !gets the C address of the fortran function f1
write(*,*)'error= ',gl10(a,b,f1_cptr)-4*atan(1.0_cl)
endprogram
function f1(x) bind(c)
use, intrinsic :: iso_c_binding,cl=>c_long_double
implicit none
real(cl),intent(in), value ::x
real(cl) f1
f1=4.0/(1.0+x**2)
endfunction
We have to use the right data types in fortran corresponding to those in C. These types are present in the module 'iso_c_binding'. We will be using two of these, i.e "c_long_double" and c_funptr. "c_long_double" is a kind of real, whereas 'c_funptr' is a derived datatype. I have just aliased the 'c_long_double' to 'cl' to reduce typing.
The interface block is used to make fortran understand how to call the function. The first function in the interface, 'gl10', is prototype for the C function that does the quadrature calculation. It is declared with the 'bind(C)' attribute along with the name of the function in C. By this interface, we give this C function, 'Gauss_Legendre_Integration_10pts', a different name in fortran, 'gl10'. Now, in fortran we can call this function as g10(...). The C function name that appears within 'bind' is case-sensitive.
There is some flexibility available in interfacing with the C functions. In general, in fortran, we can match the variables, pointers and function pointers in the argument list of a C function in two different ways. I do not know if one way is preferable over the other however.
The first two arguments in 'Gauss_Legendre_Integration_10pts' are scalar variables of type real(cl) and not pointers. As the convention in C is to pass scalars by value, the interface declares these arguments with the 'value' attribute. It is required because the C function expects a value, whereas a fortran-language call typically sends out the reference/address (usually). Now, if 'a' and 'b' in the C function definition were instead '*a' and '*b', the 'value' attribute should not be present in the interface. In our case we can edit the C code('a' to '*a' and 'b' to '*b' throughout), and remove the value attribute from 'a' and 'b' in the interface and verify that this approach also works.
Alternately, in the case with pointer arguments (like '*a','*b') for the C function, we can find the C address (using C_LOC(..)) of the corresponding variables in fortran and pass this address by value. The third argument, a C-style function pointer, is also declared with the 'value' attribute. This way the C function recieves a function address for the function pointer, f . For interoperability, for pointers that are dummy arguments, the 'value' attribute is usually required because void* matches TYPE(C_PTR), VALUE, while TYPE(C_PTR) by itself, i.e. without the 'value' attribute, matches void**. The same goes for function pointers. If the value attribute was not present in the interface for 'f', then '*f' in the C function would have to be replaced by '**f' throughout.
Here our choice of using the value attribute for all the dummy-arguments of the C function is dictated by our wish to leave the C code untouched
The function to be integrated, f1, is defined outside the scope of the main program. So, an interface has to be provided for this function as well. Since this function is being called from C code, it needs to be interoperable with C. So, this function is defined with the bind(C) attribute. The name of this function is only needed by the main fortran program. The C function access/calls this function using a pointer and so it doesnt need the name in this case.
But, in the general case this fortran function can be called from C code as well. Like before, The 'name' argument of bind can be provided if one wishes to call this function with a different name within the C code (see the end of the page).
C_FUNLOC(..) takes the name of a fortran function and returns a C function address. The type is a derived type called 'c_funptr' that is present in the iso_c_binding module. Here I have used the variable 'f1_cptr' to store the address of the fortran function 'f1'. Now that all the arguments have been defined, the C function is called in the program as 'gl10(a,b,f1_cptr)'.
To see how close the integral gets to \(\pi\) the error is calculated and printed. That wraps it up for now. The code can be downloaded from the links below.
Download Links:
To compile the code the following can be done for gcc/gfortran:
$>gcc -o gauss_legendre_10pts.o -c gauss_legendre_10pts.c
$>gfortran -o integrate.o -c integrate.f90
$>gfortran -o integrate.exe integrate.o gauss_legendre_10pts.o
Using built-in specs.
COLLECT_GCC=gfortran
COLLECT_LTO_WRAPPER=C:/TDM-GCC-64/bin/../libexec/gcc/x86_64-w64-mingw32/5.1.0/lto-wrapper.exe
Target: x86_64-w64-mingw32
Configured with: ../../../src/gcc-5.1.0/configure --build=x86_64-w64-mingw32 --enable-targets=all --enable-languages=ada,c,c++,fortran,lto,ob
jc,obj-c++ --enable-libgomp --enable-lto --enable-graphite --enable-cxx-flags=-DWINPTHREAD_STATIC --disable-build-with-cxx --disable-build-po
ststage1-with-cxx --enable-libstdcxx-debug --enable-threads=posix --enable-version-specific-runtime-libs --enable-fully-dynamic-string --enab
le-libstdcxx-threads --enable-libstdcxx-time --with-gnu-ld --disable-werror --disable-nls --disable-win32-registry --prefix=/mingw64tdm --wit
h-local-prefix=/mingw64tdm --with-pkgversion=tdm64-1 --with-bugurl=http://tdm-gcc.tdragon.net/bugs
Thread model: posix
gcc version 5.1.0 (tdm64-1)