3
* =========== DOCUMENTATION ===========
5
* Online html documentation available at
6
* http://www.netlib.org/lapack/explore-html/
11
* DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
19
*> DLAMCH 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
59
*> \date November 2011
61
*> \ingroup auxOTHERauxiliary
63
* =====================================================================
1
64
DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
3
* -- LAPACK auxiliary routine (version 2.0) --
4
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5
* Courant Institute, Argonne National Lab, and Rice University
66
* -- LAPACK auxiliary routine (version 3.4.0) --
67
* -- LAPACK is a software package provided by Univ. of Tennessee, --
68
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8
71
* .. Scalar Arguments ..
13
* $Id: dlamch.f 19697 2010-10-29 16:57:34Z d3y133 $
18
* DLAMCH determines double precision machine parameters.
23
* CMACH (input) CHARACTER*1
24
* Specifies the value to be returned by DLAMCH:
25
* = 'E' or 'e', DLAMCH := eps
26
* = 'S' or 's , DLAMCH := sfmin
27
* = 'B' or 'b', DLAMCH := base
28
* = 'P' or 'p', DLAMCH := eps*base
29
* = 'N' or 'n', DLAMCH := t
30
* = 'R' or 'r', DLAMCH := rnd
31
* = 'M' or 'm', DLAMCH := emin
32
* = 'U' or 'u', DLAMCH := rmin
33
* = 'L' or 'l', DLAMCH := emax
34
* = 'O' or 'o', DLAMCH := rmax
38
* eps = relative machine precision
39
* sfmin = safe minimum, such that 1/sfmin does not overflow
40
* base = base of the machine
42
* t = number of (base) digits in the mantissa
43
* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
44
* emin = minimum exponent before (gradual) underflow
45
* rmin = underflow threshold - base**(emin-1)
46
* emax = largest exponent before overflow
47
* rmax = overflow threshold - (base**emax)*(1-eps)
49
75
* =====================================================================
132
************************************************************************
134
SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 )
136
* -- LAPACK auxiliary routine (version 2.0) --
137
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
138
* Courant Institute, Argonne National Lab, and Rice University
141
* .. Scalar Arguments ..
149
* DLAMC1 determines the machine parameters given by BETA, T, RND, and
155
* BETA (output) INTEGER
156
* The base of the machine.
159
* The number of ( BETA ) digits in the mantissa.
161
* RND (output) LOGICAL
162
* Specifies whether proper rounding ( RND = .TRUE. ) or
163
* chopping ( RND = .FALSE. ) occurs in addition. This may not
164
* be a reliable guide to the way in which the machine performs
167
* IEEE1 (output) LOGICAL
168
* Specifies whether rounding appears to be done in the IEEE
169
* 'round to nearest' style.
174
* The routine is based on the routine ENVRON by Malcolm and
175
* incorporates suggestions by Gentleman and Marovich. See
177
* Malcolm M. A. (1972) Algorithms to reveal properties of
178
* floating-point arithmetic. Comms. of the ACM, 15, 949-951.
180
* Gentleman W. M. and Marovich S. B. (1974) More on algorithms
181
* that reveal properties of floating point arithmetic units.
182
* Comms. of the ACM, 17, 276-277.
184
* =====================================================================
186
* .. Local Scalars ..
187
LOGICAL FIRST, LIEEE1, LRND
189
DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2
191
* .. External Functions ..
192
DOUBLE PRECISION DLAMC3
195
* .. Save statement ..
196
SAVE FIRST, LIEEE1, LBETA, LRND, LT
198
* .. Data statements ..
199
DATA FIRST / .TRUE. /
201
* .. Executable Statements ..
207
* LBETA, LIEEE1, LT and LRND are the local values of BETA,
210
* Throughout this routine we use the function DLAMC3 to ensure
211
* that relevant values are stored and not held in registers, or
212
* are not affected by optimizers.
214
* Compute a = 2.0**m with the smallest positive integer m such
222
*+ WHILE( C.EQ.ONE )LOOP
232
* Now compute b = 2.0**m with the smallest positive integer m
235
* fl( a + b ) .gt. a.
240
*+ WHILE( C.EQ.A )LOOP
249
* Now compute the base. a and c are neighbouring floating point
250
* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
251
* their difference is beta. Adding 0.25 to c is to ensure that it
252
* is truncated to beta and not ( beta - 1 ).
259
* Now determine whether rounding or chopping occurs, by adding a
260
* bit less than beta/2 and a bit more than beta/2 to a.
263
F = DLAMC3( B / 2, -B / 100 )
270
F = DLAMC3( B / 2, B / 100 )
272
IF( ( LRND ) .AND. ( C.EQ.A ) )
275
* Try and decide whether rounding is done in the IEEE 'round to
276
* nearest' style. B/2 is half a unit in the last place of the two
277
* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
278
* zero, and SAVEC is odd. Thus adding B/2 to A should not change
279
* A, but adding B/2 to SAVEC should change SAVEC.
281
T1 = DLAMC3( B / 2, A )
282
T2 = DLAMC3( B / 2, SAVEC )
283
LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
285
* Now find the mantissa, t. It should be the integer part of
286
* log to the base beta of a, however it is safer to determine t
287
* by powering. So we find t as the smallest positive integer for
290
* fl( beta**t + 1.0 ) = 1.0.
296
*+ WHILE( C.EQ.ONE )LOOP
319
************************************************************************
321
SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
323
* -- LAPACK auxiliary routine (version 2.0) --
324
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
325
* Courant Institute, Argonne National Lab, and Rice University
328
* .. Scalar Arguments ..
330
INTEGER BETA, EMAX, EMIN, T
331
DOUBLE PRECISION EPS, RMAX, RMIN
337
* DLAMC2 determines the machine parameters specified in its argument
343
* BETA (output) INTEGER
344
* The base of the machine.
347
* The number of ( BETA ) digits in the mantissa.
349
* RND (output) LOGICAL
350
* Specifies whether proper rounding ( RND = .TRUE. ) or
351
* chopping ( RND = .FALSE. ) occurs in addition. This may not
352
* be a reliable guide to the way in which the machine performs
355
* EPS (output) DOUBLE PRECISION
356
* The smallest positive number such that
358
* fl( 1.0 - EPS ) .LT. 1.0,
360
* where fl denotes the computed value.
362
* EMIN (output) INTEGER
363
* The minimum exponent before (gradual) underflow occurs.
365
* RMIN (output) DOUBLE PRECISION
366
* The smallest normalized number for the machine, given by
367
* BASE**( EMIN - 1 ), where BASE is the floating point value
370
* EMAX (output) INTEGER
371
* The maximum exponent before overflow occurs.
373
* RMAX (output) DOUBLE PRECISION
374
* The largest positive number for the machine, given by
375
* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
381
* The computation of EPS is based on a routine PARANOIA by
382
* W. Kahan of the University of California at Berkeley.
384
* =====================================================================
386
* .. Local Scalars ..
387
LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND
388
INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
390
DOUBLE PRECISION A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
391
$ SIXTH, SMALL, THIRD, TWO, ZERO
393
* .. External Functions ..
394
DOUBLE PRECISION DLAMC3
397
* .. External Subroutines ..
398
EXTERNAL DLAMC1, DLAMC4, DLAMC5
400
* .. Intrinsic Functions ..
401
INTRINSIC ABS, MAX, MIN
403
* .. Save statement ..
404
SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
407
* .. Data statements ..
408
DATA FIRST / .TRUE. / , IWARN / .FALSE. /
410
* .. Executable Statements ..
418
* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
419
* BETA, T, RND, EPS, EMIN and RMIN.
421
* Throughout this routine we use the function DLAMC3 to ensure
422
* that relevant values are stored and not held in registers, or
423
* are not affected by optimizers.
425
* DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
427
CALL DLAMC1( LBETA, LT, LRND, LIEEE1 )
435
* Try some tricks to see whether or not this is the correct EPS.
439
SIXTH = DLAMC3( B, -HALF )
440
THIRD = DLAMC3( SIXTH, SIXTH )
441
B = DLAMC3( THIRD, -HALF )
442
B = DLAMC3( B, SIXTH )
449
*+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
451
IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
453
C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
454
C = DLAMC3( HALF, -C )
455
B = DLAMC3( HALF, C )
456
C = DLAMC3( HALF, -B )
457
B = DLAMC3( HALF, C )
465
* Computation of EPS complete.
467
* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
468
* Keep dividing A by BETA until (gradual) underflow occurs. This
469
* is detected when we cannot recover the previous A.
474
SMALL = DLAMC3( SMALL*RBASE, ZERO )
476
A = DLAMC3( ONE, SMALL )
477
CALL DLAMC4( NGPMIN, ONE, LBETA )
478
CALL DLAMC4( NGNMIN, -ONE, LBETA )
479
CALL DLAMC4( GPMIN, A, LBETA )
480
CALL DLAMC4( GNMIN, -A, LBETA )
483
IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
484
IF( NGPMIN.EQ.GPMIN ) THEN
486
* ( Non twos-complement machines, no gradual underflow;
488
ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
489
LEMIN = NGPMIN - 1 + LT
491
* ( Non twos-complement machines, with gradual underflow;
492
* e.g., IEEE standard followers )
494
LEMIN = MIN( NGPMIN, GPMIN )
495
* ( A guess; no known machine )
499
ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
500
IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
501
LEMIN = MAX( NGPMIN, NGNMIN )
502
* ( Twos-complement machines, no gradual underflow;
505
LEMIN = MIN( NGPMIN, NGNMIN )
506
* ( A guess; no known machine )
510
ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
511
$ ( GPMIN.EQ.GNMIN ) ) THEN
512
IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
513
LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
514
* ( Twos-complement machines with gradual underflow;
517
LEMIN = MIN( NGPMIN, NGNMIN )
518
* ( A guess; no known machine )
523
LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
524
* ( A guess; no known machine )
528
* Comment out this if block if EMIN is ok
531
WRITE( 6, FMT = 9999 )LEMIN
535
* Assume IEEE arithmetic if we found denormalised numbers above,
536
* or if arithmetic seems to round in the IEEE style, determined
537
* in routine DLAMC1. A true IEEE machine should have both things
538
* true; however, faulty machines may have one or the other.
540
IEEE = IEEE .OR. LIEEE1
542
* Compute RMIN by successive division by BETA. We could compute
543
* RMIN as BASE**( EMIN - 1 ), but some machines underflow during
547
DO 30 I = 1, 1 - LEMIN
548
LRMIN = DLAMC3( LRMIN*RBASE, ZERO )
551
* Finally, call DLAMC5 to compute EMAX and RMAX.
553
CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
567
9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
569
$ ' If, after inspection, the value EMIN looks',
570
$ ' acceptable please comment out ',
571
$ / ' the IF block as marked within the code of routine',
572
$ ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / )
578
************************************************************************
144
************************************************************************
149
*> DLAMC3 is intended to force A and B to be stored prior to doing
150
*> the addition of A and B , for use in situations where optimizers
151
*> might hold one of these in a register.
153
*> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
154
*> \date November 2011
155
*> \ingroup auxOTHERauxiliary
159
*> A is a DOUBLE PRECISION
164
*> B is a DOUBLE PRECISION
165
*> The values A and B.
580
168
DOUBLE PRECISION FUNCTION DLAMC3( A, B )
582
* -- LAPACK auxiliary routine (version 2.0) --
583
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
584
* Courant Institute, Argonne National Lab, and Rice University
170
* -- LAPACK auxiliary routine (version 3.4.0) --
171
* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
587
174
* .. Scalar Arguments ..
588
175
DOUBLE PRECISION A, B
594
* DLAMC3 is intended to force A and B to be stored prior to doing
595
* the addition of A and B , for use in situations where optimizers
596
* might hold one of these in a register.
601
* A, B (input) DOUBLE PRECISION
602
* The values A and B.
604
177
* =====================================================================
606
179
* .. Executable Statements ..
616
189
************************************************************************
618
SUBROUTINE DLAMC4( EMIN, START, BASE )
620
* -- LAPACK auxiliary routine (version 2.0) --
621
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
622
* Courant Institute, Argonne National Lab, and Rice University
625
* .. Scalar Arguments ..
627
DOUBLE PRECISION START
633
* DLAMC4 is a service routine for DLAMC2.
639
* The minimum exponent before (gradual) underflow, computed by
640
* setting A = START and dividing by BASE until the previous A
641
* can not be recovered.
643
* START (input) DOUBLE PRECISION
644
* The starting point for determining EMIN.
646
* BASE (input) INTEGER
647
* The base of the machine.
649
* =====================================================================
651
* .. Local Scalars ..
653
DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
655
* .. External Functions ..
656
DOUBLE PRECISION DLAMC3
659
* .. Executable Statements ..
666
B1 = DLAMC3( A*RBASE, ZERO )
671
*+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
672
* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
674
IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND.
678
B1 = DLAMC3( A / BASE, ZERO )
679
C1 = DLAMC3( B1*BASE, ZERO )
684
B2 = DLAMC3( A*RBASE, ZERO )
685
C2 = DLAMC3( B2 / RBASE, ZERO )
700
************************************************************************
702
SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
704
* -- LAPACK auxiliary routine (version 2.0) --
705
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
706
* Courant Institute, Argonne National Lab, and Rice University
709
* .. Scalar Arguments ..
711
INTEGER BETA, EMAX, EMIN, P
712
DOUBLE PRECISION RMAX
718
* DLAMC5 attempts to compute RMAX, the largest machine floating-point
719
* number, without overflow. It assumes that EMAX + abs(EMIN) sum
720
* approximately to a power of 2. It will fail on machines where this
721
* assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
722
* EMAX = 28718). It will also fail if the value supplied for EMIN is
723
* too large (i.e. too close to zero), probably with overflow.
728
* BETA (input) INTEGER
729
* The base of floating-point arithmetic.
732
* The number of base BETA digits in the mantissa of a
733
* floating-point value.
735
* EMIN (input) INTEGER
736
* The minimum exponent before (gradual) underflow.
738
* IEEE (input) LOGICAL
739
* A logical flag specifying whether or not the arithmetic
740
* system is thought to comply with the IEEE standard.
742
* EMAX (output) INTEGER
743
* The largest exponent before overflow
745
* RMAX (output) DOUBLE PRECISION
746
* The largest machine floating-point number.
748
* =====================================================================
751
DOUBLE PRECISION ZERO, ONE
752
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
754
* .. Local Scalars ..
755
INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
756
DOUBLE PRECISION OLDY, RECBAS, Y, Z
758
* .. External Functions ..
759
DOUBLE PRECISION DLAMC3
762
* .. Intrinsic Functions ..
765
* .. Executable Statements ..
767
* First compute LEXP and UEXP, two powers of 2 that bound
768
* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
769
* approximately to the bound that is closest to abs(EMIN).
770
* (EMAX is the exponent of the required number RMAX).
776
IF( TRY.LE.( -EMIN ) ) THEN
781
IF( LEXP.EQ.-EMIN ) THEN
788
* Now -LEXP is less than or equal to EMIN, and -UEXP is greater
789
* than or equal to EMIN. EXBITS is the number of bits needed to
790
* store the exponent.
792
IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
798
* EXPSUM is the exponent range, approximately equal to
801
EMAX = EXPSUM + EMIN - 1
802
NBITS = 1 + EXBITS + P
804
* NBITS is the total number of bits needed to store a
805
* floating-point number.
807
IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
809
* Either there are an odd number of bits used to store a
810
* floating-point number, which is unlikely, or some bits are
811
* not used in the representation of numbers, which is possible,
812
* (e.g. Cray machines) or the mantissa has an implicit bit,
813
* (e.g. IEEE machines, Dec Vax machines), which is perhaps the
814
* most likely. We have to assume the last alternative.
815
* If this is true, then we need to reduce EMAX by one because
816
* there must be some way of representing zero in an implicit-bit
817
* system. On machines like Cray, we are reducing EMAX by one
825
* Assume we are on an IEEE machine which reserves one exponent
826
* for infinity and NaN.
831
* Now create RMAX, the largest machine number, which should
832
* be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
834
* First compute 1.0 - BETA**(-P), being careful that the
835
* result is less than 1.0 .
849
* Now multiply by BETA**EMAX to get RMAX.
852
Y = DLAMC3( Y*BETA, ZERO )