Will the real epsilon please stand up?

[bartoldeman]

Hi,

this connects to the post I found here when encountering optimization problems with dlamch.f
http://icl.cs.utk.edu/lapack-forum/arch ... 00247.html

The page at
http://www.netlib.org/lapack-dev/lapack ... style.html
suggests that for DLAMCH('e') the Fortran90 intrinsic EPSILON can be used (LAPACK3E does that too). However, for IEEE754 singles and doubles it gives a different value than DLAMCH.

For double precision:
EPSILON(1.0d0) = 2.220446049250313E-016
DLAMCH('e') = 1.110223024625157E-016

This follows from the different definitions for epsilon; the Fortran (and common) definition of it being the smallest positive quantity so that 1+epsilon>1.

dlamch.f, however, defines it as (DLAMC2)
* The smallest positive number such that
* fl( 1.0 - EPS ) .LT. 1.0,

in which case epsilon can be twice as big, counting bits.

Questions:
1. Does it matter much which one is chosen? I'd imagine some subtle differences in the least significant bits of the mantissa depending on who is epsilon.
2, The code in dlamch.f actually computes EPS in DLAMC2, then ignores the result and instead computes it using:
IF( LRND ) THEN
RND = ONE
EPS = ( BASE**( 1-IT ) ) / 2
ELSE
RND = ZERO
EPS = BASE**( 1-IT )
END IF
where IT is the number of digits in the mantissa, BASE=2 on any computer I have used and LRND/RND specify the rounding mode. Is there a rationale behind this?

Bart

[sven]

Just to try to clarify what LAPACK computes for EPS.

Whilst you are correct that DLAMC2 defines EPS as

* The smallest positive number such that
* fl( 1.0 - EPS ) .LT. 1.0

that is not what dlamch returns. It computes the relative machine precision, which is a bound on the relative error in converting a real number to floating point form. (Double precision floating point form in the case of dlamch.)

If we define u = (1/2)*b^(1-t), where b is the base and t is the precision, or number of digits in the significand, usually referred to as the unit roundoff, then on arithmetics that round this is the relative machine precision returned by dlamch. On arithmetics that truncate 2*u is returned.

The Fortran 90 intrinsic EPSILON, and the MATLAB built in function EPS, return 2*u. This is the distance from 1.0 to the next larger floating point value.

I agree that the leading comment for xLAMCH could be a little clearer in what it computes.

Hope this helps,

Sven Hammarling.