f2py — prevent array reordering

Question

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.

Answer 1

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).

Answer 2

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.