Search code examples
pythonarraysnumpyfortranctypes

How to run a Fortran script with ctypes?


First of all I have no experience coding with fortran. I am trying to run a fortran code with python ctypes. I used the command gfortran -shared -g -o test.so test.f90 to convert my test.f90 file (code below) to test.so.

After reading C function called from Python via ctypes returns incorrect value and table in chapter "Fundamental data types" https://docs.python.org/3/library/ctypes.html#module-ctypes I had a clue about passing the right data type in my python code like ctypes.c_double(123) for real(kind=c_double). I received a Type Error: wrong type. I have no idea about passing the right data type to real(kind=c_double),dimension(*) it seems like an array for me. On Numpy array to ctypes with FORTRAN ordering and https://numpy.org/doc/stable/reference/generated/numpy.ndarray.ctypes.html it is well explained that numpy is providing a ctype conversion. So I tried some random arrays to pass values with np.array([[1,2],[3,4]]).ctypes.data_as(POINTER(c_double)). Now I receiving another type error: TypeError: item 2 in _argtypes_ has no from_param method.

I shared my python code also below. I could not figure out which data types I need to pass in my function. I would be grateful for an explanation.

function sticksum( anzd, w, b, a, weight, val, maxTchan, sthresh) result(spek) bind(c, name="sticksum")
    use, intrinsic :: iso_c_binding

    real(kind=c_double) :: anzd
    real(kind=c_double) :: b,a
    real(kind=c_double),dimension(*) :: w,weight
    real(kind=c_double),dimension(*) :: val
    integer(kind=c_int) :: maxTchan
    real(kind=c_double) :: sthresh
    real(kind=c_double) :: spek

    integer :: i, t, j,anz
    real    :: stick, maxw
    
    
    anz = anzd
    spek = 123
    i=1
    t = w(i) * b + a + 1.5
    if(t >= 1) THEN
        spek = anz
        stick = 0. + val(t)
        maxw = weight(i)*stick
        do i=2,anz
            t = w(i) * b + a + 1.5
            if(t > maxTchan) exit
            stick = val(t)
            maxw = max(weight(i)*stick,maxw)
            if( (w(i)*w(i)-w(i-1)*w(i-1)) > 0.5) THEN
                spek = spek + maxw
                maxw = 0
            end if
        end do
    end if
end function sticksum

from ctypes import *
import numpy as np 
so_file = "./test.so"
my_functions = CDLL(so_file)

print(type(my_functions))

my_functions.sticksum_4.argtypes = [c_double,np.ndarray.ctypes,c_double,c_double,np.ndarray.ctypes,np.ndarray.ctypes,c_int, c_double]
my_functions.restype = c_double

anzd = c_double(123) 
w = np.array([[1,2],[3,4]]).ctypes.data_as(POINTER(c_double))
b=c_double(123) 
a=c_double(123) 
weight=np.array([[1,2],[3,4]]).ctypes.data_as(POINTER(c_double))
val=np.array([[1,2],[3,4]]).ctypes.data_as(POINTER(c_double))
maxTchan=c_int(123)
sthresh=c_double(123)


sum = my_functions.sticksum_4(anzd,w,b,a,weight,val,maxTchan,sthresh)


Solution

  • I want to start by stating that Fortran is not among my skills - I came in touch with it ~5 times (other questions from SO).

    I found a bunch of problems:

    • Fortran:

      • Arguments are passed by reference to functions / subroutines ([GNU.GCC]: Argument passing conventions). That took me an hour to figure out. After I fixed all the rest, I was getting Access Violation when even printing one such argument. I added the dummy function during investigation. Anyway, the simplest version is to specify value for the ones to be passed by value (the other way around would be to pass pointers from the calling code)

      • Since I am on Win, I had to add the dllexport attribute to each function to be exported by the .dll (I was also able to do it by passing /EXPORT to the linker for each function). I had another setback, as when checking the function export, the tool that I primarily use (check [SO]: Discover missing module using command-line ("DLL load failed" error) (@CristiFati's answer)) doesn't show it:

        • Dependencies:

          Img0

        • Dependency Walker:

          Img1

      • 1st argument is the dimension, no need to be double

      • Other (minor) issues (like reorganizing, adding module stuff, ...)

    • Python:

      • Argument types didn't match the ones from Fortran definition which yields (as stated by the URL in the question) Undefined Behavior. I illustrated 2 ways of passing NumPy arrays (the uncommented one seems a bit nicer to me).
        Note that when working with multidimensional arrays, axes order is different in Fortran and C (for example for 2D, the former is column major)

    Managed to get the example working (I have no idea what all the operations from sticksum are supposed to do).

    • dll00.f90:

      module dll00
      
      use iso_c_binding
      implicit none
      
      contains
      
      #if 0
      function dummy(i0, d0) result(res) bind(C, name="dummy")
          !dec$ attributes dllexport :: dummy
          integer(kind=c_int), intent(in), value :: i0
          real(kind=c_double), intent(in), value :: d0
          integer(kind=c_int) :: res
      
          res = 1618
          print *, "xxx", res
          print *, "yyy", i0
          print *, "zzz", d0
      end function dummy
      #endif
      
      function sticksum(anz, w, b, a, weight, value, maxTchan, sthresh) result(spek) bind(C, name="sticksum")
          !dec$ attributes dllexport :: sticksum
          integer(kind=c_int), intent(in), value :: anz, maxTchan
          real(kind=c_double), intent(in), value :: b, a, sthresh
          real(kind=c_double), intent(in), dimension(*) :: w, weight, value
          real(kind=c_double) :: spek
      
          integer :: i, t, j
          real :: stick, maxw
      
          spek = 123
          i = 1
          t = w(i) * b + a + 1.5
          if (t >= 1) then
              spek = anz
              stick = 0. + value(t)
              maxw = weight(i) * stick
              do i = 2, anz
                  t = w(i) * b + a + 1.5
                  if (t > maxTchan) exit
                  stick = value(t)
                  maxw = max(weight(i) * stick, maxw)
                  if ((w(i) * w(i) - w(i - 1) * w(i - 1)) > 0.5) then
                      spek = spek + maxw
                      maxw = 0
                  end if
              end do
          end if
      end function sticksum
      
      end module dll00
      
    • code00.py:

      #!/usr/bin/env python
      
      import ctypes as cts
      import sys
      
      import numpy as np
      
      
      DLL_NAME = "./dll00.{:s}".format("dll" if sys.platform[:3].lower() == "win" else "so")
      
      ELEMENT_TYPE = cts.c_double
      
      DblPtr = cts.POINTER(cts.c_double)
      
      
      def np_mat_type(dim, element_type=ELEMENT_TYPE):
          return np.ctypeslib.ndpointer(dtype=element_type, shape=(dim,), flags="F_CONTIGUOUS")
      
      
      def main(*argv):
          dim = 5
          dll = cts.CDLL(DLL_NAME)
          sticksum = dll.sticksum
          sticksum.argtypes = (cts.c_int, np_mat_type(dim), cts.c_double, cts.c_double, np_mat_type(dim), np_mat_type(dim), cts.c_int, cts.c_double)
          #sticksum.argtypes = (cts.c_int, DblPtr, cts.c_double, cts.c_double, DblPtr, DblPtr, cts.c_int, cts.c_double)  # @TODO - cfati: alternative
          sticksum.restype = cts.c_double
      
          if 0:
              dummy = dll.dummy
              dummy.argtypes = (cts.c_int, cts.c_double)
              dummy.restype = cts.c_int
              print(dummy(cts.c_int(3141593), cts.c_double(2.718282)))
      
          w = np.arange(dim, dtype=ELEMENT_TYPE)
          a = 1
          b = 2
          weight = np.arange(1, dim + 1, dtype=ELEMENT_TYPE)
          val = np.arange(2, dim + 2, dtype=ELEMENT_TYPE)
          #print(w, weight, val)
      
          mtc = 3
          stresh = 6
      
          res = sticksum(dim, w, b, a, weight, val, mtc, stresh)
          #res = sticksum(dim, np.asfortranarray(w, dtype=ELEMENT_TYPE).ctypes.data_as(DblPtr), b, a, np.asfortranarray(weight, dtype=ELEMENT_TYPE).ctypes.data_as(DblPtr), np.asfortranarray(val, dtype=ELEMENT_TYPE).ctypes.data_as(DblPtr), mtc, stresh)  # @TODO - cfati: alternative
          print("\n{:s} returned: {:.3f}".format(sticksum.__name__, res))
      
      
      if __name__ == "__main__":
          print("Python {:s} {:03d}bit on {:s}\n".format(" ".join(elem.strip() for elem in sys.version.split("\n")),
                                                         64 if sys.maxsize > 0x100000000 else 32, sys.platform))
          rc = main(*sys.argv[1:])
          print("\nDone.\n")
          sys.exit(rc)
      

    Output:

    [cfati@CFATI-5510-0:e:\Work\Dev\StackExchange\StackOverflow\q075937942]> sopr.bat
    ### Set shorter prompt to better fit when pasted in StackOverflow (or other) pages ###
    
    [prompt]> "c:\Install\pc032\Microsoft\VisualStudioCommunity\2019\VC\Auxiliary\Build\vcvarsall.bat" x64 > nul
    
    [prompt]>
    [prompt]> dir /b
    code00.py
    dll00.f90
    
    [prompt]>
    [prompt]> "f:\Install\pc032\Intel\OneAPI\Version\compiler\2021.3.0\windows\bin\intel64\ifort.exe" /nologo /fpp /c dll00.f90
    
    [prompt]>
    [prompt]> link /NOLOGO /DLL /OUT:dll00.dll /LIBPATH:"f:\Install\pc032\Intel\OneAPI\Version\compiler\2021.3.0\windows\compiler\lib\intel64_win" dll00.obj
       Creating library dll00.lib and object dll00.exp
    
    [prompt]> dir /b
    code00.py
    dll00.dll
    dll00.exp
    dll00.f90
    dll00.lib
    dll00.mod
    dll00.obj
    
    [prompt]>
    [prompt]> "e:\Work\Dev\VEnvs\py_pc064_03.10_test0\Scripts\python.exe" ./code00.py
    Python 3.10.9 (tags/v3.10.9:1dd9be6, Dec  6 2022, 20:01:21) [MSC v.1934 64 bit (AMD64)] 064bit on win32
    
    
    sticksum returned: 5.000
    
    Done.
    

    Might also be interesting to read: