Fortran 90自动分配可分配变量的内存

Question

program test

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

    allocate (a(10))

    a = 1.d0
(*) b = a

end program

在上面的代码中，我设置了两个可分配变量 - a 和 b - 但只在代码中分配了 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

- Sangjun Lee

---

由于您的循环顺序错误，b=a可能会更快（假设您在speed2中修复了分配问题）。如果您将循环顺序反过来，我希望它们需要相同的时间，但是如果代码瓶颈是将数组清零，则这是非常奇怪的代码！ - Ian Bush

请注意，这是一个重要的时刻，需要注意它是Fortran 90还是自由格式（在F90中引入）：关于b赋值的规则在Fortran 2003中发生了改变，因此行为可能取决于编译器版本和/或编译时标志。 - francescalus

@IanBush 我应该在 speed2 中修复什么？我已经尝试了几种不同的 i,j,k 组合顺序，但它们都比 b = a 语句慢。如果是二维情况，使用 Fortran 的列主序排序，我期望 j,i 的顺序会更快。但对于三维情况，哪种 i,j,k 的顺序是 Fortran 中最优的呢？在我的测试中，i,j,k 是最快的。 - Sangjun Lee

数组在分配时会自动分配空间，因为此时数组的形状已知。无法对未分配空间的数组进行索引，这会导致段错误。在Fortran中，循环顺序应该是最左边的索引在最内层循环中，例如 do k; do j; do i; arr(i,j,k)。但是，只需使用 b=a 就可以简单明了地完成操作，而且从编译器的角度来看，这可能等效于正确排序的循环。 - Jonatan Öström

@JonatanÖström 谢谢您清晰的解释 :) - Sangjun Lee

Answer 1

我时间不多，所以这可能会比较简略。请注意，正如上面的评论所指出的（为了突出和永久性而复制在此处），代码的有效性以及本答案的有效性取决于版本（有关更多详细信息，请参见此问题）。

正如您发现的那样，语句

b = a

这将自动分配一个数组 b，使其大小（和形状）与数组a相同，并赋予其元素与数组a中的元素相同的值。这符合标准。然而，该语句

b(i) = a(i)

这里的左侧不是一个数组，而是一个数组元素，在这种情况下Fortran不会自动分配一个数组，这很不幸，但我认为这是编译器不需要按照语言标准检查错误的情况之一，这就是你在运行时发现错误的原因。你的编译器似乎与我最近使用过的其他编译器表现一致——实际上并没有元素b(1)来赋值给a(1)等。

至于带有语句的子程序中的“分配”：

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

这并不完全相同。这里的a是假定形状，仅仅假定与传递给程序的相关参数具有相同的形状（和值）。标准对于如何实现此操作保持沉默。大多数编译器不会复制数组（出于通常对Fortran程序员很重要的性能原因），但有些编译器可能会复制，并且可以构造出示例，在这些示例中，大多数编译器将复制作为参数传递的数据结构。

至于你最后两个代码片段哪个更快 - 为什么不告诉我们呢？

- High Performance Mark

Answer 2

High Performance Mark的答案是正确和有帮助的，但稍微以略微不同的方式重新表述一些要点可能会有一些优势。

如先前在那个答案中和链接的问题中提到的那样，像这样的赋值

b = a

如果未使用Fortran 2003+规则，b是否被分配取决于情况。但是，如果b未被分配，则对b元素的赋值会导致错误，如下所示：

b(i) = a(i)

一直以来，在编程语言的修订中，当一个数组元素出现在赋值的左侧（或者在变量引用中），如果下标的值超出了数组的范围，那么这种做法都是错误的。但是，一个未分配的可分配数组没有定义的边界。¹

这与Matlab等语言不同，如果对数组进行超出当前大小范围的赋值，数组会“增长”。Fortran（2003+）整个数组赋值是一个例外，因为整个数组在右侧完整地接收（重新）分配。在这样的整个数组赋值中，所需的大小是已知的；而在逐个元素赋值时，最终的分配最初是未知的。

关于子程序test_sub的虚拟参数a，有一点需要补充High Performance Mark的答案：argument association意味着子程序中的a可能与调用子程序的主程序中的某个数组共享某些方面。这确实可能包括之前给定的任何分配/存储。

为了完整性，Fortran是其中一种语言，其语义不是逐行定义的。 子程序中的...隐藏了其他人可能感兴趣的事情。 在...内部，我们可以更改对a的解释。

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

突然之间，a 不再是假定形状，而是 延迟形状。我不会在此讨论其后果，但可以在其他问题和答案中找到。

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

你的赋值（在大多数情况下）将产生相同的效果，但编译器可能会重新排序或利用向量化/并行化。同样地，它也可以对循环进行完全相同的分析（包括重新排序循环，如 Ian Bush 的评论建议可能有帮助，或者响应诸如 OpenMP 指令之类的事情）。如果你想知道在你的情况下性能的精确细节，你就必须进行测试，但如果你使用手工制作的循环进行微观优化，你最终肯定会惹恼某个人，就像你使用整个数组赋值一样。

---

¹对于深入了解的人，请参阅Fortran 2018 8.5.8.4和9.5.3.1。

- francescalus

Answer 3

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

第一段代码可以编译，因为在赋值b=a时，b采用了a的形状，这是Fortran标准特性。

第二段代码是错误的：未分配的数组无法进行索引，它没有形状或条目。因此会出现段错误。

（在这种简单情况下，编译器可以很容易地在编译时检测到此错误并给出错误。但是可分配的数组可以在子程序之间传递并在那里进行分配/取消分配，这使得在一般情况下很难检测到该错误。这就是为什么不进行此检查的原因，我想。）

第三段代码：子例程/函数参数中的“假定形状”数组a(:,:)采用传递给例程的数组的形状和值。它通常是对同一内存的引用。

第四和第五段代码：Fortran数组是按列主序排列的，这意味着最左边的索引应该是最内层的循环等等，以实现最快的执行速度。像这样：

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可能只是编译器的等效语句。

- Jonatan Öström