CODEWITHSUNDEEP
pythonnumpyfortranf2py
atbug
atbug

Reputation: 838

f2py: input not fortran contiguous

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?

Upvotes: 3

Views: 3140

Answers (3)

Jason
Jason

Reputation: 3266

I got this fortran contiguous error too and made it gone by changing some intent(intout)s to intent(in) in the .pyf signature file, e.g.:

real(kind=8) dimension(:,:),intent(in) :: kernel
integer dimension(:,:),intent(in) :: kernelmask

kernel and kernelmask are both intent(in) according to the code, but f2py generated .pyf gives them intent(inout). Correcting this somehow solved the contiguous error. It may not apply to other cases, but considering so little info out there you might want to give it a try should similar issue happen.

Upvotes: 0

g.b.
g.b.

Reputation: 347

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

Upvotes: 9

farenorth
farenorth

Reputation: 10771

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.

Upvotes: 2

Related Questions