DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) 11,1
!
! -- 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 ) 1,13
!
! -- 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 ) 1,18
!
! -- 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 ) 32
!
! -- 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 ) 4,5
!
! -- 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 ) 1,2
!
! -- 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