dolfin team mailing list archive

[Johan Hake <783380@xxxxxxxxxxxxxxxxxx>]

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

** Changed in: dolfin
       Status: New => Confirmed

** Changed in: dolfin
   Importance: Undecided => High

** Changed in: dolfin
     Assignee: (unassigned) => Johan Hake (johan-hake)

-- 
You received this bug notification because you are a member of DOLFIN
Team, which is subscribed to DOLFIN.
https://bugs.launchpad.net/bugs/783380

Title:
  Python interface ignores numpy strides

Status in DOLFIN:
  Confirmed

Bug description:
  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.