[Matrix-SIG] Re: FPIG: advice needed

Pearu Peterson pearu@ioc.ee
Tue, 13 Jul 1999 21:15:36 +0300 (EETDST)


On Tue, 13 Jul 1999, Travis Oliphant wrote:
> > 
> > Travis,
> > 
> > I have found that there are number of ways how to access fortran/c raw
> > functions using f2py.py from python. I need to choose one before I can
> > write f2py.py any further.

I have rewritten f2py.py on the part where *module.c is generated.
For those, how don't have access to CVS server (where Multipack is), I
have copied the relavant files to
  http://koer.ioc.ee/~pearu/python/f2py/

> > 
> > Say, we have a fortran subroutine:
> > subroutine f(n,x)
> > integer n
> > real*8 x(*)
> > ... (say, fill array x of length n with ones)
> > end
> > In order to access it from python, the following policies are possible:
> > 1) Currently implemented one generates the following python function:
> >   {'n':n} = _f(n,x)

Now f2py.py generates python function
	None = _f(n,x)
It means that _f will return always None on success (in unsuccess
exception is raised anyway --- not implemented yet; if f is function,
its value will be returned).
It is always assumed that x is an array (if fortran code requires that).
Period.

Now if x=Numeric.array([2,3],'d') then
n=len(x)
_f(n,x)
print x
should display [1.,1.].
Note that in this way the syntax of _f will be exactly the same as in
fortran code. It means that it is much easier to use _f as no changes is
needed in the documentation of fortran function f.

> > That is if argument is of numeric type (like 'n' is int here), its value
> > (if it is changed in f) can be accessed from the dictionary that _f
> > returns. If argument is numeric array ('x'), its contents is accessed by
> > its reference, so that _f changes x in place. But this matter makes
> > checking arguments difficult: if x is assumed to be an array, but I call
> > _f with x being, say, float, then it is not immidately clear how to make
> > the conversation.
> 
> Why does this make it difficult to check the arguments?  You have an
> "O!" in the appropriate position, if it is called with a float
> argument the routine returns with an exception.

I probably will use "O!" for argument checking and only.

> 
> Force it to be an array.  Again, this is a low-level tool.  The user will
> have to make Python wrappers (which are easier to write and maintain) if
> he/she wants adjustments to the default handling.   
> 
> > Advantage is that the interface is very direct and fast --- but again, it
> > works only if arguments are correct.
> 
> That's O.K.  I like this approach better.  I think I would prefer, though, 
> that the dictionary contain a reference to the input argument 'x' for
> completeness.  This would make it easier for the user who would know what
> to expect.
> 
> O.K. wait, I think I see the trouble.  The problem is not if x is a
> Python float but if is the wrong array type (say integer rather than
> double).
> 
> Yes, this needs to be fixed even for a low-level tool.  The initiated user
> will know that a performance hit might be taken if the routine has to copy
> an integer array to a double array but at least no one will get core-dumps
> for incorrect arrays being pased to raw fortran code that doesn't check
> bounds.
> 
> It's not hard to make sure the array is the right type.  A call
> to PyContiguous_FromObject(obj, type, min_dim, max_dim) will
> do it for you. If the array is already the right type, no copying takes
> place.   This also allows arbitrary sequences to be used as input
> arguments, which is nice.

This does not work for current implementation (I think) where f2py.
Say, if x = [2,3] (it is a list),
then calling from python
	_f(len(x),x)
first in C 
  PyArg_ParseTuple("OO",&n,&x) 
will give x containing reference to the list. 
Second a new array should be generated in C: say
  x_api=PyArray_Contiguous_FromObject(x,...);
Third, call fortran function
  f_(&(n->ob_ival),(double *)(x_api->data));
(that fulfills x_api with ones)
But now I don't see how to get the contents of x_api to x without copying:
  x=x_api; will not work (I guess).

Travis: About you concern of changing from format "iO" to "OO":
I needed that because some fortran routines may want to change, say,
flags, etc that are integers, or more generally they may change numeric
types. If I will use "iO", the code will be more complicated to get these 
changed values back (previously dictionary were returned with numeric
variables). With "OO", as it is now, it is relatively clean and simple.
For argument checking I will use "O!" construction as you suggested.

Again, f2py.py generates a very low level python functions. The user
should be responsible that input arrays are in the right length because
f2py cannot always find out the correct length of an argument array (as
in fortran some times there is statement 'REAL*8 x(*)' instead of 'REAL*8
x(n)'). Best that f2py can do, is to check only the dimensions
(1D,2D,...).
So, f2py cannot prevent code dumps in general: it is programmers
responsibility.

Pearu