Mercurial > octave

comparison libcruft/lapack/slamc2.f @ 7789:82be108cc558

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First attempt at single precision tyeps * * * corrections to qrupdate single precision routines * * * prefer demotion to single over promotion to double * * * Add single precision support to log2 function * * * Trivial PROJECT file update * * * Cache optimized hermitian/transpose methods * * * Add tests for tranpose/hermitian and ChangeLog entry for new transpose code
author David Bateman <dbateman@free.fr>
date Sun, 27 Apr 2008 22:34:17 +0200
parents
children
comparison
equal deleted inserted replaced
7788:45f5faba05a2 7789:82be108cc558
1 SUBROUTINE SLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
2 *
3 * -- LAPACK auxiliary routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Scalar Arguments ..
8 LOGICAL RND
9 INTEGER BETA, EMAX, EMIN, T
10 REAL EPS, RMAX, RMIN
11 * ..
12 *
13 * Purpose
14 * =======
15 *
16 * SLAMC2 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) REAL
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) REAL
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) REAL
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 REAL A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
70 $ SIXTH, SMALL, THIRD, TWO, ZERO
71 * ..
72 * .. External Functions ..
73 REAL SLAMC3
74 EXTERNAL SLAMC3
75 * ..
76 * .. External Subroutines ..
77 EXTERNAL SLAMC1, SLAMC4, SLAMC5
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 SLAMC3 to ensure
100 * that relevant values are stored and not held in registers, or
101 * are not affected by optimizers.
102 *
103 * SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
104 *
105 CALL SLAMC1( 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 = SLAMC3( B, -HALF )
118 THIRD = SLAMC3( SIXTH, SIXTH )
119 B = SLAMC3( THIRD, -HALF )
120 B = SLAMC3( 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 = SLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
132 C = SLAMC3( HALF, -C )
133 B = SLAMC3( HALF, C )
134 C = SLAMC3( HALF, -B )
135 B = SLAMC3( 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 = SLAMC3( SMALL*RBASE, ZERO )
153 20 CONTINUE
154 A = SLAMC3( ONE, SMALL )
155 CALL SLAMC4( NGPMIN, ONE, LBETA )
156 CALL SLAMC4( NGNMIN, -ONE, LBETA )
157 CALL SLAMC4( GPMIN, A, LBETA )
158 CALL SLAMC4( 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
205 FIRST = .FALSE.
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 SLAMC1. 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 = SLAMC3( LRMIN*RBASE, ZERO )
228 30 CONTINUE
229 *
230 * Finally, call SLAMC5 to compute EMAX and RMAX.
231 *
232 CALL SLAMC5( 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 $ ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / )
252 *
253 * End of SLAMC2
254 *
255 END