openlang
Numerical computing

f2py copies my C-ordered array on every call

The conversation continues 3 ↓Explore the sources ↓

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.

The conversation continues 3

Follow the reasoning
  1. Permalink

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

  2. Permalink

    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.

    One more case worth checking. A slice 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. Permalink
    Following an earlier reply ↗

    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.

    One more catch. 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.

Other sides of the same question

  1. Numerical computing

    sum(X) on a matrix with one row

    A function takes a matrix of observations, rows are observations and columns are features, and normalizes each column with X…

    Replies 6matlab
  2. Numerical computing

    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…

    Replies 4r
  3. Numerical computing

    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…

    Replies 4julia