REAL FUNCTION SLAMCH( CMACH )
*
* -- LAPACK auxiliary routine (preliminary version) --
* Univ. of Tennessee, Oak Ridge National Lab, Argonne National Lab,
* Courant Institute, NAG Ltd., and Rice University
* March 26, 1990
*
* .. Scalar Arguments ..
CHARACTER*1 CMACH
* ..
*
* Purpose
* =======
*
* SLAMCH determines single precision machine parameters.
*
*-----------------------------------------------------------------------
*
* The value returned by SLAMCH is determined by the parameter CMACH as
* follows:
*
* CMACH = 'E' or 'e', SLAMCH := eps
* CMACH = 'S' or 's , SLAMCH := sfmin
* CMACH = 'B' or 'b', SLAMCH := base
* CMACH = 'N' or 'n', SLAMCH := t
* CMACH = 'R' or 'r', SLAMCH := rnd
* CMACH = 'M' or 'm', SLAMCH := emin
* CMACH = 'U' or 'u', SLAMCH := rmin
* CMACH = 'L' or 'l', SLAMCH := emax
* CMACH = 'O' or 'o', SLAMCH := rmax
*
* where
*
* eps = relative machine precision
* sfmin = safe minimum, such that 1/sfmin does not overflow
* base = base of the machine
* t = number of (base) digits in the mantissa
* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
* emin = minimum exponent before (gradual) underflow
* rmin = underflow threshold - base**(emin-1)
* emax = largest exponent before overflow
* rmax = overflow threshold - (base**emax)*(1-eps)
*
*
* .. Parameters ..
REAL ONE , ZERO
PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0)
* .. Local Scalars ..
REAL EPS, EMIN, RMIN, EMAX, RMAX,
$ BASE, RND, T, SFMIN, SMALL,
$ RMACH
INTEGER BETA, IT, IMIN, IMAX
LOGICAL FIRST, LRND
* .. External Functions ..
EXTERNAL SLAMC2
* .. Save Statement ..
SAVE FIRST, EPS, SFMIN, BASE, T, RND,
$ EMIN, RMIN, EMAX, RMAX
* .. Data Statements ..
DATA FIRST/ .TRUE. /
* ..
* .. Executable Statements ..
IF( FIRST )THEN
FIRST = .FALSE.
CALL SLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
BASE = BETA
T = IT
IF( LRND )THEN
RND = ONE
EPS = ( BASE**( 1 - IT ) )/2
ELSE
RND = ZERO
EPS = BASE**( 1 - IT )
END IF
EMIN = IMIN
EMAX = IMAX
SFMIN = RMIN
SMALL = ONE/RMAX
IF( SMALL.GE.SFMIN )THEN
*
* Use SMALL plus a bit, to avoid the possibility of rounding
* causing overflow when computing 1/sfmin.
*
SFMIN = SMALL*( ONE + EPS )
END IF
END IF
*
IF( ( CMACH.EQ.'E' ).OR.( CMACH.EQ.'e' ) )THEN
RMACH = EPS
ELSE IF( ( CMACH.EQ.'S' ).OR.( CMACH.EQ.'s' ) )THEN
RMACH = SFMIN
ELSE IF( ( CMACH.EQ.'B' ).OR.( CMACH.EQ.'b' ) )THEN
RMACH = BASE
ELSE IF( ( CMACH.EQ.'N' ).OR.( CMACH.EQ.'n' ) )THEN
RMACH = T
ELSE IF( ( CMACH.EQ.'R' ).OR.( CMACH.EQ.'r' ) )THEN
RMACH = RND
ELSE IF( ( CMACH.EQ.'M' ).OR.( CMACH.EQ.'m' ) )THEN
RMACH = EMIN
ELSE IF( ( CMACH.EQ.'U' ).OR.( CMACH.EQ.'u' ) )THEN
RMACH = RMIN
ELSE IF( ( CMACH.EQ.'L' ).OR.( CMACH.EQ.'l' ) )THEN
RMACH = EMAX
ELSE IF( ( CMACH.EQ.'O' ).OR.( CMACH.EQ.'o' ) )THEN
RMACH = RMAX
END IF
*
SLAMCH = RMACH
RETURN
*
* End of SLAMCH
*
END
*
************************************************************************
*
SUBROUTINE SLAMC1( BETA, T, RND, IEEE1 )
LOGICAL RND, IEEE1
INTEGER BETA, T
*
* SLAMC1 returns the machine parameters given by:
*
* BETA - INTEGER.
* The base of the machine.
*
* T - INTEGER.
* The number of ( BETA ) digits in the mantissa.
*
* RND - LOGICAL.
* Whether proper rounding ( RND = .TRUE. ) or chopping
* ( RND = .FALSE. ) occurs in addition. This may not be a
* reliable guide to the way in which the machine perfoms
* its arithmetic.
*
* IEEE1 - LOGICAL.
* Whether rounding appears to be done in the IEEE 'round
* to nearest' style.
*
* The routine is based on the routine ENVRON by Malcolm and
* incorporates suggestions by Gentleman and Marovich. See
*
* Malcolm M. A. (1972) Algorithms to reveal properties of
* floating-point arithmetic. Comms. of the ACM, 15, 949-951.
*
* Gentleman W. M. and Marovich S. B. (1974) More on algorithms
* that reveal properties of floating point arithmetic units.
* Comms. of the ACM, 17, 276-277.
*
*
EXTERNAL SLAMC3
LOGICAL FIRST , LIEEE1, LRND
INTEGER LBETA , LT
REAL A , B , C , F , ONE , SAVEC ,
$ T1 , T2 , QTR
REAL SLAMC3
SAVE FIRST , LIEEE1 , LBETA , LRND , LT
*
DATA FIRST/ .TRUE. /
*
IF( FIRST )THEN
FIRST = .FALSE.
ONE = 1
*
* LBETA, LIEEE1, LT and LRND are the local values of BETA,
* IEEE1, T and RND.
*
* Throughout this routine we use the function SLAMC3 to ensure
* that relevant values are stored and not held in registers, or
* are not affected by optimizers.
*
* Compute a = 2.0**m with the smallest positive integer m such
* that
*
* fl( a + 1.0 ) = a.
*
A = 1
C = 1
*
*+ WHILE( C.EQ.ONE )LOOP
10 IF ( C.EQ.ONE )THEN
A = 2*A
C = SLAMC3( A, ONE )
C = SLAMC3( C, -A )
GO TO 10
END IF
*+ END WHILE
*
* Now compute b = 2.0**m with the smallest positive integer m
* such that
*
* fl( a + b ) .gt. a.
*
B = 1
C = SLAMC3( A, B )
*
*+ WHILE( C.EQ.A )LOOP
20 IF ( C.EQ.A )THEN
B = 2*B
C = SLAMC3( A, B )
GO TO 20
END IF
*+ END WHILE
*
* Now compute the base. a and c are neighbouring floating point
* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
* their difference is beta. Adding 0.25 to c is to ensure that it
* is truncated to beta and not ( beta - 1 ).
*
QTR = ONE/4
SAVEC = C
C = SLAMC3( C, -A )
LBETA = C + QTR
*
* Now determine whether rounding or chopping occurs, by adding a
* bit less than beta/2 and a bit more than beta/2 to a.
*
B = LBETA
F = SLAMC3( B/2, -B/100 )
C = SLAMC3( F , A )
IF( C.EQ.A )THEN
LRND = .TRUE.
ELSE
LRND = .FALSE.
END IF
F = SLAMC3( B/2, B/100 )
C = SLAMC3( F , A )
IF( ( LRND ).AND.( C.EQ.A ) )
$ LRND = .FALSE.
*
* Try and decide whether rounding is done in the IEEE 'round to
* nearest' style. B/2 is half a unit in the last place of the two
* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
* zero, and SAVEC is odd. Thus adding B/2 to A should not change
* A, but adding B/2 to SAVEC should change SAVEC.
*
T1 = SLAMC3( B/2, A )
T2 = SLAMC3( B/2, SAVEC )
LIEEE1 = ( T1.EQ.A ).AND.( T2.GT.SAVEC ).AND.LRND
*
* Now find the mantissa, t. It should be the integer part of
* log to the base beta of a, however it is safer to determine t
* by powering. So we find t as the smallest positive integer for
* which
*
* fl( beta**t + 1.0 ) = 1.0.
*
LT = 0
A = 1
C = 1
*
*+ WHILE( C.EQ.ONE )LOOP
30 IF ( C.EQ.ONE )THEN
LT = LT + 1
A = A*LBETA
C = SLAMC3( A, ONE )
C = SLAMC3( C, -A )
GO TO 30
END IF
*+ END WHILE
*
END IF
*
BETA = LBETA
T = LT
RND = LRND
IEEE1 = LIEEE1
RETURN
*
* End of SLAMC1
*
END
*
************************************************************************
*
SUBROUTINE SLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
LOGICAL RND
INTEGER BETA, T, EMIN, EMAX
REAL EPS, RMIN, RMAX
*
* SLAMC2 returns the machine parameters given by:
*
* BETA - INTEGER.
* The base of the machine.
*
* T - INTEGER.
* The number of ( BETA ) digits in the mantissa.
*
* RND - LOGICAL.
* Whether proper rounding ( RND = .TRUE. ) or chopping
* ( RND = .FALSE. ) occurs in addition. This may not be a
* reliable guide to the way in which the machine perfoms
* its arithmetic.
*
* EPS - REAL.
* The smallest positive number such that
*
* fl( 1.0 - EPS ) .LT. 1.0,
*
* where fl denotes the computed value.
*
* EMIN - INTEGER.
* The minimum exponent before (gradual) underflow occurs.
*
* RMIN - REAL.
* The smallest normalized number for the machine given by
* BASE**( EMIN - 1 ), where BASE is the floating point
* value of BETA.
*
* EMAX - INTEGER.
* The maximum exponent before overflow occurs.
*
* RMAX - REAL.
* The largest positive number for the machine given by
* BASE**EMAX * ( 1 - EPS ), where BASE is the floating
* point value of BETA.
*
*
* The computation of EPS is based on a routine, PARANOIA by
* W. Kahan of the University of California at Berkeley.
*
*
EXTERNAL SLAMC3
EXTERNAL SLAMC5, SLAMC1, SLAMC4
INTRINSIC ABS , MAX , MIN
LOGICAL FIRST , IEEE , IWARN , LIEEE1, LRND
INTEGER GNMIN , GPMIN , I , LBETA , LEMAX , LEMIN ,
$ LT , NGNMIN, NGPMIN
REAL A , B , C , HALF , LEPS , LRMAX ,
$ LRMIN , SLAMC3, ONE , RBASE , SIXTH , SMALL ,
$ THIRD , TWO , ZERO
SAVE FIRST , IWARN , LBETA , LEMAX , LEMIN , LEPS ,
$ LRMAX , LRMIN , LT
DATA FIRST/ .TRUE. / , IWARN/ .FALSE. /
*
IF (FIRST) THEN
FIRST = .FALSE.
ZERO = 0
ONE = 1
TWO = 2
*
* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
* BETA, T, RND, EPS, EMIN and RMIN.
*
* Throughout this routine we use the function SLAMC3 to ensure
* that relevant values are stored and not held in registers, or
* are not affected by optimizers.
*
* SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
*
CALL SLAMC1( LBETA, LT, LRND, LIEEE1 )
*
* Start to find EPS.
*
B = LBETA
A = B**( -LT )
LEPS = A
*
* Try some tricks to see whether or not this is the correct EPS.
*
B = TWO/3
HALF = ONE/2
SIXTH = SLAMC3( B , -HALF )
THIRD = SLAMC3( SIXTH, SIXTH )
B = SLAMC3( THIRD, -HALF )
B = SLAMC3( B , SIXTH )
B = ABS ( B )
IF( B.LT.LEPS )
$ B = LEPS
*
LEPS = 1
*
*+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
20 IF ( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )THEN
LEPS = B
C = SLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
C = SLAMC3( HALF , -C )
B = SLAMC3( HALF , C )
C = SLAMC3( HALF , -B )
B = SLAMC3( HALF , C )
GO TO 20
END IF
*+ END WHILE
*
IF( A.LT.LEPS )
$ LEPS = A
*
* Computation of EPS complete.
*
* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
* Keep dividing A by BETA until (gradual) underflow occurs. This
* is detected when we cannot recover the previous A.
*
RBASE = ONE/LBETA
SMALL = ONE
DO 40, I = 1, 3
SMALL = SLAMC3( SMALL*RBASE, ZERO )
40 CONTINUE
A = SLAMC3( ONE , SMALL )
CALL SLAMC4( NGPMIN, ONE, LBETA )
CALL SLAMC4( NGNMIN, -ONE, LBETA )
CALL SLAMC4( GPMIN , A , LBETA )
CALL SLAMC4( GNMIN , -A , LBETA )
IEEE = .FALSE.
*
IF( ( NGPMIN.EQ.NGNMIN ).AND.( GPMIN.EQ.GNMIN ) )THEN
IF( NGPMIN.EQ.GPMIN )THEN
LEMIN = NGPMIN
* ( Non twos-complement machines, no gradual underflow;
* e.g., VAX )
ELSE IF( ( GPMIN - NGPMIN ).EQ.3 )THEN
LEMIN = NGPMIN - 1 + LT
IEEE = .TRUE.
* ( Non twos-complement machines, with gradual underflow;
* e.g., IEEE standard followers )
ELSE
LEMIN = MIN( NGPMIN, GPMIN )
* ( A guess; no known machine )
IWARN = .TRUE.
END IF
*
ELSE IF( ( NGPMIN.EQ.GPMIN ).AND.( NGNMIN.EQ.GNMIN ) )THEN
IF( ABS( NGPMIN - NGNMIN ).EQ.1 )THEN
LEMIN = MAX( NGPMIN, NGNMIN )
* ( Twos-complement machines, no gradual underflow;
* e.g., CYBER 205 )
ELSE
LEMIN = MIN( NGPMIN , NGNMIN )
* ( A guess; no known machine )
IWARN = .TRUE.
END IF
*
ELSE IF( ( ABS( NGPMIN - NGNMIN ).EQ.1 ).AND.
$ ( GPMIN.EQ.GNMIN ) )THEN
IF( ( GPMIN - MIN( NGPMIN, NGNMIN ) ).EQ.3 )THEN
LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
* ( Twos-complement machines with gradual underflow;
* no known machine )
ELSE
LEMIN = MIN( NGPMIN, NGNMIN )
* ( A guess; no known machine )
IWARN = .TRUE.
END IF
*
ELSE
LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
* ( A guess; no known machine )
IWARN = .TRUE.
END IF
***
* Comment out this if block if EMIN is ok
IF( IWARN )THEN
FIRST = .TRUE.
WRITE( 6, FMT=99999 )LEMIN
END IF
***
*
* Assume IEEE arithmetic if we found denormalised numbers above,
* or if arithmetic seems to round in the IEEE style, determined
* in routine SLAMC1. A true IEEE machine should have both things
* true; however, faulty machines may have one or the other.
*
IEEE = IEEE .OR. LIEEE1
*
* Compute RMIN by successive division by BETA. We could compute
* RMIN as BASE**( EMIN - 1 ), but some machines underflow during
* this computation.
*
LRMIN = 1
DO 60, I = 1, 1 - LEMIN
LRMIN = SLAMC3( LRMIN*RBASE, ZERO )
60 CONTINUE
*
* Finally, call SLAMC5 to compute EMAX and RMAX.
*
CALL SLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
END IF
*
BETA = LBETA
T = LT
RND = LRND
EPS = LEPS
EMIN = LEMIN
RMIN = LRMIN
EMAX = LEMAX
RMAX = LRMAX
*
RETURN
*
99999 FORMAT( // ' WARNING. The value EMIN may be incorrect:-',
$ ' EMIN = ', I8,
$ / ' If, after inspection, the value EMIN looks',
$ ' acceptable please comment out ',
$ / ' the IF block as marked within the code of routine',
$ ' SLAMC2,'
$ / ' otherwise supply EMIN explicitly.',
$ / )
*
* End of SLAMC2
*
END
*
************************************************************************
*
REAL FUNCTION SLAMC3( A, B )
REAL A, B
*
* SLAMC3 is intended to force A and B to be stored prior to doing
* the addition of A and B. For use in situations where optimizers
* might hold one of these in a register.
*
*
SLAMC3 = A + B
*
RETURN
*
* End of SLAMC3
*
END
*
************************************************************************
*
SUBROUTINE SLAMC4( EMIN, START, BASE)
INTEGER EMIN, BASE
REAL START
*
* Service routine for SLAMC2.
*
*
INTEGER I
REAL A , B1 , B2 , C1 , C2 , D1 ,
$ D2 , SLAMC3, ONE , RBASE , ZERO
EXTERNAL SLAMC3
*
A = START
ONE = 1
RBASE = ONE/BASE
ZERO = 0
EMIN = 1
B1 = SLAMC3( A*RBASE, ZERO )
C1 = A
C2 = A
D1 = A
D2 = A
*+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
20 IF ( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
$ ( D1.EQ.A ).AND.( D2.EQ.A ) )THEN
EMIN = EMIN - 1
A = B1
B1 = SLAMC3( A /BASE, ZERO )
C1 = SLAMC3( B1*BASE, ZERO )
D1 = ZERO
DO 40, I = 1, BASE
D1 = D1 + B1
40 CONTINUE
B2 = SLAMC3( A *RBASE, ZERO )
C2 = SLAMC3( B2/RBASE, ZERO )
D2 = ZERO
DO 60, I = 1, BASE
D2 = D2 + B2
60 CONTINUE
GO TO 20
END IF
*+ END WHILE
*
RETURN
*
* End of SLAMC4
*
END
*
************************************************************************
*
SUBROUTINE SLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
*
* Given BETA, the base of floating-point arithmetic, P, the
* number of base BETA digits in the mantissa of a floating-point
* value, EMIN, the minimum exponent, and IEEE, a logical
* flag saying whether or not the arithmetic system is thought
* to comply with the IEEE standard, this routine attempts to
* compute RMAX, the largest machine floating-point number,
* without overflow. The routine assumes that EMAX + abs(EMIN)
* sum approximately to a power of 2. It will fail on machines
* where this assumption does not hold, for example the Cyber 205
* (EMIN = -28625, EMAX = 28718). It will also fail if the value
* supplied for EMIN is too large (i.e. too close to zero),
* probably with overflow.
*
*
* .. Parameters ..
REAL ZERO , ONE
PARAMETER (ZERO = 0.0E0, ONE = 1.0E0 )
* .. Scalar Arguments ..
REAL RMAX
INTEGER BETA, EMAX, EMIN, P
LOGICAL IEEE
* .. Local Scalars ..
REAL OLDY, RECBAS, Y, Z
INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY,
* UEXP
* .. External Functions ..
REAL SLAMC3
EXTERNAL SLAMC3
* .. Intrinsic Functions ..
INTRINSIC MOD
* .. Executable Statements ..
*
* First compute LEXP and UEXP, two powers of 2 that bound
* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
* approximately to the bound that is closest to abs(EMIN).
* (EMAX is the exponent of the required number RMAX).
*
LEXP = 1
EXBITS = 1
20 CONTINUE
TRY = LEXP*2
IF(TRY.LE.( -EMIN ) )THEN
LEXP = TRY
EXBITS = EXBITS + 1
GO TO 20
END IF
IF (LEXP.EQ.-EMIN) THEN
UEXP = LEXP
ELSE
UEXP = TRY
EXBITS = EXBITS + 1
END IF
*
* Now -LEXP is less than or equal to EMIN, and -UEXP is greater
* than or equal to EMIN. EXBITS is the number of bits needed to
* store the exponent.
*
IF( ( UEXP + EMIN ).GT.( -LEXP - EMIN ) )THEN
EXPSUM = 2*LEXP
ELSE
EXPSUM = 2*UEXP
END IF
*
* EXPSUM is the exponent range, approximately equal to
* EMAX - EMIN + 1 .
*
EMAX = EXPSUM + EMIN - 1
NBITS = 1 + EXBITS + P
*
* NBITS is the total number of bits needed to store a
* floating-point number.
*
IF( ( MOD( NBITS, 2 ).EQ.1 ).AND.( BETA.EQ.2 ) )THEN
*
* Either there are an odd number of bits used to store a
* floating-point number, which is unlikely, or some bits are
* not used in the representation of numbers, which is possible,
* (e.g. Cray machines) or the mantissa has an implicit bit,
* (e.g. IEEE machines, Dec Vax machines), which is perhaps the
* most likely. We have to assume the last alternative.
* If this is true, then we need to reduce EMAX by one because
* there must be some way of representing zero in an implicit-bit
* system. On machines like Cray, we are reducing EMAX by one
* unnecessarily.
*
EMAX = EMAX - 1
END IF
*
IF (IEEE) THEN
*
* Assume we are on an IEEE machine which reserves one exponent
* for infinity and NaN.
*
EMAX = EMAX - 1
END IF
*
* Now create RMAX, the largest machine number, which should
* be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
*
* First compute 1.0 - BETA**(-P), being careful that the
* result is less than 1.0 .
*
RECBAS = ONE/BETA
Z = ONE
Y = ZERO
DO 40, I = 1, P
Z = Z*RECBAS
IF( Y.LT.ONE )
$ OLDY = Y
Y = SLAMC3( Y, Z )
40 CONTINUE
IF( Y.GE.ONE )
$ Y = OLDY
*
* Now multiply by BETA**EMAX to get RMAX.
*
DO 60, I = 1, EMAX
Y = SLAMC3( Y*BETA, ZERO )
60 CONTINUE
*
RMAX = Y
RETURN
*
* End of SLAMC5
*
END