|
2329
|
1 SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, 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 IEEE |
|
|
9 INTEGER BETA, EMAX, EMIN, P |
|
|
10 DOUBLE PRECISION RMAX |
|
|
11 * .. |
|
|
12 * |
|
|
13 * Purpose |
|
|
14 * ======= |
|
|
15 * |
|
|
16 * DLAMC5 attempts to compute RMAX, the largest machine floating-point |
|
|
17 * number, without overflow. It assumes that EMAX + abs(EMIN) sum |
|
|
18 * approximately to a power of 2. It will fail on machines where this |
|
|
19 * assumption does not hold, for example, the Cyber 205 (EMIN = -28625, |
|
|
20 * EMAX = 28718). It will also fail if the value supplied for EMIN is |
|
|
21 * too large (i.e. too close to zero), probably with overflow. |
|
|
22 * |
|
|
23 * Arguments |
|
|
24 * ========= |
|
|
25 * |
|
|
26 * BETA (input) INTEGER |
|
|
27 * The base of floating-point arithmetic. |
|
|
28 * |
|
|
29 * P (input) INTEGER |
|
|
30 * The number of base BETA digits in the mantissa of a |
|
|
31 * floating-point value. |
|
|
32 * |
|
|
33 * EMIN (input) INTEGER |
|
|
34 * The minimum exponent before (gradual) underflow. |
|
|
35 * |
|
|
36 * IEEE (input) LOGICAL |
|
|
37 * A logical flag specifying whether or not the arithmetic |
|
|
38 * system is thought to comply with the IEEE standard. |
|
|
39 * |
|
|
40 * EMAX (output) INTEGER |
|
|
41 * The largest exponent before overflow |
|
|
42 * |
|
|
43 * RMAX (output) DOUBLE PRECISION |
|
|
44 * The largest machine floating-point number. |
|
|
45 * |
|
|
46 * ===================================================================== |
|
|
47 * |
|
|
48 * .. Parameters .. |
|
|
49 DOUBLE PRECISION ZERO, ONE |
|
|
50 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) |
|
|
51 * .. |
|
|
52 * .. Local Scalars .. |
|
|
53 INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP |
|
|
54 DOUBLE PRECISION OLDY, RECBAS, Y, Z |
|
|
55 * .. |
|
|
56 * .. External Functions .. |
|
|
57 DOUBLE PRECISION DLAMC3 |
|
|
58 EXTERNAL DLAMC3 |
|
|
59 * .. |
|
|
60 * .. Intrinsic Functions .. |
|
|
61 INTRINSIC MOD |
|
|
62 * .. |
|
|
63 * .. Executable Statements .. |
|
|
64 * |
|
|
65 * First compute LEXP and UEXP, two powers of 2 that bound |
|
|
66 * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum |
|
|
67 * approximately to the bound that is closest to abs(EMIN). |
|
|
68 * (EMAX is the exponent of the required number RMAX). |
|
|
69 * |
|
|
70 LEXP = 1 |
|
|
71 EXBITS = 1 |
|
|
72 10 CONTINUE |
|
|
73 TRY = LEXP*2 |
|
|
74 IF( TRY.LE.( -EMIN ) ) THEN |
|
|
75 LEXP = TRY |
|
|
76 EXBITS = EXBITS + 1 |
|
|
77 GO TO 10 |
|
|
78 END IF |
|
|
79 IF( LEXP.EQ.-EMIN ) THEN |
|
|
80 UEXP = LEXP |
|
|
81 ELSE |
|
|
82 UEXP = TRY |
|
|
83 EXBITS = EXBITS + 1 |
|
|
84 END IF |
|
|
85 * |
|
|
86 * Now -LEXP is less than or equal to EMIN, and -UEXP is greater |
|
|
87 * than or equal to EMIN. EXBITS is the number of bits needed to |
|
|
88 * store the exponent. |
|
|
89 * |
|
|
90 IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN |
|
|
91 EXPSUM = 2*LEXP |
|
|
92 ELSE |
|
|
93 EXPSUM = 2*UEXP |
|
|
94 END IF |
|
|
95 * |
|
|
96 * EXPSUM is the exponent range, approximately equal to |
|
|
97 * EMAX - EMIN + 1 . |
|
|
98 * |
|
|
99 EMAX = EXPSUM + EMIN - 1 |
|
|
100 NBITS = 1 + EXBITS + P |
|
|
101 * |
|
|
102 * NBITS is the total number of bits needed to store a |
|
|
103 * floating-point number. |
|
|
104 * |
|
|
105 IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN |
|
|
106 * |
|
|
107 * Either there are an odd number of bits used to store a |
|
|
108 * floating-point number, which is unlikely, or some bits are |
|
|
109 * not used in the representation of numbers, which is possible, |
|
|
110 * (e.g. Cray machines) or the mantissa has an implicit bit, |
|
|
111 * (e.g. IEEE machines, Dec Vax machines), which is perhaps the |
|
|
112 * most likely. We have to assume the last alternative. |
|
|
113 * If this is true, then we need to reduce EMAX by one because |
|
|
114 * there must be some way of representing zero in an implicit-bit |
|
|
115 * system. On machines like Cray, we are reducing EMAX by one |
|
|
116 * unnecessarily. |
|
|
117 * |
|
|
118 EMAX = EMAX - 1 |
|
|
119 END IF |
|
|
120 * |
|
|
121 IF( IEEE ) THEN |
|
|
122 * |
|
|
123 * Assume we are on an IEEE machine which reserves one exponent |
|
|
124 * for infinity and NaN. |
|
|
125 * |
|
|
126 EMAX = EMAX - 1 |
|
|
127 END IF |
|
|
128 * |
|
|
129 * Now create RMAX, the largest machine number, which should |
|
|
130 * be equal to (1.0 - BETA**(-P)) * BETA**EMAX . |
|
|
131 * |
|
|
132 * First compute 1.0 - BETA**(-P), being careful that the |
|
|
133 * result is less than 1.0 . |
|
|
134 * |
|
|
135 RECBAS = ONE / BETA |
|
|
136 Z = BETA - ONE |
|
|
137 Y = ZERO |
|
|
138 DO 20 I = 1, P |
|
|
139 Z = Z*RECBAS |
|
|
140 IF( Y.LT.ONE ) |
|
|
141 $ OLDY = Y |
|
|
142 Y = DLAMC3( Y, Z ) |
|
|
143 20 CONTINUE |
|
|
144 IF( Y.GE.ONE ) |
|
|
145 $ Y = OLDY |
|
|
146 * |
|
|
147 * Now multiply by BETA**EMAX to get RMAX. |
|
|
148 * |
|
|
149 DO 30 I = 1, EMAX |
|
|
150 Y = DLAMC3( Y*BETA, ZERO ) |
|
|
151 30 CONTINUE |
|
|
152 * |
|
|
153 RMAX = Y |
|
|
154 RETURN |
|
|
155 * |
|
|
156 * End of DLAMC5 |
|
|
157 * |
|
|
158 END |
