Loop optimization with f2py, Cython and Numba

Question

I tried to used f2py, Cython and Numba to make a simple loop structure be faster in python. Python implementation:

def Python(A):
    B = 0.0
    for i in range(100):
        for j in range(100):
            for k in range(100):
                for l in range(100):
                    B = B + A[i,j,k,l]
return B

Numa implementation:

@jit(double(double [:,:,:,:]), nopython=True)
def Numba(A):
    B = 0.0
    for i in range(100):
        for j in range(100):
            for k in range(100):
                for l in range(100):
                    B = B + A[i,j,k,l]
return B

Cython implementation:

cdef double Cython(double [:,:,:,:] A):
    cdef double B
    cdef int i, j, k, l
    B = 0.0
    for i in range(0,100):
        for j in range(0,100):
            for k in range(0,100):
                for l in range(0,100):
                    B = B + A[i,j,k,l]
return B

f2py implementation:

SUBROUTINE F2PY_run(A,B)
    INTEGER :: i, j, k, l
    REAL(8), DIMENSION(100,100,100,100), INTENT(in) :: A
    REAL(8), INTENT(out) :: B

    B = 0.0d0
    DO i=1, 100
        DO j=1, 100
            DO k=1, 100
                DO l=1,100
                   B = B + A(l,k,j,i)
                END DO
            END DO
        END DO
    END DO
END SUBROUTINE F2PY_run

Cython compilation:

setup(ext_modules=cythonize(Extension('runtest',['runtest.pyx'])), include_dirs=[numpy.get_include()])

f2py compilation:

f2py -c runFortran.F90 --f90flags=-O3 -m fortranrun

The functions was run by passing a np.random.rnd(100,100,100,100) into the functions. The timings were as follows:

22.998 sec; Python 
 0.114 sec; Cython 
 2.678 sec; f2py   
 0.106 sec; numba  

I know that for this specific problem I have set up, I can use np.einsum(). np.einsum() had a timing of 0.061 sec. But for my "real" code the problem cannot be solved with einsum. It can be seen that Cython and Numba executes at about the same speed, whereas f2py is much slower. The f2py function is written in Fortran90, so I would have thought it would have been atleast as fast. My question is thus, am I using f2py wrongly? Or does f2py have speed problems?

UPDATE 1 By suggestion from Stefano M, I made a F2PY perfermance report.

                  /-----------------------\ 
                 < F2PY performance report >  
                  \-----------------------/ 
Overall time spent in ...
(a) wrapped (Fortran/C) functions           :      105 msec 
(b) f2py interface,                1 calls  :     2925 msec
(c) call-back (Python) functions            :        0 msec
(d) f2py call-back interface,      0 calls  :        0 msec
(e) wrapped (Fortran/C) functions (acctual) :      105 msec

As it can be seen almost all the time is used in the interface! Here is how to function is called:

import numpy as np
import time
import fortranrun

A = np.random.rand(100,100,100,100)

start = time.time()
B = fortranrun.f2py_run(A)
print(time.time()-start, B, 'f2py')

Answer 1

I think that the problem is linked to the way in which f2py generates the fortran interface: the argument to fortranrun.f2py should be stored as a F_CONTIGUOUS array, otherwise the interface will create an internal copy with the correct storage order.

Python 3.6.2 (default, Jul 22 2017, 21:19:22) 
[GCC 7.1.1 20170516] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy as np
>>> import fortranrun
>>> A = np.random.rand(100,100,100,100)
>>> A.flags
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
>>> A.T.flags
  C_CONTIGUOUS : False
  F_CONTIGUOUS : True
  OWNDATA : False
  WRITEABLE : True
  ALIGNED : True
  UPDATEIFCOPY : False
>>> b = fortranrun.f2py(A.T)
>>> 
                      /-----------------------\
                     < F2PY performance report >
                      \-----------------------/
Overall time spent in ...
(a) wrapped (Fortran/C) functions           :      761 msec
(b) f2py interface,                1 calls  :        0 msec
(c) call-back (Python) functions            :        0 msec
(d) f2py call-back interface,      0 calls  :        0 msec
(e) wrapped (Fortran/C) functions (acctual) :      761 msec

Here you see that A is C_CONTIGUOUS, while the transpose A.T is F_CONTIGUOUS, and calling b = fortranrun.f2py(A.T) has no f2py interface overhead.

On the contrary, when you call b = fortranrun.f2py(A), since A is not F_CONTIGUOUS, the f2py interface will allocate scratch memory and copy the original array A in the correct storage order.

Memory-order is a common topic for C-FORTRAN interfaces: you should either construct A as a F_CONTIGUOUS array, or pass A.T to your fortran routine, but remember that you are operating on the transpose array.

Warning about array argument copies can be generated with compile option

-DF2PY_REPORT_ON_ARRAY_COPY=1

See also Internal organization of numpy arrays