7

I am wanting to mask a Fortran array. Here's the way I am currently doing it...

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

So then I get my masked array with:

masked = my_array * mask_array

Is there a more concise way to do this?

user14241
  • 727
  • 1
  • 8
  • 27

4 Answers4

12

Use the MERGE intrinsic function:

masked = my_array * merge(1,0,my_array<=15.0)
Fortranner
  • 2,525
  • 2
  • 22
  • 25
7

Or, sticking with where,

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

I expect that there are differences, in speed and memory consumption, between the use of where and the use of merge but off the top of my head I don't know what they are.

High Performance Mark
  • 77,191
  • 7
  • 105
  • 161
5

There are two different approaches already given here: one retaining where and one using merge. In the first, High Performance Mark mentions that there may be differences in speed and memory use (think about temporary arrays). I'll point out another potential consideration (without making a value judgment).

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

That is, the merge solution can use automatic allocation. Of course, there are times when one doesn't want this (such as when working with lots of my_arrays of the same size: there are often overheads when checking array sizes in these cases): use masked(:) = MERGE(...) after handling the allocation (which may be relevant even for the question code).

Community
  • 1
  • 1
francescalus
  • 30,576
  • 16
  • 61
  • 96
0

I find it useful to define a function where which takes an array of logicals and returns the integer indices of the .true. values, so e.g.

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

This function can be defined as

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

With this where function, your problem can be solved as

my_array(where(my_array>15.0)) = 0

This is probably not the most performant way of doing this, but I think it is very readable and concise. This where function can also be more flexible than the where intrinsic, as it can be used e.g. for specific dimensions of multi-dimensional arrays.

Limitations:

Note however that (as @francescalus points out) this will not work for arrays which are not 1-indexed. This limitation cannot easily be avoided, as performing comparison operations on such arrays drops the indexing information, e.g.

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".

For arrays which are not 1-indexed, in order to use this where function you would need the rather ugly

my_array(where(my_array>15.0)+lbound(my_array)-1) = 0
veryreverie
  • 2,871
  • 2
  • 13
  • 26
  • Afraid I'm very much against this: consider what happens when trying to use with the array `dimension my_array(-5:12)`. Also, WHERE constructs/statements are _much_ more flexible in general (I agree, for simple cases like of the question this won't matter). – francescalus Jan 06 '22 at 11:08
  • @francescalus true. Personally, I don't use arrays which aren't 1-indexed, because I think it's too easy to drop the indexing. Also, could you give an example (aside from non-1-indexed arrays) where `where` constructs are more flexible than this `where` function? – veryreverie Jan 06 '22 at 11:43
  • (My) "flexible" was probably the wrong word, but `where (f(x)) x=g(y)` would seem to be clearer than the obvious rewritings with a `where` function. As a more general rule, I'd say that if you want to do masked assignment, use masked assignment rather than assignment with a vector subscript or assignment in a loop. (Which relates to a further concern I have: one uses this `where` function, finds it pretty neat, then one day decide to use `x(where(...))` in a place outside assignment where a vector subscript really isn't allowed/wise.) – francescalus Jan 06 '22 at 12:18