FINDLOC equivalent for NAN

ereisch · ‎02-26-2021

Is there a FINDLOC() equivalent that can be used for IEEE NAN values? It is my understanding I cannot use FINDLOC with IEEE_VALUE( 1,0, IEEE_QUIET_NAN ), because the IEEE NAN values (signaling or quiet) don't equal anything, "including itself".

I can obviously just put in a DO loop across the array and check for IEEE_IS_NAN, but I'm assuming (possibly erroneously?) that FINDLOC() has a more efficient implementation under the hood.

Arjen_Markus · ‎02-26-2021

It is probably not what you are looking for, but here is a concise method:

write(*,*) findloc( ieee_is_nan(array) )

ereisch · ‎02-26-2021

While that's concise in terms of written code, it's very expensive, as IEEE_IS_NAN() will need to operate on the entirety of array before returning to FINDLOC(). So if array is very large and you have a NAN value fairly early, it will still perform the check on all the remaining members of array, whereas (theoretically) FINDLOC() will short-circuit.

jimdempseyatthecove · ‎02-26-2021

Presumably a NAN will be infrequent in your code.

IF(IEEE_IS_NAN(SUM(array))) THEN
! look for it
...
ENDIF

.OR.

the SUM(array) can be replaced by something in the prior code that contains the accumulation of what was in the array. e.g. total momentum, total energy, ...

Jim Dempsey

ereisch · ‎02-26-2021

Unfortunately, they are a little more frequent than usual, as this particular array we're initializing to NaN at startup (so as to not confuse with a "perfectly valid" 0.0). The code block in question is to find the first "unused" member of the array.

jimdempseyatthecove · ‎02-26-2021

A computed NAN is not a specific value, whereas a compiler generated NAN is a specific value.

The computed NAN will have

s111 1111 1qxx xxxx xxxx xxxx xxxx xxxx (single precision)
s111 1111 1111 qxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx xxxx (double)

s=sign (don't care)
1's exponent of all 1's
q=1=QNAN, q=0=SNAN
x's remaining bits of any value except all 0's

It would be recommended that your code should differentiate between the two.

I would suggest you generate a unique signalling NAN that is not the compiler generated signaling NAN and is not (likely) a computer signaling NAN. Such as (single precision). Such as

0xFF89ABCD (single precision)
0xFFF789ABCDEF0123 (double precision)

These are signaling NANs, if they are inadvertently used in computations, the execution will stop.
These signaling NANs are unlikely to be generated (due to trap) unless something mistakenly overwrites this with trash data (with same values).

You would then treat (cast) the array as integer(4) or integer(8) as the case may be, and then search for magic number that indicates a free location.

Jim Dempsey

ereisch · ‎03-01-2021

@jimdempseyatthecove in this case, I am not looking for a computed NAN; I simply want if a NAN is used in an operation, for all subsequent (floating-point) operations to carry-forward that NAN. I have yet to decide if I want to implement it as a signaling NAN to catch immediately, as that may have other implications I may not want to get into.

This brings up an interesting quirk though: If I have a structure with default initializers as such:

TYPE, BIND(C) :: FOO_PACKET
  INTEGER(KIND=C_INT32_T) :: ID = -1
  REAL(KIND=C_FLOAT)      :: X  = Z'FFFFFFFF'
  REAL(KIND=C_FLOAT)      :: Y  = Z'FFFFFFFF'
END TYPE

...Intel Fortran doesn't complain, and happily compiles the routine that includes this block. However, if I compile with gfortran, it complains about the BOZ literal constants as being invalid in this context. A quick check on the GNU site indicates that the use of BOZ here is not standards-conforming.

So, which is right: GNU gfortran or ifort? Is this only valid syntax after a certain Fortran standard version?

andrew_4619 · ‎03-01-2021

warning #6915: Fortran 2018 does not allow data values on a type statement. [Z'FFFFFFFF'] . You need standards checking, as otherwise it is assuming an intel language extension.

ereisch · ‎03-01-2021

What flags are you using?

andrew_4619 · ‎03-01-2021

/stand:f18 but /stand should do it

ereisch · ‎03-01-2021

Hmm... "-stand" isn't in the ifort Linux manual page, but it's in the HTML documentation; that's why I missed it. Thanks

jimdempseyatthecove · ‎03-01-2021

Look at FAST_IS_NAN in the following.

This example is not vectorized, but I have vectorized it in a different project.

!  TestNaN.f90 
module MOD_NEGATIVE_IS_NAN
    
    integer, parameter :: IEEE754_SP_SIGN = Z'80000000'
    integer, parameter :: IEEE754_SP_EXPONENT = Z'7F800000'
    ! *** Due to compiler (Intel 2020u1/u2) optimization bug we cannot use a parameter
!    integer, parameter :: IEEE754_SP_MANTISSA = Z'007FFFFF'
    ! instead using memory stored variable works
    ! at some point this can be restored to parameter
    integer :: IEEE754_SP_MANTISSA = Z'007FFFFF'
   
    integer(8), parameter :: IEEE754_DP_SIGN = Z'8000000000000000'
    integer(8), parameter :: IEEE754_DP_EXPONENT = Z'7FF0000000000000'
    ! *** Due to compiler (Intel 2020u1/u2) optimization bug we cannot use a parameter
     integer(8), parameter :: IEEE754_DP_MANTISSA =  Z'000FFFFFFFFFFFFF'
    ! instead using memory stored variable works
    ! at some point this can be restored to parameter
!    integer(8) :: IEEE754_DP_MANTISSA =  Z'000FFFFFFFFFFFFF'
   
    ! Function that returns a negative integer (of size of REAL/DOUBLE PRECISION)
    ! when input argument is a NaN  (QNaN as SNaN would have caused an abort)
    ! these functions are elemental and capable of operating on vectors
    interface NEGATIVE_IS_NAN
        MODULE PROCEDURE NEGATIVE_IS_NAN_SP
        MODULE PROCEDURE NEGATIVE_IS_NAN_DP
    end interface NEGATIVE_IS_NAN
    
    interface FAST_IS_NAN
        MODULE PROCEDURE FAST_IS_NAN_SP
        MODULE PROCEDURE FAST_IS_NAN_DP
    end interface FAST_IS_NAN
    CONTAINS
    
    INTEGER(4) ELEMENTAL FUNCTION NEGATIVE_IS_NAN_SP(X)
        REAL, INTENT(IN) :: X
        NEGATIVE_IS_NAN_SP = TRANSFER(ABS(X),0)+IEEE754_SP_MANTISSA
    END FUNCTION NEGATIVE_IS_NAN_SP
    
    INTEGER(8) ELEMENTAL FUNCTION NEGATIVE_IS_NAN_DP(X)
        DOUBLE PRECISION, INTENT(IN) :: X
        NEGATIVE_IS_NAN_DP = TRANSFER(ABS(X),0_8)+IEEE754_DP_MANTISSA
    END FUNCTION NEGATIVE_IS_NAN_DP
    
    LOGICAL FUNCTION FAST_IS_NAN_SP(X)
        REAL, INTENT(IN) :: X
        FAST_IS_NAN_SP = (NEGATIVE_IS_NAN_SP(X) < 0)
    END FUNCTION FAST_IS_NAN_SP
    
    LOGICAL FUNCTION FAST_IS_NAN_DP(X)
        DOUBLE PRECISION, INTENT(IN) :: X
        FAST_IS_NAN_DP = (NEGATIVE_IS_NAN_DP(X) < 0_8)
    END FUNCTION FAST_IS_NAN_DP
end module MOD_NEGATIVE_IS_NAN
    
program TestNaN
   use, intrinsic :: ieee_arithmetic, only : ieee_value, ieee_quiet_nan, ieee_is_nan
   use, intrinsic :: ieee_exceptions, only : ieee_set_halting_mode, ieee_divide_by_zero
    use MOD_NEGATIVE_IS_NAN
    use omp_lib
    implicit none
    
    integer, parameter :: N = 100000000     ! some number larger than Last Level Cache
    integer, parameter :: reps = 50         ! number of times to scann for NaN
!    integer, parameter :: N = 1000     ! some number larger than Last Level Cache
!    integer, parameter :: reps = 5000000         ! number of times to scann for NaN
    real, allocatable :: SP(:)
    real :: SP_NAN, SP_R, SP_DR
    double precision, allocatable :: DP(:)
    real(8) :: DP_NAN, DP_R, DP_DR
    double precision :: T0, T1
    integer :: I, J
    LOGICAL :: FOUND_NAN
    INTEGER :: COUNT_NANS
   
    ALLOCATE(SP(N), DP(N))
    CALL RANDOM_NUMBER(SP)
    CALL RANDOM_NUMBER(SP_R)
    SP_DR = SP_R / 10.0
    where ( (SP > (SP_R-SP_DR)).and.(SP < (SP_R+SP_DR)) )
        SP = ieee_value( SP_R, ieee_quiet_nan )
    end where
    
    CALL RANDOM_NUMBER(DP)
    CALL RANDOM_NUMBER(DP_R)
    DP_DR = DP_R / 10.0
    where ( (DP > (DP_R-DP_DR)).and.(DP < (DP_R+DP_DR)) )
        DP = ieee_value( DP_R, ieee_quiet_nan )
    end where
    
    ! ASSURE CONSTANTS ARE CORRECT
    SP_NAN = TRANSFER(IEEE754_SP_EXPONENT+1234, 0.0)    ! Arbitrary QNaN
    FOUND_NAN = IEEE_IS_NAN(SP_NAN)
    PRINT *,'IEEE_IS_NAN(SP_NAN)', FOUND_NAN
    FOUND_NAN = (NEGATIVE_IS_NAN(SP_NAN) < 0)
    PRINT *,'NEGATIVE_IS_NAN(SP_NAN)', FOUND_NAN
    FOUND_NAN = FAST_IS_NAN(SP_NAN)
    PRINT *,'FAST_IS_NAN(SP_NAN)', FOUND_NAN
    
    ! ASSURE CONSTANTS ARE CORRECT
    DP_NAN = TRANSFER(IEEE754_DP_EXPONENT+1234_8, 0.0_8)    ! Arbitrary QNaN
    FOUND_NAN = IEEE_IS_NAN(DP_NAN)
    PRINT *,'IEEE_IS_NAN(DP_NAN)', FOUND_NAN
    FOUND_NAN = (NEGATIVE_IS_NAN(DP_NAN) < 0_8)
    PRINT *,'NEGATIVE_IS_NAN(DP_NAN)', FOUND_NAN
    FOUND_NAN = FAST_IS_NAN(DP_NAN)
    PRINT *,'FAST_IS_NAN(DP_NAN)', FOUND_NAN
    
  ! Make pass without report to precondition cache to some stable state
    COUNT_NANS = 0
    DO I=1,N
        SP(I) = SP(I) * 1.00001
        IF(IEEE_IS_NAN(SP(I))) COUNT_NANS = COUNT_NANS + 1
    END DO
    
    ! Now timed passes
    T0 = OMP_GET_WTIME()
    DO J=1, reps
        COUNT_NANS = 0
        DO I=1,N
            SP(I) = SP(I) * 1.00001
            IF(IEEE_IS_NAN(SP(I))) COUNT_NANS = COUNT_NANS + 1
        END DO
    END DO
    T1 = OMP_GET_WTIME()
    PRINT *, "IEEE_IS_NAN(SP(I))          ", T1-T0, COUNT_NANS

    ! Make pass without report to precondition cache to some stable state
    COUNT_NANS = 0
    DO I=1,N
        SP(I) = SP(I) * 1.00001
        IF(NEGATIVE_IS_NAN(SP(I)) < 0) COUNT_NANS = COUNT_NANS + 1
    END DO
    
    ! Now timed passes
    T0 = OMP_GET_WTIME()
    DO J=1, reps
        COUNT_NANS = 0
        DO I=1,N
            SP(I) = SP(I) * 1.00001
            IF(NEGATIVE_IS_NAN(SP(I)) < 0) COUNT_NANS = COUNT_NANS + 1
        END DO
    END DO
    T1 = OMP_GET_WTIME()
    PRINT *, "NEGATIVE_IS_NAN(SP(I)) < 0  ", T1-T0, COUNT_NANS
    
    
  ! Make pass without report to precondition cache to some stable state
    COUNT_NANS = 0
    DO I=1,N
        SP(I) = SP(I) * 1.00001
        IF(IEEE_IS_NAN(SP(I))) COUNT_NANS = COUNT_NANS + 1
    END DO
    
    ! Now timed passes
    T0 = OMP_GET_WTIME()
    DO J=1, reps
        COUNT_NANS = 0
        BLOCK
            INTEGER :: COUNT_NANS_4(0:15)
            INTEGER :: J
            COUNT_NANS_4 = 0
            DO I=1,N - MOD(N,16), 16
                DO J=0,15
                    SP(I+J) = SP(I+J) * 1.00001
                    IF(NEGATIVE_IS_NAN(SP(I+J)) < 0) COUNT_NANS_4(J) = COUNT_NANS_4(J) + 1
                END DO
            END DO
            DO J=0,15
                COUNT_NANS = COUNT_NANS + COUNT_NANS_4(J)
            END DO
            DO I=N-MOD(N,16)+1,N
                SP(I) = SP(I) * 1.00001
                IF(NEGATIVE_IS_NAN(SP(I)) < 0) COUNT_NANS = COUNT_NANS + 1
            END DO
        END BLOCK
    END DO
    T1 = OMP_GET_WTIME()
    PRINT *, "NEGATIVE_IS_NAN(SP(I)) COUNT_NANS_4", T1-T0, COUNT_NANS

    ! Make pass without report to precondition cache to some stable state
    COUNT_NANS = 0
    DO I=1,N
        SP(I) = SP(I) * 1.00001
        IF(FAST_IS_NAN(SP(I))) COUNT_NANS = COUNT_NANS + 1
    END DO
    
    ! Now timed passes
    T0 = OMP_GET_WTIME()
    DO J=1, reps
        COUNT_NANS = 0
        DO I=1,N
            SP(I) = SP(I) * 1.00001
            IF(FAST_IS_NAN(SP(I))) COUNT_NANS = COUNT_NANS + 1
        END DO
    END DO
    T1 = OMP_GET_WTIME()
    PRINT *, "FAST_IS_NAN(SP(I))          ", T1-T0, COUNT_NANS
    
    ! Make pass without report to precondition cache to some stable state
    COUNT_NANS = 0
    DO I=1,N
        DP(I) = DP(I) * 1.00001_8
        IF(IEEE_IS_NAN(DP(I))) COUNT_NANS = COUNT_NANS + 1
    END DO
    
    ! Now timed passes
    T0 = OMP_GET_WTIME()
    DO J=1, reps
        COUNT_NANS = 0
        DO I=1,N
            DP(I) = DP(I) * 1.00001_8
            IF(IEEE_IS_NAN(DP(I))) COUNT_NANS = COUNT_NANS + 1
        END DO
    END DO
    T1 = OMP_GET_WTIME()
    PRINT *, "IEEE_IS_NAN(DP(I))          ", T1-T0, COUNT_NANS
    
    ! Make pass without report to precondition cache to some stable state
    COUNT_NANS = 0
    DO I=1,N
        DP(I) = DP(I) * 1.00001_8
        IF(NEGATIVE_IS_NAN(DP(I)) < 0_8) COUNT_NANS = COUNT_NANS + 1
    END DO
    
    ! Now timed passes
    T0 = OMP_GET_WTIME()
    DO J=1, reps
        COUNT_NANS = 0
        DO I=1,N
            DP(I) = DP(I) * 1.00001_8
            IF(NEGATIVE_IS_NAN(DP(I)) < 0_8) COUNT_NANS = COUNT_NANS + 1
        END DO
    END DO
    T1 = OMP_GET_WTIME()
    PRINT *, "NEGATIVE_IS_NAN(DP(I)) < 0_8", T1-T0, COUNT_NANS
    
    ! Make pass without report to precondition cache to some stable state
    COUNT_NANS = 0
    DO I=1,N
        DP(I) = DP(I) * 1.00001_8
        IF(NEGATIVE_IS_NAN(DP(I)) < 0_8) COUNT_NANS = COUNT_NANS + 1
    END DO
    
    ! Now timed passes
    T0 = OMP_GET_WTIME()
    DO J=1, reps
        COUNT_NANS = 0
        BLOCK
            INTEGER(8) :: COUNT_NANS_8(0:7)
            INTEGER :: J
            COUNT_NANS_8 = 0_8
            DO I=1,N - MOD(N,8), 8
                DO J=0,7
                    DP(I+J) = DP(I+J) * 1.00001_8
                    IF(NEGATIVE_IS_NAN(DP(I+J)) < 0_8) COUNT_NANS_8(J) = COUNT_NANS_8(J) + 1_8
                END DO
            END DO
            DO J=0,7
                COUNT_NANS = COUNT_NANS + COUNT_NANS_8(J)
            END DO
            DO I=N-MOD(N,8)+1,N
                DP(I) = DP(I) * 1.00001_8
                IF(NEGATIVE_IS_NAN(DP(I)) < 0_8) COUNT_NANS = COUNT_NANS + 1
            END DO
        END BLOCK
    END DO
    T1 = OMP_GET_WTIME()
    PRINT *, "NEGATIVE_IS_NAN(DP(I)) COUNT_NANS_8", T1-T0, COUNT_NANS
    
    ! Make pass without report to precondition cache to some stable state
    FOUND_NAN = .FALSE.
    DO I=1,N
        DP(I) = DP(I) * 1.00001_8
        IF(FAST_IS_NAN(DP(I))) FOUND_NAN = .TRUE.
    END DO
    
    ! Now timed passes
    T0 = OMP_GET_WTIME()
    DO J=1, reps
        COUNT_NANS = 0
        DO I=1,N
            DP(I) = DP(I) * 1.00001_8
            IF(FAST_IS_NAN(DP(I))) COUNT_NANS = COUNT_NANS + 1
        END DO
    END DO
    T1 = OMP_GET_WTIME()
    PRINT *, "FAST_IS_NAN(DP(I))", T1-T0, COUNT_NANS
    
end program TestNaN

The above should get you started.

Jim Dempsey

jimdempseyatthecove · ‎03-01-2021

Test NaN on Core i7 2600K

 IEEE_IS_NAN(SP_NAN) T
 NEGATIVE_IS_NAN(SP_NAN) T
 FAST_IS_NAN(SP_NAN) T
 IEEE_IS_NAN(DP_NAN) T
 NEGATIVE_IS_NAN(DP_NAN) T
 FAST_IS_NAN(DP_NAN) T
 IEEE_IS_NAN(SP(I))             42.8022058000788         16763089
 NEGATIVE_IS_NAN(SP(I)) < 0     4.22263870015740         16763089
 NEGATIVE_IS_NAN(SP(I)) COUNT_NANS_4   2.59221769985743         16763089
 FAST_IS_NAN(SP(I))             12.2822704999708         16763089
 IEEE_IS_NAN(DP(I))             45.0686576000880         16628001
 NEGATIVE_IS_NAN(DP(I)) < 0_8   9.07150449999608         16628001
 NEGATIVE_IS_NAN(DP(I)) COUNT_NANS_8   5.25665029999800         16628001
 FAST_IS_NAN(DP(I))   12.9640327000525         16628001

Jim Dempsey

jimdempseyatthecove · ‎03-01-2021

I should mention that the compiler vectorized the COUNT_NANS_4 and COUNT_NANS_8.

Orignally written on system with AVX512F, but as run above on system with AVX (AVX1).

Jim Dempsey

jimdempseyatthecove · ‎03-02-2021

I haven't experimented with this...

integer and logical are interchangeable

You might try:

    
    LOGICAL FUNCTION FAST_IS_NAN_SP(X)
        REAL, INTENT(IN) :: X
        FAST_IS_NAN_SP = ISHL(NEGATIVE_IS_NAN_SP(X), -31)
    END FUNCTION FAST_IS_NAN_SP
    
    LOGICAL FUNCTION FAST_IS_NAN_DP(X)
        DOUBLE PRECISION, INTENT(IN) :: X
        FAST_IS_NAN_DP = ISHL(NEGATIVE_IS_NAN_DP(X), -63)
    END FUNCTION FAST_IS_NAN_DP

Try that, and report back.

In the codes that I used this in, I used the return of the NEGATIVE_IS_NAN as opposed to the LOGICAL.

Jim Dempsey

JohnNichols · ‎03-01-2021

If zero is a reasonable number, but you want to know it is calculated use a weird random number that is unlikely to set your array -- extremely unlikely

1.23456789

ereisch · ‎03-01-2021

@JohnNichols no, floating point equality tests are almost universally discouraged (unless comparing to 0.0), as the stored value goes through a discretization step when it is fetched to the floating point registers, and isn't guaranteed to exactly match your compile-time constant. You could define an explicit bit value, but then you might as well just use NANs.

JohnNichols · ‎03-02-2021

@ereisch noted Unfortunately, they are a little more frequent than usual, as this particular array we're initializing to NaN at startup (so as to not confuse with a "perfectly valid" 0.0). The code block in question is to find the first "unused" member of the array.

----------------------------------------------------------------------------------------------------------------

In the types of coding I work on as an engineer there are usually a valid range of numbers for the results.

It is then reasonably easy to pick a number to initialize that is miles away from the usual range that if it pops up, it is some form of error.

The problem with the NaN is occasionally you get them in the calculations - not often but statistically it is an annoying potential problem.

Either way there are challenges, but it is fun to follow this thread.

ereisch · ‎03-02-2021

So this brings up a related question -- is there a "standards compliant" method of initializing derived type members to a NAN value?

As pointed out earlier, the following is not compliant (but works with ifort):

TYPE, BIND(C) :: FOO_PACKET
  REAL(KIND=C_FLOAT) :: X = Z'FFFFFFFF'
  REAL(KIND=C_FLOAT) :: Y = Z'FFFFFFFF'
END TYPE

Similarly, the following end-around also seems illegal (though accepted in ifort), as BOZ literals are not allowed as arguments to 'transfer':

TYPE, BIND(C) :: FOO_PACKET
  REAL(KIND=C_FLOAT) :: X = TRANSFER(Z'FFFFFFFF', 1.0)
  REAL(KIND=C_FLOAT) :: Y = TRANSFER(Z'FFFFFFFF', 1.0)
END TYPE

I read that there might be a compliant way to override the default initializer, and then BOZ literals might be able to used inside that (since it would just be a regular function), but I can't find a good documentation source on how to do this. Would the derived type then need to be wrapped up in a module to do this, or is there a less intrusive way of going about it?

JohnNichols · ‎03-02-2021

Taking you code and placing it in a standard program

program Console11
    use Base

    implicit none
    
    real X1
    logical yes
    TYPE, BIND(C) :: FOO_PACKET
      REAL(KIND=dp) :: X = Z'FFFFFFFF'
      REAL(KIND=dp) :: Y = Z'FFFFFFFF'
    END TYPE
    
    Type(FOO_PACKET) :: P1
    yes = .false.
    X1 = 1.0
    ! Variables
    ! Body of Console11
    print *, 'Hello World'
    yes = ISNAN(X1)
    print *,yes
    yes = ISNAN(P1%X)
    print *,yes
    end program Console11

The X and Y are assigned as

	Name	Value	Type
	P1%X	2.121995790471207D-314	REAL(8)

Did I make a mistake in understanding what you are saying?

if I set the z'FFF to an integer i it is equal to -1.

JohnNichols · ‎03-02-2021

If I do not assign values the type has two zero values. You can never be sure what Intel Fort will set variables to if you just create them -- it is a pain when you bring in old programs that assumed some value.