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

[Charles R Harris]

On Sat, Jul 10, 2010 at 12:30 PM, Kurt Smith <kwmsmith at gmail.com> wrote:

> On Sat, Jul 10, 2010 at 12:31 PM, Charles R Harris
> <charlesr.harris at 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
>
>
Since 'int' is 32 bits with gcc and msvc regardless of whether the operating
system is 32 or 64 bits, we are going to have a problem with all the fortran
code in scipy that uses integer types for array indexing. That is the case
for all the code in fitpack, for instance. The only way I see to get around
that is to have two fortran libraries, one modified for 64 bits and another
for 32 bits. How does your cythonized fortran deal with this?


> 2) "Should we be using npy_intp for array dimensions at all?"
>
> To interoperate with PEP-3118 buffers, the short answer is "yes."
>
>
Well, yes. The question is what to do with all that fortran code that
doesn't support it.


> 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.
>
>
Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20100710/5864e728/attachment.html>