3

I am trying to understand how F2PY works. To do so, I wrote a simple Fortran function which takes an array as input and returns the sum of the elements of the array.

I wrote three different versions of the same functions, which I expect to hold the same result:

function rsum(arr)
real, dimension(:), intent(in) :: arr
real :: rsum
rsum=sum(arr)
end function rsum

function rsum1(arr)
real(8), dimension(:), intent(in) :: arr
real(8) :: rsum1
rsum1=sum(arr)
end function rsum1

function rsum2(arr) result(s)
real, dimension(:), intent(in) :: arr
real :: s
s=sum(arr)
end function rsum2

function rsum3(arr) result(s)
real(8), dimension(:), intent(in) :: arr
real(8) :: s
s=sum(arr)
end function rsum3

My python script to test these functions is the following:

from numpy import *
import ftest as f

a=array(range(3))

print(f.rsum(a))
print(f.rsum1(a))
print(f.rsum2(a))
print(f.rsum3(a))

but the result is this:

3.0
0.0
3.0
3.0

All the results are correct, except the one of rsum1, which is 0.0. What I find even more strange is that rsum3, in which I merely change the name of the result of the function (or, at least, I think I am doing so), works perfectly!

I know this has something to do with the type conversion between Fortran and numpy, but I don't understand what the problem is.

PS: I only learned Fortran very recently.

Improve this question
6
  • 2
    If you are new, don't learn the bad habit of real(8), where do people even pick it up? It is not in any good textbook. Commented Mar 14, 2017 at 17:08
  • 2
    rsum1 and rsum3 should be completely identical. But maybe f2py interprets one of the wrong. Commented Mar 14, 2017 at 17:09
  • I can replicate it. When I call them in a Fortran subroutine which is called from f2py the result is correct. But even in Python the reported signature of the f.* functions is correct. Commented Mar 14, 2017 at 17:32
  • @VladimirF Regarding the notation real(8), someone taught it to me. Why is it so bad? Also, I don't understand: are you getting my same error when running the script? Commented Mar 14, 2017 at 23:02
  • Yes, I am getting the same error. Real(8) is not portable, it will not compile with some compilers. See stackoverflow.com/documentation/fortran/939/data-types/4390/… and stackoverflow.com/a/856243/721644 Commented Mar 14, 2017 at 23:21

1 Answer 1

Reset to default
1

Short answer and workaround

The root cause of the problem is related the use of assumed-shape dummy arguments (i.e. arr) in your functions. Fortran requires such functions to have explicit interfaces. @VladimirF gave an excellent answer to your (related?) question here about just that, indicating that the preferred solution is to put the functions in a module. Assuming your code listing with functions is saved in a file called funcs.f90, you can simply put them into a module, e.g. called mod_funcs.f90, like so:

module mod_funcs
    implicit none
    contains
        include "funcs.f90"
end module mod_funcs

Wrap this with F2PY python -m numpy.f2py -m ftest -c mod_funcs.f90, update the import statement in your python test script to from ftest import mod_funcs as f and then run it to get the expected result:

3.0
3.0
3.0
3.0

The longer answer and explanation

Fortran functions are wrapped in subroutines by F2PY. In order to support assumed-shape arrays in a Fortran standard compliant way, the F2PY created subroutine wrappers contain interfaces for user defined functions with assumed-shape dummy arguments. You can have a look at these wrappers by specifying a build directory with a --build-dir flag when wrapping with F2PY, e.g. like this:

python -m numpy.f2py --build-dir .\build -m ftest -c funcs.f90

Looking at the wrapper that is created for the problematic function rsum1 is revealing (I'm copying verbatim from ftest-f2pywrappers.f keeping F2PY's indentation):

subroutine f2pywraprsum1 (rsum1f2pywrap, arr, f2py_arr_d0)
integer f2py_arr_d0
real(8) arr(f2py_arr_d0)
real(8) rsum1f2pywrap
interface
function rsum1(arr) 
    real(8), dimension(:),intent(in) :: arr
end function rsum1
end interface
rsum1f2pywrap = rsum1(arr)
end

Note that due to implicit data typing rules, the interface for rsum1 implies a function with real data type, not real(8) as intended - so there is a data type mismatch in the interface! This explains why the seemingly identical function with an explicit result statement (rsum3) returns the correct result in your original example, its result has the correct data type. By good fortune in naming your functions, rsum has the correct interface. If you change rsum's name to e.g. isum, implicit data typing rules in its F2PY subroutine wrapper interface will imply that it has an integer result, and you will get the following output from your (modified to reflect the name change from fsum to isum) python script:

0.0
0.0
3.0
3.0

So to me it appears as if there may be a bug in how F2PY creates interfaces for functions with assumed-shape dummy arguments (which can be bypassed by putting those functions into a module directly, or by explicitly declaring return value of the function using result).

For completeness, I was using Python 3.6.3 :: Intel Corporation, with NumPy 1.14.3, and GNU Fortran (GCC) 8.2.0.

Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.