|
2329
|
1 SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 ) |
|
|
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 IEEE1, RND |
|
|
9 INTEGER BETA, T |
|
|
10 * .. |
|
|
11 * |
|
|
12 * Purpose |
|
|
13 * ======= |
|
|
14 * |
|
|
15 * DLAMC1 determines the machine parameters given by BETA, T, RND, and |
|
|
16 * IEEE1. |
|
|
17 * |
|
|
18 * Arguments |
|
|
19 * ========= |
|
|
20 * |
|
|
21 * BETA (output) INTEGER |
|
|
22 * The base of the machine. |
|
|
23 * |
|
|
24 * T (output) INTEGER |
|
|
25 * The number of ( BETA ) digits in the mantissa. |
|
|
26 * |
|
|
27 * RND (output) LOGICAL |
|
|
28 * Specifies whether proper rounding ( RND = .TRUE. ) or |
|
|
29 * chopping ( RND = .FALSE. ) occurs in addition. This may not |
|
|
30 * be a reliable guide to the way in which the machine performs |
|
|
31 * its arithmetic. |
|
|
32 * |
|
|
33 * IEEE1 (output) LOGICAL |
|
|
34 * Specifies whether rounding appears to be done in the IEEE |
|
|
35 * 'round to nearest' style. |
|
|
36 * |
|
|
37 * Further Details |
|
|
38 * =============== |
|
|
39 * |
|
|
40 * The routine is based on the routine ENVRON by Malcolm and |
|
|
41 * incorporates suggestions by Gentleman and Marovich. See |
|
|
42 * |
|
|
43 * Malcolm M. A. (1972) Algorithms to reveal properties of |
|
|
44 * floating-point arithmetic. Comms. of the ACM, 15, 949-951. |
|
|
45 * |
|
|
46 * Gentleman W. M. and Marovich S. B. (1974) More on algorithms |
|
|
47 * that reveal properties of floating point arithmetic units. |
|
|
48 * Comms. of the ACM, 17, 276-277. |
|
|
49 * |
|
|
50 * ===================================================================== |
|
|
51 * |
|
|
52 * .. Local Scalars .. |
|
|
53 LOGICAL FIRST, LIEEE1, LRND |
|
|
54 INTEGER LBETA, LT |
|
|
55 DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2 |
|
|
56 * .. |
|
|
57 * .. External Functions .. |
|
|
58 DOUBLE PRECISION DLAMC3 |
|
|
59 EXTERNAL DLAMC3 |
|
|
60 * .. |
|
|
61 * .. Save statement .. |
|
|
62 SAVE FIRST, LIEEE1, LBETA, LRND, LT |
|
|
63 * .. |
|
|
64 * .. Data statements .. |
|
|
65 DATA FIRST / .TRUE. / |
|
|
66 * .. |
|
|
67 * .. Executable Statements .. |
|
|
68 * |
|
|
69 IF( FIRST ) THEN |
|
|
70 ONE = 1 |
|
|
71 * |
|
|
72 * LBETA, LIEEE1, LT and LRND are the local values of BETA, |
|
|
73 * IEEE1, T and RND. |
|
|
74 * |
|
|
75 * Throughout this routine we use the function DLAMC3 to ensure |
|
|
76 * that relevant values are stored and not held in registers, or |
|
|
77 * are not affected by optimizers. |
|
|
78 * |
|
|
79 * Compute a = 2.0**m with the smallest positive integer m such |
|
|
80 * that |
|
|
81 * |
|
|
82 * fl( a + 1.0 ) = a. |
|
|
83 * |
|
|
84 A = 1 |
|
|
85 C = 1 |
|
|
86 * |
|
|
87 *+ WHILE( C.EQ.ONE )LOOP |
|
|
88 10 CONTINUE |
|
|
89 IF( C.EQ.ONE ) THEN |
|
|
90 A = 2*A |
|
|
91 C = DLAMC3( A, ONE ) |
|
|
92 C = DLAMC3( C, -A ) |
|
|
93 GO TO 10 |
|
|
94 END IF |
|
|
95 *+ END WHILE |
|
|
96 * |
|
|
97 * Now compute b = 2.0**m with the smallest positive integer m |
|
|
98 * such that |
|
|
99 * |
|
|
100 * fl( a + b ) .gt. a. |
|
|
101 * |
|
|
102 B = 1 |
|
|
103 C = DLAMC3( A, B ) |
|
|
104 * |
|
|
105 *+ WHILE( C.EQ.A )LOOP |
|
|
106 20 CONTINUE |
|
|
107 IF( C.EQ.A ) THEN |
|
|
108 B = 2*B |
|
|
109 C = DLAMC3( A, B ) |
|
|
110 GO TO 20 |
|
|
111 END IF |
|
|
112 *+ END WHILE |
|
|
113 * |
|
|
114 * Now compute the base. a and c are neighbouring floating point |
|
|
115 * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so |
|
|
116 * their difference is beta. Adding 0.25 to c is to ensure that it |
|
|
117 * is truncated to beta and not ( beta - 1 ). |
|
|
118 * |
|
|
119 QTR = ONE / 4 |
|
|
120 SAVEC = C |
|
|
121 C = DLAMC3( C, -A ) |
|
|
122 LBETA = C + QTR |
|
|
123 * |
|
|
124 * Now determine whether rounding or chopping occurs, by adding a |
|
|
125 * bit less than beta/2 and a bit more than beta/2 to a. |
|
|
126 * |
|
|
127 B = LBETA |
|
|
128 F = DLAMC3( B / 2, -B / 100 ) |
|
|
129 C = DLAMC3( F, A ) |
|
|
130 IF( C.EQ.A ) THEN |
|
|
131 LRND = .TRUE. |
|
|
132 ELSE |
|
|
133 LRND = .FALSE. |
|
|
134 END IF |
|
|
135 F = DLAMC3( B / 2, B / 100 ) |
|
|
136 C = DLAMC3( F, A ) |
|
|
137 IF( ( LRND ) .AND. ( C.EQ.A ) ) |
|
|
138 $ LRND = .FALSE. |
|
|
139 * |
|
|
140 * Try and decide whether rounding is done in the IEEE 'round to |
|
|
141 * nearest' style. B/2 is half a unit in the last place of the two |
|
|
142 * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit |
|
|
143 * zero, and SAVEC is odd. Thus adding B/2 to A should not change |
|
|
144 * A, but adding B/2 to SAVEC should change SAVEC. |
|
|
145 * |
|
|
146 T1 = DLAMC3( B / 2, A ) |
|
|
147 T2 = DLAMC3( B / 2, SAVEC ) |
|
|
148 LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND |
|
|
149 * |
|
|
150 * Now find the mantissa, t. It should be the integer part of |
|
|
151 * log to the base beta of a, however it is safer to determine t |
|
|
152 * by powering. So we find t as the smallest positive integer for |
|
|
153 * which |
|
|
154 * |
|
|
155 * fl( beta**t + 1.0 ) = 1.0. |
|
|
156 * |
|
|
157 LT = 0 |
|
|
158 A = 1 |
|
|
159 C = 1 |
|
|
160 * |
|
|
161 *+ WHILE( C.EQ.ONE )LOOP |
|
|
162 30 CONTINUE |
|
|
163 IF( C.EQ.ONE ) THEN |
|
|
164 LT = LT + 1 |
|
|
165 A = A*LBETA |
|
|
166 C = DLAMC3( A, ONE ) |
|
|
167 C = DLAMC3( C, -A ) |
|
|
168 GO TO 30 |
|
|
169 END IF |
|
|
170 *+ END WHILE |
|
|
171 * |
|
|
172 END IF |
|
|
173 * |
|
|
174 BETA = LBETA |
|
|
175 T = LT |
|
|
176 RND = LRND |
|
|
177 IEEE1 = LIEEE1 |
|
7034
|
178 FIRST = .FALSE. |
|
2329
|
179 RETURN |
|
|
180 * |
|
|
181 * End of DLAMC1 |
|
|
182 * |
|
|
183 END |
