Logical Indexing based on "find" in Fortran 90

Question

I am trying to make a logical array (B) to use in logical indexing based on values between .1 and .999 in an array (EP_G2) using a couple different methods 1) where loop 2) ANY.

program flux_3d

implicit none 
INTEGER :: RMAX, YMAX, ZMAZ, timesteps
DOUBLE PRECISION, PARAMETER ::  pmin=0.1
DOUBLE PRECISION, PARAMETER ::  pmax=0.999
INTEGER :: sz
DOUBLE PRECISION, ALLOCATABLE :: EP_G2(:,:), C(:)
INTEGER, DIMENSION(RMAX*ZMAX*YMAX) :: B
LOGICAL, DIMENSION(RMAX*ZMAX*YMAX) :: A
! dimensions of array, 
RMAX = 540 
YMAX = 204 
ZMAX = 54 
timesteps = 1
!Open ascii array
OPEN(100, FILE ='EP_G2', form = 'formatted')
ALLOCATE( EP_G2(RMAX*ZMAX*YMAX, timesteps))
READ(100, *) EP_G2

WHERE(pmin<EP_G2(:,timesteps)<pmax) B = 1 
    ELSEWHERE
        B = 0
ENDWHERE

PRINT*, B

! SZ # OF POINTS IN B
sz = count(B.eq.1)

!alternate way of finding points between pmin and pmax
A = ANY(pmin<EP_G2(:,t)<pmax)
print*, sz

!Then use the logical matrix B (or A) to make new array 
!Basically I want a new array that isolates just the points in array between 
!pmin and pmax
ALLOCATE(C(sz))
C = EP_G2(LOGICAL(B), 1)

The issue I am getting is that I either get the whole array or nothing and the command LOGICAL(B) gets an error that it isn't the right kind. I am new to Fortran coming from Matlab where I would just use find. Additionally, since this array is over 5,948,640 x 1 computational time is an issue. I am using the Intel Fortran compiler (15.0 I believe).

Basically, I am looking for the fastest way to find the indexes of points in an array between two numbers.

Any help would be very appreciated.

Answer 1

I'm a little confused by your code and question but I think you want to find the indices of elements in an array whose values lie between 0.1 and 0.999. It's not entirely clear that you want the indices, or just the elements themselves, but bear with me and I'll explain how to get both.

Suppose your original array is declared something like

real, dimension(10**6) :: values

and that it is populated with values. Then I might declare an index array like this

integer, dimension(10**6) :: indices = [(ix,ix=1,10**6)]  

(obviously you'll need to declare the variable ix too).

Now, the expression

pack(indices, values>pmin.and.values<pmax)

will return those elements of indices whose values lie strictly between pmin and pmax. Of course, if you only want the values themselves you can dispense with indices entirely and write

pack(values, values>pmin.and.values<pmax)

Both of these uses of pack will return an array of rank 1, and you can assign the returned array to an allocatable, Fortran will take care of the sizing for you.

I'll leave it to you to expand this approach to the rank-2 array you are actually working with, but it shouldn't be too difficult to do that. Ask for more help if you need it.

And is this the fastest approach ? I'm not sure, but it was very quick to write and if you're interested in its run-time performance I suggest you test that.

Incidentally, the Fortran intrinsic function logical is for transforming logical values from one kind to another, it's not defined on integers (or any other intrinsic type). Fortran predates the madness of regarding 0 and 1 as logical values.

Answer 2

How can the arrays be the same?

RMAX, YMAX, ZMAZ, and timesteps values are not known until after A, B, and C are declared - So they (A and B) will not likely be the size you want.

implicit none 
INTEGER :: RMAX, YMAX, ZMAZ, timesteps
DOUBLE PRECISION, PARAMETER ::  pmin=0.1
DOUBLE PRECISION, PARAMETER ::  pmax=0.999
INTEGER :: sz
DOUBLE PRECISION, ALLOCATABLE :: EP_G2(:,:), C(:)
INTEGER, DIMENSION(RMAX*ZMAX*YMAX) :: B 
LOGICAL, DIMENSION(RMAX*ZMAX*YMAX) :: A
! dimensions of array, 
RMAX = 540 
YMAX = 204 
ZMAX = 54 
timesteps = 1

You probably want either this:

implicit none 
INTEGER ,         PARAMETER :: RMAX      = 504
INTEGER ,         PARAMETER :: YMAX      = 204
INTEGER ,         PARAMETER :: ZMAZ      = 54
INTEGER ,         PARAMETER :: timesteps = 1
DOUBLE PRECISION, PARAMETER :: pmin      = 0.1
DOUBLE PRECISION, PARAMETER :: pmax      = 0.999
INTEGER :: sz
DOUBLE PRECISION, ALLOCATABLE :: EP_G2(:,:), C(:)
INTEGER, DIMENSION(RMAX*ZMAX*YMAX) :: B 
LOGICAL, DIMENSION(RMAX*ZMAX*YMAX) :: A
! dimensions of array, 
!RMAX = 540 
!YMAX = 204 
!ZMAX = 54 
!timesteps = 1

Or this:

implicit none 
! dimensions of array
INTEGER ,         PARAMETER :: RMAX      = 504
INTEGER ,         PARAMETER :: YMAX      = 204
INTEGER ,         PARAMETER :: ZMAZ      = 54
INTEGER ,         PARAMETER :: timesteps = 1

DOUBLE PRECISION, PARAMETER :: pmin      = 0.1
DOUBLE PRECISION, PARAMETER :: pmax      = 0.999
INTEGER :: sz
DOUBLE PRECISION, ALLOCATABLE :: EP_G2(:,:), C(:)
INTEGER, DIMENSION(:), ALLOCATABLE  :: B 
LOGICAL, DIMENSION(:), ALLOCATABLE  :: A
! Then allocate A and B

And you may also want to consider the use of SHAPE or SIZE to see if the rank and size of the array are correct.

IF(SHAPE(A) /= SHAPE(B) )  ... chuck and error message.
IF(SIZE(A,1) /= SIZE(B,1) ) etc