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