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