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