# Common block mystery **URL:** \ **Category:** Help\ **Created:** [August 31, 2024, 10:25pm UTC](https://fortran-lang.discourse.group/t/common-block-mystery/8545 "2024-08-31T22:25:36Z")\ **Posts on this page:** 12\ **Page:** 1
**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\ **Post date:** [August 31, 2024, 10:25pm UTC](https://fortran-lang.discourse.group/t/common-block-mystery/8545/1 "2024-08-31T22:25:37Z")
Consider the following routine (written in 1992, from [TSPACK](https://dl.acm.org/doi/abs/10.1145/151271.151277)): ```fortran 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: ```fortran 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? ---
**Author:** ![FedericoPerini](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/federicoperini/32/1750_2.png) [@FedericoPerini](https://fortran-lang.discourse.group/u/FedericoPerini)\ **Post date:** [August 31, 2024, 11:12pm UTC](https://fortran-lang.discourse.group/t/common-block-mystery/8545/2 "2024-08-31T23:12:06Z")
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 ```Fortran RTOL>epsilon(0.d0) ``` ---
**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\ **Post date:** [August 31, 2024, 11:39pm UTC](https://fortran-lang.discourse.group/t/common-block-mystery/8545/3 "2024-08-31T23:39:37Z")
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](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580) ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [September 1, 2024, 2:29am UTC](https://fortran-lang.discourse.group/t/common-block-mystery/8545/4 "2024-09-01T02:29:00Z")
> [@FedericoPerini](#): > > 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: […] 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. ---
**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\ **Post date:** [September 1, 2024, 2:53pm UTC](https://fortran-lang.discourse.group/t/common-block-mystery/8545/5 "2024-09-01T14:53:56Z")
Here’s another odd use of this function: ```fortran IF (STORE(SIG*EXP(-SIG)+.5D0) .EQ. .5D0) GO TO 5 ``` I guess this is maybe equivalent to: ```fortran if (2.0d0*sig*exp(-sig) < epsilon(0.0d0)) ... ``` ---
**Author:** ![ilayn](https://avatars.discourse-cdn.com/v4/letter/i/c89c15/32.png) [@ilayn](https://fortran-lang.discourse.group/u/ilayn)\ **Post date:** [September 1, 2024, 6:33pm UTC](https://fortran-lang.discourse.group/t/common-block-mystery/8545/6 "2024-09-01T18:33:19Z")
As you noticed previously in SLSQP, old foxes were quite self-aware of [floating-point spacing](https://www.exploringbinary.com/the-spacing-of-binary-floating-point-numbers/) 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); ```python >>> 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`. ```python >>> 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. ---
**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\ **Post date:** [September 1, 2024, 7:14pm UTC](https://fortran-lang.discourse.group/t/common-block-mystery/8545/7 "2024-09-01T19:14:52Z")
Note: Fortran also has [SPACING](https://gcc.gnu.org/onlinedocs/gfortran/SPACING.html). Don’t think I’ve ever actually used it. ---
**Author:** ![rwmsu](https://avatars.discourse-cdn.com/v4/letter/r/48db29/32.png) [@rwmsu](https://fortran-lang.discourse.group/u/rwmsu)\ **Post date:** [September 1, 2024, 7:46pm UTC](https://fortran-lang.discourse.group/t/common-block-mystery/8545/8 "2024-09-01T19:46:18Z")
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. ```auto 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 ``` ```auto ``` ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [September 2, 2024, 2:08am UTC](https://fortran-lang.discourse.group/t/common-block-mystery/8545/10 "2024-09-02T02:08:37Z")
> [@jacobwilliams](#): > > Note: Fortran also has [SPACING](https://gcc.gnu.org/onlinedocs/gfortran/SPACING.html). 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. ---
**Author:** ![rwmsu](https://avatars.discourse-cdn.com/v4/letter/r/48db29/32.png) [@rwmsu](https://fortran-lang.discourse.group/u/rwmsu)\ **Post date:** [September 3, 2024, 12:54am UTC](https://fortran-lang.discourse.group/t/common-block-mystery/8545/11 "2024-09-03T00:54:56Z")
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 ```auto 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 \< ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [September 3, 2024, 3:17am UTC](https://fortran-lang.discourse.group/t/common-block-mystery/8545/12 "2024-09-03T03:17:36Z")
> [@rwmsu](#): > > As detailed in the link referenced by @ilayn, the SPACING (or gap between representable numbers) varies with order of magnitude. 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. ---
**Author:** ![JohnCampbell](https://avatars.discourse-cdn.com/v4/letter/j/5daacb/32.png) [@JohnCampbell](https://fortran-lang.discourse.group/u/JohnCampbell)\ **Post date:** [September 3, 2024, 3:20am UTC](https://fortran-lang.discourse.group/t/common-block-mystery/8545/13 "2024-09-03T03:20:19Z")
> [@jacobwilliams](#): > > `high precision registers.` 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.