I just discovered here recently, very glad that there is a community to discuss about Fortran.

Most of time I use python for visualization or some simple calculation. Of course I sometimes encounter performance problems and I just wonder if I can call Fortran subroutines or functions in python to solve this problem.

I have some experience in using `f2py`

or `ctypes`

to call Fortran subroutines/functions in python. But I’m not sure which one is faster.

I’ve written some codes and compared them, I found that the performances are similar but `f2py`

is slightly better. But this is only one case, I’m not sure that if using `f2py`

is always better than using `ctypes`

.

Does anyone has some experience for it? Or maybe knowing the other tool which is better than `f2py`

and `ctypes`

.

I decided to use `f2py`

and `ctypes`

just because I didn’t need to install anything else.

p.s.

If needed, here is the code what I compare between `f2py`

and `ctypes`

. They calculate the vertical component of vorticity.

BTW, pure `numpy`

solution is the fastest in my test… I am not sure if the cost in calling external library slowdown the performance of calling fortran subroutine?

p.s.2

I’m not sure if this topic is suitable in here. If not I’ll delete this post as soon as possible.

```
! f2py -c solve_f2py.f90 -m solvef90f2py
subroutine vorticity_f90loop(nx, ny, u, v, dx, dy, vort)
implicit none
integer, intent(in) :: nx, ny
real*8, intent(in) :: u(ny,nx), v(ny,nx)
real*8, intent(in) :: dx, dy
real*8, intent(out) :: vort(ny,nx)
real*8 :: dx2, dy2
integer :: i, j
! ---------------------------------------------
dx2 = 2. * dx
dy2 = 2. * dy
vort(:,:) = 0.
do i = 2, nx-1
do j = 2, ny-1
vort(j,i) = (v(j,i+1)-v(j,i-1)) / dx2 - (u(j+1,i)-u(j-1,i)) / dy2
end do
end do
return
end subroutine vorticity_f90loop
```

```
! gfortran -shared -fPIC solve_ctypes.f90 -o solvef90ctypes.so
subroutine vorticity_f90loop(nx, ny, u, v, dx, dy, vort) bind(c, name='vorticity_f90loop')
use iso_c_binding, only : c_double, c_int
implicit none
integer(c_int), intent(in) :: nx, ny
real(c_double), intent(in) :: u(ny,nx), v(ny,nx)
real(c_double), intent(in) :: dx, dy
real(c_double) :: vort(ny,nx)
real(c_double) :: dx2, dy2
integer(c_int) :: i, j
! ---------------------------------------------
dx2 = 2. * dx
dy2 = 2. * dy
vort(:,:) = 0.
do i = 2, nx-1
do j = 2, ny-1
vort(j,i) = (v(j,i+1)-v(j,i-1)) / dx2 - (u(j+1,i)-u(j-1,i)) / dy2
end do
end do
return
end subroutine vorticity_f90loop
```

```
# solve.py
import time
import ctypes as ct
import numpy as np
import solvef90f2py as f90f2py
def vorticity_py(u, v, dx, dy):
ny = len(u)
nx = len(u[0])
vort =[[0 for _ in range(nx)] for _ in range(ny)]
dx2 = dx * 2
dy2 = dy * 2
for j in range(1, ny-1):
for i in range(1, nx-1):
vort[j][i] = (v[j][i+1] - v[j][i-1]) / dx2 - (u[j+1][i] - u[j-1][i]) / dy2
return vort
def vorticity_np(u, v, dx, dy):
ny, nx = u.shape
vort = np.zeros((ny, nx))
vort[1:-1,1:-1] = (v[1:-1,2:] - v[1:-1,:-2]) / (2*dx) - (u[2:,1:-1] - u[:-2,1:-1]) / (2*dy)
return vort
def solve_f90ctypes_loop(nx, ny, u, v, dx, dy):
f90lib = ct.CDLL('solvef90ctypes.so')
f90lib.vorticity_f90loop.argtypes = [
ct.POINTER(ct.c_int),
ct.POINTER(ct.c_int),
np.ctypeslib.ndpointer(dtype=np.float64, ndim=2, shape=(ny, nx)),
np.ctypeslib.ndpointer(dtype=np.float64, ndim=2, shape=(ny, nx)),
ct.POINTER(ct.c_double),
ct.POINTER(ct.c_double),
np.ctypeslib.ndpointer(dtype=np.float64, ndim=2, shape=(ny, nx))
]
f90lib.vorticity_f90loop.restype = None
u = np.asfortranarray(u)
v = np.asfortranarray(v)
vort = np.empty_like(u, dtype=np.float64, order='F')
f90lib.vorticity_f90loop(
ct.byref(ct.c_int(nx)),
ct.byref(ct.c_int(ny)),
u,
v,
ct.byref(ct.c_double(dx)),
ct.byref(ct.c_double(dy)),
vort
)
return vort
if __name__ == '__main__':
# prepare data
dx = 0.5
dy = 0.5
ny = 1000
nx = 1000
u = np.random.rand(ny, nx)
v = np.random.rand(ny, nx)
u_lst = u.tolist()
v_lst = v.tolist()
# pure python
t1 = time.time()
vort_py = vorticity_py(u_lst, v_lst, dx, dy)
t2 = time.time()
print(f'(python) time = {t2-t1:.8f}')
# numpy
t1 = time.time()
vort_np = vorticity_np(u, v, dx, dy)
t2 = time.time()
print(f'(numpy) time = {t2-t1:.8f}')
# f2py - loop
t1 = time.time()
vort_f90loop = f90f2py.vorticity_f90loop(u, v, dx, dy)
t2 = time.time()
print(f'(f90+f2py, loop) time = {t2-t1:.8f}')
# ctypes - loop
t1 = time.time()
vort_ctloop = solve_f90ctypes_loop(nx, ny, u, v, dx, dy)
t2 = time.time()
print(f'(f90+ct, loop) time = {t2-t1:.8f}')
# check result
print()
print(vort_py[2][:5])
print(vort_np[2,:5])
print(vort_f90loop[2,:5])
print(vort_ctloop[2,:5])
```

Execute result:

```
$ python solve.py
(python) time = 0.32275057
(numpy) time = 0.01257682
(f90+f2py, loop) time = 0.01558447
(f90+ct, loop) time = 0.01573086
[0, -0.5746908151001381, -0.4502448148151549, 0.4887173001423759, 0.17991335715166268]
[ 0. -0.57469082 -0.45024481 0.4887173 0.17991336]
[ 0. -0.57469082 -0.45024481 0.4887173 0.17991336]
[ 0. -0.57469082 -0.45024481 0.4887173 0.17991336]
```