Equivalence between Fortran allocatable arrays and pointers

4
votes

I have a Fortran program with an allocatable array A as follows:

real, dimension(:,:) allocatable :: A
...
allocate(A(x0:x1;y0:y1))

This array is eventually passed as argument to a subroutine which looks like

subroutine my_subroutine(arr)
   real, dimension(x0:x1,y0:y1) :: arr
   ...
end subroutine my_subroutine

I wanted to replace the allocate statement of Fortran by a custom memory allocation function my_alloc implemented in a C library. I changed the first code sample into:

type(c_ptr) :: cptr
real, pointer, dimension(:,:) :: A
...
cptr = my_alloc(...)
call c_f_pointer(cptr,A,[x1-x0+1,y1-y0+1])

This works fine, except that by specifying extents instead of lower/upper bounds in the c_f_pointer function, I lose the original shape (x0:x1,y0:y1) of the array. But this is not a big problem: the pointer is passed as argument of the subroutine, the subroutine expects an array and considers the pointer as an array, with the proper bounds.

My real problem is: when I want to also rewrite the subroutine's code to have a pointer instead of an array.

subroutine my_subroutine(arr)
   real, pointer, dimension(x0:x1,y0:y1) :: arr
   ...
end subroutine my_subroutine

The code above doesn't work; gfortran says

Array pointer 'arr' at (1) must have a deferred shape

The following code can be compiled

subroutine my_subroutine(arr)
   real, pointer, dimension(:,:) :: arr
   ...
end subroutine my_subroutine

but it doesn't provide the bounds and the program crashes when I try to perform a loop from x0 to x1 and from y0 to y1.

How can I handle this case? Within the subroutine, I need fortran to know that arr is a pointer to an array shaped (x0:x1,y0;y1).

3
arrayspointersfortran

3 Answers

5
votes

Yes, this was a problem because of the limitation of c_f_pointer. As you have found the intrinsic c_f_pointer only supports bounds starting at index 1. People frequently state that Fortran is a one-indexed language but this not true. One indexing is only the default and Fortran has long supported declaring any starting bound that the programmer wants. So it was a step backwards that c_f_pointer forced you to use one indexing. But with Fortran 2003 there is a fix: pointer bounds remapping: 

arr (0:n-1) => arr

instead of 1:n, or whatever you wish.

Then pass the array to the subroutine and it will receive the intended bounds.

EDIT: improve demo program showing the difference between allocatables and pointers. The pointer passes the bounds of the array. A regular array passes the shape ... you can declare the first dimension in a subroutine, if you wish, and let the shape control the second.

module mysubs

implicit none

contains

subroutine testsub ( ptr, alloc, start, array )

   real, pointer, dimension (:) :: ptr
   real, dimension (:), intent (in) :: alloc
   integer, intent (in) :: start
   real, dimension (start:), intent (in) :: array

   write (*, *) "pointer in sub:", lbound (ptr, 1), ubound (ptr, 1)
   write (*, *) ptr

   write (*, *) "1st array in sub:", lbound (alloc, 1), ubound (alloc, 1)
   write (*, *) alloc

   write (*, *) "2nd array in sub:", lbound (array, 1), ubound (array, 1)
   write (*, *) array

   return

end subroutine testsub

end module mysubs


program test_ptr_assignment

use mysubs

implicit none

real, pointer, dimension(:) :: test
real, allocatable, dimension(:) :: alloc1, alloc2
real, allocatable, dimension(:) :: alloc1B, alloc2B

allocate ( test (1:5), alloc1 (1:5), alloc1B (1:5) )
test = [ 1.0, 2.0, 3.0, 4.0, 5.0 ]
alloc1 = test
alloc1B = test

write (*, *) "A:", lbound (test, 1), ubound (test, 1)
write (*, *) test

call testsub (test, alloc1, 1, alloc1B )

test (0:4) => test
allocate ( alloc2 (0:4), alloc2B (0:4) )
alloc2 = test
alloc2B = test

write (*, *)
write (*, *) "B:", lbound (test, 1), ubound (test, 1)
write (*, *) test

call testsub (test, alloc2, 0, alloc2B)

stop

end program test_ptr_assignment
1
votes

You can use lbound and ubound built-in functions to find out the bounds of the array within the subroutine, e.g. 

program test
  real, dimension(:,:), pointer :: A

  allocate(A(3:5,7:8))
  A = 1
  call my_print(A)

contains

  subroutine my_print(X)
    integer :: i,ml,mu
    real, dimension(:,:), pointer :: X
    ml = lbound(X,1)
    mu = ubound(X,1)
    do i = ml,mu
       write(*,*) X(i,:)
    end do
  end subroutine my_print

end program test
0
votes

The defered shape the compiler insists on means, you can not declare upper bounds of the dummy arguments. But you can declare the lower ones, and the upper ones you will get automatically:

subroutine my_subroutine(arr)
   real, dimension(x0:,y0:) :: arr
   ...
end subroutine my_subroutine

It's true also for allocatable arrays.