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 FUNCTION DLAMCH ! !*********************************************************************** ! 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 DLAMC1 ! !*********************************************************************** ! 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 SUBROUTINE DLAMC2 ! !*********************************************************************** ! 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 FUNCTION DLAMC3 ! !*********************************************************************** ! 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 DLAMC4 ! !*********************************************************************** ! 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 SUBROUTINE DLAMC5