Fortran 可分配数组和指针之间的等价性

Question

我有一个带有可分配数组A 的 Fortran 程序，如下所示：

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

这个数组最终作为参数传递给一个看起来像

的子例程

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

我想用 C 库中实现的自定义内存分配函数 my_alloc 替换 Fortran 的 allocate 语句。我将第一个代码示例更改为：

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

这很好用，除了通过在 c_f_pointer 函数中指定范围而不是下限/上限，我失去了数组的原始形状 (x0:x1,y0:y1)。但这不是一个大问题：指针作为子程序的参数传递，子程序需要一个数组并将指针视为一个数组，具有适当的边界。

我真正的问题是：当我还想重写子程序的代码以使用指针而不是数组时。

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

上面的代码不起作用； gfortran 说

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

以下代码可以编译

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

但是当我尝试执行从 x0 到 x1 以及从 y0 到 y1 的循环时，它没有提供边界并且程序崩溃。

我该如何处理这种情况？在子例程中，我需要 fortran 知道 arr 是指向形状为 (x0:x1,y0;y1) 的数组的指针。

Answer 1

是的，由于 c_f_pointer 的限制，这是一个问题。正如您发现的那样，内在的 c_f_pointer 仅支持从索引 1 开始的边界。人们经常说 Fortran 是一种单索引语言，但事实并非如此。一个索引只是默认设置，Fortran 长期以来一直支持声明程序员想要的任何起始界限。因此，c_f_pointer 迫使您使用一个索引是一种倒退。但是对于 Fortran 2003 有一个修复：指针边界重新映射：

arr (0:n-1) => arr

而不是 1:n，或者任何你想要的。

然后将数组传递给子例程，它将接收预期的边界。

编辑：改进演示程序，显示可分配对象和指针之间的区别。指针传递数组的边界。常规数组传递形状...如果您愿意，您可以在子例程中声明第一个维度，并让形状控制第二个维度。

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

Answer 2

您可以使用lbound 和ubound 内置函数来找出子程序中数组的边界，例如

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

Answer 3

编译器坚持的延迟形状意味着，你不能声明虚拟参数的上限。但是你可以声明低的，高的你会自动得到：

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

对于可分配数组也是如此。