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