M_COMPARE_FLOAT_NUMBERS.F90 Subprograms

Last modified: Fri Feb 11 14:23:42 2005.

• Overview of Subprograms
• .EQUALTO.
• .GREATERTHAN.
• .LESSTHAN.
• COMPARE_FLOAT
• SOURCE:
• Listing of m_Compare_Float_Numbers.f90

Compare Float Numbers Module

This module defines elemental relational operators and functions for comparing floating point numbers. We all know that floating point numbers can be problematical beasties, particularly when we compare two real numbers. With apologies to everyone who is already familiar with the concept, the following is how not to do it,

  USE M_COMPARE_FLOAT_NUMBERS
  REAL( fp_kind ) :: x, y
  IF ( x == y ) THEN
    ! -- x equals y, so perform some operation
    ....
  END IF

Most people realise that real numbers should be compared within some tolerance value. However, a simple tolerance value will vary depending on the magnitudes of the numbers being compared. Without delving too much into the details, a better floating point comparison algorithm is,

  IF ( ABS( x < y ) < ( ULP * SPACING( MAX(ABS(x),ABS(y)) ) ) ) THEN
    ! -- x equals y, so perform some operation
    ....
  END IF

The Fortran 90/95 intrinsic function SPACING(x) returns the absolute spacing of numbers near the value of x. The ULP value, the "unit in the last place", is the smallest possible increment or decrement that can be made using a machine's floating point arithmetic. A 0.5 ULP maximum error is the best you could hope for, since this corresponds to always rounding to the nearest representable floating point number. The ULP value can be used to scale the comparison based on knowledge of the "numerical quality" of the two numbers being compared.

This floating point comparison capability is provided in the Compare_Float_Numbers module via the logical Compare_Float() function. In addition, simple relational operators .EqualTo., .GreaterThan., and .LessThan. are also available. Both single and double precision numbers can be compared, as can arrays from rank-1 to rank-7. Using the previous example, the floating point number comparison for a particular ULP scaling becomes,

  USE M_Compare_Float_Numbers
  REAL( fp_kind ) :: x, y
  IF ( Compare_Float( x, y, ULP=5 ) ) ) THEN
    ! -- x effectively equals y, so perform some operation
    ....
  END IF

For no ULP scaling (ULP=1), the relational operator can be used instead,

  USE M_Compare_Float_Numbers
  REAL( fp_kind ) :: x, y
  IF ( x .EqualTo. y ) THEN
    ! -- x effectively equals y, so perform some operation
    ....
  END IF

Subprogram Descriptions

.EQUALTO.

 NAME:
       .EqualTo.

 PURPOSE:
       Relational operator to test the equality of floating point numbers.

 CATEGORY:
       Utility

 LANGUAGE:
       Fortran-95

 CALLING SEQUENCE:

       IF ( x .EqualTo. y ) THEN
         .....
       END IF

 OPERANDS:
       x, y:        Two congruent floating point data objects to compare.

                     UNITS:      N/A
                     TYPE:       REAL( Single )   [ == default real]
                                   OR
                                 REAL( Double )
                     DIMENSION:  Scalar, or any allowed rank array.

 OPERATOR RESULT:
       (x .EqualTo. y)    The result is a logical value indicating whether
                          the operands are equal to within numerical precision

                           UNITS:      N/A
                           TYPE:       LOGICAL
                           DIMENSION:  Same as operands.

 PROCEDURE:
       The test performed is

         ABS( x - y ) < SPACING( MAX(ABS(x),ABS(y)) )

       If the result is .TRUE., the numbers are considered equal.

 CREATION HISTORY:
       Written by:     Paul van Delst, CIMSS/SSEC 30-Aug-2003
                       paul.vandelst@ssec.wisc.edu

.GREATERTHAN.

 NAME:
       .GreaterThan.

 PURPOSE:
       Relational operator to test if one operand is greater than another.

 CATEGORY:
       Utility

 LANGUAGE:
       Fortran-95

 CALLING SEQUENCE:
       IF ( x .GreaterThan. y ) THEN
         .....
       END IF

 OPERANDS:
       x, y:        Two congruent floating point data objects to compare.
                    UNITS:      N/A
                    TYPE:       REAL( Single )   [ == default real]
                                  OR
                                REAL( Double )
                    DIMENSION:  Scalar, or any allowed rank array.

 OPERATOR RESULT:
       (x .GreaterThan. y)    The result is a logical value indicating whether
                              the operand x is greater than y by more than
                              the spacing between representable floating point
                              numbers.
                              UNITS:      N/A
                              TYPE:       LOGICAL
                              DIMENSION:  Same as operands.

 PROCEDURE:
       The test performed is

         ( x - y ) >= SPACING( MAX(ABS(x),ABS(y)) )

       If the result is .TRUE., x is considered greater than y.

 CREATION HISTORY:
       Written by:     Paul van Delst, CIMSS/SSEC 30-Aug-2003
                       paul.vandelst@ssec.wisc.edu

.LESSTHAN.

 NAME:
       .LessThan.

 PURPOSE:
       Relational operator to test if one operand is less than another.

 CATEGORY:
       Utility

 LANGUAGE:
       Fortran-95

 CALLING SEQUENCE:
       IF ( x .LessThan. y ) THEN
         .....
       END IF

 OPERANDS:
       x, y:        Two congruent floating point data objects to compare.
                    UNITS:      N/A
                    TYPE:       REAL( Single )   [ == default real]
                                  OR
                                REAL( Double )
                    DIMENSION:  Scalar, or any allowed rank array.

 OPERATOR RESULT:
       (x .LessThan. y)    The result is a logical value indicating whether
                           the operand x is less than y by more than the
                           spacing between representable floating point
                           numbers.
                           UNITS:      N/A
                           TYPE:       LOGICAL
                           DIMENSION:  Same as operands.

 PROCEDURE:
       The test performed is

         ( y - x ) >= SPACING( MAX(ABS(x),ABS(y)) )

       If the result is .TRUE., x is considered less than y.

 CREATION HISTORY:
       Written by:     Paul van Delst, CIMSS/SSEC 30-Aug-2003
                       paul.vandelst@ssec.wisc.edu

COMPARE_FLOAT

 NAME:
       Compare_Float

 PURPOSE:
       Function to compare floating point scalars and arrays with adjustable
       precision tolerance.

 CATEGORY:
       Utility

 LANGUAGE:
       Fortran-95

 CALLING SEQUENCE:
       Result = Compare_Float( x, y,     &  ! Input
                               ULP = ULP )  ! Optional input

 INPUT ARGUMENTS:
       x, y:        Two congruent floating point data objects to compare.
                    UNITS:      N/A
                    TYPE:       REAL( Single )   [ == default real]
                                  OR
                                REAL( Double )
                    DIMENSION:  Scalar, or any allowed rank array.
                    ATTRIBUTES: INTENT( IN )

 OPTIONAL INPUT ARGUMENTS:
       ULP:         Unit of data precision. The acronym stands for "unit in
                    the last place," the smallest possible increment or decrement
                    that can be made using a machine's floating point arithmetic.
                    A 0.5 ulp maximum error is the best you could hope for, since
                    this corresponds to always rounding to the nearest representable
                    floating-point number. Value must be positive - if a negative
                    value is supplied, the absolute value is used.
                    If not specified, the default value is 1.
                    UNITS:      N/A
                    TYPE:       INTEGER
                    DIMENSION:  Scalar
                    ATTRIBUTES: OPTIONAL, INTENT( IN )

 OUTPUT ARGUMENTS:
       None.

 OPTIONAL OUTPUT ARGUMENTS:
       None.

 FUNCTION RESULT:
       Result:      The return value is a logical value indicating whether
                    the inputs are equal (to within the required precision)
                    .TRUE.  - if the floating point numbers are equal to
                              within the specified tolerance.
                    .FALSE. - if the floating point numbers are different.
                    UNITS:      N/A
                    TYPE:       LOGICAL
                    DIMENSION:  Scalar

 CALLS:
       None.

 SIDE EFFECTS:
       None.

 RESTRICTIONS:
       None.

 PROCEDURE:
       The test performed is

         ABS( x - y ) < ( ULP * SPACING( MAX(ABS(x),ABS(y)) ) )

       If the result is .TRUE., the numbers are considered equal.

       The intrinsic function SPACING(x) returns the absolute spacing of numbers
       near the value of x,

                      {     EXPONENT(x)-DIGITS(x)
                      {  2.0                        for x /= 0
         SPACING(x) = {
                      {
                      {  TINY(x)                    for x == 0

       The ULP optional argument scales the comparison.

       James Van Buskirk and James Giles suggested this method for floating
       point comparisons in the comp.lang.fortran newsgroup.

 CREATION HISTORY:
       Written by:     Paul van Delst, CIMSS/SSEC 01-Apr-2003
                       paul.vandelst@ssec.wisc.edu

Listing of m_Compare_Float_Numbers.f90