9

I have an array which is read from a fortran subroutine as a 1D array via f2py. Then in python, that array gets reshaped:

a=np.zeros(nx*ny*nz)
read_fortran_array(a)
a=a.reshape(nz,ny,nx)  #in fortran, the order is a(nx,ny,nz), C/Python it is reversed

Now I would like to pass that array back to fortran as a 3D array. 

some_data=fortran_routine(a)

The problem is that f2py keeps trying to transpose a before passing to fortran_routine. fortran routine looks like:

subroutine fortran_routine(nx,ny,nz,a,b)
real a
real b
integer nx,ny,nz
!f2py intent(hidden) nx,ny,nz
!f2py intent(in) a
!f2py intent(out) b
...
end subroutine

How do I prevent all the transposing back and forth? (I'm entirely happy to use the different array indexing conventions in the two languages).

EDIT

It seems that np.asfortranarray or np.flags.f_contiguous should have some part in the solution, I just can't seem to figure out what part that is (or maybe a ravel followed by a reshape(shape,order='F')?

EDIT

It seems this post has caused some confusion. The problem here is that f2py attempts to preserve the indexing scheme instead of the memory layout. So, If I have a numpy array (in C order) with shape (nz, ny, nx), then f2py tries to make the array have shape (nz, ny, nx) in fortran too. If f2py were preserving memory layout, the array would have shape (nz, ny, nx) in python and (nx, ny ,nz) in fortran. I want to preserve the memory layout.

Chiel
  • 6,006
  • 2
  • 32
  • 57
mgilson
  • 300,191
  • 65
  • 633
  • 696
  • hi mgilson, a quick question regarding this very problem. I've written a fortran code that takes a 3d array as input: – toylas Nov 04 '15 at 02:17

2 Answers2

4

Fortran doesn't reverse the axis order, it simply stores the data in memory differently from C/Python. You can tell numpy to store the data in Fortran order which is not the same as reversing the axes.

I would rewrite your code as this

a=np.zeros(nx*ny*nz)
read_fortran_array(a)
a=a.reshape(nx,ny,nz, order='F') # It is now in Fortran order

Now, f2py won't try to reorder the array when passing.

As a side note, this will also work

a=a.reshape(nx,ny,nz) # Store in C order

because behind the scenes, f2py performs these operations when you pass a C-order array to a Fortran routine:

a=a.flatten() # Flatten array (Make 1-D)
a=a.reshape(nx,ny,nz, order='F')  # Place into Fortran order

But of course it is more efficient to store in Fortran order from the beginning.

In general, you shouldn't have to worry about the array ordering unless you have a performance-critical section because f2py takes care of this for you.

SethMMorton
  • 45,752
  • 12
  • 65
  • 86
  • Even simpler would be to pass the order='F' paramter to the initial call to numpy.zeros. Then no reshaping or reordering should be required. – DaveP Jun 28 '12 at 22:55
  • I'm sorry about the slow feedback on this one (I was out of town). I really don't understand your solution here. `a=a.reshape(nx,ny,nz)` is *not what I want* since `nx` is supposed to be the fast index (here it is the slow index). I can't use the first solution since I've written a lot more code assuming C/Python indexing on the python side. I tried `b=a.flatten().reshape(tuple(reversed(a.shape)),order='F')`, but that didn't work (Error: `0-th dimension must be fixed to 0 but got 448`). – mgilson Jul 02 '12 at 00:36
  • 1
    sorry. ID10T error. PEBKAC, whatever. `b=a.flatten().reshape(tuple(reversed(a.shape)),order='F')` does work (Thanks for the idea). Efficiency here shouldn't be an issue since numpy isn't making a copy of the array data, only a change to the meta-data. The `reversed` does need to be there since I'm using C/Python indexing in the C/Python part and Fortran indexing in the fortran part of the code. It makes the most sense (to me) that way. (If you update your post to include this, I'll happily accept it your answer). – mgilson Jul 02 '12 at 00:56
  • Final update. It appears the `b=np.ravel(a).reshape(tuple(reversed(a.shape)),order='F')` works. `a.flatten()` seems to return a new array. – mgilson Jul 02 '12 at 12:17
  • @mgilson Reguarding fast indexing, in Fortran the first elements should be the fastest, so `a=a.reshape(nx,ny,nz)` would make `nx` the fastest index in Fortran. Also, don't confuse array storage order in memory with array indexing. – SethMMorton Jul 03 '12 at 03:14
  • @KalAurum -- what do you mean "don't confuse array storage order in memory with array indexing" -- Those two things are tied together. In a computer, an array is simply a sequence of memory addresses. the indexing is just a way of pointing to the index you want. In fortran `a(ix,iy,iz)` points to the same memory element that `a[iz+1,iy+1,ix+1]` points to in python -- If you keep that straight, there is no need for `f2py` to create temporary arrays, or to specify any special memory ordering, etc. Ultimately, it's all just a 1D array anyway... – mgilson Jul 03 '12 at 11:56
  • @mgilson I should have said array shape, not indexing. What I mean is that a numpy array is accessed by the user the same way whether it is stored in C order or Fortran order even though the data is laid out differently, and reversing the shape (transposing) is not the same as converting from C order to Fortran order. – SethMMorton Jul 03 '12 at 12:14
  • @mgilson Also, even though for the user `a(ix,iy,iz)` (in Fortran) and `a[iz,iy,ix]` (in python) access the same element (note that you shouldn't have to offset your indexing within python), they point to the same location in the array **only** if the python array is stored in Fortran order. Otherwise they actually point to different memory locations. [This link](http://and-what-happened.blogspot.com/2011/05/row-major-and-column-major-matrices.html) explains what I am talking about in more detail than a comment allows. – SethMMorton Jul 03 '12 at 12:23
  • You do need to add the offset -- unless you play tricks on the fortran side (fortran is 1 based, numpy is 0 based -- even with fortran order), and they point to the same element **only** if the python array is stored in "C" order -- which is the default, natural way to index an array in python. – mgilson Jul 03 '12 at 12:25
  • My solution works because I have a C_CONTIGUOUS array (`A`) of dimensions nz,ny,nx -- (native python indexing). When you use the transpose (`A.T`), you get a view. The memory layout is still the same, but the view has dimensions nx,ny,nz and is marked as F_CONTIGUOUS. When f2py sees the transposed array and looks at it's meta-data, it sees a fortran array that has shape nx,ny,nz -- exactly what it wants. The confusion might be because when I posted the question, I thought transposing the array returned a copy and so I was trying to avoid having f2py transpose for me. – mgilson Jul 03 '12 at 12:33
  • A agree with the first statement, I logged back on to change that (looks like I can't now, but I should have written the Fortran array as `a(ix+1,iy+1,iz+1)` to account for the 1 offset). However, did you read the link I posted? Also, check out [the Wikipedia entry](http://en.wikipedia.org/wiki/Row-major_order) for a more visual explanation. – SethMMorton Jul 03 '12 at 12:34
  • let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/13354/discussion-between-mgilson-and-kalaurum) – mgilson Jul 03 '12 at 12:34
3

It looks like the answer is reasonably simple:

b=np.ravel(a).reshape(tuple(reversed(a.shape)),order='F')

works, but apparently, this is the same thing as:

b=a.T

since transpose returns a view and a quick look at b.flags compared with a.flags shows that this is what I want. (b.flags is F_CONTIGUOUS).

mgilson
  • 300,191
  • 65
  • 633
  • 696
  • Incredibly easy solution for what seemed to be a very complex problem. Great answer. – Chiel Oct 01 '17 at 11:37
  • Yes, .T is the right answer! Because .T returns a view of the unaltered data but with the indexing flipped to Fortran order (F_CONTIGUOUS), so f2py leaves the array alone, which is what you want. – Jonatan Öström Mar 02 '22 at 09:29