1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
|
*
************************************************************************
*
SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
*
* -- LAPACK auxiliary routine (version 2.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* October 31, 1992
*
* .. 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
FIRST = .FALSE.
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
***
* 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
|
