Underflow/Overflow Errors

QUESTION: I keep getting floating underflow/overflow errors when I run my programs, but the math results seem to be correct. Can these warnings be ignored or do they mean I'm doing something wrong? If they are nonsense, can I turn the darn things off?

ANSWER: Our answer today is provided by the ever-reliable math guy, Craig Markwardt, along with George N. White III and Paul van Delst. First, Craig.

In all likelihood you can ignore the overflow and underflow messages. You can use !EXCEPT=2 to find out where in your program the exceptions are being created. You can also turn off exceptions using !EXCEPT=0, if that is crucial.

While I have argued in the past that most people don't need underflow errors, and that they should be silent, other folks on the newsgroup have argued that we should strive to avoid them, so such errors should be printed. Leaving the question of right verses wrong on error messages aside for the moment, I indeed think it is important to avoid over and underflows.

My guess is that you are getting both when you use EXP() in the Fermi distribution. To avoid this you can use some simple techniques. One idea is to use thresholding to keep all values in-bounds, like this, EXP((X) > (-1e-38) < 1e38). That is not really satisfactory though, because sometimes you *want* the effect of an underflow. That is perhaps best solved using the WHERE() command to locate extremal values and treat them specially.

Paul mostly concurs.

First test your code after setting !EXCEPT=2. This will tell you on what lines of code you are getting the overflow/underflow. Then you can alter your algorithm to avoid them depending on your needs. For example, if the numerical precision is a good enough tolerance level, you can do something like this.
  ; Set some tolerance level, with double keyword just in case.
  tolerance = ( machar(double=double) ).eps

  if ( x LT tolerance ) then $
    result = KEYWORD_SET( double ) ? 0.0d0 : 0.0 $
  else $
    result = y / x
I think a lot of people will think this is probably overkill, and maybe it is for what you want to do. My recent experience has shown, in general, that you can ignore the underflow errors if (1) you are sure your code will always be executed in the same regime, and (2) you don't care if you code crashes and burns, or produces a crappy number every now and again.

I say this because the global forecast model I use here runs up to about the 2hPa pressure level. At these high altitudes the amount of water vapour is negligible. But it's not negligible compared to numerical precision. I prototyped, in IDL, the forward, then tangent-linear (TL), then adjoint (AD) model of the radiative transfer code used to simulate the satellite measurements in the forecast model - which uses up the FOURTH power of the integrated water vapour amount. To test the adjoint model, you expect agreement between the TL and AD modeuls to numerical precision (in double precision). I was getting huge errors at the top of the atmosphere for some satellite channels in the AD model. The reason was that my test code went up to 0.005hPa where the integrated water vapour amounts are next to nothing and (next to nothing)^4 in a denominator was catastrophic. I put in a test (similar to the above code) and the errors went away, no more underflow/overflow errors. The forecast model doesn't go up as high as I tested the radiative transfer code....yet. But in the future it may (it recently got bumped up to the 0.1hPa level).

So, like I said before, it depends on what you want to do. The operational forecast model always has to run. It can't crash with an error or produce a crappy number unless it's flagged as crappy.

Just my opinion of course, but I treat all floating point errors as serious.

George weighs in on the same topic.

"Seems to be correct" only offers some hope that when you understand the reason for underflow you will realize that it doesn't affect your results. Even then, you may want to alter the program to avoid underflow for two reasons.

1. Underflow is often an expensive way to set a variable to zero. When the source of underflow is understood it may become obvious that a significant chunk of calculation isn't needed and shouldn't be used (e.g., by adding a range test to omit the section where underlow occurs when the inputs are "out of bounds").

2. Modern hardware with branch prediction and combined floating point ops is generally tuned for peak performance in typical cases, and can fall down very badly when handling exceptions. Even if your current program performs adequately, once you have analyzed the underflow you may want to document it in case you encounter performance problems in the future.

And finally Craig with one more trick.

On the other hand, the performance of IDL falls down rather badly when dealing with conditional tests on large arrays, especially when FOR loops cannot be avoided. Even using WHERE() usually makes a pretty large performance hit.

One little trick for avoiding underflows in exponentials might go like this. If you wish to compute Y = EXP(-X), for large X, then you can usually use a "mask" variable to avoid the underflow.
   MASK = (X LT UPPER_LIMIT)
   Y = MASK * EXP(-X*MASK)
Note that this doesn't catch overflows which will be more disastrous.

QUESTION: I've encountered a really interesting--though annoying--problem. I'm trying to use the command Exp(variable), where variable is greater than 88 or so. My calculator can do it, but IDL can't. It gives me a "floating underflow" error. Can anyone explain this to me?

ANSWER: This answer to the IDL newsgroup question above was given by William Clodius of Los Alamos National Lab.

"Computers store and retrieve data from memory that (except for a few eperimental machines that will never run IDL) invariably treats memory as a large linear sequence of bits, 2 valued entities. For convenience these bits are grouped into chunks called bytes that are almost invariably 8 bits long. IDL will not run on machines without this definition of bytes. Bytes can be thought of as integers whose values range from -128 to 127 or from 0 to 255 depending on context. In other contexts, on processors in countries with reasonably small character sets, these integer values are treated as mapping to individual characters. In certain contexts bytes are grouped together in sets of two and treated as integers whose values may range from -32768 to 32767. This is the default IDL integer type, called SHORT. In other contexts bytes are grouped together in sets of four and treated as integers whose values range from -2147483648 to 2147483647. This is the IDL integer type called LONG. Most computers provide only these four integer types and they are the only ones IDL provides, although a few machines now also provide integers defined in terms of eight bytes. It is possible for languages to define integers beyond those intrinsic to the machine, but they have to be emulated through software and suffer a significant performance disadvantage. IDL does not define such integers. (Note the value ranges above assume a twos complement implementation of integers).

"In many cases it is desirable to have locations in memory that are treated as representing rational numbers. The first thing to realize is that these rational numbers will at best be approximations to irrational numbers because of the finite memory of the machine. On the vast majority of machines these representations of rational numbers will be in terms of four or eight byte numbers, where the four byte numbers are IDL's FLOAT data type and the eight byte numbers are IDL's DOUBLE data type. The second thing to realize is that, because of the finite size of the storage for these numbers, only a (relatively small) fixed set of rational number can be represented by these types. The machine treats these numbers as consisting of two parts a mantisa, which can be thought of as an integer value, and an exponent, which represents a power of two that multiplies the integer. The third thing to realize is that the useage of a power of two for the exponent means that the processor can not exactly represent many numbers that novices assume are simple enought to be exactly represented, i.e., 1./3. or 0.2. The power of two is desirable because it makes the math easy to implement and ensures uniform and predictable coverage of the real number space. Again it is possible for languages to provide floating point math other than what is provided by the machine, but it requires emulation in software with a significant performance impact. IDL only provides the machine based representations. (NOTE: because calculators often deal with few numbers at a time, they can often devote more storeage space per number and hence can represent numbers with a higher degree of accuracy than most computers.)

"The vast majority of processors now define their rational number approximate representation in terms of the IEEE Standard For Binary Floating Point Arithmetic (ANSI/IEEE Std 754-1985). This defines the processor's behavior for certain conditions so that it can reliably report real or potential errors. These errors include generating an infinity (e.g., dividing a finite number by zero), generating an undefined number (e.g., dividing zero by zero), overflow (e.g., mutiplying two large numbers to generate a number too large to be represented), and underflow (e.g., dividing a large number into a small number and generating a finite number that is too small to be represented accurately). Machines with IEEE arithmetic cannot represent e^89 as a FLOAT, but can represent it (in some relatively good approximation) as a DOUBLE. Because you got an underflow error I suspect your calculation either involved e^(-89) or 1./e^89 and not just e^89.

"William Kahan's page is a useful source of additional information."

   http://HTTP.CS.Berkeley.EDU/~wkahan/