How to run a Fortran script with ctypes?

Question

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)

Answer 1

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:

      

    • Dependency Walker:

      

  • 1^st 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:

• [SO]: What is numpy.ctypeslib.as_ctypes exacty doing (@CristiFati's answer)

• [SO]: How to convert a ctypes C pointer to a Numpy array? (@CristiFati's answer)

• [SO]: How can i cast a double pointer ctype to numpy array? (@CristiFati's answer)