Difference between revisions of "Floating point math issues"
From Geos-chem
(→Safe floating-point division) |
|||
| Line 1: | Line 1: | ||
| + | == The pitfalls of floating point mathematics == | ||
| + | |||
| + | |||
| + | |||
| + | |||
| + | * [http://www.lahey.com/float.htm //The Perils of Floating Point// by Bruce M. Bush] | ||
| + | * [http://blogs.msdn.com/ericlippert/archive/2005/01/20/fun-with-floating-point-arithmetic-part-five.aspx Fun with Floating Point Arithmetic] | ||
| + | |||
== Safe floating-point division == | == Safe floating-point division == | ||
| + | |||
| + | |||
From [http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/8b367f44c419fa1d/ http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/8b367f44c419fa1d/] | From [http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/8b367f44c419fa1d/ http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/8b367f44c419fa1d/] | ||
| Line 17: | Line 27: | ||
:Thank you again, | :Thank you again, | ||
:Julio. | :Julio. | ||
| + | |||
| + | We have implemented the above algorithm into GEOS-Chem. See routine SAFE_DIV in error_mod.f: | ||
| + | |||
| + | FUNCTION SAFE_DIV( N, D, ALTV ) RESULT( Q ) | ||
| + | ! | ||
| + | !****************************************************************************** | ||
| + | ! Subroutine SAFE_DIV performs "safe division", that is to prevent overflow, | ||
| + | ! underflow, NaN, or infinity errors. An alternate value is returned if the | ||
| + | ! division cannot be performed. (bmy, 2/26/08) | ||
| + | ! | ||
| + | ! For more information, see the discussion on: http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/8b367f44c419fa1d/ | ||
| + | ! | ||
| + | ! Arguments as Input: | ||
| + | ! ============================================================================ | ||
| + | ! (1 ) N : Numerator for the division | ||
| + | ! (2 ) D : Divisor for the division | ||
| + | ! (3 ) ALTV : Alternate value to be returned if the division can't be done | ||
| + | ! | ||
| + | ! NOTES: | ||
| + | !****************************************************************************** | ||
| + | ! | ||
| + | ! Arguments | ||
| + | REAL*8, INTENT(IN) :: N, D, ALTV | ||
| + | |||
| + | ! Function value | ||
| + | REAL*8 :: Q | ||
| + | |||
| + | !================================================================== | ||
| + | ! SAFE_DIV begins here! | ||
| + | !================================================================== | ||
| + | IF ( EXPONENT(N) - EXPONENT(D) >= MAXEXPONENT(N) .or. D==0 ) THEN | ||
| + | Q = ALTV | ||
| + | ELSE | ||
| + | Q = N / D | ||
| + | ENDIF | ||
| + | |||
| + | ! Return to calling program | ||
| + | END FUNCTION SAFE_DIV | ||
Revision as of 14:32, 9 April 2008
The pitfalls of floating point mathematics
Safe floating-point division
From http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/8b367f44c419fa1d/
- Thank you to all that offered their suggestions for a "safe division" routine to prevent overflow. For those that are curious about the solution to the problem, I found useful to adopt a subroutine along the lines suggested by William Long:
if(exponent(a) - exponent (b) >= maxexponent(a) .OR. b==0)Then
q=altv
else
q=a/b
endif
- The = in the >= is to take into account the case when the fractional part of a is 1.111... and that of b is 1.
- It works very well.
- Thank you again,
- Julio.
We have implemented the above algorithm into GEOS-Chem. See routine SAFE_DIV in error_mod.f:
FUNCTION SAFE_DIV( N, D, ALTV ) RESULT( Q ) ! !****************************************************************************** ! Subroutine SAFE_DIV performs "safe division", that is to prevent overflow, ! underflow, NaN, or infinity errors. An alternate value is returned if the ! division cannot be performed. (bmy, 2/26/08) ! ! For more information, see the discussion on: http://groups.google.com/group/comp.lang.fortran/browse_thread/thread/8b367f44c419fa1d/ ! ! Arguments as Input: ! ============================================================================ ! (1 ) N : Numerator for the division ! (2 ) D : Divisor for the division ! (3 ) ALTV : Alternate value to be returned if the division can't be done ! ! NOTES: !****************************************************************************** ! ! Arguments REAL*8, INTENT(IN) :: N, D, ALTV ! Function value REAL*8 :: Q
!==================================================================
! SAFE_DIV begins here!
!==================================================================
IF ( EXPONENT(N) - EXPONENT(D) >= MAXEXPONENT(N) .or. D==0 ) THEN
Q = ALTV
ELSE
Q = N / D
ENDIF
! Return to calling program
END FUNCTION SAFE_DIV
