1 *> \brief \b DLAMCHF77 deprecated
3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
11 * DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
19 *> DLAMCHF77 determines double precision machine parameters.
27 *> Specifies the value to be returned by DLAMCH:
28 *> = 'E' or 'e', DLAMCH := eps
29 *> = 'S' or 's , DLAMCH := sfmin
30 *> = 'B' or 'b', DLAMCH := base
31 *> = 'P' or 'p', DLAMCH := eps*base
32 *> = 'N' or 'n', DLAMCH := t
33 *> = 'R' or 'r', DLAMCH := rnd
34 *> = 'M' or 'm', DLAMCH := emin
35 *> = 'U' or 'u', DLAMCH := rmin
36 *> = 'L' or 'l', DLAMCH := emax
37 *> = 'O' or 'o', DLAMCH := rmax
39 *> eps = relative machine precision
40 *> sfmin = safe minimum, such that 1/sfmin does not overflow
41 *> base = base of the machine
43 *> t = number of (base) digits in the mantissa
44 *> rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
45 *> emin = minimum exponent before (gradual) underflow
46 *> rmin = underflow threshold - base**(emin-1)
47 *> emax = largest exponent before overflow
48 *> rmax = overflow threshold - (base**emax)*(1-eps)
54 *> \author Univ. of Tennessee
55 *> \author Univ. of California Berkeley
56 *> \author Univ. of Colorado Denver
61 *> \ingroup auxOTHERauxiliary
63 * =====================================================================
64 DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
66 * -- LAPACK auxiliary routine (version 3.4.1) --
67 * -- LAPACK is a software package provided by Univ. of Tennessee, --
68 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
71 * .. Scalar Arguments ..
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
84 * .. External Functions ..
88 * .. External Subroutines ..
91 * .. Save statement ..
92 SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,
95 * .. Data statements ..
98 * .. Executable Statements ..
101 CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
106 EPS = ( BASE**( 1-IT ) ) / 2
116 IF( SMALL.GE.SFMIN ) THEN
118 * Use SMALL plus a bit, to avoid the possibility of rounding
119 * causing overflow when computing 1/sfmin.
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
155 ************************************************************************
161 *> DLAMC1 determines the machine parameters given by BETA, T, RND, and
167 *> The base of the machine.
172 *> The number of ( BETA ) digits in the mantissa.
177 *> Specifies whether proper rounding ( RND = .TRUE. ) or
178 *> chopping ( RND = .FALSE. ) occurs in addition. This may not
179 *> be a reliable guide to the way in which the machine performs
185 *> Specifies whether rounding appears to be done in the IEEE
186 *> 'round to nearest' style.
188 *> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
190 *> \ingroup auxOTHERauxiliary
192 *> \details \b Further \b Details
195 *> The routine is based on the routine ENVRON by Malcolm and
196 *> incorporates suggestions by Gentleman and Marovich. See
198 *> Malcolm M. A. (1972) Algorithms to reveal properties of
199 *> floating-point arithmetic. Comms. of the ACM, 15, 949-951.
201 *> Gentleman W. M. and Marovich S. B. (1974) More on algorithms
202 *> that reveal properties of floating point arithmetic units.
203 *> Comms. of the ACM, 17, 276-277.
206 SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 )
208 * -- LAPACK auxiliary routine (version 3.4.1) --
209 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
212 * .. Scalar Arguments ..
216 * =====================================================================
218 * .. Local Scalars ..
219 LOGICAL FIRST, LIEEE1, LRND
221 DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2
223 * .. External Functions ..
224 DOUBLE PRECISION DLAMC3
227 * .. Save statement ..
228 SAVE FIRST, LIEEE1, LBETA, LRND, LT
230 * .. Data statements ..
231 DATA FIRST / .TRUE. /
233 * .. Executable Statements ..
238 * LBETA, LIEEE1, LT and LRND are the local values of BETA,
241 * Throughout this routine we use the function DLAMC3 to ensure
242 * that relevant values are stored and not held in registers, or
243 * are not affected by optimizers.
245 * Compute a = 2.0**m with the smallest positive integer m such
253 *+ WHILE( C.EQ.ONE )LOOP
263 * Now compute b = 2.0**m with the smallest positive integer m
266 * fl( a + b ) .gt. a.
271 *+ WHILE( C.EQ.A )LOOP
280 * Now compute the base. a and c are neighbouring floating point
281 * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
282 * their difference is beta. Adding 0.25 to c is to ensure that it
283 * is truncated to beta and not ( beta - 1 ).
290 * Now determine whether rounding or chopping occurs, by adding a
291 * bit less than beta/2 and a bit more than beta/2 to a.
294 F = DLAMC3( B / 2, -B / 100 )
301 F = DLAMC3( B / 2, B / 100 )
303 IF( ( LRND ) .AND. ( C.EQ.A ) )
306 * Try and decide whether rounding is done in the IEEE 'round to
307 * nearest' style. B/2 is half a unit in the last place of the two
308 * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
309 * zero, and SAVEC is odd. Thus adding B/2 to A should not change
310 * A, but adding B/2 to SAVEC should change SAVEC.
312 T1 = DLAMC3( B / 2, A )
313 T2 = DLAMC3( B / 2, SAVEC )
314 LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
316 * Now find the mantissa, t. It should be the integer part of
317 * log to the base beta of a, however it is safer to determine t
318 * by powering. So we find t as the smallest positive integer for
321 * fl( beta**t + 1.0 ) = 1.0.
327 *+ WHILE( C.EQ.ONE )LOOP
351 ************************************************************************
357 *> DLAMC2 determines the machine parameters specified in its argument
360 *> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
362 *> \ingroup auxOTHERauxiliary
366 *> The base of the machine.
371 *> The number of ( BETA ) digits in the mantissa.
376 *> Specifies whether proper rounding ( RND = .TRUE. ) or
377 *> chopping ( RND = .FALSE. ) occurs in addition. This may not
378 *> be a reliable guide to the way in which the machine performs
384 *> The smallest positive number such that
385 *> fl( 1.0 - EPS ) .LT. 1.0,
386 *> where fl denotes the computed value.
391 *> The minimum exponent before (gradual) underflow occurs.
396 *> The smallest normalized number for the machine, given by
397 *> BASE**( EMIN - 1 ), where BASE is the floating point value
403 *> The maximum exponent before overflow occurs.
408 *> The largest positive number for the machine, given by
409 *> BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
413 *> \details \b Further \b Details
416 *> The computation of EPS is based on a routine PARANOIA by
417 *> W. Kahan of the University of California at Berkeley.
419 SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
421 * -- LAPACK auxiliary routine (version 3.4.1) --
422 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
425 * .. Scalar Arguments ..
427 INTEGER BETA, EMAX, EMIN, T
428 DOUBLE PRECISION EPS, RMAX, RMIN
430 * =====================================================================
432 * .. Local Scalars ..
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
439 * .. External Functions ..
440 DOUBLE PRECISION DLAMC3
443 * .. External Subroutines ..
444 EXTERNAL DLAMC1, DLAMC4, DLAMC5
446 * .. Intrinsic Functions ..
447 INTRINSIC ABS, MAX, MIN
449 * .. Save statement ..
450 SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
453 * .. Data statements ..
454 DATA FIRST / .TRUE. / , IWARN / .FALSE. /
456 * .. Executable Statements ..
463 * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
464 * BETA, T, RND, EPS, EMIN and RMIN.
466 * Throughout this routine we use the function DLAMC3 to ensure
467 * that relevant values are stored and not held in registers, or
468 * are not affected by optimizers.
470 * DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
472 CALL DLAMC1( LBETA, LT, LRND, LIEEE1 )
480 * Try some tricks to see whether or not this is the correct EPS.
484 SIXTH = DLAMC3( B, -HALF )
485 THIRD = DLAMC3( SIXTH, SIXTH )
486 B = DLAMC3( THIRD, -HALF )
487 B = DLAMC3( B, SIXTH )
494 *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
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 )
510 * Computation of EPS complete.
512 * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
513 * Keep dividing A by BETA until (gradual) underflow occurs. This
514 * is detected when we cannot recover the previous A.
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
531 * ( Non twos-complement machines, no gradual underflow;
533 ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
534 LEMIN = NGPMIN - 1 + LT
536 * ( Non twos-complement machines, with gradual underflow;
537 * e.g., IEEE standard followers )
539 LEMIN = MIN( NGPMIN, GPMIN )
540 * ( A guess; no known machine )
544 ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
545 IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
546 LEMIN = MAX( NGPMIN, NGNMIN )
547 * ( Twos-complement machines, no gradual underflow;
550 LEMIN = MIN( NGPMIN, NGNMIN )
551 * ( A guess; no known machine )
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
559 * ( Twos-complement machines with gradual underflow;
562 LEMIN = MIN( NGPMIN, NGNMIN )
563 * ( A guess; no known machine )
568 LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
569 * ( A guess; no known machine )
574 * Comment out this if block if EMIN is ok
577 WRITE( 6, FMT = 9999 )LEMIN
581 * Assume IEEE arithmetic if we found denormalised numbers above,
582 * or if arithmetic seems to round in the IEEE style, determined
583 * in routine DLAMC1. A true IEEE machine should have both things
584 * true; however, faulty machines may have one or the other.
586 IEEE = IEEE .OR. LIEEE1
588 * Compute RMIN by successive division by BETA. We could compute
589 * RMIN as BASE**( EMIN - 1 ), but some machines underflow during
593 DO 30 I = 1, 1 - LEMIN
594 LRMIN = DLAMC3( LRMIN*RBASE, ZERO )
597 * Finally, call DLAMC5 to compute EMAX and RMAX.
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.', / )
624 ************************************************************************
630 *> DLAMC3 is intended to force A and B to be stored prior to doing
631 *> the addition of A and B , for use in situations where optimizers
632 *> might hold one of these in a register.
639 *> The values A and B.
642 DOUBLE PRECISION FUNCTION DLAMC3( A, B )
644 * -- LAPACK auxiliary routine (version 3.4.1) --
645 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
648 * .. Scalar Arguments ..
649 DOUBLE PRECISION A, B
651 * =====================================================================
653 * .. Executable Statements ..
663 ************************************************************************
669 *> DLAMC4 is a service routine for DLAMC2.
674 *> The minimum exponent before (gradual) underflow, computed by
675 *> setting A = START and dividing by BASE until the previous A
676 *> can not be recovered.
681 *> The starting point for determining EMIN.
686 *> The base of the machine.
689 SUBROUTINE DLAMC4( EMIN, START, BASE )
691 * -- LAPACK auxiliary routine (version 3.4.1) --
692 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
695 * .. Scalar Arguments ..
697 DOUBLE PRECISION START
699 * =====================================================================
701 * .. Local Scalars ..
703 DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
705 * .. External Functions ..
706 DOUBLE PRECISION DLAMC3
709 * .. Executable Statements ..
716 B1 = DLAMC3( A*RBASE, ZERO )
721 *+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
722 * $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
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 )
750 ************************************************************************
756 *> DLAMC5 attempts to compute RMAX, the largest machine floating-point
757 *> number, without overflow. It assumes that EMAX + abs(EMIN) sum
758 *> approximately to a power of 2. It will fail on machines where this
759 *> assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
760 *> EMAX = 28718). It will also fail if the value supplied for EMIN is
761 *> too large (i.e. too close to zero), probably with overflow.
766 *> The base of floating-point arithmetic.
771 *> The number of base BETA digits in the mantissa of a
772 *> floating-point value.
777 *> The minimum exponent before (gradual) underflow.
782 *> A logical flag specifying whether or not the arithmetic
783 *> system is thought to comply with the IEEE standard.
788 *> The largest exponent before overflow
793 *> The largest machine floating-point number.
796 SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
798 * -- LAPACK auxiliary routine (version 3.4.1) --
799 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
802 * .. Scalar Arguments ..
804 INTEGER BETA, EMAX, EMIN, P
805 DOUBLE PRECISION RMAX
807 * =====================================================================
810 DOUBLE PRECISION ZERO, ONE
811 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
813 * .. Local Scalars ..
814 INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
815 DOUBLE PRECISION OLDY, RECBAS, Y, Z
817 * .. External Functions ..
818 DOUBLE PRECISION DLAMC3
821 * .. Intrinsic Functions ..
824 * .. Executable Statements ..
826 * First compute LEXP and UEXP, two powers of 2 that bound
827 * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
828 * approximately to the bound that is closest to abs(EMIN).
829 * (EMAX is the exponent of the required number RMAX).
835 IF( TRY.LE.( -EMIN ) ) THEN
840 IF( LEXP.EQ.-EMIN ) THEN
847 * Now -LEXP is less than or equal to EMIN, and -UEXP is greater
848 * than or equal to EMIN. EXBITS is the number of bits needed to
849 * store the exponent.
851 IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
857 * EXPSUM is the exponent range, approximately equal to
860 EMAX = EXPSUM + EMIN - 1
861 NBITS = 1 + EXBITS + P
863 * NBITS is the total number of bits needed to store a
864 * floating-point number.
866 IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
868 * Either there are an odd number of bits used to store a
869 * floating-point number, which is unlikely, or some bits are
870 * not used in the representation of numbers, which is possible,
871 * (e.g. Cray machines) or the mantissa has an implicit bit,
872 * (e.g. IEEE machines, Dec Vax machines), which is perhaps the
873 * most likely. We have to assume the last alternative.
874 * If this is true, then we need to reduce EMAX by one because
875 * there must be some way of representing zero in an implicit-bit
876 * system. On machines like Cray, we are reducing EMAX by one
884 * Assume we are on an IEEE machine which reserves one exponent
885 * for infinity and NaN.
890 * Now create RMAX, the largest machine number, which should
891 * be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
893 * First compute 1.0 - BETA**(-P), being careful that the
894 * result is less than 1.0 .
908 * Now multiply by BETA**EMAX to get RMAX.
911 Y = DLAMC3( Y*BETA, ZERO )