DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. CHARACTER CMACH * .. * * Purpose * ======= * * DLAMCH determines double precision machine parameters. * * Arguments * ========= * * CMACH (input) CHARACTER*1 * Specifies the value to be returned by DLAMCH: * = 'E' or 'e', DLAMCH := eps * = 'S' or 's , DLAMCH := sfmin * = 'B' or 'b', DLAMCH := base * = 'P' or 'p', DLAMCH := eps*base * = 'N' or 'n', DLAMCH := t * = 'R' or 'r', DLAMCH := rnd * = 'M' or 'm', DLAMCH := emin * = 'U' or 'u', DLAMCH := rmin * = 'L' or 'l', DLAMCH := emax * = 'O' or 'o', DLAMCH := rmax * * where * * eps = relative machine precision * sfmin = safe minimum, such that 1/sfmin does not overflow * base = base of the machine * prec = eps*base * 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 .. DOUBLE PRECISION ONE, ZERO PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) * .. * .. Local Scalars .. LOGICAL FIRST, LRND INTEGER BETA, IMAX, IMIN, IT DOUBLE PRECISION BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN, $ RND, SFMIN, SMALL, T * .. * .. External Functions .. LOGICAL LSAME EXTERNAL LSAME * .. * .. External Subroutines .. EXTERNAL DLAMC2 * .. * .. Save statement .. SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN, $ EMAX, RMAX, PREC * .. * .. Data statements .. DATA FIRST / .TRUE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN CALL DLAMC2( 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 PREC = EPS*BASE 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( LSAME( CMACH, 'E' ) ) THEN RMACH = EPS ELSE IF( LSAME( CMACH, 'S' ) ) THEN RMACH = SFMIN ELSE IF( LSAME( CMACH, 'B' ) ) THEN RMACH = BASE ELSE IF( LSAME( CMACH, 'P' ) ) THEN RMACH = PREC ELSE IF( LSAME( CMACH, 'N' ) ) THEN RMACH = T ELSE IF( LSAME( CMACH, 'R' ) ) THEN RMACH = RND ELSE IF( LSAME( CMACH, 'M' ) ) THEN RMACH = EMIN ELSE IF( LSAME( CMACH, 'U' ) ) THEN RMACH = RMIN ELSE IF( LSAME( CMACH, 'L' ) ) THEN RMACH = EMAX ELSE IF( LSAME( CMACH, 'O' ) ) THEN RMACH = RMAX END IF * DLAMCH = RMACH FIRST = .FALSE. RETURN * * End of DLAMCH * END * ************************************************************************ * SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. LOGICAL IEEE1, RND INTEGER BETA, T * .. * * Purpose * ======= * * DLAMC1 determines the machine parameters given by BETA, T, RND, and * IEEE1. * * Arguments * ========= * * BETA (output) INTEGER * The base of the machine. * * T (output) INTEGER * The number of ( BETA ) digits in the mantissa. * * RND (output) LOGICAL * Specifies 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 performs * its arithmetic. * * IEEE1 (output) LOGICAL * Specifies whether rounding appears to be done in the IEEE * 'round to nearest' style. * * Further Details * =============== * * 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. * * ===================================================================== * * .. Local Scalars .. LOGICAL FIRST, LIEEE1, LRND INTEGER LBETA, LT DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2 * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. Save statement .. SAVE FIRST, LIEEE1, LBETA, LRND, LT * .. * .. Data statements .. DATA FIRST / .TRUE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN ONE = 1 * * LBETA, LIEEE1, LT and LRND are the local values of BETA, * IEEE1, T and RND. * * Throughout this routine we use the function DLAMC3 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 CONTINUE IF( C.EQ.ONE ) THEN A = 2*A C = DLAMC3( A, ONE ) C = DLAMC3( 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 = DLAMC3( A, B ) * *+ WHILE( C.EQ.A )LOOP 20 CONTINUE IF( C.EQ.A ) THEN B = 2*B C = DLAMC3( 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 = DLAMC3( 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 = DLAMC3( B / 2, -B / 100 ) C = DLAMC3( F, A ) IF( C.EQ.A ) THEN LRND = .TRUE. ELSE LRND = .FALSE. END IF F = DLAMC3( B / 2, B / 100 ) C = DLAMC3( 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 = DLAMC3( B / 2, A ) T2 = DLAMC3( 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 CONTINUE IF( C.EQ.ONE ) THEN LT = LT + 1 A = A*LBETA C = DLAMC3( A, ONE ) C = DLAMC3( C, -A ) GO TO 30 END IF *+ END WHILE * END IF * BETA = LBETA T = LT RND = LRND IEEE1 = LIEEE1 FIRST = .FALSE. RETURN * * End of DLAMC1 * END * ************************************************************************ * SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. LOGICAL RND INTEGER BETA, EMAX, EMIN, T DOUBLE PRECISION EPS, RMAX, RMIN * .. * * Purpose * ======= * * DLAMC2 determines the machine parameters specified in its argument * list. * * Arguments * ========= * * BETA (output) INTEGER * The base of the machine. * * T (output) INTEGER * The number of ( BETA ) digits in the mantissa. * * RND (output) LOGICAL * Specifies 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 performs * its arithmetic. * * EPS (output) DOUBLE PRECISION * The smallest positive number such that * * fl( 1.0 - EPS ) .LT. 1.0, * * where fl denotes the computed value. * * EMIN (output) INTEGER * The minimum exponent before (gradual) underflow occurs. * * RMIN (output) DOUBLE PRECISION * The smallest normalized number for the machine, given by * BASE**( EMIN - 1 ), where BASE is the floating point value * of BETA. * * EMAX (output) INTEGER * The maximum exponent before overflow occurs. * * RMAX (output) DOUBLE PRECISION * The largest positive number for the machine, given by * BASE**EMAX * ( 1 - EPS ), where BASE is the floating point * value of BETA. * * Further Details * =============== * * The computation of EPS is based on a routine PARANOIA by * W. Kahan of the University of California at Berkeley. * * ===================================================================== * * .. Local Scalars .. LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT, $ NGNMIN, NGPMIN DOUBLE PRECISION A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE, $ SIXTH, SMALL, THIRD, TWO, ZERO * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. External Subroutines .. EXTERNAL DLAMC1, DLAMC4, DLAMC5 * .. * .. Intrinsic Functions .. INTRINSIC ABS, MAX, MIN * .. * .. Save statement .. SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX, $ LRMIN, LT * .. * .. Data statements .. DATA FIRST / .TRUE. / , IWARN / .FALSE. / * .. * .. Executable Statements .. * IF( FIRST ) THEN 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 DLAMC3 to ensure * that relevant values are stored and not held in registers, or * are not affected by optimizers. * * DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. * CALL DLAMC1( 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 = DLAMC3( B, -HALF ) THIRD = DLAMC3( SIXTH, SIXTH ) B = DLAMC3( THIRD, -HALF ) B = DLAMC3( B, SIXTH ) B = ABS( B ) IF( B.LT.LEPS ) $ B = LEPS * LEPS = 1 * *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP 10 CONTINUE IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN LEPS = B C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) C = DLAMC3( HALF, -C ) B = DLAMC3( HALF, C ) C = DLAMC3( HALF, -B ) B = DLAMC3( HALF, C ) GO TO 10 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 20 I = 1, 3 SMALL = DLAMC3( SMALL*RBASE, ZERO ) 20 CONTINUE A = DLAMC3( ONE, SMALL ) CALL DLAMC4( NGPMIN, ONE, LBETA ) CALL DLAMC4( NGNMIN, -ONE, LBETA ) CALL DLAMC4( GPMIN, A, LBETA ) CALL DLAMC4( 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 FIRST = .FALSE. *** * Comment out this if block if EMIN is ok IF( IWARN ) THEN FIRST = .TRUE. WRITE( 6, FMT = 9999 )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 DLAMC1. 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 30 I = 1, 1 - LEMIN LRMIN = DLAMC3( LRMIN*RBASE, ZERO ) 30 CONTINUE * * Finally, call DLAMC5 to compute EMAX and RMAX. * CALL DLAMC5( 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 * 9999 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', $ ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / ) * * End of DLAMC2 * END * ************************************************************************ * DOUBLE PRECISION FUNCTION DLAMC3( A, B ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. DOUBLE PRECISION A, B * .. * * Purpose * ======= * * DLAMC3 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. * * Arguments * ========= * * A (input) DOUBLE PRECISION * B (input) DOUBLE PRECISION * The values A and B. * * ===================================================================== * * .. Executable Statements .. * DLAMC3 = A + B * RETURN * * End of DLAMC3 * END * ************************************************************************ * SUBROUTINE DLAMC4( EMIN, START, BASE ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. INTEGER BASE, EMIN DOUBLE PRECISION START * .. * * Purpose * ======= * * DLAMC4 is a service routine for DLAMC2. * * Arguments * ========= * * EMIN (output) INTEGER * The minimum exponent before (gradual) underflow, computed by * setting A = START and dividing by BASE until the previous A * can not be recovered. * * START (input) DOUBLE PRECISION * The starting point for determining EMIN. * * BASE (input) INTEGER * The base of the machine. * * ===================================================================== * * .. Local Scalars .. INTEGER I DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. Executable Statements .. * A = START ONE = 1 RBASE = ONE / BASE ZERO = 0 EMIN = 1 B1 = DLAMC3( 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 10 CONTINUE IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND. $ ( D2.EQ.A ) ) THEN EMIN = EMIN - 1 A = B1 B1 = DLAMC3( A / BASE, ZERO ) C1 = DLAMC3( B1*BASE, ZERO ) D1 = ZERO DO 20 I = 1, BASE D1 = D1 + B1 20 CONTINUE B2 = DLAMC3( A*RBASE, ZERO ) C2 = DLAMC3( B2 / RBASE, ZERO ) D2 = ZERO DO 30 I = 1, BASE D2 = D2 + B2 30 CONTINUE GO TO 10 END IF *+ END WHILE * RETURN * * End of DLAMC4 * END * ************************************************************************ * SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX ) * * -- LAPACK auxiliary routine (version 3.1) -- * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. * November 2006 * * .. Scalar Arguments .. LOGICAL IEEE INTEGER BETA, EMAX, EMIN, P DOUBLE PRECISION RMAX * .. * * Purpose * ======= * * DLAMC5 attempts to compute RMAX, the largest machine floating-point * number, without overflow. It 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. * * Arguments * ========= * * BETA (input) INTEGER * The base of floating-point arithmetic. * * P (input) INTEGER * The number of base BETA digits in the mantissa of a * floating-point value. * * EMIN (input) INTEGER * The minimum exponent before (gradual) underflow. * * IEEE (input) LOGICAL * A logical flag specifying whether or not the arithmetic * system is thought to comply with the IEEE standard. * * EMAX (output) INTEGER * The largest exponent before overflow * * RMAX (output) DOUBLE PRECISION * The largest machine floating-point number. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ZERO, ONE PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) * .. * .. Local Scalars .. INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP DOUBLE PRECISION OLDY, RECBAS, Y, Z * .. * .. External Functions .. DOUBLE PRECISION DLAMC3 EXTERNAL DLAMC3 * .. * .. 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 10 CONTINUE TRY = LEXP*2 IF( TRY.LE.( -EMIN ) ) THEN LEXP = TRY EXBITS = EXBITS + 1 GO TO 10 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 = BETA - ONE Y = ZERO DO 20 I = 1, P Z = Z*RECBAS IF( Y.LT.ONE ) $ OLDY = Y Y = DLAMC3( Y, Z ) 20 CONTINUE IF( Y.GE.ONE ) $ Y = OLDY * * Now multiply by BETA**EMAX to get RMAX. * DO 30 I = 1, EMAX Y = DLAMC3( Y*BETA, ZERO ) 30 CONTINUE * RMAX = Y RETURN * * End of DLAMC5 * END