dolfin team mailing list archive

[Neilen Marais <launchpad@xxxxxxxxxxxxxxxxxxxxxx>]

Public bug reported:

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.

** Affects: dolfin
     Importance: Undecided
         Status: New

-- 
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:
  New

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.