Python interface ignores numpy strides

[Neilen Marais]

In the process of solving a complex valued problem, I ended up with a result in a numpy complex arrays. Trying to evaluate the resulting function gave me nonsense results, even though the matrix quantities seemed to be OK. It turns out that the dolfin C++/python interface seems to ignore numpy strides. An example:

V = dolfin.FunctionSpace(mesh, "Nedelec 1st kind H(curl)", order)
# x = numpy.complex128 complex result vector
u_re.vector()[:] = numpy.real(x)
u_im.vector()[:] = numpy.imag(x)
eval_point = [0,0,0]
# result1 will be junk
result1 = u_re(eval_pt) + 1j*(eval_pt)
# Copying ensure contiguous
u_re.vector()[:] = numpy.require(numpy.real(x), requirements='C')
u_im.vector()[:] = numpy.require(numpy.imag(x), requirements='C')
# result2 should be good
result2 = u_re(eval_pt) + 1j*(eval_pt)

The problem can be seen by looking at the strides:

numpy.real(x).strides -> (16,)
numpy.require(numpy.imag(x), requirements='C').strides -> (8,)

Suggested fix:

The dolfin C++ wrappers should do one of the following when a non-contiguous array is encountered

1) Use the strides info to properly handle strided data without copies
2) Raise an error, making it the caller's responsibility to ensure it is contiguous (using numpy.require())
3) Use the C API call PyArray_ContiguousFromAny() to ensure a contiguous array in the wrapper

Not sure which approach is preferable. 1) seems to be the "right" way, while 2) or 3) will be easier to implement. 2) lets the user of dolfin know what's going on so that they aren't surprised when copies are made behind their back, while 3) is more convenient.

[Johan Hake (johan-hake)]

Neilen!

Thanks for spotting this. Unfortunately we disregards any information about strides, and we silently assume contiguous arrays...

I think you suggestions are well good, and I think a mix of 1) and 2) will be the best solution.

Some places where NumPy arrays are used we iterate over the values and grab them. In this user case we can use the strides. This situation goes for all std::vector typemaps. However other times we really need a contiguous array, as the pointer to the data is used directly. I think the best thing is to raise an error if the array is not contiguous.

Johan