# Loop optimization with f2py, Cython and Numba

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

• Did you check the ordering of the indices for the f2py case? – nicoguaro Aug 19 '17 at 16:01
• Yes, I tried i,j,k,l and l,k,j,i. The one above is the fastests of the two cases – Erik Kjellgren Aug 19 '17 at 16:20
• I would think that l, k, j, i was faster. Might it be related to how the array is stored in numpy? docs.scipy.org/doc/numpy/f2py/python-usage.html#array-arguments – nicoguaro Aug 19 '17 at 16:35
• l,k,j,i is faster with f2py. A(l,k,j,i) is in the above f2py implementation – Erik Kjellgren Aug 19 '17 at 18:13
• I misread then, sorry about that. – nicoguaro Aug 19 '17 at 18:17

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