python - f2py: evita la reorganización de la matriz

interface fortran

Question

Tengo una matriz que se lee de una subrutina Fortran como una matriz 1D a través de f2py. Luego, en Python, esa matriz se reconfigura:

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

Ahora me gustaría pasar esa matriz a Fortran como una matriz 3D.

some_data=fortran_routine(a)

El problema es que f2py sigue tratando de transponer una antes de pasar a fortran_routine. La rutina de Fortran se parece a:

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

¿Cómo evito todas las transposiciones de ida y vuelta? (Estoy feliz de utilizar las diferentes convenciones de indexación de matrices en los dos idiomas).

EDITAR

Parece que np.asfortranarray o np.flags.f_contiguous deberían tener parte en la solución, simplemente no puedo entender qué parte es (o tal vez un ravel seguido de una reshape(shape,order=''F'') ?

EDITAR

Parece que esta publicación ha causado cierta confusión. El problema aquí es que f2py intenta preservar el esquema de indexación en lugar del diseño de la memoria . Entonces, si tengo una matriz numpy (en orden C) con forma (nz, ny, nx) , entonces f2py intenta hacer que la matriz tenga forma (nz, ny, nx) en fortran. Si f2py estuviera preservando el diseño de la memoria , la matriz tendría forma (nz, ny, nx) en python y (nx, ny ,nz) en fortran. Quiero preservar el diseño de la memoria.

Answer 1

Fortran no invierte el orden del eje, simplemente almacena los datos en la memoria de forma diferente a C / Python. Puede decir numpy para almacenar los datos en orden Fortran que no es lo mismo que invertir los ejes.

Reescribiría tu código como este

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

Ahora, f2py no intentará reordenar la matriz al pasar.

Como nota al margen, esto también funcionará

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

porque detrás de escena, f2py realiza estas operaciones cuando pasa una matriz de orden C a una rutina de Fortran:

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

Pero, por supuesto, es más eficiente almacenar en orden Fortran desde el principio.

En general, no debería tener que preocuparse por el ordenamiento de arreglos a menos que tenga una sección de rendimiento crítico porque f2py se ocupa de esto por usted.

Answer 2

Parece que la respuesta es razonablemente simple:

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

funciona, pero aparentemente, esto es lo mismo que:

b=a.T

ya que transpose devuelve una vista y una vista rápida de b.flags comparación con a.flags muestra que esto es lo que quiero. ( b.flags es F_CONTIGUO).