64 DOUBLE PRECISION FUNCTION dlamch( CMACH )
75 DOUBLE PRECISION ONE, ZERO
76 parameter ( one = 1.0d+0, zero = 0.0d+0 )
80 INTEGER BETA, IMAX, IMIN, IT
81 DOUBLE PRECISION BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
82 $ rnd, sfmin, small, t
92 SAVE first, eps, sfmin, base, t, rnd, emin, rmin,
101 CALL dlamc2( beta, it, lrnd, eps, imin, rmin, imax, rmax )
106 eps = ( base**( 1-it ) ) / 2
116 IF( small.GE.sfmin )
THEN
121 sfmin = small*( one+eps )
125 IF( lsame( cmach,
'E' ) )
THEN
127 ELSE IF( lsame( cmach,
'S' ) )
THEN
129 ELSE IF( lsame( cmach,
'B' ) )
THEN
131 ELSE IF( lsame( cmach,
'P' ) )
THEN
133 ELSE IF( lsame( cmach,
'N' ) )
THEN
135 ELSE IF( lsame( cmach,
'R' ) )
THEN
137 ELSE IF( lsame( cmach,
'M' ) )
THEN
139 ELSE IF( lsame( cmach,
'U' ) )
THEN
141 ELSE IF( lsame( cmach,
'L' ) )
THEN
143 ELSE IF( lsame( cmach,
'O' ) )
THEN
206 SUBROUTINE dlamc1( BETA, T, RND, IEEE1 )
219 LOGICAL FIRST, LIEEE1, LRND
221 DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2
224 DOUBLE PRECISION DLAMC3
228 SAVE first, lieee1, lbeta, lrnd, lt
231 DATA first / .true. /
294 f = dlamc3( b / 2, -b / 100 )
301 f = dlamc3( b / 2, b / 100 )
303 IF( ( lrnd ) .AND. ( c.EQ.a ) )
312 t1 = dlamc3( b / 2, a )
313 t2 = dlamc3( b / 2, savec )
314 lieee1 = ( t1.EQ.a ) .AND. ( t2.GT.savec ) .AND. lrnd
419 SUBROUTINE dlamc2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
427 INTEGER BETA, EMAX, EMIN, T
428 DOUBLE PRECISION EPS, RMAX, RMIN
433 LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND
434 INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
436 DOUBLE PRECISION A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
437 $ sixth, small, third, two, zero
440 DOUBLE PRECISION DLAMC3
447 INTRINSIC abs, max, min
450 SAVE first, iwarn, lbeta, lemax, lemin, leps, lrmax,
454 DATA first / .true. / , iwarn / .false. /
472 CALL dlamc1( lbeta, lt, lrnd, lieee1 )
484 sixth = dlamc3( b, -half )
485 third = dlamc3( sixth, sixth )
486 b = dlamc3( third, -half )
487 b = dlamc3( b, sixth )
496 IF( ( leps.GT.b ) .AND. ( b.GT.zero ) )
THEN
498 c = dlamc3( half*leps, ( two**5 )*( leps**2 ) )
499 c = dlamc3( half, -c )
500 b = dlamc3( half, c )
501 c = dlamc3( half, -b )
502 b = dlamc3( half, c )
519 small = dlamc3( small*rbase, zero )
521 a = dlamc3( one, small )
522 CALL dlamc4( ngpmin, one, lbeta )
523 CALL dlamc4( ngnmin, -one, lbeta )
524 CALL dlamc4( gpmin, a, lbeta )
525 CALL dlamc4( gnmin, -a, lbeta )
528 IF( ( ngpmin.EQ.ngnmin ) .AND. ( gpmin.EQ.gnmin ) )
THEN
529 IF( ngpmin.EQ.gpmin )
THEN
533 ELSE IF( ( gpmin-ngpmin ).EQ.3 )
THEN
534 lemin = ngpmin - 1 + lt
539 lemin = min( ngpmin, gpmin )
544 ELSE IF( ( ngpmin.EQ.gpmin ) .AND. ( ngnmin.EQ.gnmin ) )
THEN
545 IF( abs( ngpmin-ngnmin ).EQ.1 )
THEN
546 lemin = max( ngpmin, ngnmin )
550 lemin = min( ngpmin, ngnmin )
555 ELSE IF( ( abs( ngpmin-ngnmin ).EQ.1 ) .AND.
556 $ ( gpmin.EQ.gnmin ) )
THEN
557 IF( ( gpmin-min( ngpmin, ngnmin ) ).EQ.3 )
THEN
558 lemin = max( ngpmin, ngnmin ) - 1 + lt
562 lemin = min( ngpmin, ngnmin )
568 lemin = min( ngpmin, ngnmin, gpmin, gnmin )
577 WRITE( 6, fmt = 9999 )lemin
586 ieee = ieee .OR. lieee1
593 DO 30 i = 1, 1 - lemin
594 lrmin = dlamc3( lrmin*rbase, zero )
599 CALL dlamc5( lbeta, lt, lemin, ieee, lemax, lrmax )
613 9999
FORMAT( / /
' WARNING. The value EMIN may be incorrect:-',
615 $
' If, after inspection, the value EMIN looks',
616 $
' acceptable please comment out ',
617 $ /
' the IF block as marked within the code of routine',
618 $
' DLAMC2,', /
' otherwise supply EMIN explicitly.', / )
642 DOUBLE PRECISION FUNCTION dlamc3( A, B )
649 DOUBLE PRECISION A, B
689 SUBROUTINE dlamc4( EMIN, START, BASE )
697 DOUBLE PRECISION START
703 DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
706 DOUBLE PRECISION DLAMC3
716 b1 = dlamc3( a*rbase, zero )
724 IF( ( c1.EQ.a ) .AND. ( c2.EQ.a ) .AND. ( d1.EQ.a ) .AND.
728 b1 = dlamc3( a / base, zero )
729 c1 = dlamc3( b1*base, zero )
734 b2 = dlamc3( a*rbase, zero )
735 c2 = dlamc3( b2 / rbase, zero )
796 SUBROUTINE dlamc5( BETA, P, EMIN, IEEE, EMAX, RMAX )
804 INTEGER BETA, EMAX, EMIN, P
805 DOUBLE PRECISION RMAX
810 DOUBLE PRECISION ZERO, ONE
811 parameter ( zero = 0.0d0, one = 1.0d0 )
814 INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
815 DOUBLE PRECISION OLDY, RECBAS, Y, Z
818 DOUBLE PRECISION DLAMC3
835 IF( try.LE.( -emin ) )
THEN
840 IF( lexp.EQ.-emin )
THEN
851 IF( ( uexp+emin ).GT.( -lexp-emin ) )
THEN
860 emax = expsum + emin - 1
861 nbits = 1 + exbits + p
866 IF( ( mod( nbits, 2 ).EQ.1 ) .AND. ( beta.EQ.2 ) )
THEN
911 y = dlamc3( y*beta, zero )
double precision function dlamch(CMACH)
DLAMCH
subroutine dlamc2(BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX)
DLAMC2
subroutine dlamc1(BETA, T, RND, IEEE1)
DLAMC1
subroutine dlamc5(BETA, P, EMIN, IEEE, EMAX, RMAX)
DLAMC5
double precision function dlamc3(A, B)
DLAMC3
subroutine dlamc4(EMIN, START, BASE)
DLAMC4