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 BNuma 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 BCython 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 Bf2py 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 DOEND SUBROUTINE F2PY_runCython compilation:setup(ext_modules=cythonize(Extension('runtest',['runtest.pyx'])), include_dirs=[numpy.get_include()])f2py compilation:f2py -c runFortran.F90 --f90flags=-O3 -m fortranrunThe 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 1By 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 msecAs it can be seen almost all the time is used in the interface! Here is how to function is called:import numpy as npimport timeimport fortranrunA = np.random.rand(100,100,100,100)start = time.time()B = fortranrun.f2py_run(A)print(time.time()-start, B, 'f2py')
Loop optimization with f2py, Cython and Numba
python;fortran;performance
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 linuxType 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 msecHere 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=1See also Internal organization of numpy arrays