[SciPy-Dev] f2py, the fortran integer type, and npy_intp

[Kurt Smith]

On Sat, Jul 10, 2010 at 12:31 PM, Charles R Harris
<charlesr.harris@gmail.com> wrote:
> I note that some c programs try to call f2py generated fortran interfaces
> using the npy_intp type for the array dimension. There are two problems
> here, first the prototype for the generated wrapper uses int causing
> compiler warnings on a 64 bit os, and second, the fortran subroutines
> themselves use integer types. So some questions. What precision are fortran
> integers? Should we be using npy_intp for array dimensions at all? Is there
> any transparent way to have 64 bit support?

I can comment on some of the above (not the f2py-specific bits, though).

1) "What precision are fortran integers?"

According to a somewhat authoritative interlanguage programming site
[0], a fortran 'integer' corresponds to a C 'int'.  This can be relied
upon, but other fortran type declarations cannot (e.g. 'integer*8' !=
a C 'long' on all platforms, with all compilers).  If sizeof(npy_intp)
== sizeof(int), then everything will be fine, but it's a near
certainty there are platforms where this doesn't hold, and may lead to
problems if array sizes are large.

[0] http://www.math.utah.edu/software/c-with-fortran.html

2) "Should we be using npy_intp for array dimensions at all?"

To interoperate with PEP-3118 buffers, the short answer is "yes."

As I'm sure you know, the PyArrayObject struct uses npy_intp's for the
dimensions and the strides.  Further, the PEP-3118 buffer protocol
calls for Py_ssize_t (which maps to npy_intp inside numpy) for the
"shape", "stride" and "suboffset" arrays inside the Py_buffer struct.
As I understand it, this is primarily required for when "suboffset" is
>=0 in a specific dimension, since it indicates the buffer is arranged
in a non-contiguous-memory layout and pointer arithmetic is required
to access the array elements.

Here's the relevant PEP-3118 section:

"""
suboffsets

    address of a Py_ssize_t * variable that will be filled with a
pointer to an array of Py_ssize_t of length *ndims. If these suboffset
numbers are >=0, then the value stored along the indicated dimension
is a pointer and the suboffset value dictates how many bytes to add to
the pointer after de-referencing. A suboffset value that it negative
indicates that no de-referencing should occur (striding in a
contiguous memory block). If all suboffsets are negative (i.e. no
de-referencing is needed, then this must be NULL (the default value).
If this is not requested by the caller (PyBUF_INDIRECT is not set),
then this should be set to NULL or an PyExc_BufferError raised if this
is not possible.

    For clarity, here is a function that returns a pointer to the
element in an N-D array pointed to by an N-dimesional index when there
are both non-NULL strides and suboffsets:

    void *get_item_pointer(int ndim, void *buf, Py_ssize_t *strides,
                           Py_ssize_t *suboffsets, Py_ssize_t *indices) {
        char *pointer = (char*)buf;
        int i;
        for (i = 0; i < ndim; i++) {
            pointer += strides[i] * indices[i];
            if (suboffsets[i] >=0 ) {
                pointer = *((char**)pointer) + suboffsets[i];
            }
        }
        return (void*)pointer;
    }

    Notice the suboffset is added "after" the dereferencing occurs.
Thus slicing in the ith dimension would add to the suboffsets in the
(i-1)st dimension. Slicing in the first dimension would change the
location of the starting pointer directly (i.e. buf would be
modified).
"""

Hope this gives you some info to go on.

Kurt

>
> Example code
>
> (from scipy/interpolate/src/__fitpack.h)
>
> void
> SPLDER(double*,int*,double*,int*,int*,double*,double*,int*,int*,double*,int*);
>                                                            ^^^^
> void SPLEV(double*,int*,double*,int*,double*,double*,int*,int*,int*);
>                                                           ^^^^
>
> static char doc_spl_[] = " [y,ier] = _spl_(x,nu,t,c,k,e)";
> static PyObject *fitpack_spl_(PyObject *dummy, PyObject *args)
> {
>     int n, nu, ier, k, e=0;
>     npy_intp m;
>
>     ...
>
>     if (nu) {
>         SPLDER(t, &n, c, &k, &nu, x, y, &m, &e, wrk, &ier);
>                                         ^^
>     }
>     else {
>         SPLEV(t, &n, c, &k, x, y, &m, &e, &ier);
>
> ^^
>     }
>  ...
> }
>
> (from scipy/interpolate/src/fitpack.pyf
>
>      subroutine splev(t,n,c,k,x,y,m,e,ier)
>        ! y = splev(t,c,k,x,[e])
>        real*8 dimension(n),intent(in) :: t
>        integer intent(hide),depend(t) :: n=len(t)
>        real*8 dimension(n),depend(n,k),check(len(c)==n),intent(in) :: c
>        integer :: k
>        real*8 dimension(m),intent(in) :: x
>        real*8 dimension(m),depend(m),intent(out) :: y
>        integer intent(hide),depend(x) :: m=len(x)
>        integer check(0<=e && e<=2) :: e=0
>        integer intent(hide) :: ier
>      end subroutine splev
>
> (from scipy/interpolate/fitpack/splev.f)
>
>       subroutine splev(t,n,c,k,x,y,m,e,ier)
>
> ...
>
> c  ..scalar arguments..
>       integer n, k, m, e, ier
>
> ...
>
>
> Note also that the checks generated by f2py -- and they are generated --
> don't seem to work.
>
> Chuck
>
>
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
>