Search code examples
pythonfortranf2py

Unable to change and return value in Fortran subroutine + notebook


I have a simple fortran subroutine (just for testing python-fortran interface). It looks like this:

subroutine sum2(x,y,z)
 real(kind=8),intent(in)::x,y
 real(kind=8),intent(inout)::z
 z = x + y
 print *, "sum is ", z
end subroutine sum2

After compilation with f2py I go to python and do the following:

>>> import sum2
>>> x = 1.0
>>> y = 2.0
>>> z = 0.0
>>> sum2.sum2(x,y,z)
sum is 3
>>> z
0.0

So, even though z is specified as inout, its value is not changed by the function. I need to know why. Another question concerns notebook. If I import sum2 there and run sum2.sum2(x,y,z) I do not even see a message sum is .... So, the question is, if it is possible to call fortran subroutines in a notebook?

EDIT

There was a reasonable comment, that in my example I used an immutable data type. So, I decided to change it, but still I have the same issue. So, my new fortran subroutine looks like this:

subroutine arr(x)
    real(kind=8),dimension(1)::x
    x(1) = 2 ! new value
    print *, "x[0] = ", x(1)
end subroutine arr

Again, I compile it with f2py and go to python:

>>> import arr
>>> x = [1]
>>> arr.arr(x)
x[0] = 2.000000000
>>> x
[1]

So, even though I now use a mutable type of data, I still have the same problem - I can't pass around variables between python and fortran code (or it is better to say, that I have one way road).


Solution

  • In order to interface correctly between between Fortran and Python, the C wrapper around Fortran needs numpy arrays.

    According to the notes inside the getting started tutorial of f2py (http://docs.scipy.org/doc/numpy-dev/f2py/getting-started.html#the-quick-way), you could best pass a numpy array with a compatible data type. For real(kind=8), this is np.float64.

    For the sum2 subroutine the following code worked for me:

    >>> import sum2
    >>> import numpy as np
    >>> x=np.array(1.0,dtype=np.float64)
    >>> y=np.array(2.0,dtype=np.float64)
    >>> z=np.array(0.0,dtype=np.float64)
    >>> sum2.sum2(x,y,z)
     sum is    3.0000000000000000
    >>> z
    array(3.0)
    >>>