Fortran 90 自动(?)分配可分配项

Fortran 90 auto(?) allocation of allocatables

program test

    implicit none
    real(dp),allocatable :: a(:), b(:)

    allocate (a(10))

    a = 1.d0
(*) b = a

end program

在上面的代码中,我设置了两个可分配项 - ab - 并且在代码中只分配了 a。我预计代码无法编译,但它编译得很好并且看起来运行良好,而下面的代码,它做了类似的工作,显示 SEGFAULT。

program test_fail

   implicit none
   real(dp),allocatable :: a(:), b(:)
   integer              :: i

   allocate ( a(10) )

   a = 1.d0

   do i = 1, 10
      b(i) = a(i)
   enddo

end program

前面的代码是自动分配b的意思好理解吗?

subroutine test_sub(a)

   implicit none

   real(dp),intent(inout) :: a(:,:)

   ...

end program

并且在上面的带有整形数组输入的子例程中,是否可以理解代码自动检测输入数组的大小a,并在子例程中分配自己的数组并释放它返回到它的上框架?

最后,当我将一个数组复制到另一个数组时,哪个更快?

program speed1

  implicit none
  real(dp), allocatable :: a(:,:,:), b(:,:,:)

  allocate( a(10000,10000,10000) )

  a = 1.d0

  b = a

end program
program speed2

  implicit none
  real(dp), allocatable :: a(:,:,:), b(:,:,:)
  integer               :: i, j, k

  allocate( a(10000,10000,10000), b(10000,10000,10000) )

  a = 1.d0

  do i = 1,10000
    do j = 1,10000
      do k = 1,10000
        b(i,j,k) = a(i,j,k)
      enddo
    enddo
  enddo

end program

我没有太多时间,所以这可能会被压缩。还要注意,正如上面评论中指出的那样(为了突出和后代而在此处复制)代码的有效性以及此答案的有效性取决于版本(有关详细信息,请参见 this question)。

如您所见,语句

b = a

会自动将数组 b 分配为与 a 相同的大小(和形状),并为其元素赋予与 a 的元素相同的值。这一切都符合标准。然而,声明

b(i) = a(i)

左边没有 array,它有一个 array element,在这种情况下 Fortran 不会自动分配数组。不幸的是,但这是(我相信)根据语言标准不需要编译器来发现错误的情况之一——这就是为什么你会在 运行 时间发现错误。您的编译器的行为似乎与我最近体验过的其他编译器一样 - 实际上,没有元素 b(1) 可以将 a(1) 的值分配给 etc .

至于子程序中的'allocation'语句

real(dp),intent(inout) :: a(:,:) 

这不太一样。这里 aassumed-shape 并且简单地假定传递给例程的相应参数的形状(和值)。该标准没有说明这是如何完成的。大多数编译器不会复制数组(出于性能原因,这对 Fortran 程序员通常很重要)但有些可能,并且很容易构造大多数编译器将复制作为参数传递的数据结构的示例。

至于最后两个小码中哪个更快 - 你为什么不告诉我们?

High Performance Mark 的回答是正确且有帮助的,但以稍微不同的方式重述一些要点可能会有一些好处。

如前所述,在该答案和链接问题中,赋值如

b = a

for b 未分配取决于是否使用 Fortran 2003+ 规则。但是,如果未分配 b,则分配给 b 的元素,如

b(i) = a(i)

总是错误(直到当前语言修订)。每当数组元素出现在赋值的左侧(或变量引用中)时,下标的值必须在数组的范围内 - 但未分配的可分配数组没有定义的范围.1

这与 Matlab 等语言形成鲜明对比,在这些语言中,如果存在超出数组大小当前范围的赋值,数组就会“增长”。 Fortran (2003+) 整个数组分配是一个例外,其中(重新)分配发生在整个数组完全位于右侧的情况下。在这样的整个数组分配中,所需的大小是已知的;对于逐个元素分配,最终分配最初是未知的。

关于子例程 test_sub 的伪参数 a 有一点要添加到 High Performance Mark 的回答中: 意味着子例程中的 a 可能共享调用子例程的主程序中的数组的某些方面。这确实可能包括之前给出的任何 allocation/storage。

不过,为了完整起见:Fortran 是其中一种未逐行定义语义的语言。子例程的 ... 隐藏了其他人可能感兴趣的内容。在 ... 中,我们可以更改对 a 的解释。

subroutine  test_sub(a)
   real(dp),intent(inout) :: a(:,:)
   allocatable :: a  ! or "pointer :: a"
end program

突然a不是假定形状而是延迟形状。我不会在这里讨论后果,但可以在其他问题和答案中找到它们。

最后,关于赋值速度,我们来比较整个数组赋值和遍历数组的所有元素(左侧正确分配)。 Fortran 的规则保证,对于左侧的数组,赋值是“对相应的数组元素逐个元素执行”,但“处理器可以按任何顺序逐个元素执行赋值”。

您的赋值(在大多数情况下)会产生相同的效果,但编译器可能会重新排序或利用它认为合适的 vectorization/parallelization。然而,同样地,它可以对循环进行完全相同的分析(包括重新排序循环,正如 Ian Bush 的评论所建议的那样可能有帮助,或者响应诸如 OpenMP 指令之类的事情)。如果你想知道你的情况下性能的精确细节,你必须进行测试,但如果你用手工制作的循环进行微优化,你最终会惹恼别人,就像你肯定会惹恼其他人一样-数组赋值。


1对于非常感兴趣的一方,请参阅 Fortran 2018 8.5.8.4 和 9.5.3.1。

让我扩展我的评论作为答案。

1st 代码被截断编译,因为 b 在分配 b=a 时采用 a 的形式,这是 Fortran 标准特征

2nd 片段是错误的:无法索引未分配的数组,它没有形状或条目。因此出现段错误。

(在这种简单的情况下,编译器可以很容易地在编译时检测到这个错误并给出错误。但是可分配数组可以在子例程和 allocated/deallocated 之间传递,因此很难检测到错误一般情况。我想这就是为什么没有完成此检查的原因。)

3rd 片段:subroutine/function 参数中的 assumed-shape 数组 a(:,:),采用传递给常规。它(通常)是对同一内存的引用。

第 4 和第 5 片段:Fortran 数组是 column-major 意味着最左边的索引应该是最内层的循环等以实现最快的执行。喜欢:

do k = 1, n
    do j = 1, n
        do i = 1, n
            b(i,j,k) = a(i,j,k)

但是,在简单的情况下,编译器会在应用优化时修复此问题(例如 gfortran -O3)。

由于 b 未使用,可能是编译器在这两种情况下都删除了复制操作。

复制很少成为瓶颈,所以无论性能如何,我都会使用 b=a,因为它清晰而简短。无论如何,上面的循环和 a=b 可能只等同于编译器。