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.

francescalus
  • 30,576
  • 16
  • 61
  • 96
valerio
  • 677
  • 4
  • 12
  • 25
  • 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. – Vladimir F Героям слава Mar 14 '17 at 17:08
  • 2
    `rsum1` and `rsum3` should be completely identical. But maybe f2py interprets one of the wrong. – Vladimir F Героям слава Mar 14 '17 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. – Vladimir F Героям слава Mar 14 '17 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? – valerio Mar 14 '17 at 23:02
  • Yes, I am getting the same error. Real(8) is not portable, it will not compile with some compilers. See http://stackoverflow.com/documentation/fortran/939/data-types/4390/precision-of-floating-point-numbers#t=201703142319329773677 and http://stackoverflow.com/a/856243/721644 – Vladimir F Героям слава Mar 14 '17 at 23:21
  • I have experimented with f2py a few times and as I recall f2py doesn't play great with assumed shape arrays. Maybe pass the size of the array to the function and check. – Vikram Mar 15 '17 at 11:29

1 Answers1

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.

jbdv
  • 1,263
  • 1
  • 11
  • 18