[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