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