openlang
数值计算

f2py copies my C-ordered array on every call

后续讨论 3 ↓查看参考资料 ↓

Passing a NumPy array into a Fortran subroutine through f2py. The array is created the usual way, np.zeros((n, m)), so it's C-ordered. Fortran wants the columns contiguous. f2py knows that, and for an ordinary input argument it makes an element-wise copy in the right order without saying anything. No warning. It works, the numbers are right. Every call in the loop just copies the matrix.

intent(inout) is worse but more honest. A copy there would swallow the changes, so f2py demands a properly laid out array and raises an exception otherwise. The cure is order='F' at creation, or one np.asfortranarray call before the loop.

It's the old column-order story, only this time it sits where two languages meet, and the bill is a copy per call.

Annoying. I'd rather have had the warning.

后续讨论 3

顺着彼此的依据读下去

  1. One-dimensional arrays are exempt. A contiguous vector is C- and F-contiguous at once. No copy. The trouble starts at rank two.

  2. Not a Fortran person, but I'm curious how many copies the loop really makes. The f2py docs build their example with -DF2PY_REPORT_ON_ARRAY_COPY=1, and then the wrapper prints a message whenever it copies an array. Flip the order and nothing else, then count the messages again.

    A slice is worth checking too. Something like a[:, ::2] taken after np.asfortranarray isn't contiguous anymore, so my guess is the copy comes back. Better to see the messages than guess, though.

  3. 接着前面的回复 ↗

    mtsuda's flag is the right tool for this. With -DF2PY_REPORT_ON_ARRAY_COPY=1, a C-ordered array reports a copy on every call, one np.asfortranarray before the loop silences it, and the a[:, ::2] slice copies again, as predicted above.

    And it gets worse. That same slice passed to intent(inout) raises, even though the parent array is in Fortran order, because intent(inout) needs a contiguous array and a strided slice isn't one.

同一个问题的其他角度

  1. 数值计算

    Ported a short MATLAB script to R and one result came back as a matrix

    Trying to redo a short MATLAB script in R, and I've barely touched MATLAB before. Most of it came over fine, except for a sin…

    回复 7rmatlab
  2. 数值计算

    na.rm = TRUE on a group that is all NA

    Per-group totals with tapply(x, g, sum, na.rm = TRUE), and one group came back as 0. Every value in that group is NA. Not a s…

    回复 5r
  3. 数值计算

    A loop that allocated a column at a time

    Slicing in Julia copies, which I knew and then forgot inside a loop. Take a function that walks the columns of a matrix and d…

    回复 4julia