SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlamc1()

subroutine dlamc1 ( integer  beta,
integer  t,
logical  rnd,
logical  ieee1 
)

Definition at line 139 of file tools.f.

140*
141* -- LAPACK auxiliary routine (version 2.0) --
142* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
143* Courant Institute, Argonne National Lab, and Rice University
144* October 31, 1992
145*
146* .. Scalar Arguments ..
147 LOGICAL IEEE1, RND
148 INTEGER BETA, T
149* ..
150*
151* Purpose
152* =======
153*
154* DLAMC1 determines the machine parameters given by BETA, T, RND, and
155* IEEE1.
156*
157* Arguments
158* =========
159*
160* BETA (output) INTEGER
161* The base of the machine.
162*
163* T (output) INTEGER
164* The number of ( BETA ) digits in the mantissa.
165*
166* RND (output) LOGICAL
167* Specifies whether proper rounding ( RND = .TRUE. ) or
168* chopping ( RND = .FALSE. ) occurs in addition. This may not
169* be a reliable guide to the way in which the machine performs
170* its arithmetic.
171*
172* IEEE1 (output) LOGICAL
173* Specifies whether rounding appears to be done in the IEEE
174* 'round to nearest' style.
175*
176* Further Details
177* ===============
178*
179* The routine is based on the routine ENVRON by Malcolm and
180* incorporates suggestions by Gentleman and Marovich. See
181*
182* Malcolm M. A. (1972) Algorithms to reveal properties of
183* floating-point arithmetic. Comms. of the ACM, 15, 949-951.
184*
185* Gentleman W. M. and Marovich S. B. (1974) More on algorithms
186* that reveal properties of floating point arithmetic units.
187* Comms. of the ACM, 17, 276-277.
188*
189* =====================================================================
190*
191* .. Local Scalars ..
192 LOGICAL FIRST, LIEEE1, LRND
193 INTEGER LBETA, LT
194 DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2
195* ..
196* .. External Functions ..
197 DOUBLE PRECISION DLAMC3
198 EXTERNAL dlamc3
199* ..
200* .. Save statement ..
201 SAVE first, lieee1, lbeta, lrnd, lt
202* ..
203* .. Data statements ..
204 DATA first / .true. /
205* ..
206* .. Executable Statements ..
207*
208 IF( first ) THEN
209 first = .false.
210 one = 1
211*
212* LBETA, LIEEE1, LT and LRND are the local values of BETA,
213* IEEE1, T and RND.
214*
215* Throughout this routine we use the function DLAMC3 to ensure
216* that relevant values are stored and not held in registers, or
217* are not affected by optimizers.
218*
219* Compute a = 2.0**m with the smallest positive integer m such
220* that
221*
222* fl( a + 1.0 ) = a.
223*
224 a = 1
225 c = 1
226*
227*+ WHILE( C.EQ.ONE )LOOP
228 10 CONTINUE
229 IF( c.EQ.one ) THEN
230 a = 2*a
231 c = dlamc3( a, one )
232 c = dlamc3( c, -a )
233 GO TO 10
234 END IF
235*+ END WHILE
236*
237* Now compute b = 2.0**m with the smallest positive integer m
238* such that
239*
240* fl( a + b ) .gt. a.
241*
242 b = 1
243 c = dlamc3( a, b )
244*
245*+ WHILE( C.EQ.A )LOOP
246 20 CONTINUE
247 IF( c.EQ.a ) THEN
248 b = 2*b
249 c = dlamc3( a, b )
250 GO TO 20
251 END IF
252*+ END WHILE
253*
254* Now compute the base. a and c are neighbouring floating point
255* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
256* their difference is beta. Adding 0.25 to c is to ensure that it
257* is truncated to beta and not ( beta - 1 ).
258*
259 qtr = one / 4
260 savec = c
261 c = dlamc3( c, -a )
262 lbeta = c + qtr
263*
264* Now determine whether rounding or chopping occurs, by adding a
265* bit less than beta/2 and a bit more than beta/2 to a.
266*
267 b = lbeta
268 f = dlamc3( b / 2, -b / 100 )
269 c = dlamc3( f, a )
270 IF( c.EQ.a ) THEN
271 lrnd = .true.
272 ELSE
273 lrnd = .false.
274 END IF
275 f = dlamc3( b / 2, b / 100 )
276 c = dlamc3( f, a )
277 IF( ( lrnd ) .AND. ( c.EQ.a ) )
278 $ lrnd = .false.
279*
280* Try and decide whether rounding is done in the IEEE 'round to
281* nearest' style. B/2 is half a unit in the last place of the two
282* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
283* zero, and SAVEC is odd. Thus adding B/2 to A should not change
284* A, but adding B/2 to SAVEC should change SAVEC.
285*
286 t1 = dlamc3( b / 2, a )
287 t2 = dlamc3( b / 2, savec )
288 lieee1 = ( t1.EQ.a ) .AND. ( t2.GT.savec ) .AND. lrnd
289*
290* Now find the mantissa, t. It should be the integer part of
291* log to the base beta of a, however it is safer to determine t
292* by powering. So we find t as the smallest positive integer for
293* which
294*
295* fl( beta**t + 1.0 ) = 1.0.
296*
297 lt = 0
298 a = 1
299 c = 1
300*
301*+ WHILE( C.EQ.ONE )LOOP
302 30 CONTINUE
303 IF( c.EQ.one ) THEN
304 lt = lt + 1
305 a = a*lbeta
306 c = dlamc3( a, one )
307 c = dlamc3( c, -a )
308 GO TO 30
309 END IF
310*+ END WHILE
311*
312 END IF
313*
314 beta = lbeta
315 t = lt
316 rnd = lrnd
317 ieee1 = lieee1
318 RETURN
319*
320* End of DLAMC1
321*
double precision function dlamc3(a, b)
Definition tools.f:586
Here is the caller graph for this function: