simplified cython wrappers to fortran
Ondrej Certik
I have figured out how to use Cython directly to wrap almost any
fortran code, so I am posting it here to get some feedback.
You don't need any fortran interface nor C header files, just
cythonize this single file (at the end of my email), compile and link
with the fortran .o, that's it. Very slick, as you can see, it wraps
any fortran array (float double/single, int) as long as it's not
assume shape (assume size is ok), it works for in/doubles and also for
strings.
Obviously, all of these things are compiler dependent, but this simple
cython script essentially contains the user's definitions (types
mappings), that at some point have to be there anyway (in fwrap). In
my case, I usually just need to wrap couple functions, so it's not a
big deal to fix this one cython file at different platforms by hand.
For larger projects it will need some automation (fwrap), that's
clear. Also I have thorough tests to test all kinds of arguments to
make sure the wrappers are right, full code is here:
http://github.com/certik/qsnake/tree/master/fortran_test/
Another tricky part was how to call Python from fortran, so that I can
construct numpy arrays on the fly (without a-priori knowing their
size!). To do that, I just call this from fortran:
integer(kind=8) :: ctx
call push_array(ctx, "lam"//CHAR(0), lam, size(lam))
and here is the definition of "push_array":
cdef extern void sle1d_(long *)
cdef extern void push_array_(long *ctx, char *name, double *lam,
int *lam_len, int name_len):
"""
Creates a numpy array from "lam" and copies it to the namespace in "ctx".
"""
from numpy import array
cdef object namespace = <object>(<void *>ctx)
cdef ndarray _lam = array([1, 2], dtype="f8")
cdef n = lam_len[0]
memcpy(_lam.data, lam, sizeof(double)*n)
namespace[name] = _lam
def sle1d():
from numpy import array
cdef object namespace = {}
sle1d_(<long *>(<void *>namespace))
return namespace
So the user calls sle1d(), which constructs a python namespace, passes
the pointer ctx to fortran's sle1d(ctx), the pointer has to be typed
to a long (if anyone knows a better solution, let me know), and in the
push_array, I just recover the python namespace and push this in. So
the workflow is:
* write a simple fortran driver for the current fortran program, using
the program's (library) functionality, but don't read/write any files
* accept input as arguments
* construct output using push_array to dynamically construct numpy arrays
* in the cython wrappers, construct this namespace, fill it with
push_array, and return it to the user
Here is how it looks like in practice:
In [1]: import sle1d
In [2]: d = sle1d.sle1d()
Eigenvector orthogonality:
c_i Am c_j:
1 1: 7.73967964260E-01
1 2: 6.32011804003E-14
2 2: 3.15490838047E+00
c_i Bm c_j:
1 1: 1.00000000000E+00
1 2: 1.64798730218E-16
2 2: 1.00000000000E+00
Eigenvector residuals ||Am c - lam Bm c||/||Am c||:
1: 1.90903827911E-12
2: 6.48740175577E-13
In [3]: d
Out[3]: {'lam': array([ 0.77396796, 3.15490838])}
In [4]: lam = d["lam"]
In [5]: lam
Out[5]: array([ 0.77396796, 3.15490838])
In [6]: lam[0]
Out[6]: 0.77396796426050707
In [7]: lam[1]
Out[7]: 3.1549083804681266
I am posting it here, because it was all a bit tricky to get
everything right (string args, pointer handling in fortran, ...).
Future improvements:
* improve push_array to accept fortran-formatted strings so that in
fortran, I can just do "call push_array(ctx, "lam", lam, size(lam))",
this will then look like a native fortran subroutine
* possibly figure out a better handling of the "ctx" pointer. But I
actually really like that it's just a fortran integer(kind=8), because
then in fortran I really don't have to worry about this at all and
fortran uses integers to handle writing/reading to files anyway, like
in:
open(1,file='input.in',status='old')
read(1,*) ! header
so this just follows a fortran way of doing things.
* I recently wrote a Python() C++ class to handle nice and easy way to
call Python from C++:
Python *p = new Python();
p->push("i", c2py_int(5));
p->exec("i = i*2");
int i = py2c_int(p->pull("i"));
_assert(i == 10);
delete p;
(which automatically handles creating/deletion of the Python
interpreter, reference counting, ..., simply everything)
so it seems to me that I can write the exact same thing for fortran, e.g.:
integer(kind=8) ctx
integer i
ctx = new_Python()
push(ctx, "i", f2py_int(5))
exec(ctx, "i = i*2")
i = py2f_int(pull(ctx, "i"))
delete_Python(ctx)
and you can hook this up into any fortran program, push some things
into the Python namespace, launch ipython (let's say) in the middle of
some calculation and play with things a bit. I am very excited about
this and especially that Cython can directly be used to wrap fortran.
Ondrej
-----------------
cdef extern from "arrayobject.h":
ctypedef int intp
ctypedef extern class numpy.ndarray [object PyArrayObject]:
cdef char *data
cdef int nd
cdef intp *dimensions
cdef intp *strides
cdef int flags
cdef extern void args_subr1_(signed char *, short *, int *, long *, long *)
cdef extern long args_func1_(signed char *, short *, int *, long *)
cdef extern int int_arg_(int *, int *, int *)
cdef extern float single_arg_(float *, float *, float *)
cdef extern double double_arg_(double *, double *, double *)
cdef extern void string_arg1_(char *, int)
cdef extern void string_arg2_(char *, char *, char *, int, int, int)
cdef extern void string_arg3_(char *, char *, int *, char *, int *,
int, int, int)
cdef extern void single_array1_(float *)
cdef extern void double_array1_(double *)
cdef extern void double_array2_(double *, int *)
cdef extern void int4_array1_(int *)
cdef extern void int8_array1_(long *)
def args_subr1(signed char a, short b, int c, long d):
cdef long o
args_subr1_(&a, &b, &c, &d, &o)
return o
def args_func1(signed char a, short b, int c, long d):
return args_func1_(&a, &b, &c, &d)
def int_arg(int a, int b, int c):
return int_arg_(&a, &b, &c)
def float_arg(float a, float b, float c):
return single_arg_(&a, &b, &c)
def double_arg(double a, double b, double c):
return double_arg_(&a, &b, &c)
def string_arg1(char *a):
string_arg1_(a, len(a))
return a
def string_arg2(char *a, char *b, char *c):
string_arg2_(a, b, c, len(a), len(b), len(c))
return a, b, c
def string_arg3(char *a, char *b, int n1, char *c, int n2):
string_arg3_(a, b, &n1, c, &n2, len(a), len(b), len(c))
return a, b, n1, c, n2
def float_array1(ndarray a):
single_array1_(<float *>(a.data))
def double_array1(ndarray a):
double_array1_(<double *>(a.data))
def double_array2(ndarray a):
cdef int n = len(a)
double_array2_(<double *>(a.data), &n)
def int4_array1(ndarray a):
int4_array1_(<int *>(a.data))
def int8_array1(ndarray a):
int8_array1_(<long *>(a.data))
Kurt Smith
> Hi,
>
> I have figured out how to use Cython directly to wrap almost any
> fortran code, so I am posting it here to get some feedback.
>
> You don't need any fortran interface nor C header files, just
> cythonize this single file (at the end of my email), compile and link
> with the fortran .o, that's it. Very slick, as you can see, it wraps
> any fortran array (float double/single, int) as long as it's not
> assume shape (assume size is ok), it works for in/doubles and also for
> strings.
This is great. Very nice.
One benefit fwrap will allow is total portability -- between different
platforms & different fortran compilers. All the name mangling
issues, datatype sizes, etc. is handled for you. What you have here
is very nice.
Very nice. I like how it parallels the functionality in the C++ & Python stuff.
>
> So the user calls sle1d(), which constructs a python namespace, passes
> the pointer ctx to fortran's sle1d(ctx), the pointer has to be typed
> to a long (if anyone knows a better solution, let me know), and in the
> push_array, I just recover the python namespace and push this in. So
> the workflow is:
>
There is a 'c_ptr' type provided in the ISO C BINDING intrinsic module
(all modern fortran compilers have this module; gfortran >= 4.3.3.)
You'd put this at the top of your module or procedure:
use, intrinsic :: ISO_C_BINDING
and then the ctx declaration would be:
type(c_ptr) :: ctx
You'd have to explicitly declare the push_array's interface like this:
interface
subroutine push_array(ctx, name, lam, lam_len, name_len) bind(c,
name='push_array')
use, intrinsic :: ISO_C_BINDING
implicit none
type(c_ptr), value :: ctx ! note the 'value' attribute
integer(c_int), intent(in) :: name_len
character(kind=c_char), dimension(name_len) :: name
integer(c_int), intent(in) :: lam_len
real(c_double), dimension(lam_len) :: lam
end subroutine push_array
end interface
Yes it's really verbose, but that's fortran for you. The 'bind' stuff
on the declaration line ensures that when fortran calls 'push_array'
it calls exactly that function in C -- without the underscore. This
takes care of one portability issue between fortran compilers.
You could do this easily with null-terminated strings -- just pass in
"lam"//C_NULL_CHAR as the name & it will act like a normal char *
argument.
> * possibly figure out a better handling of the "ctx" pointer. But I
> actually really like that it's just a fortran integer(kind=8), because
> then in fortran I really don't have to worry about this at all and
> fortran uses integers to handle writing/reading to files anyway, like
> in:
You'd have to use a type(c_ptr) declaration as shown above.
Thanks for posting this -- very nice.
> --
> You received this message because you are subscribed to the Google Groups "Fwrap Users" group.
> To post to this group, send email to fwrap...@googlegroups.com.
> To unsubscribe from this group, send email to fwrap-users...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/fwrap-users?hl=en.
>
>
Dag Sverre Seljebotn
> Hi,
>
> I have figured out how to use Cython directly to wrap almost any
> fortran code, so I am posting it here to get some feedback.
>
> You don't need any fortran interface nor C header files, just
> cythonize this single file (at the end of my email), compile and link
> with the fortran .o, that's it. Very slick, as you can see, it wraps
> any fortran array (float double/single, int) as long as it's not
> assume shape (assume size is ok), it works for in/doubles and also for
> strings.
>
> Obviously, all of these things are compiler dependent, but this simple
> cython script essentially contains the user's definitions (types
> mappings), that at some point have to be there anyway (in fwrap). In
>
already. That's what fwrap uses, so there's no compiler dependency as
long as the compiler supports that standard (which all major ones does
in current versions, though older compilers might get problems). And the
purpose of fwrap's wrapper Fortran modules is to wrap non-iso-c-binding
Fortran code into Fortran code using the ISO C bindings (so that
everything becomes standards compliant).
Of course, what you are doing works great (and I've been doing the same
for a library I'm using myself for some time) -- this is just to explain
the additional challenges fwrap deals with, and how the traditional
approach used by you and me is different from fwrap.
Dag Sverre