Fixing NaN Problem in f90 Do Loop

[fun]

Hi, I encountered a problem in f90 today. The problem is that somewhere in middle of the do loop, a variable goes larger than exp(-307). I have to invert that variable, then I am getting NaN. So f90 is working until the variable reaches exp(-307), but after that it is becoming infinity, so I am not able to invert it. Could you please able to help in fixing this issue.

 

[jedishrfu]

usually when you hit something like this you jump to the next higher precision datatype like from real to double but in this case it looks like you've hit that limit leaving three choices:

1) shift to logs

2) rework your calculations some that your numbers are in a more reasonable range such as +/-7 ie factor out 10^300 and put it back in later at the end.

3) use a big number library if one is available, but this is usually much slower than working with doubles.

 

[fun]

I am already using the double precision. The third way, I do not know how to implement it. Could you please give some source links or examples. The second way I tried, but its a scientific calculation, so i do not have any control over that. Actually i need to multiply a matrix 20 times. Thats whre i got this error.

Regards,
Fun.

 

[jedishrfu]

http://en.wikipedia.org/wiki/IMSL_Numerical_Libraries

the IMSL may have something on it although I couldn't find a direct reference.

 

[jedishrfu]

another thought is to intersperse the multiply and divide operations in your calculations if possible ie

x * y * z / (a * b * c) = x / a * y / b * z / c

this prevents any intermediate result from getting too large.

 

[chiro]

Hey Fun and welcome to the forums.

If you can read C++, you will find plenty of examples of big integer libraries where you can get the source code.

Ultimately it will depend on what operations you wish to use. If you are using only arithmetic then you will find plenty, but if you are using exponentials, logarithms or something else then you might have to end up buying a decent commercial library like MIRACL if you want to use it commercially.

What operations (arithmetic, exponentiation, logarithms, etc) are you using?

 

[jedishrfu]

I found some info on arbitrary precision computing but no specific link to fortran usage:

http://en.wikipedia.org/wiki/Arbitrary-precision_arithmetic

 

[jedishrfu]

since you're doing matrix multiplies can't you factor out a 10^300 from all elements of the matrices and then factor it back in at the end or wherever convenient?

 

[gsal]

Well, nobody seemed to have bothered asked the very first question!

fun:

Before I asked, though, will you please clarify something for me...in your posting, you say that a variable goes "larger than exp(-307)"...can you write this number in scientific notation...I am not 100 percent sure what it is.

Here are my questions:

Question 1.- Are the numbers supposed to get this large?

and then we continue:

Do you know the solution to your problem?
Have you tried your algorithm with a smaller problem you know the solution to?

before you go out on a wild goose chase looking for libraries, etc.

Although the solution suggested about reducing your number in the first place is not a bad idea...maybe, if you do not want to reduce your number right up front, you can do it once in a while and keep track how many times you have divided by 10...maybe you can divide by 10 (or 100 or whatever necessary) after every multiplication.

But, first, verify your algorithm is working correctly!

my 2 cents

 

[jedishrfu]

In fortran double precision limits are in the order of 1 e -308 to 1 e 308

Double-precision real floating-point values in IEEE T_float format ranging from 2.2250738585072013D-308 to 1.7976931348623158D308.

One other thing to check in your program is where a variable is the result of a sqrt or other math function that can't accept a negative number as an argument as fortran will set the result to nan.

 

[gsal]

Good to know.

But I just wanted fun to tell us more.

First, I thought fun's statement sounded contradictory in saying that something 'goes larger than exp(-307)'...shouldn't it be smaller than?

Also, I was not 100% sure of what fun meant with his choice of 'words'; after all, exp() is a Fortran intrinsic function and, so, when fun said exp(-307)...I did not know if he meant 2.71**-307...which, I guess, is still pretty bad.

Anyway, I am still looking forward to the answers to my questions. fun?

 

[gsal]

By the way, taking one step back...having something like 1.0e-307 ...shouldn't this simply taken as a zero?

 

[gsal]

Or did fun meant that something went so large (larger than 1.0e308) that it went negative?

 

[jedishrfu]

i looked at some f90 comments on the web and they said you get a NaN when trying to do the square root of a negative number and that he should backtrack thru his code to see where that could happen.

With respect to the 1E-308, I think there would have to be an implicit check in each calculation to readjust it to zero and that that would be unacceptable to many programmers both from efficiency and from accuracy of computation.

Instead, the burden is on the programmer to use the math rounding functions to explicitly change the number to zero when it is so close within some tolerance to zero.