Fortran 流程控制（二）：forall和do concurrent孰优孰劣

自学气象人

发布于 2023-06-20 16:43:42

4870

发布于 2023-06-20 16:43:42

文章被收录于专栏：自学气象人自学气象人

在《Fortran 流程控制（一）：where》一文中，我们介绍了一种面向数组的条件判断结构，类似于面向标量的if结构。对于数组，同样有类似于标量里的do循环类似的结构：forall与do concurrent。

FORALL 结构

forall结构可以看作是隐式循环的一种拓展，可以实现通过条件判断是否给数组赋值的功能。

基本格式

forall (triplet1[, triplet2[, triplet3...]], scalar_mask_expression)
    statement1 
    statement2
    ...
end forall

triplet是用于赋值的数组坐标范围，scalar_mask_expression为条件判断值，只有scalar_mask_expression成立才会运行forall中的语句。

Example

program forall_construct
  implicit none
  integer, parameter :: N = 5
  character(len=100) :: outFormat
  integer :: i, j
  real :: a(N,N) = 0, b(N,N) = 0, threshold = 0.5, &
       c(N,N) = 0, d(N,N) = 0 ! used in next examples

  ! write some values in a
  call random_number( a )

  ! Create dynamic format, with internal-file(=string) outFormat.
  ! This way, the format is adjusted automatically if N changes.
  write(outFormat, *) "(", N, "(x, f8.2))"

  write(*, '(a)') "a = "
  write(*, fmt=outFormat) &
       ( (a(i,j), j=1,N), i=1,N )

  ! ** Forall construct **
  forall( i = 2:N:1, j = 1:N:1)
     b(i, j) = a(i, j)
  end forall

  forall( i = 1:N, j = 1:N, i==j)
     c(i, j) = a(i, j)
  end forall

  forall( i = 1:N, j = 1:N, a(i, j) > threshold)
     d(i, j) = a(i, j)
  end forall

  write(*, '(/,a)') "b = "
  write(*, fmt=outFormat) ( (b(i,j), j=1,N), i=1,N )

  write(*, '(/,a)') "c = "
  write(*, fmt=outFormat) ( (c(i,j), j=1,N), i=1,N )

  write(*, '(/,a)') "d = "
  write(*, fmt=outFormat) ( (d(i,j), j=1,N), i=1,N )

end program forall_construct

输出：

a = 
     0.00     0.84     0.35     0.73     0.04
     0.03     0.34     0.87     0.30     0.09
     0.35     0.92     0.09     0.05     0.56
     0.67     0.80     0.89     0.91     0.93
     0.96     0.83     0.70     0.10     0.08

b = 
     0.00     0.00     0.00     0.00     0.00
     0.03     0.34     0.87     0.30     0.09
     0.35     0.92     0.09     0.05     0.56
     0.67     0.80     0.89     0.91     0.93
     0.96     0.83     0.70     0.10     0.08

c = 
     0.00     0.00     0.00     0.00     0.00
     0.00     0.34     0.00     0.00     0.00
     0.00     0.00     0.09     0.00     0.00
     0.00     0.00     0.00     0.91     0.00
     0.00     0.00     0.00     0.00     0.08

d = 
     0.00     0.84     0.00     0.73     0.00
     0.00     0.00     0.87     0.00     0.00
     0.00     0.92     0.00     0.00     0.56
     0.67     0.80     0.89     0.91     0.93
     0.96     0.83     0.70     0.00     0.00

do concurrent 结构

do concurrent结构在Fortran 2008中引入，可以用来改善数组操作的性能和简洁性。严格来讲，这种结构更加通用，因为也可以用它来处理标量数据。然而，我们在这里讨论它，因为它对数组特别有用，并且还因为它有效地取代了另一个面向数组的结构，即前文所提到的forall。

基本格式

do concurrent( [type_spec ::] list_of_indices_with_ranges & [, scalar_mask_expression] ) 
    statement1 
    statement2 
    ... 
end do

其中，list_of_indices_with_ranges可以是索引范围规范（如在正常do循环后出现那样），也可以是此类规范的逗号分隔列表（在这种情况下，构造等同于一组嵌套循环）。我们在这部分的结尾讨论了可选的type_spec。scalar_mask_expression（如果存在）用于将语句应用程序限制为索引范围内表达式求值为.true.的部分片段。

Example

以下示例对do concurrent结构进行了说明，其中属于棋盘图案的矩阵a的元素被复制到矩阵b中：

program do_concurrent_checkerboard_selection
  implicit none
  integer, parameter :: DOUBLE_REAL = selected_real_kind(15, 307)
  integer, parameter :: N = 5 ! side-length of the matrices
  integer :: i, j ! dummy-indices
  real(kind=DOUBLE_REAL), dimension(N,N) :: a, b ! the matrices
  character(len=100) :: outFormat

  ! Create dynamic format, using internal file
  write(outFormat, *) "(", N, "(x, f8.2))"
  ! Initialize matrix a to some random values
  call random_number( a )

  ! Pattern-selection with do concurrent
  do concurrent( i=1:N, j=1:N, mod(i+j, 2)==1 )
     b(i,j) = a(i,j)
  end do

  ! Print matrix a
  write(*, '(/,a)') "a ="
  write(*, fmt=outFormat) ( (a(i,j), j=1,N), i=1,N )

  ! Print matrix b
  write(*, '(/,a)') "b ="
  write(*, fmt=outFormat) ( (b(i,j), j=1,N), i=1,N )
end program do_concurrent_checkerboard_selection

输出：

a =
     0.00     0.84     0.35     0.73     0.04
     0.03     0.34     0.87     0.30     0.09
     0.35     0.92     0.09     0.05     0.56
     0.67     0.80     0.89     0.91     0.93
     0.96     0.83     0.70     0.10     0.08

b =
     0.00     0.84     0.00     0.73     0.00
     0.03     0.00     0.87     0.00     0.09
     0.00     0.92     0.00     0.05     0.00
     0.67     0.00     0.89     0.00     0.93
     0.00     0.83     0.00     0.10     0.00

在语法构成上，上述示例中第15-17行的代码可以用do循环和if语句实现其功能，如下所示：

do i=1,N
    do j=1,N
        if( mod(i+j, 2)==1 ) then
            b(i,j) = a(i,j)
        end if
    end do
end do

可以看出，使用do concurrent的代码结构更加紧凑。更重要的是，该结构还允许使用嵌套do循环对版本进行一些编译器优化。

优势与限制

不过，do concurrent也有其使用上的限制。其中一些限制是编译器可以检查的（如果违反了这些限制，则会发出编译时错误），而另一些则无法自动检查，需要程序员自行保证满足这些限制。举个例子：

• 大多数限制与防止程序员在do concurrent结构之外分支有关。造成这种分支的示例包括return、go to 、exit、cycle或者err=（用于错误处理）。安全的处理方法是避免使用这些语句。
• 允许在结构体中调用其他程序，只要这些程序是纯的，这意味着该程序没有副作用。会使程序变得不纯的副作用的例子有：
- • 将全局或局地实体中的程序状态更改为下次调用该过程时可能使用的过程。
- • 在一次迭代期间产生输出，在另一次迭代中读取输出。
• 程序员也要保证在迭代过程中编译器没有数据依赖（比如，通过共享变量，在一次迭代过程中分配数据内存，在另一个迭代过程中释放数据内存，或者在不同的迭代过程中从外部渠道读取和写入数据）。由此可见，虽然do concurrent结构有着更加紧凑的优势，但也存在结构上的使用限制，会使得代码脚本难以更加通用。这就意味着在使用do concurrent结构时需要权衡利弊。然而，对于一些优先考虑性能的应用程序，可以考虑使用do concurrent结构，因为它迫使程序员以有利于后期并行化的方式重新构造算法。 " "type_spec" 选项 关于do concurrent结构，一个有趣的注意事项是：标准还允许指定结构中索引的类型（类型总是integer，但kind参数可以自定义）。这非常方便，因为它让类型定义紧挨着变量使用的地方（否则，这些索引需要在（子）程序的开头声明，如先前的示例一般）。例如，先前示例程序中的模式选择部分可以写成：do concurrent( integer :: l=1:N, m=1:N, mod(l+m, 2) == 1) b(l,m) = a(l,m) end do 参考
• Chirila, D.B., Lohmann, G., 2015. Introduction to Modern Fortran for the Earth System Sciences. Springer Berlin Heidelberg, Berlin, Heidelberg.