UP | HOME

Calling Fortran from C

[2013-04-17 Wed]

Eventually, I want to be able to call Fortran from Python. However, as a stepping stone, I first call Fortran from C and then use that for the next step. As an example I'll be using a fractal calculating program. The code can be accessed with 

git clone https://maurow@bitbucket.org/maurow/mauro_learning_fortran.git
cd mauro_learning_fortran
git checkout v0.1

and is in the folder mandel/.

The Fortran side

The function calc_num_iter in Fortran module mandel computes whether the points in the grid (re,im) are part of the Mandelbrot set.

mandel.f90: 

module mandel
  implicit none

  integer, parameter :: dp=kind(0.d0) ! double precision

contains

  pure function mandel_frac(z, c) result(out)
    ! The Mandelbrot function z -> z^2 + c
    complex(dp), intent(in):: z, c
    complex(dp):: out
    out = z**2 + c
  end function mandel_frac

  pure function calc_num_iter(re, im, itermax, escape) result(out)
    ! Iterates on mandel_frac
    real(dp), intent(in):: re(:), im(:), escape
    integer, intent(in):: itermax
    integer:: out(size(re), size(im))
    integer:: ii, jj, kk, itt
    complex(dp):: zz, cc

    do ii=1,size(re)
       do jj=1,size(im)
          zz = 0
          cc = cmplx(re(ii), im(jj), dp)
          itt = 0
          do kk=1,itermax
             zz = mandel_frac(zz, cc)
             itt = kk
             if (abs(zz)>escape) then
                exit
             end if
          end do
          out(ii,jj) = itt
       end do
    end do

  end function calc_num_iter

end module mandel

This can, of course, be used in Fortran like so:

run_mandel_from_fortran.f90: 

program mandelbrot
use mandel
implicit none

integer, parameter :: nre=31, nim=21
real(dp), parameter :: rer(2)=[-2._dp, 1._dp], imr(2)=[-1._dp,1._dp]
integer :: ii, jj, kk, out(nre,nim), itermax=99
real(dp) :: re(nre), im(nim), escape=2._dp
complex(dp) :: zz

re = [ (dble(ii)/(nre-1)*(rer(2)-rer(1)) + rer(1) , ii=0, nre-1, 1) ]
im = [ (dble(ii)/(nim-1)*(imr(2)-imr(1)) + imr(1) , ii=0, nim-1, 1) ]

out = calc_num_iter(re, im, itermax, escape)

! plot array as x/y coordinates
do jj=nim,1,-1
   print '(30000I3.2)', out(:,jj)
end do

end program mandelbrot

Compile it with 

gfortran mandel.f90 run_mandel_from_fortran.f90 -o fout

The wrapper

To make the function calc_num_iter available in C it needs to be wrapped using the iso_c_binding module (in the Fortran 2003 standard). Resources for the dumb Fortran programmer, like myself, are a bit sparse (e.g. gfortran docs). The fortran wrapper is

mandel_wrap.f90: 

module mandel_wrap
  ! to wrap calc_num_iter for use in C

  use iso_c_binding, only: c_double, c_int
  use mandel, only:  calc_num_iter

  implicit none

contains

  ! need to make a subroutine as only scalars can be returned
  subroutine c_calc_num_iter(nre, re, nim, im, itermax, escape, out) bind(c)
    real(c_double), intent(in):: re(nre), im(nim)
    real(c_double), intent(in), value:: escape
    integer(c_int), intent(in), value:: itermax, nre, nim
    ! note that in C the indices will be reversed!:
    integer(c_int), intent(out):: out(nim, nre)
    ! thus the transpose here:
    out = transpose(calc_num_iter(re, im, itermax, escape))
  end subroutine c_calc_num_iter

end module mandel_wrap

There are a few things to note:

  • bind(c) needs to be appended to the wrapper function
  • the c_ variable types are needed to make them C interoperable
  • The attribute value specifies that this variable uses call by value, which is the standard for C (but not Fortran).
  • I define out with the indices reversed compared to calc_num_iter as Frotran uses column-major and C row-major array storage. I then transpose the result of calc_num_iter to make it fit into out.

The C side

Now it is possible to call c_calc_num_iter from C.

run_mandel_from_c.c: 

#include <stdio.h>
#include <math.h>
#define NRE 31
#define NIM 21

int main ( void ) {
  int i, j, itermax=99;
  // NOTE: index not reversed!
  int res[NRE][NIM];
  double re[NRE], im[NIM],  escape=2., dx, dy;

  // make real and imaginary vectors
  dx = 3.0/(NRE-1);
  re[0] = -2.;
  for (i=1; i<NRE; i++)
    re[i] = re[i-1] + dx;
  dy = 2.0/(NIM-1);
  im[0] = -1.;
  for (i=1; i<NIM; i++){
    im[i] = im[i-1] + dy;
    if (im[i]<1e-6 && im[i]>-1e-6)
      im[i] = 0.; // otherwise result is not equal to Fortran's
  }
  //printvec(NRE, re);
  //printvec(NIM, im);

  // call to fortran
  c_calc_num_iter(NRE, re, NIM, im, itermax, escape, res);

  // print result
  printarr(NRE, NIM, res, re, im);
  return 0;
}

This needs to be compiled: 

gcc -c run_mandel_from_c.c mandel.f90 mandel_wrap.f90

Produces the object files *.o, which need to be linked into an executable with 

gfortran run_mandel_from_c.o mandel.f90 mandel_wrap.f90 -o cout

Note the use of gfortran here. Running the same command using gcc does not work as the linker does not find the Fortran functions. Probably some flags for gcc could solve this?

Next I call c_calc_num_iter from python through a cython interface.

Links