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