Ctypes
======
Ctypes is a very interesting python package which lets you import
shared object libraries into python and call them directly. I
should say that even though this is called ctypes, it can be used
just as well to call functions from libraries written in fortran.
The only complication is you need to know what a fortran function
looks like to C. This is simple everything is a pointer, so if your
fortran function would be called as foo(A,N) where A is an array
and N is its length, then to call it from C it takes a pointer to
an array of doubles and a pointer to an int. The other thing to be
aware of is that from C, fortran functions usually have an
underscore appended. That is, a fortran function foo would appear
as foo
from C (this is usually the case but is compiler dependent). Having
said this, the following examples are in C.
As an example suppose you write the following simple C program
.. code-block:: c
#include
int sum(double *x,int n)
{
int i;
double counter;
counter = 0;
for(i=0;i
#include
double timestep(double *u,int nx,int ny,double dx,double dy)
{
double tmp, err, diff,dx2,dy2,dnr_inv;
dx2=dx*dx;
dy2=dy*dy;
dnr_inv=0.5/(dx2+dy2);
err = 0.0;
int i,j;
for (i=1; i 1e-6)
{
err=timestep(u,nx,ny,dx,dy);
iter++;
}
return err;
}
We can compile it by running at the command line
::
gcc -c laplace.c
gcc -shared -o laplace.so laplace.o
Now in sage (notebook or command line) execute
::
from ctypes import *
laplace=CDLL('/home/jkantor/laplace.so')
laplace.timestep.restype=c_double
laplace.solve_in_C.restype=c_double
import numpy
u=numpy.zeros((51,51),dtype=float)
pi_c=float(pi)
x=numpy.arange(0,pi_c+pi_c/50,pi_c/50,dtype=float)
u[0,:]=numpy.sin(x)
u[50,:]=numpy.sin(x)
def solve(u):
iter =0
err = 2
n=c_int(int(51))
pi_c=float(pi/50)
dx=c_double(pi_c)
while(iter <5000 and err>1e-6):
err=laplace.timestep(u.ctypes.data_as(c_void_p),n,n,dx,dx)
iter+=1
if(iter %50==0):
print((err,iter))
return (u,err,iter)
Note the line laplace.timestep.restype=c
double. By default ctypes assumes the return values are ints. If
they are not you need to tell it by setting restype to the correct
return type. If you execute the above code, then solve(u) will
solve the system. It is comparable to the weave or fortran
solutions taking around .2 seconds. Alternatively you could do
::
n=c_int(int(51))
dx=c_double(float(pi/50))
laplace.solve_in_C(n.ctypes.data_as(c_void_p),n,n,dx,dx)
which computes the solution entirely in C. This is very fast.
Admittedly we could have had our fortran or weave routines do the
entire solution at the C/Fortran level and we would have the same
speed.
As I said earlier you can just as easily call a shared object
library that is written in Fortran using ctypes. The key point is
it must be a shared object library and all fortran arguments are
passed by reference, that is as pointers or using byref. Also even
though we used very simple data types, it is possible to deal with
more complicated C structures. For this and more about ctypes see
http://python.net/crew/theller/ctypes/