Common block mystery

Help
1

Consider the following routine (written in 1992, from TSPACK):

C   This function forces its argument X to be stored in a
C memory location, thus providing a means of determining
C floating point number characteristics (such as the machine
C precision) when it is necessary to avoid computation in
C high precision registers.
C
C On input:
C
C       X = Value to be stored.
C
C X is not altered by this function.
C
C On output:
C
C       STORE = Value of X after it has been stored and
C               possibly truncated or rounded to the single
C               precision word length.
      DOUBLE PRECISION FUNCTION STORE (X)
      DOUBLE PRECISION X
      DOUBLE PRECISION Y
      COMMON/STCOM/Y
      Y = X
      STORE = Y
      RETURN
      END

It’s used in various places like this:

IF (STORE(RTOL+1.D0) .GT. 1.D0) GO TO 1

It seems to be some kind of numerical trickery or epsilon test. I don’t understand why stuffing the input into a common block makes any difference. What is going on here?

1 Like
2

Renka put this also in TRIPACK. In my understanding, saving the variable to a common is a trick to ensure that the arithmetic expression does not get simplified out by the compiler: it could decide to just check for RTOL>0, that is not what the author is trying to do here.

In modern Fortran, you should just check for

RTOL>epsilon(0.d0)
2 Likes
3

Yep, sounds reasonable! See this earlier post for another old school way to trick the compilers to do something like that: Looking at some old code

2 Likes
4

This was one way to force the intermediate to be represented as a 64-bit real rather than an 80-bit real. There were also compiler options that could be used to accomplish this same thing, but this approach would work regardless of the compiler options.

3 Likes
5

Here’s another odd use of this function:

IF (STORE(SIG*EXP(-SIG)+.5D0) .EQ. .5D0) GO TO 5

I guess this is maybe equivalent to:

if (2.0d0*sig*exp(-sig) < epsilon(0.0d0)) ...
1 Like
6

As you noticed previously in SLSQP, old foxes were quite self-aware of floating-point spacing

Hence without the necessary tooling back then, they adopted some -isms that you see a lot in old code. Another common one is res = 0.5D0 + (0.5D0 - x) from the 70s.

The case they are trying to avoid is the following (sorry for the NumPy code);

>>> import numpy as np
>>> zzz = np.float64(1234567891011.123456789)

>>> zzz + 1e-6 == zzz
True

>>> zzz + 1e-5 == zzz
True

>>> zzz + 1e-1 == zzz
False

In other words, when is the added value cannot be distinguished in that particular dtype for that magnitude? The switching point from True to False is the zzz*del and del is the distance to the next number in that particular dtype for that magnitude (between two consecutive powers of 2). NumPy also has a convenience function with the same name numpy.spacing.

>>> s = np.spacing(zzz)
>>> s
np.float64(0.000244140625)

>>> s / zzz
np.float64(1.9775390788759812e-16)

>>> np.log2(s / zzz)
np.float64(-52.1671433135915)

The test number is always a power of two to get the relevant spacing in that particular magnitude hence you would see 0.5, 1.0 or other powers of two depending on what they are testing.

2 Likes
7

Note: Fortran also has SPACING. Don’t think I’ve ever actually used it.

1 Like
8

I have. The following is an almostEqual function patterned after a post on Intel Fortran Forum that uses SPACING. You could also use EPSILON (which is what I originally had). I use this in most places where I need to test for equlity (to at least a tolerance) of two floating point numbers. This version is for REAL32 but I also have one for REAL64. The SPACING_FACTOR is set to a constant 5 but thats just what the original post had. You could use 2 if you want a tighter tolerance.



  Elemental Function almostEqualR32(a, b) Result(aeqb)

! Function to test if two WP floating point numbers 
! almost equal to close to machine precision 

! Code taken from post on Intel Fortran Forum 11/15/2016

! Argument variables

    Real(REAL32), Intent(IN) :: a, b
    Logical              :: aeqb

    aeqb = ABS(a-b) < SPACING_FACTOR32*SPACING(MAX(ABS(a), ABS(b)))

  End Function almostEqualR32
2 Likes
10

I think fortran has had the SPACING() intrinsic since f90. If so, that is 15 years or so before numpy existed. I don’t know where the idea originated before that.

1 Like
11

Here is a little program that demonstrates the difference between SPACING and EPSILON. As detailed in the link referenced by @ilayn, the SPACING (or gap between representable numbers) varies with order of magnitude. EPSILON represents the gap in numbers between 1 and 2 and is therefore constant irrespective of the number you give it. To me the results from this little exercise demonstrates why you have to be careful if you want to use default 32 bit precision in floating point intensive calculations. If you have data that varies several orders of magnitude (ie powers of 10), then you should consider normalizing your data back into a range of probably 0 to 100.

This program

Program spacing_v_eps

  USE ISO_FORTRAN_ENV, SP=>REAL32, DP=>REAL64

  Implicit NONE

  Real(SP) :: a4
  Real(DP) :: a8

  Integer :: i

  Print *,''
  Print *,' SPACING - 10**i'
  Print *,''
  Do i = 0,10
    a4 = 10.0_SP**i
    a8 = 10.0_DP**i
    Write(*,'(" i = ",i2," r4 = ", ES0.7, " r8 = ", ES0.16)') i, SPACING(a4), &
              SPACING(a8)
  End Do
  Print *,''
  Print *,' EPSILON - 10**i'
  Print *,''
  Do i = 0,10
    a4 = 10.0_SP**i
    a8 = 10.0_DP**i
    Write(*,'(" i = ",i2," r4 = ", ES0.7, " r8 = ", ES0.16)') i, EPSILON(a4), &
              EPSILON(a8)
  End Do

End Program spacing_v_eps

gives these results when compiled with ifx 2024.1

SPACING - 10**i

i = 0 r4 = 1.1920929-07 r8 = 2.2204460492503131-16
i = 1 r4 = 9.5367432-07 r8 = 1.7763568394002505-15
i = 2 r4 = 7.6293945-06 r8 = 1.4210854715202004-14
i = 3 r4 = 6.1035156-05 r8 = 1.1368683772161603-13
i = 4 r4 = 9.7656250-04 r8 = 1.8189894035458565-12
i = 5 r4 = 7.8125000-03 r8 = 1.4551915228366852-11
i = 6 r4 = 6.2500000-02 r8 = 1.1641532182693481-10
i = 7 r4 = 1.0000000+00 r8 = 1.8626451492309570-09
i = 8 r4 = 8.0000000+00 r8 = 1.4901161193847656-08
i = 9 r4 = 6.4000000+01 r8 = 1.1920928955078125-07
i = 10 r4 = 1.0240000+03 r8 = 1.9073486328125000-06

EPSILON - 10**i

i = 0 r4 = 1.1920929-07 r8 = 2.2204460492503131-16
i = 1 r4 = 1.1920929-07 r8 = 2.2204460492503131-16
i = 2 r4 = 1.1920929-07 r8 = 2.2204460492503131-16
i = 3 r4 = 1.1920929-07 r8 = 2.2204460492503131-16
i = 4 r4 = 1.1920929-07 r8 = 2.2204460492503131-16
i = 5 r4 = 1.1920929-07 r8 = 2.2204460492503131-16
i = 6 r4 = 1.1920929-07 r8 = 2.2204460492503131-16
i = 7 r4 = 1.1920929-07 r8 = 2.2204460492503131-16
i = 8 r4 = 1.1920929-07 r8 = 2.2204460492503131-16
i = 9 r4 = 1.1920929-07 r8 = 2.2204460492503131-16
i = 10 r4 = 1.1920929-07 r8 = 2.2204460492503131-16 <

12

All of the floating point numbers with the same binary exponent will have the same SPACING() value, so you can think of the spacing as epsilon*2**exponent while ignoring the mantissa. In contrast, most numerical analysis is based on the idea that x*(1+epsilon) is the floating point result with x being the exact mathematical value. This obviously does depend on the mantissa bits. These two values are slightly different perspectives of floating point error.

13

If 1992, this may be a reference to 80-bit registers.
I recall Lahey Fortran provided compiler options to limit real*8 calculations to 64-bit precision, rather than the default of re-use of the values stored in the 80-bit registers.

1 Like