f2py: input not fortran contiguous

Question

I wrapped some fortran code with f2py. Here's the fortran code:

MODULE iteration
implicit none
contains

SUBROUTINE iterate(alpha, beta, e, es, rank, omega, smearing, prec, max_step)
    REAL(kind=8), INTENT(in) :: omega, smearing, prec
    INTEGER :: max_step, step, rank, cnt
    COMPLEX(kind=16) :: alpha(rank,rank), beta(rank,rank), omega_mat(rank, rank), green(rank, rank)
    COMPLEX(kind=16), INTENT(inout) ::  e(rank,rank), es(rank,rank)
    step = 0
    omega_mat = 0
    DO cnt=1, rank
        omega_mat(cnt, cnt) = 1.0
    ENDDO
    omega_mat = omega_mat * (omega + (0.0, 1.0) * smearing)
    DO WHILE (minval(abs(alpha)) .gt. prec .or.  minval(abs(beta)) .gt. prec .and. step .lt. max_step)
        green = zInverse(rank, omega_mat - e)
        e = e + matmul(alpha, matmul(green, beta)) + matmul(beta, matmul(green, alpha))
        es = es + matmul(alpha, matmul(green, beta))
        alpha = matmul(alpha, matmul(green, alpha))
        beta = matmul(beta, matmul(green, beta))
        step = step + 1
    ENDDO
END SUBROUTINE iterate

FUNCTION zInverse(n, a)  result(ra)
    INTEGER :: n,lda,ipiv(n),info,lwork
    COMPLEX(kind=16)::a(n,n),ra(n,n),work(n)
    ra=a
    lwork=n
    lda=n
    CALL zgetrf(n, n, ra, lda, ipiv, info)
    IF(info/=0) WRITE(0,*) 'Error occured in zgetrf!'
    CALL zgetri(n, ra, lda, ipiv, work, lwork, info)
    IF(info/=0) WRITE(0,*) 'Error occured in zgetri!'
END FUNCTION zInverse
END MODULE iteration

and then I compiled the code with f2py -L/usr/lib -llapack -m pyiteration -c iteration.F90, and tested with

import numpy as np
import pyiteration
alpha = np.array([[1,0],[0,1]], dtype='complex')
beta = np.array([[1,0],[0,1]], dtype='complex')
e = np.array([[1,0],[0,1]], dtype='complex')
es = np.array([[1,0],[0,1]], dtype='complex')
# f2py is automatically generating rank for me
pyiteration.iteration.iterate(alpha,beta, e, es, 1.0, 0.001, 0.001, 100)

However, I got the following error: ValueError: failed to initialize intent(inout) array -- input not fortran contiguous.

I googled and found f2py should automatically make the array fortran contiguous. Then what is happening here?

Answer 1

You have to create the arrays as Fortran-contiguous ordered using order='F', so that:

alpha = np.array([[1,0],[0,1]], dtype='complex', order='F')
beta = np.array([[1,0],[0,1]], dtype='complex', order='F')
e = np.array([[1,0],[0,1]], dtype='complex', order='F')
es = np.array([[1,0],[0,1]], dtype='complex', order='F')

Answer 2

Arrays that are created within Fortran (i.e. out) will return as Fortran contiguous. However, my understanding is that all arrays that are passed in (i.e. in and inout) must be specified to be Fortran contiguous in Python. If f2py were changing them from C to F contiguous, that would probably mean creating a copy, which would take extra memory, which isn't very efficient.

To fix this, all you should need to do is add a order='F' kwarg to each of your np.array calls.