f2py, problems passing a Python function to Fortran

Question

I have a this simple Fortran code (stack.f90):

      subroutine fortran_sum(f,xs,nf,nxs)
      integer nf,nxs
      double precision xs,result
      dimension xs(nxs),result(nf)
      external f 
      result = 0.0
      do I = 1,nxs
         result = result + f(xs(I))
         print *,xs(I),f(xs(I))
      enddo
      return
      end 

Which I am compiling using:

f2py -c --compiler=mingw32 -m stack2 stack2.f90

And then testing using this Python script (stack.py):

import numpy as np
from numpy import cos, sin , exp
from stack import fortran_sum
def func(x):
    return x**2

if __name__ == '__main__':
    xs = np.linspace(0.,10.,10)
    ans =  fortran_sum(func,xs,1)
    print 'Fortran:',ans
    print 'Python:',func(xs).sum()

When I run using "python stack.py" it gives:

   0.0000000000000000        0.00000000
   1.1111111111111112              Infinity
   2.2222222222222223              Infinity
   3.3333333333333335        9.19089638E-26
   4.4444444444444446              Infinity
   5.5555555555555554        0.00000000
   6.6666666666666670        9.19089638E-26
   7.7777777777777786        1.60398502E+09
   8.8888888888888893              Infinity
   10.000000000000000        0.00000000
Fortran: None
Python: 351.851851852

My questions are:

• why the function is not being evaluated correctly?

• how to return result to Python?

• is it possible to evaluate the array xs at once in Fortran?

Thank you!

---

EDIT: With the great tips from @SethMMorton I came ended up with the following:

      subroutine fortran_sum(f,xs,nxs,result)
      implicit none
      integer :: I
      integer, intent(in) :: nxs
      double precision, intent(in) :: xs(nxs)
      double precision, intent(out) :: result
      double precision :: f
      external :: f 
      ! "FIX" will be added here
      result = 0.0
      do I = 1,nxs
        result = result + f(xs(I))
        print *,xs(I),f(xs(I))
      enddo
      return
      end 

Running stack.py with this command modified: ans = fortran_sum(func,xs); gives:

   0.0000000000000000        0.0000000000000000
   1.1111111111111112        3.8883934247189009E+060
   2.2222222222222223        3.8883934247189009E+060
   3.3333333333333335        9.1908962428537221E-026
   4.4444444444444446        3.8883934247189009E+060
   5.5555555555555554        5.1935286092977251E-060
   6.6666666666666670        9.1908962428537221E-026
   7.7777777777777786        1603984978.1728516
   8.8888888888888893        3.8883934247189009E+060
   10.000000000000000        0.0000000000000000
Fortran: 1.55535736989e+61
Python: 351.851851852

Which is wrong. This weird behavior does not happen if I add an intermediate variable x=x(I) AND call the function using this variable f(x). The funny thing is that if I call f(x) once, the desired call f(x(I)) also works. After applying this "FIX":

double precision :: x, f_dummy
x = xs(I)
f_dummy = f(x)

then compiling and running, the right result is obtained:

   0.0000000000000000        0.0000000000000000
   1.1111111111111112        1.2345679012345681
   2.2222222222222223        4.9382716049382722
   3.3333333333333335        11.111111111111112
   4.4444444444444446        19.753086419753089
   5.5555555555555554        30.864197530864196
   6.6666666666666670        44.444444444444450
   7.7777777777777786        60.493827160493836
   8.8888888888888893        79.012345679012356
   10.000000000000000        100.00000000000000
Fortran: 351.851851852
Python: 351.851851852

It would be nice if someone could explain why?

Answer 1

Why is the function not being evaluated correctly?

You are sending the result of the callback function f directly to the print method. Because Fortran has no idea what data type f will returned (since it is not declared and it is not a Fortran function), print cannot correctly interpret the data to print correctly. if you assign the result to a variable, then print the variable you should see the expected result:

double precision x
....
x = f(xs(1))
print *, x

How to return result to Python?

You need to inform f2py of the intent of the variables. This is actually possible to do with standard Fortran 90 syntax, which I HIGHLY recommend doing since it makes your code clearer (I know you are using Fortran 90 because you are able to print an entire array without looping over it.).

subroutine stack2(f,xs,result,nf,nxs)
    implicit none
    integer,          intent(in)  :: nf, nxs
    double precision, intent(in)  :: xs(nxs)
    double precision, intent(out) :: result(nf)
    external f
    ...
end subroutine stack2

Now, it is clear which variables come into the routine, and which are going out. Note that result has to be an argument to the subroutine in order for it to be passed out.

f2py will now understand that result should be returned from the stack2 function.

Is it possible to evaluate the array xs at once in Fortran

I'm not sure what you mean exactly, but I assume you want to know if you can do this on the entire array at a time, instead of in a do loop.

Possibly

Since xs is an array, you should be able to just pass the whole array to f and it should operate on the entire array. This might not work, because Fortran requires a function to be ELEMENTAL to do this, and this might not be within the scope of f2py. If f2py does not make the function ELEMENTAL, then you must do this in a loop.

An issue you might want to address

The line result = result + f(xs(I)) is assigning the same value to every element of result. It is not clear if this is what you want, but it might be something to address if it is not.

A code-based summary

Below is an example of a version of your code using these new suggestions.

subroutine stack2(f,xs,result,nf,nxs)
    implicit none
    integer,          intent(in)  :: nf, nxs
    double precision, intent(in)  :: xs(nxs)
    double precision, intent(out) :: result(nf)
    double precision :: x
    integer          :: I
    external f
    result = 0.0d0 ! Make this a double constant since result is double
    do I = 1,nxs
        x = f(xs(I))
        result = result + x
        print *, x
    end do
    print *, xs
    return ! Not needed in Fortran 90
end subroutine stack2

Note that f2py will properly figure out the length of your input array, so when you call it you only need to define the length of the result:

import numpy as np
from stack2 import stack2
def func(x):
    return x**2

if __name__ == '__main__':
    xs = np.linspace(0.,10.,10)
    ans =  stack2(func,xs,5) # This was changed
    print 'ans:',ans