SWIG and Ctypes¶
Interfacing Fortran code through C interface
I found this reference with good examples.
Problem: Collatz conjecture¶
Choose \(u_0\)
If \(u_k\) even, \(u_{k+1} \rightarrow \frac{u_k}{2}\) ;
If \(u_k\) odd, \(u_{k+1} \rightarrow 3 u_k+1\)
The Collatz conjecture is: For all \(u_0>0\) , the process will eventually reach \(u_k=1\).
The programs below compute number of steps (named flight) to reach \(f(u_0)\) for \(1\leq u_0 \leq N\), \(N\) given.
C program¶
%%file syracuse.c
#include <stdlib.h>
#include <stdio.h>
long syracuse(long n) {
long count = 0L ;
while (n > 1) {
if ((n&1)==0)
n /= 2;
else
n = 3*n+1;
count++;
}
return count ;
}
int main() {
const long N = 1000000;
double t1, t2;
long i , *flights ;
flights = (long*)malloc(N*sizeof(long));
for (i = 0; i <N; i++) flights[i] = syracuse(i+1);
return EXIT_SUCCESS;
}
Writing syracuse.c
%%bash
gcc -O3 syracuse.c
time ./a.out
real 0m0.102s
user 0m0.001s
sys 0m0.002s
Python program¶
%%time
from itertools import count
def syracuse(n):
x = n
for steps in count() :
if x & 1 :
x = 3*x+1
else:
x = x // 2
if x == 1:
return steps
N = 1000000
flights = [syracuse(i) for i in range(1,N+1)]
CPU times: user 21.1 s, sys: 262 ms, total: 21.4 s
Wall time: 21.5 s
Performances¶
The python syntax is simpler.
100 times slower
Solution : call the C function from python.
Ctypes¶
This is the C function we will call from python
%%file syrac.c
long syracuse(long n)
{
long count = 0L ;
while (n > 1)
{
if ((n&1)==0)
n /= 2;
else
n = 3*n+1;
count++;
}
return count ;
}
Writing syrac.c
Build the shared library
%%bash
gcc -fPIC -shared -O3 \
-o syrac.so syrac.c
%%time
import time
from ctypes import *
syracDLL = CDLL("./syrac.so")
syracuse = syracDLL.syracuse
flights = [syracuse(i) for i in range(1,N+1)]
CPU times: user 532 ms, sys: 8.51 ms, total: 540 ms
Wall time: 652 ms
Ctypes with Fortran module¶
If you change the fortran file you have to restart the kernel
%%file syrac.F90
module syrac_f90
use iso_c_binding
implicit none
contains
function f_syrac(n) bind(c, name='c_syrac') result(f)
integer(c_long) :: f
integer(c_long), intent(in), value :: n
integer(c_long) :: x
x = n
f = 0_8
do while(x>1)
if (iand(x,1_8) == 0) then
x = x / 2
else
x = 3*x+1
end if
f = f + 1_8
end do
end function f_syrac
end module syrac_f90
Writing syrac.F90
%%bash
rm -f *.o *.so *.dylib
gfortran -fPIC -shared -O3 -o syrac.dylib syrac.F90
bash: line 2: gfortran: command not found
---------------------------------------------------------------------------
CalledProcessError Traceback (most recent call last)
/var/folders/24/8k48jl6d249_n_qfxwsl6xvm0000gn/T/ipykernel_28719/4175552032.py in <module>
----> 1 get_ipython().run_cell_magic('bash', '', 'rm -f *.o *.so *.dylib\ngfortran -fPIC -shared -O3 -o syrac.dylib syrac.F90\n')
/usr/local/lib/python3.8/site-packages/IPython/core/interactiveshell.py in run_cell_magic(self, magic_name, line, cell)
2417 with self.builtin_trap:
2418 args = (magic_arg_s, cell)
-> 2419 result = fn(*args, **kwargs)
2420 return result
2421
/usr/local/lib/python3.8/site-packages/IPython/core/magics/script.py in named_script_magic(line, cell)
140 else:
141 line = script
--> 142 return self.shebang(line, cell)
143
144 # write a basic docstring:
/usr/local/lib/python3.8/site-packages/decorator.py in fun(*args, **kw)
230 if not kwsyntax:
231 args, kw = fix(args, kw, sig)
--> 232 return caller(func, *(extras + args), **kw)
233 fun.__name__ = func.__name__
234 fun.__doc__ = func.__doc__
/usr/local/lib/python3.8/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
185 # but it's overkill for just that one bit of state.
186 def magic_deco(arg):
--> 187 call = lambda f, *a, **k: f(*a, **k)
188
189 if callable(arg):
/usr/local/lib/python3.8/site-packages/IPython/core/magics/script.py in shebang(self, line, cell)
243 sys.stderr.flush()
244 if args.raise_error and p.returncode!=0:
--> 245 raise CalledProcessError(p.returncode, cell, output=out, stderr=err)
246
247 def _run_script(self, p, cell, to_close):
CalledProcessError: Command 'b'rm -f *.o *.so *.dylib\ngfortran -fPIC -shared -O3 -o syrac.dylib syrac.F90\n'' returned non-zero exit status 127.
from ctypes import *
syrac_f90 = CDLL('./syrac.dylib')
syrac_f90.c_syrac.restype = c_long
syrac_f90.c_syrac(1000)
%%time
N = 1000000
flights = [syrac_f90.c_syrac(i) for i in range(1,N+1)]
Faster than pure Python
We can call function from DLL windows libraries.
Unfortunately you need to adapt the syntax to the operating system.
SWIG¶
Interface file syrac.i for C function in syrac.c
%%file syrac.i
%module syracuseC
%{
extern long syracuse(long n);
%}
extern long syracuse(long n);
%%bash
swig -python syrac.i
Build the python module¶
Using command line
swig -python syrac.i
gcc `python3-config --cflags` -fPIC \
-shared -O3 -o _syracuseC.so syrac_wrap.c syrac.c `python3-config --ldflags`
With distutils
%%file setup.py
from numpy.distutils.core import Extension, setup
module_swig = Extension('_syracuseC', sources=['syrac_wrap.c', 'syrac.c'])
setup( name='Syracuse',
version = '0.1.0',
author = "Pierre Navaro",
description = """Simple C Fortran interface example """,
ext_modules = [module_swig],
)
import sys, os
if sys.platform == "darwin":
os.environ["CC"] = "gcc-10"
!{sys.executable} setup.py build_ext --inplace --quiet
import _syracuseC
syracuse = _syracuseC.syracuse
syracuse(1000)
%%time
N=1000000
flights = [syracuse(i) for i in range(1,N+1)]
References¶
Python Scripting for Computational Science de Hans Petter Langtangen chez Springer