F2Py with allocatable and assumed shape arrays

Question

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.

Answer 1

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.

8. Instead of

   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.

Answer 2

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.]

Answer 3

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.]

Answer 4

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