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.
_______________________________________________
Mailing list: https://launchpad.net/~dolfin
Post to : dolfin@lists.launchpad.net
Unsubscribe : https://launchpad.net/~dolfin
More help : https://help.launchpad.net/ListHelp
