topic Re: need zero / x equals zero in Intel® Fortran Compiler

need zero / x equals zero

jim_cox — Mon, 19 Oct 2009 04:12:18 GMT

Apparently IFort thinks zero / x = NaN

Our program logic expects zero / x = zero

How do we persuade IFort - or is writing our own exception handler the only only option ?

advTHANKXance

Jimmy

Re: need zero / x equals zero

DavidWhite — Mon, 19 Oct 2009 04:31:46 GMT

Quoting - jim.cox

Apparently IFort thinks zero / x = NaN

Our program logic expects zero / x = zero

How do we persuade IFort - or is writing our own exception handler the only only option ?

advTHANKXance

Jimmy

Jimmy,

Can you please supply a sample of code demonstrating this issue? 0 / x should be zero unless x is zero, when NaN is the appropriate response.

David

Re: need zero / x equals zero

TimP — Mon, 19 Oct 2009 04:33:38 GMT

If you have 0./0. that produces NaN according to IEEE standard. I don't see how an exception handler is relevant. If you have an expression involving divide, and you wish to set a zero result when the numerator is zero, regardless of divisor, that should be easily accomplished by testing before divide.

Re: need zero / x equals zero

jimdempseyatthecove — Mon, 19 Oct 2009 15:37:11 GMT

Even if you compile with "/assume:protect_parens" and change your line to

result =numerator * (1.0/denominator)

When denominator is NaN, the result will be NaN. When denominator = 0.0 you will either get a div 0 error or NaN for (1.0/denominator). However

result = numerator
if(result .ne. 0.0) result = result / denominator

Will produce 0.0 when numerator=0.0 regardless of value of denominator.

Jim Dempsey

Re: need zero / x equals zero

TimP — Mon, 19 Oct 2009 16:46:13 GMT

Quoting - jimdempseyatthecove

Depending on your detailed requirements, you might prefer Jim's suggestion (as we would have 20 years ago), or
result = MERGE(numerator,numerator/denominator,abs(numerator < TINY(numerator)))

Re: need zero / x equals zero

jim_cox — Mon, 19 Oct 2009 19:38:02 GMT

Quoting - David White

Jimmy,

Can you please supply a sample of code demonstrating this issue? 0 / x should be zero unless x is zero, when NaN is the appropriate response.

David

we are getting the slope of a line

when

VV(2,1) = 2.0
VV(1,1)= 2.0
QQ(2,1)=500.0
QQ(1,1)=0.0

SLOPEV(JP,I) = (VV(JP+1,I)-VV(JP,I)) / (QQ(JP+1,I)-QQ(JP,I))

= NaN

:(

Re: need zero / x equals zero

jim_cox — Mon, 19 Oct 2009 20:39:07 GMT

Quoting - jimdempseyatthecove

However

result = numerator
if(result .ne. 0.0) result = result / denominator

Will produce 0.0 when numerator=0.0 regardless of value of denominator.

That is elegant code thank you Jim :)

But will it work if the denominator is zero?

Better perhaps with...

result = numerator
if(result .ne. 0.0.and.denominator.ne.0.0) result = result / denominator

Wonder how many places I have to change my code?

There are only 2060 divide statements in this one program :(

Jim

=mjc=

.

Re: need zero / x equals zero

TimP — Mon, 19 Oct 2009 20:56:40 GMT

For what it's worth, behavior of your program may depend on whether you set /Qftz- for compilation of main program, either directly, or indirectly by /fp:source

Re: need zero / x equals zero

jim_cox — Mon, 19 Oct 2009 23:35:43 GMT

Quoting - tim18

For what it's worth, behavior of your program may depend on whether you set /Qftz for compilation of main program, either directly, or indirectly by /fp:source

Thankx,

Had a quick play with this setting

But in my case that doesnt seem to make any difference

Re: need zero / x equals zero

DavidWhite — Mon, 19 Oct 2009 23:54:14 GMT

Quoting - jim.cox

we are getting the slope of a line

when

VV(2,1) = 2.0
VV(1,1)= 2.0
QQ(2,1)=500.0
QQ(1,1)=0.0

SLOPEV(JP,I) = (VV(JP+1,I)-VV(JP,I)) / (QQ(JP+1,I)-QQ(JP,I))

= NaN

:(

Jim,

Assuming that JP and I are set correctly for the values you quote here, the only reason I can think that you could be getting NaN is that perhaps you are exceeding your array bounds, and accessing other memory.

David

Re: need zero / x equals zero

jim_cox — Tue, 20 Oct 2009 00:25:25 GMT

Quoting - David White

Both JP = 1 and I = 1

and SLOPEV is declared with

INTEGER, PARAMETER :: NP=10

COMMON/COEF/QQ(NP+1,99),VV(NP+1,99),SLOPEV(NP+1,99),NPT(99)

REAL QQ, VV, SLOPEV

Slopev is happily initialised to zero elsewhere

Bounds checking is enabled and not complaining

There is so much "weirdness" going on with this code that I figure it has to be some sort of memory issue.

And I'm ripping my hair out trying to work out what / where it is

Ah the joys of an old old old codebase...

Re: need zero / x equals zero

DavidWhite — Tue, 20 Oct 2009 00:58:37 GMT

Jim,

You said "There is so much "weirdness" going on with this code that I figure it has to be some sort of memory issue."

Perhaps it is notthese arrays causing the problem, but they are being trampled on from elsewhere.

Are youusing IMPLICIT NONE and ensuringthat all variables are declared and initialized properly -- may be agood place to start -- could identify variables which are not declared in the way you think they are.

Happy hunting.

David

Re: need zero / x equals zero

jim_cox — Tue, 20 Oct 2009 01:24:01 GMT

Quoting - David White

The codebase is such that much of the old code was written expecting implicit typing, and often with no declaration at all.

Personally I'm fussier, and the newer code is 'better'

I'm also now compiling with warnings for undeclared, unused, unalligned and interfaces etc all on - and fixing them as I go.

The routines displaying this beghaviour are compiling without complaint

more....

I tried IMPLICIT NONE

But that stopped compilation at the line

L_DEL = GDLY(.TRUE.,.5,ANGL,AS*1000/3600.,DS*1000/3600.,R)

Where GDLY is one of our library functions

Go figure...

Thankx for your input

Re: need zero / x equals zero

DavidWhite — Tue, 20 Oct 2009 01:43:16 GMT

Quoting - jim.cox

Jim,

if you add GDLY to a type statement, I think this should be resolved.

D

Re: need zero / x equals zero

jimdempseyatthecove — Tue, 20 Oct 2009 04:21:11 GMT

When numerator .ne. 0.0 and denominator .eq. 0.0 just what sort of number do you want produced?

The code I listed will produce either a div by 0 exception or a NaN. Either of which are commonly accepted results. Returning 0.0 is not (not even for 0.0/0.0). The code I produces returns 0.0 for 0.0/x where x is anyting including 0.0 and NaN. This is the behavior you asked for. Y/0.0 where Y is not 0.0 produces what in your number scheme?

Wherever your program will produce a divide by 0, you must protect that section of code. This is likely a very small set of your 2060 statements. And if the NaN's are produced in places they ought not to be produced then this is an indication of an error elsewhere in your program (i.e. something else bunged up your data).

Jim Dempsey.

Re: need zero / x equals zero

portoon — Tue, 20 Oct 2009 08:18:56 GMT

Quoting - jim.cox

Apparently IFort thinks zero / x = NaN

Our program logic expects zero / x = zero

How do we persuade IFort - or is writing our own exception handler the only only option ?

advTHANKXance

Jimmy

I was concerning for this issue, always.

When I use following code, I got a good results.

It's idea came from tim18 andJim Dempsey.

Thanks tim18 and Jim Dempsey

Whatever a, b it's result seems to be reasonable.

a=0.0, huge, tiny, -huge, -tiny; b=0.0, huge, tiny, -huge, -tiny...etc.

____________________________________________________________

implicit none
real*8 :: a, b, c

a = 0.0d0 ! 0.0d0, tiny(a), huge(a)...
b = 0.0d0 ! 0.0d0, tiny(b), huge(b)...

c = merge( a ,a/b,dabs(a)< tiny(a))
c = merge(dsign(huge(c),c),c ,dabs(c)>huge(c))

write(*,*) a/b, c

stop
end

____________________________________________________________

Re: need zero / x equals zero

portoon — Tue, 20 Oct 2009 08:44:11 GMT

Quoting - David White

Jim,

if you add GDLY to a type statement, I think this should be resolved.

D

I was concerning for this issue, always.
When I use following code, I got a good results.
It's idea came from tim18 and Jim Dempsey.
Thanks tim18 and Jim Dempsey

If a or b is not NaN, it's result seems to be reasonable.
a=0.0, huge, tiny; b=0.0, huge, tiny...etc.

____________________________________________________________
implicit none
real*8 :: a, b, c

a = 0.0d0 ! 0.0d0, tiny(a), huge(a)...
b = 0.0d0 ! 0.0d0, tiny(b), huge(b)...

c = merge( a ,a/b,dabs(a)< tiny(a))
c = merge(dsign(huge(c),c),c ,dabs(c)>huge(c))

write(*,*) a/b, c

stop
end
____________________________________________________________

Re: need zero / x equals zero

jimdempseyatthecove — Tue, 20 Oct 2009 14:00:21 GMT

Jim,

I forgot to mention that in both my suggestion and in David's (even if merge is inlined) is that compiler optimizations may perform the division regardless of if the if test passes or fails. The newer versions of the instruction set have a conditional move. It is often faster to perform the division speculatively then perform a conditional move without use of branch instruction. In this situation (conditional move) you may/will experience divide by zero.

/Qfp-speculation:off

might help (for the source files with the questionable statements).

And use /fpe:1 or /fpe:3

Note, use of /fpe:1 or /fpe:3 will generate NaN's and asumes it is your responsibility for detection of floating point errors.

Also,

If you are the programmer and NOT the research scientist using the progrma, then on every instance where you assume you want this test, you must verify with the scientist as to what result to generate. Do not assume 0.0 is the desired value. The correct value could be any of

error/abort
0.0
TINY (+/-)
NaN
denormal number

Each case may (generaly does) have a different set of rules.

Jim Dempsey

Re: need zero / x equals zero

jim_cox — Wed, 21 Oct 2009 19:27:14 GMT

Quoting - jimdempseyatthecove

Yes Jim - I'm the current mug programer tending this simple old codebase (sigh)

As stated in the title we need 0.0 / X = 0.0

If X = 0 crashes, thats ok

Thankx for the suggestiions re the switches. I'm wondering if some of the floating point options might influence. I'll keep playing...

We are trying to move to the IFort compiler so that we can take advantage of the multi-porcessor support, so any suggestions that improved reliability would be helpful.

Cheers

Jim

Re: need zero / x equals zero

jimdempseyatthecove — Wed, 21 Oct 2009 22:21:44 GMT

When I converted a 700 file project from F77 to F95 (asside from some wierdness with DATA statements) conversion went rather well excepting for...

When converting from COMMON to MODULE some problems were exposed relating to EQUIVILENCE statements. The conversion to modules caught several problems relating to this but in the end was well worth the effort. The biggest initial improvement was in incorporationg allocatable arrays and use of pointers. The allocatable arrays meant the internal storage requirements could adjust according to input file parameters as opposed to being fixed at compile time to worst case. The pointers eliminated a lot of block moves of data in/out of working sets. The initial conversion attained 10x performance boost (principaly due to the appropriate use of pointers). Following that, compiler optimizations providing vectorization, then adding OpenMP. End result was 4 core system showed 40x+ improvement and provided for input file sizing of internal arrays. Your Milage May Vary.

One of the "hacks" I used to aid in the conversion was to extensively use the FPP. While keeping the majority of the code untouched, I could have the FPP redirect the once COMMON variables to use allocatable array and pointers. e.g.

#define SEGLEN pTether.pFiniteSolution.rSEGLEN

Placed in a xxx.INC include file permitted the code to continue to use SEGLEN(JSEG) without edits. Debugging was a little strange seeing that the InteliSence saw SEGLEN but the compiler saw pTether.pFiniteSolution.rSEGLEN. However, debugging was minimalized due to virtually no edits to working code.

There were several other conversion helpers the FPP provide in the conversion project.

Try to use all the tricks available - make life easy.

Jim Dempsey