F2PY: Passing single-element arrays to Fortran

The following Fortran code fills a 2D array x with value v

      subroutine fill(x,v,m,n)
      real*8 x(m,n),v
      integer m,n,i
cf2py intent(in) :: x,v,m,n
      forall(i=1:m,j开发者_开发问答=1:n) x(i,j) = v
      end

When calling this function from Python:

x = numpy.array([[[0.0]]],order='F')
fill(x[:,:,0],2.0)
assert(x[0,0,0]==2.0) # Assertion failed

Why is this assertion failing ?

x should be declared as intent(inout) if you want it to pass values back to the caller.

However this causes an additional problem because passing array slices doesn't work for intent(inout) arrays. In this simple example you can get around it by calling from python:

fill(x, 2.0).

If you really wanted to pass a slice then you need to declare x as intent(in,out), and call from python: x[:,:,0] = fill(x[:,:,0],2.0)

The description of the different attributes can be found at:

http://cens.ioc.ee/projects/f2py2e/usersguide/index.html#attributes

I just had this issue. Neither intent(inout) nor intent(inplace) fixed it. The problem is obviously in the array checking routine array_from_pyobj() in fortranobject.c, which ships with f2py and is linked with every module one builds. array_from_pyobj() makes every effort that any input is converted to a properly shaped contiguous array, doing many many checks. One of them does not work properly with single-element arrays, such that a copy is made instead of working on the original array.

One could fix that, but ...well... I don't want any of this polymorphism stuff in an interface to a performance library anyway... I have a Python class wrapping the library calls that already guarantees that all the arguments are passed on correctly.

So my fix was to make my own copy of fortranobject.c and simply replace array_from_pyobj() by the following dummy version:

extern PyArrayObject * 
array_from_pyobj(const int type_num, npy_intp *dims,
                 const int rank, const int intent, PyObject *obj) {

    int i;
    PyArrayObject *arr = (PyArrayObject *)obj;

    for (i=0; i<arr->nd; ++i)
        dims[i] = arr->dimensions[i];
    return arr;
}

My plan is to use the original fortranobject.c during development of my Python wrapper class, where admittedly it is an advantage to get gentle error messages instead of hard crashes every time one makes a mistake calling a library function. Once I'm sure that all the library calls work, I'll use my custom fortranobject.c, which also works with single-element arrays.