[Matrix-SIG] Re: FPIG: advice needed

[Travis Oliphant]

> 
> 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.
> 
> 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)
> 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.

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.


> 
> 2) More secure would be, if the result of f2py is python function having
> form:
>   n,x = _f(n,x)
> that is, if needed, type conversation can be implemented easily and in a
> general way. Disadvantage is that using it in a python, results a little
> longer code. For example, dswap from the blas should be always called as
> follows:
>   n,x,incx,y,incy = blas._dswap(n,x,incx,y,incy)
> instead of current shorter version:
>   blas._dswap(len(x),x,1,y,1)

I don't understand why the "longer-code" here.  Maybe I don't understand
the point you are making.  But, why couldn't you have the same calling
syntax as 1.

> Final implementation of 1) may be complex (I don't see right know how to
> avoid its problems. These could be ignored but then only correct
> arguments should be assumed). Implementation of 2) means more code but it
> seems easier to carry out and should be safer. But performance will
> probably be lower, I think.
> 
> What do you think? Which ways would you prefer as a f2py.py user?
> Or should be both ways present trough an optional key in f2py.py?

I'm not sure exactly what the concern is.  If the concern is over type
casting, then we should use PyArray_ContiguousFromObject()

Otherwise I'm not sure what the concern is here.