# F2Py with allocatable and assumed shape arrays

I would like to use f2py with modern Fortran. In particular I'm trying to get the following basic example to work. This is the smallest useful example I could generate.

! alloc_test.f90
subroutine f(x, z)
implicit none

! Argument Declarations !
real*8, intent(in) ::  x(:)
real*8, intent(out) :: z(:)

! Variable Declarations !
real*8, allocatable :: y(:)
integer :: n

! Variable Initializations !
n = size(x)
allocate(y(n))

! Statements !
y(:) = 1.0
z = x + y

deallocate(y)
return
end subroutine f


Note that n is inferred from the shape of input parameter x. Note that y is allocated and deallocated within the body of the subroutine.

When I compile this with f2py

f2py -c alloc_test.f90 -m alloc


And then run in Python

from alloc import f
from numpy import ones
x = ones(5)
print f(x)


I get the following error

ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (-1,)


So I go and create and edit the pyf file manually

f2py -h alloc_test.pyf -m alloc alloc_test.f90


Original

python module alloc ! in
interface  ! in :alloc
subroutine f(x,z) ! in :alloc:alloc_test.f90
real*8 dimension(:),intent(in) :: x
real*8 dimension(:),intent(out) :: z
end subroutine f
end interface
end python module alloc


Modified

python module alloc ! in
interface  ! in :alloc
subroutine f(x,z,n) ! in :alloc:alloc_test.f90
integer, intent(in) :: n
real*8 dimension(n),intent(in) :: x
real*8 dimension(n),intent(out) :: z
end subroutine f
end interface
end python module alloc


Now it runs but the values of the output z are always 0. Some debug printing reveals that n has the value 0 within the subroutine f. I assume that I'm missing some f2py header magic to manage this situation properly.

More generally what is the best way to link the above subroutine into Python? I'd strongly prefer not to have to modify the subroutine itself.

• Matt, are you familiar with Ondrej Certik's best practices guide, specifically, the Interfacing with Python section? We've been discussing a similar interfacing issue for PyClaw and haven't resolved it yet at this point either :) – Aron Ahmadia Apr 24 '13 at 19:41
• One just needs to specify the shape of z. All else is fine. See mine or @Prometheous answer below. The currently accepted answer is an interesting read but doesn't relate much to the issue here. – Jonatan Öström Jul 5 '20 at 21:12

I am not super familiar with f2py internals, but I am very familiar with wrapping Fortran. F2py just automates some or all of the things below.

1. You first need to export to C using the iso_c_binding module, as described for example here:

http://fortran90.org/src/best-practices.html#interfacing-with-c

Disclaimer: I am the main author of the fortran90.org pages. This is the only platform and compiler independent way of calling Fortran from C. This is F2003, so these days there is no reason to use any other way.

2. You can only export/call arrays with full length specified (explicit shape), that is:

integer(c_int), intent(in) :: N
real(c_double), intent(out) :: mesh(N)


but not assume shape:

real(c_double), intent(out) :: mesh(:)


That is because the C language does not support such arrays itself. There is talk to include such support in either F2008 or later (I am not sure), and the way it would work is through some supporting C data structures, as you need to carry shape information about the array.

In Fortran, you should mainly use the assume shape, only in special cases you should use explicit shape, as described here:

http://fortran90.org/src/best-practices.html#arrays

That means, that you need to write a simple wrapper around your assume shape subroutine, that will wrap things into explicit shape arrays, per my first link above.

3. Once you have a C signature, just call it from Python in any way you like, I use Cython, but you can use ctype, or C/API by hand.

4. The deallocate(y) statement is not needed, Fortran deallocates automatically.

http://fortran90.org/src/best-practices.html#allocatable-arrays

5. real*8 should not be used, but rather real(dp):

http://fortran90.org/src/best-practices.html#floating-point-numbers

6. The statement y(:) = 1.0 is assigning 1.0 in single precision, so the rest of digits will be random! This is a common pitfall:

http://fortran90.org/src/gotchas.html#floating-point-numbers

You need to use y(:) = 1.0_dp.

7. Instead of writing y(:) = 1.0_dp, you can just write y = 1, that's it. You can assign integer to a floating point number without losing accuracy, and you don't need to put the redundant (:) in there. Much simpler.

y = 1
z = x + y


just use

z = x + 1


and don't bother with the y array at all.

9. You don't need the "return" statement at the end of the subroutine.

10. Finally, you should probably be using modules, and just put implicit none on the module level and you don't need to repeat it in each subroutine.

Otherwise it looks good to me. Here is the code following the suggestions 1-10 above::

module test
use iso_c_binding, only: c_double, c_int
implicit none
integer, parameter :: dp=kind(0.d0)

contains

subroutine f(x, z)
real(dp), intent(in) ::  x(:)
real(dp), intent(out) :: z(:)
z = x + 1
end subroutine

subroutine c_f(n, x, z) bind(c)
integer(c_int), intent(in) :: n
real(c_double), intent(in) ::  x(n)
real(c_double), intent(out) :: z(n)
call f(x, z)
end subroutine

end module


It shows the simplified subroutine as well as a C wrapper.

As far as f2py, it probably tries to write this wrapper for you and fails. I am also not sure whether it is using the iso_c_binding module. So for all these reasons, I prefer to wrap things by hand. Then it's exactly clear what is happening.

• As far as I know, f2py does not rely on ISO C bindings (its primary target is Fortran 77 and Fortran 90 code). – Aron Ahmadia Apr 24 '13 at 21:31
• I knew I was being a bit dumb with y but I wanted to make something was allocated (my actual code has non-trivial allocations). I did not know about many of the other points though. Looks like I should go look into some sort of Fortran90 best practices guide.... Thanks for the thorough answer! – MRocklin Apr 24 '13 at 21:49
• Note that using today's Fortran compilers, you wrap F77 in exactly the same way --- by writing a simple iso_c_binding wrapper and call the legacy f77 subroutine from it. – Ondřej Čertík Apr 24 '13 at 23:28

All you have to do is the following:

!alloc_test.f90
subroutine f(x, z, n)
implicit none

! Argument Declarations !
integer :: n
real*8, intent(in) ::  x(n)
real*8, intent(out) :: z(n)

! Variable Declarations !
real*8, allocatable :: y(:)

! Variable Initializations !
allocate(y(n))

! Statements !
y(:) = 1.0
z = x + y

deallocate(y)
return
end subroutine f


Though the size of array x and z is now passed as an explicit argument, f2py makes the argument n optional. Following is the docstring of the function as it appears to python:

Type:       fortran
String Form:<fortran object>
Docstring:
f - Function signature:
z = f(x,[n])
Required arguments:
x : input rank-1 array('d') with bounds (n)
Optional arguments:
n := len(x) input int
Return objects:
z : rank-1 array('d') with bounds (n)


Importing and calling it from python:

from alloc import f
from numpy import ones
x = ones(5)
print f(x)


gives the following output:

[ 2.  2.  2.  2.  2.]

• Is there a way to use some non-trivial expression as size? For instance, I pass n and want to get an array of size 2 ** n. So far I have to pass also 2 ** n as a separate argument. – Alleo Mar 14 '16 at 13:32

I experienced the same issue with assumed-shape dummy arrays, and the approach by Jonatan Ostrom (here) worked for me also. Because f2py creates a fresh variable for dummy arguments with intent(out), an explicit shape information seems necessary for such arguments. This is in contrast to dummy arrays with intent(in) and intent(inout) or with no intent, for which Numpy ndarrays (created on the Python side) are passed to Fortran, so f2py can probably create wrappers with sufficient shape information for assumed-shape arrays.

mycode.f90:

module fmod
implicit none
contains

!! Use assumed-shape dummy args.
subroutine fsub( a_in, a_inout, a_out, a_none )

real(8), intent(in)    :: a_in(:)
real(8), intent(inout) :: a_inout(:)

real(8), intent(out)   :: a_out( size(a_in) )  !! OK
!real(8), intent(out)  :: a_out(:)             !! fails

real(8)                :: a_none(:)   !! no intent

print *
print *, "fort| setting 777 to a_out(:)"
a_out(:) = 777

print *
print *, "fort| a_in    = ", a_in
print *, "fort| a_inout = ", a_inout
print *, "fort| a_out   = ", a_out
print *, "fort| a_none  = ", a_none

print *
print *, "fort| multiplying a_none(:) by 100"
a_none(:) = a_none(:) * 100
end

end module


main.py:

import numpy as np

# Load a Python module created by f2py.
import pymod

# Show info.
print( pymod.__doc__ )
print( pymod.fmod.__doc__ )

# Prepare arrays on the Python side.
b_in    = np.ones( 3, dtype='d' ) * 1.0
b_inout = np.ones( 3, dtype='d' ) * 2.0
b_none  = np.ones( 3, dtype='d' ) * 3.0

print( "\nBefore:" )
print( "py| b_in    = ", b_in )
print( "py| b_inout = ", b_inout )
print( "py| b_none  = ", b_none )

# Call Fortran.
b_out = pymod.fmod.fsub( a_in = b_in,
a_inout = b_inout,
a_none = b_none )

print( "\nAfter:" )
print( "py| b_in    = ", b_in )
print( "py| b_inout = ", b_inout )
print( "py| b_out   = ", b_out )
print( "py| b_none  = ", b_none )


Build (f2py):

python3.8 -m numpy.f2py -c mycode.f90 -m pymod


Run:

\$ py main.py

This module 'pymod' is auto-generated with f2py (version:2).
Functions:
Fortran 90/95 modules:
fmod --- fsub().
a_out = fsub(a_in,a_inout,a_none)

Wrapper for fsub.

Parameters
----------
a_in : input rank-1 array('d') with bounds (f2py_a_in_d0)
a_inout : in/output rank-1 array('d') with bounds (f2py_a_inout_d0)
a_none : input rank-1 array('d') with bounds (f2py_a_none_d0)

Returns
-------
a_out : rank-1 array('d') with bounds (size(a_in))

Before:
py| b_in    =  [1. 1. 1.]
py| b_inout =  [2. 2. 2.]
py| b_none  =  [3. 3. 3.]

fort| setting 777 to a_out(:)

fort| a_in    =    1.0000000000000000        1.0000000000000000        1.0000000000000000
fort| a_inout =    2.0000000000000000        2.0000000000000000        2.0000000000000000
fort| a_out   =    777.00000000000000        777.00000000000000        777.00000000000000
fort| a_none  =    3.0000000000000000        3.0000000000000000        3.0000000000000000

fort| multiplying a_none(:) by 100

After:
py| b_in    =  [1. 1. 1.]
py| b_inout =  [2. 2. 2.]
py| b_out   =  [777. 777. 777.]
py| b_none  =  [300. 300. 300.]


I don't think the accepted answer actually answers the question. The allocations are fine. The single problem is that the output is assumed shape. This is not possible since there will be no output-array passed from Python from which the Fortran routine can get the shape. So change

real*8, intent(out) :: z(:)


into

real*8, intent(out) :: z(size(x))