1      REAL             FUNCTION SLAMCH( CMACH )
2*
3*  -- LAPACK auxiliary routine (version 1.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 1.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 1.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 1.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 1.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 1.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
858