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