Dear Colleagues,

Quote:

> >> Actually, a tool to create python wrappers for fortran libraries

> >> automatically would do a lot to popularize Python in science. Like it

> >> or not, most good scientific libraries are written in Fortran.

> Andrew> Indeed; there's also piles of code in repositories like Netlib

> Andrew> that could be quickly added to Numeric's box of tools.

> Here's an idea for doing this quick-n-dirty to start with, then more cleanly

> later on.

> 1. If it doesn't already do it, modify f2c to generate C .h files that

> match the Fortran function signatures. (I know this was done one

> upon a time and can ask the person who did it if his mods are still

> available.)

> 2. Use SWIG to wrap those generated wrappers. That should allow you to

> effectively use SWIG to wrap Fortran code without teaching SWIG how

> to parse Fortran.

> Short term, the .h files may only reflect the C signatures, thus SWIG would

> be wrapping the wrapper code. With a little effort, I think SWIG's output

> phase could be trained to realize it was actually wrapping f2c output, and

> generate wrappers for the Fortran code directly...

It is a good idea to involve f2c at the first stage.

The next step depends on what one want to do with this in Python.

Suppose, one will want to pass matrices and vectors produced by NumPy to

those subroutines.

Consider the following example:

We have a matrix A (M x N), two vectors X and Y of dimension N and a vector C

of dimension M. symbolically:

real*8 A(M,N), X(N), Y(N), C(M)

and we want to pass them to a subroutine as follows:

CALL PROCESS1(A,N,X,IFLAG,Y,M,C,IMODE)

where the variables IFLAG and IMODE control the behaviour of that 'PROCESS1'.

Now, in Python, we obtain four multiarray objects a, x, y and c which

dimensions could be extracted via NumPy C-interface functions. How can

we automatically distribute these arrays and their dimensions among the

parameters of this subroutine?

Two possible solutions:

1. When modifying f2c, make sure that arrays are distinct from scalars

when converted to 'double *'. E.g. 'Double *' for true arrays.

Then build a SWIG typemap to handle arrays in a primitive manner.

See examples in SWIG tutorial.

All the conversion from NumPy objects to those 'Double *'s is then

left to Python.

For example, a transformation of a symmetric matrix to a linearized

triangular form to be used in a symmetric matrix inverse,

is done here explicitly in Python.

It seems inefficient.

Even if not, I don't like this. :-)

2. Semi-manual clean-up with more sophisticated typemaps to deal with

multiple groups of dimensions. A primitive example is given:

------------------------------------------------------------------------

%module sum

%typemap(python, ignore) int NCOLS(int *__ncols__)

{

__ncols__ = &$target;

Quote:

}

%typemap(python, ignore) int NROWS(int *__nrows__)

{

__nrows__ = &$target;

Quote:

}

%typemap(python, in) double *MATRIX

{

/* It may be a sophisticated slice. We need a plain pointer */

PyArrayObject *data_obj = (PyArrayObject *)

PyArray_ContiguousFromObject($source, PyArray_DOUBLE, 2, 2);

*__ncols__ = data_obj->dimensions[1];

*__nrows__ = data_obj->dimensions[0];

$target = (double *)data_obj->data;

Quote:

}

/*

..........

the same for MATRIX2, NCOLS2, NROWS2, ...

to describe sets of parms with different dimensions

...........

*/

extern double sum(double *MATRIX, int NCOLS, int NROWS); /* sum all the elems */

------------------------------------------------------------------------

import sum # this is the module above

from Matrix import * # this is NumPy

def test():

a = Matrix( [ [ 100.0, 108.0, 2011.0 ],

[ 107.0, 106.0, 2081.1 ],

.............................

[ 108.0, 100.0, 2091.0 ] ])

nrows, ncols = a.shape

s0 = 0.0

rcols = range(0, ncols)

for i in range(0, nrows):

for j in rcols:

s0 = s0 + a[i,j]

# We call then directly the wrapped function:

s = sum.sum(a)

print s0, s, s0 - s

test()

------------------------------------------------------------------------

In this way, the prototype for PROCESS1 may look like:

void process1_(double *MATRIX, int NROWS, double *COLUMN, int iflag,

double *COLUMN, int NCOLS, double *ROW, int imode);

where pattern ROWi is supposed to set (, check and use) __ncols__i,

COLUMNi -- __nrows__i. (i = , 1, 2, ... say, 9)

Patterns for SYMMETRIC and DISTANCE (zero diagonal) could also be

introduced with conversion code supplied.

Finally, recall that Fortran matrices are transposed w.r.t. C ones...

(The above example was intended to a C library -- that is my case.

Note that most of really *cool* libraries are translated to C long

ago, or designed in a multilingual manner -- e.g. Numeric Recipes.)

I do not mention arrays with more dimensions -- I do not use them

practically -- and many other issues.

This example is 2 days old, I have begun to use SWIG only 3-4 days ago.

So, there may be something naive in that message.

Anyway, both Python and SWIG are brilliant things!

Thanks to their authors and contributors!

Regards

Alexander