Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

how to write wrapper for 'allocate'

Tags:

fortran

I am trying to write a wrapper for 'allocate' function, i.e. function which receives an array and dimensions, allocates memory and returns allocated array. The most important thing is that the function must work with arrays of different rank. But I have to explicitly state rank of array in function interface, and in this case code only compiles if I pass arrays of certain rank as a parameter. For example, this code does not compile:

module memory_allocator
contains 

  subroutine memory(array, length)
    implicit none

    real(8), allocatable, intent(out), dimension(:) :: array
    integer, intent(in) :: length

    integer :: ierr

    print *, "memory: before: ", allocated(array)

    allocate(array(length), stat=ierr)
    if (ierr /= 0) then
      print *, "error allocating memory: ierr=", ierr
    end if

    print *, "memory: after: ", allocated(array)

  end subroutine memory

  subroutine freem(array)
    implicit none

    real(8), allocatable, dimension(:) :: array

    print *, "freem: before: ", allocated(array)
    deallocate(array)
    print *, "freem: after: ", allocated(array)

  end subroutine freem

end module memory_allocator

program alloc
  use memory_allocator
  implicit none

  integer, parameter :: n = 3
  real(8), allocatable, dimension(:,:,:) :: foo
  integer :: i, j, k

  print *, "main: before memory: ", allocated(foo)
  call memory(foo, n*n*n)
  print *, "main: after memory: ", allocated(foo)

  do i = 1,n
    do j = 1,n
      do k = 1, n
        foo(i, j, k) = real(i*j*k)
      end do
    end do
  end do

  print *, foo

  print *, "main: before freem: ", allocated(foo)
  call freem(foo)  
  print *, "main: after freem: ", allocated(foo)

end program alloc

Compilation error:

gfortran -o alloc alloc.f90 -std=f2003
alloc.f90:46.14:

  call memory(foo, n*n*n)
              1
Error: Rank mismatch in argument 'array' at (1) (1 and 3)
alloc.f90:60.13:

  call freem(foo)  
             1
Error: Rank mismatch in argument 'array' at (1) (1 and 3)  

Is there any way of implementing such wrapper?..

Thanks!

like image 500
robusta Avatar asked Feb 13 '10 11:02

robusta


2 Answers

This can be done via a generic interface block. You have to create procedures for each rank that you want to handle, e.g., memory_1d, memory_2d, ... memory_4d. (Obviously a lot of cut & pasting.) Then you write a generic interface block that gives all of these procedures the alternative name memory as a generic procedure name. When you call memory, the compiler distinguishes which memory_Xd should be called based on the rank of the argument. The same for your freem functions.

This is how intrinsic functions such as sin have long worked -- you can call sin with a real arguments of various previsions, or with a complex argument, and the compiler figures out with actual sin function to call. In really old FORTRAN you had to use different names for the different sin functions. Now modern Fortran you can setup the same thing with your own routines.

Edit: adding a code example demonstrating the method & syntax:

module double_array_mod

   implicit none

   interface double_array
      module procedure double_vector
      module procedure double_array_2D
   end interface double_array

   private  ! hides items not listed on public statement 
   public :: double_array

contains

   subroutine double_vector (vector)
      integer, dimension (:), intent (inout) :: vector
      vector = 2 * vector
   end subroutine double_vector

   subroutine double_array_2D (array)
      integer, dimension (:,:), intent (inout) :: array
      array = 2 * array
   end subroutine double_array_2D

end module double_array_mod


program demo_user_generic

   use double_array_mod

   implicit none

   integer, dimension (2) :: A = [1, 2]
   integer, dimension (2,2) :: B = reshape ( [11, 12, 13, 14], [2,2] )
   integer :: i

   write (*, '( / "vector before:", / 2(2X, I3) )' )  A
   call double_array (A)
   write (*, '( / "vector after:", / 2(2X, I3) )' )  A

   write (*, '( / "2D array before:" )' )
   do i=1, 2
      write (*, '( 2(2X, I3) )' )  B (i, :)
   end do
   call double_array (B)
   write (*, '( / "2D array after:" )' )
   do i=1, 2
      write (*, '( 2(2X, I3) )' )  B (i, :)
   end do   

   stop
end program demo_user_generic
like image 196
M. S. B. Avatar answered Nov 01 '22 16:11

M. S. B.


subroutine memory(array, length) has as it first dummy parameter 1-dimensional array (real(8), allocatable, intent(out), dimension(:) :: array).

Calling this subroutine from your main program with 3-dimensional array foo (real(8), allocatable, dimension(:,:,:) :: foo) is error obviously. And this is what compiler actually said.

If you really need such subroutines write one pair memory/freem subroutines for each array of different dimension - one subroutines pair for 1-dimensional array, another for 2-dimensional array, etc.

By the way, memory subroutines will be different in general because in order to allocate n-dimensional array you need to pass n extents to above-mentioned subroutine.

like image 1
Wildcat Avatar answered Nov 01 '22 15:11

Wildcat