I enabled floating point exception trapping using
the ieee_exceptions module in a code I'm working with,
and started getting some unexpected test suite failures.
I tracked them down to a merge statement of the form
real :: dx(:)
dx = merge(0.0, 1.0/dx, abs(dx) < 1.0e-3)
This particular compiler (pathscale) apparently
fully evaluates the array expression 1.0/dx, which
is throwing an exception when a 0-valued component
of dx is encountered.
I'm not actually surprised that the argument expression
was fully evaluated; that would be required if
merge were a user-function. Do intrinsic functions
have special rules in this regard?
That got me to wondering why this bit of code doesn't
cause any problems with the other compilers we use,
which do trap on fp exceptions by default. Poking
around with the nag compiler, I've found that it
compiles the merge to a loop that only evaluates those
elements of the array expression that are used, and
so never tries to do a 1.0/0.0. I'd guess that the
other compilers are doing the same thing.
I actually prefer this behavior (in this case at least),
but I'm suspicious that this may not be strictly kosher.
-Neil
--
Neil Carlson
Computational Physics and Methods (CCS-2)
Los Alamos National Laboratory
PO Box 1663, MS D413 · Los Alamos, NM 87545
505.665.6386 · 505.665.4972 (fax)
|