Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Fortran nested WHERE statement

I have a Fortran 90 source code with a nested WHERE statement. There is a problem but it seems difficult to understand what exactly happens. I would like to transform it into DO-IF structure in order to debug. What it is not clear to me is how to translate the nested WHERE.

All the arrays have the same size.

WHERE (arrayA(:) > 0)
    diff_frac(:) = 1.5 * arrayA(:)
    WHERE (diff_frac(:) > 2)
        arrayC(:) = arrayC(:) + diff_frac(:)
    ENDWHERE
ENDWHERE

My option A:

DO i=1, SIZE(arrayA)
  IF (arrayA(i) > 0) THEN
    diff_frac(i) = 1.5 * arrayA(i)
    DO j=1, SIZE(diff_frac)
      IF (diff_frac(j) > 2) THEN
          arrayC(j) = arrayC(j) + diff_frac(j)
      ENDIF
    ENDDO
  ENDIF
ENDDO

My option B:

DO i=1, SIZE(arrayA)
  IF (arrayA(i) > 0) THEN
    diff_frac(i) = 1.5 * arrayA(i)        
    IF (diff_frac(i) > 2) THEN
        arrayC(i) = arrayC(i) + diff_frac(i)
    ENDIF        
  ENDIF
ENDDO

Thank you

like image 632
rgrun Avatar asked Mar 04 '16 12:03

rgrun


People also ask

What is nested IF statement in Fortran?

Fortran - Nested If Construct You can use one if or else if statement inside another if or else if statement(s).

What is the function of the Fortran where statement?

The where statement was introduced in FORTRAN 90 to aid in operations involving arrays. It provides a way to mask the assignment of arrays or the evaluations of arrays.


2 Answers

According to the thread "Nested WHERE constructs" in comp.lang.fortran (particularly Ian's reply), it seems that the first code in the Question translates to the following:

do i = 1, size( arrayA )
    if ( arrayA( i ) > 0 ) then
        diff_frac( i ) = 1.5 * arrayA( i )
    endif
enddo

do i = 1, size( arrayA )
    if ( arrayA( i ) > 0 ) then
        if ( diff_frac( i ) > 2 ) then
            arrayC( i ) = arrayC( i ) + diff_frac( i )
        endif
    endif
enddo

This is almost the same as that in Mark's answer except for the second mask part (see below). Key excerpts from the F2008 documents are something like this:

7.2.3 Masked array assignment – WHERE (page 161)

7.2.3.2 Interpretation of masked array assignments (page 162)

... 2. Each statement in a WHERE construct is executed in sequence.

... 4. The mask-expr is evaluated at most once.

... 8. Upon execution of a WHERE statement that is part of a where-body-construct, the control mask is established to have the value m_c .AND. mask-expr.

... 10. If an elemental operation or function reference occurs in the expr or variable of a where-assignment-stmt or in a mask-expr, and is not within the argument list of a nonelemental function reference, the operation is performed or the function is evaluated only for the elements corresponding to true values of the control mask.

If I understand the above thread/documents correctly, the conditional diff_frac( i ) > 2 is evaluated after arrayA( i ) > 0, so corresponding to double IF blocks (if I assume that A .and. B in Fortran does not specify the order of evaluation).


However, as noted in the above thread, the actual behavior may depend on compilers... For example, if we compile the following code with gfortran5.2, ifort14.0, or Oracle fortran 12.4 (with no options)

integer, dimension(4) :: x, y, z
integer :: i

x = [1,2,3,4]
y = 0 ; z = 0

where ( 2 <= x )
   y = x
   where ( 3.0 / y < 1.001 )  !! possible division by zero
      z = -10
   end where
end where

print *, "x = ", x
print *, "y = ", y
print *, "z = ", z

they all give the expected result:

x =            1           2           3           4
y =            0           2           3           4
z =            0           0         -10         -10

But if we compile with debugging options

gfortran -ffpe-trap=zero
ifort -fpe0
f95 -ftrap=division  (or with -fnonstd)

gfortran and ifort abort with floating-point exception by evaluating y(i) = 0 in the mask expression, while f95 runs with no complaints. (According to the linked thread, Cray behaves similarly to gfortran/ifort, while NAG/PGI/XLF are similar to f95.)


As a side note, when we use "nonelemental" functions in WHERE constructs, the control mask does not apply and all the elements are used in the function evaluation (according to Sec. 7.2.3.2, sentence 9 of the draft above). For example, the following code

integer, dimension(4) :: a, b, c

a = [ 1, 2, 3, 4 ]
b = -1 ; c = -1

where ( 3 <= a )
    b = a * 100
    c = sum( b )
endwhere

gives

a =  1 2 3 4
b =  -1 -1 300 400
c =  -1 -1 698 698

which means that sum( b ) = 698 is obtained from all the elements of b, with the two statements evaluated in sequence.

like image 177
roygvib Avatar answered Jan 02 '23 12:01

roygvib


Why not

WHERE (arrayA(:) > 0)
    diff_frac(:) = 1.5 * arrayA(:)
ENDWHERE

WHERE (diff_frac(:) > 2 .and. arrayA(:) > 0)
    arrayC(:) = arrayC(:) + diff_frac(:)
ENDWHERE

?

I won't say it can't be done with nested wheres, but I don't see why it has to be. Then, if you must translate to do loops, the translation is very straightforward.

Your own attempts suggest you think of where as a kind of looping construct, I think it's better to think of it as a masked assignment (which is how it's explained in the language standard) in which each individual assignment happens at the same time. These days you might consider translating into do concurrent constructs.

like image 29
High Performance Mark Avatar answered Jan 02 '23 12:01

High Performance Mark