How to deal with "missing points" in arrays
in [Fortran]
|
From: robert.corbett on 9 Jun 2010 23:36 On Jun 9, 2:41 pm, Gib Bogle <g.bo...(a)auckland.no.spam.ac.nz> wrote: > glen herrmannsfeldt wrote: > > I don't recommend the LOGICAL variable method, unless it is > > necessary to have all REAL values legal. If you need portability > > to older compilers, you could do conditional compilation on a > > test for NaN or 9.9999e30 (fits in single precision on most > > machines), or 9.9999e99 (in double precision on many machines). > > Is 9.9999e30 exactly representable? Anyway, why not use 1.0e30, which > definitely is? I like your idea of incrementally modifying true values that hit > the "missing point" value. 1.0E30 is not exactly representable as an IEEE Std. 754 single- precision value. It is not exactly representable as an IEEE Std. 754 double- precision value. IEEE Std. 754's single-precision format has 24 significand bits. The value 1.0E30 is 2^30 x 5^30. The representation of the 2^30 portion requires no significand bits. The representation of 5^30 requires 70 bits. The IEEE Std. 754 single-precision value that is the closest to 1.0E30 is 1000000015047466219876688855040. Bob Corbett
From: William Clodius on 10 Jun 2010 00:14 deltaquattro <deltaquattro(a)gmail.com> wrote: > Hi, > > this is really more of a "numerical computing" question, so I cross- > post to sci.math.num.analysis too. I decided to post on > comp.lang.fortran, anyway, because here is full of computational > scientists and anyway there are some sides of the issue specifically > related to Fortran language. > > The problem is this: I am modifying a legacy code, and I need to > compute some REAL values which I then store in large arrays. Sometimes > it's impossible to compute these values: for example, think of > interpolating a table to a given abscissa, it may happen that the > abscissa falls outside the curve boundaries. I have code which checks > for this possibility, and if this happens the interpolation is not > performed. However, now I must "store" somewhere the information that > interpolation was not possible for that array element, and inform the > user of it. Since the values can be either positive or negative, I > cannot use tricks like initializing the array element to a negative > values. > <snip> Is there more than one way the interpolation can fail and shoult the user be aware of these distinctions? If so then I suggest a corresponding flag array of chararacters of length 1. Define character parameters of the same length, e.g., TOO_LOW, TOO_HIGH, VALID, ... with distinct values and assign to the duplicate array as the proper state is recognized. Consider initializing the character array element to a special distinct value, say UNPROCESSED. Also consider assignining an NaN to invalid elements so that if you fail to to use the flag array to determine the processing flow the failure will be obvious. -- Bill Clodius los the lost and net the pet to email
From: Gib Bogle on 10 Jun 2010 03:57 robert.corbett(a)oracle.com wrote: > > 1.0E30 is not exactly representable as an IEEE Std. 754 single- > precision > value. It is not exactly representable as an IEEE Std. 754 double- > precision > value. IEEE Std. 754's single-precision format has 24 significand > bits. > The value 1.0E30 is 2^30 x 5^30. The representation of the 2^30 > portion > requires no significand bits. The representation of 5^30 requires 70 > bits. > The IEEE Std. 754 single-precision value that is the closest to 1.0E30 > is > 1000000015047466219876688855040. > > Bob Corbett Right, I wasn't thinking straight. Very silly of me.
From: Gib Bogle on 10 Jun 2010 04:00 robert.corbett(a)oracle.com wrote: > 1.0E30 is not exactly representable as an IEEE Std. 754 single- > precision > value. ... In any case I was thinking that it would make sense to choose as the "missing point" value a number that is exactly representable - but does it really make any difference? Maybe not.
From: Tobias Burnus on 10 Jun 2010 04:22
On 06/10/2010 01:55 AM, Richard Maine wrote: > glen herrmannsfeldt <gah(a)ugcs.caltech.edu> wrote: >> The usual >> NaN test of (x.ne.x) is likely to be optimized out by some compilers. Well, it also depends on the options you use. For GCC (gcc, g++, gfortran, but also g95), if you compile with -ffast-math (or -Ofast), you imply the option "-ffinite-math-only Allow optimizations for floating-point arithmetic that assume that arguments and results are not NaNs or +-Infs." Thus, if you rely on NaN, make sure you don't tell the compiler that it can assume that you don't have any. On some systems, one also needs to enable IEEE support (which adds some software checks as the hardware does not have full IEEE support); in case of GCC that's -mieee which is seemingly used by Alphas and SH. > I don't consider that the "usual" test. In fact, I strongly recommend > against it. I consider the "usual" test to involve invoking a separate > function for the purpose. For today's compilers, that would be the > appropriate IEEE function. (Don't most of them do that by now?) I think most do, but gfortran does not yet do so.* Cf. also http://fortranwiki.org/fortran/show/Fortran+2003+status Tobias * For gfortran implementing the Fortran 2003 IEEE support is not immediately planned - currently the focus is on object-oriented programming [fixing bugs, adding FINAL and CLASS(*), etc.]. gfortran supports "(x /= x)" and also the somewhat common vendor intrinsic "isnan" - but not yet F2003's "ieee_is_nan". |