Fortran equivalent to matlab find - application to slicing matrix without memory duplication

Question

I use the command find quite a lot in matlab, and I wonder how to translate this smartly in fortran to extract a slice of an array. In matlab you can slice with either logicals or indexes, but in fortran you need indexes to slice. I'm aware of the intrinsic subroutines pack et al, but have never used them. Also, since I'm dealing with big matrices, I would like to avoid duplicating memory. I want the sliced matrix to be manipulated within a subroutine. I've read somewhere that slices of array were not duplicated. I don't know how this the slicing in done in matlab though. I'm puzzled also because in matlab some allocations are transparent to you.

I'd like to know how to reproduce the examples below, and make sure I'm not duplicating stuff in memory and that it's actually elegant to do so. Otherwise, I would forget about slicing and just send the whole matrix(since it's by reference) and loop through an index array I...

Matlab example 1: simply reproducing find

  v=[1 2 3 4];
  I=find(v==3);

Matlab example 2:

m=rand(4,4);
bools=logical([ 1 0 0 1]);
I=find(bools==1); 
% which I could also do like: 
I=1:size(m,1); 
I=I(bools);

  for i=1:length(I)
      % here dealing with m(I(i)),:)  and do some computation
      % etc.

Example 3: just call a subroutine on m(I,:) , but using directly booleans for slicing

   foo( m(bools, :) , arg2, arg3 )

In advance thank you for your help!

Solution 1

Fortran doesn't have an exact match for Matlab's find but you can generally use either where or forall, and sometimes both, to replace its functionality.

For example, given an array v such as you have in your first example, the Fortran statement

where (v==3) do_stuff

will operate only on the elements of v which are equal to 3. This doesn't get you the indices of those elements as find does, but much of the use of find is for selecting elements for having stuff done to them, and in most of those cases the where construct is applicable.

Given v as before, and an index array ix which, in Fortran, is an array of logicals like this:

[.true., .false., .false., .true.]

you can use ix, so long as it is the same shape as v, in masked array assignments such as

where (ix) v = some_value

ix doesn't have to be an array of logicals, it can be an array of any type, if it were an array of reals you might have an expression such as

where (ix>=0.0) v = some_value

I don't think that any of the current Fortran compilers make copies of arrays to implement where constructs. I'll leave you to read about the forall construct.

Don't forget, either, that you can use arrays as indices for Fortran arrays, so the expression

v([1,3]) = 0

sets elements 1 and 3 of v to 0. You can, of course, use multiple array indices if your array has rank greater than 1.

When you start using this sort of indexing to pass non-contiguous sections of an array to a sub-program, then you have to start worrying about copying into temporary arrays (if that's the sort of thing that you want to worry about). I believe that compilers may make temporary copies if you do something like

call my_subroutine(array(1:12:3, 2:12:4))

to enable the subroutine, which does not know the indices of the elements of the array section at run-time, to operate on what it 'sees' as a contiguous array.

Solution 2

You can use "pack" with an implied do loop:

I = pack([(j,j=1,size(v))],v==3)

Solution 3

Bellow, in FORTRAN CODE, is an example of a subroutine that makes the equivalent of find in matlab or scilab. In the examples below, we want to (a) find the position of the vector where the vector is equal to 22 (b) find the positions of even numbers in a vector

     PROGRAM Principal
     REAL*8 A(8)
     INTEGER n, npos, pos(8)
     n=8
     A = (/ 19, 20, 21, 22, 23, 24, 25, 26 /)
     ! Find the positions of vector A that is equal to 22 
     CALL FindInVector(n,A==22,npos,pos)
     WRITE(*,*) pos(1:npos)

     ! Find the positions of vector A that contains even numbers 
     CALL FindInVector(n,ABS(A/2.d0-INT(A/2.d0))<1.d-2,npos,pos)
     WRITE(*,*) pos(1:npos)

     END PROGRAM Principal

!________________________________________________________________

! Subroutine that find positions of a vector of logicals TF (True or False). The first npos elements contains the response.

     SUBROUTINE FindInVector(n,TF,npos,pos)
    ! Inlet variables
    INTEGER,INTENT(IN):: n      ! Dimension of logical vector
    LOGICAL,INTENT(IN):: TF(n)  ! Logical vector (True or False)
    ! Outlet variables
    INTEGER npos                ! number of "true" conditions
    INTEGER pos(n)              ! position of "true" conditions
    ! Internal variables
    INTEGER i                   ! counter
    INTEGER v(n)                ! vector of all positions

    pos = 0                     ! Initialize pos
    FORALL(i=1:n)   v(i) = i    ! Enumerate all positions
    npos  = COUNT(TF)           ! Count the elements of TF that are .True.
    pos(1:npos)= pack(v, TF)    ! With Pack function, verify position of true conditions

    ENDSUBROUTINE FindInVector