Equivalence of Fortran Distributed Arrays and Pointers

I have a Fortran program with allocated array A as follows:

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

This array is eventually passed as an argument to the routine, which looks like

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

I wanted to replace the allocate for Fortran with the special memory allocation function my_alloc implemented in the C library. I changed the first code example to:

 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 great, except that by specifying extents instead of the 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 an argument to the subroutine, the subroutine expects an array and considers the pointer as an array with the corresponding boundaries.

My real problem: when I want to also rewrite the routine code with 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 does not 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 does not provide boundaries, and the program crashes when I try to execute a cycle from x0 to x1 and from y0 to y1.

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

+4
source share
3 answers

Yes, it was a problem due to the c_f_pointer restriction. As you have discovered, the built-in c_f_pointer only supports borders starting with index 1. People often claim that Fortran is an indexed language, but this is not so. One indexing is only the default, and Fortran has long supported the declaration of any initial boundary that the programmer wants. So, it was a step backward that c_f_pointer made you use one indexing. But with Fortran 2003 there is a fix: overriding pointers:

 arr (0:n-1) => arr 

instead of 1: n or whatever.

Then pass the array to the subroutine and it will get the assigned boundaries.

EDIT: Improve the demo showing the difference between allocatables and pointers. The pointer passes the boundaries of the array. A regular array passes the form ... you can declare the first dimension in the routine, if you want, and let the shape control be 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 
+5
source

You can use the built-in functions lbound and ubound to determine the boundaries of an array in a routine, for example.

 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 
+1
source

The deferred form, which the compiler insists on means, cannot declare the upper bounds of dummy arguments. But you can declare the lower ones, and the upper ones automatically:

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

This is also true for distributed arrays.

0
source

Source: https://habr.com/ru/post/1397467/


All Articles