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?
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