# Indexing arrays by an array of logicals

Imagine that you have two arrays `A` and `B` of the same shape, and you would like to set `A` to `0` wherever `B` is positive. This can be done by

``````where (B > 0)
A = 0
end where
``````

This is great, but in MATLAB, it can be done in a more concise way, which is

``````A(B > 0) = 0
``````

Here, A is indexed by an array of logicals (masks). Operations are affecting only the positions where the mask is `.true.`. This seems to me a very convenient feature. It might be something worth considering in Fortran as well.

See here for more information about logical indexing in MATLAB: Matrix Indexing in MATLAB - MATLAB & Simulink

Python has a similar feature:
https://www.pythonlikeyoumeanit.com/Module3_IntroducingNumpy/AdvancedIndexing.html

4 Likes

There are a few ways more ways to achieve the same:

``````do concurrent (i = 1:size(b), b(i) > 0)
a(i) = 0
end do
``````
``````a(findloc(b > 0,.true.)) = 0
``````

I can imagine a small helper function would make the second version shorter, but perhaps less performant. I guess it doesnâ€™t generalize to arbitrary rank either.

5 Likes

If your only purpose is to make the code look nicer and you donâ€™t need `elsewhere`, nothing prevent you from using

``````where (b > 0) a = 0
``````

This is standard approval (F2008 7.2.3.1), and to me this is closer to natural language than `a(b > 0) = 0`.

8 Likes

In Fortran, array indices have to be integers; therefore, â€śan array of logicalsâ€ť cannot be used to index an array (of any type).

In the proposed assignment

`A(B > 0) = 0`

what is the type of the expression B > 0 ? Is it an array of integers or is it an array of logical values?

The following short program illustrates the issues. Many ideas may seem â€śsmartâ€ť in a narrow context or to offer convenience, but (i) will they work in all contexts allowed by the language, and (ii) will the suggested changes in semantics break existing code?

``````program WhereWhat
implicit none
integer, parameter :: N = 5
integer :: ix(N), A(N), B(N)
logical :: gt(N)
!
B = [2, 3, 4, 5, 7]
A = B
where (B > 4) A = -A
print 10,'A : ',A
print *,'(B > 4) :      ',(B > 4) ! what is the type of the expression?
!ix = B > 4     ! compile time error, logical array value placed into integer array
gt = B > 4
print 10,'ix = ',ix
print 20,'gt = ',gt
print 10,'A(ix) = ',A(ix)  ! illegal subscripts used, results meaningless
!print *,'A(gt) = ',A(gt) ! compile time error: logical values used as subscripts
10 format(1x,A,T15,5I4)
20 format(1x,A,T15,5L4)
end program

``````
1 Like

In MATLAB this would be an array of logical values (represented as zeros and ones).

MATLAB also has a function `find` which can be used to map logical values (or nonzero elements more generally) to linear indexes. The syntax is:

``````k = find(X)         % or
[row,col] = find(X)
``````

This could be emulated to some extent in Fortran, but would be more verbose than the `where` construct, or the use of a logical array as an index, like in MATLAB.

2 Likes

Yes, the value of the logical expression B > 0 as represented in Matlab is [1, 1, 1, 0, 0], so we cannot use that expression as an index expression. The 0 indices are less than the lower bound of the array, and the only other value is 1, and the result is not at all what we want. That is the problem: we need to implicitly map the positions that have 1 as the value into a new integer valued index array, [1, 2, 3], and use those values as the indices into A. Furthermore, this would have to be done consistently in all contexts, since in Fortran the value of an expression has an associated type.

1 Like

I have written a function that returns a vector index from a logical array

``````function true_pos(tf) result(ipos)
! return in ipos the positions of the true elements in tf(:)
logical, intent(in) :: tf(:)
integer             :: ipos(count(tf))
integer             :: i,j
j = 0
do i=1,size(tf)
if (tf(i)) then
j = j + 1
ipos(j) = i
end if
end do
end function true_pos
``````

and would write

`A(true_pos(B > 0)) = 0`

Logical indexing would be convenient, but given a simple alternative Iâ€™m not sure it is worth changing the language. If you are referring to array elements but not setting them you can use the expression

`pack(A,B>0)`

3 Likes

@zaikunzhang thanks for bringing this up. Indeed, Python/NumPy also allows this.

There are always usually some reasons why it wasnâ€™t done in Fortran, as others have pointed out above. It might also be very difficult to make such changes now and keep everything backwards compatible.

However, I still like to have these discussions and to brainstorm such possibilities. If nothing else, one could for example use these to create better compiler error messages, say if you write `A(B>0) = 0`, the compiler can tell you to write it as `where (B>0) A = 0`. But I think these discussions can also be useful either to try to improve Fortran, or if somebody is designing a new language, how it could be done.

3 Likes

@Beliavsky, I think this would make a nice addition to the standard library. Could you open a new issue for it?

Some care is required to make the proposed function work correctly when its argument is not an array with a lower bound of 1, as the following test code demonstrates (it uses Beliavskyâ€™s true_pos function from #7).

``````program xtpos
implicit none
integer :: ia(2:6) = [1,2,3,4,5], i
!
write(*,'(5(1x,I3))')(i,i=lbound(ia,1),ubound(ia,1))
do i=lbound(ia,1),ubound(ia,1)
write(*,'(1x,L3)',advance='no')ia(i) > 2
end do
print *
print 10, 'True_Pos',true_pos(ia > 2)
print 10, 'Expected',[4,5,6]
10 format(A8,2x,3(1x,I3))
contains

function true_pos(tf) result(ipos)
! return in ipos the positions of the true elements in tf(:)
logical, intent(in) :: tf(:)
integer             :: ipos(count(tf))
integer             :: i,j
j = 0
do i=1,size(tf)
if (tf(i)) then
j = j + 1
ipos(j) = i
end if
end do
end function true_pos

end program
``````

The output:

``````   2   3   4   5   6
F   F   T   T   T
True_Pos     3   4   5
Expected     4   5   6
``````

A related issue: what should the function return when an array section or an array reference with a vector subscript is the input argument?

2 Likes

What? None of you suggested:

`a = merge (a,0,b>0)`

???

6 Likes

Nice! I donâ€™t use merge in my codes, so I havenâ€™t thought of that.

When you spend too much time merging pull requests, you forget about the Fortran `merge`.

1 Like

I generally would not recommend `MERGE` for any instructions except for those involving scalar arguments. My experience is both the practitioners and compiler implementations donâ€™t quite â€śgetâ€ť this intrinsic.

Exhibit A: going by the original post, it should be

``````a = merge (0, a, b>0)
``````

Exhibit B: compiler implementations can do needless operations.

C:\Temp>a.exe
Block 1: use MERGE
In assign_w_t: performing lhs%n(currently= 1 ) = 1 assignment
In assign_w_t: performing lhs%n(currently= 2 ) = 0 assignment
In assign_w_t: performing lhs%n(currently= 3 ) = 3 assignment
In assign_w_t: performing lhs%n(currently= 4 ) = 4 assignment
In assign_w_t: performing lhs%n(currently= 5 ) = 5 assignment
In assign_w_t: performing lhs%n(currently= 6 ) = 6 assignment
In assign_w_t: performing lhs%n(currently= 7 ) = 7 assignment
In assign_w_t: performing lhs%n(currently= 8 ) = 0 assignment
In assign_w_t: performing lhs%n(currently= 9 ) = 9 assignment
In assign_w_t: performing lhs%n(currently= 10 ) = 10 assignment
a = 1 0 3 4 5 6
7 0 9 10

Block 2: use WHERE
In assign_w_t: performing lhs%n(currently= 2 ) = 0 assignment
In assign_w_t: performing lhs%n(currently= 8 ) = 0 assignment
a = 1 0 3 4 5 6
7 0 9 10

C:\Temp>

Click for code!

## Illustrate array operations using a derived type

``````module w_m
type :: w_t
integer :: n = 0
end type
generic :: operator(>) => compare_w_t
generic :: assignment(=) => assign_w_t
type(w_t), parameter :: ZERO = w_t()
contains
elemental function compare_w_t( lhs, rhs ) result(r)
type(w_t), intent(in) :: lhs
type(w_t), intent(in) :: rhs
logical :: r
r = lhs%n > rhs%n
end function
impure elemental subroutine assign_w_t( lhs, rhs )
type(w_t), intent(inout) :: lhs
type(w_t), intent(in)    :: rhs
print *, "In assign_w_t: performing lhs%n(currently=", lhs%n, ") = ", rhs%n, " assignment"
lhs%n = rhs%n
end subroutine
end module
use w_m
type(w_t), dimension(10) :: a, b
blk1: block
print *, "Block 1: use MERGE"
call init()
a = merge( ZERO, a, b > ZERO )  !<-- use of MERGE intrinsic
print *, "a = ", a%n
end block blk1
print *
blk2: block
print *, "Block 2: use WHERE"
call init()
where ( b > ZERO ) a = ZERO   !<-- WHERE statement
print *, "a = ", a%n
end block blk2
contains
subroutine init()
b(2)%n = 1 ; b(8)%n = 1
do concurrent ( integer :: i=1:size(a) )
a(i)%n = i
end do
end subroutine
end
``````

C:\Temp>ifort /standard-semantics a.f90
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.3.0 Build 20210609_000000
Copyright (C) 1985-2021 Intel Corporation. All rights reserved.

Microsoft (R) Incremental Linker Version 14.29.30038.1
Copyright (C) Microsoft Corporation. All rights reserved.

-out:a.exe
-subsystem:console
a.obj

C:\Temp>

3 Likes

It is true that MERGE is not well optimized, but it IS the language feature being asked for and if more people used it, compiler writers would work harder to optimize it. I think Cray Fortran already does this to some degree.

One could have said similar things about array assignments in the early 90s.

2 Likes

No, in the context of this thread, what is being asked for is some sort of â€śprecis writingâ€ť for what the language offers with `where ( b > 0 ) a = 0` and that ainâ€™t `MERGE`.

2 Likes

To be fair, `merge` is one way it could be done (and letâ€™s assume compilers implement this well in the future), but as beautifully shown by you above, the `merge` syntax is so confusing that anybody can get it wrong. I personally do not use it for the same reason, I just never found it clear enough what it does to warrant using it. `where` is pretty good. So is the NumPy syntax `a[b>0] = 0`.

1 Like

Thank everyone for the ideas and insights, from which I have learned a lot, particularly the usage of `merge` proposed by Dr. Fortran @sblionel.

It seems that the idea of logical indexing has gained little popularity during the discussion. There do exist a few languages for numerical computing, namely Python (as mentioned by @certik), Julia, R, and MATLAB, that support logical indexing. So I do not think it is just yet another â€śsmart ideaâ€ť.

However, I do understand that including logical indexing necessities significant augmentation to the existing indexing mechanism of Fortran, and I do not assume by any means that such augmentation is trivial to implement.

Meanwhile, it will be very interesting to see a snippet illustrating that logical indexing can hurt backward compatibility, namely a piece of code that works now but would become broken if logical indexing were introduced. Mathematically, it is a bit hard for me to imagine how an augmentation/generalization (rather than a modification) can harm backward compatibility.

It may be worth mentioning that `A(B>0) = 0` is only a simple illustration for logical indexing. I totally agree that `where` as suggested by @han190 is a better implementation for this job, because `where` is closer to the natural language. However, logical indexing can be applied in much broader contexts, for instance, `A(B>0) = C(D<0)`. The `true_pos` proposed by @Beliavsky is very nice, and it is essentially the same as logical indexing except that it is implemented as a function.

Anyway, even if logical indexing were to be included in the language by a future standard, I would not expect my projects to make use of it. I am quite conservative regarding new features in order to keep my code as portable as possible. I am still refraining from some convenient F2008 features because they are not supported yet by some recent compilers; for example, Absoft 21.0 does not support automatic allocation upon intrinsic assignments. One may complain that the compiler vendors are too slow to support the standards, but one can never blame users for using such compilers.

2 Likes

I do not find merge at all confusing but then I use it with high frequency. I find many things confusing and needing a visit to help for the sole reason of infrequent use. Many simple if /else/ endif blocks become one line instead five lines! It saves typing and makes the code easier to scan read IMO.

3 Likes

In MoA,
(
(((rav A) = (rav B))/(iota (x red ( rho A)))) Omega (<0 1>) <2 4>
) psi A = 0

So when psi reduced, we just do the assignment.
The rav A = rav B Creates a Boolean vector, when compressed with iota product rho the shape, finds the index, then using Omega, applies gamma prime to get the indices in the two by 4 array , eg. Now we assign 0 to them.