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?
A masked array is the combination of a standard numpy. ndarray and a mask. A mask is either nomask , indicating that no value of the associated array is invalid, or an array of booleans that determines for each element of the associated array whether the value is valid or not.
C arrays always start at zero, but by default Fortran arrays start at 1. There are two usual ways of approaching indexing. You can use the Fortran default, as in the preceding example. Then the Fortran element B(2) is equivalent to the C element b[1] .
Fortran stores higher dimensional arrays as a contiguous sequence of elements. It is important to know that 2-dimensional arrays are stored by column. So in the above example, array element (1,2) will follow element (3,1). Then follows the rest of the second column, thereafter the third column, and so on.
All arrays consist of contiguous memory locations. The lowest address corresponds to the first element and the highest address to the last element. Arrays can be one- dimensional (like vectors), two-dimensional (like matrices) and Fortran allows you to create up to 7-dimensional arrays.
Use the MERGE intrinsic function:
masked = my_array * merge(1,0,my_array<=15.0)
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.
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_array
s 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).
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With