F2py Linking Quadmath Libraries? Use Ctypes For Fortran Wrapper Instead?

Update 11/23/2019: This started out as a question about why I could not get f2py to work for a simple fortran wrapper. My 'answer' (below) is to use ctypes instead. Original post:

Solution 1:

I did not exactly figure out the problem with f2py, so please don't ding me for this nonanswer. All I wanted was to put a simple wrapper around some fortran. I solved that problem though, so this post is my answer. I am posting it in hopes it may help someone else with a similar problem.

After the post, I decided I needed to drop down to a lower level to figure out what was going on. Since f2py is built on the C-api, I decided to work with it directly. C - fortran interoperability is well defined. I got some simple code to work, but I was scared off when I started looking at passing arrays and read about reference counting. So, I decided the C-api it was too low level for my needs. Cython appears oriented toward writing Python-like code that can be compiled. I finally settled on ctypes, since all I really need is to convert the internal python types and then call (what it thinks is) c code. There are examples around but none seemed to fit exactly, so I will post some test code in hopes it may help others. The example shows passing ints, doubles and arrays back and forth by both functions and subroutines. The same code works in windows if .so is replaced with .dll.This information was helpful - How to dereference a memory location from python ctypes?

The fortran is:

Integer Function intfunc(n)
   integer, intent(in) :: n
   intfunc = n*n
end Function intfunc
! ---------------------------------------------------------------------------------------
subroutine Asub(acf,af)
   include 'defs.fi' !  integer, parameter :: float = selected_real_kind(12)
   Real(float), intent(in) :: acf
   Real(float), intent(out) :: af
   if(acf < 0.49)then
      af = 0.37464 + acf*(1.54226 - 0.26992*acf)
   else
      af = 0.379642 + acf*(1.48503 + acf*(-0.164423 + 0.016666*acf))
   endif
end subroutine Asub
! ---------------------------------------------------------------------------------------
Function Afactor(acf) result(af)
   include 'defs.fi'
   Real(float), intent(in) :: acf
   Real(float) :: af
   if(acf < 0.49)then
      af = 0.37464 + acf*(1.54226 - 0.26992*acf)
   else
      af = 0.379642 + acf*(1.48503 + acf*(-0.164423 + 0.016666*acf))
   endif
End Function Afactor
! ---------------------------------------------------------------------------------------
Function MultXY(x,y) result(z)
   include 'defs.fi'
   Real(float), intent(in) :: x,y
   Real(float)::z
   z = x*y
End Function MultXY
! ---------------------------------------------------------------------------------------
subroutine arrays(n,m,a)
   include 'defs.fi'
   integer, intent(in) :: n,m
   real(float), intent(out) :: a(n,m)
   integer :: i,j
   do j=1,m
      do i=1,n
         a(i,j) = real(100*i+j,float)
      enddo
   enddo
end subroutine arrays

The code is compiled into a shared library by:

gfortran -shared -o code.so code.f90 -fPIC -Xlinker -Map=code.map

It can be run with the following python code. A number of the steps in this code could be combined, but I thought it was more illustrative to break them down. I don't completely understand the function prototypes (argtypes & restypes). They seem to allow some shortcuts, such that Python performs the appropriate cast. Can someone elaborate?

import ctypes as ct
import numpy as np
from numpy.ctypeslib import ndpointer
flib = ct.CDLL('code.so')
flib.asub_.argtypes = [ct.POINTER(ct.c_double),ct.POINTER(ct.c_double)]
flib.afactor_.argtypes = [ct.POINTER(ct.c_double)]
flib.afactor_.restype = ct.c_double
flib.multxy_.argtypes = [ct.POINTER(ct.c_double),ct.POINTER(ct.c_double)]
flib.multxy_.restype = ct.c_double
flib.arrays_.argtypes = [ct.POINTER(ct.c_int),ct.POINTER(ct.c_int),ndpointer(ct.c_double)]
def main():
   nval = 4
   nptr = ct.pointer(ct.c_int(nval))
   # integer function intfunc(n)
   nx = flib.intfunc_(nptr)
   print('intfunc(',nval,') =',nx)
   x = 0.35
   xc = ct.c_double(x)
   yc = ct.c_double(0.50)
   zc = ct.c_double(0.0)
   # subroutine Asub(acf,af)
   flib.asub_(ct.byref(xc),ct.byref(zc))
   z = zc.value
   x = xc.value
   y = yc.value
   print(' x,z =',x,z)
   # Function Afactor(acf) result(af)
   z = flib.afactor_(ct.byref(xc))
   print(' x,z =',x,z)
   # Function MultXY(x,y) result(z)
   z = flib.multxy_(ct.byref(xc),ct.byref(yc))
   print('multxy(',x,y,') =',z)
   n = 3
   m = 4
   nc = ct.c_int(n)
   mc = ct.c_int(m)
   a = np.empty((m,n),dtype=np.double) # m,n not n,m
   # subroutine arrays(n,m,a)
   flib.arrays_(ct.byref(nc),ct.byref(mc),a)
   for i in range(n):   # note reversal of indices
      for j in range(m):
         print('i,j,a(j,i) =',i+1,j+1,a[j,i])
main()

The output is:

intfunc( 4 ) = 16
 x,z = 0.35 0.8813658067584037
 x,z = 0.35 0.8813658067584037
multxy( 0.35 0.5 ) = 0.175
i,j,a(j,i) = 1 1 101.0
i,j,a(j,i) = 1 2 102.0
i,j,a(j,i) = 1 3 103.0
i,j,a(j,i) = 1 4 104.0
i,j,a(j,i) = 2 1 201.0
i,j,a(j,i) = 2 2 202.0
i,j,a(j,i) = 2 3 203.0
i,j,a(j,i) = 2 4 204.0
i,j,a(j,i) = 3 1 301.0
i,j,a(j,i) = 3 2 302.0
i,j,a(j,i) = 3 3 303.0
i,j,a(j,i) = 3 4 304.0