18
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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 :) $\endgroup$ – Aron Ahmadia Apr 24 '13 at 19:41
23
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ As far as I know, f2py does not rely on ISO C bindings (its primary target is Fortran 77 and Fortran 90 code). $\endgroup$ – Aron Ahmadia Apr 24 '13 at 21:31
  • $\begingroup$ 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! $\endgroup$ – MRocklin Apr 24 '13 at 21:49
  • $\begingroup$ 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. $\endgroup$ – Ondřej Čertík Apr 24 '13 at 23:28
6
$\begingroup$

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.]
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Alleo Mar 14 '16 at 13:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.