[Numpy-discussion] Unnecessarily bad performance of elementwise operators with Fortran-arrays

[Sebastian Haase]

Hi,
this comment might be counterproductive to this discussion:
However.
I'm also using numpy as a basis for putting together an image analysis
"platform"/library.
This includes both 2D and higher dimensional images.
Since all my code, if not n Python, is written in C or C++ and not
Fortran, I decided early on that I had to get used to "invese
indexing", as in
image[y,x] or image[z,y,x]

"Inverse" only refers to the common (alphabetecally oriented)
intuition to say "at point x,y,z" rather than "z,y,x".
I also argueed that mathematical matrices (mostly / often) use an
index order as "row , column" which essentially is "y, x" (if the
matrix is seen as an image)

In any case, getting used to a "z,y,x" order saved me lots of
technical troubles, and after all, most code is or will be C/C++ and
not Fortran.

Hope I was not entirey misguided,
Sebastian Haase



On Nov 9, 2007 12:37 PM, Hans Meine <meine at informatik.uni-hamburg.de> wrote:
> Am Freitag, 09. November 2007 00:16:12 schrieb Travis E. Oliphant:
> > C-order is "special" in NumPy due to the history.  I agree that it
> > doesn't need to be and we have taken significant steps in that
> > direction.
>
> Thanks for this clarifying statement.
>
> > Right now, the fundamental problem is probably due to the
> > fact that the output arrays are being created as C-order arrays when the
> > input is a Fortran order array.  Once that is changed then we can cause
> > Fortran-order inputs to use the simplest path through the code.
>
> That looks like a plan.  I have read about __array_wrap__ in your book, and it
> sounds exactly like what we need for our imaging library; this way we can
> ensure that the result of operations on images is again an image.  Here, it
> is crucial that the resulting image has Fortran order, too (for C/C++
> algorithms to agree with the Python idea of the images' orientation).
>
> > Fortran order arrays can be preserved but it takes a little extra work
> > because backward compatible expectations had to be met.  See for example
> > the order argument to the copy method of arrays.
>
> What do you mean exactly (if you have something specific in mind at all)?
> Should the order be configurable for all operations, and default to "C" for
> backwards compatibility?
> (BTW: copy() has no arguments yet in my numpybook, page 57.)
>
> At last, I had a look at the relevant code for adding arrays; this is what I
> found: In core/src/umathmodule.c.src, the code templates for the relevant
> ufuncs (e.g. add) are defined, which will get expanded to e.g. DOUBLE_add.
> This will get wrapped into an ufunc in 'InitOperators' using
>
>  f = PyUFunc_FromFuncAndData(add_functions, add_data, add_signatures, 18,
>                              2, 1, PyUFunc_Zero, "add",
>                              "adds the arguments elementwise.", 0);
>
> and inserted as "add" into a dictionary which is later passed to
> PyArray_SetNumericOps.  This will install the ufunc in the static
> NumericOps struct 'n_ops', which is used in the 'array_add' function
> installed in the PyTypeObject.
>
> Thus, it is the ufunc code (not suprisingly), which seems to call
> DOUBLE_add for the inner loop:
>
> | static void
> | DOUBLE_add(char **args, intp *dimensions, intp *steps, void *func)
> | {
> |     register intp i;
> |     intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
> |     char *i1=args[0], *i2=args[1], *op=args[2];
> |     for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
> |         *((double *)op)=*((double *)i1) + *((double *)i2);
> |     }
> | }
>
> If I understood David correctly, he suggested an unstrided variant of
> this inner loop, which simply uses pointer increments.  While this is
> a good idea (also probably quite some work), the real thing bugging me is
> that the above DOUBLE_add could (and should!) be called by the ufunc
> framework in such a way that it is equally efficient for C and Fortran
> arrays.  (I.e. with steps set to sizeof(double), without prior
> copying/conversion.)  It could even be used with negative step sizes when the
> input and output arrays have different orders.  Disclaimer: So far, I did not
> look at the ufunc code yet, so I do not know which code paths exist.
>
> BTW: What are the "void *func" in core/src/umathmodule.c.src arguments
> for?  Shouldn't that be "void * /*func*/", since the parameters are
> never used?
>
> --
> Ciao, /  /
>     /--/
>    /  / ANS
>
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion at scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>