Annotation of rpl/lapack/lapack/dlamch.f, revision 1.2
1.1 bertrand 1: DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
2: *
3: * -- LAPACK auxiliary routine (version 3.2) --
4: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5: * November 2006
6: *
7: * .. Scalar Arguments ..
8: CHARACTER CMACH
9: * ..
10: *
11: * Purpose
12: * =======
13: *
14: * DLAMCH determines double precision machine parameters.
15: *
16: * Arguments
17: * =========
18: *
19: * CMACH (input) CHARACTER*1
20: * Specifies the value to be returned by DLAMCH:
21: * = 'E' or 'e', DLAMCH := eps
22: * = 'S' or 's , DLAMCH := sfmin
23: * = 'B' or 'b', DLAMCH := base
24: * = 'P' or 'p', DLAMCH := eps*base
25: * = 'N' or 'n', DLAMCH := t
26: * = 'R' or 'r', DLAMCH := rnd
27: * = 'M' or 'm', DLAMCH := emin
28: * = 'U' or 'u', DLAMCH := rmin
29: * = 'L' or 'l', DLAMCH := emax
30: * = 'O' or 'o', DLAMCH := rmax
31: *
32: * where
33: *
34: * eps = relative machine precision
35: * sfmin = safe minimum, such that 1/sfmin does not overflow
36: * base = base of the machine
37: * prec = eps*base
38: * t = number of (base) digits in the mantissa
39: * rnd = 1.0 when rounding occurs in addition, 0.0 otherwise
40: * emin = minimum exponent before (gradual) underflow
41: * rmin = underflow threshold - base**(emin-1)
42: * emax = largest exponent before overflow
43: * rmax = overflow threshold - (base**emax)*(1-eps)
44: *
45: * =====================================================================
46: *
47: * .. Parameters ..
48: DOUBLE PRECISION ONE, ZERO
49: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
50: * ..
51: * .. Local Scalars ..
52: LOGICAL FIRST, LRND
53: INTEGER BETA, IMAX, IMIN, IT
54: DOUBLE PRECISION BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,
55: $ RND, SFMIN, SMALL, T
56: * ..
57: * .. External Functions ..
58: LOGICAL LSAME
59: EXTERNAL LSAME
60: * ..
61: * .. External Subroutines ..
62: EXTERNAL DLAMC2
63: * ..
64: * .. Save statement ..
65: SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,
66: $ EMAX, RMAX, PREC
67: * ..
68: * .. Data statements ..
69: DATA FIRST / .TRUE. /
70: * ..
71: * .. Executable Statements ..
72: *
73: IF( FIRST ) THEN
74: CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )
75: BASE = BETA
76: T = IT
77: IF( LRND ) THEN
78: RND = ONE
79: EPS = ( BASE**( 1-IT ) ) / 2
80: ELSE
81: RND = ZERO
82: EPS = BASE**( 1-IT )
83: END IF
84: PREC = EPS*BASE
85: EMIN = IMIN
86: EMAX = IMAX
87: SFMIN = RMIN
88: SMALL = ONE / RMAX
89: IF( SMALL.GE.SFMIN ) THEN
90: *
91: * Use SMALL plus a bit, to avoid the possibility of rounding
92: * causing overflow when computing 1/sfmin.
93: *
94: SFMIN = SMALL*( ONE+EPS )
95: END IF
96: END IF
97: *
98: IF( LSAME( CMACH, 'E' ) ) THEN
99: RMACH = EPS
100: ELSE IF( LSAME( CMACH, 'S' ) ) THEN
101: RMACH = SFMIN
102: ELSE IF( LSAME( CMACH, 'B' ) ) THEN
103: RMACH = BASE
104: ELSE IF( LSAME( CMACH, 'P' ) ) THEN
105: RMACH = PREC
106: ELSE IF( LSAME( CMACH, 'N' ) ) THEN
107: RMACH = T
108: ELSE IF( LSAME( CMACH, 'R' ) ) THEN
109: RMACH = RND
110: ELSE IF( LSAME( CMACH, 'M' ) ) THEN
111: RMACH = EMIN
112: ELSE IF( LSAME( CMACH, 'U' ) ) THEN
113: RMACH = RMIN
114: ELSE IF( LSAME( CMACH, 'L' ) ) THEN
115: RMACH = EMAX
116: ELSE IF( LSAME( CMACH, 'O' ) ) THEN
117: RMACH = RMAX
118: END IF
119: *
120: DLAMCH = RMACH
121: FIRST = .FALSE.
122: RETURN
123: *
124: * End of DLAMCH
125: *
126: END
127: *
128: ************************************************************************
129: *
130: SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 )
131: *
132: * -- LAPACK auxiliary routine (version 3.2) --
133: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
134: * November 2006
135: *
136: * .. Scalar Arguments ..
137: LOGICAL IEEE1, RND
138: INTEGER BETA, T
139: * ..
140: *
141: * Purpose
142: * =======
143: *
144: * DLAMC1 determines the machine parameters given by BETA, T, RND, and
145: * IEEE1.
146: *
147: * Arguments
148: * =========
149: *
150: * BETA (output) INTEGER
151: * The base of the machine.
152: *
153: * T (output) INTEGER
154: * The number of ( BETA ) digits in the mantissa.
155: *
156: * RND (output) LOGICAL
157: * Specifies whether proper rounding ( RND = .TRUE. ) or
158: * chopping ( RND = .FALSE. ) occurs in addition. This may not
159: * be a reliable guide to the way in which the machine performs
160: * its arithmetic.
161: *
162: * IEEE1 (output) LOGICAL
163: * Specifies whether rounding appears to be done in the IEEE
164: * 'round to nearest' style.
165: *
166: * Further Details
167: * ===============
168: *
169: * The routine is based on the routine ENVRON by Malcolm and
170: * incorporates suggestions by Gentleman and Marovich. See
171: *
172: * Malcolm M. A. (1972) Algorithms to reveal properties of
173: * floating-point arithmetic. Comms. of the ACM, 15, 949-951.
174: *
175: * Gentleman W. M. and Marovich S. B. (1974) More on algorithms
176: * that reveal properties of floating point arithmetic units.
177: * Comms. of the ACM, 17, 276-277.
178: *
179: * =====================================================================
180: *
181: * .. Local Scalars ..
182: LOGICAL FIRST, LIEEE1, LRND
183: INTEGER LBETA, LT
184: DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2
185: * ..
186: * .. External Functions ..
187: DOUBLE PRECISION DLAMC3
188: EXTERNAL DLAMC3
189: * ..
190: * .. Save statement ..
191: SAVE FIRST, LIEEE1, LBETA, LRND, LT
192: * ..
193: * .. Data statements ..
194: DATA FIRST / .TRUE. /
195: * ..
196: * .. Executable Statements ..
197: *
198: IF( FIRST ) THEN
199: ONE = 1
200: *
201: * LBETA, LIEEE1, LT and LRND are the local values of BETA,
202: * IEEE1, T and RND.
203: *
204: * Throughout this routine we use the function DLAMC3 to ensure
205: * that relevant values are stored and not held in registers, or
206: * are not affected by optimizers.
207: *
208: * Compute a = 2.0**m with the smallest positive integer m such
209: * that
210: *
211: * fl( a + 1.0 ) = a.
212: *
213: A = 1
214: C = 1
215: *
216: *+ WHILE( C.EQ.ONE )LOOP
217: 10 CONTINUE
218: IF( C.EQ.ONE ) THEN
219: A = 2*A
220: C = DLAMC3( A, ONE )
221: C = DLAMC3( C, -A )
222: GO TO 10
223: END IF
224: *+ END WHILE
225: *
226: * Now compute b = 2.0**m with the smallest positive integer m
227: * such that
228: *
229: * fl( a + b ) .gt. a.
230: *
231: B = 1
232: C = DLAMC3( A, B )
233: *
234: *+ WHILE( C.EQ.A )LOOP
235: 20 CONTINUE
236: IF( C.EQ.A ) THEN
237: B = 2*B
238: C = DLAMC3( A, B )
239: GO TO 20
240: END IF
241: *+ END WHILE
242: *
243: * Now compute the base. a and c are neighbouring floating point
244: * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
245: * their difference is beta. Adding 0.25 to c is to ensure that it
246: * is truncated to beta and not ( beta - 1 ).
247: *
248: QTR = ONE / 4
249: SAVEC = C
250: C = DLAMC3( C, -A )
251: LBETA = C + QTR
252: *
253: * Now determine whether rounding or chopping occurs, by adding a
254: * bit less than beta/2 and a bit more than beta/2 to a.
255: *
256: B = LBETA
257: F = DLAMC3( B / 2, -B / 100 )
258: C = DLAMC3( F, A )
259: IF( C.EQ.A ) THEN
260: LRND = .TRUE.
261: ELSE
262: LRND = .FALSE.
263: END IF
264: F = DLAMC3( B / 2, B / 100 )
265: C = DLAMC3( F, A )
266: IF( ( LRND ) .AND. ( C.EQ.A ) )
267: $ LRND = .FALSE.
268: *
269: * Try and decide whether rounding is done in the IEEE 'round to
270: * nearest' style. B/2 is half a unit in the last place of the two
271: * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
272: * zero, and SAVEC is odd. Thus adding B/2 to A should not change
273: * A, but adding B/2 to SAVEC should change SAVEC.
274: *
275: T1 = DLAMC3( B / 2, A )
276: T2 = DLAMC3( B / 2, SAVEC )
277: LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
278: *
279: * Now find the mantissa, t. It should be the integer part of
280: * log to the base beta of a, however it is safer to determine t
281: * by powering. So we find t as the smallest positive integer for
282: * which
283: *
284: * fl( beta**t + 1.0 ) = 1.0.
285: *
286: LT = 0
287: A = 1
288: C = 1
289: *
290: *+ WHILE( C.EQ.ONE )LOOP
291: 30 CONTINUE
292: IF( C.EQ.ONE ) THEN
293: LT = LT + 1
294: A = A*LBETA
295: C = DLAMC3( A, ONE )
296: C = DLAMC3( C, -A )
297: GO TO 30
298: END IF
299: *+ END WHILE
300: *
301: END IF
302: *
303: BETA = LBETA
304: T = LT
305: RND = LRND
306: IEEE1 = LIEEE1
307: FIRST = .FALSE.
308: RETURN
309: *
310: * End of DLAMC1
311: *
312: END
313: *
314: ************************************************************************
315: *
316: SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
317: *
318: * -- LAPACK auxiliary routine (version 3.2) --
319: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
320: * November 2006
321: *
322: * .. Scalar Arguments ..
323: LOGICAL RND
324: INTEGER BETA, EMAX, EMIN, T
325: DOUBLE PRECISION EPS, RMAX, RMIN
326: * ..
327: *
328: * Purpose
329: * =======
330: *
331: * DLAMC2 determines the machine parameters specified in its argument
332: * list.
333: *
334: * Arguments
335: * =========
336: *
337: * BETA (output) INTEGER
338: * The base of the machine.
339: *
340: * T (output) INTEGER
341: * The number of ( BETA ) digits in the mantissa.
342: *
343: * RND (output) LOGICAL
344: * Specifies whether proper rounding ( RND = .TRUE. ) or
345: * chopping ( RND = .FALSE. ) occurs in addition. This may not
346: * be a reliable guide to the way in which the machine performs
347: * its arithmetic.
348: *
349: * EPS (output) DOUBLE PRECISION
350: * The smallest positive number such that
351: *
352: * fl( 1.0 - EPS ) .LT. 1.0,
353: *
354: * where fl denotes the computed value.
355: *
356: * EMIN (output) INTEGER
357: * The minimum exponent before (gradual) underflow occurs.
358: *
359: * RMIN (output) DOUBLE PRECISION
360: * The smallest normalized number for the machine, given by
361: * BASE**( EMIN - 1 ), where BASE is the floating point value
362: * of BETA.
363: *
364: * EMAX (output) INTEGER
365: * The maximum exponent before overflow occurs.
366: *
367: * RMAX (output) DOUBLE PRECISION
368: * The largest positive number for the machine, given by
369: * BASE**EMAX * ( 1 - EPS ), where BASE is the floating point
370: * value of BETA.
371: *
372: * Further Details
373: * ===============
374: *
375: * The computation of EPS is based on a routine PARANOIA by
376: * W. Kahan of the University of California at Berkeley.
377: *
378: * =====================================================================
379: *
380: * .. Local Scalars ..
381: LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND
382: INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
383: $ NGNMIN, NGPMIN
384: DOUBLE PRECISION A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
385: $ SIXTH, SMALL, THIRD, TWO, ZERO
386: * ..
387: * .. External Functions ..
388: DOUBLE PRECISION DLAMC3
389: EXTERNAL DLAMC3
390: * ..
391: * .. External Subroutines ..
392: EXTERNAL DLAMC1, DLAMC4, DLAMC5
393: * ..
394: * .. Intrinsic Functions ..
395: INTRINSIC ABS, MAX, MIN
396: * ..
397: * .. Save statement ..
398: SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
399: $ LRMIN, LT
400: * ..
401: * .. Data statements ..
402: DATA FIRST / .TRUE. / , IWARN / .FALSE. /
403: * ..
404: * .. Executable Statements ..
405: *
406: IF( FIRST ) THEN
407: ZERO = 0
408: ONE = 1
409: TWO = 2
410: *
411: * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of
412: * BETA, T, RND, EPS, EMIN and RMIN.
413: *
414: * Throughout this routine we use the function DLAMC3 to ensure
415: * that relevant values are stored and not held in registers, or
416: * are not affected by optimizers.
417: *
418: * DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1.
419: *
420: CALL DLAMC1( LBETA, LT, LRND, LIEEE1 )
421: *
422: * Start to find EPS.
423: *
424: B = LBETA
425: A = B**( -LT )
426: LEPS = A
427: *
428: * Try some tricks to see whether or not this is the correct EPS.
429: *
430: B = TWO / 3
431: HALF = ONE / 2
432: SIXTH = DLAMC3( B, -HALF )
433: THIRD = DLAMC3( SIXTH, SIXTH )
434: B = DLAMC3( THIRD, -HALF )
435: B = DLAMC3( B, SIXTH )
436: B = ABS( B )
437: IF( B.LT.LEPS )
438: $ B = LEPS
439: *
440: LEPS = 1
441: *
442: *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
443: 10 CONTINUE
444: IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
445: LEPS = B
446: C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
447: C = DLAMC3( HALF, -C )
448: B = DLAMC3( HALF, C )
449: C = DLAMC3( HALF, -B )
450: B = DLAMC3( HALF, C )
451: GO TO 10
452: END IF
453: *+ END WHILE
454: *
455: IF( A.LT.LEPS )
456: $ LEPS = A
457: *
458: * Computation of EPS complete.
459: *
460: * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)).
461: * Keep dividing A by BETA until (gradual) underflow occurs. This
462: * is detected when we cannot recover the previous A.
463: *
464: RBASE = ONE / LBETA
465: SMALL = ONE
466: DO 20 I = 1, 3
467: SMALL = DLAMC3( SMALL*RBASE, ZERO )
468: 20 CONTINUE
469: A = DLAMC3( ONE, SMALL )
470: CALL DLAMC4( NGPMIN, ONE, LBETA )
471: CALL DLAMC4( NGNMIN, -ONE, LBETA )
472: CALL DLAMC4( GPMIN, A, LBETA )
473: CALL DLAMC4( GNMIN, -A, LBETA )
474: IEEE = .FALSE.
475: *
476: IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
477: IF( NGPMIN.EQ.GPMIN ) THEN
478: LEMIN = NGPMIN
479: * ( Non twos-complement machines, no gradual underflow;
480: * e.g., VAX )
481: ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
482: LEMIN = NGPMIN - 1 + LT
483: IEEE = .TRUE.
484: * ( Non twos-complement machines, with gradual underflow;
485: * e.g., IEEE standard followers )
486: ELSE
487: LEMIN = MIN( NGPMIN, GPMIN )
488: * ( A guess; no known machine )
489: IWARN = .TRUE.
490: END IF
491: *
492: ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
493: IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
494: LEMIN = MAX( NGPMIN, NGNMIN )
495: * ( Twos-complement machines, no gradual underflow;
496: * e.g., CYBER 205 )
497: ELSE
498: LEMIN = MIN( NGPMIN, NGNMIN )
499: * ( A guess; no known machine )
500: IWARN = .TRUE.
501: END IF
502: *
503: ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
504: $ ( GPMIN.EQ.GNMIN ) ) THEN
505: IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
506: LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
507: * ( Twos-complement machines with gradual underflow;
508: * no known machine )
509: ELSE
510: LEMIN = MIN( NGPMIN, NGNMIN )
511: * ( A guess; no known machine )
512: IWARN = .TRUE.
513: END IF
514: *
515: ELSE
516: LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
517: * ( A guess; no known machine )
518: IWARN = .TRUE.
519: END IF
520: FIRST = .FALSE.
521: ***
522: * Comment out this if block if EMIN is ok
523: IF( IWARN ) THEN
524: FIRST = .TRUE.
525: WRITE( 6, FMT = 9999 )LEMIN
526: END IF
527: ***
528: *
529: * Assume IEEE arithmetic if we found denormalised numbers above,
530: * or if arithmetic seems to round in the IEEE style, determined
531: * in routine DLAMC1. A true IEEE machine should have both things
532: * true; however, faulty machines may have one or the other.
533: *
534: IEEE = IEEE .OR. LIEEE1
535: *
536: * Compute RMIN by successive division by BETA. We could compute
537: * RMIN as BASE**( EMIN - 1 ), but some machines underflow during
538: * this computation.
539: *
540: LRMIN = 1
541: DO 30 I = 1, 1 - LEMIN
542: LRMIN = DLAMC3( LRMIN*RBASE, ZERO )
543: 30 CONTINUE
544: *
545: * Finally, call DLAMC5 to compute EMAX and RMAX.
546: *
547: CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
548: END IF
549: *
550: BETA = LBETA
551: T = LT
552: RND = LRND
553: EPS = LEPS
554: EMIN = LEMIN
555: RMIN = LRMIN
556: EMAX = LEMAX
557: RMAX = LRMAX
558: *
559: RETURN
560: *
561: 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
562: $ ' EMIN = ', I8, /
563: $ ' If, after inspection, the value EMIN looks',
564: $ ' acceptable please comment out ',
565: $ / ' the IF block as marked within the code of routine',
566: $ ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / )
567: *
568: * End of DLAMC2
569: *
570: END
571: *
572: ************************************************************************
573: *
574: DOUBLE PRECISION FUNCTION DLAMC3( A, B )
575: *
576: * -- LAPACK auxiliary routine (version 3.2) --
577: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
578: * November 2006
579: *
580: * .. Scalar Arguments ..
581: DOUBLE PRECISION A, B
582: * ..
583: *
584: * Purpose
585: * =======
586: *
587: * DLAMC3 is intended to force A and B to be stored prior to doing
588: * the addition of A and B , for use in situations where optimizers
589: * might hold one of these in a register.
590: *
591: * Arguments
592: * =========
593: *
594: * A (input) DOUBLE PRECISION
595: * B (input) DOUBLE PRECISION
596: * The values A and B.
597: *
598: * =====================================================================
599: *
600: * .. Executable Statements ..
601: *
602: DLAMC3 = A + B
603: *
604: RETURN
605: *
606: * End of DLAMC3
607: *
608: END
609: *
610: ************************************************************************
611: *
612: SUBROUTINE DLAMC4( EMIN, START, BASE )
613: *
614: * -- LAPACK auxiliary routine (version 3.2) --
615: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
616: * November 2006
617: *
618: * .. Scalar Arguments ..
619: INTEGER BASE, EMIN
620: DOUBLE PRECISION START
621: * ..
622: *
623: * Purpose
624: * =======
625: *
626: * DLAMC4 is a service routine for DLAMC2.
627: *
628: * Arguments
629: * =========
630: *
631: * EMIN (output) INTEGER
632: * The minimum exponent before (gradual) underflow, computed by
633: * setting A = START and dividing by BASE until the previous A
634: * can not be recovered.
635: *
636: * START (input) DOUBLE PRECISION
637: * The starting point for determining EMIN.
638: *
639: * BASE (input) INTEGER
640: * The base of the machine.
641: *
642: * =====================================================================
643: *
644: * .. Local Scalars ..
645: INTEGER I
646: DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO
647: * ..
648: * .. External Functions ..
649: DOUBLE PRECISION DLAMC3
650: EXTERNAL DLAMC3
651: * ..
652: * .. Executable Statements ..
653: *
654: A = START
655: ONE = 1
656: RBASE = ONE / BASE
657: ZERO = 0
658: EMIN = 1
659: B1 = DLAMC3( A*RBASE, ZERO )
660: C1 = A
661: C2 = A
662: D1 = A
663: D2 = A
664: *+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.
665: * $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP
666: 10 CONTINUE
667: IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND.
668: $ ( D2.EQ.A ) ) THEN
669: EMIN = EMIN - 1
670: A = B1
671: B1 = DLAMC3( A / BASE, ZERO )
672: C1 = DLAMC3( B1*BASE, ZERO )
673: D1 = ZERO
674: DO 20 I = 1, BASE
675: D1 = D1 + B1
676: 20 CONTINUE
677: B2 = DLAMC3( A*RBASE, ZERO )
678: C2 = DLAMC3( B2 / RBASE, ZERO )
679: D2 = ZERO
680: DO 30 I = 1, BASE
681: D2 = D2 + B2
682: 30 CONTINUE
683: GO TO 10
684: END IF
685: *+ END WHILE
686: *
687: RETURN
688: *
689: * End of DLAMC4
690: *
691: END
692: *
693: ************************************************************************
694: *
695: SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )
696: *
697: * -- LAPACK auxiliary routine (version 3.2) --
698: * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
699: * November 2006
700: *
701: * .. Scalar Arguments ..
702: LOGICAL IEEE
703: INTEGER BETA, EMAX, EMIN, P
704: DOUBLE PRECISION RMAX
705: * ..
706: *
707: * Purpose
708: * =======
709: *
710: * DLAMC5 attempts to compute RMAX, the largest machine floating-point
711: * number, without overflow. It assumes that EMAX + abs(EMIN) sum
712: * approximately to a power of 2. It will fail on machines where this
713: * assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
714: * EMAX = 28718). It will also fail if the value supplied for EMIN is
715: * too large (i.e. too close to zero), probably with overflow.
716: *
717: * Arguments
718: * =========
719: *
720: * BETA (input) INTEGER
721: * The base of floating-point arithmetic.
722: *
723: * P (input) INTEGER
724: * The number of base BETA digits in the mantissa of a
725: * floating-point value.
726: *
727: * EMIN (input) INTEGER
728: * The minimum exponent before (gradual) underflow.
729: *
730: * IEEE (input) LOGICAL
731: * A logical flag specifying whether or not the arithmetic
732: * system is thought to comply with the IEEE standard.
733: *
734: * EMAX (output) INTEGER
735: * The largest exponent before overflow
736: *
737: * RMAX (output) DOUBLE PRECISION
738: * The largest machine floating-point number.
739: *
740: * =====================================================================
741: *
742: * .. Parameters ..
743: DOUBLE PRECISION ZERO, ONE
744: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
745: * ..
746: * .. Local Scalars ..
747: INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
748: DOUBLE PRECISION OLDY, RECBAS, Y, Z
749: * ..
750: * .. External Functions ..
751: DOUBLE PRECISION DLAMC3
752: EXTERNAL DLAMC3
753: * ..
754: * .. Intrinsic Functions ..
755: INTRINSIC MOD
756: * ..
757: * .. Executable Statements ..
758: *
759: * First compute LEXP and UEXP, two powers of 2 that bound
760: * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
761: * approximately to the bound that is closest to abs(EMIN).
762: * (EMAX is the exponent of the required number RMAX).
763: *
764: LEXP = 1
765: EXBITS = 1
766: 10 CONTINUE
767: TRY = LEXP*2
768: IF( TRY.LE.( -EMIN ) ) THEN
769: LEXP = TRY
770: EXBITS = EXBITS + 1
771: GO TO 10
772: END IF
773: IF( LEXP.EQ.-EMIN ) THEN
774: UEXP = LEXP
775: ELSE
776: UEXP = TRY
777: EXBITS = EXBITS + 1
778: END IF
779: *
780: * Now -LEXP is less than or equal to EMIN, and -UEXP is greater
781: * than or equal to EMIN. EXBITS is the number of bits needed to
782: * store the exponent.
783: *
784: IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
785: EXPSUM = 2*LEXP
786: ELSE
787: EXPSUM = 2*UEXP
788: END IF
789: *
790: * EXPSUM is the exponent range, approximately equal to
791: * EMAX - EMIN + 1 .
792: *
793: EMAX = EXPSUM + EMIN - 1
794: NBITS = 1 + EXBITS + P
795: *
796: * NBITS is the total number of bits needed to store a
797: * floating-point number.
798: *
799: IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
800: *
801: * Either there are an odd number of bits used to store a
802: * floating-point number, which is unlikely, or some bits are
803: * not used in the representation of numbers, which is possible,
804: * (e.g. Cray machines) or the mantissa has an implicit bit,
805: * (e.g. IEEE machines, Dec Vax machines), which is perhaps the
806: * most likely. We have to assume the last alternative.
807: * If this is true, then we need to reduce EMAX by one because
808: * there must be some way of representing zero in an implicit-bit
809: * system. On machines like Cray, we are reducing EMAX by one
810: * unnecessarily.
811: *
812: EMAX = EMAX - 1
813: END IF
814: *
815: IF( IEEE ) THEN
816: *
817: * Assume we are on an IEEE machine which reserves one exponent
818: * for infinity and NaN.
819: *
820: EMAX = EMAX - 1
821: END IF
822: *
823: * Now create RMAX, the largest machine number, which should
824: * be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
825: *
826: * First compute 1.0 - BETA**(-P), being careful that the
827: * result is less than 1.0 .
828: *
829: RECBAS = ONE / BETA
830: Z = BETA - ONE
831: Y = ZERO
832: DO 20 I = 1, P
833: Z = Z*RECBAS
834: IF( Y.LT.ONE )
835: $ OLDY = Y
836: Y = DLAMC3( Y, Z )
837: 20 CONTINUE
838: IF( Y.GE.ONE )
839: $ Y = OLDY
840: *
841: * Now multiply by BETA**EMAX to get RMAX.
842: *
843: DO 30 I = 1, EMAX
844: Y = DLAMC3( Y*BETA, ZERO )
845: 30 CONTINUE
846: *
847: RMAX = Y
848: RETURN
849: *
850: * End of DLAMC5
851: *
852: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlamch.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack