CODEWITHSUNDEEP
pythonnumpyctypes
BenT
BenT

Reputation: 3200

Issues Transferring Numpy array into Fortran using Ctypes

I am trying to wrap some fortran code into python using the ctypes library but have am having major issues with the data transfer. When I print the data from my python script, it looks substantially different from when I print it within the fortran code. Can someone please help me figure out what is going on here? I tried playing around with the datatypes which did not fix the solution and so far I have not found any other SO questions that have addressed my issue. Also please note that using f2py will not work for the end product that this example refers too.

Below is an example of my code:

temp.f90

       subroutine test(num_mod_sel,goodval,x,y,coeff,coeff_flag)
       implicit none


! Input/Output
       integer num_mod_sel, goodval
       real, dimension(goodval,num_mod_sel) :: x
       real, dimension(goodval)             :: y
       real, dimension(num_mod_sel)         :: coeff
       integer, dimension(num_mod_sel)      :: coeff_flag


      print*, num_mod_sel,goodval,x


      return
      end subroutine test

!===================================================================================================!

The above f90 code is compiled with:

gfortran -shared -fPIC -o temp.so temp.f90

test.py

from ctypes import CDLL,cdll, POINTER, c_int, c_float
import numpy as np


def test_fcode(num_mod_sel,goodval,x,y,coeff,coeff_flag):

    fortran = CDLL('./temp.so')

    fortran.test_.argtypes = [ POINTER(c_int),
                                     POINTER(c_int),
                                     POINTER(c_float),
                                     POINTER(c_float),
                                     POINTER(c_float),
                                     POINTER(c_int) ]

    fortran.test_.restype = None

    num_mod_sel_ = c_int(num_mod_sel)
    goodval_     = c_int(goodval)
    x_           = x.ctypes.data_as(POINTER(c_float))
    y_           = y.ctypes.data_as(POINTER(c_float))
    coeff_       = coeff.ctypes.data_as(POINTER(c_float))
    coeff_flag_  = coeff_flag.ctypes.data_as(POINTER(c_int))

    fortran.test_(num_mod_sel_,goodval_,x_,y_,coeff_,coeff_flag_)


    
 #Create some test data   
    
num_mod_sel = 4
goodval = 10
x = np.full((num_mod_sel,goodval),999.,dtype=float)
x[:] = np.random.rand(num_mod_sel,goodval)
y = np.full(goodval,999.,dtype=float)
y[:] = np.random.rand(goodval)

coeff = np.empty(num_mod_sel,dtype=float)
coeff_flag = np.empty(num_mod_sel,dtype=int)


#Run the fortran code
test_fcode(num_mod_sel,goodval,x,y,coeff,coeff_flag)

print(x) from the python code:

[[0.36677304 0.8734628 0.72076823 0.20234787 0.91754331 0.26591916 0.46325577 0.00334941 0.98890871 0.3284262 ] [0.15428096 0.24979671 0.97374747 0.83996786 0.59849493 0.55188578 0.9668523 0.98441142 0.50954678 0.22003844] [0.54362548 0.42636074 0.65118397 0.69455346 0.30531619 0.88668116 0.97278714 0.29046492 0.64851937 0.64885967] [0.31798739 0.37279389 0.88855305 0.38754276 0.94985151 0.56566525 0.99488508 0.13812829 0.0940132 0.07921261]]

print*, x from f90:

-2.91465824E-17 1.68338645 13.0443134 1.84336567 -7.44153724E-34 1.80519199 -2.87629426E+27 1.57734776 -1297264.38 1.85438573 -236487.531 1.63295949 -1.66118658E-33 1.73162782 -6.73423983E-09 0.919681191 -1.09687280E+21 1.87222707 5.50313165E+09 1.66421306 8.38275158E+34 1.52928090 -2.15154066E-13 1.62479663 3.88800366E+30 1.86843681 127759.977 1.83499193 -3.55062879E+15 1.77462363 2.43241945E+19 1.76297140 3.16150975E-03 1.86671305 1.35183692E+21 1.87110281 1.74403865E-31 1.75238669 9.85857248E-02 1.59503841 -2.33541620E+30 1.79045486 -1.86185171E+11 1.78229403 4.23132255E-20 1.81525886 2.96771497E-04 1.82888138 -4.55096013E-26 1.86097753 0.00000000 3.68934881E+19 -7.37626273E+15 1.58494916E+29 0 -1064355840 -646470284 -536868869

Upvotes: 1

Views: 214

Answers (1)

Michael Butscher
Michael Butscher

Reputation: 10959

The problem is a mismatch of datatypes.

The Fortran real is usually a 32 bit float (C float) while numpy interprets the Python datatype float as numpy.float_ which is an alias of numpy.float64, the C double with 64 bits.

Solution: In Python use numpy.float32 as dtype for numpy array creation.

Upvotes: 2

Related Questions