|
2329
|
1 SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX ) |
|
|
2 * |
|
3339
|
3 * -- LAPACK auxiliary routine (version 3.0) -- |
|
2329
|
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
|
|
5 * Courant Institute, Argonne National Lab, and Rice University |
|
|
6 * October 31, 1992 |
|
|
7 * |
|
|
8 * .. Scalar Arguments .. |
|
|
9 LOGICAL RND |
|
|
10 INTEGER BETA, EMAX, EMIN, T |
|
|
11 DOUBLE PRECISION EPS, RMAX, RMIN |
|
|
12 * .. |
|
|
13 * |
|
|
14 * Purpose |
|
|
15 * ======= |
|
|
16 * |
|
|
17 * DLAMC2 determines the machine parameters specified in its argument |
|
|
18 * list. |
|
|
19 * |
|
|
20 * Arguments |
|
|
21 * ========= |
|
|
22 * |
|
|
23 * BETA (output) INTEGER |
|
|
24 * The base of the machine. |
|
|
25 * |
|
|
26 * T (output) INTEGER |
|
|
27 * The number of ( BETA ) digits in the mantissa. |
|
|
28 * |
|
|
29 * RND (output) LOGICAL |
|
|
30 * Specifies whether proper rounding ( RND = .TRUE. ) or |
|
|
31 * chopping ( RND = .FALSE. ) occurs in addition. This may not |
|
|
32 * be a reliable guide to the way in which the machine performs |
|
|
33 * its arithmetic. |
|
|
34 * |
|
|
35 * EPS (output) DOUBLE PRECISION |
|
|
36 * The smallest positive number such that |
|
|
37 * |
|
|
38 * fl( 1.0 - EPS ) .LT. 1.0, |
|
|
39 * |
|
|
40 * where fl denotes the computed value. |
|
|
41 * |
|
|
42 * EMIN (output) INTEGER |
|
|
43 * The minimum exponent before (gradual) underflow occurs. |
|
|
44 * |
|
|
45 * RMIN (output) DOUBLE PRECISION |
|
|
46 * The smallest normalized number for the machine, given by |
|
|
47 * BASE**( EMIN - 1 ), where BASE is the floating point value |
|
|
48 * of BETA. |
|
|
49 * |
|
|
50 * EMAX (output) INTEGER |
|
|
51 * The maximum exponent before overflow occurs. |
|
|
52 * |
|
|
53 * RMAX (output) DOUBLE PRECISION |
|
|
54 * The largest positive number for the machine, given by |
|
|
55 * BASE**EMAX * ( 1 - EPS ), where BASE is the floating point |
|
|
56 * value of BETA. |
|
|
57 * |
|
|
58 * Further Details |
|
|
59 * =============== |
|
|
60 * |
|
|
61 * The computation of EPS is based on a routine PARANOIA by |
|
|
62 * W. Kahan of the University of California at Berkeley. |
|
|
63 * |
|
|
64 * ===================================================================== |
|
|
65 * |
|
|
66 * .. Local Scalars .. |
|
|
67 LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND |
|
|
68 INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT, |
|
|
69 $ NGNMIN, NGPMIN |
|
|
70 DOUBLE PRECISION A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE, |
|
|
71 $ SIXTH, SMALL, THIRD, TWO, ZERO |
|
|
72 * .. |
|
|
73 * .. External Functions .. |
|
|
74 DOUBLE PRECISION DLAMC3 |
|
|
75 EXTERNAL DLAMC3 |
|
|
76 * .. |
|
|
77 * .. External Subroutines .. |
|
|
78 EXTERNAL DLAMC1, DLAMC4, DLAMC5 |
|
|
79 * .. |
|
|
80 * .. Intrinsic Functions .. |
|
|
81 INTRINSIC ABS, MAX, MIN |
|
|
82 * .. |
|
|
83 * .. Save statement .. |
|
|
84 SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX, |
|
|
85 $ LRMIN, LT |
|
|
86 * .. |
|
|
87 * .. Data statements .. |
|
|
88 DATA FIRST / .TRUE. / , IWARN / .FALSE. / |
|
|
89 * .. |
|
|
90 * .. Executable Statements .. |
|
|
91 * |
|
|
92 IF( FIRST ) THEN |
|
|
93 FIRST = .FALSE. |
|
|
94 ZERO = 0 |
|
|
95 ONE = 1 |
|
|
96 TWO = 2 |
|
|
97 * |
|
|
98 * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of |
|
|
99 * BETA, T, RND, EPS, EMIN and RMIN. |
|
|
100 * |
|
|
101 * Throughout this routine we use the function DLAMC3 to ensure |
|
|
102 * that relevant values are stored and not held in registers, or |
|
|
103 * are not affected by optimizers. |
|
|
104 * |
|
|
105 * DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. |
|
|
106 * |
|
|
107 CALL DLAMC1( LBETA, LT, LRND, LIEEE1 ) |
|
|
108 * |
|
|
109 * Start to find EPS. |
|
|
110 * |
|
|
111 B = LBETA |
|
|
112 A = B**( -LT ) |
|
|
113 LEPS = A |
|
|
114 * |
|
|
115 * Try some tricks to see whether or not this is the correct EPS. |
|
|
116 * |
|
|
117 B = TWO / 3 |
|
|
118 HALF = ONE / 2 |
|
|
119 SIXTH = DLAMC3( B, -HALF ) |
|
|
120 THIRD = DLAMC3( SIXTH, SIXTH ) |
|
|
121 B = DLAMC3( THIRD, -HALF ) |
|
|
122 B = DLAMC3( B, SIXTH ) |
|
|
123 B = ABS( B ) |
|
|
124 IF( B.LT.LEPS ) |
|
|
125 $ B = LEPS |
|
|
126 * |
|
|
127 LEPS = 1 |
|
|
128 * |
|
|
129 *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP |
|
|
130 10 CONTINUE |
|
|
131 IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN |
|
|
132 LEPS = B |
|
|
133 C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) |
|
|
134 C = DLAMC3( HALF, -C ) |
|
|
135 B = DLAMC3( HALF, C ) |
|
|
136 C = DLAMC3( HALF, -B ) |
|
|
137 B = DLAMC3( HALF, C ) |
|
|
138 GO TO 10 |
|
|
139 END IF |
|
|
140 *+ END WHILE |
|
|
141 * |
|
|
142 IF( A.LT.LEPS ) |
|
|
143 $ LEPS = A |
|
|
144 * |
|
|
145 * Computation of EPS complete. |
|
|
146 * |
|
|
147 * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). |
|
|
148 * Keep dividing A by BETA until (gradual) underflow occurs. This |
|
|
149 * is detected when we cannot recover the previous A. |
|
|
150 * |
|
|
151 RBASE = ONE / LBETA |
|
|
152 SMALL = ONE |
|
|
153 DO 20 I = 1, 3 |
|
|
154 SMALL = DLAMC3( SMALL*RBASE, ZERO ) |
|
|
155 20 CONTINUE |
|
|
156 A = DLAMC3( ONE, SMALL ) |
|
|
157 CALL DLAMC4( NGPMIN, ONE, LBETA ) |
|
|
158 CALL DLAMC4( NGNMIN, -ONE, LBETA ) |
|
|
159 CALL DLAMC4( GPMIN, A, LBETA ) |
|
|
160 CALL DLAMC4( GNMIN, -A, LBETA ) |
|
|
161 IEEE = .FALSE. |
|
|
162 * |
|
|
163 IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN |
|
|
164 IF( NGPMIN.EQ.GPMIN ) THEN |
|
|
165 LEMIN = NGPMIN |
|
|
166 * ( Non twos-complement machines, no gradual underflow; |
|
|
167 * e.g., VAX ) |
|
|
168 ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN |
|
|
169 LEMIN = NGPMIN - 1 + LT |
|
|
170 IEEE = .TRUE. |
|
|
171 * ( Non twos-complement machines, with gradual underflow; |
|
|
172 * e.g., IEEE standard followers ) |
|
|
173 ELSE |
|
|
174 LEMIN = MIN( NGPMIN, GPMIN ) |
|
|
175 * ( A guess; no known machine ) |
|
|
176 IWARN = .TRUE. |
|
|
177 END IF |
|
|
178 * |
|
|
179 ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN |
|
|
180 IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN |
|
|
181 LEMIN = MAX( NGPMIN, NGNMIN ) |
|
|
182 * ( Twos-complement machines, no gradual underflow; |
|
|
183 * e.g., CYBER 205 ) |
|
|
184 ELSE |
|
|
185 LEMIN = MIN( NGPMIN, NGNMIN ) |
|
|
186 * ( A guess; no known machine ) |
|
|
187 IWARN = .TRUE. |
|
|
188 END IF |
|
|
189 * |
|
|
190 ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND. |
|
|
191 $ ( GPMIN.EQ.GNMIN ) ) THEN |
|
|
192 IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN |
|
|
193 LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT |
|
|
194 * ( Twos-complement machines with gradual underflow; |
|
|
195 * no known machine ) |
|
|
196 ELSE |
|
|
197 LEMIN = MIN( NGPMIN, NGNMIN ) |
|
|
198 * ( A guess; no known machine ) |
|
|
199 IWARN = .TRUE. |
|
|
200 END IF |
|
|
201 * |
|
|
202 ELSE |
|
|
203 LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN ) |
|
|
204 * ( A guess; no known machine ) |
|
|
205 IWARN = .TRUE. |
|
|
206 END IF |
|
|
207 *** |
|
|
208 * Comment out this if block if EMIN is ok |
|
|
209 IF( IWARN ) THEN |
|
|
210 FIRST = .TRUE. |
|
|
211 WRITE( 6, FMT = 9999 )LEMIN |
|
|
212 END IF |
|
|
213 *** |
|
|
214 * |
|
|
215 * Assume IEEE arithmetic if we found denormalised numbers above, |
|
|
216 * or if arithmetic seems to round in the IEEE style, determined |
|
|
217 * in routine DLAMC1. A true IEEE machine should have both things |
|
|
218 * true; however, faulty machines may have one or the other. |
|
|
219 * |
|
|
220 IEEE = IEEE .OR. LIEEE1 |
|
|
221 * |
|
|
222 * Compute RMIN by successive division by BETA. We could compute |
|
|
223 * RMIN as BASE**( EMIN - 1 ), but some machines underflow during |
|
|
224 * this computation. |
|
|
225 * |
|
|
226 LRMIN = 1 |
|
|
227 DO 30 I = 1, 1 - LEMIN |
|
|
228 LRMIN = DLAMC3( LRMIN*RBASE, ZERO ) |
|
|
229 30 CONTINUE |
|
|
230 * |
|
|
231 * Finally, call DLAMC5 to compute EMAX and RMAX. |
|
|
232 * |
|
|
233 CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX ) |
|
|
234 END IF |
|
|
235 * |
|
|
236 BETA = LBETA |
|
|
237 T = LT |
|
|
238 RND = LRND |
|
|
239 EPS = LEPS |
|
|
240 EMIN = LEMIN |
|
|
241 RMIN = LRMIN |
|
|
242 EMAX = LEMAX |
|
|
243 RMAX = LRMAX |
|
|
244 * |
|
|
245 RETURN |
|
|
246 * |
|
|
247 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-', |
|
|
248 $ ' EMIN = ', I8, / |
|
|
249 $ ' If, after inspection, the value EMIN looks', |
|
|
250 $ ' acceptable please comment out ', |
|
|
251 $ / ' the IF block as marked within the code of routine', |
|
|
252 $ ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / ) |
|
|
253 * |
|
|
254 * End of DLAMC2 |
|
|
255 * |
|
|
256 END |
