0

I would like the create an fortran-contiguous view of a 3-D C-contiguous array named a. For a 2-D matrix, this is as simple as a.T, but transpose is not defined for 3-D arrays. Happily the reshape function takes an order argument, which allows the output array to be given F-contiguous flag.

Unfortunately, naively reshaping the array seem to be very slow unless I ravel() it first.

What is going on with the speed difference? Is something else happening?

Is there a better way to view the array in fortran order without making a copy (asfortranarray makes a copy in order to preserve the shape).

Ultimately I'm trying to provide this array to a Fortran routine using f2py. Is there an easy way to suggest f2py should auto-transpose the array?

In [1]: a = np.empty([1000,2000,3000],dtype=np.float32)

In [2]: %timeit b = np.reshape(a.ravel(),[3000,2000,1000],order='F')
1.43 µs ± 5.42 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

In [3]: %timeit c = np.reshape(a,[3000,2000,1000],order='F')
24 s ± 46.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Both the b and c matrix seem to be equivalent views of a

Update: As @hpaulj points out, the c version is slower because it makes a copy. Relevant documentation:

The order keyword gives the index ordering both for fetching the values from a, and then placing the values into the output array

Which is why the c array is not the same ordering, and is a copy.

Additionally transpose (a.T) is defined on n-D arrays, so that's the simplest solution.

L. Robison
  • 151
  • 1
  • 9
  • 2
    The order of elements in `b` and `c` is not the same. `b` is a `view` of `a`; `c` is not. – hpaulj Dec 02 '21 at 06:37
  • The expensive cost comes from the implicit required transposition. See [How to avoid huge overhead of single-threaded NumPy's transpose?](https://stackoverflow.com/questions/67431966) for more information. Still, transposing a 24 GiB array will take some time (a good server machine can do that in less than a second with a good algorithm)... – Jérôme Richard Dec 02 '21 at 10:30
  • @hpaulj Hm, you are right, `c` is not a `view` of `a`, verified by changing values in `a` not reflected in `c`. Somewhat puzzling because I didn't realize reshape would ever make a copy, and also `c` has the `OWNDATA` flag as `False`. – L. Robison Dec 02 '21 at 18:07
  • @JérômeRichard I disagree. There is no intent to transpose happening here, I'm simply asking Python to provide the same data as-is to a Fortran function, which requires python to consider the dimension ordering to be reversed without moving data. – L. Robison Dec 02 '21 at 18:10

0 Answers0