! This file: http://ftp.aset.psu.edu/pub/ger/fortran/hdk/testdumach.f90
!
!===Two ways to compute exactly the machine unit roundoff.
! (sometimes called machine precision;see *Description below).
! Epsilon is a Fortran 90 Intrinsic Procedure.
! Dumach is a Fortran algorithm whose result is identically
! equal to Epsilon (using Double Precision in this case).
double precision dumach, x
print *, " Dumach=",Dumach()
print *, "Epsilon(x)=",Epsilon(x)
if (Dumach() .eq. Epsilon(x)) then
print *, "Dumach() = Epsilon(x)"
else
print *, "Dumach() /= Epsilon(x)"
endif
end
DOUBLE PRECISION FUNCTION DUMACH ()
!***BEGIN PROLOGUE DUMACH
!***PURPOSE Compute the unit roundoff of the machine.
!***CATEGORY R1
!***TYPE DOUBLE PRECISION (RUMACH-S, DUMACH-D)
!***KEYWORDS MACHINE CONSTANTS
!***AUTHOR Hindmarsh, Alan C., (LLNL)
!***DESCRIPTION
! *Usage:
! DOUBLE PRECISION A, DUMACH
! A = DUMACH()
!
! *Function Return Values:
! A : the unit roundoff of the machine.
!
! *Description:
! "for real x, Epsilon(x) returns a real result with the same
! type parameter as x that is almost negligible compared with
! the value one.in the model that includes x." Section 8.7.2
! in 'Fortran 90/2003 Explained' by Michael Metcalf, et al.
!
! "The unit roundoff is defined as the smallest positive machine
! number u such that 1.0 + u .ne. 1.0 "for *all* rounding modes.
! (as pointed out by James Giles at comp.lang.fortran on
! 11 July 2006).
!
! "Note that the definition in the Standard makes *NO* reference to
! the operation of addition. It is thus independent of any issues
! of rounding mode, optimizations, or other machine-dependent
! properties of the addition operation,,,
! If you see any alleged definition of epsilon() that involves
! propertiesof addition, then that wording is an approximation
! or a deduction based on specific implementation assumptions.
! It is not definitive." (as pointed out by Richard Maine,
! at comp.lang.fortran on 11 July 2006).
!
!
! This is computed by DUMACH
! in a machine-independent manner. The result matches the
! F90 Intrinsic Epsilon(x) on the IEEE model.
!
!***REFERENCES (NONE)
!***ROUTINES CALLED DUMSUM
!***REVISION HISTORY (YYYYMMDD)
! 19930216 DATE WRITTEN
! 19930818 Added SLATEC-format prologue. (FNF)
! 20030707 Added DUMSUM to force normal storage of COMP. (ACH)
!***END PROLOGUE DUMACH
!
DOUBLE PRECISION U, COMP
!***FIRST EXECUTABLE STATEMENT DUMACH
U = 1.0D0
10 U = U*0.5D0
CALL DUMSUM(1.0D0, U, COMP)
IF (COMP .NE. 1.0D0) GO TO 10
DUMACH = U*2.0D0
RETURN
!----------------------- End of Function DUMACH ------------------------
END
SUBROUTINE DUMSUM(A,B,C)
! Routine to force normal storing of A + B, for DUMACH.
DOUBLE PRECISION A, B, C
C = A + B
RETURN
END