[Numpy-discussion] f2py with allocatable arrays

Jim Vickroy jim.vickroy@noaa....
Tue Jul 3 06:09:24 CDT 2012


On 7/2/2012 7:17 PM, Casey W. Stark wrote:
> Hi numpy.
>
> Does anyone know if f2py supports allocatable arrays, allocated inside 
> fortran subroutines? The old f2py docs seem to indicate that the 
> allocatable array must be created with numpy, and dropped in the 
> module. Here's more background to explain...
>
> I have a fortran subroutine that returns allocatable positions and 
> velocities arrays. I wish I could get rid of the allocatable part, but 
> you don't know how many particles it will create until the subroutine 
> does some work (it checks if each particle it perturbs ends up in the 
> domain).
>
> module zp
>   implicit none
>   contains
>   subroutine ics(..., num_particles, particle_mass, positions, velocities)
>     use data_types, only : dp
>     implicit none
>     ... inputs ...
>     integer, intent(out) :: num_particles
>     real (kind=dp), intent(out) :: particle_mass
>     real (kind=dp), intent(out), dimension(:, :), allocatable :: 
> positions, velocities
>     ...
>   end subroutine
> end module
>
> I tested this with a fortran driver program and it looks good, but 
> when I try with f2py, it cannot compile. It throws the error "Error: 
> Actual argument for 'positions' must be ALLOCATABLE at (1)". I figure 
> this has something to do with the auto-generated "*-f2pywrappers2.f90" 
> file, but the build deletes the file.
>
> If anyone knows an f2py friendly way to rework this, I would be happy 
> to try. I'm also fine with using ctypes if it can handle this case.
>
> Best,
> Casey
>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
Hi,

... not sure what version of numpy you are using but it can be done with 
the f2py included with numpy 1.6.x.  Here is a (hopefully intelligible) 
code fragment from a working application:



module Grid_
    ! a representation of a two-dimensional, rectangular, (image) sensor 
grid
    ! -- this component is used in the "smoothing" portion of the 
Thematic MaP solution algorithm (smoothness prior probabilities)

    use runtime_
    implicit none
    integer, private, parameter                  :: begin_ = 1, end_ = 
2, row_ = 1, column_ = 2
    integer, private,     dimension(begin_:end_) :: rows__ = (/0,0/), 
columns__ = (/0,0/) ! grid rows and columns extent
    integer, allocatable, dimension(:,:)         :: pixel_neighbors
       ! the neighbors of a specified pixel -- ((row,column), ..., 
(row,column))
       !    + this array is managed (allocated, deallocated, populated) 
by the neighbors_of() procedure
       !    + when allocated, the first dimension is always 2 -- pixel 
(row,column)
       !    + it would be preferable for this array to be the return 
result, of procedure neighbors_of(),
       !      but F2PY does not seem to support an allocatable array as 
a function result

    contains

       [snip]

       function neighbors_of (pixel) result(completed)
          ! determines all second-order (grid) neighbors of *pixel*
          !    -- returns .true. if successful
          !    -- use the error_diagnostic() procedure if .false. is 
returned
          !    -- upon successful completion of this procedure, 
Grid_.pixel_neighbors contains the neighbors of *pixel*
          !    -- *pixel* may not be on the grid; in this case, 
Grid_.pixel_neighbors is not allocated and .false. is returned
          ! integer, dimension(row_:column_), intent(in) :: pixel ! 
pixel to be considered ... f2py does not support this
          integer, dimension(2), intent(in) :: pixel  ! pixel under 
consideration (row,column)
          logical                           :: completed
          ! integer, dimension(begin_:end_) :: neighbor_rows, 
neighbor_columns ... f2py does not support this
          ! integer, dimension(row_:column_)   :: neighbor ... f2py does 
not support this
          integer, dimension(2)             :: neighbor_rows, 
neighbor_columns ! each is: (begin,end)
          integer, dimension(2)             :: neighbor ! (row,column)
          integer                           :: row, column, code, count
          character (len=100)               :: diagnostic
          completed = .false.
          count     = 0 ! *pixel* has no neighbors
          if (allocated (pixel_neighbors)) deallocate (pixel_neighbors)
          if (is_interior (pixel)) then
             count                     = 8 ! interior pixels have eight, 
second-order neighbors
             neighbor_rows    (begin_) = pixel(row_)-1
             neighbor_rows    (end_)   = pixel(row_)+1
             neighbor_columns (begin_) = pixel(column_)-1
             neighbor_columns (end_)   = pixel(column_)+1
          else if (is_border (pixel)) then
             count = 5 ! non-corner, border pixels have five, 
second-order neighbors, but ...
             if (is_corner (pixel)) count = 3 ! corner pixels have 
three, second-order neighbors
             neighbor_rows    (begin_) = max (pixel(row_)-1, rows__(begin_))
             neighbor_rows    (end_)   = min (pixel(row_)+1, rows__(end_))
             neighbor_columns (begin_) = max (pixel(column_)-1, 
columns__(begin_))
             neighbor_columns (end_)   = min (pixel(column_)+1, 
columns__(end_))
             end if
          if (count > 0) then
             allocate (pixel_neighbors(row_:column_,count), stat=code, 
errmsg=diagnostic)
                if (code /= 0) then
                   call set_error (code,diagnostic)
                   return
                   end if
             count = 0
             do row = neighbor_rows(begin_), neighbor_rows(end_)
                do column = neighbor_columns(begin_), neighbor_columns(end_)
                   neighbor(row_) = row;  neighbor(column_) = column
                   if (neighbor(row_) == pixel(row_) .and. 
neighbor(column_) == pixel(column_)) cycle ! neighbor is pixel
                   count = count + 1;   pixel_neighbors(:,count) = neighbor
                   end do ! columns
                end do ! rows
             completed = .true.
             end if
          end function neighbors_of

    end module Grid_


-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/numpy-discussion/attachments/20120703/1c655ce2/attachment-0001.html 


More information about the NumPy-Discussion mailing list