1
REAL FUNCTION SLAMCH( 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
8
* .. Scalar Arguments ..
15
* SLAMCH determines single precision machine parameters.
20
* CMACH (input) CHARACTER*1
21
* Specifies the value to be returned by SLAMCH:
22
* = 'E' or 'e', SLAMCH := eps
23
* = 'S' or 's , SLAMCH := sfmin
24
* = 'B' or 'b', SLAMCH := base
25
* = 'P' or 'p', SLAMCH := eps*base
26
* = 'N' or 'n', SLAMCH := t
27
* = 'R' or 'r', SLAMCH := rnd
28
* = 'M' or 'm', SLAMCH := emin
29
* = 'U' or 'u', SLAMCH := rmin
30
* = 'L' or 'l', SLAMCH := emax
31
* = 'O' or 'o', SLAMCH := rmax
35
* eps = relative machine precision
36
* sfmin = safe minimum, such that 1/sfmin does not overflow
37
* base = base of the machine
39
* t = number of (base) digits in the mantissa
40
* rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
41
* emin = minimum exponent before (gradual) underflow
42
* rmin = underflow threshold - base**(emin-1)
43
* emax = largest exponent before overflow
44
* rmax = overflow threshold - (base**emax)*(1-eps)
46
* =====================================================================
50
PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
54
INTEGER BETA, IMAX, IMIN, IT
55
REAL BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
56
$ RND, SFMIN, SMALL, T
58
* .. External Functions ..
62
* .. External Subroutines ..
65
* .. Save statement ..
66
SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,
69
* .. Data statements ..
72
* .. Executable Statements ..
76
CALL SLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
81
EPS = ( BASE**( 1-IT ) ) / 2
91
IF( SMALL.GE.SFMIN ) THEN
93
* Use SMALL plus a bit, to avoid the possibility of rounding
94
* causing overflow when computing 1/sfmin.
96
SFMIN = SMALL*( ONE+EPS )
100
IF( LSAME( CMACH, 'E' ) ) THEN
102
ELSE IF( LSAME( CMACH, 'S' ) ) THEN
104
ELSE IF( LSAME( CMACH, 'B' ) ) THEN
106
ELSE IF( LSAME( CMACH, 'P' ) ) THEN
108
ELSE IF( LSAME( CMACH, 'N' ) ) THEN
110
ELSE IF( LSAME( CMACH, 'R' ) ) THEN
112
ELSE IF( LSAME( CMACH, 'M' ) ) THEN
114
ELSE IF( LSAME( CMACH, 'U' ) ) THEN
116
ELSE IF( LSAME( CMACH, 'L' ) ) THEN
118
ELSE IF( LSAME( CMACH, 'O' ) ) THEN
129
************************************************************************
131
SUBROUTINE SLAMC1( BETA, T, RND, IEEE1 )
133
* -- LAPACK auxiliary routine (version 2.0) --
134
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
135
* Courant Institute, Argonne National Lab, and Rice University
138
* .. Scalar Arguments ..
146
* SLAMC1 determines the machine parameters given by BETA, T, RND, and
152
* BETA (output) INTEGER
153
* The base of the machine.
156
* The number of ( BETA ) digits in the mantissa.
158
* RND (output) LOGICAL
159
* Specifies whether proper rounding ( RND = .TRUE. ) or
160
* chopping ( RND = .FALSE. ) occurs in addition. This may not
161
* be a reliable guide to the way in which the machine performs
164
* IEEE1 (output) LOGICAL
165
* Specifies whether rounding appears to be done in the IEEE
166
* 'round to nearest' style.
171
* The routine is based on the routine ENVRON by Malcolm and
172
* incorporates suggestions by Gentleman and Marovich. See
174
* Malcolm M. A. (1972) Algorithms to reveal properties of
175
* floating-point arithmetic. Comms. of the ACM, 15, 949-951.
177
* Gentleman W. M. and Marovich S. B. (1974) More on algorithms
178
* that reveal properties of floating point arithmetic units.
179
* Comms. of the ACM, 17, 276-277.
181
* =====================================================================
183
* .. Local Scalars ..
184
LOGICAL FIRST, LIEEE1, LRND
186
REAL A, B, C, F, ONE, QTR, SAVEC, T1, T2
188
* .. External Functions ..
192
* .. Save statement ..
193
SAVE FIRST, LIEEE1, LBETA, LRND, LT
195
* .. Data statements ..
196
DATA FIRST / .TRUE. /
198
* .. Executable Statements ..
204
* LBETA, LIEEE1, LT and LRND are the local values of BETA,
207
* Throughout this routine we use the function SLAMC3 to ensure
208
* that relevant values are stored and not held in registers, or
209
* are not affected by optimizers.
211
* Compute a = 2.0**m with the smallest positive integer m such
219
*+ WHILE( C.EQ.ONE )LOOP
229
* Now compute b = 2.0**m with the smallest positive integer m
232
* fl( a + b ) .gt. a.
237
*+ WHILE( C.EQ.A )LOOP
246
* Now compute the base. a and c are neighbouring floating point
247
* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
248
* their difference is beta. Adding 0.25 to c is to ensure that it
249
* is truncated to beta and not ( beta - 1 ).
256
* Now determine whether rounding or chopping occurs, by adding a
257
* bit less than beta/2 and a bit more than beta/2 to a.
260
F = SLAMC3( B / 2, -B / 100 )
267
F = SLAMC3( B / 2, B / 100 )
269
IF( ( LRND ) .AND. ( C.EQ.A ) )
272
* Try and decide whether rounding is done in the IEEE 'round to
273
* nearest' style. B/2 is half a unit in the last place of the two
274
* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
275
* zero, and SAVEC is odd. Thus adding B/2 to A should not change
276
* A, but adding B/2 to SAVEC should change SAVEC.
278
T1 = SLAMC3( B / 2, A )
279
T2 = SLAMC3( B / 2, SAVEC )
280
LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
282
* Now find the mantissa, t. It should be the integer part of
283
* log to the base beta of a, however it is safer to determine t
284
* by powering. So we find t as the smallest positive integer for
287
* fl( beta**t + 1.0 ) = 1.0.
293
*+ WHILE( C.EQ.ONE )LOOP
316
************************************************************************
318
SUBROUTINE SLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
320
* -- LAPACK auxiliary routine (version 2.0) --
321
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
322
* Courant Institute, Argonne National Lab, and Rice University
325
* .. Scalar Arguments ..
327
INTEGER BETA, EMAX, EMIN, T
334
* SLAMC2 determines the machine parameters specified in its argument
340
* BETA (output) INTEGER
341
* The base of the machine.
344
* The number of ( BETA ) digits in the mantissa.
346
* RND (output) LOGICAL
347
* Specifies whether proper rounding ( RND = .TRUE. ) or
348
* chopping ( RND = .FALSE. ) occurs in addition. This may not
349
* be a reliable guide to the way in which the machine performs
353
* The smallest positive number such that
355
* fl( 1.0 - EPS ) .LT. 1.0,
357
* where fl denotes the computed value.
359
* EMIN (output) INTEGER
360
* The minimum exponent before (gradual) underflow occurs.
363
* The smallest normalized number for the machine, given by
364
* BASE**( EMIN - 1 ), where BASE is the floating point value
367
* EMAX (output) INTEGER
368
* The maximum exponent before overflow occurs.
371
* The largest positive number for the machine, given by
372
* BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
378
* The computation of EPS is based on a routine PARANOIA by
379
* W. Kahan of the University of California at Berkeley.
381
* =====================================================================
383
* .. Local Scalars ..
384
LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND
385
INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
387
REAL A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
388
$ SIXTH, SMALL, THIRD, TWO, ZERO
390
* .. External Functions ..
394
* .. External Subroutines ..
395
EXTERNAL SLAMC1, SLAMC4, SLAMC5
397
* .. Intrinsic Functions ..
398
INTRINSIC ABS, MAX, MIN
400
* .. Save statement ..
401
SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
404
* .. Data statements ..
405
DATA FIRST / .TRUE. / , IWARN / .FALSE. /
407
* .. Executable Statements ..
415
* LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
416
* BETA, T, RND, EPS, EMIN and RMIN.
418
* Throughout this routine we use the function SLAMC3 to ensure
419
* that relevant values are stored and not held in registers, or
420
* are not affected by optimizers.
422
* SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
424
CALL SLAMC1( LBETA, LT, LRND, LIEEE1 )
432
* Try some tricks to see whether or not this is the correct EPS.
436
SIXTH = SLAMC3( B, -HALF )
437
THIRD = SLAMC3( SIXTH, SIXTH )
438
B = SLAMC3( THIRD, -HALF )
439
B = SLAMC3( B, SIXTH )
446
*+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
448
IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
450
C = SLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
451
C = SLAMC3( HALF, -C )
452
B = SLAMC3( HALF, C )
453
C = SLAMC3( HALF, -B )
454
B = SLAMC3( HALF, C )
462
* Computation of EPS complete.
464
* Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
465
* Keep dividing A by BETA until (gradual) underflow occurs. This
466
* is detected when we cannot recover the previous A.
471
SMALL = SLAMC3( SMALL*RBASE, ZERO )
473
A = SLAMC3( ONE, SMALL )
474
CALL SLAMC4( NGPMIN, ONE, LBETA )
475
CALL SLAMC4( NGNMIN, -ONE, LBETA )
476
CALL SLAMC4( GPMIN, A, LBETA )
477
CALL SLAMC4( GNMIN, -A, LBETA )
480
IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
481
IF( NGPMIN.EQ.GPMIN ) THEN
483
* ( Non twos-complement machines, no gradual underflow;
485
ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
486
LEMIN = NGPMIN - 1 + LT
488
* ( Non twos-complement machines, with gradual underflow;
489
* e.g., IEEE standard followers )
491
LEMIN = MIN( NGPMIN, GPMIN )
492
* ( A guess; no known machine )
496
ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
497
IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
498
LEMIN = MAX( NGPMIN, NGNMIN )
499
* ( Twos-complement machines, no gradual underflow;
502
LEMIN = MIN( NGPMIN, NGNMIN )
503
* ( A guess; no known machine )
507
ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
508
$ ( GPMIN.EQ.GNMIN ) ) THEN
509
IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
510
LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
511
* ( Twos-complement machines with gradual underflow;
514
LEMIN = MIN( NGPMIN, NGNMIN )
515
* ( A guess; no known machine )
520
LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
521
* ( A guess; no known machine )
525
* Comment out this if block if EMIN is ok
528
WRITE( 6, FMT = 9999 )LEMIN
532
* Assume IEEE arithmetic if we found denormalised numbers above,
533
* or if arithmetic seems to round in the IEEE style, determined
534
* in routine SLAMC1. A true IEEE machine should have both things
535
* true; however, faulty machines may have one or the other.
537
IEEE = IEEE .OR. LIEEE1
539
* Compute RMIN by successive division by BETA. We could compute
540
* RMIN as BASE**( EMIN - 1 ), but some machines underflow during
544
DO 30 I = 1, 1 - LEMIN
545
LRMIN = SLAMC3( LRMIN*RBASE, ZERO )
548
* Finally, call SLAMC5 to compute EMAX and RMAX.
550
CALL SLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
564
9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
566
$ ' If, after inspection, the value EMIN looks',
567
$ ' acceptable please comment out ',
568
$ / ' the IF block as marked within the code of routine',
569
$ ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / )
575
************************************************************************
577
REAL FUNCTION SLAMC3( A, B )
579
* -- LAPACK auxiliary routine (version 2.0) --
580
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
581
* Courant Institute, Argonne National Lab, and Rice University
584
* .. Scalar Arguments ..
591
* SLAMC3 is intended to force A and B to be stored prior to doing
592
* the addition of A and B , for use in situations where optimizers
593
* might hold one of these in a register.
599
* The values A and B.
601
* =====================================================================
603
* .. Executable Statements ..
613
************************************************************************
615
SUBROUTINE SLAMC4( EMIN, START, BASE )
617
* -- LAPACK auxiliary routine (version 2.0) --
618
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
619
* Courant Institute, Argonne National Lab, and Rice University
622
* .. Scalar Arguments ..
630
* SLAMC4 is a service routine for SLAMC2.
636
* The minimum exponent before (gradual) underflow, computed by
637
* setting A = START and dividing by BASE until the previous A
638
* can not be recovered.
641
* The starting point for determining EMIN.
643
* BASE (input) INTEGER
644
* The base of the machine.
646
* =====================================================================
648
* .. Local Scalars ..
650
REAL A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
652
* .. External Functions ..
656
* .. Executable Statements ..
663
B1 = SLAMC3( A*RBASE, ZERO )
668
*+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
669
* $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
671
IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND.
675
B1 = SLAMC3( A / BASE, ZERO )
676
C1 = SLAMC3( B1*BASE, ZERO )
681
B2 = SLAMC3( A*RBASE, ZERO )
682
C2 = SLAMC3( B2 / RBASE, ZERO )
697
************************************************************************
699
SUBROUTINE SLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
701
* -- LAPACK auxiliary routine (version 2.0) --
702
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
703
* Courant Institute, Argonne National Lab, and Rice University
706
* .. Scalar Arguments ..
708
INTEGER BETA, EMAX, EMIN, P
715
* SLAMC5 attempts to compute RMAX, the largest machine floating-point
716
* number, without overflow. It assumes that EMAX + abs(EMIN) sum
717
* approximately to a power of 2. It will fail on machines where this
718
* assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
719
* EMAX = 28718). It will also fail if the value supplied for EMIN is
720
* too large (i.e. too close to zero), probably with overflow.
725
* BETA (input) INTEGER
726
* The base of floating-point arithmetic.
729
* The number of base BETA digits in the mantissa of a
730
* floating-point value.
732
* EMIN (input) INTEGER
733
* The minimum exponent before (gradual) underflow.
735
* IEEE (input) LOGICAL
736
* A logical flag specifying whether or not the arithmetic
737
* system is thought to comply with the IEEE standard.
739
* EMAX (output) INTEGER
740
* The largest exponent before overflow
743
* The largest machine floating-point number.
745
* =====================================================================
749
PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 )
751
* .. Local Scalars ..
752
INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
753
REAL OLDY, RECBAS, Y, Z
755
* .. External Functions ..
759
* .. Intrinsic Functions ..
762
* .. Executable Statements ..
764
* First compute LEXP and UEXP, two powers of 2 that bound
765
* abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
766
* approximately to the bound that is closest to abs(EMIN).
767
* (EMAX is the exponent of the required number RMAX).
773
IF( TRY.LE.( -EMIN ) ) THEN
778
IF( LEXP.EQ.-EMIN ) THEN
785
* Now -LEXP is less than or equal to EMIN, and -UEXP is greater
786
* than or equal to EMIN. EXBITS is the number of bits needed to
787
* store the exponent.
789
IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
795
* EXPSUM is the exponent range, approximately equal to
798
EMAX = EXPSUM + EMIN - 1
799
NBITS = 1 + EXBITS + P
801
* NBITS is the total number of bits needed to store a
802
* floating-point number.
804
IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
806
* Either there are an odd number of bits used to store a
807
* floating-point number, which is unlikely, or some bits are
808
* not used in the representation of numbers, which is possible,
809
* (e.g. Cray machines) or the mantissa has an implicit bit,
810
* (e.g. IEEE machines, Dec Vax machines), which is perhaps the
811
* most likely. We have to assume the last alternative.
812
* If this is true, then we need to reduce EMAX by one because
813
* there must be some way of representing zero in an implicit-bit
814
* system. On machines like Cray, we are reducing EMAX by one
822
* Assume we are on an IEEE machine which reserves one exponent
823
* for infinity and NaN.
828
* Now create RMAX, the largest machine number, which should
829
* be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
831
* First compute 1.0 - BETA**(-P), being careful that the
832
* result is less than 1.0 .
846
* Now multiply by BETA**EMAX to get RMAX.
849
Y = SLAMC3( Y*BETA, ZERO )
added
removed
routines/lapack/slamch.f
