屏蔽 Fortran 数组的更好方法？答案

【问题标题】：Better way to mask a Fortran array?屏蔽 Fortran 数组的更好方法？
【发布时间】：2015-04-22 14:40:43
【问题描述】：

我想屏蔽一个 Fortran 数组。这是我目前的做法...

where (my_array <=15.0)
    mask_array = 1
elsewhere
    mask_array = 0
end where

然后我得到我的掩码数组：

masked = my_array * mask_array

有没有更简洁的方法来做到这一点？

【问题讨论】：

标签： arrays fortran fortran90 masking

【解决方案1】：

我发现定义一个函数 where 很有用，它接受 logicals 的数组并返回 .true. 值的 integer 索引，例如

x = where([.true., .false., .false., .true.]) ! sets `x` to [1, 4].

这个函数可以定义为

function where(input) result(output)
  logical, intent(in) :: input(:)
  integer, allocatable :: output(:)
  
  integer :: i
  
  output = pack([(i, i=1, size(input))], input)
end function

有了这个where函数，你的问题就可以解决了

my_array(where(my_array>15.0)) = 0

这可能不是最高效的方式，但我认为它非常易读和简洁。这个where 函数也可以比where 内在函数更灵活，因为它可以用于例如用于多维数组的特定维度。

限制：

但是请注意（正如@francescalus 指出的那样）这不适用于非 1 索引的数组。这种限制无法轻易避免，因为对此类数组执行比较操作会丢弃索引信息，例如

real :: my_array(-2,2)
integer, allocatable :: indices(:)

my_array(-2:2) = [1,2,3,4,5]
indices = my_array>3
write(*,*) lbound(indices), ubound(indices) ! Prints "1 5".

对于不是 1-indexed 的数组，为了使用这个 where 函数，你需要相当丑陋的

my_array(where(my_array>15.0)+lbound(my_array)-1) = 0

【讨论】：

恐怕我非常反对这一点：考虑一下尝试使用数组 dimension my_array(-5:12) 时会发生什么。此外，WHERE 构造/语句通常更加更灵活（我同意，对于像问题这样的简单情况，这无关紧要）。
@francescalus 是的。就个人而言，我不使用不是 1 索引的数组，因为我认为删除索引太容易了。另外，你能举一个例子（除了非 1 索引数组）where 构造比这个where 函数更灵活吗？
(My)“灵活”可能是错误的词，但where (f(x)) x=g(y) 似乎比使用where 函数的明显重写更清晰。作为更一般的规则，我会说，如果您想做屏蔽赋值，请使用屏蔽赋值，而不是使用向量下标进行赋值或循环赋值。（这与我的另一个担忧有关：一个人使用了这个 where 函数，发现它非常整洁，然后有一天决定在赋值之外的地方使用 x(where(...)) 确实不允许/明智地使用向量下标。）

【解决方案2】：

这里已经给出了两种不同的方法：一种保留where，另一种使用merge。首先，高性能标记提到速度和内存使用可能存在差异（考虑临时数组）。我会指出另一个潜在的考虑因素（不做价值判断）。

subroutine work_with_masked_where(my_array)
  real, intent(in) :: my_array(:)
  real, allocatable :: masked(:)

  allocate(masked(SIZE(my_array)), source=0.)
  where (my_array <=15.0) masked = my_array
  ! ...
 end subroutine

subroutine work_with_masked_merge(my_array)
  real, intent(in) :: my_array(:)
  real, allocatable :: masked(:)

  masked = MERGE(my_array, 0., my_array<=15.)
  ! ...
end subroutine

即merge解决方案可以使用自动分配。当然，有时人们不希望这样做（例如，在处理大量相同大小的my_arrays 时：在这些情况下检查数组大小时通常会产生开销）：在处理完之后使用masked(:) = MERGE(...)分配（甚至可能与问题代码相关）。

【讨论】：

【解决方案3】：

或者，坚持where，

masked = 0
where (my_array <=15.0) masked = my_array

我预计where 和merge 的使用在速度和内存消耗方面存在差异，但我不知道它们是什么。

【讨论】：

【解决方案4】：

使用 MERGE 内在函数：

masked = my_array * merge(1,0,my_array<=15.0)

【讨论】：

完美！这正是我想要的。
同一个主题masked = MERGE(my_array, 0._blah, my_array<=15).